Quantum simulation in aeroacoustics

The linearized acoustic Euler equations with prescribed sources are reduced to a finite system of ordinary differential equations by Fourier–Galerkin projection. This system is then lifted to a KvN wavefunction evolution on phase space. The generator is discretized with periodic summation-by-parts (SBP) operators of order 4, 6, or 8 under Weyl symmetrization. These operators are exactly skew-symmetric, so the discrete generator is exactly anti-Hermitian and the propagator is exactly unitary at every resolution. Discretization can therefore degrade accuracy, but it can never violate probability conservation.

The analytic core is an exact operator identity for the commutator of the position operator with the discrete generator. The identity isolates the entire grid-spacing dependence of the moment error into a single residual operator that is computable from the stencil alone. From this identity we prove closed-form bounds on the first- and second-moment errors. The errors decay at the stencil order in the ratio of wavepacket width to grid spacing, and the prefactors distinguish the three Lighthill source types through their modal-amplitude scaling. Inverting the bounds turns a target moment accuracy into the required grid resolution, the qubit count per axis, and the best stencil order, with explicit crossover thresholds between orders.

The circuit construction exploits a structural property of the acoustic drift. No drift component depends on its own coordinate, so the drift and derivative operators act on different tensor factors and commute exactly. The propagator then factors by Strang splitting into one term per phase-space axis. Each term is implemented exactly by three operations, a quantum Fourier transform, a diagonal phase built from the known eigenvalues of the SBP stencil, and the inverse transform. The coupling to the other register reduces to one local diagonal plus one singly controlled diagonal per qubit, because the position value is an affine function of the bits. The only approximation in the evolution is the splitting between the axes. It plays the same role a time integrator plays in a classical solver, and its error is kept below the discretization error. The Gaussian initial state is prepared exactly by the Grover–Rudolph method. The moments are read out by measuring the position registers and averaging, with amplitude estimation available as the large-scale route.

Figure: SBP-KvN propagator circuit at M = 8. Panel (a) shows the Strang ordering of the axis factors, panel (b) the anatomy of one factor, Fourier-conjugated diagonal phases with one controlled diagonal per bit of the coupled register.

The demonstration is self-contained and corresponds to categories (b) and (c) of the call. We implemented the full circuit in Qiskit, executed it on the Aer simulator, and verified it at five levels against exact references.

The first level checks the construction itself. A single splitting step on a grid of 8 points per axis was compared with the exact matrix exponential. The error decreased by factors of 7.87 and 7.97 when the time step was halved, in agreement with the third-order prediction of 8.0. The Grover–Rudolph preparation matched the target Gaussian amplitudes to 4e-16, and the two equivalent forms of the axis factor agreed to 2e-15, so the structured circuit and the fast executable form are the same unitary operation.

The second level reproduces the discretization defect that the bounds control. On a grid of 32 points per axis with 128 splitting steps, the overlap between the circuit wavefunction and the exact discrete evolution differed from one by 1.9e-11. The circuit gave a first-moment defect of +1.303e-3 against the exact value of +1.304e-3, an agreement of 0.07 percent, and both lie below the closed-form bound. We then measured the convergence law from the circuit itself. With the number of splitting steps increased so that the splitting error stays below one percent of the defect, doubling the grid reduced the circuit defect by a factor of 15.31. The exact discrete value is 15.32, and finer exact references continue the sequence through 15.84 and 16.13 toward the limit of 16 for the order-4 stencil. Measuring the position register one million times recovered the moment to 0.07 standard errors. This error accounting is not decoration. An early run in which the splitting error contaminated the measured ratio was caught by exactly this budgeting.

The circuit results rest on a broader classical campaign of 276 exact state-vector runs over the same discretization. The campaign spans one- and two-mode systems at three stencil orders, a 216-case study of all three source types with two spatial profiles and three amplitudes over two decades, and a 36-case study with weakly stratified mean flow. Measured convergence rates matched the predicted exponents within six percent over three decades of error. The structural source-type predictions held to four to six significant digits. The conservatism of the bounds was quantified as well, with a factor of two to ten of slack from the Cauchy–Schwarz step of the proof.

Figure: Classical reference campaign: measured first-moment errors versus grid resolution for the two-mode system at SBP orders 4, 6, and 8, with the predicted slopes.

The scaling argument is a complexity analysis of memory, the dominant constraint for statistical aeroacoustic computations. Classical storage grows as 16 M^d bytes for M grid points per axis in d dimensions, while the register grows as d log2(M) qubits. The bounds make this comparison concrete. At a target accuracy of 1e-3 the order-8 design rule gives seven qubits per axis, so a six-dimensional problem needs about 42 qubits and an eight-dimensional one about 56, against classical state vectors at gigabyte and exabyte scale. The crossover thresholds between stencil orders are explicit, and order 8 is best for any accuracy target below eight percent of the drift norm.

The circuit costs are measured rather than estimated. After compilation to a standard two-qubit gate set, one splitting step costs 1446 entangling gates at 8 qubits and 3666 at 10. This is the M log M cost of synthesizing a general diagonal phase. A route with cost polynomial in the number of qubits is available. It evaluates the stencil eigenvalues into an arithmetic ancilla register and applies the phase from there, and the structure of our diagonals, with the nonlinearity confined to one register, is what makes that route applicable.

The hybridization strategy follows the Lighthill analogy itself. The turbulent source is computed classically by large-eddy or direct numerical simulation, exactly as industry does today, and enters the quantum evolution as a precomputed time-dependent gate schedule with no feedback loop. Classical post-processing converts the measured moments into engineering observables. The quantum register carries only the classically intractable part of the computation, the high-dimensional statistical propagation.

We state plainly what is not yet established. A runtime advantage requires the gate-count analysis of the polynomial route, and that analysis has not been performed. What we claim is a proven register-size requirement, a demonstrated separation of the discretization, splitting, and measurement errors with the discretization side dominant, and direct evidence that a circuit realization inherits exactly the error the analysis predicts. The bounds therefore fix the qubit budget for any future implementation, whichever simulation algorithm it uses.Per-axis qubit count required to enforce the first-moment error bound versus target accuracy, at SBP orders 4, 6, and 8.

Figure: Per-axis qubit count required to enforce the first-moment error bound versus target accuracy, at SBP orders 4, 6, and 8.