Classical simulation of non-Gaussian bosonic circuits

TL;DR

This work establishes a scalable classical framework for simulating non-Gaussian bosonic circuits by augmenting the standard covariance-matrix formalism to track relative phases in a superposition of Gaussian states. It delivers an exact strong-simulation algorithm with runtime $O(s\chi^2 n^3)$ and a faster approximate algorithm with runtime $O\left(s\chi n^3 N e^{2z_{tot}}/(\varepsilon^3 p_f)\right)$, where $\chi$ is the number of Gaussian terms and $z_{tot}$ sums the squeezing parameters; both rely on a triple-overlap formula and a fast norm-estimation method. The framework is particularly relevant for analyzing bosonic error-correction schemes (e.g., GKP and cat codes) and oscillator-to-oscillator encodings, where non-Gaussian resources are interleaved with Gaussian operations. By providing explicit procedures for extending Gaussian descriptions, computing overlaps including phases, evolving under Gaussian unitaries, and handling heterodyne measurements, the approach makes a broad class of non-Gaussian CV circuits tractable for classical simulation and study. This advances the understanding of when and how classical computation can simulate continuous-variable quantum processes, with practical implications for evaluating CV error-correction architectures and near-term simulability of CV quantum information tasks.

Abstract

We propose efficient classical algorithms which (strongly) simulate the action of bosonic linear optics circuits applied to superpositions of Gaussian states. Our approach relies on an augmented covariance matrix formalism to keep track of relative phases between individual terms in a linear combination. This yields an exact simulation algorithm whose runtime is polynomial in the number of modes and the size of the circuit, and quadratic in the number of terms in the superposition. We also present a faster approximate randomized algorithm whose runtime is linear in this number. Our main building blocks are a formula for the triple overlap of three Gaussian states and a fast algorithm for estimating the norm of a superposition of Gaussian states up to a multiplicative error. Our construction borrows from earlier work on simulating quantum circuits in finite-dimensional settings, including, in particular, fermionic linear optics with non-Gaussian initial states and Clifford computations with non-stabilizer initial states. It provides algorithmic access to a practically relevant family of non-Gaussian bosonic circuits.

Figures (7)

• Figure 1: The amplitude $|\langle x,\overline{0}_{\mathrm{GKP}}\rangle|^2$ of the (approximate) GKP state $|\overline{0}\rangle_{\mathrm{GKP}}$ in the position-basis $\{|x\rangle\}_{x\in\mathbb{R}}$.
• Figure 2: The $3$-mode GKP-squeezed repetition encoder by Noh et al. nohoscillatortooscillator
• Figure 3: A graphical representation of the algorithm $\mathsf{overlaptriple}$. Solid lines represent inner products between the states at the vertices that are given / have been computed and are non-zero. Inner products of the form $\langle \psi_3, D(\lambda) \psi_1 \rangle$ are represented by arrows. a) The input to the algorithm $\mathsf{overlaptriple}$ consists of covariance matrices and displacements of three Gaussian states $\psi_1,\psi_2,\psi_3$ together with non-zero overlaps $u = \langle \psi_3, D(\lambda) \psi_1 \rangle, v = \langle \psi_1, \psi_2 \rangle$. b) Applying $\mathsf{overlaptriple}$ provides the inner product $\langle\psi_2,\psi_3\rangle$. In this diagrammatic representation, this completes the triangle with vertices $\psi_1,\psi_2,\psi_3$.
• Figure 4: An illustration of the algorithm $\mathsf{overlap}$. a) The input to $\mathsf{overlap}$ consists of descriptions of two Gaussian states $\psi(\Delta_1),\psi(\Delta_2)$ which includes $\alpha_1,\alpha_2\in\mathbb{C}^n$ and the overlaps $r_j=\langle \alpha_j,\psi(\Delta_j)\rangle, j\in [2]$. b) The algorithm computes $\lambda\in\mathbb{C}^n$ and $\vartheta\in\mathbb{R}$ such that $D(\lambda) |\alpha_1 \rangle = e^{i\vartheta} |\alpha_2\rangle$. This implies that $\langle \alpha_2, D(\lambda) \alpha_1 \rangle = e^{i\vartheta}$ and $\langle \psi(\Delta_2), D(\lambda)\alpha_1 \rangle = e^{i\vartheta} \overline{r_2}$ are known. c) Lastly, the subroutine $\mathsf{overlaptriple}$ is applied to complete a triangle. This amounts to computing the inner product $w=\langle \psi(\Delta_1),\psi(\Delta_2) \rangle$ which the algorithm returns.
• Figure 5: An illustration of the algorithm $\mathsf{squeezing}$. Dashed lines correspond to inner products between the states at the vertices which are known to be non-zero but have not been computed. a) The input to the algorithm $\mathsf{squeezing}$ is the description of a Gaussian state $\psi$ which includes the overlap $r = \langle \alpha, \psi\rangle$, as well as $z\in\mathbb{R}$ and an index $j\in[n]$ associated to a single-mode squeezing operation $S_j(z)$. b) By unitarity of $S_j(z)$, the input data provides the inner product $\langle S_j(z) \alpha, S_j(z) \psi \rangle = \langle \alpha, \psi \rangle = r$. c) The reference state of the evolved state $S_j(z)|\psi\rangle$ is the coherent state $|\alpha'\rangle$ with $\alpha'\in\mathbb{C}^n$ such that $S(z) |\alpha\rangle = D(\alpha') S(z) | 0 \rangle$. This is due to $S_j(z)|\psi\rangle$ having the same displacement vector $d(\alpha')$ as $S_j(z)|\alpha\rangle$. d) The inner product $\langle \alpha', S_j(z) \alpha \rangle$ is computed using \ref{['eq:squeezedcoherentstateCoherentStOverlap']}. e) In the last step, the subroutine $\mathsf{overlaptriple}$ is used to compute $r'=\langle\alpha', S_j(z) \psi\rangle$. Thus, $(\mathop{\mathrm{\Gamma}}\nolimits', \alpha', r')$ is the description of $S_j(z) |\psi\rangle$.

Theorems & Definitions (34)

• Theorem 1.1: Exact strong simulation
• Theorem 1.2: Approximate strong simulation
• Definition 3.1: Extended description of a Gaussian state
• Lemma 3.2
• proof
• Corollary 3.3
• proof
• Lemma 3.4
• proof
• Theorem 3.5