LCHS: Linear Combination of Hamiltonian Simulation

• LCHS is a quantum framework that decomposes non-unitary dynamics into an integral of unitary Hamiltonian simulations.
• It employs kernel truncation and Gaussian quadrature to achieve near-optimal query scaling and resource efficiency in simulating dissipative systems.
• The method leverages block-encoding, amplitude amplification, and advanced techniques to extend its applicability to PDEs and eigenvalue transformations.

Linear Combination of Hamiltonian Simulation (LCHS) is a quantum-algorithmic framework for simulating linear non-unitary dynamics by rewriting the target propagator as a linear combination of unitary evolutions, each of which is a Hamiltonian simulation problem. In its standard form, LCHS addresses linear ordinary differential equations of the form $\partial_t u(t)=-A(t)u(t)$ or $\partial_t u(t)=-A(t)u(t)+b(t)$, with the generator decomposed as $A(t)=L(t)+iH(t)$, where $L(t)$ and $H(t)$ are Hermitian and $L(t)\succeq 0$ after a possible identity shift. The central idea is that the non-unitary propagator $\mathcal T e^{-\int_0^t A(s)\,ds}$ can be expressed as an integral over unitaries of the form $\mathcal T e^{-i\int_0^t(kL(s)+H(s))\,ds}$, after which truncation, quadrature, block-encoding, linear-combination-of-unitaries, and amplitude amplification yield a quantum implementation (An et al., 2023, An et al., 2023).

1. Foundational identity and problem setting

The modern LCHS formalism emerged from a theorem showing that a general dissipative generator can be represented as a continuous linear combination of Hamiltonian evolutions. For the time-dependent homogeneous ODE

$\partial_t u(t)=-A(t)u(t),\qquad A(t)=L(t)+iH(t),$

with $L(t)\succeq 0$, one has

$\partial_t u(t)=-A(t)u(t)+b(t)$0

for a suitable kernel $\partial_t u(t)=-A(t)u(t)+b(t)$1 satisfying normalization and decay conditions (An et al., 2023). The same framework extends to the inhomogeneous case via Duhamel’s principle, and in the time-independent case it reduces to a direct integral formula for $\partial_t u(t)=-A(t)u(t)+b(t)$2 (An et al., 2023).

A widely used one-parameter family of kernels is

$\partial_t u(t)=-A(t)u(t)+b(t)$3

which leads to the identity

$\partial_t u(t)=-A(t)u(t)+b(t)$4

for the propagator $\partial_t u(t)=-A(t)u(t)+b(t)$5 (An et al., 2023). In the time-independent specialization summarized in later constant-factor analyses,

$\partial_t u(t)=-A(t)u(t)+b(t)$6

with $\partial_t u(t)=-A(t)u(t)+b(t)$7 (Pocrnic et al., 25 Jun 2025).

This identity places LCHS conceptually between dissipative dynamics and standard Hamiltonian simulation. A plausible implication is that the method transfers advances in unitary simulation—such as qubitization, generalized quantum signal processing, and improved block-encoding analyses—into the non-unitary setting with comparatively little change in the outer algorithmic structure.

2. Kernel design, truncation, and quadrature

An LCHS algorithm becomes finite only after truncating the infinite $\partial_t u(t)=-A(t)u(t)+b(t)$8-integral and discretizing it. Early formulations used a cutoff $\partial_t u(t)=-A(t)u(t)+b(t)$9 and trapezoidal discretization for error $A(t)=L(t)+iH(t)$0 (An et al., 2023). The subsequent near-optimal formulation replaced this with a sharper kernel family and composite Gaussian quadrature, obtaining

$A(t)=L(t)+iH(t)$1

and hence

$A(t)=L(t)+iH(t)$2

LCU terms after discretization (An et al., 2023). This is the source of the widely cited logarithmic-to-sub-polynomial error dependence of LCHS.

A later constant-factor study derived explicit non-asymptotic formulas for both truncation and Gaussian quadrature. For truncation error $A(t)=L(t)+iH(t)$3, the cutoff may be chosen as

$A(t)=L(t)+iH(t)$4

where $A(t)=L(t)+iH(t)$5 is the principal Lambert $A(t)=L(t)+iH(t)$6-function. For discretization error $A(t)=L(t)+iH(t)$7 on $A(t)=L(t)+iH(t)$8, the order-$A(t)=L(t)+iH(t)$9 Gauss–Legendre rule can be chosen as

$L(t)$0

with total LCU size

$L(t)$1

(Pocrnic et al., 25 Jun 2025). These expressions recover the asymptotic estimates but, more importantly for resource estimation, fix the constants that dominate realistic fault-tolerant costs.

Kernel construction itself has also diversified. A Fourier-transform-based reformulation shows that one may choose any $L(t)$2 satisfying $L(t)$3 for $L(t)$4, and then define

$L(t)$5

This removes the technical requirements beyond the real axis present in the original half-plane formulation and yields a new kernel family. For truncation error $L(t)$6, that formalism reports a $L(t)$7 times reduction in the quantum differential equation algorithms based on LCHS and an $L(t)$8 times reduction in its quantum circuit depth (Huang et al., 27 Aug 2025).

3. Coherent implementation: block-encoding, SELECT, and amplification

The coherent realization of LCHS follows the linear-combination-of-unitaries pattern. In the time-dependent oracle model, one assumes access to HAM-T oracles that block-encode time samples of $L(t)$9 and $H(t)$0 on short time slices, prepares an ancilla superposition over quadrature nodes, and implements

$H(t)$1

after which

$H(t)$2

block-encodes the discrete approximation to the non-unitary propagator (An et al., 2023).

For the time-independent block-encoding model, later work makes the SELECT step substantially more explicit. If $H(t)$3 is an $H(t)$4 block encoding of $H(t)$5, define the effective Hamiltonian

$H(t)$6

This $H(t)$7 can be $H(t)$8-encoded using just one call to $H(t)$9 and one to $L(t)\succeq 0$0 plus $L(t)\succeq 0$1 controlled-$L(t)\succeq 0$2 gates. Qubitization or generalized quantum signal processing then implements

$L(t)\succeq 0$3

as a $L(t)\succeq 0$4 block encoding using

$L(t)\succeq 0$5

calls to controlled-$L(t)\succeq 0$6 or its inverse, with $L(t)\succeq 0$7 (Pocrnic et al., 25 Jun 2025).

The outer LCU circuit is then

$L(t)\succeq 0$8

which realizes the approximate propagator $L(t)\succeq 0$9 as a block encoding. Applying $\mathcal T e^{-\int_0^t A(s)\,ds}$0 to an initial-state block encoding $\mathcal T e^{-\int_0^t A(s)\,ds}$1 yields a block encoding of $\mathcal T e^{-\int_0^t A(s)\,ds}$2 with

$\mathcal T e^{-\int_0^t A(s)\,ds}$3

(Pocrnic et al., 25 Jun 2025).

Because postselection amplitude depends on the norm of the evolved state, fixed-point oblivious amplitude amplification is used to boost success probability to near unity. If $\mathcal T e^{-\int_0^t A(s)\,ds}$4 is a known lower bound on the overlap, the amplification cost is

$\mathcal T e^{-\int_0^t A(s)\,ds}$5

with the explicit expression for $\mathcal T e^{-\int_0^t A(s)\,ds}$6 given in the constant-factor analysis, and the total number of calls to $\mathcal T e^{-\int_0^t A(s)\,ds}$7 becomes

$\mathcal T e^{-\int_0^t A(s)\,ds}$8

subject to the sub-error budget

$\mathcal T e^{-\int_0^t A(s)\,ds}$9

(Pocrnic et al., 25 Jun 2025).

4. Complexity, optimality, and constant factors

A defining property of LCHS is its near-optimal parameter dependence. In the 2023 near-optimal analysis, letting $\mathcal T e^{-i\int_0^t(kL(s)+H(s))\,ds}$0, the total number of Hamiltonian-simulation queries is

$\mathcal T e^{-i\int_0^t(kL(s)+H(s))\,ds}$1

while the number of initial-state-oracle queries is $\mathcal T e^{-i\int_0^t(kL(s)+H(s))\,ds}$2, which matches the $\mathcal T e^{-i\int_0^t(kL(s)+H(s))\,ds}$3 lower bound for state preparation cost (An et al., 2023). The earlier formulation stated the same optimality more explicitly: any quantum ODE solver must use $\mathcal T e^{-i\int_0^t(kL(s)+H(s))\,ds}$4 queries to the initial-state oracle, whereas LCHS uses $\mathcal T e^{-i\int_0^t(kL(s)+H(s))\,ds}$5 (An et al., 2023).

In the time-independent special case, replacing truncated Dyson simulation by QSP or QSVT yields overall matrix-query complexity

$\mathcal T e^{-i\int_0^t(kL(s)+H(s))\,ds}$6

which was described as fully recovering the “fast-forwarded” performance asymptotically up to the $\mathcal T e^{-i\int_0^t(kL(s)+H(s))\,ds}$7 factor (An et al., 2023). Later work sharpened the error dependence further. For bounded, time-independent $\mathcal T e^{-i\int_0^t(kL(s)+H(s))\,ds}$8, a generalized LCHS construction gives a block-encoding of $\mathcal T e^{-i\int_0^t(kL(s)+H(s))\,ds}$9 with

$\partial_t u(t)=-A(t)u(t),\qquad A(t)=L(t)+iH(t),$0

queries to the block-encoding oracle for $\partial_t u(t)=-A(t)u(t),\qquad A(t)=L(t)+iH(t),$1, and the normalized evolved state can be prepared with

$\partial_t u(t)=-A(t)u(t),\qquad A(t)=L(t)+iH(t),$2

queries to the initial-state oracle; both complexities are stated to be optimal in all parameters (Low et al., 26 Aug 2025). The same work also states that any improvement exceeding a constant factor of approximately $\partial_t u(t)=-A(t)u(t),\qquad A(t)=L(t)+iH(t),$3 is infeasible (Low et al., 26 Aug 2025).

Constant factors, however, are decisive in practice. The detailed 2025 resource analysis identifies several improvements: tight truncation via a Lambert $\partial_t u(t)=-A(t)u(t),\qquad A(t)=L(t)+iH(t),$4 bound, tighter Gaussian quadrature bounds, a diagonal block-encoding of $\partial_t u(t)=-A(t)u(t),\qquad A(t)=L(t)+iH(t),$5 that avoids $\partial_t u(t)=-A(t)u(t),\qquad A(t)=L(t)+iH(t),$6 separate controlled calls to $\partial_t u(t)=-A(t)u(t),\qquad A(t)=L(t)+iH(t),$7, generalized quantum signal processing “doubling,” and constant-factor robust-block-encoding analyses for qubitization and fixed-point amplitude amplification (Pocrnic et al., 25 Jun 2025). Under the specialization

$\partial_t u(t)=-A(t)u(t),\qquad A(t)=L(t)+iH(t),$8

the resulting query counts are reported as two orders of magnitude smaller than those of the “randomized linear solver” method, with approximately $\partial_t u(t)=-A(t)u(t),\qquad A(t)=L(t)+iH(t),$9 fewer block-encoding calls under identical assumptions (Pocrnic et al., 25 Jun 2025). That same analysis summarizes the resulting behavior as linear scaling in $L(t)\succeq 0$0, logarithmic-to-sub-polynomial scaling in $L(t)\succeq 0$1, and markedly reduced block-encoding query complexity in the “no-fast-forwarding” regime (Pocrnic et al., 25 Jun 2025).

5. Variants, generalizations, and adjacent formalisms

LCHS has rapidly expanded beyond its original role as an ODE solver. One direction generalizes from matrix exponentials to broader eigenvalue transformations. A Laplace-transform-based extension represents

$L(t)\succeq 0$2

and then substitutes the LCHS representation for $L(t)\succeq 0$3, enabling transformations such as $L(t)\succeq 0$4 and $L(t)\succeq 0$5 without explicitly inverting $L(t)\succeq 0$6 (An et al., 2024). A later Weyl-calculus-based theory goes further, constructing LCHS formulas for general matrix functions $L(t)\succeq 0$7 analytic on the numerical range of $L(t)\succeq 0$8, including non-normal $L(t)\succeq 0$9, and reports optimal $\partial_t u(t)=-A(t)u(t)+b(t)$00 query complexity scaling for quantum eigenvalue transformation (Ni et al., 29 Jun 2026).

A second direction extends the domain of operators. Inf-LCHS generalizes the finite-dimensional theorem to infinite-dimensional Hilbert spaces and possibly unbounded operators, under domain, kernel, and common-domain assumptions. It introduces Inf-LCHS-Gaussian and Inf-LCHS-MC as quadrature-based and Monte-Carlo-based discretizations, respectively, and applies the method to linear parabolic PDEs, birth–death processes, Schrödinger equations with complex potentials, Lindblad equations, and black hole thermal field equations (Lu et al., 27 Feb 2025).

A third direction concerns implementation architecture. Hybrid oscillator–qubit LCHS replaces the discrete quadrature register by a continuous-variable ancilla mode, using

$\partial_t u(t)=-A(t)u(t)+b(t)$01

to encode the kernel directly. This removes the explicit $\partial_t u(t)=-A(t)u(t)+b(t)$02 ancilla-qubit overhead of qubit-only LCHS. In heat-equation benchmarks, the Law–Eberly protocol achieves end-to-end solution fidelity at least $\partial_t u(t)=-A(t)u(t)+b(t)$03, and a comparison with a matrix-product-state-based discrete-variable implementation reports a drastically reduced CNOT count, stated as $\partial_t u(t)=-A(t)u(t)+b(t)$04 versus $\partial_t u(t)=-A(t)u(t)+b(t)$05 in a Dirichlet benchmark (Das et al., 11 May 2026).

A fourth direction targets circuit structure rather than asymptotics. Random-LCHS randomizes both the outer LCU layer and the inner Hamiltonian-simulation layer, with observable-driven and symmetry-aware variants for expectation estimation and physically symmetric models (Yang et al., 9 Sep 2025). Another recent refinement shows that implementing the Hamiltonian-simulation subroutines with Multi-Product Formulas yields commutator-sensitive complexity bounds, and that the chosen quadrature rule affects not only discretization error but also commutator structure and query complexity (Aftab et al., 9 Jun 2026).

6. Applications, implementations, and interpretive cautions

LCHS has been applied most visibly to linear PDEs. One line of work maps second-order PDEs with spatially varying parameters to first-order ODE systems by spatial discretization, then uses logic minimization and matrix-product-state compression to reduce the number of Pauli strings and compile the coefficient oracle. In a room-acoustics example on a $\partial_t u(t)=-A(t)u(t)+b(t)$06 grid, a naive $\partial_t u(t)=-A(t)u(t)+b(t)$07 terms for one coefficient collapsed to just $\partial_t u(t)=-A(t)u(t)+b(t)$08 terms, and an MPS of bond $\partial_t u(t)=-A(t)u(t)+b(t)$09 for the LCU coefficient state achieved $\partial_t u(t)=-A(t)u(t)+b(t)$10 fidelity for $\partial_t u(t)=-A(t)u(t)+b(t)$11 quadrature points (Sato et al., 2024). Related work gives explicit near-optimal circuits for the advection-diffusion equation, with $\partial_t u(t)=-A(t)u(t)+b(t)$12, $\partial_t u(t)=-A(t)u(t)+b(t)$13, total qubits about $\partial_t u(t)=-A(t)u(t)+b(t)$14, and success probability $\partial_t u(t)=-A(t)u(t)+b(t)$15 after amplitude amplification in one representative run (Novikau et al., 19 Jan 2025). Subsequent papers provide boundary-aware and structure-preserving PDE frameworks, including mixed Dirichlet, Neumann, Robin, and periodic boundary conditions, explicit block-encodings, and end-to-end complexity balances (Liang et al., 2 Jun 2026, Chen et al., 17 May 2026).

LCHS has also been explored in hardware-oriented settings. A variational algorithm on the superconducting processor Wukong uses LCHS only at the loss-function level, replacing coherent long-depth LCU by a fixed-depth parameterized ansatz and simplified Hadamard tests; the reported circuit depth is independent of total simulation time, and experiments on dissipative transverse Ising and interacting Hatano–Nelson models tracked theoretical curves accurately (Liu et al., 23 Oct 2025). At the fault-tolerant end, an end-to-end algorithm for rapidly distorted turbulence formulates the non-unitary generator through LCHS and estimates that a minimal non-trivial instance on a $\partial_t u(t)=-A(t)u(t)+b(t)$16 grid can be encoded in about $\partial_t u(t)=-A(t)u(t)+b(t)$17 qubits with circuit depth approximately $\partial_t u(t)=-A(t)u(t)+b(t)$18 controlled-unitary layers (Meng et al., 24 Nov 2025).

Two interpretive cautions recur in the literature. First, LCHS should not be conflated with the broader linear-combination-of-unitaries literature for unitary Hamiltonian simulation. Earlier work on linear combinations of unitary operations and multi-product formulas established the ancilla-based block-encoding and amplitude-amplification toolkit on which LCHS relies (Childs et al., 2012), and photonic experiments demonstrated such linear-combination methods for unitary dynamics (Yu et al., 2022); the later LCHS literature is distinguished by its specific integral representation of non-unitary evolution in terms of Hamiltonian simulations (An et al., 2023). Second, the condition $\partial_t u(t)=-A(t)u(t)+b(t)$19 is not a cosmetic assumption. It encodes dissipativity or non-positive logarithmic norm, guarantees $\partial_t u(t)=-A(t)u(t)+b(t)$20 in the standard setting, and underlies both the correctness of the kernel identity and the controlled success probability of postselected implementations (An et al., 2023).

Within quantum algorithms for non-unitary dynamics, LCHS has therefore become both a specific algorithmic family and a broader design pattern: represent the target contraction as a superposition of unitary evolutions, discretize that superposition carefully, and spend algorithmic effort on the constants hidden inside truncation, quadrature, block-encoding, and amplification. The progression from the original identity to constant-factor analyses, optimal-scaling variants, continuous-variable encodings, and PDE-specific compilers suggests that the central research questions have shifted from mere asymptotic existence to resource-tight realizations and domain-specific structure exploitation.