Efficient simulation of low-entanglement bosonic Gaussian states in polynomial time

TL;DR

This paper introduces a polynomial-time classical algorithm to efficiently simulate low-entanglement bosonic Gaussian states by converting them into matrix product states (MPS) without computing computationally hard hafnians. The method combines a Gaussian singular value decomposition (GSVD) to compress the mode space with a projected-creation-operator (PCO) mapping to build finite-dimensional MPS tensors, with complexity governed by entanglement rather than mode count. Benchmarking on Jiuzhang 2.0 and 4.0 covariance data demonstrates substantial speedups over hafnian-based tensor-network approaches in the low-entanglement regime, while maintaining high fidelity to exact representations. The framework broadens the applicability of MPS techniques to bosonic systems, offering a scalable path for simulating open or lossy Gaussian states and potentially informing broader quantum many-body problems and Gaussian-state-based simulations.

Abstract

Bosonic Gaussian states appear ubiquitously in quantum optics and condensed matter physics but remain difficult to simulate classically due to the hafnian bottleneck. We present an efficient algorithm that converts pure bosonic Gaussian states into matrix product states (MPSs), with a computational cost governed solely by the entanglement and not by the number of bosonic modes. Our method combines a Gaussian singular value decomposition with a projected-creation-operator mapping that constructs local MPS tensors without computing hafnians. Benchmarking on covariance matrices from the Jiuzhang 2.0 and Jiuzhang 4.0 Gaussian boson sampling experiments demonstrates substantial speedups over previous tensor-network approaches in the low-entanglement regime relevant to lossy devices. The method provides a scalable classical simulation framework for bosonic Gaussian states with limited entanglement and extends the applicability of MPS-based methods to a broad range of bosonic systems.

Figures (5)

• Figure 1: (a) Schematic of GSVD on the $m$-th bond. A pure BGS $|\phi\rangle$ is decomposed as $|\phi\rangle=\langle I_m|(|\phi_L\rangle\otimes|\phi_R\rangle)$, where $|\phi_L\rangle$ ($|\phi_R\rangle$) is a BGS of physical and virtual modes in the subsystem $L$ ($R$) and $|I_m\rangle$ is a maximally entangled BGS of virtual modes. (b) Iterative GSVD steps. At the $m$-th step, the GSVD gives $|\phi^{m-1}_R\rangle=\langle I_m|(|A^m\rangle\otimes|\phi^m_R\rangle)$, where $|\phi^{m-1}_R\rangle$ becomes the input in the $(m+1)$-th GSVD step. For notational simplicity, we display the situation with the number of virtual modes being the same here, i.e., $n^{m-1}_e=n^m_e=n_e$. (c) Conversion of a local BGS $|A^m\rangle$ to an MPS tensor, involving truncation on the virtual and physical fermion occupation numbers. This truncation is performed using a projection operator $P_{\mathbbm{V}_m}$, where $\mathbbm{V}_m$ represents the local kept subspace.
• Figure 2: Performance of the PCO mapping algorithm for the Jiuzhang 2.0 J2-P65-5 system. The wall timed (blue line) and truncation error (red line) as a function of $D^2$ for constructing the $72$ mode. The error bars for the wall time are contained within the markers. All computations were performed on on a laptop laptop.
• Figure 3: Comparison between the hafnian-based and the PCO mapping methods. The performance is benchmarked by constructing an MPS tensor for the 7th mode in Jiuzhang 2.0 J2-P65-5 system. (a, b) Data collapse of the measured wall times for (a) the hafnian-based method, shown as a function of $d^{\alpha}\log(D)$, and (b) the PCO mapping method, as a function of $d^{\beta}D^2$. (c, d) The same wall time results replotted to show their dependence on the local dimension $d$, with separate curves for each value of $D$. The wall times for the hafnian-based method (a, c) are presented on a logarithmic scale. The error bars for the wall time are smaller than the marker size.
• Figure 4: Performance of the PCO mapping algorithm for the Jiuzhang 4.0 S65 system. The wall time (blue line) and the truncation error (red line) as a function of $D^2$ for constructing the MPS tensor on the $2168$th mode. All computations were performed on a workstation with 32 Intel(R) Xeon(R) Gold 6326 CPUs and 512 GB of RAM.
• Figure 5: Evolution of the loss $\mathcal{L}$ [Eq. \ref{['eq:loss']}] as a function of the effective mode number $N_{\text{eff}}$ during the optimization process for S64 and M256. Convergence is achieved after approximately 3800 steps for S64 and 24000 steps for M256. For visual clarity, data points are shown at intervals of every 50 steps for S64 and every 400 steps for M256.