Gaussian Boson Sampling with Strawberry Fields
Run a Gaussian boson sampling circuit in Strawberry Fields: squeezed states, beamsplitters, and the connection to molecular vibronic spectra and graph similarity.
What Is Gaussian Boson Sampling?
Gaussian boson sampling (GBS) is a model of photonic quantum computation designed to produce output distributions that are classically hard to simulate. The protocol starts with a set of optical modes prepared in squeezed vacuum states, passes them through a linear optical interferometer (a network of beamsplitters and phase shifters), and measures the number of photons at each output mode.
The classical hardness comes from the mathematics of the output. The probability of observing a particular photon-number pattern at the output is related to the permanent or hafnian of a matrix derived from the interferometer. Computing the matrix permanent is #P-hard: no polynomial-time classical algorithm is known, and the best known algorithm is exponential.
Unlike qubit-based quantum computing, GBS does not require universal quantum gates or error correction to demonstrate a computational advantage. It simply requires preparing good squeezed states and a low-loss interferometer.
Photonic Basics in Strawberry Fields
Strawberry Fields (by Xanadu) describes quantum optics computations using continuous-variable (CV) quantum mechanics. The key concepts are:
Modes: Each optical mode is like a qubit, but its state is described over an infinite-dimensional Hilbert space (photon number states 0, 1, 2, …).
Squeezed states: A squeezed vacuum state S(r)|0> has reduced uncertainty in one quadrature at the cost of increased uncertainty in the other. The squeezing parameter r controls how much. Larger r means more photons on average and a higher probability of multi-photon events at the output.
Beamsplitter: The standard two-mode linear optical element. A beamsplitter BS(theta, phi) mixes two modes according to a 2x2 unitary. Chains of beamsplitters implement arbitrary interferometers.
Fock measurement: Measuring the photon number in a mode collapses the state to a Fock state |n> and returns the integer n.
The Strawberry Fields Programming Model
A Strawberry Fields program uses three objects:
sf.Program(n_modes): defines the circuitsf.Engine(backend, backend_options): selects the simulator- Operations from
strawberryfields.ops: gates and measurements
The "gaussian" backend simulates Gaussian states efficiently in polynomial time, but can only compute probabilities, not sample photon patterns. The "fock" backend simulates the full Fock space but scales exponentially. For GBS, the Gaussian backend is used to compute output probabilities; sampling from real hardware or the Fock backend gives actual photon patterns.
pip install strawberryfields
A Small GBS Circuit
This example uses 4 modes with squeezing parameter r=1.0 and a fixed interferometer.
import strawberryfields as sf
from strawberryfields.ops import Sgate, BSgate, MeasureFock
import numpy as np
n_modes = 4
r = 1.0 # squeezing parameter
prog = sf.Program(n_modes)
with prog.context as q:
# Prepare all modes in squeezed vacuum states
for i in range(n_modes):
Sgate(r) | q[i]
# Interferometer: a fixed sequence of beamsplitters
# A real GBS device uses a Haar-random interferometer
BSgate(np.pi / 4, 0.0) | (q[0], q[1])
BSgate(np.pi / 3, 0.5) | (q[2], q[3])
BSgate(np.pi / 5, 0.2) | (q[1], q[2])
BSgate(np.pi / 4, 0.0) | (q[0], q[1])
BSgate(np.pi / 4, 0.0) | (q[2], q[3])
# Measure photon number at each output mode
MeasureFock() | q
# Run on the Gaussian simulator
eng = sf.Engine("gaussian")
result = eng.run(prog)
print("State object:", result.state)
print("Mean photon number per mode:", result.state.mean_photon(0), result.state.mean_photon(1))
Note: MeasureFock() in the Gaussian backend does not produce actual samples; it computes the marginal photon-number distribution. To generate samples, use the "fock" backend (small mode number) or Xanadu’s cloud hardware.
Sampling Photon Patterns
To generate actual photon-number samples from the circuit, use the Fock backend with a truncation dimension:
import strawberryfields as sf
from strawberryfields.ops import Sgate, BSgate, MeasureFock
import numpy as np
n_modes = 4
r = 0.5 # smaller r for Fock backend tractability
prog = sf.Program(n_modes)
with prog.context as q:
for i in range(n_modes):
Sgate(r) | q[i]
BSgate(np.pi / 4, 0.0) | (q[0], q[1])
BSgate(np.pi / 4, 0.0) | (q[2], q[3])
BSgate(np.pi / 4, 0.0) | (q[1], q[2])
MeasureFock() | q
eng = sf.Engine("fock", backend_options={"cutoff_dim": 7})
results = []
for _ in range(200):
result = eng.run(prog)
sample = result.samples[0] # photon counts per mode
results.append(sample)
# Show the 10 most common photon patterns
from collections import Counter
pattern_counts = Counter(tuple(s) for s in results)
for pattern, freq in pattern_counts.most_common(10):
total_photons = sum(pattern)
print(f" {pattern} n_photons={total_photons} freq={freq}")
The most frequent patterns tend to have low total photon counts (vacuum or single-photon events). Patterns with 3 or more total photons are rarer but carry the computational weight of GBS.
Application 1: Molecular Vibronic Spectra
One of the most striking natural applications of GBS is computing the vibronic spectrum of a molecule. When a molecule undergoes an electronic transition, its nuclear coordinates shift to a new equilibrium geometry. The resulting vibrational states are related by a Bogoliubov transformation, which is exactly what a GBS circuit computes.
The Franck-Condon factors, which determine the intensity of each spectral line, are proportional to the squared overlap between vibrational states before and after the transition. These overlaps are hafnians of matrices derived from the Doktorov transformation.
Strawberry Fields provides a dedicated module for this:
from strawberryfields.apps import vibronic
# Duschinsky matrix and displacement vector for a model molecule
# (from quantum chemistry software such as Gaussian or Psi4)
U = np.array([[0.9950, 0.0998], [-0.0998, 0.9950]])
delta = np.array([0.1, 0.0])
T = 0 # temperature in Kelvin (0 = ground state)
# Compute the GBS encoding for this molecule
Ud, r, alpha, phi = vibronic.gbs_params(U, delta, T)
print("Squeezing parameters:", r)
print("Interferometer:", Ud)
The circuit parameters produced by gbs_params directly encode the molecular transition. Running the GBS device produces samples whose photon-number patterns correspond to the vibrational quanta excited, and their frequencies reproduce the vibronic spectrum.
Application 2: Graph Similarity
GBS also provides a natural encoding for graphs. Given an adjacency matrix A of an undirected graph, set the interferometer and squeezing parameters so that the GBS output distribution concentrates probability on dense subgraphs of A.
The mathematical connection: the probability of observing a photon pattern is related to the hafnian of a submatrix of the GBS covariance matrix. When the covariance encodes A, this hafnian equals the number of perfect matchings in the induced subgraph, which concentrates on dense parts of the graph.
from strawberryfields.apps import graph, sample
# A simple 4-node graph with adjacency matrix
A = np.array([
[0, 1, 0, 1],
[1, 0, 1, 1],
[0, 1, 0, 1],
[1, 1, 1, 0],
], dtype=float)
# Generate GBS samples from the graph encoding
n_samples = 500
samples = sample.sample(A, n_mean=4, n_samples=n_samples)
# Dense subgraph candidates are patterns with many photons
dense_candidates = [s for s in samples if sum(s) >= 3]
print(f"High-photon samples: {len(dense_candidates)} of {n_samples}")
# The most frequent high-photon patterns identify dense subgraphs
counts = Counter(tuple(s) for s in dense_candidates)
print("Top dense-subgraph candidates:")
for pattern, freq in counts.most_common(5):
nodes = [i for i, n in enumerate(pattern) if n > 0]
print(f" pattern={pattern} nodes={nodes} freq={freq}")
Graph similarity can then be computed by comparing the GBS output distributions of two graphs, using classical kernel methods on the resulting sample histograms.
The Borealis Quantum Advantage Experiment
In 2022, Xanadu demonstrated quantum computational advantage using a device called Borealis. Borealis is a photonic chip with 216 squeezed modes, a programmable loop-based interferometer, and photon-number-resolving detectors.
The task was GBS on a specific interferometer. Xanadu estimated that simulating the Borealis output distribution classically would require approximately 9000 years on the best available supercomputer. The device completed the same task in 36 microseconds.
Unlike earlier quantum advantage demonstrations (Google’s Sycamore, USTC’s Jiuzhang), Borealis was fully programmable: users could specify the interferometer settings, making the device more broadly applicable than a fixed-circuit experiment.
Limitations and Caveats
GBS’s classical hardness relies on several assumptions that break down in realistic devices.
Photon loss: Every optical component introduces some loss. Lost photons effectively reduce the squeezing and collapse the output toward a classical distribution. At high enough loss rates, efficient classical simulation becomes possible again.
Photon-number resolution: True GBS hardness requires detectors that can distinguish 0, 1, 2, 3, … photons. Threshold detectors, which only distinguish “zero” from “one or more,” reduce the problem to a different (and easier) variant.
Distinguishable photons: If photons are not perfectly identical (due to timing jitter or spectral impurity), the output distribution is a mixture of GBS distributions, which may be easier to simulate classically.
These constraints mean that a convincing quantum advantage from GBS requires very high optical circuit fidelity, low-loss components, and photon-number-resolving detection across all modes simultaneously.
Summary
| Concept | Detail |
|---|---|
| Core operation | Squeezed states through a linear interferometer |
| Output | Photon-number patterns per mode |
| Classical hardness | Matrix hafnian computation |
| Key application 1 | Molecular vibronic spectra (Franck-Condon factors) |
| Key application 2 | Dense subgraph detection via graph encoding |
| Xanadu hardware | Borealis (216 modes, 2022 advantage claim) |
| Main limitations | Photon loss, threshold vs. PNR detectors, photon distinguishability |
GBS sits at the intersection of quantum optics, complexity theory, and graph theory. Its applications to molecular simulation and combinatorial optimization are both concrete and theoretically grounded, making it one of the most compelling near-term photonic quantum computing targets.
Was this tutorial helpful?