Quantum simulation of real-world nonlinear dynamics via Koopman method

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Quantum simulation of real-world nonlinear dynamics via Koopman method

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

artificialintelligence knowledgerepresentation reasoning

artificial intelligence, reasoning, knowledge representation

Quantum simulation of real-world nonlinear dynamics via Koopman method Baoyang Zhang,1 Dong An,2 Zhaoyuan Meng,3 Yefei Yu,4 Xiaoxiao Xiao,4 Zhen Lu,1, ∗ and Yue Yang1, 5, †

arXiv:2607.07338v1 [quant-ph] 8 Jul 2026

1 State Key Laboratory for Turbulence and Complex Systems, School of Mechanics and Engineering Science, Peking University, Beijing 100871, China 2 Beijing International Center for Mathematical Research, Peking University, Beijing 100871, China 3 Institute of Mechanics, State Key Laboratory of Nonlinear Mechanics, Chinese Academy of Sciences, Beijing 100190, China 4 Beijing Academy of Quantum Information Sciences, Beijing 100193, China 5 HEDPS-CAPT, Peking University, Beijing 100871, China (Dated: July 9, 2026)

Nonlinear dynamics is ubiquitous in nature, ranging from chemical pattern formation to ocean circulation, yet its simulation on quantum computers is fundamentally limited by the unitary nature of quantum evolution. We propose the quantum Koopman method, a data-driven framework that embeds nonlinear dynamics into a learned linear representation and implements the resulting evolution using shallow quantum circuits. This method learns Koopman observables from trajectory data, projects the lifted dynamics onto a finite-dimensional subspace, and decomposes the corresponding non-unitary propagator into parallel spectral channels. We utilize the Koopman method on a superconducting processor to simulate three distinct nonlinear systems, comprising reaction-diffusion dynamics, fluid motion on a sphere, and satellite-derived observations of Gulf Stream currents, employing up to 32 parallel circuits of 10 qubits. These quantum simulations capture the dominant multiscale patterns and statistical signatures of the underlying dynamics, and reveal a transition from performance limited by hardware noise in weakly nonlinear systems to performance limited by finite-dimensional Koopman representations as nonlinear scale interactions increase. This transition identifies a practical boundary for quantum-amenable nonlinear dynamics, establishing a hardware-validated route for simulating moderately nonlinear dynamics on near-term quantum hardware.

Quantum computing promises exponential speedups for the scientific simulation of complex systems [1, 2]. However, quantum state evolution is governed by unitary operators, rendering these dynamics strictly linear and reversible [3]. This intrinsic linearity stands in sharp contrast to the nonlinear dynamics ubiquitous in real-world phenomena, ranging from planetary-scale ocean circulation [4] to microscale kinetics in engine combustion [5]. Bridging this divide is essential for extending quantum utility to nonlinear regimes [6–8], yet any viable framework must reconcile mathematical rigor with the resource constraints of current quantum hardware [9, 10]. Efficiently embedding nonlinear, dissipative dynamics within the unitary operations of quantum processors therefore remains a fundamental open challenge for the quantum simulation of real-world physical problems. Existing quantum approaches to nonlinear dynamics generally fall into two broad families. The first follows a three-stage pipeline consisting of analytical linearization [11–15], mapping onto a quantum linear-system or Hamiltonian-simulation algorithms [16–21], and compilation into executable circuits [22], with several alternative methods bypassing explicit linearization [23–26]. However, this family remains constrained by two limitations. First, the three stages are typically developed in isolation, yielding circuits whose depth exceeds what current noisy intermediate-scale quantum (NISQ) [27] devices can reliably execute [10, 28]. Second, the underlying analytical linearizations typically converge only under weak nonlinearity, leaving moderately nonlinear dynamics of realworld systems out of reach. The second family comprises vari-

∗ [email protected] † [email protected]

ational quantum algorithms [29–33], which treat the quantum circuits as trainable, high-dimensional regressors. Without a structural connection to the underlying physics, these problemagnostic circuit ansatzes lack the inductive bias needed to navigate complex parameter landscapes and frequently encounter barren plateaus [34] and local minima [35]. Consequently, hardware demonstrations of these methods have been limited to 1D and 2D benchmarks using no more than 10 qubits [36– 40]. Although a recent Koopman-inspired approach [41] has explored embedding physical priors into variational circuits, a unified framework that can simultaneously access moderately nonlinear regimes and compile into hardware-feasible circuits remains absent. Here, we introduce the quantum Koopman method (QKM), which provides both a theoretical foundation and a hardwareefficient implementation for simulating nonlinear dynamics on an 𝑛-qubit quantum processor, as illustrated in Fig. 1. Specifically, the Koopman representation lifts the nonlinear dynamics (Fig. 1a) into an infinite-dimensional linear space (Fig. 1c), which is subsequently projected onto a finite, 2𝑛 -dimensional subspace spanned by learned observables (Figs. 1b,d). A suite of three theorems connects these non-unitary linearized dynamics to hardware-executable quantum circuits, yielding a topology-native ansatz derived from the linear combination of Hamiltonian simulation (LCHS) [18] rather than heuristic design. To prepare quantum states, a classical neural network (NN) encoder maps physical initial conditions to circuit parameters, achieving a poly(𝑛)-depth encoding by leveraging the regularity of physical fields (Fig. 1d). The Koopman observables and operator parameters are learned jointly from data, enabling the direct execution of the resulting circuits on a superconducting quantum processor. We validate the QKM across three progressively challeng-

2 a

Nonlinear dynamics ẋ = f (x)

Infinite-dimensional linear space

Observable functions Linear evolution in observable space ˙ ffi(x(t)) = Affi(x(t))

g1 (x)

t0 g2 (x) g3 (x)

Koopman lifting

g4 (x)

gN (x)

N-dimensional subspace

ffi(x(t0 )) ffi(x(t1 )) ffi(x(t2 )) ffi(x(tm )) g1 g2 g3 A A g4 gN

Learned finite Koopman embedding

Non-unitary evolution via h unitary channels

u(0)

Input x(0) exp(At)

Output x(t) Linear combination

Encoder Shared operator A across time

u(t)

Decoder

Quantum processor

FIG. 1. Schematic of the QKM framework for simulating nonlinear dynamics on a superconducting quantum processor. a, A nonlinear dynamical system 𝒙¤ = 𝒇 (𝒙) generates discrete spatiotemporal snapshots 𝒙(𝑡 0 ), 𝒙(𝑡1 ), · · · , 𝒙(𝑡 𝑚 ). b, The nonlinear dynamics are mapped to an infinite-dimensional linear observable space via Koopman lifting in c, and are subsequently projected onto a finite 𝑁-dimensional subspace spanned by learned observables in d, yielding the linear evolution 𝒖¤ = 𝐴𝒖. c, In full Koopman theory, this infinite-dimensional space comprises all observables 𝑔 𝑗 (𝒙) acting on the physical state, such that 𝝓(𝒙) = (𝑔1 (𝒙), 𝑔2 (𝒙), · · · ) T . Although the physical field 𝒙(𝑡) evolves nonlinearly, ¤ the lifted observables evolve linearly under the Koopman generator according to the equation 𝝓(𝒙(𝑡)) = 𝐴𝝓(𝒙(𝑡)). d, A finite Koopman embedding is learned from the data in a, yielding a low-dimensional latent state 𝒖(𝑡) that preserves the dominant dynamics. These learned observables define a finite-dimensional representation in which the same linear generator 𝐴 propagates the state through time. In practice, the physical input 𝒙(0) is encoded into the finite Koopman state 𝒖(0), evolved under the non-unitary propagator e 𝐴𝑡 , and decoded to reconstruct the physical state 𝒙(𝑡). The non-unitary propagator e 𝐴𝑡 is implemented via the LCHS, which decomposes the propagator e 𝐴𝑡 into ℎ weighted unitary channels (exemplified by the yellow, red, and purple components). These ℎ independent circuits, characterized by distinct ring-topology qubit layouts, are executed in parallel on a superconducting quantum processor.

ing regimes. These benchmarks comprise a reaction-diffusion system on a 3D cubic grid, shallow-water dynamics on a curvilinear sphere, and satellite-derived observations of the Gulf Stream. Our experiments utilize up to 32 parallel circuits of 10 physical qubits each, corresponding to 105 computational grid points. This effort represents a significant expansion in the scale of quantum resource integration for classical nonlinear dynamics, extending prior demonstrations that have been largely restricted to 1D and 2D benchmarks on fewer qubits. Together, these theoretical and experimental advances delineate the practical boundary of quantum utility for simulating nonlinear dynamics. The QKM operates effectively in the regime where the finite-dimensional Koopman representation and the shallow LCHS-derived ansatz remain faithful, attaining a speedup ratio of S ∼ O (2𝑛 /𝑛3 ) over classical Koopman propagation. Our experiments demonstrate a transition from hardware-noise-limited to theory-limited accuracy, which operationally defines the boundary of this regime.

Results Quantum-circuit realization of Koopman dynamics The Koopman representation lifts the original nonlinear dy-

namics into an 𝑁-dimensional linear system d𝑢 (𝑡) /d𝑡 = 𝐴𝑢 (𝑡), where 𝐴 ∈ C 𝑁 × 𝑁 is the finite-dimensional approximation of the Koopman operator and 𝑢 denotes the vector of learned observables (see Methods). The resulting propagator e 𝐴𝑡 is generically non-unitary, as the eigenvalues of 𝐴 can exhibit non-zero real parts associated with dissipative or unstable dynamics. To bridge the gap between this non-unitary evolution and hardware-executable quantum circuits, we establish a chain of three theorems that progressively transform e 𝐴𝑡 into an ensemble of shallow quantum circuits (see Methods). As illustrated in Fig. 2, Theorems 1 (diagonalized LCHS), 2 (spectral sampling convergence), and 3 (universal approximation bound for diagonal unitaries) establish a framework for the efficient simulation of nonlinear dynamics via ℎ independent quantum circuits, where each circuit requires only a single layer of 𝑅 𝑧 for time evolution. With the 𝑁-dimensional Koopman state encoded into 𝑛 = log2 𝑁 qubits, explicit classical propagation of e 𝐴𝑡 𝑢 (0) requires at least O (𝑁) = O (2𝑛 ) operations. The QKM replaces this dense-vector propagation by ℎ shallow circuits with O (𝑛) gates per channel, yielding an evolution-step cost of O (ℎ𝑛) and a corresponding speedup over classical Koopman propagation.

3 a

Unitary operator (Theorem 1) †

V1 eiΛ1 t V1

eAt

† V2 eiΛ2 t V2

†

Vh eiΛh t Vh

‘=1

Q1 : |0

Q2 : |0

Q3 : |0

Qn : |0

l=1 U3 (¸1;l;‘ ; ˛1;l;‘ ; ‚1;l;‘ )

l=1 U3 (¸2;l;‘ ; ˛2;l;‘ ; ‚2;l;‘ )

l=1 U3 (¸3;l;‘ ; ˛3;l;‘ ; ‚3;l;‘ )

l=1 U3 (¸n;l ;‘ ; ˛n;l ;‘ ; ‚n;l ;‘ )

h circuits (Theorem 2)

Rz („1 t)

Rz („2 t)

Rz („3 t)

Rz („n t)

Time evolution (Theorem 3)

State preparation

b “Yudu”

U U3 gate

Layer 1

U U U

U U

9 U U U

U U

8 U

U U U U

CZ gate

U U

FIG. 2. Hardware-efficient quantum circuit of the QKM. a, The non-unitary propagator e 𝐴𝑡 for nonlinear dynamics is first reformulated as a spectral integral of diagonal unitary operators (Theorem 1). Each of the ℎ independent PQCs implements a single spectral component 𝑉 (𝑘)eiΛ(𝑘 )𝑡 𝑉 † (𝑘) of the diagonalized LCHS decomposition (Theorem 2). Each 𝑛-qubit circuit comprises two functional blocks. The statepreparation block utilizes 𝑅 alternating layers of parameterized 𝑈3 gates with rotation angles {𝛼 𝑗,𝑙,ℓ , 𝛽 𝑗,𝑙,ℓ , 𝛾 𝑗,𝑙,ℓ } generated by the classical encoder, and CZ pairs that alternate between odd layers {(𝑛, 1), (2, 3), . . . , (4, 5)} and even layers {(1, 2), (3, 4), . . . , (𝑛 − 1, 𝑛)}, forming a parity-based ring topology. The time-evolution block implements a sandwich ansatz that nests a central 𝑅 𝑧 layer (Theorem 3), which encodes the diagonal eigenvalues Λ(𝑘), between 𝑈3 blocks realizing the basis transformation 𝑉 (𝑘) and its adjoint 𝑉 † (𝑘). b, Hardware-native compilation of the QKM circuit on the superconducting processor “Yudu”. The circuit ansatz in a is mapped onto a connected 10-qubit subgraph of the device such that all two-qubit gates respect the native nearest-neighbour couplings of the chip. Starting from the ground state |0⟩ ⊗10 , we apply parallel single- and two-qubit gates layer by layer up to a circuit depth of nine. In practice, each layer of parallel single-qubit (two-qubit) gates requires 40 (90) ns. Consequently, the total execution time is 560 ns, which is substantially shorter than the median qubit lifetime 48𝜇s.

To complete the simulation pipeline, the quantum evolution must be preceded by state preparation. Preparing an arbitrary 2𝑛 -dimensional quantum state requires O (2𝑛 ) quantum resources [42, 43], which would eliminate the advantage in the evolution step. Fortunately, most physically relevant initial conditions are not arbitrary and possess additional spatial structure, such as smoothness and symmetry [22]. We exploit this structure by employing a classical NN encoder that maps physical states to the rotation angles of a parameterized quantum circuit (PQC) of depth O (poly(𝑛)). The universal approximation capability of NNs [44] ensures that such a mapping is learnable, while the universality of the gate set [45] guarantees sufficient expressivity. A complementary oracle-based analysis in Sec. 4 in Supplementary Information (SI) [46] justifies the efficient state-preparation assumption for structured physical initial conditions. Numerical verification of this expressibility is provided in Sec. 8 B in SI [46]. Hardware-efficient end-to-end implementation The QKM framework is realized through a hybrid procedure comprising a classical encoder-decoder pair and ℎ parallel PQCs in Fig. 1d. Given an initial condition 𝑥(0), the encoder maps it to ℎ sets of rotation angles, with each set preparing an 𝑛-qubit quantum state on its respective PQC. Each PQC then

applies a time-evolution block to implement a single spectral component 𝑉 (𝑘)eiΛ(𝑘 )𝑡 𝑉 † (𝑘) of the diagonalized LCHS in Eq. (1). Following computational-basis measurements, the shot-count distributions from the ℎ circuits are processed by the decoder to reconstruct the evolved state 𝑥(𝑡). This workflow directly instantiates our theoretical framework: the encoder leverages physical input regularity to ensure efficient state preparation; the time-evolution block realizes the unitary representation of e 𝐴𝑡 guaranteed by Theorems 1 and 2; and the single-layer 𝑅 𝑧 configuration yields the shallow circuit prescribed by Theorem 3. Comprehensive details regarding the encoder, decoder, and circuit construction are provided in Methods and Secs. 2 and 5 in SI [46]. Crucially, the circuit structure is dictated by the LCHS decomposition rather than an empirical hardware-efficient ansatz. This systematic construction provides interpretability and error traceability while facilitating the training of the QKM. The state-preparation block (Fig. 2a) comprises 𝑅 alternating layers, each consisting of 𝑟 single-qubit gates followed by CZ entangling gates. Each time-evolution block (Fig. 2a) features a sandwich architecture, wherein a central layer of 𝑅 𝑧 rotations encoding the diagonal Λ(𝑘) is flanked by two 𝑈3 blocks realizing the basis transformation 𝑉 (𝑘) and its ad-

4 TABLE I. Quantum simulation parameters. Summary of simulation settings across the three cases, detailing the classical grid resolution D and the structural parameters of the ℎ PQCs, including the qubit count 𝑛, circuit depth 𝑅, and 𝑈3 gate density 𝑟. The theoretical evolution speedup scales as Sevo = 2𝑛 /(ℎ𝑛), while the number of measurement repetitions is determined by the shot count 𝑀. Experimental metrics for the total circuit depth 𝐿 and gate count 𝐺 are utilized to evaluate the quantum speedup S = 2𝑛 /(ℎ𝐺). In experiments, the relative theoretical error is quantified by the training loss 𝜀th ≈ ℓtrain .

Case 3D reaction-diffusion Spherical fluid Ocean currents

D 64 × 64 × 64 512 × 256 256 × 256

𝑛 6 10 10

𝑅 4 3 3

joint. In the transpiled circuit, the CZ gates alternate between odd and even layers across adjacent qubits, forming a ring topology that maps directly onto the connectivity of the superconducting processor “Yudu” in Fig. 2b. This hardware-native entangling layout eliminates SWAP insertions during transpilation and compresses the post-transpilation circuit depth by up to 50% (see Methods), which is critical for preserving simulation fidelity on NISQ devices. Furthermore, the parallelcircuit structure guaranteed by Theorem 2 can be preserved at the hardware level, with the ℎ PQCs executed concurrently on disjoint qubit subsets of the processor. To accommodate NISQ-era constraints, we set 𝑛 = O (log D), ℎ = O (𝑛), and 𝑅𝑟 = O (𝑛) for a nonlinear system with D grid points, yielding O (𝑛)-depth state preparation. The choice of ℎ balances the spectral approximation bound of Theorem 2 against resource overhead, whereas 𝑅𝑟 balances state-preparation expressibility against gate-error accumulation. Both trade-offs are detailed in Secs. 10 A and 8 B in SI [46]. The specific values of 𝑛, ℎ, 𝑅, and 𝑟 employed in each benchmark are summarized in Tab. I. Simulation case setup and hardware overview We evaluate the QKM across three physically distinct and progressively challenging regimes, each executed on the superconducting processor “Yudu” [47]. The first is a 3D reactiondiffusion system governed by the Gray-Scott equations on a cubic domain, where the competition between diffusion and nonlinear reaction yields complex, self-organizing spatiotemporal patterns. The second involves the shallow-water equations on a sphere, a canonical model in geophysical fluid dynamics, where multiscale rotational dynamics impose stringent demands on simulation accuracy. The third moves beyond prescribed equations and applies the QKM to satellite-derived Gulf Stream observations, converting real geophysical data into a quantum-executable model of ocean-current dynamics. The superconducting processor “Yudu” comprises 72 frequency-tunable transmon qubits arranged in a lattice geometry, with comprehensive device specifications provided in Sec. 6 in SI [46]. In this work, we utilize up to 10 qubits per subcircuit and up to 32 parallel subcircuits to perform the experiments. The target experimental circuits are implemented using the native gate set {𝑈3 , CZ}. Through optimized control procedures, we achieve parallel single- and two-qubit gate fidelities of 99.92% and 98.82%, respectively. Experimental results are reconstructed by the NN decoder from the computational-basis measurements on all qubits of the ℎ

𝑟 1 3 3

ℎ 8 32 8

𝑀 6144 10240 10240

𝐿 12 17 17

𝐺 58 146 145

𝜀 th 0.001 0.015 0.020

Sevo 1.33 3.2 12.8

S 0.14 0.22 0.88

parallel circuits, with each circuit executed for 𝑀 shots. The reference solutions are obtained from direct numerical simulations for the 3D reaction–diffusion and spherical fluid dynamics cases, and from satellite-derived altimetry observations for the Gulf Stream case. Tab. I summarizes the problem sizes, circuit configurations and experimental parameters used in the three benchmarks. 3D reaction-diffusion systems The first test case targets the 3D Gray-Scott reaction-diffusion system [48] (see Sec. 7 A in SI [46]). In this system, the competition between chemical reaction and spatial diffusion triggers the emergence of complex, self-organizing dissipative structures. This benchmark assesses the capacity of the QKM to encode high-dimensional (D = 643 ) nonlinearities into a compact quantum representation using 𝑛 = 6 qubits per circuit. The QKM successfully simulates the 3D reaction-diffusion dynamics, characterized by the emergence, growth, and coalescence of self-organized dissipative structures of the reactant-concentration field 𝑢, as shown in Fig. 3a. Quantitatively, the total energy tracks the dissipative relaxation of the system (Fig. 3b), and the relative error 𝜀 𝐿2 remains below 0.05 and the median Hellinger fidelity is maintained near 0.93 throughout the evolution (Fig. 3c). These results demonstrate that the diagonalized unitary dynamics established in Theorems 1 to 3 enable physically consistent 3D spatiotemporal simulations on a NISQ device. Spherical fluid dynamics The second benchmark extends the evaluation to spherical fluid dynamics governed by the shallow-water equations [49] (see Sec. 7 B in SI [46]). This case introduces topological constraints absent in Euclidean domains, specifically periodic longitudinal closure and polar coordinate singularities. The resulting dynamics, which include mid-latitude jet instability, planetary-scale vortex shedding, and turbulent spectral cascades, exhibit a heightened sensitivity to phase errors and discretization artifacts compared to their flat-domain counterparts. At a grid resolution of D = 512 × 256 encoded using 𝑛 = 10 qubits per circuit and ℎ = 32 parallel PQCs, this case constitutes the largest circuit deployment and the most geometrically complex test within this study. The QKM simulates the long-time vorticity evolution under spherical shallow-water dynamics in Fig. 3d. The quantum simulation captures the growth of the mid-latitude instability and the emergence of large-scale coherent vortices on the sphere, while quantum noise manifests primarily as localized

100

b t

300

High Reactant concentration

Ref.

Exp.

0.5 0.4 0.3 0.2 0.1

Exp. Ref.

1.00

0.08

0.93

0.06

0.86 Relative error Hellinger fidelity

Spherical fluid dynamics 40 60 80

e 0.1 100 t

Exp. Vorticity

0:3

0.72 0

0.8

0.72 0.64

0.06

Ref.

−0:3

Relative error Hellinger fidelity

0.08

0.04

0.56

0.02

0.48

0.4

0.79

0.02 0

100 150 200 250 300 Time

0.1

0.04 Low

40 80 Time

100 150 200 250 300 Time

106 E! 102 10−2

0.65

Exp. Ref.

100

101 k

103

KDE

3D reaction-diffusion system 200

Total energy

102 Exp. Ref.

101 10−1 10−3

10−2

10−1

100

FIG. 3. Assessment of the QKM framework on benchmark nonlinear dynamical systems. a, Quantum simulation of a 3D reaction-diffusion system using the QKM on the superconducting processor “Yudu”. Volume renderings compare the experimental reactant-concentration fields 𝑢 with the reference solutions at 𝑡 = 100, 200, and 300. The QKM result accurately captures the temporal evolution of the total energy ⟨𝑢 2 ⟩/2 in b, where ⟨·⟩ denotes volume average. The relative error 𝜀 𝐿2 remains small, and the median Hellinger fidelity stays close to unity in c, demonstrating quantitative agreement with the nonlinear dynamics. d, Quantum simulation of spherical fluid dynamics using the QKM on the superconducting processor “Yudu”. The experimental vorticity fields 𝜔 are compared with the reference solutions at 𝑡 = 0, 20, 40, 60, 80, and 100. The relative error remains below 10%, and the median Hellinger fidelity stays near 0.65 throughout the simulation in e. The comparison of the experimental enstrophy spectrum 𝐸 𝜔 and vorticity KDE with the reference solutions at 𝑡 = 80 in f confirms the accuracy of the QKM.

fluctuations in the vorticity snapshots. As illustrated in Fig. 3e, the median Hellinger fidelity stabilizes near 0.6, whereas the relative error grows at a moderate rate. The enstrophy spectrum 𝐸 𝜔 closely matches the multiscale distribution of vorticity across wavenumbers, while the kernel density estimate (KDE) of the vorticity recovers consistent one-point statistics, as shown in Fig. 3f. These diagnostic assessments demonstrate that the QKM maintains the core vortical, spectral, and statistical properties of spherical fluid dynamics. Moreover, the close agreement of the enstrophy spectrum over a broad range of scales indicates that the quantum simulation preferentially preserves the core dynamics of the system. The principal low-wavenumber modes containing the bulk of the enstrophy are encoded with high fidelity, whereas discrepancies are restricted to the high-wavenumber tail where both Koopman truncation and quantum noise accumulate. This spectral selectivity is consistent with the complexity analysis, which assumes that many physically relevant systems are dominated by low-frequency components that admit efficient Koopman truncation. A comparison with a noiseless simulation indicates that the theoretical error 𝜀th and hardware noise

𝜀 noise are of comparable magnitude in this case (detailed in Sec. 9 B in SI [46]). Real-world ocean currents The final benchmark evaluates the QKM on the complex mesoscale dynamics of the Gulf Stream using real-world observational data [50] (see Sec. 7 C in SI [46]). Unlike previous benchmarks governed by closed-form partial differential equations under controlled initial conditions, this scenario utilizes satellite-derived altimetry products characterized by observational gaps and irregular land-sea boundaries. Simulating this system tests the capacity of the QKM to function not merely as an equation solver, but as a data-driven framework that assimilates observational data into a quantum-executable form. The Gulf Stream region (20◦ N to 52◦ N, 33◦ W to 65◦ W) is discretized at a resolution of D = 256 × 256 and encoded using 𝑛 = 10 qubits with ℎ = 8 parallel PQCs. A comparison between the satellite-derived data in Fig. 4a and the QKM results in Fig. 4b demonstrates that the QKM successfully simulates the principal mesoscale features of the surface geostrophic velocity field 𝑣, including jet and eddy

6 2024-01-07

1010

2024-01-13 v

1.8

52◦ N

0.6

Exp.

0.0 20◦ N

0.4

65◦ W

33◦ W

65◦ W

33◦ W

65◦ W

1.0

0.2

0.6

0.0

0.2

01-01 01-03 01-05 01-07 01-09 01-11 01-13 Date

FIG. 4. Quantum simulation of real-world Gulf Stream ocean currents on the superconducting processor “Yudu”. a, Surface geostrophic velocity 𝑣 derived from satellite observations on three dates spanning January 1 to 13, 2024. b, Corresponding QKM results simulated by the “Yudu” processor. c, Temporal evolution of the relative error 𝜀 𝐿2 and the median Hellinger fidelity, with box plots showing the fidelity distribution across circuit instances on January 1, 7, and 13, 2024.

structures, over a 13-day evaluation period in January 2024. As shown in Fig. 4c, the error 𝜀 𝐿2 exhibits a gradual growth characteristic of chaotic dynamical systems, yet remains bounded at approximately 0.2 after nearly two weeks of evolution. Throughout the evaluation window, the median Hellinger fidelity 𝐹 remains stable at approximately 0.6, indicating robust statistical agreement between the hardware-sampled and noiseless quantum distributions. These findings demonstrate that the QKM possesses the expressive capacity and noise resilience required to manage the inherent uncertainties of real-world geophysical observations. The growth of the 𝐿 2 error, in contrast with the temporally stable Hellinger fidelity, suggests that while the quantum processor continues to faithfully execute the learned circuits, the finite-dimensional Koopman observable space gradually loses its capacity to track the chaotic trajectory over extended horizons. Consequently, expanding the dimension 𝑁 of the observable space represents the primary avenue for improving long-term predictive accuracy. Together, these three benchmarks demonstrate that the QKM provides a unified and hardware-validated framework for the digital quantum simulation of nonlinear dynamics. The QKM maintains physically consistent accuracy on a superconducting processor without requiring modification to the circuit topology or the LCHS formulation, confirming the portability and practical viability of the underlying framework. Pathway to quantum utility with QKM Although the QKM speedup S ∼ O (2𝑛 /𝑛3 ) derived in Methods (blue line in Fig. 5) is attainable for any nonlinear system

(iii) QKM-prohibitive (i)

∼

(ii)

102

ent

rk wo

rr Cu

(iii) 10−2

st Be

se ca

(ii) QKM-intermediate

104

100

33◦ W

Hellinger fidelity Relative error

106

(i) QKM-amenable

nonlinearity sing rea Inc

1.2

20◦ N

108 Quantum speedup, S

2024-01-01

Ref.

a 52◦ N

Worst case S ∼ 1/n

20 30 Number of qubits, n

FIG. 5. Quantum-speedup regimes of the QKM. The quantum speedup S as a function of the number of qubits 𝑛. The upper blue curve represents the best-case scaling, S ∼ 2𝑛 /𝑛3 , achieved when a first-order circuit ansatz corresponding to a single 𝑅 𝑧 layer already provides sufficient accuracy. Problems in this limit are classified as QKM-amenable. With increasing nonlinearity of the target dynamics, higher-order QKM representations and more retained terms are required, thereby increasing the quantum cost and lowering the achievable speedup. The region between the best-case curve and the threshold S = 1 defines the QKM-intermediate regime, where the first-order approximation is insufficient, yet the quantum implementation remains faster than the corresponding classical Koopman propagation. The dashed line S = 1 marks the crossover between quantum advantage and the absence of advantage. Below this threshold, problems enter the QKM-prohibitive regime, where the required accuracy approaches the full expansion limit, the worst-case scaling (red curve) decreases to S ∼ 1/𝑛. The star indicates the scale of the present experiments with 𝑛 = 10, near the S = 1 utility threshold.

using the single-layer 𝑅 𝑧 ansatz, the accuracy at which this ceiling is reached varies across systems. The performance metric is the relative theoretical error 𝜀 th , which combines the approximation error of the finite-dimensional Koopman representation and that of the single-layer 𝑅 𝑧 ansatz. In our experiments, the ansatz approximation error 𝜀 ansatz due to the limited dimension 𝑁 dominates 𝜀 th (see Sec. 10 B in SI [46]). Theorem 3 therefore establishes the scalability limit, bounding 𝜀 ansatz by the high-order Pauli-𝑍 string coefficients of the target Hamiltonian. Retaining higher-order terms enriches ansatz and reduces 𝜀th by capturing additional spectral content, albeit at the expense of a reduced S. The choice of ansatz therefore defines an accuracy-speedup trade-off curve, ranging from the shallow limit at S ∼ O (2𝑛 /𝑛3 ) to the deep limit of exact representation at S ∼ O (1/𝑛). While any system can be simulated using shallow circuits, complex systems must trade speedup S for accuracy 𝜀 th . We adopt a threshold of 𝜀 ∗ = 10−3 on 𝜀 th as the operational criterion, quantified in practice by the training loss ℓtrain . This trade-off partitions nonlinear dynamical systems into three distinct regimes based on their spectral compressibility under a learned Koopman embedding. For QKM-amenable systems, the single-layer 𝑅 𝑧 ansatz achieves an error of 𝜀 th ≲ 𝜀 ∗ , enabling the speedup S ∼ O (2𝑛 /𝑛3 ) to be realized

7 at scientific-computing accuracy. For QKM-intermediate systems, the shallow ansatz yields 𝜀 th > 𝜀 ∗ , but augmenting the ansatz with 𝑅 𝑧𝑧 and higher-order gates reduces 𝜀th at a polynomial cost in 𝑛. In this regime, quantum advantage remains theoretically viable at a reduced speedup and becomes practically accessible once gate fidelities improve sufficiently to support the deeper circuits. Conversely, QKM-prohibitive systems are characterized by intrinsically incompressible spectral dynamics, so that no finite-dimensional Koopman embedding can concentrate their spectral weight into low-order modes. A canonical example is fully-developed turbulence [41], in which the energy cascade distributes spectral weight across all wavenumbers. Capturing such dynamics demands an ansatz whose cost scales exponentially with 𝑛, degrading the worstcase performance to S ∼ O (1/𝑛) (red line in Fig. 5) and eliminating the quantum advantage. Representing these systems instead requires tailored methods focused on extracting dominant coherent structures [22]. The three benchmarks calibrate this partition. With measured training losses ℓtrain ranging from 10−3 to 2 × 10−2 (see Tab. I), the 3D reaction-diffusion benchmark falls within the QKM-amenable regime, whereas the spherical shallow-water and Gulf Stream benchmarks lie in the QKM-intermediate regime. Comparing hardware execution against a noiseless simulation (Sec. 9 in SI [46]) resolves the dominant error source in each case, revealing a progression from hardwarenoise-limited to theory-limited performance as nonlinearity grows. Critically, by learning observables directly from data rather than relying on an analytical expansion, the QKM maps dynamics with broadband spectral content into the reach of shallow circuits, including mesoscale geophysical flows that are beyond the scope of truncation-based linearization. The framework thereby establishes a broader operational envelope for simulating general nonlinear problems.

Discussion The QKM unifies the Koopman operator theory, Hamiltonian simulation, and hardware-native compilation into a single jointly trained pipeline, recasting non-unitary Koopman dynamics into shallow parallel circuits on a superconducting processor. The scaling of speedup ratio S ∼ O (2𝑛 /𝑛3 ) establishes the performance ceiling of the QKM; the regime criterion 𝜀 th ≲ 𝜀 ∗ specifies when a given system achieves this bound under the single-layer 𝑅 𝑧 ansatz. Implemented on up to 32 parallel 10-qubit circuits, the 3D reaction-diffusion benchmark falls within the QKM-amenable regime, establishing that this performance ceiling is not merely asymptotic but practically reachable on NISQ devices. The spherical shallow-water and Gulf Stream benchmarks exceed this threshold under the shallow ansatz, thereby mapping the boundary behavior predicted by the QKM and indicating where richer ansatzes are required to further reduce 𝜀 th . These results are achieved using a structured circuit ansatz derived from the LCHS decomposition rather than generic variational heuristics. The global representational capacity of the QKM stems from integrating Koopman lifting and an NN encoder, which jointly learn a finite-dimensional observable subspace that linearizes the dynamics. The LCHS decomposition then translates this representation into physically mo-

tivated circuits, while a topology-native layout preserves this structure on hardware without requiring SWAP insertions. Together, these features provide interpretability and error traceability while encoding underlying physical structures beyond those of empirical ansatz. Analytical embeddings, such as Carleman [11] and Koopman-von Neumann [12, 13] formulations, are typically restricted to regimes of weak nonlinearity and strong dissipation. By contrast, learning observables from data identifies an invariant subspace adapted to the dynamics, extending the accessible range to the moderately nonlinear regimes demonstrated in this work. However, several limitations remain. First, the quantumclassical input/output bandwidth is the primary obstacle to scaling up the practical speedup. The input/output latency and noise accumulation preclude efficient on-device training, so the QKM is trained classically in a one-time offline preprocessing phase before deployment on quantum hardware. This restriction constrains the accessible observable dimension, limiting the present implementation to reduced-dimensional demonstrations and deferring the validation of practical advantages at larger qubit scales to future studies. Second, strongly multi-scale dynamics, such as turbulence [51, 52], challenge both the shallow ansatz and the encoder, because high-order Pauli-𝑍 interactions increase the ansatz error while the smooth bias of the NN encoder can attenuate fine-scale information. Finally, establishing rigorous bounds on the projection error remains an open challenge for data-driven Koopman methods. Extending the operational reach of the QKM requires the coordinated development of both hardware and algorithms. Specifically, improved gate fidelities enable the reliable execution of deeper circuits, while broader qubit connectivity facilitates native, multi-qubit entangling operations. These hardware advancements pave the way for augmenting the timeevolution block with higher-order entangling gates, rendering the QKM-intermediate regime practically accessible. On the algorithmic side, learning observables that concentrate spectral weight into lower-order Pauli-𝑍 interactions can transition a QKM-intermediate system into the QKM-amenable regime, providing quantum utility without requiring deeper circuits. Finally, refining the theoretical bounds and exploring early fault-tolerant architectures will further clarify the boundary between QKM-amenable and QKM-prohibitive systems, ultimately enabling quantum utility in simulating complex, realworld physical phenomena.

Methods Koopman operator theory The Koopman operator theory reformulates nonlinear dynamics as linear time evolution in an infinite-dimensional function space [53]. Consider a smooth dynamical system d𝑥 (𝑡) /d𝑡 = 𝑓 [𝑥 (𝑡)], where 𝑥 ∈ X represents the state on a manifold X ⊂ R D embedded in the D-dimensional real space R D , and 𝑓 is a nonlinear vector field. Its solution defines a flow map Φ𝑡 : X → X, such that 𝑥(𝑡) = Φ𝑡 [𝑥(0)] for an initial condition 𝑥(0). Instead of directly evolving the state 𝑥, the Koopman framework considers the evolution of an observable 𝑔 : X → C belonging to a function space G(X), where C denotes the complex space. The Koopman operator K 𝑡 acts on this ob-

8 servable via K 𝑡 𝑔 (𝑥) = 𝑔[Φ𝑡 (𝑥)]. Although the flow map Φ𝑡 may be nonlinear, K 𝑡 is linear by construction. Consequently, the infinitesimal generator A of K 𝑡 induces a linear dynamical system d𝑔/d𝑡 = A𝑔 [54, 55], where A represents the Lie derivative of 𝑔 along 𝑓 (𝑥), defined as A𝑔 := lim𝑡→0 (K 𝑡 𝑔 − 𝑔)/𝑡 = ∇𝑔 · 𝑓 . Although the observable space G(X) is inherently infinitedimensional, practical computation necessitates projecting the dynamics onto a finite-dimensional subspace. Identifying the optimal dimensionality for such a projection remains an open challenge in Koopman analysis [56]. In the QKM, the state 𝑥 is mapped to a vector of 𝑁 basis functions 𝑢 = [𝑔1 (𝑥) , · · · , 𝑔 𝑁 (𝑥)] T ∈ C 𝑁 (see Fig. 1c), yielding an approximate linear dynamical system d𝑢 (𝑡) /d𝑡 = 𝐴𝑢 (𝑡) , 𝑢 (0) = 𝑢 0 , where 𝐴 ∈ C 𝑁 × 𝑁 represents the finite-dimensional approximation of A. Both the observable functions {𝑔𝑖 } and the operator 𝐴 are learned jointly from data via end-to-end training of the NN encoder and the quantum circuit parameters. Theorems for unitary mapping Theorem 1 (Diagonalized LCHS). Let 𝐴 ∈ C 𝑁 × 𝑁 be decomposed into Hermitian and anti-Hermitian parts, 𝐴 = 𝐿 + i𝐻, with 𝐿 = ( 𝐴 + 𝐴† )/2 and 𝐻 = ( 𝐴 − 𝐴† )/2i. If 𝐿 ⪯ 0, then for 𝑡 ⩾ 0, the propagator e 𝐴𝑡 admits the spectral integral representation ∫ 1 𝑉 (𝑘)eiΛ(𝑘 )𝑡 𝑉 † (𝑘) d𝑘, (1) e 𝐴𝑡 = 2 R 𝜋(1 + 𝑘 ) where 𝐻+𝑘 𝐿 = 𝑉 (𝑘)Λ(𝑘)𝑉 † (𝑘) is the spectral decomposition of the Hermitian operator 𝐻 + 𝑘 𝐿 for each 𝑘 ∈ R. Theorem 1 reformulates the non-unitary propagator as a spectral integral of diagonal unitary operators, demonstrating the feasibility of a unitary representation. The diagonal matrix Λ(𝑘) reduces the time evolution within each spectral component to phase rotations on individual basis states, a simplification exploited in Theorem 3 to guarantee a shallow circuit depth. The condition 𝐿 ⪯ 0 is satisfied without loss of generality via a rescaling of variables (see Sec. 1 A in SI [46]). Theorem 2 (Spectral sampling convergence). Let 𝑝(𝑘) = [𝜋(1 + 𝑘 2 )] −1 be the Cauchy-Lorentz density and 𝑈 𝑘 (𝑡) = ei(𝐻+𝑘 𝐿)𝑡 = 𝑉 (𝑘)eiΛ(𝑘 )𝑡 𝑉 † (𝑘) be the unitary operator defined in Theorem 1. For an ℎ-point uniform approximation over the truncated interval [−𝐾, 𝐾] with spacing Δ𝑘 = 2𝐾/ℎ, the error is bounded by ∥e

𝐴𝑡

−

ℎ ∑︁ 𝑗=1

√ 1/3 3 3 3 + 8𝑡 ∥𝐿 ∥ 2 Δ𝑘 · 𝑝(𝑘 𝑗 )𝑈 𝑘 𝑗 (𝑡) ∥ 2 ⩽ , (2) 𝜋 4ℎ

where the optimal truncation radius 𝐾 = minimizes this upper bound.

√ 3 3+8𝑡 ∥ 𝐿 ∥ 2 −1/3 4ℎ

Theorem 2 replaces the spectral integral in Eq. (1) with a finite sum over ℎ points (see Sec. 1 B in SI [46]), where ℎ serves as a hyperparameter controlling the approximation fidelity. This decomposition of a deep circuit into ℎ shallow ones is essential for near-term hardware. Although the ℎ spectral components could be encoded into a single log2 (ℎ𝑁)-qubit

circuit through entanglement, doing so would demand significantly deeper circuits with controlled multi-qubit operations. Distributing the computation across ℎ independent 𝑛 = log2 𝑁 qubits keeps each circuit shallow. Since these circuits share no quantum state, they can be executed concurrently to disjoint regions of the processors, as illustrated in Fig. 1e. Theorem 3 (Universal approximation bound for diagonal uni𝑛 𝑛 taries). Let 𝐻 ∗ ∈ C2 ×2 be a diagonal Hamiltonian with the ∗ i𝐻 ∗ . The single-layer 𝑅 subspace induced unitary 𝑈 = eË 𝑧 𝑛 i𝜙 is defined as F = {e 𝑗=1 𝑅 𝑧 (𝜃 𝑗 ) | 𝜃 𝑗 , 𝜙 ∈ R}, where 𝑅 𝑧 (𝜃) = e−i𝜎𝑧 𝜃/2 and 𝜎𝑧 is the Pauli-𝑍 operator. For any 𝑈 ∈ F , the minimum approximation error satisfies ∑︁ 1 ∥𝑈 − 𝑈 ∗ ∥ 2𝐹 = 𝛼𝑠2 + O (max 𝛼𝑠4 ) 𝑛 𝑈∈ F 2 |𝑠 |⩾2 min

(3)

|𝑠 |⩾2

where 𝑠 = 𝑠1 . . . 𝑠 𝑛 ∈ {0, 1} 𝑛 is a bitstring with Hamming weight |𝑠|, ∥ · ∥ 𝐹 denotes the Frobenius norm, and 𝛼𝑠 = 2−𝑛 Tr(𝐻 ∗ 𝑃𝑠 ) are the Pauli coefficients ofË 𝐻 ∗ associ𝑛 𝑠𝑘 ated with the multi-qubit Pauli-𝑍 operators 𝑃𝑠 = 𝑘=1 𝜎𝑧 . Although the diagonal unitaries yielded by Theorem 1 are compatible with quantum architectures, their exact implementation requires multi-qubit Pauli-𝑍 interactions of all orders, rendering them prohibitively error-prone on NISQ devices. Theorem 3 establishes that any diagonal Hamiltonian admits a unique decomposition into multi-body Pauli-𝑍 operators (see Sec. 1 C in SI [46]), and demonstrates that a single layer of 𝑅 𝑧 rotations captures the zeroth- and first-order terms exactly, leaving a residual error 𝜖 ansatz determined by the second- and higher-order coefficients. This bound implies that systems with spectral content concentrated in low-order interactions admit highly faithful approximations. Although Theorem 3 ∗ addresses the static unitary 𝑈 ∗ = ei𝐻 , the bound extends to ∗ the time-dependent evolution operator ei𝐻 𝑡 . Consequently, 𝜖 ansatz grows with 𝑡, a trend consistent with the experimental observations presented below. Experimental setup All quantum experiments are executed on the superconducting processor “Yudu”, accessed via the Quafu cloud platform [47], deploying 8 × 6-qubit subcircuits for 3D reaction-diffusion, 32 × 10-qubit for spherical fluid dynamics, and 8 × 10-qubit for real-world ocean current simulations. Because the ringstructured entangling topology of the circuit ansatz maps directly onto the processor (see Sec. 6 in SI [46]), transpilation requires no SWAP gate insertions. Consequently, minor practical adjustments reduce the circuit depth from 12 to 9 layers for 3D reaction-diffusion systems, and from 17 to 9 for spherical fluid and ocean currents simulations, achieving a reduction of nearly 50% in the most complex case. This demonstrates that the ansatz is hardware-friendly and well suited to NISQ-era execution. The ℎ independent PQCs are mapped to spatially disjoint qubit subsets, enabling concurrent execution within a single job submission, thereby maximizing hardware throughput. Quantum resource allocation and processor specifications are detailed in Tab. I and Sec. 6 in SI [46], respectively. Throughout this work, we distinguish relative errors, denoted by 𝜀, from absolute errors, denoted by 𝜖. A detailed

9 error analysis is provided in Sec. 8 in SI [46]. We assess simulation performance using two metrics. The relative 𝐿 2 error 𝜀 𝐿2 = (∥ 𝑥ˆ 𝑘 − 𝑥 𝑘 ∥ 22 )/∥𝑥 𝑘 ∥ 22 quantifies the deviation between the QKM prediction 𝑥ˆ 𝑘 and the reference 𝑥 𝑘 at a discrete time step 𝑘. The error induced by quantum is evaluated Pnoise √ using the Hellinger fidelity 𝐹 (𝑄, 𝑊) = ( 𝑖 𝑞 𝑖 𝑤 𝑖 ) 2 , which measures the statistical overlap between the experimentally sampled distribution 𝑄 and the noise-free ideal distribution 𝑊. Fidelity values close to unity indicate that the hardware execution more closely reproduces the ideal noise-free quantum state. Complexity analysis The 𝑁-dimensional Koopman observable space is encoded into 𝑛 = log2 𝑁 qubits. For each of the ℎ components in the discretized LCHS integral of Eq. (1), the quantum circuit comprises a state-preparation block and a time-evolution block. The gate complexity for preparing the initial state of ℎ PQCs is Cprep = O (ℎ𝑛𝑅𝑟 + ℎ𝑛𝑅/2 + 4ℎ𝑛) = O (ℎ𝑛𝑅𝑟). The gate complexity for the time evolution of ℎ PQCs is Cevo = O (ℎ𝑛). The total gate complexity is therefore Cqc = Cprep + Cevo = O (ℎ𝑛𝑅𝑟). We compare this cost against classical propagation of the same Koopman dynamics. Evaluating the matrix exponential e 𝐴𝑡 𝑢(0) classically incurs a cost Ccc that is optimally O (𝑁) for sparse 𝐴, and ranges from O (𝑁 2 ) to O (𝑁 3 ) for dense exponentiation. Taking the sparse case as the conservative reference, the QKM enables a theoretical evolution speedup of Sevo = Ccc /Cevo = O (2𝑛 /(ℎ𝑛)). Under the scaling 𝑅𝑟 = O (𝑛) and ℎ = O (𝑛) (see Sec. 8 B and 10 A in SI [46]), the resulting quantum speedup is S = Ccc /Cqc = O (2𝑛 /𝑛3 ). Standard projective measurements over 𝑀 shots incur a statistical error √ 𝜖meas = O (1/ 𝑀), yielding an end-to-end speedup of Stotal = 2 Ccc /(𝑀 Cqc ) = O (2𝑛 𝜖meas /𝑛3 ). This scaling indicates that a practical quantum speedup is achievable on a sufficiently large 2 observable space of dimension 𝑁 ≳ O (1/𝜖meas ). Although quantum execution is efficient, the framework requires a one-time classical pre-processing phase to train NNs and the unitary operator parameters. This training stage, performed on a classical computer, optimizes these parameters to encode the nonlinear dynamics into the circuit rotation angles.

Data Availability The data presented in the figures and that support the other findings of this study will be publicly available upon its publication.

Code Availability The source code has been deposited (https://github.com/YYgroup/QKM) [57].

QKM

Acknowledgments This work has been supported by the National Natural Science Foundation of China (Grant Nos. 12525201, 12432010, 12588201, and 52306126), and the Beijing Natural Science Foundation (Grant No. F261001).

Author Contributions B.Z., Z.L., and Y.Yang conceived the theoretical idea. B.Z.

and Z.L. developed the quantum Koopman method. B.Z. conducted the quantum simulations. D.A. and Z.M. contributed to the theoretical analysis. Y.Yu and X.X. provided the quantum hardware support. Y.Yang supervised the project. All authors contributed to data analysis, discussion of the results, and writing of the manuscript.

REFERENCES [1] R. P. Feynman, Simulating physics with computers, Int. J. Theor. Phys. 21, 467 (1982). [2] A. J. Daley, I. Bloch, C. Kokail, S. Flannigan, N. Pearson, M. Troyer, and P. Zoller, Practical quantum advantage in quantum simulation, Nature 607, 667 (2022). [3] M. A. Nielsen and I. L. Chuang, Quantum computation and quantum information (Cambridge University Press, 2010). [4] N. Beech, T. Rackow, T. Semmler, S. Danilov, Q. Wang, and T. Jung, Long-term evolution of ocean eddy activity in a warming world, Nat. Clim. Chang. 12, 910 (2022). [5] P. Givi, A. J. Daley, D. Mavriplis, and M. Malik, Quantum speedup for aeroscience and engineering, AIAA J. 58, 3715 (2020). [6] S. Succi, W. Itani, K. Sreenivasan, and R. Steijl, Quantum computing for fluids: where do we stand?, Europhys. Lett. 144, 10001 (2023). [7] Z. Meng, C. Song, and Y. Yang, Challenges of simulating fluid flows on near-term quantum computer, Sci. China Phys., Mech. Astron. 68, 104705 (2025). [8] F. Tennie, S. Laizet, S. Lloyd, and L. Magri, Quantum computing for nonlinear differential equations and turbulence, Nat. Rev. Phys. 7, 220 (2025). [9] T. Hoefler, T. Häner, and M. Troyer, Disentangling hype from practicality: on realistically achieving quantum advantage, Commun. ACM 66, 82 (2023). [10] S. Aaronson, A. M. Childs, E. Farhi, A. W. Harrow, and B. C. Sanders, Future of quantum computing, Quantum Mach. Intell. 8, 3 (2026). [11] 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, Proc. Natl. Acad. Sci. U. S. A. 118, e2026805118 (2021). [12] I. Joseph, Koopman-von Neumann approach to quantum simulation of nonlinear classical dynamics, Phys. Rev. Res. 2, 043102 (2020). [13] S. Jin, N. Liu, and Y. Yu, Time complexity analysis of quantum algorithms via linear representations for nonlinear ordinary and partial differential equations, J. Comput. Phys. 487, 112149 (2023). [14] S. Succi, W. Itani, C. Sanavio, K. R. Sreenivasan, and R. Steijl, Ensemble fluid simulations on quantum computers, Comput. Fluids 270, 106148 (2024). [15] H. Alipanah, F. Zhang, Y.-X. Yao, R. Thompson, N. Nguyen, J. Liu, P. Givi, B. J. McDermott, and J. J. Mendoza-Arenas, Quantum dynamics simulation of the advection-diffusion equation, Phys. Rev. Res. 7, 043318 (2025). [16] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum algorithm for linear systems of equations, Phys. Rev. Lett. 103, 150502 (2009). [17] A. M. Childs, R. Kothari, and R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision, SIAM J. Comput. 46, 1920 (2017). [18] D. An, J.-P. Liu, and L. Lin, Linear combination of Hamil-

10 tonian simulation for nonunitary dynamics with optimal state preparation cost, Phys. Rev. Lett. 131, 150603 (2023). [19] S. Jin, N. Liu, and Y. Yu, Quantum simulation of partial differential equations via Schrödingerization, Phys. Rev. Lett. 133, 230602 (2024). [20] Z. Lu and Y. Yang, Quantum computing of reacting flows via Hamiltonian simulation, Proc. Combust. Inst. 40, 105440 (2024). [21] P. Brearley and S. Laizet, Quantum algorithm for solving the advection equation using Hamiltonian simulation, Phys. Rev. A 110, 12430 (2024). [22] Z. Meng, X. Zhang, X. Yuan, and Y. Yang, Geometric encoding of turbulence for end-to-end quantum simulation, arXiv preprint arXiv:2508.05346 (2025). [23] Z. Meng and Y. Yang, Quantum computing of fluid dynamics using the hydrodynamic Schrödinger equation, Phys. Rev. Res. 5, 033182 (2023). [24] Z. Meng and Y. Yang, Quantum spin representation for the Navier-Stokes equation, Phys. Rev. Res. 6, 043130 (2024). [25] B. Wang, Z. Meng, Y. Zhao, and Y. Yang, Quantum lattice Boltzmann method for simulating nonlinear fluid dynamics, npj Quantum Inf. 11, 196 (2025). [26] M. Lee, Z. Song, S. Kocherla, A. Adams, A. Alexeev, and S. H. Bryngelson, A multiple-circuit approach to quantum resource reduction with application to the quantum lattice Boltzmann method, Future Gener. Comput. Syst. 174, 107975 (2026). [27] J. Preskill, Quantum computing in the NISQ era and beyond, Quantum 2, 79 (2018). [28] A. A. Mele, A. Angrisani, S. Ghosh, S. Khatri, J. Eisert, D. Stilck França, and Y. Quek, Noise-induced shallow circuits and the absence of barren plateaus, Nat. Phys. , 1 (2026). [29] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Quantum machine learning, Nature 549, 195 (2017). [30] M. Lubasch, J. Joo, P. Moinier, M. Kiffner, and D. Jaksch, Variational quantum algorithms for nonlinear problems, Phys. Rev. A 101, 010301 (2020). [31] M. Cerezo, G. Verdon, H.-Y. Huang, L. Cincio, and P. J. Coles, Challenges and opportunities in quantum machine learning, Nat. Comput. Sci. 2, 567 (2022). [32] P. Pfeffer, F. Heyder, and J. Schumacher, Reduced-order modeling of two-dimensional turbulent Rayleigh-Bénard flow by hybrid quantum-classical reservoir computing, Phys. Rev. Res. 5, 043242 (2023). [33] D. Jaksch, P. Givi, A. J. Daley, and T. Rung, Variational quantum algorithms for computational fluid dynamics, AIAA J. 61, 1885 (2023). [34] M. Larocca, S. Thanasilp, S. Wang, K. Sharma, J. Biamonte, P. J. Coles, L. Cincio, J. R. McClean, Z. Holmes, and M. Cerezo, Barren plateaus in variational quantum computing, Nat. Rev. Phys. , 1 (2025). [35] E. R. Anschuetz and B. T. Kiani, Quantum variational algorithms are swamped with traps, Nat. Commun. 13, 7760 (2022). [36] S. S. Bharadwaj and K. R. Sreenivasan, Hybrid quantum algorithms for flow problems, Proc. Natl. Acad. Sci. U. S. A. 120, e2311014120 (2023). [37] L. Wright, C. Mc Keever, J. T. First, R. Johnston, J. Tillay, S. Chaney, M. Rosenkranz, and M. Lubasch, Noisy intermediate-scale quantum simulation of the one-dimensional wave equation, Phys. Rev. Res. 6, 043169 (2024). [38] Z. Meng, J. Zhong, S. Xu, K. Wang, J. Chen, F. Jin, X. Zhu, Y. Gao, Y. Wu, C. Zhang, et al., Simulating unsteady flows on a superconducting quantum processor, Commun. Phys. 7, 349 (2024). [39] Z.-Y. Chen, T.-Y. Ma, C.-C. Ye, L. Xu, W. Bai, L. Zhou, M.-

Y. Tan, X.-N. Zhuang, X.-F. Xu, Y.-J. Wang, et al., Enabling large-scale and high-precision fluid simulations on near-term quantum computers, Comput. Methods Appl. Mech. Eng. 432 (2024). [40] Z. Wang, J. Zhong, K. Wang, Z. Zhu, Z. Bao, C. Zhu, W. Zhao, Y. Zhao, Y. Yang, C. Song, et al., Simulating fluid vortex interactions on a superconducting quantum processor, Nat. Commun. 17, 2602 (2026). [41] B. Zhang, Z. Lu, Y. Zhao, and Y. Yang, Data-driven quantum Koopman method for simulating nonlinear dynamics, preprint arXiv:2507.21890 (2025). [42] X.-M. Zhang, T. Li, and X. Yuan, Quantum state preparation with optimal circuit depth: implementations and applications, Phys. Rev. Lett. 129, 230504 (2022). [43] X. Sun, G. Tian, S. Yang, P. Yuan, and S. Zhang, Asymptotically optimal circuit depth for quantum state preparation and general unitary synthesis, IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst. 42, 3301 (2023). [44] K. Hornik, M. Stinchcombe, and H. White, Multilayer feedforward networks are universal approximators, Neural Netw. 2, 359 (1989). [45] A. Barenco, C. H. Bennett, R. Cleve, D. P. DiVincenzo, N. Margolus, P. Shor, T. Sleator, J. A. Smolin, and H. Weinfurter, Elementary gates for quantum computation, Phys. Rev. A 52, 3457 (1995). [46] See supplementary information for details about the QKM algorithm, theoretical foundations, method comparisons, device information, benchmarks, ablations and error analysis. [47] BAQIS, Quafu superconducting quantum computing, https://quafu-sqc.baqis.ac.cn (2024). [48] P. Gray and S. Scott, Autocatalytic reactions in the isothermal, continuous stirred tank reactor: isolas and other forms of multistability, Chem. Eng. Sci. 38, 29 (1983). [49] J. Galewsky, R. K. Scott, and L. M. Polvani, An initial-value problem for testing numerical models of the global shallowwater equations, Tellus Ser. A-Dyn. Meteorol. Oceanogr. 56, 429 (2004). [50] E.U. Copernicus Marine Service (CMEMS), Global ocean gridded L4 sea surface heights and derived variables reprocessed 1993 ongoing, https://doi.org/10.48670/moi-00148 (2024). [51] N. Gourianov, M. Lubasch, S. Dolgov, Q. Y. van den Berg, H. Babaee, P. Givi, M. Kiffner, and D. Jaksch, A quantuminspired approach to exploit turbulence structures, Nat. Comput. Sci. 2, 30 (2022). [52] Z. Li, W. Han, Y. Zhang, Q. Fu, J. Li, L. Qin, R. Dong, H. Sun, Y. Deng, and L. Yang, Learning spatiotemporal dynamics with a pretrained generative model, Nat. Mach. Intell. 6, 1566 (2024). [53] S. L. Brunton, M. Budišić, E. Kaiser, and J. N. Kutz, Modern Koopman theory for dynamical systems, SIAM Rev. 64, 229 (2022). [54] M. Wang, X. Xue, M. Gao, and P. V. Coveney, Quantuminformed machine learning for predicting spatiotemporal chaos with practical quantum advantage, Sci. Adv. 12, eaec5049 (2026). [55] D. Jennings, K. Korzekwa, M. Lostaglio, and G. Wang, Quantum Koopman Algorithms, arXiv preprint arXiv:2605.19054 (2026). [56] Y. T. Lin, Y. Tian, D. Livescu, and M. Anghel, Data-driven learning for the Mori-Zwanzig formalism: a generalization of the Koopman learning framework, SIAM J. Appl. Dyn. Syst. 20, 2558 (2021). [57] The code is available at github.com/YYgroup/QKM.

Supplementary Information for “Quantum simulation of real-world nonlinear dynamics via Koopman method” Baoyang Zhang1 , Dong An2 , Zhaoyuan Meng3 , Yefei Yu4 , Xiaoxiao Xiao4 , Zhen Lu1,∗ and Yue Yang1,5,† 1 State Key Laboratory for Turbulence and Complex Systems, School of Mechanics and Engineering Science, Peking University, Beijing 100871, China 2 Beijing International Center for Mathematical Research, Peking University, Beijing 100871, China 3 Institute of Mechanics, State Key Laboratory of Nonlinear Mechanics, Chinese Academy of Sciences, Beijing 100190, China 4 Beijing Academy of Quantum Information Sciences, Beijing 100193, China 5 HEDPS-CAPT, Peking University, Beijing 100871, China ∗ [email protected]

† [email protected]

CONTENTS

1. Theorem proofs A. Proof of Theorem 1 B. Proof of Theorem 2 C. Proof of Theorem 3

S2 S2 S2 S4

2. QKM algorithm

3. Comparison with existing methods

4. State preparation

5. Autoencoder of QKM

6. Device information

S10

7. Description of benchmarks A. 3D reaction-diffusion systems B. Spherical fluid dynamics C. Real-world ocean currents

S12 S12 S12 S13

8. Detailed error analysis A. Theoretical estimation of measurement error B. Quantitative analysis of state preparation error C. Quantitative analysis of theoretical and optimization errors

S13 S14 S15 S16

9. Ideal noiseless simulation A. 3D reaction-diffusion systems B. Spherical fluid dynamics C. Real-world ocean currents

S17 S17 S17 S17

10. Ablation study A. About Theorem 2 B. About Theorem 3

S18 S18 S19

References

S21

S2 1.

THEOREM PROOFS

Proof of Theorem 1

˜ with We begin by recalling the LCHS theorem [S1] for time-independent matrices. For a general matrix 𝐴˜ = 𝐿˜ + i𝐻, 𝐿˜ = ( 𝐴˜ + 𝐴˜ † )/2, 𝐻˜ = ( 𝐴˜ − 𝐴˜ † )/2i and 𝐿˜ ⪰ 0, the non-unitary time evolution under − 𝐴˜ is given by ∫ 1 ˜ ˜ ˜ − 𝐴𝑡 e = e−i( 𝐻+𝑘 𝐿)𝑡 𝑑𝑘. (S1) 2 R 𝜋(1 + 𝑘 ) To map the original LCHS formulation to our context, we define 𝐴˜ := −𝐴, where 𝐴 = 𝐿 + i𝐻 with Hermitian matrices 𝐿 and 𝐻 subject to the stability condition 𝐿 ⪯ 0. It follows that the corresponding components are 𝐿˜ = −𝐿 and 𝐻˜ = −𝐻. The condition 𝐿 ⪯ 0 trivially guarantees that 𝐿˜ ⪰ 0, thereby satisfying the requirement of the original theorem. Substituting these components into Eq. (S1) yields ∫ ∫ 1 1 ˜ −i(−𝐻 −𝑘 𝐿)𝑡 e 𝑑𝑘 = ei(𝐻+𝑘 𝐿)𝑡 𝑑𝑘. (S2) e 𝐴𝑡 = e− 𝐴𝑡 = 2 2 R 𝜋(1 + 𝑘 ) R 𝜋(1 + 𝑘 ) Next, we define the parameterized matrices 𝐺 (𝑘) := 𝐻 + 𝑘 𝐿. Since both 𝐻 and 𝐿 are Hermitian, 𝐺 (𝑘) remains strictly Hermitian for all 𝑘 ∈ R. Consequently, by the spectral theorem, 𝐺 (𝑘) admits a unitary diagonalization 𝐺 (𝑘) = 𝑉 (𝑘)Λ(𝑘)𝑉 † (𝑘),

(S3)

where Λ(𝑘) is a real diagonal matrix containing the eigenvalues of 𝐺 (𝑘), and 𝑉 (𝑘) is the unitary matrix composed of its corresponding eigenvectors. Using the standard property of matrix exponentials for diagonalizable matrices, the unitary evolution operator associated with 𝐺 (𝑘) can be expanded as ei𝐺 (𝑘 )𝑡 = 𝑉 (𝑘)eiΛ(𝑘 )𝑡 𝑉 † (𝑘).

(S4)

Substituting Eq. (S4) into Eq. (S2) directly yields the desired result ∫ 1 e 𝐴𝑡 = 𝑉 (𝑘)eiΛ(𝑘 )𝑡 𝑉 † (𝑘)𝑑𝑘. 2 R 𝜋(1 + 𝑘 )

(S5)

The condition 𝐿 ⪯ 0 is satisfied without loss of generality. For any 𝐴, the substitution 𝑢 (𝑡) = e𝑏𝑡 𝑐 (𝑡) with the shift parameter 𝑏 ⩾ 𝜆max (𝐿) yields a modified system whose Hermitian part 𝐿 − 𝑏𝐼 is negative semi-definite.

Proof of Theorem 2

The continuous diagonal LCHS formulation expresses the evolution operator e 𝐴𝑡 as ∫ 𝐴𝑡 𝑈 (𝑡) = e = 𝑝(𝑘)𝑈 𝑘 (𝑡)𝑑𝑘,

(S6)

R 1 where 𝑝(𝑘) = 𝜋 (1+𝑘 2 ) is the density of the Cauchy-Lorentz distribution with

∫

𝑝(𝑘)𝑑𝑘 = 1, and 𝑈 𝑘 (𝑡) = 𝑉 (𝑘)eiΛ(𝑘 )𝑡 𝑉 † (𝑘) =

ei(𝐻+𝑘 𝐿)𝑡 denotes parameterized unitary operators. In practice, the integral in Eq. (S6) is approximated by a discrete quadrature over ℎ points. We first truncate the infinite integration domain to a symmetric interval [−𝐾, 𝐾]. Assuming a uniform sampling strategy with grid spacing Δ𝑘 = 2𝐾/ℎ, the total approximation error 𝜖spec is bounded by the sum of the truncation error 𝜖 trunc and the quadrature error 𝜖 quad . By the triangle

S3 inequality, we have ℎ ∑︁

∫ 𝜖spec :=

𝑝(𝑘)𝑈 𝑘 (𝑡)𝑑𝑘 − R

∫

∫ 𝐾

𝑝(𝑘)𝑈 𝑘 (𝑡)𝑑𝑘 −

𝑝(𝑘)𝑈 𝑘 (𝑡)𝑑𝑘 +

−𝐾

∫

∫ 𝐾 𝑝(𝑘)𝑈 𝑘 (𝑡)𝑑𝑘 −

⩽

Δ𝑘 · 𝑝(𝑘 𝑗 )𝑈 𝑘 𝑗 (𝑡)

𝑗=1

𝑝(𝑘)𝑈 𝑘 (𝑡)𝑑𝑘 − −𝐾

𝑝(𝑘)𝑈 𝑘 (𝑡)𝑑𝑘 −

𝑝(𝑘)𝑈 𝑘 (𝑡)𝑑𝑘

−𝐾

∫ 𝐾

∫ | 𝑘 |>𝐾

𝑝(𝑘)𝑈 𝑘 (𝑡)𝑑𝑘 −

𝑝(𝑘)𝑈 𝑘 (𝑡)𝑑𝑘

−𝐾

} |

𝜖trunc

Δ𝑘 · 𝑝(𝑘 𝑗 )𝑈 𝑘 𝑗 (𝑡)

𝑗=1

∫ 𝐾

−𝐾

ℎ ∑︁

2 ℎ ∑︁

Δ𝑘 · 𝑝(𝑘 𝑗 )𝑈 𝑘 𝑗 (𝑡)

𝑗=1

Δ𝑘 · 𝑝(𝑘 𝑗 )𝑈 𝑘 𝑗 (𝑡) ,

𝑗=1

}

(S7)

𝜖quad

where 𝑘 𝑗 = −𝐾 + 𝑗Δ𝑘 denotes the 𝑗-th quadrature point and 𝑘 0 = −𝐾. Given that 𝑈 𝑘 (𝑡) is a unitary operator, its spectral norm is exactly ∥𝑈 𝑘 (𝑡) ∥ 2 = 1. Exploiting the standard inequality arctan(𝑥) ⩽ 𝑥 for 𝑥 ⩾ 0, the truncation error 𝜖trunc is bounded as follows ∫ 𝜖trunc ⩽ 𝑝(𝑘) ∥𝑈 𝑘 (𝑡) ∥ 2 𝑑𝑘 | 𝑘 |>𝐾 2 1 = arctan 𝜋 𝐾 2 ⩽ . (S8) 𝜋𝐾 Define 𝑓 (𝑘) = 𝑝(𝑘)𝑈 𝑘 (𝑡). The quadrature error scales as 𝜖 quad =

ℎ ∫ 𝑘𝑗 ∑︁ 𝑗=1

⩽

𝑘 𝑗 −1

ℎ ∫ 𝑘𝑗 ∑︁

𝑓 (𝑘) − 𝑓 (𝑘 𝑗 ) 2 𝑑𝑘

𝑘 𝑗 −1

𝑗=1

⩽

𝑓 (𝑘) − 𝑓 (𝑘 𝑗 ) 𝑑𝑘

ℎ ∫ 𝑘𝑗 ∑︁

|𝑘 − 𝑘 𝑗 |𝑑𝑘 ·

𝑗=1

max

𝑘 ∈ [−𝐾 ,𝐾 ]

𝑘 𝑗 −1

𝑓 ′ (𝑘) 2

2𝐾 2 max 𝑓 ′ (𝑘) 2 ℎ 𝑘 ∈ [ −𝐾 ,𝐾 ] " # 2𝐾 2 𝜕𝑈 𝑘 (𝑡) ′ ⩽ max | 𝑝 (𝑘)| + | 𝑝(𝑘)| ℎ 𝑘 ∈ [ −𝐾 ,𝐾 ] 𝜕𝑘 2 =

⩽

2𝐾 2 2𝐾 2 𝜕𝑈 𝑘 (𝑡) max | 𝑝 ′ (𝑘)| + max | 𝑝(𝑘)| · max . ℎ 𝑘 ∈ [ −𝐾 ,𝐾 ] ℎ 𝑘 ∈ [−𝐾 ,𝐾 ] 𝜕𝑘 𝑘 ∈ [−𝐾 ,𝐾 ] 2

(S9)

Note that 1 max | 𝑝(𝑘)| = 𝜋 𝑘 ∈ [ −𝐾 ,𝐾 ]

and

√ √ 3 3 3 max | 𝑝 (𝑘)| ⩽ | 𝑝 ( )| = . 3 8𝜋 𝑘 ∈ [−𝐾 ,𝐾 ] ′

′

(S10)

⩽ 𝑡 ∥𝐿∥ 2 .

(S11)

By Duhamel’s Formula [S2], we have 𝜕𝑈 𝑘 (𝑡) = 𝑖 𝜕𝑘 2

∫ 𝑡 0

ei𝐺 (𝑘 ) (𝑡 − 𝜏 ) 𝐿ei𝐺 (𝑘 ) 𝜏 𝑑𝜏 2

S4 Substituting Eq. (S10) and Eq. (S11) into Eq. (S9), we obtain ! √ 2𝐾 2 3 3 𝑡 ∥𝐿 ∥ 2 𝜖quad ⩽ . + ℎ 8𝜋 𝜋

(S12)

Substituting Eq. (S8) and Eq. (S12) into Eq. (S7) yields ! √ 2 2𝐾 2 3 3 𝑡 ∥𝐿 ∥ 2 𝜖 spec ⩽ . + + 𝜋𝐾 ℎ 8𝜋 𝜋

(S13)

By optimizing with the geometric mean inequality, the overall error is bounded by ! 1/3 √ 3 3 3 + 8𝑡 ∥𝐿 ∥ 2 𝜖spec ⩽ , 𝜋 4ℎ √ where the optimal truncation radius 𝐾 =

3 3+8𝑡 ∥ 𝐿 ∥ 2 4ℎ

(S14)

−1/3 minimizes this upper bound.

Proof of Theorem 3

Consider the functional space on the 𝑛-dimensional Boolean hypercube, H = { 𝑓 | 𝑓 : {0, 1} 𝑛 → R}. For any 𝑓 , 𝑔 ∈ H , the normalized inner product is defined as the expectation over all inputs, ⟨ 𝑓 , 𝑔⟩ =

1 2𝑛

∑︁

𝑓 (𝑥)𝑔(𝑥).

(S15)

𝑥 ∈ {0,1} 𝑛

This functional space H is spanned by the orthonormal Walsh basis {𝑤 𝑠 } 𝑠∈ {0,1} 𝑛 [S3], with each element 𝑤 𝑠 (𝑥) = (−1) 𝑠·𝑥 where 𝑠 · 𝑥 =

P𝑛

𝑘=1 𝑠 𝑘 𝑥 𝑘

for

𝑥 ∈ {0, 1} 𝑛

(S16)

(mod 2) represents the bitwise inner product of 𝑠 and 𝑥. ∗

Because the action of a diagonal unitary 𝑈 ∗ = ei𝐻 on any computational basis |𝑥⟩ = |𝑥 1 ⟩ ⊗ · · · ⊗ |𝑥 𝑛 ⟩ yields a specific phase distribution 𝑈 ∗ |𝑥⟩ = ei 𝑓 ( 𝑥 ) |𝑥⟩, any diagonal Hamiltonian 𝐻 ∗ can be rigorously identified as an element of the space H , establishing a mapping between the operator’s eigenvalues and the Boolean scalar field 𝐻 ∗ |𝑥⟩ = 𝑓 (𝑥)|𝑥⟩ ∑︁ = 𝛼𝑠 𝑤 𝑠 (𝑥)|𝑥⟩ 𝑠∈ {0,1} 𝑛

∑︁

𝛼𝑠 (−1)

P𝑛

𝑘=1 𝑠𝑘 𝑥 𝑘

|𝑥 1 ⟩ ⊗ |𝑥 2 ⟩ ⊗ · · · ⊗ |𝑥 𝑛 ⟩

𝑠∈ {0,1} 𝑛

∑︁

𝛼𝑠 (−1) 𝑠1 𝑥1 |𝑥 1 ⟩ ⊗ (−1) 𝑠2 𝑥2 |𝑥2 ⟩ ⊗ · · · ⊗ (−1) 𝑠𝑛 𝑥𝑛 |𝑥 𝑛 ⟩

𝑠∈ {0,1} 𝑛

∑︁

𝛼𝑠 𝜎𝑧𝑠1 |𝑥1 ⟩ ⊗ 𝜎𝑧𝑠2 |𝑥 2 ⟩ · · · ⊗ 𝜎𝑧𝑠𝑛 |𝑥 𝑛 ⟩

𝑠∈ {0,1} 𝑛

∑︁

𝛼𝑠 𝑃𝑠 |𝑥⟩.

(S17)

𝑠∈ {0,1} 𝑛

Ë𝑛 𝑠𝑘 In this decomposition, 𝑃𝑠 = 𝑘=1 𝜎𝑧 denote the multi-qubit Pauli-𝑍 operators, and the Pauli coefficients 𝛼𝑠 are precisely the Walsh-Fourier transform of the eigenvalue spectrum 𝑓 (𝑥), given by 𝛼𝑠 = ⟨ 𝑓 , 𝑤 𝑠 ⟩ = 2−𝑛 Tr(𝐻 ∗ 𝑃𝑠 ). From Eq. (S17), we obtain an isomorphism between the functional Walsh basis and the multi-qubit Pauli-𝑍 operators via the relation 𝑃𝑠 |𝑥⟩ = 𝑤 𝑠 (𝑥)|𝑥⟩

for

𝑠, 𝑥 ∈ {0, 1} 𝑛 .

(S18)

S5 Since Eq. (S18) holds for all computational basis, the target diagonal Hamiltonian ∑︁ 𝐻∗ = 𝛼𝑠 𝑃 𝑠

(S19)

𝑠∈ {0,1} 𝑛

possesses a unique spectral decomposition, where the Hamming weight |𝑠| specifies the |𝑠|-body interaction order, representing the number of qubits actively participating in each multi-body Pauli operator 𝑃𝑠 . We now establish the equivalent representation of the subspace F for an 𝑛-qubit single-layer 𝑅 𝑧 circuit. Using the mutual commutativity of the Pauli-𝑍 operators, we have ei𝜙

  𝑛 ∑︁  𝜃 𝑗 ⊗ ( 𝑗 −1) ⊗ (𝑛− 𝑗 )  ⊗ 𝜎𝑧 ⊗ 𝐼1 𝑅 𝑧 (𝜃 𝑗 ) = exp i𝜙𝐼𝑛 − i 𝐼1 , 2   𝑗=1 𝑗=1  

𝑛 Ì

(S20)

where 𝐼1 denotes the 2 × 2 identity matrix and 𝐼𝑛 denotes the 2𝑛 × 2𝑛 identity matrix. By identifying 𝛽0 = 𝜙 and 𝛽𝑠 = −𝜃 𝑗 /2 for strings with Hamming weight |𝑠| = 1 in Eq. (S20), the ansatz subspace can be equivalently expressed as ∑︁ F = {ei𝐻 | 𝐻 = 𝛽0 𝐼𝑛 + 𝛽𝑠 𝑃𝑠 , 𝛽𝑠 ∈ R, 𝑠 ∈ {0, 1} 𝑛 }. (S21) |𝑠 |=1

This equivalence demonstrates that F is the manifold of unitaries whose generators are restricted to the 0-local and 1-local spectral components of the functional space H . Given by the normalized Frobenius norm as the distance metric, the approximation error for any 𝑈 ∈ F is h † i 1 1 ∥𝑈 − 𝑈 ∗ ∥ 2𝐹 = 𝑛 Tr 𝑈 − 𝑈 ∗ 𝑈 − 𝑈 ∗ 𝑛 2 2

= Tr 𝑈 †𝑈 + 𝑈 ∗†𝑈 ∗ − 𝑈 †𝑈 ∗ − 𝑈 ∗†𝑈 1 = 𝑛 2𝑛 + 2𝑛 − 2Re 𝑈 †𝑈 ∗ 2 1 = 2 − 𝑛−1 ReTr 𝑈 †𝑈 ∗ . 2

(S22)

Since all 𝑃𝑠 are diagonal matrices in the computational basis, they form a mutually commuting set [𝑃𝑠 , 𝑃𝑠′ ] = 0. From Eqs. (S19) and (S21), this commutativity allows the unitary overlap to be expressed as ∑︁ © ª © ∑︁ ª 𝑈 †𝑈 ∗ = exp −i𝛽0 𝐼𝑛 − i𝛽𝑠 𝑃𝑠 ® · exp i𝛼𝑠 𝑃𝑠 ® |𝑠 |=1 « ¬ «𝑠∈ {0,1} 𝑛 ¬   ∑︁ ∑︁   = exp i 𝛼0 − 𝛽0 𝐼𝑛 + i 𝛼𝑠 − 𝛽 𝑠 𝑃 𝑠 + i𝛼𝑠 𝑃𝑠  .   |𝑠 |=1 |𝑠 |⩾2  

(S23)

Substituting Eq. (S23) into Eq. (S22), and utilizing Eq. (S18), we obtain   ∑︁ ∑︁    cos  (𝛼0 − 𝛽0 ) + (𝛼𝑠 − 𝛽𝑠 )𝑤 𝑠 (𝑥) + 𝛼𝑠 𝑤 𝑠 (𝑥)    𝑥 ∈ {0,1} 𝑛 |𝑠 |=1 |𝑠 |⩾2          ∑︁  ∑︁ ∑︁     1   1 − cos  (𝛼0 − 𝛽0 ) + (𝛼𝑠 − 𝛽𝑠 )𝑤 𝑠 (𝑥) + 𝛼𝑠 𝑤 𝑠 (𝑥)  = 𝑛−1  2     𝑥 ∈ {0,1} 𝑛  |𝑠 |=1 |𝑠 |⩾2   " #  P P (𝛼0 − 𝛽0 ) + |𝑠 |=1 (𝛼𝑠 − 𝛽𝑠 )𝑤 𝑠 (𝑥) + |𝑠 |⩾2 𝛼𝑠 𝑤 𝑠 (𝑥) 4 ∑︁ = 𝑛 sin2 . 2 2 𝑛

1 1 ∥𝑈 − 𝑈 ∗ ∥ 2𝐹 = 2 − 𝑛−1 𝑛 2 2

∑︁

(S24)

𝑥 ∈ {0,1}

To establish the lower bound, we apply the Taylor expansion sin2 (𝜃) = 𝜃 2 + O (𝜃 4 ) and the orthonormality of the Walsh

S6 functions ⟨𝑤 𝑠 , 𝑤 𝑠′ ⟩ = 𝛿 𝑠𝑠′ (with the Kronecker delta function 𝛿 𝑠𝑠′ ) to Eq. (S24), yielding 2  ∑︁ ∑︁    + O 𝜁4  (𝛼0 − 𝛽0 ) + (𝛼 − 𝛽 )𝑤 (𝑥) + 𝛼 𝑤 (𝑥) 𝑠 𝑠 𝑠 𝑠 𝑠    |𝑠 |=1 |𝑠 |⩾2 𝑥 ∈ {0,1} 𝑛   ∑︁ ∑︁ 2 2 2 4 = (𝛼0 − 𝛽0 ) + (𝛼𝑠 − 𝛽𝑠 ) + 𝛼𝑠 + O 𝜁

1 1 ∥𝑈 − 𝑈 ∗ ∥ 2𝐹 = 𝑛 𝑛 2 2

∑︁

|𝑠 |=1

(S25)

|𝑠 |⩾2

P P with 𝜁 = |𝛼0 − 𝛽0 | + |𝑠 |=1 |𝛼𝑠 − 𝛽𝑠 | + |𝑠 |⩾2 |𝛼𝑠 |. To minimize Eq. (S25), the variational parameters must be chosen such that the 0-local and 1-local residuals vanish, i.e., 𝛽0 = 𝛼0 and 𝛽𝑠 = 𝛼𝑠 for all |𝑠| = 1. Consequently, the minimum approximation error is given by ∑︁ 1 ∗ 2 2 4 min ∥𝑈 − 𝑈 ∥ 𝐹 = 𝛼𝑠 + O max 𝛼𝑠 . (S26) 𝑈 ∈ F 2𝑛 |𝑠 |⩾2 |𝑠 |⩾2

This confirms that the single-layer 𝑅 𝑧 ansatz is intrinsically limited by the high-order Walsh components of the target Hamiltonian.

QKM ALGORITHM

Training the QKM requires joint optimization of the autoencoder and Koopman operator parameters [S4]. The training dataset comprises multiple time-series trajectories of nonlinear dynamics, each containing 𝑇 + 1 consecutive states {𝑥 𝑘 }𝑇𝑘=0 sampled at intervals Δ𝑡, where 𝑥 𝑘 = 𝑥 (𝑘Δ𝑡). We implement the supervised learning with data organized as pairs (𝑥 𝑘 , Δ𝑘), in which Δ𝑘 ∈ [0, 𝑇] specifies the number of time steps to be predicted ahead. Given such a pair, the QKM maps 𝑥 𝑘 through the encoder E, evolves the resulting quantum state for a duration Δ𝑘Δ𝑡 via the parameterized unitary U, and reconstructs the prediction through the decoder D. The composite loss function minimizes the normalized prediction error over all training pairs,  2  D ◦ U Δ𝑘Δ𝑡 ◦ E (𝑥 𝑘 ) − 𝑥 𝑘+Δ𝑘 2  . L = E ( 𝑥𝑘 ,Δ𝑘 )   ∥𝑥 𝑘+Δ𝑘 ∥ 22    

(S27)

The expectation in Eq. (S27) is taken over the index set I consists of two complementary subsets I = {(ℓ, 0) | ℓ ∈ {0, . . . , 𝑇 }} ∪ {(0, ℓ) | ℓ ∈ {1, . . . , 𝑇 }} . | {z } | {z } Zero-step transitions

(S28)

Initial rollouts

The zero-step enforces autoencoder reconstruction fidelity for Δ𝑘 = 0, while the initial rollouts enforce the predictive accuracy of the learned Koopman dynamics from the initial condition 𝑥 (0). All variational parameters, including the NN encoder–decoder weights and the quantum circuit rotation angles, are optimized jointly on NVIDIA H100 GPUs. Detailed architectural specifications of the autoencoder are provided in Sec. 5.

S7 Algorithm 1 Training of the QKM Require: Dataset {𝑥 𝑘 , Δ𝑘, 𝑥 𝑘+Δ𝑘 }, state preparation circuit structure 1: Set circuit configuration 𝑛, ℎ, 𝑅, 𝑟 and training hyperparameters 2: Initialize classical encoder E 𝜃 , decoder D 𝜙 , and quantum evolution parameters Θ 3: for each training epoch do 4: for each mini-batch (𝑥 𝑘 , Δ𝑘, 𝑥 𝑘+Δ𝑘 ) do 5: Encode the input field: {𝑃 𝑗 } ℎ𝑗=1 ← E 𝜃 (𝑥 𝑘 ) 6:

Prepare the initial quantum state: 𝑗

|𝑢 0 ⟩ ← 𝑃 𝑗 7:

for

Evolve through quantum circuits: 𝑗

𝑗

|𝑢 Δ𝑘 ⟩ ← U (|𝑢 0 ⟩, Θ, Δ𝑘) 8:

for

𝑗 = 1, . . . , ℎ

Compute the probability distribution: 𝑗

𝑗

𝑢 Δ𝑘 ← |𝑢 Δ𝑘 ⟩ 9:

𝑗 = 1, . . . , ℎ

for

𝑗 = 1, . . . , ℎ

Decode to obtain the output field: ℎ ) 𝑥ˆ 𝑘+Δ𝑘 ← D 𝜙 (𝑢 1Δ𝑘 , . . . , 𝑢 Δ𝑘

10:

Compute loss: ∥ 𝑥ˆ 𝑘+Δ𝑘 − 𝑥 𝑘+Δ𝑘 ∥ 2 ∥𝑥 𝑘+Δ𝑘 ∥ 2

L←

11: Backpropagate and update 𝜃, 𝜙, Θ via AdamW optimizer 12: end for 13: end for 14: return Trained model E 𝜃 , D 𝜙 and quantum evolution parameters Θ

Algorithm 2 End-to-end implementation of the QKM Require: Physical initial condition 𝑥(0), target evolution time 𝑡, Measurement shot count 𝑀 Require: Trained model E 𝜃 , D 𝜙 and parameters Θ for circuits 1: Classical pre-processing: {𝑃 𝑗 } ℎ𝑗=1 ← E 𝜃 𝑥 (0) 2: Run circuits on superconducting quantum processors: 𝑗

|𝑢 𝑡 ⟩ ← Circuit(𝑃 𝑗 , Θ, 𝑡)

for

𝑗 = 1, . . . , ℎ

3: Perform projective measurements in the Z-basis for 𝑀 shots: 𝑗

𝑗

𝑢ˆ 𝑡 ← |𝑢 𝑡 ⟩

for

𝑗 = 1, . . . , ℎ

4: Classical post-processing:

𝑥(𝑡) ˆ ← D 𝜙 ( 𝑢ˆ 1𝑡 , . . . , 𝑢ˆ 𝑡ℎ ) 5: return Target-time physical fields 𝑥(𝑡) ˆ

Algorithm 1 outlines the training process used to learn the unitary Koopman operators and autoencoder parameters from

S8 nonlinear trajectories. Algorithm 2 specifies the hardware deployment pipeline, describing the parallel execution of the ℎ discretized spectral components and the subsequent reconstruction of the evolved physical state. 3.

COMPARISON WITH EXISTING METHODS

To situate the QKM within the current landscape of quantum algorithms for nonlinear dynamics, we provide a comparative analysis with two leading paradigms: the Carleman method [S9, S10] and the Koopman-von Neumann (KvN) method [S11, S12], as summarized in Table S1. TABLE S1. Comparison of quantum algorithms for nonlinear classical dynamics. The table contrasts the proposed QKM with the established Carleman method and KvN method across theoretical foundations, algorithmic primitives, and hardware requirements.

Aspect Theoretical basis Nonlinearity requirement Quantum primitive State preparation Hardware requirement Drawback

Carleman method Carleman linearization Weak QLSP Oracle required Fault-tolerant Limited convergence radius

KvN method Liouville equation Moderate (phase space) Hamiltonian simulation Oracle required Fault-tolerant Curse of dimensionality

Our QKM Koopman theory Moderate (physical space) Hamiltonian simulation End-to-end encoding NISQ-era Data-driven optimization required

Existing frameworks predominantly rely on strict analytical expansion or rigorous phase-space mapping. The Carleman method employs a Taylor-like truncation to embed nonlinear ordinary differential equations (ODEs) into an infinite-dimensional linear system. While mathematically rigorous and capable of exponential quantum speedup, its convergence fundamentally restricts the method to weakly nonlinear regimes. Although recent advances accommodate broader stable systems [S10], the limited convergence radius remains a primary drawback. The KvN method maps the classical Liouville equation to a unitary Schrödinger-like evolution. However, this continuous mapping requires discretizing the phase space of all possible system states, incurring an exponential overhead due to the curse of dimensionality. The QKM departs from these rigid analytical constraints by adopting a data-driven strategy based on Koopman theory. By jointly training the Koopman propagator and the circuit ansatz, the QKM dynamically identifies an invariant subspace to effectively capture moderately nonlinear dynamics in physical space, such as reaction-diffusion systems and fluid flows, that remain beyond the efficient reach of fixed-truncation techniques, albeit at the cost of requiring data-driven optimization. Furthermore, the algorithmic primitives and hardware prerequisites differ fundamentally across these paradigms (Table S1). Algorithms based on Carleman embeddings typically rely on quantum linear system problem (QLSP) solvers, while advanced implementations of the KvN method necessitate LCHS or quantum singular value transformation [S12]. Both pathways demand the block-encoding of non-unitary operators, requiring deep quantum circuits, substantial ancilla overhead, and complex statepreparation oracles that are exclusively viable on fault-tolerant quantum computers. Conversely, the QKM leverages Hamiltonian simulation primitives to circumvent these prohibitive overheads, compiling the propagator directly into an ensemble of shallow, topology-native circuits. By optimizing the state encoding within the basis transformation, the QKM eliminates the need for complex block-encoding and oracle queries. This yields an end-to-end simulation pipeline that is adaptive for NISQ and early fault-tolerant hardware. The capability of the QKM to reach moderately nonlinear dynamics rests on an alignment between the locality structure of physical Koopman generators and the native cost structure of quantum circuits. A general Koopman generator on a 2𝑛 dimensional observable space would require all 2𝑛 Pauli-𝑍 strings and thus exponentially many parameters. In contrast, many physical systems concentrate on low-order interactions among small subsets of observables. Theorem 3 establishes that the single-layer 𝑅 𝑧 ansatz represents zeroth- and first-order Pauli-𝑍 terms exactly, so the captured modes coincide with those favored by physics. The NN encoder reinforces this concentration by learning an observable basis that further shifts spectral weight onto the captured modes. The QKM-amenable regime therefore consists of systems whose Koopman generators are few-body in the learned observable basis, allowing an exponentially large observable space to be propagated with polynomial quantum resources. Systems with stronger high-order couplings can be addressed by adding multi-qubit gates, at the cost of deeper circuits and reduced speedup. 4.

STATE PREPARATION

We justify the efficient state-preparation assumption. In the hardware implementation, the state-preparation circuit has depth O (𝑛) under the scaling 𝑅𝑟 = O (𝑛). Here we provide a complementary oracle-based construction showing that, for structured physical initial conditions, preparing the corresponding state does not require exponential resources.

S9 The initial state of the linearized system is |𝑢 0 ⟩ = 𝑢 0 /∥𝑢 0 ∥ where 𝑢 0 = [𝑔1 (𝑥(0)); · · · , 𝑔 𝑁 (𝑥(0))]. We start with |0⟩ and apply the Hadamard gate on each qubit to obtain the uniform superposition 𝑁 −1

1 ∑︁ | 𝑗⟩. √ 𝑁 𝑗=0

(S29)

To proceed, we append two ancilla registers initialized with all zeros, denoted as |0⟩𝑔 and |0⟩𝑟 . Here |0⟩𝑔 is used to binary encode 𝑔 𝑗 (𝑥(0)) and |0⟩𝑟 is a single qubit for rotations. We then apply the operator 𝑂 𝑔 : | 𝑗⟩|0⟩𝑔 ↦→ | 𝑗⟩|𝑔 𝑗 (𝑥(0))/𝑔max ⟩, where 𝑔max = max{|𝑔 𝑗 (𝑥(0))|}, |𝑔 𝑗 (𝑥(0))/𝑔max ⟩ represents a binary encoding of 𝑔 𝑗 (𝑥(0))/𝑔max , and the cost of constructing 𝑂 𝑔 will be discussed later. This gives rise to the state 𝑁 −1

1 ∑︁ | 𝑗⟩|𝑔 𝑗 (𝑥(0))/𝑔max ⟩𝑔 |0⟩𝑟 . √ 𝑁 𝑗=0

(S30)

√︁ Applying the controlled rotation c-R : |𝜃⟩𝑔 |0⟩𝑟 ↦→ |𝜃⟩𝑔 (𝜃|0⟩𝑟 + 1 − |𝜃| 2 |1⟩𝑟 ) yields 𝑁 −1

1 ∑︁ © 𝑔 𝑗 (𝑥(0)) | 𝑗⟩|𝑔 𝑗 (𝑥(0))/𝑔max ⟩𝑔 |0⟩𝑟 + √ 𝑔max 𝑁 𝑗=0 «

√︄ 1−

|𝑔 𝑗 (𝑥(0))| 2 2 𝑔max

ª |1⟩𝑟 ® . ¬

(S31)

Uncompute the register with subscript 𝑔 by applying 𝑂 †𝑔 , and we obtain 𝑁 −1

1 ∑︁ © 𝑔 𝑗 (𝑥(0)) | 𝑗⟩|0⟩𝑔 |0⟩𝑟 + √ 𝑔max 𝑁 𝑗=0 «

√︄ 1−

|𝑔 𝑗 (𝑥(0))| 2 2 𝑔max

𝑁 −1 ∥𝑢 0 ∥ ∑︁ ª |1⟩𝑟 ® = √ |𝑢 0 ⟩|0⟩𝑔 |0⟩𝑟 + | ⊥⟩, 𝑁𝑔max 𝑗=0 ¬

(S32)

where | ⊥⟩ represents a possibly unnormalized state with |1⟩𝑟 in this ancilla qubit. Therefore, the desired state is encoded after projecting both ancilla registers to zero, which can be achieved by measurements under the computational basis or more efficiently amplitude amplification. Later on we will focus on the amplitude amplification procedure as the final step. Now we estimate the overall gate complexity of the above approach. In each round of the method, we apply 𝑂 𝑔 and c-R for O (1) times. We assume that 𝑔 𝑗 (𝑥(0))’s are classically efficiently computable (as a function of 𝑗 and 𝑥(0) which can be classically encoded with O (𝑛 + 𝑑) bits under the constant machine precision), which means that the classical cost of computing 𝑔 𝑗 (𝑥(0)) is O (poly(𝑛𝑑)). Therefore, the quantum operation 𝑂 𝑔 can be implemented with O (poly(𝑛𝑑)) complexity following the standard reversible circuit construction [S3]. The operation c-R and its variant have been widely used in many existing quantum algorithms such as the Harrow–Hassidim–Lloyd algorithm [S16] and can also be efficiently implemented with constant complexity when the number of ancilla qubits in |0⟩𝑔 is constant, which is the case when the binary encoding precision is constant (but can be sufficiently small). Therefore, the number of gates required in each round of the approach is O (poly(𝑛𝑑)), so the overall complexity after amplitude amplification is ! √ 𝑁𝑔max O poly(𝑛𝑑) . (S33) ∥𝑢 0 ∥ √

√︃P 𝑁 −1 2 𝑁. However, the norm ∥𝑢 0 ∥ = 𝑗=0 |𝑔 𝑗 (𝑥(0))| on the √ denominator is the square root of 𝑁 many terms, which also scales as 𝑁 and leads to an overall O (poly(𝑛𝑑)) complexity, under the condition that a constant portion of 𝑔 𝑗 (𝑥(0))’s are on the same scale as 𝑔max . Notice that such a condition can be satisfied as long as the distribution of 𝑔 𝑗 (𝑥(0)) does not concentrate too much around its peak, in the cases such as 𝑔max /𝑔min = O (1) where 𝑔min = min 𝑗 {|𝑔 𝑗 (𝑥(0))|} or the median of 𝑔 𝑗 (𝑥(0)) is on the same scale as 𝑔max . Therefore, this theoretical construction supports the efficient state-preparation assumption in the main text, while the practical implementation further realizes this step through the NN-generated PQC with O (𝑛) depth. At a first glance, Eq. (S33) appears to be on a huge scale ∼

AUTOENCODER OF QKM

The QKM autoencoder architecture serves as the interface between the original state space and quantum representation, learning to encode input fields into quantum circuits and reconstruct evolved fields from the measurement of circuits. ,𝑅 The encoder E maps the input state 𝑥 (0) into ℎ sets of parameters 𝑃 𝑘 = {𝛼 𝑗,𝑙,ℓ , 𝛽 𝑗,𝑙,ℓ , 𝛾 𝑗,𝑙,ℓ } 𝑛,𝑟 𝑗=1,𝑙=1,ℓ=1 for 𝑘 = 1, . . . , ℎ,

S10 which define the rotation angles of the 𝑈3 layers that prepare the initial quantum state on ℎ quantum circuits. The encoding proceeds as 𝑦 1 = R 𝑥 (0) , n o 𝑦 2 = C M (𝑦 1 ) + R ◦ C 𝑥 (0) , n o 𝑦 3 = C M 𝑦 2 + T ◦ R ◦ C 𝑥 (0) , (S34) n o 𝜙1 (0) = C M 𝑦 3 + T ◦ T ◦ R ◦ C 𝑥 (0) , n o 𝜙2 (0) = C M 𝑦 4 + T ◦ T ◦ T ◦ R ◦ C 𝑥 (0) , where R (·), C(·), and T (·) denote the residual, convolution, and attention blocks, respectively; M (·) represents a convolution operator with a kernel size of three and stride of two for spatial downsampling and ◦ denotes function composition. The gate parameters {𝑃 𝑘 } ℎ𝑘=1 correspond directly to the hierarchical feature maps {𝜙1 (0) , 𝜙2 (0)}, where the dimensions are regulated through the convolutional channel capacity to align with the PQC configurations. The decoder D reconstructs the evolved state from the quantum-processed observables 𝜙1 (𝑡) and 𝜙2 (𝑡), mapping the measurement outcomes of the evolved quantum states back to the state space X through the learned inverse transformation n o 𝑥 (𝑡) = Q ◦ C 𝑤 𝑦 1 , P (𝑧2 ) , n o 𝑧2 = R ◦ C 𝑤 𝑦 2 , P (𝑧3 ) , n o 𝑧3 = R ◦ C 𝑤 𝑦 3 , P (𝑧4 ) , (S35) n o 𝑧4 = R ◦ C 𝑤 𝜙1 (𝑡) , P (𝑧5 ) , 𝑧5 = 𝜙2 (𝑡) , where P (·) denotes a 3 × 3 transposed convolution with stride two for upsampling, Q (·) is a 1 × 1 convolution with stride one for output refinement, and 𝑤(·, ·) handles tensor reshaping and channel-wise concatenation. The autoencoder is realized through a modified U-Net [S5] architecture enhanced with vision transformer components [S6], as illustrated in Fig. S1. The system employs three distinct computational blocks for encoding and their corresponding inverse operations for decoding, with skip connections facilitating information flow between corresponding encoder-decoder levels. It combines convolutional neural networks for spatial processing and transformer mechanisms for long-range dependencies, forming a framework for learning the global linearized observable space of dynamical systems. The residual blocks R (·) (red block in Fig. S1) integrate dual convolution (Conv)–batch normalization (BatchNorm)–sigmoid linear unit (SiLU) layers with skip connections to facilitate gradient flow and preserve fine-grained features. The convolution blocks C (·) (blue block in Fig. S1) perform spatial downsampling during encoding and feature processing during decoding. The attention blocks T (·) (green block in Fig. S1) utilize overlap embedding [S7] with 4 attention units to capture long-range spatial dependencies in the feature maps. Each attention unit utilizes multi-head self-attention through multi-layer perceptron (MLP)-based transformations to capture spatial relationships within the encoded features. The integration of cascaded Conv– SiLU–MLP operations further facilitates channel mixing and feature transformation to improve the representational capacity of the model. The decoder mirrors the encoder’s hierarchical structure through transposed convolutions P (·) (yellow block in Fig. S1) for progressive upsampling and physical field recovery. The integration of time-invariant quantities {𝑦 1 , 𝑦 2 , 𝑦 3 } and quantum evolution results {𝜙1 (𝑡), 𝜙2 (𝑡)} ensures the accurate reconstruction of the dynamical state. Concatenation across the latent space maintains structural integrity and enables seamless feature fusion throughout the autoencoder.

DEVICE INFORMATION

All quantum experiments were performed on the superconducting quantum processor “Yudu” [S8]. “Yudu” consists of frequency-tunable transmon qubits connected by tunable couplers, with a lattice connectivity that supports nearest-neighbour two-qubit operations. The native gate set used in our experiments is composed of single-qubit rotations and controlled-𝑍 gates, summarized in Table S2. The QKM circuits are mapped to hardware-native subgraphs of the processor, as shown in Fig. S2. For the 3D reaction– diffusion benchmark, each spectral component is implemented on a 6-qubit subcircuit. For the spherical fluid dynamics and real-world ocean-current benchmarks, each component is implemented on a 10-qubit subcircuit. Different spectral components

S11 Input

Encoder

Decoder

Output

Attention block

Residual block

SiLU

Multi-head self-attention

Conv-Norm

MLP

SiLU

Conv-Norm

h quantum circuits

Overlap embedding

Addition C

Concatenation Deconvolution

MLP Conv 3x3(x3)

1x1(x1) Conv kernel

SiLU

3x3(x3) Conv kernel

Reshape

SiLU

Conv-Norm

MLP

FIG. S1. Architecture of the QKM autoencoder. The encoder E transforms the physical initial condition 𝑥 (0) into ℎ parameter sets to define the 𝑈3 rotation angles for the initial quantum state |𝑢 (0)⟩ within each circuit. The ℎ parallel quantum circuits implement the unitary evolution |𝑢 (𝑡)⟩ = U 𝑡 |𝑢 (0)⟩ to represent the spectral components of the diagonalized LCHS. The classical decoder D reconstructs the physical state 𝑥 (𝑡) from the measurement outcomes of the evolved quantum states |𝑢 (𝑡)⟩. The insets detail the attention blocks, residual blocks, and the computational symbol legend. The attention block utilizes multi-head self-attention mechanisms and MLP-based transformations. The residual block architecture incorporates dual convolution and batch normalization layers with SiLU activation functions. The architectural symbol legend defines computational operations including convolution, deconvolution, batch normalization, activation functions, and tensor operations.

TABLE S2. Calibration data for the qubits and CZ couplings used in the hardware experiments on the superconducting quantum processor “Yudu”. Gate fidelities 𝐹1Q and 𝐹CZ are reported as percentages (%), while coherence times 𝑇1 and 𝑇2∗ are reported in 𝜇s. Qubits shared by different subcircuits are listed only once. CZ labels omit the leading “Q”.

Qubit properties

CZ couplings 𝑇2∗

Qubit

𝐹1Q

Q25 Q30 Q31 Q36 Q37

99.93 44.6 3.8 Q38 99.85 48.3 2.2 Q42 99.90 53.0 10.1 Q43 99.89 45.0 7.4 99.92 42.8 6.8

𝑇1

Qubit

𝐹1Q

𝑇1

𝑇2∗

Qubit

99.93 57.1 7.6 Q44 99.93 48.7 3.9 Q49 99.93 37.8 14.6 Q50

𝐹1Q

𝑇1

𝑇2∗

99.95 54.8 8.2 99.94 52.0 4.8 99.95 46.3 7.3

𝐹CZ

25–31 98.77 30–25 97.69 31–38 99.59 36–30 98.64 42–37 98.99

𝐹CZ

37–31 99.22 38–43 98.42 38–44 99.17 42–36 99.34

𝐹CZ

43–49 99.50 44–50 98.55 49–42 98.87 50–43 97.85

are assigned to disjoint qubit subsets, enabling parallel execution on the same processor. The subcircuit topology is chosen to match the ring-like entangling layout of the QKM ansatz, thereby reducing compilation overhead and avoiding unnecessary SWAP operations. We also evaluate the resource scaling before and after hardware compilation. The design ansatz is used for algorithmic complexity analysis, whereas the compiled circuits quantify the practical resources required after mapping to the “Yudu” topology. As shown in Figs. S2b,c, topology-native compilation reduces both circuit depth and gate count while preserving the intended circuit structure.

S12 c

Design ansatz Compiled circuit

150

∼

×103

∼

Circuit depth

Gate count

5 3/2

1/2

∼n

0 0

25 50 75 100 Number of qubits, n

FIG. S2. Quantum processor architecture, topology-native subcircuit mapping, and resource scaling. (a) Qubit layout and lattice connectivity of the superconducting quantum processor “Yudu”, together with representative 6- and 10-qubit subcircuit mappings used for the 3D reaction– diffusion benchmark and for the spherical fluid and ocean-current benchmarks, respectively. The subcircuits are assigned to disjoint qubit subsets to enable parallel hardware execution. (b,c) Scaling of circuit depth (b) and gate count (c) with the number of qubits 𝑛, comparing the design ansatz before compilation with the compiled circuits after hardware mapping. The design ansatz exhibits approximately linear depth scaling and quadratic gate-count scaling, while topology-native compilation reduces the implemented depth and gate count. All complexity analyses are based on design ansatz (pre-compilation) for fairness.

DESCRIPTION OF BENCHMARKS A.

3D reaction-diffusion systems

The Gray-Scott equation [S13] 𝜕𝑢 = 𝐷 𝑢 ∇2 𝑢 − 𝑢𝑣 2 + 𝐹 (1 − 𝑢), 𝜕𝑡 𝜕𝑣 = 𝐷 𝐵 ∇2 𝑣 + 𝑢𝑣 2 − (𝐹 + 𝐾)𝑣 𝜕𝑡

(S36) (S37)

provides a quintessential representation of complex pattern formation in reaction-diffusion systems. This system describes the interaction between two chemical species with concentrations 𝑢 and 𝑣, where the respective diffusion coefficients are 𝐷 𝑢 = 2 × 10−5 and 𝐷 𝑣 = 1 × 10−5 . The feed rate 𝐹 and kill rate 𝐾 collectively determine the morphological evolution of the system, governing the transition between stationary states and complex Turing patterns. The dataset was generated via direct numerical simulation (DNS) using the spectral method [S14], with parameters set to 𝐹 = 0.018 and 𝐾 = 0.051. The computational domain spans [−1, 1] 3 under periodic boundary conditions, discretized with 643 grid points. We generated a dataset comprising 1200 independent trajectories initiated from random Gaussian fields. Each trajectory consists of 61 temporal snapshots, corresponding to a final index 𝑇 = 60 as defined in Eq. (S28), with a sampling interval of Δ𝑡 = 5 covering the total evolution period 𝑡 ∈ [0, 300]. The datasets were partitioned into training, validation, and test sets with a ratio of 8:1:1. Model parameters were optimized using the training set, employing a cosine decay learning rate schedule with logarithmic warm-up, while the validation set guided early stopping. Final performance was assessed on the test set.

Spherical fluid dynamics

Spherical shallow-water equations are instrumental in modeling global atmospheric circulation and geophysical flows. We adopted the barotropic instability test case [S15], which involves a potent zonal jet subject to small perturbations that trigger complex vortex shedding and turbulent cascades. The governing equations 𝜕𝑢 = −(𝑢 · ∇)𝑢 − 𝜈∇4 𝑢 − 𝑔∇ℎ − 2𝛀 × 𝑢, 𝜕𝑡 𝜕ℎ = −∇ · (ℎ𝑢) − 𝜈∇4 ℎ − 𝐻∇ · 𝑢 𝜕𝑡

(S38) (S39)

are formulated in spherical coordinates to account for the Coriolis effect and the intrinsic curvature of the Earth. Physical parameters are scaled to Earth’s characteristics, where the unit length corresponds to the mean planetary radius 6.37 × 106 meters

S13 and the temporal unit is one hour. In this dimensionless representation, the planetary radius is 𝑅 = 1 and the magnitude of the rotation rate |𝛀| is 0.26. The angular velocity vector 𝛀 is oriented parallel to the planetary rotation axis, which determines the meridional variation of the Coriolis effect. The gravitational acceleration 𝑔 and the mean fluid depth 𝐻 are set to 19.95 and 1.57 × 10−3 respectively, while the hyperviscosity coefficient 𝜈 is approximately 8.66 × 10−9 to ensure stability. The system is discretized using a spherical harmonic basis with 512 × 256 grid points across the longitudinal and latitudinal dimensions. To mitigate singularities at the poles, the latitudinal domain is extended through the concatenation of the original grid with a version flipped meridionally and shifted by 180 degrees zonally. The resulting 512 × 512 grid representation ensures numerical continuity across the polar axes and preserves the physical integrity of the spherical manifold during feature extraction. The dataset was generated through DNS using a pseudo-spectral code [S16]. The initial velocity field is constructed as a mid-latitude jet with a random peak velocity. To ensure physical consistency, the initial height field is obtained by solving a linear boundary value problem to satisfy the geostrophic balance condition. This stable equilibrium is subsequently perturbed by a localized Gaussian height anomaly with a random amplitude to trigger the barotropic instability. We generated a dataset comprising 1200 independent trajectories. Each trajectory consists of 51 temporal snapshots, corresponding to a final index 𝑇 = 50 as defined in Eq. (S28), with a sampling interval of Δ𝑡 = 2 covering the fully developed turbulent regime 𝑡 ∈ [400, 500]. The resulting data were partitioned into training, validation, and test sets according to an 8:1:1 ratio. Model parameters were optimized using the training set with a cosine decay learning rate schedule preceded by a logarithmic warm-up phase. The validation set was employed to guide the early stopping criterion to prevent overfitting, while the final predictive performance was assessed on the independent test set.

Real-world ocean currents

Real-world ocean currents exhibit complex nonlinear dynamics, motivating the selection of the Gulf Stream as a representative benchmark to validate the performance of the QKM. The spatial domain covers latitudes from 20◦ N to 52◦ N and longitudes from 33◦ W to 65◦ W, discretized into 256 × 256 grid points with a spatial resolution of 0.125◦ . We utilize absolute geostrophic velocity fields derived from Level-4 satellite reanalysis products to represent the surface flow field [S17]. To incorporate geographical constraints and boundary conditions, these velocity fields are integrated with static environmental variables [S18], specifically the land-sea mask and the geopotential representing the regional topographic height. Spanning the period from 1993 to 2024, the dataset is reorganized into independent trajectories with a final index 𝑇 = 12 and a 13-day temporal duration. The resulting data are partitioned into disjoint sets for training, validation, and independent testing. Model optimization is conducted using the training sequences from 1993 to 2022 with a cosine decay learning rate schedule preceded by a logarithmic warm-up phase. The validation data from 2022 serve to guide the early stopping criterion to prevent overfitting, while the final predictive robustness is assessed on the independent test set from 2024.

DETAILED ERROR ANALYSIS

Throughout this work, 𝜀 denotes relative errors and 𝜖 denotes absolute errors. The cumulative error 𝜖 of the QKM is partitioned into theoretical and experimental constituents. 𝜖 = 𝜖th + 𝜖 exp , 𝜖 th = 𝜖proj + 𝜖 spec + 𝜖 ansatz , 𝜖 exp = 𝜖meas + 𝜖 prep + 𝜖 opt + 𝜖 noise .

(S40)

Each term corresponds to a stage of the simulation pipeline. In the theoretical error 𝜖th , the projection error 𝜖 proj arises from truncating the infinite-dimensional Koopman operator to a finite 𝑁-dimensional subspace, which depends on the expressivity of the learned observable functions. The spectral sampling error 𝜖spec originates from approximating the integral in Eq. (1) by a discrete sum over ℎ points. The ansatz approximation error 𝜖ansatz reflects the structural bias introduced by representing each diagonalized evolution operator with a single 𝑅 𝑧 layer. This design prioritizes circuit shallowness over universal expressivity, limiting the fidelity of the unitary synthesis to the range achievable by the chosen parametric form. In the experimental error 𝜖exp , the measurement error 𝜖 meas originates from the reconstruction of the state probabilities |⟨𝑖|𝑢(𝑡)⟩| 2 (𝑖 = 1, . . . , 𝑁) based on outcomes from the ℎ PQCs. The state preparation error 𝜖 prep stems from the {𝑈3 , CZ} PQCs used to encode the initial state |𝑢 (0)⟩. The optimization error 𝜖opt includes the convergence gap of the variational optimization and the generalization error. The hardware noise error 𝜖 noise encapsulates the inaccuracies introduced by gate infidelities, decoherence, and measurement errors during execution on physical quantum processors. While the total error 𝜖 and hardware noise 𝜖noise are empirically evaluated in the Results section and the bounds for 𝜖 spec and 𝜖ansatz are provided by Theorems 2 and 3, respectively, this section offers analysis of the remaining terms. Specifically, we present

S14 theoretical estimation for 𝜖meas , and quantitative discussion of 𝜖prep and 𝜖opt based on empirical observations. Complementing the theoretical analyses in Theorems 2 and 3, we further provide an overall quantitative estimation of 𝜖 th .

Theoretical estimation of measurement error

The decoding stage of the QKM framework requires reconstructing the probability distribution 𝑢 with components 𝑢 𝑖 = |⟨𝑖|𝑢(𝑡)⟩| 2 from measurements. The reconstructed 𝑢 is not derived from a single quantum state, but rather from a linear combination of the outputs from the ℎ parallel circuits, in accordance with the discretized LCHS integral. For each of the ℎ circuits, standard projective measurements in the computational basis are performed with 𝑀 shots, leading to a statistical reconstruction error. Let the exact target probability distribution vector of the 𝑗-th circuit be 𝑢 ( 𝑗 ) ∈ R 𝑁 (for 𝑗 = 1, . . . , ℎ), where its 𝑖-th component

P𝑁 ( 𝑗 ) ( 𝑗) ( 𝑗) 𝑢 𝑖 = 𝑝 𝑖 represents the theoretical probability of observing the computational basis state |𝑖⟩, satisfying 𝑖=1 𝑝 𝑖 = 1. The Pℎ ( 𝑗 ) overall reconstructed state vector is defined by the weighted sum 𝑢 = 𝑗=1 𝑤 𝑗 𝑢 . To ensure 𝑢 is a valid probability distribution, P the weights 𝑤 𝑗 > 0 are strictly positive and normalized ℎ𝑗=1 𝑤 𝑗 = 1. The expected absolute error in the 𝐿 2 norm, which defines 𝜖meas , follows directly from Jensen’s inequality

𝜖meas := E ∥ 𝑢ˆ − 𝑢∥ 2

v u √︂ h 𝑁 i t∑︁ 2 ⩽ E ∥ 𝑢ˆ − 𝑢∥ 2 = Var( 𝑢ˆ 𝑖 ).

(S41)

𝑖=1

P ( 𝑗) The combined empirical estimator for the 𝑖-th component of the total state is 𝑢ˆ 𝑖 = ℎ𝑗=1 𝑤 𝑗 𝑢ˆ 𝑖 . Because the ℎ parallel circuits are executed and measured independently, the variance of the combined estimator is the weighted sum of the individual variances Var( 𝑢ˆ 𝑖 ) =

ℎ ∑︁

( 𝑗)

𝑤 2𝑗 Var( 𝑢ˆ 𝑖 ).

(S42)

𝑗=1

A single shot collapses the 𝑛-qubit circuit into one of the 𝑁 = 2𝑛 possible computational basis states. Therefore, measuring the ( 𝑗) circuit 𝑀 times yields an empirical probability distribution that follows a multinomial distribution. The empirical estimator 𝑢ˆ 𝑖 , ( 𝑗) ( 𝑗) ( 𝑗) as the observed frequency of state |𝑖⟩, is an unbiased estimator of the true probability E[𝑢ˆ 𝑖 ] = 𝑝 𝑖 = 𝑢 𝑖 , and its variance is given by ( 𝑗)

( 𝑗)

Var( 𝑢ˆ 𝑖 ) =

( 𝑗)

𝑝 𝑖 (1 − 𝑝 𝑖 ) . 𝑀

(S43)

Substituting Eq. (S42) and Eq. (S43) into Eq. (S41) yields an upper bound v v v u u u u u u t∑︁ t ∑︁ t∑︁ ( 𝑗) ( 𝑗) ( 𝑗) ℎ ∑︁ 𝑁 ℎ ∑︁ 𝑁 ℎ 𝑝 (1 − 𝑝 ) 𝑝 1 𝑖 𝑖 𝑖 2 2 𝑤𝑗 ⩽ 𝑤𝑗 = 𝑤 2𝑗 . 𝜖meas ⩽ 𝑀 𝑀 𝑀 𝑗=1 𝑖=1

𝑗=1 𝑖=1

(S44)

𝑗=1

For the uniform LCHS quadrature over the truncated interval [−𝐾, 𝐾], the grid spacing is Δ𝑘 = 2𝐾/ℎ, and the Cauchy-Lorentz density is universally bounded by 𝑝(𝑘) ⩽ 1/𝜋. Therefore, the weight is bounded by 𝑤 𝑗 ≈ 𝑝(𝑘 𝑗 )Δ𝑘 𝑗 ≲ 2𝐾/(𝜋ℎ). With the optimal 𝐾 = O (ℎ1/3 ) for 𝜖spec in Eq. (S14), we obtain ℎ ∑︁

𝑤 2𝑗 ⩽

max 𝑤 𝑗

∑︁ ℎ

1⩽ 𝑗⩽ℎ

𝑗=1

2 𝑤 𝑗 = max 𝑤 𝑗 = O ℎ − 3 . 1⩽ 𝑗⩽ℎ

(S45)

𝑗=1

Substituting Eq. (S45) into Eq. (S44) yields √︂ 𝜖 meas ⩽

max1⩽ 𝑗⩽ℎ 𝑤 𝑗 1 1 = O 𝑀 − 2 ℎ− 3 . 𝑀

(S46)

S15 a

b Circuit Haar-random

1.5

10−1 10−2

0.5

10−3 1.0

10−2

100

1.0

0.5 Fidelity

n=8

101

2.0

0.0 0.0

Circuit Haar-random

102

n=2 Probability

Probability

2.5

103

DKL

3.0

10−3

10−4 0.00

0.05 Fidelity

0.10

10−4

FIG. S3. Expressibility analysis of the state preparation ansatz. Empirical fidelity distributions (blue histograms) generated by the topologynative {𝑈3 , CZ} circuits for (a) 𝑛 = 2 and (b) 𝑛 = 8 qubits, compared against the analytical Haar-random probability density function (red dashed lines). The vertical axis in (b) is shown on a logarithmic scale to account for the concentration of measure in high-dimensional Hilbert spaces. (c) KL divergence 𝐷 KL between the circuit-generated fidelity distribution and the Haar-random ensemble as a function of qubit number 𝑛. Markers represent the mean 𝐷 KL value obtained from 5 independent numerical experiments, with error bars denoting the maximum and minimum values. The consistently low 𝐷 KL values indicate that the circuit effectively generates states representative of the Hilbert space, with minimal information loss compared to a Haar-random distribution.

Quantitative analysis of state preparation error

In the absence of ancillary qubits, preparing an arbitrary 𝑛-qubit quantum state requires an exponential circuit depth [S19], with the optimal scaling O (2𝑛 /𝑛) [S20]. Such complexity restricts high-dimensional quantum simulations to theoretical frameworks, rather than real-device execution. By exploiting the intrinsic regularities of real-world physical systems, such as continuity [S21, S22], our approach effectively circumvents the exponential circuit depth typically required for arbitrary quantum state synthesis. The state preparation error 𝜖prep in Eq. (S40) is fundamentally governed by the expressibility of the PQCs used to encode the initial state |𝑢(0)⟩. To quantify the resulting representational capacity of our topology-native {𝑈3 , CZ} ansatz, we evaluate its expressibility by comparing the ensemble of states it generates against the ensemble of Haar-random states [S23]. The evaluation relies on the distribution of state fidelities 𝐹 = |⟨𝜓 𝜃 |𝜓 𝜙 ⟩| 2 , obtained by uniformly sampling parameter sets 𝜃 and 𝜙 to generate state pairs. A perfectly expressive PQC is characterized by a fidelity distribution that replicates the Haar-random ensemble [S23], whose analytical probability density function is given by [S24] 𝑃Haar (𝐹) = (𝑁 − 1) (1 − 𝐹) 𝑁 −2

(S47)

with 𝑁 = 2𝑛 . The expressibility is quantified by the Kullback-Leibler (KL) divergence between the empirical fidelity distribution 𝑃PQC (𝐹) and the reference Haar-random distribution 𝑃Haar (𝐹) , evaluated as a discrete estimator partitioned into 𝑁bins uniformly spaced bins ! 𝑁 bins ∑︁ 𝑃ˆPQC (𝐹𝑖 ) . (S48) 𝐷 KL (𝑃PQC ∥𝑃Haar ) = 𝑃ˆPQC (𝐹𝑖 ) ln 𝑃ˆHaar (𝐹𝑖 ) 𝑖=1

Under this statistical framework, a lower 𝐷 KL value signifies a more expressive circuit that more effectively explores the Hilbert space [S23]. Since the representational capacity of PQCs fundamentally constrains the ability to approximate target states, superior expressivity implies a systematic reduction in the state preparation error 𝜖 prep . Consequently, the KL divergence 𝐷 KL serves as an empirical proxy to estimate 𝜖 prep within the QKM. √ Here we configure the depth parameters 𝑅 and 𝑟 as 𝑅 = 𝑟 = ⌈ 𝑛⌉. This ensures that the total circuit depth scales as 𝑅𝑟 ≈ 𝑛, providing a necessary balance between the circuit’s representational capacity and the coherence constraints of near-term quantum hardware. For a given 𝑛, we uniformly sample pairs of parameter vectors to generate the corresponding quantum states and compute their pairwise fidelities. Because the empirical estimation of the KL divergence is affected by finite-sampling effects and histogram resolution, we set the sample size to 𝑁samples = 5000 and discretize the continuous fidelity range [0, 1] into 𝑁bins = 75 uniformly spaced bins. These parameters are chosen to match the benchmarking in Ref. [S23], ensuring a consistent and valid evaluation. Fig. S3(a) and Fig. S3(b) illustrate the resulting fidelity distributions for 𝑛 = 2 and 𝑛 = 8 qubits, respectively. The empirical histograms generated by our topology-native circuit exhibit alignment with the analytical Haar-random probability density functions. For the larger system size (𝑛 = 8), the fidelity distribution heavily concentrates near zero, reflecting the concentration of measure

S16 ×10−4

102

Learning rate Training loss Validation loss

Loss

101 100 10−1

val = 0.003

10−2 10−3

b 2.4

102

1.9 1.5

×10−4

102

101

100

1.8

100

1.0 10−1

100 Epoch

200

0.1 10−3

30 Epoch

×10−4

2.2 2 1.8

val = 0.057

val = 0.044

0.6 10−2

c 2.2

1.6 10−1

1.6

1.4 10−2

1.4

1.2 10−3

50 Epoch

100

Learning rate

1.2

FIG. S4. Learning dynamics and loss of the QKM framework. The temporal evolution of the composite loss function L and the learning rate schedule for (a) the 3D reaction-diffusion system, (b) the spherical fluid dynamics, and (c) the real-world ocean currents. In each panel, the training loss ℓtrain (blue solid line) and validation loss ℓval (red solid line) are plotted against the left axis on a logarithmic scale. The learning rate schedule (green dashed line), featuring a warm-up phase followed by a cosine decay, is shown on the right axis. The annotated ℓval indicates the terminal validation error achieved at the conclusion of training via an early stopping criterion to prevent overfitting. The minimal gap |ℓval − ℓtrain | between the training and validation curves indicates that the relative optimization error relative 𝐿 2 -𝜖 opt is small, while the terminal training loss ℓtrain quantitatively estimates the theoretical error relative 𝐿 2 -𝜖 th .

typical of high-dimensional Hilbert spaces. In this regime, the circuit accurately reproduces the theoretical Haar distribution over a logarithmic scale, indicating that the parameterization explores the state space uniformly and with minimal bias. A systematic quantitative assessment of this expressibility across system sizes from 𝑛 = 2 to 𝑛 = 9 is presented in Fig. S3(c). The calculated 𝐷 KL values remain consistently low, generally on the order of 10−2 to 10−4 , across all evaluated qubit counts. For 𝑛 = 10, we the concentration of measure in high-dimensional Hilbert spaces requires finer sampling. Using 𝑁samples = 50000 and 𝑁bins = 750, we obtain 𝐷 KL = 1.04 × 10−4 . The KL divergence does not exhibit growth as the system scales; instead, it maintains a stable, low magnitude. This confirms that the representational capacity of the circuit does not degrade for larger systems, ensuring that the PQC can generate states highly representative of the full Hilbert space. By evaluating Eq. (S48) across varying qubit scales 𝑛, we numerically demonstrate that the adopted scaling of 𝑅𝑟 = O (𝑛) provides sufficient expressibility to represent the physical initial conditions, thereby establishing an empirical estimate on 𝜖 prep while circumventing the exponential depth typically required for arbitrary state synthesis.

Quantitative analysis of theoretical and optimization errors

Within the QKM framework, the pre-training phase involves the joint optimization of the neural network encoder-decoder and the unitary Koopman operators. The optimization error 𝜖 opt in Eq. (S40) encapsulates the inaccuracies introduced during this variational training process. To quantify this constituent within the context of physical prediction, we assess it using the relative metric relative 𝐿 2 -𝜖 opt , which provides a direct evaluation of the relative deviation arising from optimization. Because our composite loss function L defined in Eq. (S27) directly computes the expectation of the relative 𝐿 2 error, the optimization error relative 𝐿 2 -𝜖 opt shares the same scale and physical interpretation as relative 𝐿 2 -𝜖. Following statistical learning theory [S25], the expected error can be decomposed into the empirical training risk and the generalization gap. Therefore, without the need for scaling factors, we estimate the relative theoretical and optimization error as 𝜀th := relative 𝐿 2 -𝜖 th ≈ Ltrain ≈ ℓtrain , 𝜀 opt := relative 𝐿 2 -𝜖opt ≈ |Lval − Ltrain | ≈ |ℓval − ℓtrain |,

(S49) (S50)

where Ltrain and Lval represent the asymptotic values of the loss function evaluated on the training and unseen validation datasets, respectively, while ℓtrain and ℓval denote the terminal training and validation losses achieved at the conclusion of the training process. Figure S4 monitors the training and validation losses across the three benchmarks. For the 3D reaction-diffusion system (Fig. S4a), the training terminates at epoch 195 with ℓtrain = 0.001 and ℓval = 0.003. For the spherical fluid dynamics (Fig. S4b), the process ends at epoch 41, yielding ℓtrain = 0.015 and ℓval = 0.044. The real-world ocean current benchmark (Fig. S4c) concludes at epoch 75 with ℓtrain = 0.020 and ℓval = 0.057. In all three cases, the validation loss closely tracks the training curve without divergence, resulting in a minimal generalization gap |ℓval − ℓtrain | throughout the optimization. Substituting these terminal values into Eqs. (S49) and (S50) provides quantitative measures for both 𝜀 th and 𝜀opt . Across all

S17 benchmarks, the minimal generalization gap confirms that the optimization error 𝜖opt remains a subdominant factor. Instead, the empirical training loss ℓtrain serves as the experimental manifestation of the cumulative theoretical error 𝜖 th , capturing the fundamental constraints imposed by the Koopman projection 𝜖proj , the spectral sampling bias 𝜖spec , and the ansatz structural bias 𝜖ansatz . Comparing these estimates with the total relative 𝐿 2 -𝜖 trajectories reported in the Results section allows us to distinguish the primary error sources across different regimes. For the 3D reaction-diffusion case, the minimal ℓtrain indicates that 𝜖th is well-suppressed, leaving the hardware noise 𝜖noise as the primary source of the total simulation error. Conversely, the higher ℓtrain plateaus observed in the spherical and ocean current benchmarks demonstrate that the underlying theoretical complexity, specifically the multi-scale dynamics, dominates the error budget. In these regimes, the theoretical limit 𝜖th becomes the primary constraint on fidelity, rendering the contribution of 𝜖 noise a secondary factor. Furthermore, the statistical measurement error 𝜖meas can be systematically suppressed by increasing the sampling shot count 𝑀. Consequently, improving simulation accuracy requires a dual approach: mitigating hardware noise for weakly nonlinear systems, and refining the theoretical representation and training convergence for complex, multi-scale dynamics.

IDEAL NOISELESS SIMULATION

This section provides ideal noiseless simulations to serve as a reference for the hardware experiments discussed in the main text. By employing classical emulation of the QKM, we isolate the algorithmic performance from measurement errors 𝜖meas and hardware-induced noise 𝜖 noise .

3D reaction-diffusion systems

The classical emulation of the QKM reproduces the pattern formation of the species concentration 𝑢. As shown in Fig. S5(a), the simulated structures match the ground truth (GT) at 𝑡 = 100, 200, and 300. Quantitative assessment in Fig. S5(b) shows that the relative 𝐿 2 error remains below 0.015 throughout the evolution. This is substantially smaller than the error observed during hardware execution (∼ 0.08), indicating that the theoretical error 𝜖 th is small and hardware noise 𝜖 noise is the dominant error source in the physical experiment. Furthermore, the average energy ⟨𝑢 2 ⟩/2 (Fig. S5c) and kernel density estimations (Fig. S5d) closely align with the GT results, confirming that the theoretical framework accurately captures the dissipative dynamics and statistical properties of the system.

Spherical fluid dynamics

The classical emulation of the QKM reproduces the dynamics of the spherical fluid system. As shown in Fig. S6(a), the simulated vorticity field 𝜔 matches the GT from 𝑡 = 0 to 𝑡 = 100. The relative 𝐿 2 error (Fig. S6b) accumulates gradually but remains below 0.06 throughout the evolution. Additionally, the statistical properties at 𝑡 = 80, including the enstrophy spectrum and the kernel density estimation of 𝜔 (Fig. S6c), align with the GT, confirming that the theoretical framework preserves the spectral backbone of the system. Comparing this noiseless baseline with the hardware execution error reported in the main text, the theoretical error 𝜖th and hardware noise 𝜖noise are of similar magnitude. This indicates that for complex multiscale dynamics, both the theoretical constraints and hardware infidelities contribute comparably to the total simulation error.

Real-world ocean currents

The classical emulation of the QKM reconstructs the mesoscale features of the surface geostrophic velocity field 𝑣 over a 13-day period. Comparing the satellite-derived ground truth in Fig. S7(a) with the emulation results in Fig. S7(b), the simulated fields align with the observational data. The pointwise absolute error (Fig. S7c) shows that deviations are primarily localized along the high-gradient regions of the current. As shown in Fig. S7(d), the 𝜀 𝐿2 exhibits a gradual accumulation typical of chaotic dynamics, reaching approximately 0.21 after 13 days. Because this noiseless baseline error closely matches the total error observed during hardware execution, it indicates that the theoretical error 𝜖th is large and dominates the overall error budget for this real-world benchmark. This systematic bias is primarily attributed to the structural simplification of the evolution operator, where each of the ℎ independent circuits utilizes a single-layer 𝑅 𝑧 ansatz to approximate the diagonalized spectral components. Improving simulation fidelity for complex geophysical flows thus necessitates refining the ansatz expressivity alongside the finite-dimensional Koopman representation.

S18

t = 100

t = 200

t = 300

Sim. t=0

0.020

Error

d Sim. GT

0.5

101

0.4

0.012

KDE

0.3

0.008

10−1

0.2

0.004 0.000

100

u 2 /2

relative L2 -

0.016

100

150 t

200

250

300

0.1

t = 100 t = 200

10−2 0

100 150 200 250 300 t

t = 300 0

0.25

0.5 u

0.75

FIG. S5. Classical emulation of 3D reaction-diffusion systems, comparing with GT. (a) Contours of species concentration 𝑢 at 𝑡 = 0 (initial condition, left) and at 𝑡 = 100, 200, and 300 for GT (top row) and classical emulation. (b) Temporal evolution of the relative 𝐿 2 error. (c) Average energy ⟨𝑢 2 ⟩/2 as a function of time for QKM (blue circles) and GT (red solid line). (d) KDEs of 𝑢 at 𝑡 = 100 (circles/solid), 𝑡 = 200 (triangle/dashed), and 𝑡 = 300 (squares/dash–dot), comparing emulation (blue markers) with GT (red lines).

10.

ABLATION STUDY

To substantiate the theoretical scaling laws and error decomposition of the QKM, we perform ablation studies using the spherical fluid dynamics benchmark through classical simulations. The multi-scale interaction and non-trivial topology of this system provide a representative environment to evaluate how the spectral sampling resolution ℎ affects resource efficiency and how the circuit ansatz structure determines the theoretical accuracy floor. In the following, we provide empirical evidence for the ℎ = O (𝑛) scaling while demonstrating the dominance of 𝜖ansatz within the cumulative theoretical error 𝜖 th .

About Theorem 2

We evaluate the impact of the spectral sampling resolution ℎ by testing various values in the vicinity of 𝑛. According to our error analysis, the selection of ℎ must balance the theoretical spectral sampling error 𝜖 spec established by Theorem 2 against computational resource overhead. Our experimental results demonstrate that the linear scaling ℎ = O (𝑛) is sufficient to suppress 𝜖spec and minimize the cumulative theoretical error 𝜖 th while maintaining a minimal generalization gap. As shown in Fig. S8, increasing ℎ explicitly accelerates the convergence of the training process, with the early stopping criterion triggered progressively earlier at epochs 78, 64, and 41 for ℎ = 8, 16, and 32, respectively. However, the terminal validation loss (ℓval ) plateaus at approximately 0.044 to 0.045 across all cases, indicating no substantial gain in final predictive accuracy despite the increased circuit ensemble size. These results demonstrate that the linear scaling ℎ = O (𝑛) is already sufficient to suppress 𝜖spec and capture the essential spectral components without incurring exponential resource costs. This empirical evidence supports that the QKM achieves physically meaningful accuracy with a polynomial number of parallel circuits, thereby

S19 a GT

t=0

t = 20

= 0.010

= 0.017

t = 40

t = 60

t = 80

t = 100

= 0.026

= 0.030

= 0.034

= 0.048

Sim.

b 0.10

Error

106

0.04

100

0.02

10−2 0

100

10−4

Sim. GT

KDE

0.06

0.00

103

104 Ek

relative L2 -

0.08

108

100

10−1

Sim. GT

100

101

101 k

102

10−2 10−3

10−2

10−1

100

FIG. S6. Classical emulation of spherical fluid dynamics, comparing with GT. (a) Vorticity field 𝜔 at 𝑡 =0, 20, 40, 60, 80, and 100 for GT (top row) and classical emulation, with the relative 𝐿 2 error ℓ indicated beneath each snapshot. (b) Temporal evolution of the relative 𝐿 2 error. (c) Enstrophy spectrum 𝐸 𝜔 as a function of wavenumber 𝑘 (left) and KDE of 𝜔 (right), comparing emulation (blue circles) with GT (red solid line) at 𝑡 = 80.

securing the algorithmic speedup S = O (2𝑛 /𝑛3 ). Furthermore, as we continue to scale up the framework with a larger ℎ, the representational capacity of the model expands. Consequently, larger and more diverse datasets will be required in future large-scale implementations to effectively constrain the optimization landscape.

About Theorem 3

We further decompose the theoretical error 𝜀th to identify the primary bottleneck in simulation accuracy for multi-scale systems. Theorem 3 suggests that the structural bias of the single-layer 𝑅 𝑧 ansatz, 𝜀 ansatz , arises from the high-order Pauli-𝑍 coefficients in the Koopman generator’s spectral decomposition. To empirically evaluate this, we augment the time-evolution block by introducing 𝑅 𝑧𝑧 gates between adjacent qubit pairs {(1, 2), (2, 3), . . . , (𝑛, 1)}, thereby incorporating partial second-order interactions while strictly preserving all other hyperparameters. As illustrated in Fig. S9, this augmentation explicitly drives the terminal training loss ℓtrain , which quantitatively approximates the theoretical error 𝜀th , down from 0.015 to 0.008. The visualized field snapshots are evaluated on the training dataset. Because the training process is governed by an early stopping criterion, this significant reduction represents a genuine suppression of theoretical error rather than an artifact of overfitting. This finding confirms that 𝜀 ansatz is the dominant constituent of 𝜀th in our experiments. While the spectral sampling error 𝜀spec can be theoretically bounded and systematically suppressed by increasing ℎ according to Theorem 2, and the projection error 𝜀 proj is empirically presumed to be small based on established principles of classical reduced-order modeling, the structural bias of the shallow circuit ultimately governs the theoretical accuracy floor. These results highlight that capturing higher-order mode interactions is the critical pathway for extending the QKM toward more strongly nonlinear regimes. However, it is worth noting that while the augmented ansatz successfully minimizes ℓtrain , the validation loss remains largely unchanged. This plateau implies that merely enhancing the expressivity of the quantum circuit is insufficient on its own. To effectively translate the reduced theoretical error into improved generalization and further lower the optimization error 𝜀 opt , it is imperative to scale up the training datasets to adequately constrain the expanded parameter space of the augmented higher-order

S20 a

2024.01.01

52◦ N

2024.01.04

2024.01.07

2024.01.10

2024.01.13

1.75

1.50 20◦ N

1.25

52◦ N

Sim.

1.00

0.75 20◦ N

0.50

52◦ N

0.25

0.00 20 N

= 0.001

◦

65◦ W

= 0.045

33◦ W

65◦ W

33◦ W

2024.01.03

2024.01.05

= 0.107 65◦ W

33◦ W

= 0.168 65◦ W

33◦ W

= 0.211 65◦ W

33◦ W

relative L2 -

d 0.2 0.1 0.0

2024.01.01

2024.01.07 Date

2024.01.09

2024.01.11

2024.01.13

FIG. S7. Classical emulation of real-world Gulf Stream ocean currents, comparing with the observational data. (a) Surface geostrophic velocity 𝑣 from satellite-derived GT at five dates spanning January 1-13 2024. (b) Corresponding emulation results. (c) Pointwise absolute error 𝜖 (𝑥, 𝑦; 𝑡) = 𝑣 Sim (𝑥, 𝑦; 𝑡) − 𝑣 GT (𝑥, 𝑦; 𝑡) , annotated with the relative 𝐿 2 error ℓ indicated beneath each snapshot. (d) Relative 𝐿 2 error as a function of date over the 13-day evaluation period.

102

×10−4

h=8

Learning rate Training loss Validation loss

Loss

101 100

b 2.2

102

1.9 1.6

102

101

2.0

101

100

1.7

100

10−1

1.2 10−1

10−2

0.9 10−2

50 Epoch

c 2.2

val = 0.045

10−3

×10−4

0.6 10−3 100 0

h = 16

h = 32

2.2 2 1.8

val = 0.044

val = 0.045

40 Epoch

×10−4

1.5 10−1

1.6

1.2 10−2

1.4

1.0 10−3

30 41 Epoch

Learning rate

1.2

FIG. S8. Ablation study of circuit count ℎ for the spherical fluid dynamics benchmark. The system is configured with (𝑛, 𝑅, 𝑟) = (10, 3, 3). The temporal evolution of the composite loss function L and the learning rate schedule are shown for (a) ℎ = 8, (b) ℎ = 16, and (c) ℎ = 32. In each panel, the training loss ℓtrain (blue solid line) and validation loss ℓval (red solid line) are plotted against the left axis on a logarithmic scale. The learning rate schedule (green dashed line), featuring a warm-up phase followed by a cosine decay, is shown on the right axis. The annotated ℓval indicates the terminal validation error achieved at the conclusion of training via an early stopping criterion.

S21 a

t=0

t = 20

t = 40

t = 60

t = 80

= 0.002

= 0.025

= 0.030

= 0.028

= 0.026

= 0.000

= 0.012

= 0.015

= 0.018

Rz Rz

Sim.

Rz train = 0.015

c Rz

Rzz Rz Rzz Rz

Rz train = 0.008

FIG. S9. Ablation study on the time-evolution ansatz for the spherical fluid dynamics benchmark. The system is configured with (𝑛, ℎ, 𝑅, 𝑟) = (10, 32, 3, 3) and the visualized field snapshots are evaluated on samples from the training dataset. (a) GT snapshots of the vorticity field 𝜔 at 𝑡 = 0, 20, 40, 60 and 80. (b, c) Classical emulation results with the pointwise relative 𝐿 2 error ℓ indicated beneath each snapshot. Specifically, (b) displays the evolution using the baseline single-layer 𝑅 𝑧 ansatz (first-order approximation), achieving a terminal training loss ℓtrain = 0.015; (c) displays the evolution using an augmented ansatz that incorporates adjacent two-qubit 𝑅 𝑧𝑧 gate pairs {(1, 2), (2, 3), . . . , (𝑛, 1)} (capturing partial second-order interactions). Without altering other hyperparameters, the augmented ansatz further reduces ℓtrain to 0.008, directly demonstrating the mitigation of the structural bias 𝜖ansatz . Governed by an early stopping criterion, this reduction represents a genuine suppression of theoretical error rather than overfitting.

ansatz. REFERENCES [S1] D. An, J.-P. Liu, and L. Lin, Linear combination of Hamiltonian Simulation for Nonunitary Dynamics with Optimal State Preparation Cost, Phys. Rev. Lett. 131, 150603 (2023). [S2] R. M. Wilcox, Exponential operators and parameter differentiation in quantum physics, J. Math. Phys. 8, 962 (1967). [S3] R. O’Donnell, Analysis of boolean functions (Cambridge University Press, 2014). [S4] B. Zhang, Z. Lu, Y. Zhao, and Y. Yang, Data-driven quantum Koopman method for simulating nonlinear dynamics, preprint arXiv:2507.21890 (2025). [S5] O. Ronneberger, P. Fischer, and T. Brox, U-net: convolutional networks for biomedical image segmentation, in Med. Image Comput. Comput.-Assist. Interv. (MICCAI) (Springer, 2015) pp. 234–241. [S6] H. Li, J. Xie, C. Zhang, Y. Zhang, and Y. Zhao, A transformer-based convolutional method to model inverse cascade in forced two-dimensional turbulence, J. Comput. Phys. 520, 113475 (2025). [S7] E. Xie, W. Wang, Z. Yu, A. Anandkumar, J. M. Alvarez, and P. Luo, SegFormer: simple and efficient design for semantic segmentation with transformers, Adv. Neural Inf. Process. Syst. 34, 12077 (2021). [S8] BAQIS, Quafu superconducting quantum computing, https://quafu-sqc.baqis.ac.cn (2024). [S9] 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, Proc. Natl. Acad. Sci. U. S. A. 118, e2026805118 (2021). [S10] D. Jennings, K. Korzekwa, M. Lostaglio, A. T. Sornborger, Y. Subasi, and G. Wang, Quantum algorithms for general nonlinear dynamics based on the Carleman embedding, arXiv preprint arXiv:2509.07155 (2025). [S11] I. Joseph, Koopman–von Neumann approach to quantum simulation of nonlinear classical dynamics, Phys. Rev. Res. 2, 043102 (2020). [S12] I. Novikau and I. Joseph, Quantum algorithm for the advection-diffusion equation and the Koopman-von Neumann approach to nonlinear dynamical systems, Comput. Phys. Commun. 309, 109498 (2025). [S13] P. Gray and S. Scott, Autocatalytic reactions in the isothermal, continuous stirred tank reactor: isolas and other forms of multistability, Chem. Eng. Sci. 38, 29 (1983). [S14] The code is available at https://github.com/YYgroup/QKM.

S22 [S15] J. Galewsky, R. K. Scott, and L. M. Polvani, An initial-value problem for testing numerical models of the global shallow-water equations, Tellus Ser. A-Dyn. Meteorol. Oceanogr. 56, 429 (2004). [S16] K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, Dedalus: a flexible framework for numerical simulations with spectral methods, Phys. Rev. Res. 2, 023068 (2020). [S17] E.U. Copernicus Marine Service (CMEMS), Global ocean gridded L4 sea surface heights and derived variables reprocessed 1993 ongoing, https://doi.org/10.48670/moi-00148 (2024). [S18] H. Hersbach, B. Bell, P. Berrisford, S. Hirahara, A. Horányi, J. Muñoz-Sabater, J. Nicolas, C. Peubey, R. Radu, D. Schepers, et al., The ERA5 global reanalysis, Q. J. R. Meteorol. Soc. 146, 1999 (2020). [S19] X.-M. Zhang, T. Li, and X. Yuan, Quantum state preparation with optimal circuit depth: implementations and applications, Phys. Rev. Lett. 129, 230504 (2022). [S20] X. Sun, G. Tian, S. Yang, P. Yuan, and S. Zhang, Asymptotically optimal circuit depth for quantum state preparation and general unitary synthesis, IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst. 42, 3301 (2023). [S21] Z. Meng, X. Zhang, X. Yuan, and Y. Yang, Geometric encoding of turbulence for end-to-end quantum simulation, arXiv preprint arXiv:2508.05346 (2025). [S22] Z. Meng, L. Chen, J.-P. Liu, and G. He, Toward end-to-end quantum simulation of rapidly distorted turbulence, J. Comput. Phys. 558, 114888 (2026). [S23] S. Sim, P. D. Johnson, and A. Aspuru-Guzik, Expressibility and entangling capability of parameterized quantum circuits for hybrid quantum-classical algorithms, Adv. Quantum Technol. 2, 1900070 (2019). [S24] K. Życzkowski and H.-J. Sommers, Average fidelity between random quantum states, Phys. Rev. A 71, 032313 (2005). [S25] I. Goodfellow, Y. Bengio, A. Courville, and Y. Bengio, Deep learning, Vol. 1 (MIT press Cambridge, 2016).

Record · ID 349680 · SHA-256 b26bcf9ed3e73359
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.