Boson sampling

“The extended Church-Turing thesis states that ‘Any algorithmic process can be simulated efficiently using a probabilistic Turing machine’. This ad-hoc modification of the Church-Turing thesis should leave you feeling rather queasy.” - Nielsen and Chuang [1].

Introduced by Aaronson and Arkhipov [2], boson sampling presented a slight deviation from the general approach in quantum computation. Rather than presenting a theoretical model of universal quantum computation (i.e., a framework that enables quantum simulation of any arbitrary Hamiltonian [1]), boson sampling-based devices are instead an example of an intermediate quantum computer, designed to experimentally implement a computation that is thought to be intractable classically [3].

Boson sampling proposes the following quantum linear optics scheme. An array of single-photon sources is set up, designed to simultaneously emit single photon states into a multimode linear interferometer; the results are then generated by sampling from the probability of single photon measurements from the output of the linear interferometer.

For example, consider \(N\) single photon Fock states, \(\ket{\psi}=\ket{m_1,m_2,\dots,m_N}\), composed of \(b=\sum_i m_i\) photons, incident on an \(N\)-mode linear interferometer, which performs the following linear transformation of the input mode creation and annihilations operators:

\[\begin{split}&\hat{a}^\dagger_{out_k} = \sum_{j=0}^N U_{kj}\hat{a}^\dagger_{in_j}\\ &\hat{a}_{out_k} = \sum_{j=0}^N U_{jk}^\dagger\hat{a}_{in_j}\end{split}\]

Here, the unitary \(U\) completely describes the interferometer. Thus, the probability of detecting \(n_j\) photons at the \(j\text{th}\) output mode is given by


where \(W\) represents the action of \(U\) on the Fock basis (\(W\) is simply a homomorphism of \(U\)). The remarkable nature of the boson sampling problem to challenge the extended Church-Turing thesis lies in the fact that [2]

\[\left|\braketT{n_1,n_2,\dots,n_N}{W}{\psi}\right|^2 = \frac{\left|\text{Per}(U_{st})\right|^2}{m_1!m_2!\cdots m_N!n_1!n_2!\cdots n_N!}\]

i.e., the sampled single photon probability distribution is proportional to the permanent of \(U_{st}\), a submatrix of the interferometer unitary, dependent upon the input and output Fock states.


The permanent of a matrix, defined by

\[\text{Per}(A) = \sum_{\sigma=S_N}\prod_{i=1}^N A_{i\sigma(i)}\]

where \(S_N\) is the set of all permutations of \(N\) elements, is calculated in a similar fashion to the determinant, but unlike the determinant, the signatures of the permutations are not taken into account - every permutation is taken as a positive quantity.

In graph theory, the permanent calculates the number of perfect matchings in a bipartite graph with adjacency matrix \(A\).

Whilst the determinant can be calculated efficiently on classical computers, computing the permanent belongs to the computational complexity class of #P-Hard problems [4], which are strongly believed to be classically hard to calculate. (Surprisingly, even calculating the permanent in an approximate manner is a member of #P and intractable classically).

This implies that simulating boson sampling cannot be done efficiently on a classical computer, providing a potential challenge to the extended Church-Turing thesis, and demonstrating the power of (non-universal) quantum computation.

CV implementation

In quantum linear optics, the multimode linear interferometer is commonly decomposed into two-mode beamsplitters (BSgate) and single-mode phase shifters (Rgate) [5], allowing for a straightforward translation into a CV quantum circuit.

For example, in the case of a 4 mode interferometer, with arbitrary \(4\times 4\) unitary \(U\), the quantum optics circuit is given by


In the above, the detectors perform Fock state measurements, and the parameters of the beamsplitters and the rotation gates determines the unitary \(U\). Note that, in order to allow for arbitrary linear unitaries for \(m\) imput modes, we must have a minimum of \(m+1\) columns in the beamsplitter array [6].

Blackbird code

The boson sampling circuit displayed above, with randomly chosen rotation angles and beamsplitter parameters, can be implemented using the Blackbird quantum circuit language:

# prepare the input fock states
Fock(1) | q[0]
Fock(1) | q[1]
Vac     | q[2]
Fock(1) | q[3]

# rotation gates

# beamsplitter array
BSgate(0.7804, 0.8578)  | (q[0], q[1])
BSgate(0.06406, 0.5165) | (q[2], q[3])
BSgate(0.473, 0.1176)   | (q[1], q[2])
BSgate(0.563, 0.1517)   | (q[0], q[1])
BSgate(0.1323, 0.9946)  | (q[2], q[3])
BSgate(0.311, 0.3231)   | (q[1], q[2])
BSgate(0.4348, 0.0798)  | (q[0], q[1])
BSgate(0.4368, 0.6157)  | (q[2], q[3])


Currently, only the Fock backends (i.e., fock and tensorflow) support boson sampling, as the Gaussian backend cannot initialise Fock states.

If we wish to simulate Fock measurements, we can additionally include

Measure | q

after the beamsplitter array. After constructing the circuit and running the engine, the values of the Fock state measurements will be available within the samples attribute of the Result object returned by the engine.

Alternatively, you may omit the measurements, and extract the resulting Fock state probabilities directly via the state method all_fock_probs().


A fully functional Strawberry Fields simulation containing the above Blackbird code is included at examples/

For more details on running the above Blackbird code in Strawberry Fields, including calculations of how to determine the output Fock state probabilities using the matrix permanent and comparisons to the returned state, refer to the in-depth boson sampling tutorial.