Classical simulability is ultimately determined by both the dynamics of a quantum system and the observables being evaluated. Lie-algebraic simulation exploits the latter to make exact polynomial-time classical simulations by propagating observables through low-dimensional invariant operator spaces. However, its conventional formulation in terms of polynomial-dimensional dynamical Lie algebras does not directly accommodate bosonic systems as their algebras are neither compact nor semisimple. In this contribution, we overcome this limitation, making bosonic systems accessible to the Lie-algebraic programme of exact polynomial-time classical simulation. We prove that expectation values, fixed-order correlation functions, including multi-time correlators and out-of-time-ordered correlators, and gradients are efficiently computable whenever their operator modules have polynomial dimension. This recovers Gaussian quantum optics and extends it to non-Gaussian input states, while identifying exact polynomial regimes of interacting non-Gaussian dynamics including bounded-photon Kerr and pair-hopping Hamiltonians and nilpotent polynomial phase dynamics. We show that unlike in the finite-dimensional spin and fermionic setting treated previously, a finite-dimensional bosonic generator algebra alone does not guarantee finite observable dynamics. We further derive a controlled perturbative hierarchy for squeezing beyond exact sector confinement and confirm the predicted error orders numerically. We also evaluate operator spreading on interacting chains of up to $400$ modes and connect a topological doublon band with flux-reversed edge motion. These results provide a unified formalism for classifying, discovering, and systematically approximating tractable bosonic quantum dynamics with classical polynomial-time simulation.
Comments
20 pages main text, 25 pages appendices, 10 figures
Classical simulability is ultimately determined by both the dynamics of a quantum system and the observables being evaluated. Lie-algebraic simulation exploits the latter to make exact polynomial-time classical simulations by propagating observables through low-dimensional invariant operator spaces. However, its conventional formulation in terms of polynomial-dimensional dynamical Lie algebras does not directly accommodate bosonic systems as their algebras are neither compact nor semisimple. In this contribution, we overcome this limitation, making bosonic systems accessible to the Lie-algebraic programme of exact polynomial-time classical simulation. We prove that expectation values, fixed-order correlation functions, including multi-time correlators and out-of-time-ordered correlators, and gradients are efficiently computable whenever their operator modules have polynomial dimension. This recovers Gaussian quantum optics and extends it to non-Gaussian input states, while identifying exact polynomial regimes of interacting non-Gaussian dynamics including bounded-photon Kerr and pair-hopping Hamiltonians and nilpotent polynomial phase dynamics. We show that unlike in the finite-dimensional spin and fermionic setting treated previously, a finite-dimensional bosonic generator algebra alone does not guarantee finite observable dynamics. We further derive a controlled perturbative hierarchy for squeezing beyond exact sector confinement and confirm the predicted error orders numerically. We also evaluate operator spreading on interacting chains of up to 400 modes and connect a topological doublon band with flux-reversed edge motion. These results provide a unified formalism for classifying, discovering, and systematically approximating tractable bosonic quantum dynamics with classical polynomial-time simulation.
Efficient classical simulation has matured to a critical component of the quantum computing stack, driving hardware validation, algorithm design, benchmarking, and the study of structured quantum dynamics. Lie-algebraic simulation (g-sim) offers a compelling approach: it represents Heisenberg-picture dynamics in the adjoint space whose dimension is set by the dynamical Lie algebra (DLA) governing the circuit, enabling efficient simulation of expectation values whenever the DLA grows only polynomially with system size. Despite this promise, existing applications of g-sim have been confined to free-fermionic settings. It has therefore remained unclear if the method can be applied to other structured circuit families, especially when their generators have large Pauli expansions, and hence whether Lie-algebraic simulability presents a genuinely broader paradigm than free fermions. In this work, we resolve this question by identifying additional non-trivial families of polynomial-dimensional DLAs and introducing symmetry-adapted bases that make the required adjoint-space preprocessing tractable. In particular, we develop an explicit Pauli orbit representation for permutation-equivariant dynamics, enabling efficient processing of cubic-dimensional algebras despite exponential Pauli support, and a modified generalized Gell--Mann representation for bounded Hamming-weight (U(1)-equivariant) dynamics, yielding polynomial simulation costs on fixed excitation sectors. Together with streamlined routines for free-fermionic algebras, these constructions significantly broaden the practical scope of g-sim as a unifying simulation tool for structured quantum circuits. Numerical benchmarks confirm favorable preprocessing scaling and validate large-scale proof-of-concept simulations beyond the reach of state-vector simulation.
We study, classify, and explore the mathematical properties of finite-dimensional Lie algebras occurring in the quantum dynamics of single-mode and self-interacting bosonic systems. These Lie algebras are contained in the real skew-hermitian Weyl algebra A^1, defined as the real subalgebra of the Weyl algebra A1 consisting of all skew-hermitian polynomials. A central aspect of our analysis is the choice of basis for A^1, which is composed of skew-symmetric combinations of two elements of the Weyl algebra called monomials, namely strings of creation and annihilation operators combined with their hermitian conjugate. Motivated by the quest for analytical solutions in quantum optimal control and dynamics, we aim at answering the following three fundamental questions: (i) What are the finite-dimensional Lie subalgebras in A^1 generated by monomials alone? (ii)~What are the finite-dimensional Lie subalgebras in A^1 that contain the free Hamiltonian? (iii) What are the non-abelian and finite-dimensional Lie subalgebras that can be faithfully realized in A^1? We answer the first question by providing all possible realizations of all finite-dimensional non-abelian Lie algebras that are generated by monomials alone. We answer the second question by proving that any non-abelian and finite-dimensional subalgebra of A^1 that contains a free Hamiltonian term must be a subalgebra of the Schrödinger algebra. We partially answer the third question by classifying all nilpotent and non-solvable Lie algebras that can be realized in A^1, and comment on the remaining cases. Finally, we also discuss the implications of our results for quantum control theory. Our work constitutes an important stepping stone to understanding quantum dynamics of bosonic systems in full generality.
We consider the task of simulating the real-time dynamics of weakly interacting fermionic systems. In particular, we focus on computing the expectation value of a local observable A at time t. By analyzing the convergence of the perturbative expansion in the interaction strength λ for the Heisenberg-picture observable, we propose a polynomial-time algorithm for estimating this expectation value in the weakly interacting regime λ∣t∣2D+1=O(1), when the Hamiltonian is geometrically local on a D-dimensional lattice. Importantly, this condition is independent of the system size. If the goal is instead to approximate the time-evolved observable in normalized Frobenius norm, we extend the convergence regime to λ∣t∣=O(1) with quasi-polynomial runtime. When the non-interacting part exhibits Anderson localization, our polynomial-time algorithm can be extended up to λ∣t∣=O(1), modulo polylogarithmic factors. Our algorithm brings together ideas from continuous-time QMC, diagrammatic QMC, and Majorana Propagation, but with a new Heisenberg-picture operator-growth analysis that makes the sampling complexity rigorously controllable. This leads to provably efficient classical algorithms in regimes where the interaction is weak enough that the sampling variance remains bounded independently of system size. Together, these results identify broad regimes in which weak interactions, locality, and localization can be leveraged to make real-time fermionic dynamics classically tractable.
We introduce coherent-state propagation, a computational framework for simulating bosonic systems. We focus on bosonic circuits composed of displaced linear optics augmented by Kerr nonlinearities, a universal model of bosonic quantum computation that is also physically motivated by driven Bose-Hubbard dynamics. The method works in the Schrödinger picture representing the evolving state as a sparse superposition of coherent states. We develop approximation strategies that keep the simulation cost tractable in physically relevant regimes, notably when the number of Kerr gates is small or the Kerr nonlinearities are weak, and prove rigorous guarantees for both observable estimation and sampling. In particular, bosonic circuits with logarithmically many Kerr gates admit quasi-polynomial-time classical simulation at exponentially small error in trace distance. We further identify a weak-nonlinearity regime in which the runtime is polynomial for arbitrarily small constant precision. We complement these results with numerical benchmarks on the Bose-Hubbard model with all-to-all connectivity. The method reproduces Fock-basis and matrix-product-state reference data, suggesting that it offers a useful route to the classical simulation of bosonic systems.
Abstract Simulating non-Gaussian quantum dynamics in open continuous-variable systems is notoriously challenging: exact methods scale exponentially, while approximations sacrifice non-classical features such as Wigner negativity. Here we introduce an exact finite-dimensional reduction that overcomes this curse of dimensionality. We prove that the algebra of extended Gaussian quasi-probability densities (EGQPDs) is closed under Lindblad evolution with linear jump operators, reducing the infinite-dimensional Wigner PDE to a closed system of ODEs that scales polynomially with the number of modes. A necessary and sufficient condition for closure is established: the algebra remains closed iff every jump operator is at most linear; otherwise the extended algebra with polynomial prefactors is required. We further show that the minimal number of Gaussian components defines a non-Gaussian rank—a discrete complexity measure monotonically non-increasing under the dynamics. The framework yields explicit ODEs for the component parameters, admits a rigorous density theorem and error propagation bound, and provides closed-form algebraic criteria for Gaussianity, purity, and entanglement. Seven numerical experiments validate the method across Gaussian and non-Gaussian states, entanglement decay, PT-symmetric exceptional-point dynamics, and complexity scaling, demonstrating machine-precision accuracy and exponential speedups over grid and Fock methods. The EGQPD framework provides a systematic, exact, and efficient toolbox for non-Gaussian open quantum dynamics.
Boson sampling demonstrates quantum advantage through the interference of indistinguishable particles, with output probabilities governed by matrix permanents. Realizing it on deterministic, matter-based platforms requires encoding the bosonic modes in finite-dimensional local Hilbert spaces, which introduces a leakage channel absent in linear optics: multi-particle bunching beyond the local truncation d. We develop a unified framework for non-interacting sampling on the irreducible representations of compact Lie groups, in which the transition amplitude is the immanant of a submatrix of the single-particle transition matrix, recovering the permanent in the bosonic case. Within this framework we bound the bunching leakage through a Dyson-series analysis: decomposing the correlated many-body leakage operator into independent random matrices and applying non-commutative concentration inequalities, we prove, in a Gaussian model of the transition matrix, that its spectral norm concentrates at O~(n) rather than the O(n) worst-case of prior spin-based emulations; the passage to the physical Haar ensemble is reduced to a single submatrix-comparison input, verified at leading order. Exact numerics across local dimensions d=2--5 indicate that the bound is tight, the Haar-ensemble norm matching the closed form d(n−d+1) to sub-percent accuracy. This tightens the required mode number from m=Ω(n4) to the near-optimal m=Ω~(n1+2/(d−1)); for a spin-1 representation (d=3) the overhead falls to m=Ω~(n2), matching the collision-free threshold. The result is independent of particle statistics and applies across finite-dimensional Lie-symmetric architectures, quantifying the spatial resources needed to preserve sampling hardness.