# A complexity transition in displaced Gaussian Boson sampling

Zhenghao Li, Naomi R Solomons, Jacob F F Bulmer, Raj B Patel, Ian A Walmsley

► To cite this version:

Zhenghao Li, Naomi R Solomons, Jacob F F Bulmer, Raj B Patel, Ian A Walmsley. A complexity transition in displaced Gaussian Boson sampling. npj Quantum Information, 2025, 11 (1), pp.119. 10.1038/s41534-025-01062-5 . hal-05175875

HAL Id: hal-05175875

https://hal.science/hal-05175875v1

Submitted on 22 Jul 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# A complexity transition in displaced Gaussian Boson sampling

<figure><img src="image_3.png" /></figure>

Check for updates

Zhenghao Li<fnref n="1" />, Naomi R. Solomons²,³,⁴, Jacob F. F. Bulmer³, Raj B. Patel<fnref n="1" /> & Ian A. Walmsley<fnref n="1" />

Gaussian Boson Sampling (GBS) is the problem of sampling from the output of photon-number-resolving measurements of squeezed states input to a linear optical interferometer. For purposes of demonstrating quantum computational advantage as well as practical applications, a large photon number is often desirable. However, producing squeezed states with high photon numbers is experimentally challenging. In this work, we examine the computational complexity implications of increasing the photon number by introducing coherent states. This displaces the state in phase space and as such we call this modified problem *Displaced GBS*. By utilising a connection to the matching polynomial in graph theory, we first describe an efficient classical algorithm for Displaced GBS when displacement is high or when the output state is represented by a non-negative graph. Then we provide complexity theoretic arguments for the quantum advantage of the problem in the low-displacement regime and numerically quantify where the complexity transition occurs.

An important milestone in quantum computing is the demonstration of quantum advantage, the goal of performing a task on a quantum device while providing unequivocal evidence that the same task cannot be performed on a classical device within a reasonable amount of time. Advances in the coherent control of large quantum systems have led to experiments that are difficult or perhaps even intractable for conventional computing methods to emulate<fnref n="1" />⁻⁷.

Gaussian Boson Sampling (GBS), as a variant of the Boson Sampling problem, has emerged as a route to demonstrating quantum advantage that is particularly suited to near-term optical implementations⁸⁻¹¹. Here, one measures the photon number statistics of a Gaussian state, prepared by passing squeezed states through a Haar random, multi-mode interferometer. Plausible conjectures in complexity theory provide reasonable evidence for the classical hardness of the sampling task that is further supported by algorithmic studies¹¹,¹². Impressive experimental implementations of the scheme have already been demonstrated on the scale of 100–255 photons<fnref n="1" />⁻⁴. They represent major milestones towards demonstrating quantum advantage, while also stimulating the development of new classical algorithms that have significantly raised the bar for exponential quantum advantage in GBS¹²,¹³.

A major challenge faced by these experiments is scaling. The largest experimental demonstrations of GBS in refs. 1–4 have so far only utilised squeezed vacuum states as input to the interferometer. Preparing squeezed vacuum states with large mean number of photons is experimentally challenging and also sensitive to loss. This

poses an obstacle to further scaling up the total photon number in these experiments.

Coherent states are an under-explored resource in the context of GBS. They are incorporated in the original theoretical framework⁹,¹⁰, but their capability has not yet been exploited. This class of states presents a means by which the total number of photons in the system can be easily increased, as they are readily generated by lasers and their amplitude and phase can be easily controlled using conventional optical components. They are also necessary for certain practical applications that map onto the GBS problem, such as molecular vibronic spectra simulations¹⁴⁻¹⁸ and graph similarity measurements¹⁹. Importantly, the preparation of coherent states is tolerant to loss, as a coherent state with photon loss remains a coherent state, albeit with a smaller amplitude.

In contrast to squeezed states, coherent states input into a passive linear-optical networks do not generate entanglement at the output, and a GBS problem that only takes coherent states as input can be exactly solved efficiently on a classical computer²⁰. Therefore, adding coherent states to a GBS problem might seem to make the problem more ‘classical’, diluting the quantum character of the output state. Despite this, a small-scale experimental demonstration of GBS using both squeezed and coherent light has shown a clear distinction in its output photon distribution from a fully classical model, even when there is a considerable ‘classical’ component to the state²¹.

In this work, we analyse the computational complexity of the 'Displaced GBS' (D-GBS) problem. The name derives from the effect of coherent light

<footnote n="1">¹Department of Physics, Imperial College London, London, UK. ²Quantum Engineering Centre for Doctoral Training, Centre for Nanoscience and Quantum Information, University of Bristol, Bristol, UK. ³Quantum Engineering Technology Labs, H. H. Wills Physics Laboratory and Department of Electrical and Electronic Engineering, University of Bristol, Bristol, UK. ⁴Present address: LIP6 CNRS, Sorbonne Université, Paris, France. ☒e-mail: zhenghao.li21@imperial.ac.uk</footnote>
on the phase-space distribution of the optical quantum state, which is to displace it by an amount proportional to the amplitude of the coherent state. The sampling probability function, which is given by the loop-Hafnian function$^{10,22}$, is connected to the matching polynomial that is widely studied in graph theory$^{23-26}$. In particular, we propose a Uniform D-GBS scheme that takes identical displaced squeezed states as input to a Haar random interferometer and analyse its average-case complexity.

In complexity theory, one way to prove that a problem is 'easy' is by designing an efficient classical algorithm to solve it. Since quantum advantage is usually assumed to be an exponential separation in running time, we consider an algorithm to be efficient as long as its running time is sub-exponential. Based on the Taylor approximation method with quasi-polynomial runtime by Barvinok$^{25}$, we describe a new loop-Hafnian estimation algorithm for a D-GBS problem that is either displacement dominant or when it encodes a non-negative weighted graph.

Second, to show that the D-GBS problem can still be 'hard' in certain regimes, we apply complexity-theoretic arguments. We prove the worst-case complexity for estimating the loop-Hafnian at non-zero displacement, which rules out exact solutions to D-GBS by efficient classical methods. To study the approximate hardness of the problem, we derive the matrix distribution in Uniform D-GBS, over which the loop-Hafnians are evaluated. We conjecture that when displacement is low, this distribution leaves little structure to separate the worst- and average-case complexity of the loop-Hafnian, while also not concentrating the loop-Hafnian values around specific values. We provide evidence for the conjectures, and by using a standard technique developed in refs. 8,11,27, we prove, up to the conjectures, that there is a regime where no efficient classical algorithm is expected to solve the Uniform D-GBS problem even approximately. Numerically, we also provide evidence for how the average-case complexity of estimating the loop-Hafnian changes upon adding more displacement, based on the Taylor approximation method.

A second major obstacle to experimentally implementing GBS is photon loss. A recently developed tensor network algorithm has challenged the quantum advantage claims of the latest GBS experiments by exploiting their high levels of loss$^{28}$. While the addition of displacement cannot directly mitigate the effects of loss—indeed, adding displacement to an already lossy GBS implementation does not make it harder to simulate using the method of ref. 28—it also does not exacerbate them. Despite adding more 'classical' light, the addition of displacement does not lower the loss-threshold for the simulability of GBS, while the preparation of coherent states themselves is loss-tolerant, as lost photons in a coherent state can be compensated by preparing it with a larger amplitude. Comparing regular GBS to D-GBS, our work shows a regime where the latter can retain the former's quantum advantage without significantly adding to the experimental difficulty all while presenting further advantages in achieving large photon numbers and mapping to practical applications.

# Results

## D-GBS and the matching polynomial

The D-GBS problem can be visualised as the sampling of subgraphs from a graph representation of the Gaussian state, as shown in Fig. 1. Compared to regular non-displaced GBS, the addition of displacement introduces loops to the graph, which are edges that connect vertices to themselves. We will first describe this graph representation before connecting it to the matching polynomial that is widely studied in graph theory.

Any M-mode Gaussian state can be uniquely represented by a 2M × 2M covariance matrix, $\Sigma$, and a length-2M vector of means, $\mu$. Measuring all M modes with photon-number-resolving detectors (PNRDs) generates a sample $n = (n_1 \dots n_M)^T$. The total photon number is denoted as $N = \sum_{i=1}^M n_i$. In the absence of loss and noise, the pure-state probability of generating sample $n$ is$^{10,22}$:

$$P_n = \frac{\exp\left(-\frac{1}{2}\mu^\dagger\Sigma^{-1}\mu\right)}{\sqrt{\det(\Sigma_Q) \prod_{j=1}^M n_j!}} \times |\text{Haf}(\text{fd}(\mathbf{B}_n, \gamma_n))|^2. \quad (1)$$

where $\Sigma_Q$ is the covariance matrix of Q function, the $\text{Haf}(\cdot)$ denotes the loop-Hafnian function, and the symbol $\text{fd}(\cdot)$ denotes replacing the diagonal of matrix $\mathbf{B}_n$ with the vector $\gamma_n$. In the rest of this paper, we will occasionally drop the $\text{fd}(\cdot)$ for notation simplicity but its use is implicit in the loop-Hafnian function.

The complex symmetric matrix $\mathbf{B}$ and vector $\gamma$ are constructed from $\Sigma$ and $\mu$, and the subscript $n$ denotes repeating the i-th row and column of $\mathbf{B}$ (or the i-th element of $\gamma$) by $n_i$ times. The $N \times N$ matrix $\text{fd}(\mathbf{B}_n, \gamma_n)$ defines the adjacency matrix of an undirected, complex-weighted graph, denoted $G_n$. The loop-Hafnian function, $\text{Haf}(\mathbf{B}_n, \gamma_n)$, is defined as$^{10,22,29}$

$$\text{Haf}(\mathbf{B}_n, \gamma_n) = \sum_{\Phi \in \text{SPM}(G_n)} \prod_{(i,j) \in \Phi} (\mathbf{B}_n)_{ij} \prod_{(k,k) \in \Phi} (\gamma_n)_k \quad (2)$$

where $\text{SPM}(G_n)$ denotes the set of all single-pair matchings on $G_n$. A single pair matching is a subset of edges, including loops, that is incident on every vertex once and only once. A loop is considered to only connect to its incident vertex once.

In this paper, we primarily study D-GBS where every mode has a non-zero displacement, which results in non-zero loop weights on every vertex of the graph representation. In this case, we can define a new matrix $\tilde{\mathbf{B}}$ by rewriting the loop-Hafnian as the same function but on a matrix of diagonal of ones:

$$\text{Haf}(\mathbf{B}_n, \gamma_n) = \left( \prod_{i=1}^{M} \gamma_i^{n_i} \right) \text{Haf}(\tilde{\mathbf{B}}_n, \mathbf{1}_N), \quad (3)$$

where $\mathbf{1}_N$ is a length-N vector $\mathbf{1}_N = (1 \dots 1)^T$, the matrix $\tilde{\mathbf{B}}$ is constructed by

$$\tilde{\mathbf{B}}_{ij} = \frac{B_{ij}}{\gamma_i \gamma_j}, \quad (4)$$

and the subscript $n$ denotes the same repetition of rows and columns.

The inclusion of loops widens the range of graph problems that map to the quantum sampling problem. Furthermore, the problem of counting single-pair matchings in D-GBS is exactly the monodimer problem in statistical physics$^{23}$, and Equation (3) allows for the loop-Hafnian to be rewritten as the matching polynomial in graph theory, which enjoys a long line of research in computational complexity$^{24-26,30,31}$.

On a complex-weighted loopless graph $\tilde{\mathbf{G}}$ with an $N \times N$ adjacency matrix $\tilde{\mathbf{A}}$, the matching polynomial is defined as

$$g(z; \tilde{\mathbf{A}}) = \sum_{\Theta \in \mathcal{M}(\tilde{\mathbf{G}})} \left( \prod_{(i,j) \in \Theta} \tilde{\mathbf{A}}_{ij} \right) z^{|\Theta|}, \quad (5)$$

where $\mathcal{M}(\tilde{\mathbf{G}})$ is the set of all matchings on $\tilde{\mathbf{G}}$. A matching, $\Theta$, is a subset of edges that is incident on every vertex at most once. The cardinality $|\Theta|$ is the number of edges inside $\Theta$, which is upper bounded by $|\Theta| \le \lfloor \frac{N}{2} \rfloor$. As such, the maximum degree of the matching polynomial is $\lfloor \frac{N}{2} \rfloor$.

If we add a loop of weight one to every vertex in $\tilde{\mathbf{G}}$, to construct a new graph $\mathbf{G}$ with adjacency matrix $\text{fd}(\tilde{\mathbf{A}}, \mathbf{1}_N)$, then it is easy to see that summing over the matchings of $\tilde{\mathbf{G}}$ is equivalent to summing over the single-pair matchings of $\mathbf{G}$. Every matching, $\Theta$, in graph $\tilde{\mathbf{G}}$ maps to a unique single-pair matching, $\Phi$, in graph $\mathbf{G}$ by connecting the vertices outside the matching to themselves via loops. This is illustrated with an example in Fig. 2a, b. In other words, the matching polynomial $g(z; \tilde{\mathbf{A}})$ can be equivalently defined as:

$$g(z; \tilde{\mathbf{A}}) = \text{Haf}(z\tilde{\mathbf{A}}, \mathbf{1}_N). \quad (6)$$
A comparison between GBS and D-GBS

<table><caption>**Fig. 1 | A comparison between GBS and D-GBS with M = 8 modes and K = 2 squeezed vacuum states.** The interferometer is described by the unitary matrix U and multi-mode operator ̂U(U). An example detection sample of n = (20030001)<sup>T</sup> is shown. In D-GBS, the problem is described by a complex symmetric matrix fd(B, γ), which defines an undirected, weighted graph, G. The edges are weighted by B<sub>ij</sub>, while the loop-weights by γ<sub>k</sub>. The detection of n selects a subgraph and repeats the i-th vertex by n<sub>i</sub> times to construct G<sub>n</sub>. Repeated vertices</caption><thead><tr><th></th><th>GBS</th><th>D-GBS</th></tr></thead><tbody><tr><th>Scheme</th><td></td><td></td></tr><tr><th>Graph representation</th><td></td><td></td></tr><tr><th>Output state</th><td>|ψ<sub>out</sub>⟩ = ̂U(U) (∏<sup>K</sup><sub>i=1</sub> ̂S<sub>i</sub>(r<sub>i</sub>)) |0⟩</td><td>|ψ<sub>out</sub>⟩ = (∏<sup>M</sup><sub>j=1</sub> ̂D<sub>j</sub>(α<sub>j</sub>)) ̂U(U) (∏<sup>K</sup><sub>i=1</sub> ̂S<sub>i</sub>(r<sub>i</sub>)) |0⟩</td></tr><tr><th>Gaussian state description</th><td>Σ; μ = 0</td><td>Σ; μ = (α<sub>1</sub> ... α<sub>8</sub> α<sub>1</sub><sup>*</sup> ... α<sub>8</sub><sup>*</sup>)<sup>T</sup></td></tr><tr><th>Pure state probability</th><td>p<sub>n</sub> ∝ |Haf(B<sub>n</sub>)|^2</td><td>p<sub>n</sub> ∝ |Haf(Haf(B<sub>n</sub>, γ<sub>n</sub>))|^2</td></tr><tr><th>Counting problem</th><td>Perfect matchings in G<sub>n</sub></td><td>Single-pair matchings in G<sub>n</sub></td></tr></tbody></table>

are connected to each other by the diagonal terms of the B matrix, B<sub>ii</sub>. The probability of sampling n is proportional to |Haf(Haf(B<sub>n</sub>, γ<sub>n</sub>)|<sup>2</sup>, which sums over the single-pair matchings of G<sub>n</sub>. An example of a single-pair matching is highlighted in red. In the absence of displacement, γ = 0, the graph reduces to a loopless graph. The loop-Hafnian becomes the Hafnian, Haf(Haf(B<sub>n</sub>, 0<sub>N</sub>) = Haf(B<sub>n</sub>), which sums over the perfect matchings of G<sub>n</sub><sup>9,10</sup>, an example of which also highlighted in red.

To put this in the context of D-GBS, evaluation of Equation (3) can be rewritten as evaluating a matching polynomial at z = 1:

$$|\text{Haf}(\tilde{B}_n, \mathbf{1}_N) = g(1; \tilde{B}_n). \qquad (7)$$

Multiplicative-error approximation of the matching polynomial is a #P-hard problem<sup>24,26</sup>. By Equation (6), multiplicative-error approximation of the loop-Hafnian would then also be #P-hard. We define multiplicative-error approximation as follows: an estimator q approximates some value p to within multiplicative-error η if |q - p| ≤ η|p|. On the other hand, an additive-error approximation is defined as giving an estimator q, such that |q - p| ≤ ε, where ε is the additive error.

The complexity class of a problem is a statement of its worst-case complexity. There can, however, be classes of B̃n instances with a special structure that renders the loop-Hafnian, Haf(B̃n, 1<sub>N</sub>), efficient to estimate. In later sections we will see such complexity transitions for two special structures in the B̃ matrix: when its off-diagonal magnitudes are small and when its entries are non-negative.

## Hiding in uniform D-GBS

For the general D-GBS problem, the matrix B̃ in Equation (3) can be constructed to be any arbitrary symmetric matrix by suitable choice of squeezing, displacement and interferometer parameters<sup>32-34</sup>. In this section,

we propose a scheme, named Uniform D-GBS, following the ensuing requirements. The submatrices B̃n will then 'hide' a distribution related to independent and identically distributed (i.i.d.) Gaussian random variables, over which we will study the average-case complexity of the loop-Hafnian. We will discuss how to relax these constraints in an experimental implementation in the Methods section.

1. The input states are K identical displaced squeezed vacuum states, $\prod_{j=1}^{K} \tilde{D}_j(\beta) \hat{S}_j(r)|0\rangle$, with displacement parameter β and squeezing parameter r;

2. The M-mode interferometer is characterised by a Haar random unitary matrix U;

3. The number of non-vacuum input modes is upper-bounded by K ≤ √M; and

4. The total mean photon number is set at one photon per non-vacuum input mode: N̄ = K, such that |β|^2 + sinh²r = 1.

An illustration of the scheme is shown in Fig. 3. By Condition 1, the pair (B̃n, γn) for some sampling outcome n is derived to be:

$$B_n = \tanh(r) U_{n,1_K} U_{n,1_K}^T, \qquad (8)$$

$$\gamma_n = (\beta^* - \beta \tanh(r)) U_{n,1_K} \mathbf{1}_N, \qquad (9)$$
<figure><img src="image_4.png" /><figcaption>**Fig. 2** | Matchings and single-pair matchings. a Enumerating the matchings on some graph $\tilde{G}$ without loops. b If loops are added to every vertex of $\tilde{G}$, every single-pair matching on the updated graph $G$ has a one-to-one mapping to a unique matching in $\tilde{G}$. If all loops in $G$ are weighted by one, then the loop-Hafnian on $G$ is equivalent to the matching polynomial on $\tilde{G}$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3** | Uniform D-GBS example. Identical displaced squeezed states are input into $K$ input modes. The input states are interfered on an $M$-mode interferometer with $M \ge K^2$ and is characterised by a Haar random unitary $U$. After the interferometer the output photons are measured by PNRDs.</figcaption></figure>

where, without loss of generality, the squeezing parameter $r$ is assumed to be positive real, as its phase can be absorbed into the interferometer. The matrix $U_{n,1_K}$ is constructed by repeating the $i$-th row of $U$ by $n_i$ times and only keeping the first $K$ columns. Matrix multiplication with the column vector $1_N$ returns sums of each row. More explicitly, this is:

$$ \gamma_i = (\beta^* - \beta \tanh(r)) \sum_{j=1}^{K} U_{ij}. \qquad (10) $$

Since $U$ is Haar random, in general we have $\gamma_i \neq 0$ for all $i$.

Equations (8) and (9) allow us to separate the squeezing and displacement parameters from the interferometer transfer matrix. We go one step further and define a complex parameter $w$:

$$ w = \frac{\beta^* - \beta \tanh(r)}{\sqrt{\tanh(r)}}. \qquad (11) $$

The magnitude of $w$ characterises the ratio of displacement to squeezing in the Uniform D-GBS scheme. For a fixed $w$, one can solve for the values of $r$ and $\beta$ using Eq. (11) and Condition 4.

Using Eqs. (8), (9) and (11), the probability of sampling $n$ can be rewritten as

$$ p_n \propto |\mathrm{Haf}(U_{n,1_K} U_{n,1_K}^T, wU_{n,1_K} 1_N)|^2 \qquad (12) $$

where $w$ becomes a scaling factor on the diagonal of the loop-Hafnian. In other words, changing the displacement and squeezing in Uniform D-GBS only changes the relative weights of the diagonal and off-diagonal elements in the loop-Hafnian. The distribution of these elements are determined by the matrices $U_{n,1_K}$.

Reference 8 showed that, when $\bar{N} \le \sqrt{M}$, the likelihood of collision events, i.e. detecting more than one photon in any output mode, is small – this is the ‘collision-less’ condition. This means we can treat $U_{n,1_K}$ as an $N \times K$ random submatrix of some Haar random $U$.

For output samples with total photon number upper-bounded by $N \le \sqrt{M}$, which is a substantial subset of all outcomes given Condition 4, Ref. 8 further conjectured that the submatrix $U_{n,1_K}$ is approximately an independent and identically distributed (i.i.d.) complex Gaussian matrix $X$ with mean zero and variance $1/M$, denoted as

$$ U_{n,1_K} \sim X \in \mathcal{G}_{N,K}(0, 1/M). \qquad (13) $$

This is the ‘hiding’ conjecture.

Together, the collision-less condition and the hiding conjecture give the loop-Hafnian function in Eq. (12) an underlying Gaussian distribution. Using the same loop-Hafnian rescaling technique in Eq. (3), we can rescale
this loop-Hafnian to be:

$$ 
\begin{aligned}
& \mathrm{Haf}(\mathbf{X}\mathbf{X}^T, w\mathbf{X}\mathbf{1}_N) \\
& = \left(\prod_{i=1}^{M} (\gamma_i)^{n_i}\right) \mathrm{Haf}\left(\frac{1}{w^2}\tilde{\mathbf{X}}, \mathbf{1}_N\right) \\
& = \left(\prod_{i=1}^{M} (\gamma_i)^{n_i}\right) g\left(\frac{1}{w^2}; \tilde{\mathbf{X}}\right),
\end{aligned} 
\quad (14) $$

where matrix $\tilde{\mathbf{X}}$ is constructed from i.i.d. Gaussians $\mathbf{X} \in \mathcal{G}_{N,K}(0, 1/M)$:

$$ \tilde{\mathbf{X}}_{ij} = \frac{(\mathbf{X}\mathbf{X}^T)_{ij}}{(\mathbf{X}\mathbf{1}_N)_i (\mathbf{X}\mathbf{1}_N)_j}. \quad (15) $$

The distribution of $\tilde{\mathbf{X}}$ is preserved if the Gaussian distribution is rescaled to $\mathbf{X} \in \mathcal{G}_{N,K}(0, 1)$. We denote the matrix ensemble defined in Equation (15) as $\tilde{\mathbf{X}} \in \tilde{\mathcal{G}}_{N,K}(0)$.

### Diagonally-dominant loop-Hafnians in general D-GBS

Physically, the D-GBS problem is expected to become simulable when there is a large amount of displacement relative to squeezing$^{21}$. Mathematically, this corresponds to a small $|\tilde{\mathbf{B}}_{ij}|$ for $\tilde{\mathbf{B}}$ matrix in Eq. (3), or a large $|w|$ for the $w$ parameter in Eq. (14), both of which result in a diagonally-dominant matrix on which the loop-Hafnian is calculated.

In this special case, the loop-Hafnian can be approximated in quasi-polynomial time by the Taylor approximation method developed in refs. 25,30,35-37. The method is applicable to any polynomial that is non-zero for some domain in the complex plane and works by truncating the Taylor expansion of its logarithm in that domain.

More specifically, for an $N \times N$ matrix $\tilde{\mathbf{A}}$, if there is some disc in the complex plane with radius $R > 0$ such that the matching polynomial $g(z; \tilde{\mathbf{A}}) \neq 0$ for all $|z| < R$, then inside the non-zero disc the Taylor truncation can approximate $\ln g(z; \tilde{\mathbf{A}})$ to within arbitrary additive error $\epsilon > 0$ and can be computed in quasi-polynomial time $N^{O(\ln N - \ln \epsilon)}$. The detailed derivations are given in the 'Methods' section. If one can approximate $\ln g(z)$ to any additive error $\epsilon > 0$, then one can approximate $g(z; \tilde{\mathbf{A}})$ to any multiplicative error $\eta = \epsilon^\eta - 1$.

To establish a connection with diagonally-dominant loop-Hafnians, let's consider an $N \times N$ matrix $\tilde{\mathbf{A}}$ that satisfies $|\tilde{\mathbf{A}}_{ij}| \ll 1$, for all $i, j$. Using Equation (6), the loop-Hafnian, $\mathrm{Haf}(\tilde{\mathbf{A}}, \mathbf{1}_N)$, can be expressed as:

$$ \mathrm{Haf}(\tilde{\mathbf{A}}, \mathbf{1}_N) = 1 + \sum_{\{(i,j)\} \in \mathcal{M}} \tilde{\mathbf{A}}_{ij} + \sum_{\{(i,j),(k,l)\} \in \mathcal{M}} \tilde{\mathbf{A}}_{ij} \tilde{\mathbf{A}}_{kl} + \dots \quad (16) $$

If the magnitude of the matrix elements, $|\tilde{\mathbf{A}}_{ij}|$, are sufficiently small, the magnitude of the higher-order terms decay to zero and $\mathrm{Haf}(\tilde{\mathbf{A}}, \mathbf{1}_N)$ is non-zero.

Now suppose we can find a bound $0 < \zeta < 1$, such that Eq. (16) is non-zero if $|\tilde{\mathbf{A}}_{ij}| < \zeta$ for all $i, j$. In order to estimate $\mathrm{Haf}(\tilde{\mathbf{A}}, \mathbf{1}_N)$, we can fix a real number $\lambda$ satisfying $0 < \lambda < 1$, such that $|\tilde{\mathbf{A}}_{ij}| \le \lambda\zeta$ for all $i, j$. Then, the matching polynomial $g(z; \tilde{\mathbf{A}})$ is non-zero for any $|z| < 1/\lambda$. This enables the Taylor approximation method to approximate the matching polynomial at $z = 1$, which is precisely the loop-Hafnian value we wanted to estimate: $\mathrm{Haf}(\tilde{\mathbf{A}}, \mathbf{1}_N) = g(1; \tilde{\mathbf{A}})$. The associated runtime is $N^{O_\lambda(\ln N - \ln \epsilon)}$, where the subscript in the big-O notation denotes that the implied constant depends only on $\lambda$.

The outstanding task is therefore to identify the bound on $|\tilde{\mathbf{A}}_{ij}|$, such that $\mathrm{Haf}(\tilde{\mathbf{A}}, \mathbf{1}_N)$ is non-zero and can be efficiently estimated. We prove Theorems 1 and 2.

**Theorem 1.** For an $N \times N$ complex symmetric matrix $\tilde{\mathbf{A}}$ that satisfies

$$ |\tilde{\mathbf{A}}_{ij}| < \frac{1}{e(2N - 3)} \quad \text{for all } i \neq j, \quad (17) $$

given any $\epsilon > 0$, the function $\mathrm{Haf}(\tilde{\mathbf{A}}, \mathbf{1}_N)$ can be approximated to within multiplicative error $\epsilon$ in time $N^{O(\ln N - \ln \epsilon)}$.

**Theorem 2.** For an $N \times N$ complex symmetric matrix $\tilde{\mathbf{A}}$ that satisfies

$$ \sum_{j \neq i} |\tilde{\mathbf{A}}_{ij}| < \frac{1}{N-1} \quad \text{forall } i \in [1, N], \quad (18) $$

given any $\epsilon > 0$, the function $\mathrm{Haf}(\tilde{\mathbf{A}}, \mathbf{1}_N)$ can be approximated to within multiplicative error $\epsilon$ in time $N^{O(\ln N - \ln \epsilon)}$.

The proof for Theorem 1 uses results from refs. 38,39, and the proof for Theorem 2 adapts a method from ref. 37. The details of both proofs are provided in the Supplemental MaterialsSI. Equation (18) is typically less applicable than Eq. (17), except for a special matrix $\tilde{\mathbf{A}}$ where a row might have one off-diagonal element that violates latter, but the sums of matrix elements still satisfy the former.

For a general D-GBS problem, if matrix $\tilde{\mathbf{B}}_n$ in Eq. (3) satisfies either condition, then $\mathrm{Haf}(\tilde{\mathbf{B}}_n, \mathbf{1}_N)$ can be estimated within quasi-polynomial time, and one can approximately sample from the target distribution in Eq. (1) using a chain rule method$^{40}$. The probability of these conditions being satisfied increases as the magnitude $|\tilde{\mathbf{B}}_{ij}|$ decreases, which corresponds to a higher ratio of displacement to squeezing, as shown in the Methods section. However, the required ratio could be very significant, if $\tilde{\mathbf{B}}_n$ were to satisfy either condition deterministically for all outcomes $n$.

### Diagonally-dominant loop-Hafnians in uniform D-GBS

If we allow the algorithm to fail for a fraction of outcomes $n$, the requirement on displacement can be made less stringent. For this, we consider the underlying distribution of $\tilde{\mathbf{B}}_n \sim \frac{1}{w^2} \tilde{\mathbf{X}}$ for $\tilde{\mathbf{X}} \in \tilde{\mathcal{G}}_{N,K}(0)$ in our proposed Uniform D-GBS scheme, where $N$ is the total photon number in outcome $n$ and $K$ is the number of non-vacuum input modes. We use the hiding conjecture and consider where the loop-Hafnian, $\mathrm{Haf}(\frac{1}{w^2} \tilde{\mathbf{X}}, \mathbf{1}_N)$, or equivalently the matching polynomial, $g(\frac{1}{w^2}; \tilde{\mathbf{X}})$, is on average efficient to estimate.

We approach the problem numerically. Fig. 4(a), (b) plot the distribution of zeroes of $g(z; \tilde{\mathbf{X}})$ for $\tilde{\mathbf{X}} \in \tilde{\mathcal{G}}_{N,16}(0)$ for $N \le 16$. Clearly, there is a region near the origin of the complex plane where the zeroes are less dense, and the area of this region increases with decreasing $N$. The larger $|w|$ is, the more likely it is that, for a random $\tilde{\mathbf{X}}$, the matching polynomial satisfies $g(z; \tilde{\mathbf{X}}) \neq 0$ for $|z| \le \frac{1}{|w|^2}$. The Taylor approximation method can then approximate this to arbitrary multiplicative error.

The lack of zeroes near the origin of the complex plane is true for all $\tilde{\mathbf{X}} \in \tilde{\mathcal{G}}_{N \le K,K}(0)$ up to $K = 16$, whose matching polynomials we numerically computed. Fig. 5(a) plots the probability that $g(z; \tilde{\mathbf{X}}) \neq 0$ for $|z| < \frac{1}{|w|^2}$ and $N = K$ as a colour plot. The black crosses indicate 50% probability and the red 75%. The probabilities would be even higher for $N < K$, as illustrated by Fig. 4(b). Given $w$, by solving Eq. (11) under Condition 4, we can solve for the displacement parameter $\beta$ and squeezing $r$. The photon number ratios $\frac{|\beta|^2}{\sinh^2(r)}$ are plotted in Fig. 5(b).

In other words, above the level of displacement indicated by the black (or red) crosses in Fig. 5, the Taylor approximation method can efficiently estimate $\mathrm{Haf}(\frac{1}{w^2} \tilde{\mathbf{X}}, \mathbf{1}_N)$ for at least $1/2$ (or $3/4$) of the matrices $\tilde{\mathbf{X}} \in \tilde{\mathcal{G}}_{N \le K,K}(0)$. This reflects a transition in the average-case complexity of estimating the outcome probabilities of Uniform D-GBS, away from #P-hard, which would invalidate the proofs of hardness we present in later sections. Fig. 5 therefore, maps out a displacement regime which experiments should avoid if the goal is to demonstrate quantum advantage.

Unlike the well-defined conditions in Theorems 1 and 2, the efficient approximability of $\mathrm{Haf}(\frac{1}{w^2} \tilde{\mathbf{X}}, \mathbf{1}_N)$ in this section is conditioned on a gradual change of algorithmic success probability over matrices
**Fig. 4 | The distribution of zeroes of the matching polynomial.** a The distribution of the smallest-magnitude roots of $g(z; \tilde{X})$ for 10,000 $\tilde{X}$ matrices drawn from $\tilde{G}_{16,16}(0)$ (Eq. (15)). There is a region about the origin of the complex plane where the roots are less dense. The black box encompasses 50% of all plotted roots and the red box 25%. In other words, the matching polynomial $g(z; \tilde{X})$ is empirically non-zero with 50% (75%) probability if the complex number $z$ is inside the black (red) box.
b The magnitude distribution of the smallest-magnitude roots of $g(z; \tilde{X})$ for $N \le K = 16$. The magnitudes are higher for smaller $N$. The black and red lines correspond to the black and red boxes in (a), and the line-widths equal the error bars for 95% confidence level. In the Uniform D-GBS scheme, $N$ is the total number of photons and $K$ is the number of non-vacuum input modes.

<figure><img src="image_6.png" /></figure>

$\tilde{X} \in \tilde{G}_{N \le K, K}(0)$. It may be possible to build this into an efficient but probabilistic algorithm for simulating Uniform D-GBS, such as by the chain rule method proposed in ref. 40. But further techniques need to be developed to efficiently estimate the sampling probability of outcomes with photon number $N > \sqrt{M}$, as these cannot be straightforwardly expressed in terms of the $\tilde{X}$ matrix. One possible method is by expanding the larger loop-Hafnians in terms of smaller ones—the loop-Hafnian of any $N \times N$ matrix can be expanded as a sum of $N$ loop-Hafnians on submatrices of dimensions $(N-1) \times (N-1)$.

On the other hand, below the levels of displacement shown in Fig. 5, the increasing probability of success of the Taylor approximation method over matrices $\tilde{X} \in \tilde{G}_{N \le K, K}(0)$ does not directly translate into a proportionate algorithmic speedup in simulating instances of D-GBS. Given a Uniform D-GBS problem, due to the randomness of the matrix distribution, one cannot speed up an existing classical sampling algorithm by simply assigning the ‘easy’ instances to the efficient algorithm. In fact, if the average-case complexity of the loop-Hafnian is #P-hard, then the corresponding sampling problem plausibly remains an asymptotically hard problem, as we discuss in later sections.

One assumption we make is that the complexity transition only depends on the magnitude of $|w|$. An open question is whether the complexity transition depends on the complex phase of $w$ as well, which is determined by the relative phase between the displacement and squeezing parameters.

## Non-negative loop-Hafnians

A second special case we consider is $\text{IHaf}(\tilde{A}, \mathbf{1})$ when the elements of $\tilde{A}$ are non-negative. This is a special structure that has been previously exploited by classical algorithms for the efficient approximation of permanents in Fock-state Boson Sampling and Hafnians in regular GBS$^{41-43}$.

Ref. 25 gives a method that approximates matching polynomials on non-negative weighted graphs based on the same Taylor approximation method. Based on this, we have the following result:

**Theorem 3.** For an $N \times N$ complex symmetric matrix $\tilde{A}$ that satisfies $\tilde{A}_{ij} \ge 0 \ \forall i, j \in [1, N]$, given any $\epsilon > 0$, the function $\text{IHaf}(\tilde{A}, \mathbf{1}_N)$ can be approximated to within multiplicative error $\epsilon$ in time $N^{O(\ln N + ||\tilde{A}|| - \ln \epsilon)}$.

We define the norm of $\tilde{A}$ to be the maximum absolute row sum of the matrix:

$$||\tilde{A}|| = \max_i \sum_{j \ne i} |\tilde{A}_{ij}|. \qquad (19)$$

The runtime complexity $N^{O(\ln N + ||\tilde{A}|| - \ln \epsilon)}$ grows quasi-polynomially if $||\tilde{A}||$ grows only logarithmically. However, if $||\tilde{A}||$ grows faster than polynomially, the runtime complexity will grow exponentially, which would render the algorithm inefficient.

Theorem 3 is particularly relevant in the search for practical applications of D-GBS. Regular GBS without displacement has been found to solve
<figure><img src="image_7.png" /><figcaption>Fig. 5 | Efficient estimation of the matching polynomial. a Colour plot of the probability that g(z; X̃)≠0 for |z| < 1/w for different values of |w| and X̃ ∈ G̃_{K,K}(0). Probabilities estimated by bootstrapping from 10,000 randomly drawn X̃ matrices. Black crosses indicate 50% probability, and red 75%. Error bars are included for 95% confidence level, though for many points not visible on this scale. b The photon number ratio 1/(β|β|²) calculated from Equation (11) for the w values corresponding to the black and red crosses in (a).</figcaption></figure>

certain graph problems, such as dense-subgraph search⁴⁴⁻⁴⁶ and molecular docking⁴⁷,⁴⁸. These graph problems all involve a search for subgraphs with high connectivity and weights, which are preferentially sampled by GBS. Addition of loops may then be used to encode self-interaction in the graph problems, or to boost the sampling probability in a certain subregion by assigning loops to only a subset of vertices. Thus far, these application proposals have been limited to non-negative weighted graphs and as such show little advantage over classical algorithms⁴³. Similarly, Theorem 3 rules out exponential quantum advantage for solving non-negative weighted graph problems by D-GBS.

We note that the runtime complexity in Theorem 3, while asymptotically efficient, can still be very large for a finite-size problem. If the ability of a quantum device that runs D-GBS to solve graph problems is compared to the performance of the Taylor approximation algorithm, the quantum device may still demonstrate a sub-exponential advantage.

Quantum advantage via uniform D-GBS

In addition to studying where D-GBS becomes classically tractable, we also try to answer the question of where D-GBS remains classically intractable, i.e. where no efficient classical algorithm exists to sample from the same distribution. A trivial answer to this question is that D-GBS is classically intractable when it has zero displacement—the loop-Hafnian is worst-case

\#P-hard to estimate, and the worst case occurs when it has a zero diagonal and reduces to the Hafnian.

By connecting the loop-Hafnian to the matching polynomial in Eq. (6), however, we can already provide a better answer. Since multiplicative-error approximation of a matching polynomial is worst-case \#P-hard²⁴,²⁶, by Eq. (6), so must also be a loop-Hafnian with non-zero diagonal. One can ‘hide’ the worst-case instance in the outcomes of D-GBS, since one can encode any arbitrary symmetric matrix into the $\tilde{B}_n$ matrices of Eq. (3). Using Stockmeyer’s theorem, one can then prove that no efficient classical algorithm can solve the D-GBS problem exactly, unless the Polynomial Hierarchy collapses to the third level, which is widely believed to be implausible in complexity theory⁸,¹¹,⁴⁹⁻⁵¹.

However, this argument is not appealing for the purpose of proving quantum advantage, as the noisy quantum devices today are themselves prevented from achieving exact precision. A more common precision target in quantum sampling schemes is defined in terms of the total-variation distance (TVD). For a target distribution $\tilde{p}_n$, consider an efficient classical algorithm C that, if given some error bound $\epsilon$, can sample from some $q_n$ such that

$$ \text{TVD}(\tilde{p}_n, q_n) = \frac{1}{2} \sum_{n \in \Omega} |\tilde{p}_n - q_n| \le \epsilon, \quad (20) $$

where $\Omega$ denotes the sample space. This additive-error approximation becomes multiplicative-error if the $\tilde{p}_n$ distribution satisfies an anti-concentration property over $n$ in the form of

$$ \Pr_{n \in \Omega} \left( \tilde{p}_n \ge \frac{1}{\alpha |\Omega|} \right) > 1 - \eta \quad (21) $$

for some $\alpha, \eta > 0$. Furthermore, if the probability function $\tilde{p}_n$ is on average \#P-hard to approximate to within multiplicative error $\epsilon = O(\epsilon\delta)$ for at least a $(1 - \delta)$-fraction of outcomes $n \in \Omega$, then the approximate sampler C would also collapse the Polynomial Hierarchy⁸,⁵⁰.

For the complexity proof of D-GBS, we utilise the Uniform D-GBS framework, which puts both the anti-concentration and average-case complexity of $\tilde{p}_n$ in context of the Gaussian matrix distribution $X \in \mathcal{G}_{N,K}(0, 1)$ by Eqs. (12) and (13). Interestingly, we observe a transition of both properties with the amount of displacement.

The relevant probability distribution $\tilde{p}_n$ we consider is that of outcomes $n$ post-selected on $N = K \le \sqrt{M}$ photons. The post-selected probability of measuring $n$ is derived to be SI

$$ \tilde{p}_n = \frac{2^N}{F_N(w)} |\text{Haf} (U_{n,1_K} U_{n,1_K}^T, w U_{n,1_K} \mathbf{1}_N)|^2, \quad (22) $$

where $\frac{2^N}{F_N(w)}$ is the normalisation factor and the submatrix $U_{n,1_K}$ is conjectured to hide a Gaussian distribution according to Eq. (13). The function $F_N(w)$ is the convolution of absolute squared Hermite polynomial over $n$:

$$ F_N(w) = \sum_{\substack{m \ge 0 \\ \sum m_j = N}} \prod_{j=1}^{N} \frac{1}{m_j!} \left| H_{m_j} \left( \frac{iw}{\sqrt{2}} \right) \right|^2. \quad (23) $$

Since the average total photon number of Uniform D-GBS is designed to be $\bar{N} = K$, the post-selection of $N = K$ photon-coincidence outcomes succeeds with at least inverse polynomial probability (see Supplemental MaterialsSI). Therefore, any classical algorithm that can approximately sample from the non-post-selected distribution, $p_n$, can also approximately sample from $\tilde{p}_n$ within polynomial time. But the latter implies the unlikely collapse of the Polynomial Hierarchy, if the loop-Hafnian in Eq. (22) is average-case \#P-hard and adequately anti-concentrated. A definite proof of either property has remained elusive for various Boson Sampling proposals so far. In what follows, we make two conjectures on how these two properties for the loop-Hafnian changes with the $w$ parameter from Eq. (11).
Conjecture 4. (Average-case complexity of the loop-Hafnian) There is some positive real $\tilde{w} > 0$, such that for all complex $w$ that satisfies $|w| \le \tilde{w}$, given any constants $\epsilon, \delta > 0$, it is #P-hard to estimate $\mathrm{Haf}(XX^T, wX1_N)$ to within multiplicative error $\epsilon$ with probability at least $1 - \delta$ over the choice of $X \in \mathcal{G}_{N,N}(0, 1)$.

Conjecture 5. (Anti-concentration of the loop-Hafnian) There is some positive real $\tilde{w} > 0$, such that for all complex $w$ that satisfies $|w| \le \tilde{w}$, for all positive integers $N$ and real $\eta > 0$, there exists a polynomial $\alpha = \mathrm{poly}(N, 1/\eta)$ such that

$$ \mathrm{Pr}_{X} \left[ |\mathrm{Haf}(XX^T, wX1_N)|^2 < \frac{1}{\alpha} \frac{N! F_N(w)}{2^N} \right] < \eta \quad (24) $$

over $X \in \mathcal{G}_{N,N}(0, 1)$.

We consider the evidence for Conjectures 4 and 5, in particular noting that for smaller values of $\tilde{w}$, they are more likely to be true. More detailed discussion is provided in Methods, with further support given in Supplemental MaterialsSI.

Essentially, Conjecture 4 suggests that for some $w$, the worst-case and average-case complexity of the loop-Hafnian are equivalent and both #P-hard over $X \in \mathcal{G}_{N,N}(0, 1)$. Proof of average-case complexity conjectures akin to Conjecture 4 is a well-known open problem across all quantum random sampling schemes$^{8,11,27,50,52,53}$. Established worst-to-average-case reduction techniques based on polynomial extrapolation$^{8,54}$ do not easily extend to complex fields and have failed to converge to sufficient precision for robust quantum advantage$^{8,50,52}$. In the Supplemental MaterialsSI, we describe how these methods can be adapted to prove exact computation of the loop-Hafnian is average-case #P-hard over finite fields, and multiplicative-error estimation is average-case #P-hard over some distributions in real fields. However, significant gaps still remain in the full proof of Conjecture 4.

Despite the lack of rigorous proof, we still believe Conjecture 4 is plausibly true for some $\tilde{w}$ value. One reason is the apparent lack of structure in the Gaussian distribution, $X \in \mathcal{G}_{N,N}(0, 1)$, that could enable a classical algorithm to solve a random instance more efficiently than the worst-case instance. In similar contexts, average-case and worst-case complexity of $\mathrm{Per}(X)$ (or $\mathrm{Haf}(XX^T)$) were also conjectured to be equivalent in Boson Sampling (or Gaussian Boson Sampling)$^{8,11}$.

A second reason is the reduction of $\mathrm{Haf}(XX^T, wX1_N)$ to the Hafnian in GBS when $w$ approaches 0, whose average-case approximate complexity is similarly conjectured and supported by complexity-theoretic$^{11}$ and algorithmic evidence$^{12,29}$.

Extending the Hafnian to the loop-Hafnian, the most obvious structure that can be exploited to separate their complexities are the diagonal weights, which is precisely what is considered previously in this work. We developed an efficient Taylor approximation method to estimate diagonally-dominant loop-Hafnians by exploiting the distribution of their zeroes. Then, in Fig. 5, we observed a transition in the algorithm's average success probability over $X \in \mathcal{G}_{N,N}(0, 1)$ as $w$ decreases.

We conjecture that this observation is not isolated to the specific algorithm, but that the average-case complexity of the loop-Hafnian is linked with its distribution of zeroes: the loop-Hafnian is on average easy to estimate when $w$ lies in domains in the complex plane where the loop-Hafnian is on average non-zero, and #P-hard to estimate in domains where on average it contains zeroes. This is motivated by a similar transition in the matching polynomial over non-negative matrices. The zeroes of $g(z; \tilde{A})$ for some non-negative real matrix $\tilde{A}$ are distributed along the negative-real axis. The result is that multiplicative-error approximation to $g(z; \tilde{A})$ is efficient everywhere else on the complex plane, except on the negative-real axis where the problem is #P-hard for all negative real $z$ that is upper bounded by a value related to the degree of the graph$^{23,25,26}$.

The numerical evidence based on existing methods may then suggest a possible bound of $\tilde{w} = 1$ for Conjecture 4. For at least a fraction $\delta = \frac{1}{2}$ of $\tilde{X}$, when $|w| \le 1$, the polynomial $g(z = \frac{1}{w}; \tilde{X})$ would contain at least one zero.

Thus estimation of the loop-Hafnian escapes the applicability of the Taylor approximation method. At $w = 1$, the displacement to squeezing ratio is calculated to be $\frac{|\beta|^2}{\sinh^2(r)} = 5.9$, at which the loop-Hafnian also escapes the applicability of the k-th order approximation method proposed in ref. 21.

For Conjecture 5, the physical interpretation of Eq. (24) is that a non-negligible fraction of outcomes should occur with probability larger than that of a uniform distribution, and this fraction does not decay super-polynomially with $N$. This property is absent in the limiting case of $|w| = +\infty$, which is when squeezing is reduced to zero and the input only consists of coherent states. In this limit, the probability distribution $\tilde{p}_n$ is given by $|\mathrm{Haf}(0, X1_N)|^2 = \prod_{i=1}^N |\sum_{j=1}^K X_{ij}|^2$ for i.i.d. Gaussian random variables $X_{ij}$ whose magnitudes squared follow the exponential distribution. This results in a counter-example to anti-concentration, where only an exponentially small fraction of outcomes occurring with higher than uniform probability.

When $w$ is finite, the proof or disproof of Conjecture 5 becomes more difficult. When $|w|$ is small, we expect the loop-Hafnian distribution to be at least as anti-concentrated as the Hafnian. In the Supplemental Materials, we present numerical evidence to support this up to $|w| = 0.4$, which corresponds to having equal mean photon number from squeezing and displacement: $\frac{|\beta|^2}{\sinh^2(r)} = 1$SI, though the conjecture for higher but finite $|w|$ values cannot be ruled out.

## Discussion

In this work, we have studied the change of computational complexity when a D-GBS problem becomes increasingly displaced. Exact solution by an efficient classical algorithm is not possible, because estimating the loop-Hafnian is worst-case #P-hard. With some displacement, we argue that there is a regime where even additive-precision sampling is not possible for any efficient classical algorithm, when the parameter $w$, which characterises the relative importance of displacement to that of squeezing, is small.

As $|w|$ increases, we observe a transition in the complexity of D-GBS. Beyond $|w| = 0.4$, that is, as the ratio between the average photon number from displacement and from squeezing exceeds $|\beta|^2/\sinh^2(r) = 1$, we note that the outcomes become increasingly concentrated in a small subset, which may enable classical spoofing strategies. Beyond $|w| = 1$ or $|\beta|^2/\sinh^2(r) = 5.9$, we observe a transition in the average-case complexity of the loop-Hafnian as well. In this regime, a Taylor approximation algorithm has increasing success in estimating the sampling probability of a random outcome, until it can do so for majority of the outcomes. This invalidates the proof of hardness based on Polynomial Hierarchy collapse, and the algorithm could be used to construct an efficient classical sampler. The quoted $|w|$ values are numerical estimates based on existing methods, subject to change if tighter bounds are found on the transition in Conjectures 4 and 5. The complexity transition could also depend on other factors such as mode number and number of squeezers. Eventually, when the displacement is sufficiently high such that the condition in Theorem 1 or 2 is satisfied, all outcome probabilities can be approximated to arbitrary multiplicative-error in quasi-polynomial time.

An efficient approximate solution is also possible when the D-GBS problem encodes a non-negative graph. This constrains the types of graph problems that a D-GBS scheme can solve while maintaining an exponential quantum advantage.

One important hurdle in photonic quantum computing is overcoming loss. Tensor network methods exploiting photon loss have had great success in simulating large-scale GBS demonstrations$^{28}$. In tensor network methods, displacement is a local operation and would not change the runtime complexity. Therefore, by adding displacement without reducing their current levels of loss, experimental efforts on GBS will neither increase nor decrease the runtime separation between quantum and classical methods.

However, if loss levels can be sufficiently reduced below the simulability threshold, tensor network methods may lose out to classical algorithms that simulate the D-GBS process by loop-Hafnian calculations$^{12,40}$. The algorithm runtime of the latter scale exponentially with the number of photons. As long as optimal classical algorithms still rely on loop-Hafnian
calculations, D-GBS does present a scheme where the quantum-classical separation is increased by a factor exponential in the number of photons added by displacement.

It is also possible to simulate photon loss in GBS schemes by displacement. To sample from a mixed Gaussian state is equivalent to sampling from a pure Gaussian state with smaller mean photon number but with a random displacement<sup>28,40</sup>. It is therefore an interesting question for future work whether our framework of D-GBS can be used to study the complexity of lossy, un-displaced GBS. In the Supplemental Materials, we consider the link between the amount of loss required and the w parameter which indicates the displacement/squeezing ratio and therefore the suitability of the Taylor approximation method to simulating lossy GBS experiments SI.

Viewing displacement by a coherent state as a classical resource, our work shows that a quantum device can tolerate a certain amount of classicality without destroying its computational advantage. More importantly, since coherent states are readily generated by lasers, displacement presents a relatively easy way of scaling up the photon number—a current bottleneck in GBS experiments.

Beyond the ongoing goal of demonstrating quantum advantage, extending GBS to include displacement may support progress towards practical quantum advantage—solving classically intractable problems with real-world relevance. Displacement enhances the expressivity of graph-based GBS applications<sup>19,44-48</sup> by enabling the encoding of loops, which models self-interaction at graph vertices. Our results suggest that such applications must go beyond non-negative graphs for exponential quantum advantage. Displacement also enables chemistry applications such as molecular vibronic spectra simulation<sup>14-18</sup>. While the quantum advantage of this specific problem is under active research<sup>55</sup>, extending the Gaussian operations to include displacement unlocks a larger space to explore for complexity-theoretically grounded quantum advantage in quantum optics.

## Methods

### Photon statistics of D-GBS

The B matrix and y vector in Eq. (1) are given by refs. 10,22:

$$ \begin{pmatrix} \mathbf{B} & 0 \\ 0 & \mathbf{B}^* \end{pmatrix} = \begin{pmatrix} 0 & \mathbb{1}_M \\ \mathbb{1}_M & 0 \end{pmatrix} (\mathbb{1}_{2M} - \mathbf{\Sigma}_Q^{-1}), \quad (25) $$

$$ (y^T \quad y^\dagger) = \mu^\dagger \Sigma_Q^{-1}. \quad (26) $$

For an M-mode pure state constructed by applying an interferometer transformation, U on displaced squeezed states, $|\psi\rangle = \tilde{U}(U) \bigoplus_{j=1}^M \hat{D}_j(\beta_j) \hat{S}_j(r_j)|0\rangle$, the pair (B, y) is given by refs. 10,22 (detailed derivations also presented in the Supplemental MaterialsSI):

$$ \mathbf{B} = U \bigoplus_{j=1}^{M} \tanh(r_j) U^T, \quad (27) $$

$$ \gamma = U\beta^* - U\left(\bigoplus_{j=1}^{M} \tanh(r_j)\right)\beta. \quad (28) $$

In a general D-GBS problem, the matrix $\tilde{B}$, defined in Eq. (4) for easier mapping to the matching polynomial, is the B matrix rescaled by the an outer product of the y vector with itself. By Eqs. (27) and (28), we see that the magnitude $|\tilde{B}_{ij}|$ increases (decreases) monotonically with $|r_j|$ ($|\beta_j|$) and thus represents a ratio of squeezing to displacement.

This ratio is further simplified to the w parameter in Eq. (11) by Eqs. (8) and (9) in the Uniform D-GBS model. It is straightforward to derive Eqs. (8) and (9) from Eqs. (27) and (28), if we make use of the following requirements:

$$ r_j = \begin{cases} r_j, & \text{for } 1 \le j \le K, \\ 0, & \text{for } j > K; \end{cases} \quad (29) $$

$$ \beta_j = \begin{cases} \beta_j, & \text{for } 1 \le j \le K, \\ 0, & \text{for } j > K. \end{cases} \quad (30) $$

The photon number distribution of a single-mode displaced squeezed state, $|\psi_j\rangle = \hat{D}_j(\beta_j) \hat{S}_j(r_j)|0\rangle$, is given by refs. 56,57:

$$ p^{(j)}(n) = e^{-|\beta|^2 + j(\beta^2) \tanh^2(r)} \frac{\tanh^4(r)}{n! 2^n \cosh(r)} \left| H_n \left( \frac{\beta^* \cosh(r) - \beta \sinh(r)}{-i\sqrt{\sinh(2r)}} \right) \right|^2, \quad (31) $$

where $H_n$ is the n-th Hermite polynomial in the physicist's definition and the variable inside the Hermite polynomial can be rewritten in terms of the w parameter from Eq. (11). The probability of measuring N photons across K identically displaced squeezed states is then given by a discrete convolution over the K modes, with vacuum probability given by

$$ p_0 = \frac{1}{\cosh^K(r)} e^{-K|\beta|^2 + K \tanh(r) j(\beta^2)} \quad (32) $$

which derives the post-selected probability, $\tilde{p}_n$, in Equation (22).

### Taylor approximation to the loop-Hafnian

Here, we give an overview of the Taylor approximation method used in the main text for estimating diagonally dominant and non-negative loop-Hafnians, with a derivation of its runtime complexity. Detailed explanation can also be found in refs. 25,30. The scheme works for any general polynomial of degree d:

$$ g(z) = a_0 + a_1 z + \dots + a_d z^d. \quad (33) $$

The matching polynomial, $g(z; \tilde{A}) = \text{Haf}(z\tilde{A}, \mathbf{1}_N)$, on an $N \times N$ matrix $\tilde{A}$, has degree $d = \lfloor \frac{N}{2} \rfloor$.

If the roots of $g(z)$ are $\zeta_1, \dots, \zeta_d \in \mathbb{C}$, then $g(z)$ can be rewritten as

$$ g(z) = a_d(z - \zeta_1) \dots (z - \zeta_d). \quad (34) $$

If there is a disc in complex plane with radius R where the polynomial is non-zero, i.e. $|\zeta_i| > R$ for all $1 \le i \le d$, then we can define a function f to be the logarithm of $g(z)$:

$$ f(z) = \ln g(z) = \ln(a_d) + \ln(z - \zeta_1) + \dots + \ln(z - \zeta_d). \quad (35) $$

Function f can be Taylor expanded as

$$ f(z) = \sum_{j=1}^{\infty} \frac{z^j}{j!} f^{(j)}(0) = f(0) - \sum_{j=1}^{\infty} \frac{t_j z^j}{j}, \quad (36) $$

where $t_j = \sum_{i=1}^d \zeta_i^{-j}$.

If we truncate the Taylor expansion at degree m:

$$ T_m(z) = f(0) - \sum_{j=1}^{m} \frac{t_j z^j}{j}, \quad (37) $$

then we can approximate $f(z)$ to additive error:

$$ |f(z) - T_m(z)| = \left| \sum_{j=m+1}^{\infty} \frac{t_j z^j}{j} \right| \le \frac{1}{m+1} \sum_{j=m+1}^{\infty} |t_j z_j|. \quad (38) $$

Since the roots of the polynomial lie outside radius R in complex plane, $|\zeta_i| > R, \forall j \in [1, d]$, the absolute value of the j-th inverse power sum is bounded by $|t_j| = |\sum_{i=1}^d \zeta_i^{-j}| \le d/R^l$. The additive error is therefore bounded by the
sum of a geometric progression,

$$|f(z) - T_m(z)| = \frac{d}{m+1} \frac{q^{m+1}}{1-q}, \quad (39)$$

where the ratio $q$ is given by $q = |z|/R < 1$.

Given any error bound $\epsilon > 0$, we can choose $m = C(\ln d - \ln \epsilon)$, where the constant $C$ can be chosen as $C = (1-q)^{-1}$. Then

$$\frac{d}{m+1} \frac{q^{m+1}}{1-q} < dq^m = d \left(\frac{d}{\epsilon}\right)^{C \ln q} d \left(\frac{d}{\epsilon}\right)^{-1} = \epsilon. \quad (40)$$

To derive the runtime complexity of the algorithm for the loop-Hafnian, we make use of the loop-Hafnian expansion formula derived in refs. 10,21. For some $N \times N$ matrix $\tilde{A}$, one can express $\text{IHaf}(\tilde{A}, \mathbf{1}_N)$ as

$$\text{IHaf}(\tilde{A}, \mathbf{1}_N) = \sum_{m=0}^{[N/2]} \sum_{\substack{S \subseteq \{1, \dots, N\} \\ |S| = 2m}} \text{Haf}(\tilde{A}_S), \quad (41)$$

where $S$ is a set of integer numbers with cardinality $2m < N$, and $\tilde{A}_S$ is a submatrix of $\tilde{A}$ constructed by keeping the $i$-th row and column of $\tilde{A}$ if and only if $i \in S$. If $S = \emptyset$, then $\text{IHaf}(\tilde{A}, \emptyset) = 1$. The matching polynomial can therefore be rewritten in a similar form,

$$g(z; \tilde{A}) = \sum_{m=0}^{[N/2]} \left( \sum_{\substack{S \subseteq \{1, \dots, N\} \\ |S| = 2m}} \text{Haf}(\tilde{A}_S) \right) z^m, \quad (42)$$

whose $m$-th order derivative is:

$$g^{(m)}(0; \tilde{A}) = m! \sum_{\substack{S \subseteq \{1, \dots, N\} \\ |S| = 2m}} \text{Haf}(\tilde{A}_S). \quad (43)$$

The runtime of computing $T_m(z)$ scales with the coefficient of its leading term, $\frac{d}{dz^2} \ln g(0)$, which can be computed in polynomial number of calls to Eq. (43)$^{25}$. Since Eq. (43) can be computed by direct enumeration in time $N^{O(m)} = N^{O(\ln N - \ln \epsilon)}$, the function $T_m(z)$ can also be computed in quasi-polynomial time $N^{O(\ln N - \ln \epsilon)}$.

## Complexity proof of uniform D-GBS

Consider the classical algorithm $C$ with precision that satisfies Eq. (20). Then, by Stockmeyer's theorem, one can compute some $\tilde{q}_n$ that approximates $q_n$ to a multiplicative error, e.g. $|\tilde{q}_n - q_n| \le 2\epsilon q_n$, in BPP$^{NP,49}$. By using Markov's inequality,

$$\Pr_{n \in \Omega} \left( |q_n - \tilde{p}_n| \ge \frac{4\epsilon}{k|\Omega|} \right) \le \frac{k}{2}, \quad (44)$$

$$\Pr_{n \in \Omega} \left( q_n \ge \frac{2}{k|\Omega|} \right) \le \frac{k}{2}, \quad (45)$$

one derives that $\tilde{q}_n$ is also an approximation to $\tilde{p}_n$, but to an additive error,

$$|\tilde{q}_n - \tilde{p}_n| \le |\tilde{q}_n - q_n| + |q_n - \tilde{p}_n| \le \frac{8\epsilon}{k} \frac{1}{|\Omega|}, \quad (46)$$

with probability at least $1 - k$ over $n \in \Omega$.

Equation (46) becomes a multiplicative-error approximation to $\tilde{p}_n$, if the target distribution $\tilde{p}_n$ further satisfies an anti-concentration property

over $n$ in the form of Eq. (21). In this case, Eq. (46) can be rewritten as

$$|\tilde{q}_n - \tilde{p}_n| \le \frac{8\epsilon\alpha}{k} \tilde{p}_n \quad (47)$$

with probability at least $1 - k - \eta$ over $n \in \Omega$. If the probability function $\tilde{p}_n$ is on average #P-hard to approximate to within multiplicative error $\epsilon = O(\frac{\epsilon\alpha}{\delta})$ for at least a $(1 - \delta)$-fraction of outcomes $n \in \Omega$, then the approximate sampler $C$ would collapse the Polynomial Hierarchy to the third level.

In some complexity proofs for quantum random sampling schemes, the constants $\epsilon$ and $\delta$ in Conjecture 4 are given precise values$^{27}$, while we keep them in their generic form. To some extent, they can be traded off against each other and against the $\tilde{w}$ value. More importantly, they can also be traded off against $\epsilon$ and $\alpha$ in Eq. (47), which defines the precision of the adversarial classical algorithm and how anti-concentrated $\tilde{p}_n$ needs to be to rule it out. If the loop-Hafnian is average-case #P-hard for polynomial multiplicative error, $\epsilon = \text{poly}(N)$, then classical samplers with even polynomial TVD precision, $\epsilon = \text{poly}(N)$, are ruled out, and the $\alpha$ in Eq. (21) can be polynomial in $N$. In the Supplemental Materials, we provide a proof that multiplicative-error approximation of the loop-Hafnian is indeed worst-case #P-hard for polynomial error bounds, $\epsilon \le \text{poly}(N)$ SI. But we leave the question open on whether this could be conjectured for the average-case as well.

A similar trade-off exists for the anti-concentration conjecture. If the loop-Hafnian lacks anti-concentration entirely, e.g. if $\alpha = \Omega(\exp(N))$ in Eq. (24), then the adversarial classical algorithm needs to be exponentially precise in TVD $\epsilon = \frac{1}{\exp(N)}$ (Eq. (47)), in order to collapse the Polynomial Hierarchy. In other words, the quantum experiment needs to be exponentially close to the ground-truth distribution to guarantee an escape of classical simulability. On the other hand, if a stronger version of Conjecture 5 can hold true for some D-GBS variant, e.g. if $\alpha$ is constant in Eq. (24), then any adversarial classical sampler with constant TVD precision could be ruled out by assuming the Polynomial Hierarchy doesn't collapse.

A recent study on the Hafnian has also revealed a loss of anti-concentration in GBS if the number of initially-squeezed modes, $K$, scales too slowly with total mode number, $M^{58}$. A similar transition could also exist in D-GBS, but if the anti-concentration $K$ scale more quickly than $K = \sqrt{M}$, the output distribution may no longer 'hide' a Gaussian distribution of matrices. The average-case complexity of the loop-Hafnian needs then be studied over a different distribution of matrices. However, for a sufficiently random distribution, we still expect equivalence between the worst-case and the average-case complexity of the loop-Hafnian when displacement is small.

## Towards experimental implementations

Displacement of a single mode can be experimentally implemented by interfering the single-mode quantum state with a strong local oscillator on a highly asymmetric beamsplitter$^{59,60}$. The proposed Uniform D-GBS scheme requires identical displacement operations on all non-vacuum input modes, so that we can separate the $w$ parameter from the Gaussian submatrices of the Haar random unitary. These displacement operations can be implemented by an extra $M$ beamsplitters and phase-locked coherent states with individual control over their reflectivities.

Identical displacement on all populated input modes is an assumption made for mathematical convenience and is not a necessary condition for the D-GBS problem to be hard. If this requirement is dropped, ref. 21 showed that a single-mode coherent state is sufficient to apply displacements to a multi-mode state provided both are inputted into a fully connected interferometer, in which every input mode is able to interfere with every other input mode. Alternatively, displacement with a single coherent state may also be achieved by seeding the non-linear process that generates the squeezed vacuum state with a coherent state in the target wavelength$^{61,62}$. In either case, the action of the interferometer applies displacements to each mode depending on its transfer matrix, which in general can be non-zero for every mode. This then enables the same rewriting of the loop-Hafnian as the
matching polynomial, albeit on a different distribution of matrices from the Uniform D-GBS scheme.

In a similar vein, the requirement for the number of non-vacuum input modes to be upper bounded by $K \le \sqrt{M}$ was also made for mathematical convenience and to target a Gaussian distribution in the loop-Hafnian matrices. Increasing $K$ beyond this bound could be experimentally beneficial when the squeezing operations are marred by photon loss$^{28}$. Provided the new distribution is sufficiently random and is absent of special structures, we expect a similar reduction of the average-case to worst-case complexity of estimating the loop-Hafnian.

In our paper, we discussed two special matrix structures that would allow efficient estimation of the loop-Hafnian: diagonal-dominance and non-negativity. Another example is the sparsity. An interferometer with only limited connectivity, instead of the full connectivity that quantum advantage would stipulate, reduces the loop-Hafnian calculations to sparser matrices that allow for more efficient classical algorithms$^{13,63,64}$.

Applying displacement via the interferometer couples the covariance matrix $\Sigma$ and the mean vector $\mu$ defining the output state. To generate arbitrary displacement independent of the interferometer, one needs to revert back to the use of $M$ additional individually controlled beamsplitters. The ability of individually controlling squeezing, displacement and interferometer enables encoding of any arbitrary complex symmetric matrix into the D-GBS scheme. Therefore, if there is a specific distribution of matrices, on which stronger evidence is found for average-case complexity and the anti-concentration of the loop-Hafnian, one can design a quantum advantage experiment with D-GBS from this distribution.

Our complexity proof required the total average photon number to be upper bounded by square root of the number of modes, $\bar{N} \le \sqrt{M}$, such that the photon statistics remain collision-less. If this requirement were relaxed, collisions would introduce repetitions to the rows and columns of $fd(B_m, \gamma_n)$, which increases the dimensions of the matrix without increasing its rank. However, with a moderate amount of collisions, there is no evidence that the sampling problem becomes efficient to solve. The collision-less condition is also required in regular GBS proposals$^{11}$, although recent work suggests it may not be required for the hardness of Fock-basis Boson Sampling$^{65}$. We leave as open question for future work the implications of relaxing the collision-less requirement for D-GBS.

## Data availability

The numerical data underlying the results presented in this paper are available from the corresponding author on reasonable request.

Received: 8 November 2024; Accepted: 9 June 2025; Published online: 09 July 2025

## References

1. Zhong, H.-S. et al. Quantum computational advantage using photons. *Science* **370**, 1460 (2020).

2. Zhong, H.-S. et al. Phase-programmable Gaussian boson sampling using stimulated squeezed light. *Phys. Rev. Lett.* **127**, 180502 (2021).

3. Deng, Y.-H. et al. Gaussian boson sampling with pseudo-photon-number-resolving detectors and quantum computational advantage. *Phys. Rev. Lett.* **131**, 150601 (2023).

4. Madsen, L. S. et al. Quantum computational advantage with a programmable photonic processor. *Nature* **606**, 75 (2022).

5. Arute, F. et al. Quantum supremacy using a programmable superconducting processor. *Nature* **574**, 505 (2019).

6. Wu, Y. et al. Strong quantum computational advantage using a superconducting quantum processor. *Phys. Rev. Lett.* **127**, 180501 (2021).

7. Zhu, Q. et al. Quantum computational advantage via 60-qubit 24-cycle random circuit sampling. *Sci. Bull.* **67**, 240 (2022).

8. Aaronson, S. & Arkhipov, A. The computational complexity of linear optics, in *Proceedings of the Forty-Third Annual ACM Symposium on*

*Theory of Computing*, STOC '11 pp. 333-342 (Association for Computing Machinery, 2011).

9. Hamilton, C. S. et al. Gaussian boson sampling. *Phys. Rev. Lett.* **119**, 170501 (2017).

10. Kruse, R. et al. Detailed study of Gaussian boson sampling. *Phys. Rev. A* **100**, 032326 (2019).

11. Deshpande, A. et al. Quantum computational advantage via high-dimensional Gaussian boson sampling. *Sci. Adv.* **8**, eabi7894 (2022).

12. Bulmer, J. F. F. et al. The boundary for quantum advantage in Gaussian boson sampling. *Sci. Adv.* **8**, eabl9236 (2022).

13. Oh, C., Lim, Y., Fefferman, B. & Jiang, L. Classical simulation of boson sampling based on graph structure. *Phys. Rev. Lett.* **128**, 190501 (2022).

14. Huh, J., Guerreschi, G. G., Peropadre, B., McClean, J. R. & Aspuru-Guzik, A. Boson sampling for molecular vibronic spectra. *Nat. Photonics* **9**, 615 (2015).

15. Huh, J. & Yung, M.-H. Vibronic boson sampling: generalized Gaussian boson sampling for molecular vibronic spectra at finite temperature. *Sci. Rep.* **7**, 1 (2017).

16. Clements, W. R. et al. Approximating vibronic spectroscopy with imperfect quantum optics. *J. Phys. B At. Mol. Opt. Phys.* **51**, 245503 (2018).

17. Quesada, N. Franck-Condon factors by counting perfect matchings of graphs with loops. *J. Chem. Phys.* **150**, 164113 (2019).

18. Arrazola, J. M. et al. Quantum circuits with many photons on a programmable nanophotonic chip. *Nature* **591**, 54 (2021).

19. Schuld, M., Brádler, K., Israel, R., Su, D. & Gupt, B. Measuring the similarity of graphs with a Gaussian boson sampler. *Phys. Rev. A* **101**, 032314 (2020).

20. Rahimi-Keshari, S., Ralph, T. C. & Caves, C. M. Sufficient conditions for efficient classical simulation of quantum optics. *Phys. Rev. X* **6**, 021039 (2016).

21. Thekkadath, G. S. et al. Experimental demonstration of Gaussian boson sampling with displacement. *PRX Quantum* **3**, 020336 (2022).

22. Quesada, N. et al. Simulating realistic non-Gaussian state preparation. *Phys. Rev. A* **100**, 022341 (2019).

23. Heilmann, O. J. & Lieb, E. H. Theory of monomer-dimer systems. *Commun. Math. Phys.* **25**, 190 (1972).

24. Jerrum, M. Two-dimensional monomer-dimer systems are computationally intractable. *J. Stat. Phys.* **48**, 121 (1987).

25. Barvinok, A. *Combinatorics and Complexity of Partition Functions* (Springer International Publishing, 2016).

26. Bezáková, I., Galanis, A., Goldberg, L. A. & Stefankovic, D. The complexity of approximating the matching polynomial in the complex plane. *ACM Trans. Comput. Theory* **13**, 37 (2021).

27. Bremner, M. J., Montanaro, A. & Shepherd, D. J. Average-case complexity versus approximate simulation of commuting quantum computations. *Phys. Rev. Lett.* **117**, 080501 (2016).

28. Oh, C., Liu, M., Alexeev, Y., Fefferman, B. & Jiang, L. Classical algorithm for simulating experimental Gaussian boson sampling. *Nat Phys* **20**, 1461–1468 (2024).

29. Björklund, A., Gupt, B. & Quesada, N. A faster hafnian formula for complex matrices and its benchmarking on a supercomputer. *ACM J. Exp. Algorithmics* **24**, 17 (2019).

30. Patel, V. & Regts, G. Deterministic polynomial-time approximation algorithms for partition functions and graph polynomials. *SIAM J. Comput.* **46**, 1893 (2017).

31. Bayati, M., Gamarnik, D., Katz, D., Nair, C. & Tetali, P. Simple deterministic approximation algorithms for counting matchings, in *STOC '07: Proceedings of the thirty-ninth annual ACM symposium on Theory of computing* pp.122–127 (Association for Computing Machinery, 2007).

32. Brádler, K., Dallaire-Demers, P.-L., Rebentrost, P., Su, D. & Weedbrook, C. Gaussian boson sampling for perfect matchings of arbitrary graphs. *Phys. Rev. A* **98**, 032310 (2018).
33. Cariolaro, G. & Pierobon, G. Bloch-Messiah reduction of Gaussian unitaries by Takagi factorization. *Phys. Rev. A* **94**, 062109 (2016).

34. Houde, M., McCutcheon, W. & Quesada, N. Matrix decompositions in quantum optics: Takagi/Autonne, Bloch-Messiah/Euler, Iwasawa, and Williamson, *Can. J. Phys.* **102**, 497-507 (2024).

35. Barvinok, A. Computing the permanent of (some) complex matrices. *Found. Comput. Math.* **16**, 329 (2016).

36. Barvinok, A. Approximating permanents and hafnians. *Discrete Analysis* **2**, 34 (2017).

37. Barvinok, A. Computing permanents of complex diagonally dominant matrices and tensors. *Isr. J. Math.* **232**, 931-945 (2019).

38. Kotecký, R. & Preiss, D. Cluster expansion for abstract polymer models. *Commun. Math. Phys.* **103**, 491 (1986).

39. Mann, R. L. & Minko, R. M. Algorithmic cluster expansions for quantum problems. *PRX Quantum* **5**, 010305 (2024).

40. Quesada, N. et al. Quadratic speed-up for simulating Gaussian boson sampling. *PRX Quantum* **3**, 010306 (2022).

41. Jerrum, M., Sinclair, A. & Vigoda, E. A polynomial-time approximation algorithm for the permanent of a matrix with nonnegative entries. *J. ACM* **51**, 671-697 (2004).

42. Rudelson, M., Samorodnitsky, A. & Zeitouni, O. Hafnians, perfect matchings and Gaussian matrices. *Ann. Probab.* **44**, 2858 (2016).

43. Oh, C., Fefferman, B., Jiang, L. & Quesada, N. Quantum-inspired classical algorithm for graph problems by Gaussian boson sampling. *PRX Quantum* **5**, 020341 (2024).

44. Arrazola, J. M. & Bromley, T. R. Using Gaussian boson sampling to find dense subgraphs. *Phys. Rev. Lett.* **121**, 030503 (2018).

45. Sempere-Llagostera, S., Patel, R. B., Walmsley, I. A. & Kolthammer, W. S. Experimentally finding dense subgraphs using a time-bin encoded Gaussian boson sampling device. *Phys. Rev. X* **12**, 031045 (2022).

46. Deng, Y.-H. et al. Solving graph problems using Gaussian boson sampling. *Phys. Rev. Lett.* **130**, 190601 (2023).

47. Banchi, L., Fingerhuth, M., Babej, T., Ing, C. & Arrazola, J. M. Molecular docking with Gaussian boson sampling. *Sci. Adv.* **6**, eaax1950 (2020).

48. Yu, S. et al. A universal programmable Gaussian boson sampler for drug discovery. *Nat. Comput. Sci.* **3**, 839-848 (2023).

49. Stockmeyer, L. The complexity of approximate counting, in *Proceedings of the Fifteenth Annual ACM Symposium on Theory of Computing*, STOC '83 pp. 118-126 (Association for Computing Machinery, 1983).

50. Hangleiter, D. & Eisert, J. Computational advantage of quantum random sampling. *Rev. Mod. Phys.* **95**, 035001 (2023).

51. Toda, S. PP is as hard as the polynomial-time hierarchy. *SIAM J. Comput.* **20**, 865 (1991).

52. Bouland, A., Fefferman, B., Landau, Z. & Liu, Y. Noise and the frontier of quantum supremacy, in *2021 IEEE 62nd Annual Symposium on Foundations of Computer Science (FOCS)* pp. 1308-1317 (2022).

53. Movassagh, R. The hardness of random quantum circuits. *Nat. Phys.* **19**, 1719 (2023).

54. Lipton, R. J. New directions in testing, in *Distributed computing and cryptography*, Vol. 2, edited by Feigenbaum, J. and Merritt, M. J. pp. 191-202 (American Mathematical Society, 1991).

55. Oh, C., Lim, Y., Wong, Y., Fefferman, B. & Jiang, L. Quantum-inspired classical algorithms for molecular vibronic spectra. *Nat. Phys.* **20**, 225 (2024).

56. Wang, S., Fan, H.-Y. & Hu, L.-Y. Photon-number distributions of non-Gaussian states generated by photon subtraction and addition. *J. Opt. Soc. Am. B* **29**, 1020 (2012).

57. Fitzke, E., Niederschuh, F. & Walther, T. Simulating the photon statistics of multimode Gaussian states by automatic differentiation of generating functions. *APL Photonics* **8**, 026106 (2023).

58. Ehrenberg, A., Iosue, J. T., Deshpande, A., Hangleiter, D. & Gorshkov, A.V. Transition of Anticoncentration in Gaussian Boson Sampling. *Phys. Rev. Lett.* **134**, 140601 (2025).

59. Paris, M. G. Displacement operator by beam splitter. *Phys. Lett. A* **217**, 78 (1996).

60. Laiho, K., Avenhaus, M., Cassemiro, K. N. & Silberhorn, C. Direct probing of the Wigner function by time-multiplexed detection of photon statistics. *N. J. Phys.* **11**, 043012 (2009).

61. Roux, F. S. Stimulated parametric down-conversion for spatiotemporal metrology. *Phys. Rev. A* **104**, 043514 (2021).

62. Keller, A. et al. Reconstructing the full modal structure of photonic states by stimulated-emission tomography in the low-gain regime. *Phys. Rev. A* **106**, 063709 (2022).

63. Qi, H. et al. Efficient sampling from shallow Gaussian quantum-optical circuits with local interactions. *Phys. Rev. A* **105**, 052412 (2022).

64. Yang, T.-Y. & Wang, X.-B. Speeding up the classical simulation of Gaussian boson sampling with limited connectivity. *Sci. Rep.* **14**, 7680 (2024).

65. Bouland, A. et al. Complexity-theoretic foundations of BosonSampling with a linear number of modes https://arxiv.org/abs/2312.00286 (2023).

## Acknowledgements

The authors would like to thank Ryan L. Mann for informative discussions around the matching polynomial, complexity proof techniques and relevant literature. The authors would also like to thank Guillaume Thekkadath, Santiago Sempere-Llagostera, Nicolás Quesada, Nicolas Menicucci, Anthony Laing, Dara McCutcheon, Oliver Thomas, and Dominik Hangleiter for useful discussions. This work is supported by the Quantum Computing & Simulation Hub under the Engineering and Physical Sciences Research Council (EPSRC) UK Quantum Technologies Programme (project number: EP/T001062/1) and UK Research & Innovation (project number: MR/ W011794/1). Z.L. acknowledges partial funding support from ORCA Computing. N.R.S. is supported by the European Union's Horizon Europe research and innovation program under grant agreement No. 101102140, and acknowledges partial support from the Quantum Engineering Centre for Doctoral Training, EPSRC (project number: EP/SO23607/1).

## Author contributions

Z.L. developed the theory and performed the numerical calculations with contributions from N.R.S., J.F.F.B., R.B.P. and I.A.W. Z.L. and N.R.S. completed the complexity proof for Uniform D-GBS. N.R.S. detailed the link between D-GBS and lossy GBS with contributions from J.F.F.B. and Z.L. I.A.W., R.B.P. and J.F.F.B. conceived and managed the project. All authors contributed to the discussion and writing.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41534-025-01062-5.

**Correspondence** and requests for materials should be addressed to Zhenghao Li.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2025