Arbitrary Boundary Conditions and Constraints in Quantum Algorithms for Differential Equations via Penalty Projections

Al\'an Aspuru-Guzik; Jin-Peng Liu; Nathan Wiebe; Philipp Schleich; Tyler Kharazi; Xiangyu Li

REVIEW 3 major objections 4 minor 1 cited by

Arbitrary Boundary Conditions and Constraints in Quantum Algorithms for Differential Equations via Penalty Projections

T0 review · 3 major / 4 minor · reviewed 2026-08-06 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read The paper claims that adding a penalty projection $-i\lambda P_c$ to the generator of a differential equation enforces arbitrary boundary conditions in quantum algorithms, with the gate-cost overhead growing only as $O(\log \lambda)$ when…

desk verdict A useful penalty-projection framework for constraints in quantum DE solvers with an appealing O(log λ) overhead, but the time-dependent error proofs—especially Lemma 13—need a rigorous rewrite before the central guarantees are accepted. read the letter →

arxiv 2506.21751 v1 pith:RN3ZFNBB submitted 2025-06-26 quant-ph

Philipp Schleich , Tyler Kharazi , Xiangyu Li , Jin-Peng Liu , Alán Aspuru-Guzik , Nathan Wiebe This is my paper

classification quant-ph MSC 81P6868Q1265M06 PACS 03.67.Ac

keywords penaltyprojectionboundaryconditionsquantumdifferentialequationslinearcombinationofHamiltoniansimulationinteractionpicturefast-forwardingDirichletandNeumannconstraineddynamics

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Quantum algorithms for differential equations have mostly sidestepped boundary conditions, assuming periodic or trivial boundaries. This paper proposes to enforce arbitrary constraints—Dirichlet values, Neumann derivatives, Robin combinations, and interface conditions—by adding a penalty term $-i\lambda P_c$ to the dynamics that damps every component of the solution living on the infeasible subspace. It proves that the leftover constraint error shrinks like $O(1/\lambda)$, so any target accuracy $\varepsilon$ can be reached by choosing $\lambda$ large enough, with explicit lower bounds such as $\lambda \ge 2 v_{\max}^2 \|A_0\|/\varepsilon$ for stable homogeneous systems. The decisive claim is that the gate-cost overhead of the penalty is at most $O(\log \lambda)$, not linear in $\lambda$, provided the projection can be fast-forwarded—which the paper demonstrates for Dirichlet and Neumann projections built from phase gates and swaps. Numerical experiments on the heat and wave equations confirm the predicted scaling, and the heat-equation example puts the total overhead at $\widetilde{O}(d \log n + \log t)$.

What carries the argument

The load-bearing object is the fast-forwardable penalty projection $P_c$, added to the generator as $-i\lambda P_c$. Dirichlet conditions are enforced by $P_c^D = \sum_{j \in \Gamma_D} |j\rangle\langle j|$, a phase gate on boundary basis states; Neumann conditions use the swap-based projector $P_c^N = \frac{1}{2}(I - S)$ built from a finite-difference stencil, whose Hamiltonian simulation $e^{i\lambda t P} = P^\perp + e^{i\lambda t} P$ has cost independent of $\lambda t$. The error analysis runs through a non-Hermitian, time-dependent, inhomogeneous generalization of the Kubo formula that bounds the expectation $\langle P_c \rangle$, i.e. the infeasible-space norm, in terms of the modified anticommutator $\{P_c, A_0\}_{\sim}$ and the $L^1$ norm of the forcing term. The quantum algorithm then solves the interaction-picture ODE $v_I(t) = e^{i\lambda P_c t} v(t)$ with the Linear Combination of Hamiltonian Simulation solver, whose query complexity is governed by the norm of the conjugated generator (independent of $\lambda$) and whose discretization steps grow only logarithmically with $\lambda$.

What would settle it

Measure the gate count of the proposed fast-forwarded projection circuits as a function of $\lambda$: if the cost of $e^{-i\lambda P_c t}$ grows with $\lambda$ for any boundary geometry satisfying the paper's assumptions, the central claim fails. More directly, solve the penalized ODE for a fixed stable $A_0$ and boundary, and verify the predicted scaling $\|P_c v(t)\|^2 = O(\|A_0\| v_{\max}^2 / \lambda)$ and the threshold $\lambda \ge 2 v_{\max}^2 \|A_0\|/\varepsilon$; a violation of the $O(1/\lambda)$ decay for large $\lambda$ would falsify the error analysis. A third concrete check: construct a boundary set where constraint points have overlapping neighbor sets, violating conditions (105) and (106), and test whether the simulation cost still grows only logarithmically in $\lambda$.

Watch

Extended reading notes

Core claim

The central claim is that the constrained evolution $\frac{\mathrm{d}}{\mathrm{d}t} v(t) = (A_0(t) - i\lambda P_c) v(t) + b(t)$ with a Hermitian projection $P_c$ onto the constraint-violating subspace is a faithful surrogate for the constrained problem $P_c v = 0$: the norm of the solution on the infeasible subspace obeys $\|P_c v(t)\|^2 \le \varepsilon$ provided $\lambda$ is chosen at least $\sim 2\|A_0\|(v_{\max}^2 + 2 v_{\max} t B + t^2 B^2)/\varepsilon$ for stable dynamics with forcing bounded by $B$, and the solution on the feasible subspace is exactly unchanged by the penalty. The algorithmic consequence is that simulating the penalized system in the interaction picture of $P_c$ makes the cost of enforcing the constraint grow only as $O(\log \lambda)$ in the penalty strength, because the interaction-picture generator $e^{i\lambda P_c t} A_0 e^{-i\lambda P_c t}$ has norm independent of $\lambda$, and the transformation $e^{-i\lambda P_c t}$ itself is fast-forwardable for the Dirichlet (phase-gate) and Neumann (swap-based) projections the paper constructs. For the heat equation this yields a gate-complexity overhead of $\widetilde{O}(d \log n + \log t)$ over the unconstrained solver.

Load-bearing premise

The entire $O(\log \lambda)$ efficiency claim rests on the assumption that the evolution $e^{-i\lambda P_c t}$ of the constraint projection can be fast-forwarded, meaning its circuit cost does not grow with $\lambda t$; the paper proves this only for the specific Dirichlet and Neumann projections it constructs, and for a truly arbitrary constraint a user must supply a fast-forwardable projection or the logarithmic overhead disappears.

Editorial extensions

If this is right

Any constraint representable by an orthogonal, fast-forwardable projection—Dirichlet values, Neumann derivatives, Robin combinations, interface conditions—can be imposed inside an LCHS-based quantum DE solver at $O(\log \lambda)$ gate overhead, instead of the linear-in-$\lambda$ cost a direct encoding of the penalized generator would incur.
For the discretized heat equation, the gate-cost overhead of enforcing boundaries is $\widetilde{O}(d \log n + \log t)$, so boundary treatment does not change the solver's scaling with system size.
The explicit lower bounds on $\lambda$ (e.g. $\lambda \ge 2 v_{\max}^2 \|A_0\|/\varepsilon$ for stable homogeneous dynamics) let a user fix the penalty strength in advance from the initial data, the operator norm, and the target accuracy $\varepsilon$.
The feasible-subspace solution is exactly invariant under the penalty (Lemma 8), so enforcing constraints never perturbs the unconstrained dynamics on the admissible part of the state.
Non-zero Dirichlet and Neumann data are reduced to homogeneous constraints by shifting the solution, so the method covers realistic boundary data, not only zero boundaries.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

The practical reach of the method is the class of constraints whose projector has a fast-forwardable circuit; for any constraint where the cost of $e^{-i\lambda P_c t}$ grows with $\lambda$, the advertised logarithmic overhead reverts to polynomial in $\lambda$, so a fast-forwardability check should accompany any deployment outside the constructed Dirichlet and Neumann cases.
Reading the penalty as a high-frequency measurement ties the $O(1/\lambda)$ error bound to quantum-Zeno-type suppression; this suggests a measurement-based variant of the scheme in which repeated projections replace the coherent interaction-picture evolution, at the price of a success probability that scales with the constraint accuracy.
The numerical results show the actual constraint error several orders of magnitude below the first-order bound; if a tighter (e.g. second-order) analysis of the $1/\lambda$ expansion were available, the required $\lambda$—and therefore the practical gate count—could drop substantially.
Because the Neumann projector is built from swaps, the machinery should extend to constraints on non-local subspaces such as global conservation laws or symmetry sectors, provided the violating-subspace projector admits a fast-forwardable decomposition—an extension the paper only gestures at.

Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

Short version: worth engaging for the penalty-projection idea and the interaction-picture cost analysis; the numerics are honest and the code is out. But the analytic core, especially Lemma 13, has proof gaps that need fixing before the time-dependent inhomogeneous claims are trusted.

The genuinely new thing here is treating boundary conditions as a penalty projection -iλP_c in the generator and simulating in the interaction frame of that projection, with fast-forwardable constructions for Dirichlet, Neumann, and Robin conditions. That combination is new relative to complex absorbing potentials and previous LCHS work, and the O(log λ) overhead story is plausible and interesting. The construction of P_c via SWAP networks for derivative constraints is concrete, and the numerical experiments back the qualitative behavior, with error decaying roughly as λ^{-1/2} or a bit faster. Code is available, which is real evidence.

Soft spots: the proofs in Section II B are rougher than the rest. Lemma 9 contains an inequality as written that is false: it claims exp(2||Re(A0)||t) ≤ 1. The intended bound actually follows from the logarithmic norm since Re(A) ≤ 0, so this is fixable, but as printed it needs correction. Lemma 12 mixes commutator and anticommutator notation and needs a careful rewrite. More serious: Lemma 13 applies the generalized Kubo formula, which is only first-order with an O(ζ²) remainder, and the proof bounds the first-order response without estimating that remainder. Since ζ = 1/λ and λ is taken large, the remainder should be suppressible, but a rigorous λ lower bound needs it controlled explicitly. As it stands, the central ε-guarantee for time-dependent A0 with forcing is not fully proven, and it feeds Corollary 14 and the complexity theorem. I don't think the approach is wrong—the simpler lemmas and the numerics support the mechanism—but the strongest claim in the abstract rests on unproven asymptotics.

Who it's for: people working on quantum DE solvers, LCHS, and boundary-condition input models. It deserves a serious referee, with the expectation of major revision: make Lemma 13's remainder explicit, clean up Lemma 12, and fix the Lemma 9 inequality. I would cite this after those repairs.

Referee Report

3 major / 4 minor

Summary. The paper proposes enforcing boundary conditions and general constraints in quantum algorithms for ODEs/PDEs by adding a penalty projection -iλP_c to the dynamical generator. It derives lower bounds on λ for the constraint error (the norm of the solution on the infeasible subspace) under constant, time-dependent, and inhomogeneous dynamics; constructs fast-forwardable circuit implementations for Dirichlet, Neumann, Robin, and interface projections; presents numerical experiments for the heat and wave equations; and analyzes the gate complexity of an LCHS-based interaction-picture solver, claiming at most O(log λ) overhead for enforcing constraints.

Significance. If the analytical guarantees held, the proposed method would be a useful and modular way to impose boundary conditions in LCHS-based quantum differential-equation solvers, with a clearly identified assumption (fast-forwardability of the projection) and explicit circuit constructions. The paper also provides numerical evidence and a public code repository, which are assets. However, the central mathematical claims are not currently established: the most general error bound (Lemma 13) ignores a remainder term, and the feasible-space equivalence (Lemma 8) appears false as stated. The practical value of the approach therefore depends on repairing these proofs.

major comments (3)

[II.B.2, Lemma 13 and Eq. (72)] The proof of Lemma 13 bounds only the first-order term coming from Proposition 6 and never estimates the O(ζ²)=O(λ⁻²) remainder in Eq. (C1)/(59). Consequently the stated lower bound on λ does not by itself imply ∥P_c v(t)∥² ≤ ε for time-dependent inhomogeneous dynamics. This is load-bearing because Corollary 14, Table Ia, and Theorem 17 use Lemma 13 to state the ε-guarantee and complexity for non-zero, time-dependent, and inhomogeneous cases. In addition, the proof of Proposition 6 in Appendix C relies on the identity T_t T_{-t}^† = I for non-Hermitian time-dependent generators, which is not true in general. Please provide a complete remainder estimate, or replace the Kubo-formula argument with a direct Dyson-series bound that controls all orders.
[II.B, Lemma 8 and Eq. (18)] Lemma 8 is false as stated. The Taylor-expansion step P⊥_c e^{(A-iλP_c)t} = Σ_k t^k/k! (P⊥_c A - iλ P⊥_c P_c)^k is invalid because P⊥_c (A-iλP_c)^k is not equal to (P⊥_c(A-iλP_c))^k. For example, with A0 = [[-1,1],[1,-1]], P_c = diag(0,1), and v(0)=e1, the feasible components of e^{A0t}v(0) and e^{(A0-iλP_c)t}v(0) differ for finite λ and t>0. The later lemmas only bound ∥P_c v(t)∥, so the paper does not currently establish that the solution inside the feasible space is accurate. Since the method's usefulness requires both small constraint violation and small interior error, this lemma needs to be replaced by a valid convergence or error estimate.
[II.B.1, Lemma 9 proof, Eq. (23)] The proof of Lemma 9 contains a false inequality: it states ∥exp(At)^†∥∥exp(At)∥ ≤ exp(2∥Re(A0)∥t) ≤ 1, but ∥Re(A0)∥ is nonnegative, so the first exponential is at least 1 and the second inequality is generally false. The intended conclusion ∥e^{At}∥≤1 does follow from Re(A0)⪯0 via the log-norm bound, so the lemma is salvageable, but the written argument is invalid and should be corrected.

minor comments (4)

[Definitions and Problem 5] Definition 3 says 'P_c projects onto the constrained space and P_c to its complement'; the second occurrence should presumably be P⊥_c. Also, Eq. (6)-(7) contain a doubled comma after 't≥0'.
[II.D.3, Eq. (126)] The object P_R^c defined in Eq. (126) is a Hermitian linear combination of a value projector and a derivative projector, not generally an orthogonal projection. The paper's definitions and Lemma 15 are stated for orthogonal projections, so the Robin construction needs an explicit clarification of which properties of P_R^c are used.
[Appendix B, Eqs. (160)-(161)] The binomial identity for derivatives of e^{Bt}A(t)e^{-Bt} assumes B is constant in time; this is true in the paper's setting, but the assumption should be stated explicitly in the appendix.
[References] References [31] and [34] appear to cite the same work (Mangin-Brinet et al., Quantum 8, 1311 (2024)); please consolidate to avoid duplication.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity: the penalty error bounds are derived analytically, the fast-forwardable projection is an explicit assumption, and the complexity overhead is a substitution of the derived lambda; remaining proof gaps are correctness issues, not circular reductions.

full rationale

The paper's derivation chain is self-contained against its stated inputs. The constraint error is defined as the norm of the solution on the infeasible subspace, ||P_c v(t)||, and the bounds in Lemmas 9, 10, 12, and 13 are analytic inequalities obtained from the penalized ODE dynamics, with lambda chosen to make the derived upper bound smaller than epsilon. The numerical experiments in Section III validate these bounds but do not fit them, so there is no fitted-input-called-prediction pattern. The O(log lambda) gate overhead is obtained by substituting the derived lambda thresholds into the LCHS smoothness parameters Lambda_I and Xi_I in Section IV B 2; this is a substitution of an analytically derived quantity, not a renaming of the desired result. The fast-forwardability of the projection is explicitly assumed: Section IV states, 'A necessary condition of this to work efficiently is that the cost of the interaction picture transformation itself, i.e., e^{-i lambda P_c t}, does not depend on lambda,' and the Dirichlet and Neumann constructions in Section II D demonstrate that the assumption is satisfiable for those cases. The self-citations to [11] (co-author Liu) and [46] (co-author Wiebe) are load-bearing for the LCHS solver and interaction-picture simulation, but they are published, parameter-free theorems with stated assumptions that do not include the present paper's target result, so they are independent evidence rather than circularity. Two mathematical gaps should be flagged as correctness risks, not circularity: Lemma 13 in Section II B 2 drops the O(zeta^2) remainder of Proposition 6 when passing from Eq. (60) to Eq. (72), so the epsilon-guarantee for time-dependent inhomogeneous dynamics is not fully proven; and Lemma 9's estimate ||exp(At)^dagger exp(At)|| <= exp(2||Re(A0)||t) <= 1 is unjustified as written. Neither gap makes a claimed prediction equivalent to an input by construction, so the circularity score is 0.

Assumptions & free parameters 1 free parameters · 6 assumptions · 0 invented entities

The central claim depends on standard quantum algorithm assumptions (oracle access, amplitude encoding), the stability of the unconstrained dynamics, and fast-forwardability of the constraint projection. No new physical entities are introduced. The only tunable parameter is λ, which is set by the desired accuracy and is not fit to data.

free parameters (1)

penalty strength λ = not fitted; chosen as λ ≥ 2 v_max^2 ||A0|| / ε (Lemma 9), with modifications in Lemma 10-13
λ is a tunable algorithmic parameter introduced to enforce constraints. The paper provides lower bounds for target error ε; it is not fitted to data but is a user-chosen parameter that enters the gate complexity.

assumptions (6)

domain assumption The constraint projection P_c is an orthogonal (Hermitian) projection onto the infeasible subspace.
Definition 3 and throughout; used in all error lemmas. Non-Hermitian or non-orthogonal projections are explicitly left to future work (Section II A).
domain assumption The unconstrained dynamics are stable: Re(A0) ⪯ 0.
Assumed in Lemmas 9, 10, 12, 13 and in the complexity analysis. Lemma 11 relaxes this for finite time at the cost of an exponential factor.
domain assumption The constraint projection P_c is fast-forwardable, meaning e^{-iλP_c t} can be simulated with cost independent of λt.
Stated in Section IV and Section II D. This is required for the O(log λ) overhead; if it fails, the main complexity claim collapses.
domain assumption Initial data satisfies the constraints exactly, or can be prepared with small overlap ∥P_c v(0)∥ ≤ ζ (Lemma 16).
Problem 5 and Remark 4; the error bounds assume P_c v(0) = 0. Lemma 16 gives a procedure to prepare such states with O(1/ζ) cost.
domain assumption The Dirichlet and Neumann constraint projections commute where needed (I_ΓD ∩ I_ΓN = ∅).
Remark 2 and Definition 1. This ensures orthogonality and simplifies the combined Hamiltonian simulation in Eq. (129).
domain assumption The dynamics and forcing term are smooth and have bounded higher derivatives, and the forcing is bounded in L1 norm.
Used in the smoothness factor analysis in Section IV B and in Lemmas 10 and 13.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Arbitrary Boundary Conditions and Constraints in Quantum Algorithms for Differential Equations via Penalty Projections." pith.science (2026). https://pith.science/paper/RN3ZFNBB

@misc{pith2026250621751,
  author       = {Pith},
  title        = {Pith review of: Arbitrary Boundary Conditions and Constraints in Quantum Algorithms for Differential Equations via Penalty Projections},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/RN3ZFNBB}},
  note         = {Machine review of arXiv:2506.21751}
}

abstract

Complicated boundary conditions are essential to accurately describe phenomena arising in nature and engineering. Recently, the investigation of a potential speedup through quantum algorithms in simulating the governing ordinary and partial differential equations of such phenomena has gained increasing attention. We design an efficient quantum algorithms for solving differential equations with arbitrary boundary conditions. Specifically, we propose an approach to enforce arbitrary boundary conditions and constraints through adding a penalty projection to the governing equations. Assuming a fast-forwardable representation of the projection to ensure an efficient interaction picture imulation, the cost of to enforce the constraints is at most $O(\log\lambda)$ in the strength of the penalty $\lambda$ in the gate complexity; in the worst case, this goes as $O([\|v(0)\|^2\|A_0\| + \|b\|_{L^1[0;t]}^2)]t^2/\varepsilon)$, for precision $\varepsilon$ and a dynamical system $\frac{\rm d}{{\rm d}t} v(t) = A_0(t) v(t) + b(t)$ with negative semidefinite $A_0(t)$ of size $n^d\times n^d$. E.g., for the heat equation, this leads to a gate complexity overhead of $\widetilde O(d\log n + \log t)$. We show constraint error bounds for the penalty approach and provide validating numerical experiments, and estimate the circuit complexity using the Linear Combination of Hamiltonian Simulation.

Figures

Figures reproduced from arXiv: 2506.21751 by the authors.

**Figure 2.** FIG. 2: Visualization of relevant index sets as a graph depicted as a uniform grid. Recall that “inner” points [PITH_FULL_IMAGE:figures/full_fig_p007_2.png] view at source ↗

**Figure 3.** FIG. 3: Inducing a constraint subspace via a penalty projection [PITH_FULL_IMAGE:figures/full_fig_p008_3.png] view at source ↗

Figures from the paper (6 more)

**Figure 4.** Figure 4: FIG. 4: Conflicts arising in derivative constraints on a regular grid (left). Mitigation is possible by either [PITH_FULL_IMAGE:figures/full_fig_p028_4.png]

**Figure 5.** Figure 5: FIG. 5: Numerical simulations of heat equation with [PITH_FULL_IMAGE:figures/full_fig_p035_5.png]

**Figure 6.** Figure 6: FIG. 6: Numerical simulations of heat equation with [PITH_FULL_IMAGE:figures/full_fig_p036_6.png]

**Figure 7.** Figure 7: FIG. 7: Numerical simulation of heat equation with [PITH_FULL_IMAGE:figures/full_fig_p037_7.png]

**Figure 8.** Figure 8: FIG. 8: Numerical simulations of heat equation with [PITH_FULL_IMAGE:figures/full_fig_p038_8.png]

**Figure 9.** Figure 9: FIG. 9: Numerical experiments for wave equation with zero boundary conditions at radius [PITH_FULL_IMAGE:figures/full_fig_p039_9.png]

Discussion (0). Sign in to comment.

Forward citations

Cited by 1 Pith paper

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score. Full citation record

Circuit-Efficient Randomized Quantum Simulation of Non-Unitary Dynamics with Observable-Driven and Symmetry-Aware Designs
quant-ph 2025-09 reject novelty 5.0 of 10

A randomized compilation of LCHS for non-unitary dynamics, with an observable-driven variant and a symmetry-aware sampler, claims reduced ancilla and circuit depth at the cost of more repetitions.

Reference graph

Works this paper leans on

76 extracted references · 61 canonical work pages · cited by 1 Pith paper

[1]

On a grid, this naturally induces a projection matrix by projecting directly onto basis elements

Projections for Dirichlet conditions and value constraints Dirichlet boundary conditions in PDEs are point-wise value constraints. On a grid, this naturally induces a projection matrix by projecting directly onto basis elements. We can define P𝐷 𝑐 to be P𝐷 𝑐 = ∑︁ 𝒋∈IΓ𝐷 | 𝒋⟩⟨ 𝒋|, (97) which is equivalent to conditioning on the boundary flag from the bounda...

work page
[2]

Thus, we will continue only with bounding the error in representing a difference formula, denoted by D, using our projection method

Neumann boundary conditions Within this work, we will not worry about the approximation error with respect to the continuum limit and assume that we are provided with an ODE that follows from a sensible discretization, achieving a target discretization error𝜀 >0. Thus, we will continue only with bounding the error in representing a difference formula, den...

work page
[3]

Projectors for Neumann conditions and derivative constraints To obtain a finite-dimensional approximation to a Neumann projection, we employ a discrete finite- difference derivative D. As observed in the statement of Problem 5, a crucial ingredient for our method to work is that the projection P𝑁 𝑐 for the approximate derivative condition is orthogonal to...

work page
[4]

Input Model for Robin Boundary Conditions The ability to implement value constraints and derivative constraints allows us to extend the applicability to another class of boundary conditions. Constraints described by a superposition of Dirichlet and Neumann conditions are called Robin boundary conditions [51]; they can be expressed as 𝛼𝑣(𝑥;𝑡)+ 𝛽𝜕𝑛·𝑣(𝑥;𝑡)= ...

work page
[5]

(133) with the temperature distribution 𝑢(𝑥,𝑦 ;𝑡) and diffusion coefficient 𝐷 > 0; for the sake of simplicity, we choose isotropic diffusion here

Heat equation with Dirichlet and Neumann boundary conditions We consider the following 2D-problem for the heat equation: 𝜕𝑡𝑢(𝑥,𝑦 ;𝑡)=𝐷Δ𝑢(𝑥,𝑦 ;𝑡)+ 𝑓(𝑥,𝑦 ;𝑡) (130) 𝑢(𝑥,𝑦 ;𝑡)|𝑥,𝑦∈Γ𝐷 =𝑔(𝑥,𝑦 ;𝑡) (131) 𝜕𝑛·𝑢(𝑥,𝑦 ;𝑡)|𝑥,𝑦∈Γ𝑁 =ℎ(𝑥,𝑦 ;𝑡) (132) 𝑢(𝑥,𝑦 ;𝑡= 0)=𝑢0(𝑥,𝑦). (133) with the temperature distribution 𝑢(𝑥,𝑦 ;𝑡) and diffusion coefficient 𝐷 > 0; for the sake of sim...

work page
[6]

In order to solve the wave equation, let 𝑤(𝑥,𝑦 ;𝑡) = 𝜕𝑡𝑢(𝑥,𝑦 ;𝑡)

Wave equation with Dirichlet boundary conditions The 2-D wave equation is defined as 𝜕𝑡𝑡𝑢(𝑥,𝑦 ;𝑡)=𝑐2Δ𝑢(𝑥,𝑦 ;𝑡), (136) with a speed-of-sound parameter 𝑐2 > 0. In order to solve the wave equation, let 𝑤(𝑥,𝑦 ;𝑡) = 𝜕𝑡𝑢(𝑥,𝑦 ;𝑡). Then, Eq. (136) is equivalent to the system 𝜕𝑡v(𝑡)=𝜕𝑡 𝑢(𝑡) 𝑤(𝑡) = 0 I 𝑐2Δ 0 𝑢(𝑡) 𝑤(𝑡) = Av(𝑡), (137) with initial data v(𝑡) = h 𝑢(𝑥,𝑦...

work page
[7]

𝑁 = 25, Δ𝑡= 10−4,𝑡= 1,𝑐2= 1

work page
[8]

5 to 8 and Fig

Discussion The results to our simulations are shown in Figs. 5 to 8 and Fig. 9. For all simulations, we can observe that the errors that we witness decrease quadratically, or slightly faster, with 𝜆, i.e. ∥u∥𝑆𝑐 ∈ 𝑂(𝜆−1/2). This aligns very well with the upper bounds derived previously. The measured errors are several orders of magnitude smaller than the u...

work page

Show all 76 references

[9]

Schr¨odinger picture

Oracles HAM-T oracles The HAM-T-oracles follow the convention from [11, 12]. The truncated Dyson series is implemented across short-time integrators with lengthΔ𝑡 using𝑀𝐷 time-steps. Let 𝑞∈[ 𝑀𝐷] and|0⟩𝑎 a register of ancillas, then ⟨0|𝑎 HAM-T𝐴𝐼,𝑞|0⟩𝑎= 𝑀𝐷−1∑︁ 𝑙=0 |𝑙⟩⟨𝑙|⊗ 𝐴𝐼(𝑞Δ𝑡...
[10]

Final Result In Table Ia, we summarize the necessary strengths of the penalty through lower bounds on 𝜆. Then, as we show in the subsequent section after stating the main complexity result, this allows us to derive the smoothness parametersΛ𝐼,Ξ𝐼 as collected in Table Ib, which...
[11]

Intuitively, simulating 𝐴𝐼(𝑡) in the interaction picture with simulation parameter𝜆𝑡 means that we simulate a highly oscillatory system, as𝜆 is large

Analysis of smoothness factors Before we can do so, we look at the impact of the penalty projection. Intuitively, simulating 𝐴𝐼(𝑡) in the interaction picture with simulation parameter𝜆𝑡 means that we simulate a highly oscillatory system, as𝜆 is large. Therefore, this will have...
[12]

In our work, this ‘refinement’ is expressed by a number of grid points𝑛𝑑

Implications for the solution of discretized PDEs with boundary conditions Typically, when quantifying the cost of numerical methods, we are interested in the effect of the particular choice of discretization given a target error on the final cost ( with a ‘finer’ discretizati...

2022
[13]

D. W. Berry, High-order quantum algorithm for solving linear differential equations, Journal of Physics A: Mathe- matical and Theoretical 47, 105301 (2014)

2014
[14]

D. W. Berry, A. M. Childs, A. Ostrander, and G. Wang, Quantum algorithm for linear differential equations with exponentially improved dependence on precision, Communications in Mathematical Physics356, 1057 (2017)

2017
[15]

D. W. Berry and P. Costa, Quantum algorithm for time-dependent differential equations using dyson series, arXiv preprint arXiv:2212.03544 (2022)

2022 arXiv
[16]

Krovi, Improved quantum algorithms for linear and nonlinear differential equations, Quantum 7, 913 (2023)

H. Krovi, Improved quantum algorithms for linear and nonlinear differential equations, Quantum 7, 913 (2023)

2023
[17]

A. M. Childs and J.-P. Liu, Quantum spectral methods for differential equations, Communications in Mathematical Physics 375, 1427–1457 (2020)

2020
[18]

A. M. Childs, J.-P. Liu, and A. Ostrander, High-precision quantum algorithms for partial differential equations, Quantum 5, 574 (2021)

2021
[19]

D. Fang, L. Lin, and Y. Tong, Time-marching based quantum solvers for time-dependent linear differential equations, Quantum 7, 955 (2023)

2023
[20]

S. S. Bharadwaj and K. R. Sreenivasan, Compact quantum algorithms for time-dependent differential equations, arXiv preprint arXiv:2405.09767 (2024)

2024 arXiv
[21]

An, J.-P

D. An, J.-P. Liu, D. Wang, and Q. Zhao, A theory of quantum differential equation solvers: limitations and fast-forwarding, arXiv preprint arXiv:2211.05246 (2022)

2022 arXiv
[22]

An and K

D. An and K. Trivisa, Quantum algorithms for linear and non-linear fractional reaction-diffusion equations, preprint (2023), arXiv:2310.18900 [quant-ph]

2023
[23]

An, J.-P

D. An, J.-P. Liu, and L. Lin, Linear combination of hamiltonian simulation for non-unitary dynamics with optimal state preparation cost, Physical Review Letters131, 150603 (2023)

2023
[24]

D. An, A. M. Childs, and L. Lin, Quantum algorithm for linear non-unitary dynamics with near-optimal dependence on all parameters, arXiv preprint arXiv:2312.03916 (2023)

2023
[25]

S. Jin, N. Liu, and Y. Yu, Quantum simulation of partial differential equations via schr ¨odingerization, Physical Review Letters 133, 230602 (2024)

2024
[26]

Lu, H.-E

R. Lu, H.-E. Li, Z. Liu, and J.-P. Liu, Infinite-dimensional extension of the linear combination of hamiltonian simulation: Theorems and applications, arXiv preprint arXiv:2502.19688 (2025)

2025 arXiv
[27]

Shang, N

Z.-X. Shang, N. Guo, D. An, and Q. Zhao, Design nearly optimal quantum algorithm for linear differential equations via lindbladians, arXiv preprint arXiv:2410.19628 (2024)

2024
[28]

Y. Cao, A. Papageorgiou, I. Petras, J. Traub, and S. Kais, Quantum algorithm and circuit design solving the poisson equation, New Journal of Physics 15, 013021 (2013)

2013
[29]

Bagherimehrab, K

M. Bagherimehrab, K. Nakaji, N. Wiebe, and A. Aspuru-Guzik, Fast quantum algorithm for differential equations, arXiv preprint arXiv:2306.11802 (2023)

2023
[30]

Novikau and I

I. Novikau and I. Joseph, Quantum algorithm for the advection-diffusion equation and the koopman-von neumann approach to nonlinear dynamical systems, Computer Physics Communications 309, 109498 (2025)

2025
[31]

P. C. Costa, S. Jordan, and A. Ostrander, Quantum algorithm for simulating the wave equation, Physical Review A 99, 012323 (2019)

2019
[32]

S. Jin, N. Liu, and C. Ma, Quantum simulation of maxwell’s equations via schr¨odingerisation, ESAIM: Mathematical Modelling and Numerical Analysis 58, 1853 (2024)

2024
[33]

J.-P. Liu, H. Ø. Kolden, H. K. Krovi, N. F. Loureiro, K. Trivisa, and A. M. Childs, Efficient quantum algorithm for dissipative nonlinear differential equations, Proceedings of the National Academy of Sciences 118, e2026805118 (2021). 51

2021
[34]

J.-P. Liu, D. An, D. Fang, J. Wang, G. H. Low, and S. Jordan, Efficient quantum algorithm for nonlinear reaction– diffusion equations and energy estimation, Communications in Mathematical Physics 404, 963 (2023)

2023
[35]

Costa, P

P. Costa, P. Schleich, M. E. Morales, and D. W. Berry, Further improving quantum algorithms for nonlinear differential equations via higher-order methods and rescaling, arXiv preprint arXiv:2312.09518 (2023)

2023 arXiv
[36]

X. Li, X. Yin, N. Wiebe, J. Chun, G. K. Schenter, M. S. Cheung, and J. M¨ ulmenst¨adt, Potential quantum advantage for simulation of fluid dynamics, Phys. Rev. Res.7, 013036 (2025)

2025
[37]

S. S. Bharadwaj and K. R. Sreenivasan, Towards simulating fluid flows with quantum computing, S ¯adhan¯a 50, 1 (2025)

2025
[38]

Penuel, A

J. Penuel, A. Katabarwa, P. D. Johnson, C. Farquhar, Y. Cao, and M. C. Garrett, Feasibility of accelerating incompressible computational fluid dynamics simulations with fault-tolerant quantum computers, arXiv preprint arXiv:2406.06323 (2024)

2024
[39]

Novikau, I

I. Novikau, I. Dodin, and E. Startsev, Simulation of linear non-hermitian boundary-value problems with quantum singular-value transformation, Physical Review Applied19 (2023)

2023
[40]

S. Jin, X. Li, N. Liu, and Y. Yu, Quantum simulation for partial differential equations with physical boundary or interface conditions, Journal of Computational Physics 498, 112707 (2024)

2024
[41]

S. Jin, X. Li, N. Liu, and Y. Yu, Quantum simulation for quantum dynamics with artificial boundary conditions, SIAM Journal on Scientific Computing 46, B403 (2024)

2024
[42]

Kharazi, A

T. Kharazi, A. M. Alkadri, J.-P. Liu, K. K. Mandadapu, and K. B. Whaley, Explicit block encodings of boundary value problems for many-body elliptic operators, Quantum 9, 1764 (2025)

2025
[43]

Mangin-Brinet, J

M. Mangin-Brinet, J. Zhang, D. Lacroix, and E. A. R. Guzman, Efficient solution of the non-unitary time-dependent schrodinger equation on a quantum computer with complex absorbing potential, Quantum 8, 1311 (2024)

2024
[44]

D. G. Luenberger, Y. Ye, D. G. Luenberger, and Y. Ye, Penalty and barrier methods, Linear and Nonlinear Programming , 397 (2016)

2016
[45]

J. Muga, J. Palao, B. Navarro, and I. Egusquiza, Complex absorbing potentials, Physics Reports 395, 357 (2004)

2004
[46]

Mangin-Brinet, J

M. Mangin-Brinet, J. Zhang, D. Lacroix, and E. A. Ruiz Guzman, Efficient solution of the non-unitary time- dependent Schrodinger equation on a quantum computer with complex absorbing potential, Quantum 8, 1311 (2024)

2024
[47]

S. V. Tsynkov, Artificial boundary conditions for the numerical simulation of unsteady acoustic waves, Journal of Computational Physics 189, 626 (2003)

2003
[48]

Burgarth, P

D. Burgarth, P. Facchi, H. Nakazato, S. Pascazio, and K. Yuasa, Generalized adiabatic theorem and strong-coupling limits, Quantum 3, 152 (2019)

2019
[49]

Burgarth, P

D. Burgarth, P. Facchi, G. Gramegna, and K. Yuasa, One bound to rule them all: from adiabatic to zeno, Quantum 6, 737 (2022)

2022
[50]

Van der Sluis, Perturbations of eigenvalues of non-normal matrices, Communications of the ACM18, 30 (1975)

A. Van der Sluis, Perturbations of eigenvalues of non-normal matrices, Communications of the ACM18, 30 (1975)

1975
[51]

A. B. Catli, S. Simon, and N. Wiebe, Exponentially better bounds for quantum optimization via dynamical simula- tion, arXiv preprint arXiv:2502.04285 (2025)

2025
[52]

Kubo, Statistical-mechanical theory of irreversible processes

R. Kubo, Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems, Journal of the physical society of Japan 12, 570 (1957)

1957
[53]

Sticlet, B

D. Sticlet, B. D ´ora, and C. P. Moca, Kubo formula for non-hermitian systems and tachyon optical conductivity, Physical Review Letters128, 016802 (2022)

2022
[54]

K. T. Geier and P. Hauke, From non-hermitian linear response to dynamical correlations and fluctuation-dissipation relations in quantum many-body systems, PRX Quantum 3, 030308 (2022)

2022
[55]

L. Pan, X. Chen, Y. Chen, and H. Zhai, Non-hermitian linear response theory, Nature Physics 16, 767 (2020)

2020
[56]

Ames, Numerical Methods for Partial Differential Equations , Computer Science and Scientific Computing (Elsevier Science, 2014)

W. Ames, Numerical Methods for Partial Differential Equations , Computer Science and Scientific Computing (Elsevier Science, 2014)

2014
[57]

Evans, J

G. Evans, J. Blackledge, and P. Yardley, Numerical Methods for Partial Differential Equations , Springer Under- graduate Mathematics Series (Springer London, 2012)

2012
[58]

G. H. Low and N. Wiebe, Hamiltonian simulation in the interaction picture, arXiv preprint arXiv:1805.00675 (2018)

2018 arXiv
[59]

Atia and D

Y. Atia and D. Aharonov, Fast-forwarding of hamiltonians and exponentially precise measurements, Nature com- munications 8, 1572 (2017)

2017
[60]

Gidney, Halving the cost of quantum addition, Quantum 2, 74 (2018)

C. Gidney, Halving the cost of quantum addition, Quantum 2, 74 (2018)

2018
[61]

S. Gu, R. D. Somma, and B. S ¸ahino˘glu, Fast-forwarding quantum evolution, Quantum5, 577 (2021). 52

2021
[62]

Cen and Q

S. Cen and Q. Zhang, Multiphysics Modeling: Numerical Methods and Engineering Applications(Academic Press, 2015)

2015
[63]

Gustafson and T

K. Gustafson and T. Abe, The third boundary condition–was it robin’s?, Mathematical Intelligencer 20 (1998)

1998
[64]

Linden, A

N. Linden, A. Montanaro, and C. Shao, Quantum vs. classical algorithms for solving the heat equation, Communi- cations in Mathematical Physics 395, 601 (2022)

2022
[65]

Babbush, D

R. Babbush, D. W. Berry, R. Kothari, R. D. Somma, and N. Wiebe, Exponential quantum speedup in simulating coupled classical oscillators, Physical Review X13, 041041 (2023)

2023
[66]

Bogacki and L

P. Bogacki and L. F. Shampine, A 3 (2) pair of runge-kutta formulas, Applied Mathematics Letters 2, 321 (1989)

1989
[67]

S. Jin, N. Liu, C. Ma, and Y. Yu, On the schr\” odingerization method for linear non-unitary dynamics with optimal dependence on matrix queries, arXiv preprint arXiv:2505.00370 (2025)

2025
[68]

Gily ´en, Y

A. Gily ´en, Y. Su, G. H. Low, and N. Wiebe, Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics, in Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing (2019) pp. 193–204

2019
[69]

J. D. Lambert, Numerical methods for ordinary differential systems: the initial value problem(John Wiley & Sons, Inc., 1991)

1991
[70]

I. D. Kivlichan, N. Wiebe, R. Babbush, and A. Aspuru-Guzik, Bounding the costs of quantum simulation of many- body physics in real space, Journal of Physics A: Mathematical and Theoretical 50, 305301 (2017), publisher: {IOP} Publishing

2017
[71]

Riss and H.-D

U. Riss and H.-D. Meyer, Calculation of resonance energies and widths using the complex absorbing potential method, Journal of Physics B: Atomic, Molecular and Optical Physics 26, 4503 (1993)

1993
[72]

Baer, findiff software package (2018), https://github.com/maroba/findiff

M. Baer, findiff software package (2018), https://github.com/maroba/findiff

2018
[73]

Dattoli, A

G. Dattoli, A. Torre, and R. Mignani, Non-hermitian evolution of two-level quantum systems, Physical Review A 42, 1467 (1990)

1990
[74]

D. C. Brody, Biorthogonal quantum mechanics, Journal of Physics A: Mathematical and Theoretical 47, 035305 (2013)

2013
[75]

A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Theory of trotter error with commutator scaling, Physical Review X 11, 011020 (2021)

2021
[76]

G. H. Hardy, J. E. Littlewood, and G. P´olya, Inequalities (Cambridge mathematical library) (Cambridge university press, 1934). 53 Appendix A: Proof of Lemma 15 Lemma (Parametrized exponential of an orthogonal projection; Lemma 15 in main text) . Let P, Q be orthogonal project...

1934

Pith tools

Reviewed August 6, 2026 · model on record in the stance chip above.

[1] [1]

On a grid, this naturally induces a projection matrix by projecting directly onto basis elements

Projections for Dirichlet conditions and value constraints Dirichlet boundary conditions in PDEs are point-wise value constraints. On a grid, this naturally induces a projection matrix by projecting directly onto basis elements. We can define P𝐷 𝑐 to be P𝐷 𝑐 = ∑︁ 𝒋∈IΓ𝐷 | 𝒋⟩⟨ 𝒋|, (97) which is equivalent to conditioning on the boundary flag from the bounda...

work page

[2] [2]

Thus, we will continue only with bounding the error in representing a difference formula, denoted by D, using our projection method

Neumann boundary conditions Within this work, we will not worry about the approximation error with respect to the continuum limit and assume that we are provided with an ODE that follows from a sensible discretization, achieving a target discretization error𝜀 >0. Thus, we will continue only with bounding the error in representing a difference formula, den...

work page

[3] [3]

Projectors for Neumann conditions and derivative constraints To obtain a finite-dimensional approximation to a Neumann projection, we employ a discrete finite- difference derivative D. As observed in the statement of Problem 5, a crucial ingredient for our method to work is that the projection P𝑁 𝑐 for the approximate derivative condition is orthogonal to...

work page

[4] [4]

Input Model for Robin Boundary Conditions The ability to implement value constraints and derivative constraints allows us to extend the applicability to another class of boundary conditions. Constraints described by a superposition of Dirichlet and Neumann conditions are called Robin boundary conditions [51]; they can be expressed as 𝛼𝑣(𝑥;𝑡)+ 𝛽𝜕𝑛·𝑣(𝑥;𝑡)= ...

work page

[5] [5]

(133) with the temperature distribution 𝑢(𝑥,𝑦 ;𝑡) and diffusion coefficient 𝐷 > 0; for the sake of simplicity, we choose isotropic diffusion here

Heat equation with Dirichlet and Neumann boundary conditions We consider the following 2D-problem for the heat equation: 𝜕𝑡𝑢(𝑥,𝑦 ;𝑡)=𝐷Δ𝑢(𝑥,𝑦 ;𝑡)+ 𝑓(𝑥,𝑦 ;𝑡) (130) 𝑢(𝑥,𝑦 ;𝑡)|𝑥,𝑦∈Γ𝐷 =𝑔(𝑥,𝑦 ;𝑡) (131) 𝜕𝑛·𝑢(𝑥,𝑦 ;𝑡)|𝑥,𝑦∈Γ𝑁 =ℎ(𝑥,𝑦 ;𝑡) (132) 𝑢(𝑥,𝑦 ;𝑡= 0)=𝑢0(𝑥,𝑦). (133) with the temperature distribution 𝑢(𝑥,𝑦 ;𝑡) and diffusion coefficient 𝐷 > 0; for the sake of sim...

work page

[6] [6]

In order to solve the wave equation, let 𝑤(𝑥,𝑦 ;𝑡) = 𝜕𝑡𝑢(𝑥,𝑦 ;𝑡)

Wave equation with Dirichlet boundary conditions The 2-D wave equation is defined as 𝜕𝑡𝑡𝑢(𝑥,𝑦 ;𝑡)=𝑐2Δ𝑢(𝑥,𝑦 ;𝑡), (136) with a speed-of-sound parameter 𝑐2 > 0. In order to solve the wave equation, let 𝑤(𝑥,𝑦 ;𝑡) = 𝜕𝑡𝑢(𝑥,𝑦 ;𝑡). Then, Eq. (136) is equivalent to the system 𝜕𝑡v(𝑡)=𝜕𝑡 𝑢(𝑡) 𝑤(𝑡) = 0 I 𝑐2Δ 0 𝑢(𝑡) 𝑤(𝑡) = Av(𝑡), (137) with initial data v(𝑡) = h 𝑢(𝑥,𝑦...

work page

[7] [7]

𝑁 = 25, Δ𝑡= 10−4,𝑡= 1,𝑐2= 1

work page

[8] [8]

5 to 8 and Fig

Discussion The results to our simulations are shown in Figs. 5 to 8 and Fig. 9. For all simulations, we can observe that the errors that we witness decrease quadratically, or slightly faster, with 𝜆, i.e. ∥u∥𝑆𝑐 ∈ 𝑂(𝜆−1/2). This aligns very well with the upper bounds derived previously. The measured errors are several orders of magnitude smaller than the u...

work page