Getting Started with Strawberry Fields

Strawberry Fields is Xanadu’s Python library for photonic quantum computing. Unlike qubit-based systems that work with discrete two-level states, photonic systems use continuous degrees of freedom: the amplitude and phase of light. This is called the continuous-variable (CV) model.

Continuous Variable vs Qubit Computing

In the qubit model, the basic unit is a two-state system and gates are rotations on the Bloch sphere. In the CV model:

• The basic unit is an optical mode (a beam of light)
• States are described by quadrature operators X (position-like) and P (momentum-like)
• Superposition and entanglement still exist, but over a continuous infinite-dimensional Hilbert space
• The vacuum state |0⟩ is the starting point: a mode with no photons

Photonic systems have natural advantages: they operate at room temperature, photons are the natural carrier of quantum information over fiber, and Gaussian operations are highly efficient to simulate classically.

Phase Space Intuition: The Wigner Function

Before diving into code, it helps to build intuition about how CV states look in phase space. In classical mechanics, a particle has a definite position x and momentum p. In quantum mechanics, the Heisenberg uncertainty principle forbids simultaneous precision in both. The Wigner function W(x, p) is the quantum analog of a classical phase-space probability distribution. It maps every point (x, p) to a quasi-probability value that encodes the full quantum state.

Key properties of the Wigner function:

• For a vacuum state, W(x, p) is a round Gaussian centered at the origin with width 1/2 in both directions. This width represents the vacuum fluctuations, the minimum uncertainty allowed by quantum mechanics.
• For a coherent state |alpha⟩, W(x, p) is the same round Gaussian but shifted to the point (Re(alpha), Im(alpha)). The shape is identical to vacuum; only the center moves.
• For a squeezed state, W(x, p) is an elongated Gaussian. One quadrature is compressed below the vacuum level while the other expands, keeping the total area constant (as required by the uncertainty principle).
• Negative values of W(x, p) signal genuinely nonclassical states. Gaussian states never go negative, but Fock states and other non-Gaussian states can.

Strawberry Fields provides a direct way to compute the Wigner function from the Fock backend using state.wigner(mode, xvec, pvec). Here is how to visualize vacuum, displaced, and squeezed states side by side:

import strawberryfields as sf
from strawberryfields import ops
import numpy as np
import matplotlib.pyplot as plt

# Phase space grid
xvec = np.linspace(-4, 4, 200)
pvec = np.linspace(-4, 4, 200)

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# --- Vacuum state ---
prog_vac = sf.Program(1)
eng_vac = sf.Engine("fock", backend_options={"cutoff_dim": 15})
with prog_vac.context as q:
    pass  # No gates: mode stays in vacuum
result_vac = eng_vac.run(prog_vac)
W_vac = result_vac.state.wigner(0, xvec, pvec)
axes[0].contourf(xvec, pvec, W_vac, levels=60, cmap="RdBu_r")
axes[0].set_title("Vacuum |0⟩")
axes[0].set_xlabel("x")
axes[0].set_ylabel("p")
axes[0].set_aspect("equal")

# --- Coherent (displaced) state ---
prog_coh = sf.Program(1)
eng_coh = sf.Engine("fock", backend_options={"cutoff_dim": 15})
with prog_coh.context as q:
    ops.Dgate(1.5, 0.5) | q[0]  # displacement with amplitude 1.5, phase 0.5 rad
result_coh = eng_coh.run(prog_coh)
W_coh = result_coh.state.wigner(0, xvec, pvec)
axes[1].contourf(xvec, pvec, W_coh, levels=60, cmap="RdBu_r")
axes[1].set_title("Coherent state (displaced)")
axes[1].set_xlabel("x")
axes[1].set_ylabel("p")
axes[1].set_aspect("equal")

# --- Squeezed state ---
prog_sq = sf.Program(1)
eng_sq = sf.Engine("fock", backend_options={"cutoff_dim": 15})
with prog_sq.context as q:
    ops.Sgate(1.0) | q[0]  # squeezing parameter r=1.0
result_sq = eng_sq.run(prog_sq)
W_sq = result_sq.state.wigner(0, xvec, pvec)
axes[2].contourf(xvec, pvec, W_sq, levels=60, cmap="RdBu_r")
axes[2].set_title("Squeezed state (r=1.0)")
axes[2].set_xlabel("x")
axes[2].set_ylabel("p")
axes[2].set_aspect("equal")

plt.tight_layout()
plt.savefig("wigner_comparison.png", dpi=150)
plt.show()

The vacuum plot shows a symmetric circle. The coherent state is the same circle shifted away from the origin. The squeezed state is an ellipse: narrow in x, wide in p (or the reverse, depending on the squeezing phase).

Installation

pip install strawberryfields

For cloud access to Xanadu’s Borealis hardware:

pip install strawberryfields xanadu-cloud-client

For computing Hafnians and verifying Gaussian boson sampling results:

pip install thewalrus

Core Concepts: Programs and Engines

In Strawberry Fields, a circuit is a Program. An Engine runs the program on a chosen backend. This separation keeps the circuit description backend-agnostic.

import strawberryfields as sf
from strawberryfields import ops

# A program on 2 modes (optical modes, not qubits)
prog = sf.Program(2)

The two available local backends are:

• Fock backend: represents the state as a tensor in the photon number (Fock) basis, truncated at cutoff_dim. It can handle arbitrary states including non-Gaussian ones, but memory grows as cutoff_dim^n where n is the number of modes.
• Gaussian backend: represents the state via a covariance matrix and mean vector. It handles only Gaussian states and Gaussian measurements, but scales polynomially in the number of modes.

Key Gates

Gate	Symbol	Effect
Displacement	D(alpha)	Shifts the state in phase space
Squeezing	S(r)	Compresses one quadrature, expands the other
Rotation	R(phi)	Rotates in phase space
Beamsplitter	BS(theta, phi)	Mixes two modes (like a 50/50 mirror)
Two-mode squeezing	S2(r)	Generates entanglement between two modes
Kerr gate	K(kappa)	Non-Gaussian, needed for universal CV computing

Gaussian gates (D, S, R, BS, S2) form an efficiently simulable subset. Adding any non-Gaussian gate (like K or cubic phase) makes the model universal but also harder to simulate.

Coherent States and Displacement

Coherent states |alpha⟩ are the closest quantum analog of classical laser fields. They are eigenstates of the annihilation operator: a|alpha⟩ = alpha|alpha⟩. Every laser produces light in a coherent state.

The displacement gate creates a coherent state from vacuum:

ops.Dgate(r, phi) | q[0]  # alpha = r * exp(i * phi)

Here r is the displacement amplitude and phi is the displacement phase. The complex number alpha = r * exp(i * phi) determines the location in phase space: Re(alpha) sets the X displacement, Im(alpha) sets the P displacement.

A coherent state has a Poisson photon number distribution:

P(n) = exp(-|alpha|^2) * |alpha|^(2n) / n!

The mean photon number is |alpha|^2 and the variance equals the mean (the defining property of a Poisson distribution). Let’s verify this:

import strawberryfields as sf
from strawberryfields import ops
import numpy as np
from scipy.stats import poisson

# Create a coherent state with alpha = 1.0 (real displacement)
prog = sf.Program(1)
eng = sf.Engine("fock", backend_options={"cutoff_dim": 20})

with prog.context as q:
    ops.Dgate(1.0, 0.0) | q[0]  # alpha = 1.0, mean photon number = |alpha|^2 = 1.0

result = eng.run(prog)
state = result.state

# Get Fock state probabilities from the simulation
probs_sf = state.all_fock_probs()

# Compare with the Poisson distribution (mean = |alpha|^2 = 1.0)
mean_n = 1.0  # |alpha|^2
probs_poisson = [poisson.pmf(n, mean_n) for n in range(len(probs_sf))]

print("n | SF prob   | Poisson prob")
print("-" * 35)
for n in range(10):
    print(f"{n} | {probs_sf[n]:.6f}  | {probs_poisson[n]:.6f}")

# The two columns match to high precision, confirming the Poisson distribution
# Mean photon number from SF:
mean_from_sf = state.mean_photon(0)[0]
print(f"\nMean photon number: {mean_from_sf:.4f} (expected: 1.0000)")

The Poisson distribution means that measuring the photon number of a coherent state yields a random result, most likely near the mean |alpha|^2, with fluctuations proportional to sqrt(|alpha|^2).

Example 1: Squeezed State with the Fock Backend

The Fock backend works in the photon number (Fock) basis. It is exact up to a chosen cutoff dimension and can represent arbitrary states including non-Gaussian ones.

import strawberryfields as sf
from strawberryfields import ops
import numpy as np

# Program on 1 mode, Fock backend with photon number cutoff of 10
prog = sf.Program(1)
eng = sf.Engine("fock", backend_options={"cutoff_dim": 10})

with prog.context as q:
    # Start in vacuum, apply squeezing with parameter r=0.5
    ops.Sgate(0.5) | q[0]

result = eng.run(prog)

# Get the state object
state = result.state

# Fock state probabilities: P(n photons)
probs = state.all_fock_probs()
print("Photon number probabilities:")
for n, p in enumerate(probs):
    if p > 0.01:
        print(f"  |{n}⟩: {p:.4f}")

# Squeezed states only have even photon number components
# |0⟩: ~0.8868, |2⟩: ~0.1049, |4⟩: ~0.0077

Squeezing compresses the uncertainty in one quadrature below the vacuum level while expanding the other. The resulting state has photon number components only at even n, which is a signature of single-mode squeezing. The physical reason: the squeezing operator S(r) = exp(r/2 * (a^2 - a†^2)) creates and destroys photons in pairs, so it can only change the photon number by 0, 2, 4, and so on. Starting from vacuum (n=0), only even photon numbers are reachable.

Homodyne Detection in Depth

Homodyne detection measures a single quadrature of the optical field. The measurement angle determines which quadrature is measured: angle=0 measures the X quadrature, angle=pi/2 measures the P quadrature, and intermediate angles measure linear combinations.

The measurement outcome is a continuous real number (unlike Fock measurement, which yields a non-negative integer). For the vacuum state, X measurement outcomes follow a Gaussian distribution N(0, 1/2), where the variance of 1/2 represents the vacuum noise floor.

Squeezing reduces the variance in one quadrature below this vacuum level. For squeezing parameter r, the X-quadrature variance becomes (1/2) * exp(-2r) and the P-quadrature variance becomes (1/2) * exp(2r). Let’s verify this by running many shots:

import strawberryfields as sf
from strawberryfields import ops
import numpy as np
import matplotlib.pyplot as plt

n_shots = 1000
r = 1.0  # squeezing parameter

# --- Vacuum state homodyne ---
samples_vacuum = []
for _ in range(n_shots):
    prog = sf.Program(1)
    eng = sf.Engine("fock", backend_options={"cutoff_dim": 15})
    with prog.context as q:
        ops.MeasureHomodyne(0) | q[0]  # measure X on vacuum
    result = eng.run(prog)
    samples_vacuum.append(result.samples[0, 0])

# --- Squeezed state homodyne (X quadrature, squeezed direction) ---
samples_squeezed = []
for _ in range(n_shots):
    prog = sf.Program(1)
    eng = sf.Engine("fock", backend_options={"cutoff_dim": 15})
    with prog.context as q:
        ops.Sgate(r) | q[0]
        ops.MeasureHomodyne(0) | q[0]  # measure X (squeezed direction)
    result = eng.run(prog)
    samples_squeezed.append(result.samples[0, 0])

# --- Squeezed state homodyne (P quadrature, anti-squeezed direction) ---
samples_antisqueezed = []
for _ in range(n_shots):
    prog = sf.Program(1)
    eng = sf.Engine("fock", backend_options={"cutoff_dim": 15})
    with prog.context as q:
        ops.Sgate(r) | q[0]
        ops.MeasureHomodyne(np.pi / 2) | q[0]  # measure P (anti-squeezed)
    result = eng.run(prog)
    samples_antisqueezed.append(result.samples[0, 0])

samples_vacuum = np.array(samples_vacuum)
samples_squeezed = np.array(samples_squeezed)
samples_antisqueezed = np.array(samples_antisqueezed)

# Compute variances
var_vacuum = np.var(samples_vacuum)
var_squeezed = np.var(samples_squeezed)
var_antisqueezed = np.var(samples_antisqueezed)

# Expected values
expected_vacuum = 0.5
expected_squeezed = 0.5 * np.exp(-2 * r)     # ~0.068
expected_antisqueezed = 0.5 * np.exp(2 * r)  # ~3.695

print(f"Vacuum X variance:        {var_vacuum:.3f} (expected: {expected_vacuum:.3f})")
print(f"Squeezed X variance:      {var_squeezed:.3f} (expected: {expected_squeezed:.3f})")
print(f"Anti-squeezed P variance: {var_antisqueezed:.3f} (expected: {expected_antisqueezed:.3f})")

# Plot histograms
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

axes[0].hist(samples_vacuum, bins=40, density=True, alpha=0.7, color="gray")
axes[0].set_title(f"Vacuum X (var={var_vacuum:.3f})")
axes[0].set_xlabel("X")

axes[1].hist(samples_squeezed, bins=40, density=True, alpha=0.7, color="blue")
axes[1].set_title(f"Squeezed X (var={var_squeezed:.3f})")
axes[1].set_xlabel("X")

axes[2].hist(samples_antisqueezed, bins=40, density=True, alpha=0.7, color="red")
axes[2].set_title(f"Anti-squeezed P (var={var_antisqueezed:.3f})")
axes[2].set_xlabel("P")

plt.tight_layout()
plt.savefig("homodyne_histograms.png", dpi=150)
plt.show()

The squeezed histogram is visibly narrower than vacuum, while the anti-squeezed histogram is visibly wider. With 1000 shots and r=1.0, the squeezed variance should be close to 0.068 (about 7x below vacuum noise). This sub-vacuum noise is the resource that makes CV quantum computing and CV quantum key distribution possible.

Example 2: Entanglement via Two-Mode Squeezing

The two-mode squeezed vacuum (TMSV) state is the CV analog of a Bell state. It is the primary entanglement resource in photonic quantum computing and quantum communication.

import strawberryfields as sf
from strawberryfields import ops
import numpy as np

prog = sf.Program(2)
eng = sf.Engine("gaussian")

with prog.context as q:
    # Two-mode squeezing creates entanglement between modes 0 and 1
    ops.S2gate(1.0) | (q[0], q[1])

result = eng.run(prog)
state = result.state

# The covariance matrix describes all quadrature correlations
# For 2 modes, it is 4x4: rows/columns are [x0, x1, p0, p1]
cov = state.cov()
print("Covariance matrix:")
print(np.round(cov, 4))
print()

# Compute the symplectic eigenvalues to check entanglement
# For a separable state, all symplectic eigenvalues >= 1
# For an entangled state, at least one < 1 (after partial transpose)
print("Symplectic eigenvalues:", state.symplectic_eigenvalues())

The entanglement in a TMSV state manifests as EPR-type correlations between the two modes. For the quadrature operators X0, X1, P0, P1, the TMSV state satisfies:

• Var(X0 - X1) = exp(-2r), which approaches zero for large r
• Var(P0 + P1) = exp(-2r), which approaches zero for large r

In the limit r → infinity, measuring X on mode 0 perfectly predicts X on mode 1 (they become identical), and measuring P on mode 0 perfectly predicts -P on mode 1. This is exactly the situation Einstein, Podolsky, and Rosen described in their 1935 paper. Let’s verify these correlations from the covariance matrix:

import strawberryfields as sf
from strawberryfields import ops
import numpy as np

for r in [0.5, 1.0, 1.5, 2.0]:
    prog = sf.Program(2)
    eng = sf.Engine("gaussian")

    with prog.context as q:
        ops.S2gate(r) | (q[0], q[1])

    result = eng.run(prog)
    cov = result.state.cov()

    # Strawberry Fields uses the ordering [x0, x1, p0, p1]
    # Var(X0 - X1) = Var(X0) + Var(X1) - 2*Cov(X0, X1)
    var_x_diff = cov[0, 0] + cov[1, 1] - 2 * cov[0, 1]

    # Var(P0 + P1) = Var(P0) + Var(P1) + 2*Cov(P0, P1)
    var_p_sum = cov[2, 2] + cov[3, 3] + 2 * cov[2, 3]

    expected = np.exp(-2 * r)

    print(f"r={r:.1f}: Var(X0-X1)={var_x_diff:.4f}, "
          f"Var(P0+P1)={var_p_sum:.4f}, expected={expected:.4f}")

As r increases, both variances decrease exponentially. At r=2.0, the correlations are strong enough that the outcomes of mode 0 and mode 1 are nearly perfectly correlated. This is the resource that enables CV quantum teleportation and CV quantum key distribution.

The Beamsplitter as a Photonic Hadamard Gate

The beamsplitter gate BSgate(theta, phi) is the fundamental two-mode operation in photonic circuits. It transforms the annihilation operators of two modes as:

a_out = cos(theta) * a_in  +  exp(i*phi) * sin(theta) * b_in
b_out = -exp(-i*phi) * sin(theta) * a_in  +  cos(theta) * b_in

When theta = pi/4 and phi = pi/2, this gives a 50/50 beamsplitter (equal transmission and reflection). The action on Fock states is particularly instructive: if one photon enters mode 0 and vacuum enters mode 1, the output is a superposition:

|1,0⟩ → (|1,0⟩ + i|0,1⟩) / sqrt(2)

This is the photonic equivalent of a Hadamard gate acting on a qubit encoded in dual-rail (one photon in two modes). Let’s verify this with the Fock backend:

import strawberryfields as sf
from strawberryfields import ops
import numpy as np

# Prepare |1, 0⟩ and apply a 50/50 beamsplitter
prog = sf.Program(2)
eng = sf.Engine("fock", backend_options={"cutoff_dim": 5})

with prog.context as q:
    ops.Fock(1) | q[0]   # one photon in mode 0
    ops.Vac() | q[1]     # vacuum in mode 1 (default, but explicit for clarity)
    ops.BSgate(np.pi / 4, np.pi / 2) | (q[0], q[1])  # 50/50 beamsplitter

result = eng.run(prog)
state = result.state

# Check Fock probabilities
p10 = state.fock_prob([1, 0])  # probability of |1,0⟩
p01 = state.fock_prob([0, 1])  # probability of |0,1⟩
p00 = state.fock_prob([0, 0])  # probability of |0,0⟩ (should be ~0)
p11 = state.fock_prob([1, 1])  # probability of |1,1⟩ (should be ~0)

print(f"P(1,0) = {p10:.4f}")  # should be ~0.5
print(f"P(0,1) = {p01:.4f}")  # should be ~0.5
print(f"P(0,0) = {p00:.4f}")  # should be ~0.0
print(f"P(1,1) = {p11:.4f}")  # should be ~0.0

The photon exits through one port or the other with 50/50 probability. This single-photon interference is the foundation of linear optical quantum computing (the KLM protocol). Notice that the total photon number is conserved: the beamsplitter cannot create or destroy photons, only redistribute them between modes.

Gaussian Simulation Efficiency

The Gaussian backend represents the state of n modes using:

• A 2n x 2n covariance matrix (quadrature correlations)
• A 2n-length mean vector (quadrature displacements)

Applying a Gaussian gate updates these matrices with O(n^2) matrix multiplications. This is enormously more efficient than the Fock backend, which stores a rank-n tensor of dimension cutoff_dim in each index. The memory requirement is cutoff_dim^n complex amplitudes.

For concrete numbers:

Modes	Fock (cutoff=10)	Gaussian
1	10 amplitudes	2x2 + 2 = 6 numbers
2	100 amplitudes	4x4 + 4 = 20 numbers
5	100,000 amplitudes	10x10 + 10 = 110 numbers
10	10^10 amplitudes	20x20 + 20 = 420 numbers
50	impossible	100x100 + 100 = 10,100 numbers

At 10 modes, the Fock backend needs 10 billion complex numbers (about 160 GB of memory). The Gaussian backend needs 420 real numbers. This exponential gap is why the Gaussian backend exists and why Gaussian states are considered “easy” from a computational perspective.

Let’s time both backends to see the difference in practice:

import strawberryfields as sf
from strawberryfields import ops
import numpy as np
import time

cutoff = 10

for n_modes in [1, 2, 3, 4, 5]:
    # Fock backend
    prog_fock = sf.Program(n_modes)
    with prog_fock.context as q:
        for i in range(n_modes):
            ops.Sgate(0.3) | q[i]
        for i in range(n_modes - 1):
            ops.BSgate(np.pi / 4, 0) | (q[i], q[i + 1])

    t0 = time.time()
    eng_fock = sf.Engine("fock", backend_options={"cutoff_dim": cutoff})
    result_fock = eng_fock.run(prog_fock)
    t_fock = time.time() - t0

    # Gaussian backend
    prog_gauss = sf.Program(n_modes)
    with prog_gauss.context as q:
        for i in range(n_modes):
            ops.Sgate(0.3) | q[i]
        for i in range(n_modes - 1):
            ops.BSgate(np.pi / 4, 0) | (q[i], q[i + 1])

    t0 = time.time()
    eng_gauss = sf.Engine("gaussian")
    result_gauss = eng_gauss.run(prog_gauss)
    t_gauss = time.time() - t0

    print(f"Modes={n_modes}: Fock={t_fock:.4f}s, Gaussian={t_gauss:.4f}s, "
          f"ratio={t_fock / max(t_gauss, 1e-9):.1f}x")

You will see the Fock backend time grow rapidly while the Gaussian backend stays nearly constant. Beyond 5 or 6 modes with cutoff=10, the Fock backend becomes impractical. For circuits that only use Gaussian operations, always prefer the Gaussian backend.

Measurements

Strawberry Fields supports several measurement types:

prog = sf.Program(2)
eng = sf.Engine("fock", backend_options={"cutoff_dim": 8})

with prog.context as q:
    ops.Sgate(0.5) | q[0]
    ops.BSgate(np.pi / 4, 0) | (q[0], q[1])

    # Homodyne detection: measures a quadrature (like X or P)
    ops.MeasureHomodyne(0) | q[0]   # angle=0 measures X quadrature

    # Fock measurement: measures photon number
    ops.MeasureFock() | q[1]

result = eng.run(prog, shots=5)
print("Homodyne samples (X quadrature):", result.samples[:, 0])
print("Fock samples (photon number):", result.samples[:, 1])

The available measurement operations are:

• MeasureHomodyne(phi): measures the quadrature at angle phi. Angle 0 gives X, angle pi/2 gives P. The result is a continuous real number.
• MeasureHeterodyne(): simultaneously measures both X and P, at the cost of added noise (1 unit of vacuum noise in each quadrature). The result is a complex number.
• MeasureFock(): measures the photon number. The result is a non-negative integer. This measurement is non-Gaussian and is only available in the Fock backend.

The Blackbird Language

Strawberry Fields circuits can also be written in Blackbird, a domain-specific language for photonic circuits. This is useful for saving and sharing circuits in a hardware-independent format:

prog = sf.Program(2)
with prog.context as q:
    ops.Sgate(0.6) | q[0]
    ops.BSgate(0.43, 0.1) | (q[0], q[1])
    ops.MeasureFock() | q[0]
    ops.MeasureFock() | q[1]

# Export to Blackbird
bb = prog.serialize()
print(bb)

Gaussian Boson Sampling

Gaussian Boson Sampling (GBS) is the task that Xanadu’s Borealis processor was designed for. The circuit is conceptually simple:

1. Prepare n single-mode squeezed states
2. Apply a random interferometer (a network of beamsplitters and phase shifters that implements a unitary transformation on the modes)
3. Measure the photon number at each output mode

Classically simulating the output probability distribution of this circuit requires computing the Hafnian of a matrix related to the interferometer. The Hafnian is similar to the permanent (which appears in standard boson sampling) but applies to symmetric matrices. Computing the Hafnian of an n x n matrix is a #P-hard problem, which is the computational hardness that makes GBS a candidate for quantum advantage.

GBS has practical applications beyond demonstrating quantum advantage:

• Graph optimization: GBS samples are biased toward dense subgraphs of a graph encoded in the interferometer
• Molecular simulation: vibronic spectra of molecules can be computed from GBS samples
• Point processes: GBS generates samples from a determinantal point process, useful in machine learning

Here is a complete 4-mode GBS circuit:

import strawberryfields as sf
from strawberryfields import ops
import numpy as np

# 4-mode GBS circuit
prog = sf.Program(4)
eng = sf.Engine("fock", backend_options={"cutoff_dim": 7})

with prog.context as q:
    # Step 1: Squeeze all modes
    for i in range(4):
        ops.Sgate(0.4) | q[i]

    # Step 2: Random interferometer (beamsplitters + phases)
    # Layer 1
    ops.BSgate(0.6, 0.1) | (q[0], q[1])
    ops.BSgate(0.3, 0.7) | (q[2], q[3])
    # Layer 2
    ops.BSgate(0.5, 0.2) | (q[1], q[2])
    # Layer 3
    ops.BSgate(0.4, 0.9) | (q[0], q[1])
    ops.BSgate(0.7, 0.3) | (q[2], q[3])

    # Phase rotations
    ops.Rgate(0.3) | q[0]
    ops.Rgate(1.1) | q[1]
    ops.Rgate(0.8) | q[2]
    ops.Rgate(2.0) | q[3]

result = eng.run(prog)
state = result.state

# Probability of specific output patterns
patterns = [[1, 1, 0, 0], [0, 1, 1, 0], [1, 0, 0, 1], [0, 0, 1, 1],
            [2, 0, 0, 0], [0, 0, 0, 2], [1, 1, 1, 1]]

print("Output pattern probabilities:")
for pattern in patterns:
    prob = state.fock_prob(pattern)
    if prob > 1e-6:
        print(f"  |{pattern}⟩: {prob:.6f}")

The key insight: these probabilities are proportional to the Hafnian of a submatrix of the circuit’s adjacency matrix. For a detection pattern S = (s0, s1, …, sn-1), the probability involves Haf(A_S), where A_S is the submatrix of A formed by repeating row/column i exactly s_i times. This is the connection between GBS and computational complexity.

The Walrus Library for Hafnians

The Walrus is Xanadu’s companion library for computing Hafnians, loop Hafnians, permanents, and other matrix functions that arise in quantum optics. It provides efficient C++ implementations of these #P-hard computations for small matrices.

The Hafnian of a 2n x 2n symmetric matrix A is defined as:

Haf(A) = sum over all perfect matchings M of prod_{(i,j) in M} A_ij

For example, a 4x4 matrix has 3 perfect matchings: {(0,1),(2,3)}, {(0,2),(1,3)}, {(0,3),(1,2)}.

from thewalrus import hafnian
import numpy as np

# A small 4x4 symmetric matrix
A = np.array([
    [0.0, 0.5, 0.3, 0.1],
    [0.5, 0.0, 0.2, 0.4],
    [0.3, 0.2, 0.0, 0.6],
    [0.1, 0.4, 0.6, 0.0]
])

haf = hafnian(A)
print(f"Hafnian of A: {haf:.6f}")

# Verify by hand: the 3 perfect matchings of a 4x4 matrix are:
# {(0,1),(2,3)}: A[0,1]*A[2,3] = 0.5 * 0.6 = 0.30
# {(0,2),(1,3)}: A[0,2]*A[1,3] = 0.3 * 0.4 = 0.12
# {(0,3),(1,2)}: A[0,3]*A[1,2] = 0.1 * 0.2 = 0.02
# Sum = 0.30 + 0.12 + 0.02 = 0.44
print(f"Manual computation: {0.30 + 0.12 + 0.02:.6f}")

In a GBS experiment, the probability of detecting photons in pattern S = (n_0, n_1, …, n_{m-1}) is:

P(S) = |Haf(A_S)|^2 / (sqrt(det(Q)) * n_0! * n_1! * … * n_{m-1}!)

where A_S is the submatrix of the adjacency matrix corresponding to the output pattern, and Q = sigma + I/2 is the Husimi covariance matrix. This formula is how you verify GBS circuits classically for small instances, and why large instances are intractable: the number of perfect matchings grows combinatorially.

Non-Gaussian Gates and Universality

Gaussian operations alone are not enough for universal quantum computing. They can be efficiently simulated classically (as we saw with the Gaussian backend), which means they cannot provide a computational advantage over classical computers for arbitrary tasks. To reach universality, you need at least one non-Gaussian element.

The Kerr gate K(kappa) applies a photon-number-dependent phase shift:

K(kappa)|n⟩ = exp(i * kappa * n * (n-1)) |n⟩

This means the phase applied to each Fock component depends quadratically on the photon number. A coherent state, which is a superposition of many Fock states, gets its components rotated by different amounts. The result is a distorted state whose Wigner function develops non-Gaussian features.

import strawberryfields as sf
from strawberryfields import ops
import numpy as np
import matplotlib.pyplot as plt

xvec = np.linspace(-4, 4, 200)
pvec = np.linspace(-4, 4, 200)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

# --- Coherent state (Gaussian) ---
prog_coh = sf.Program(1)
eng_coh = sf.Engine("fock", backend_options={"cutoff_dim": 25})
with prog_coh.context as q:
    ops.Dgate(1.5, 0.0) | q[0]
result_coh = eng_coh.run(prog_coh)
W_coh = result_coh.state.wigner(0, xvec, pvec)
axes[0].contourf(xvec, pvec, W_coh, levels=60, cmap="RdBu_r")
axes[0].set_title("Coherent state (Gaussian)")
axes[0].set_xlabel("x")
axes[0].set_ylabel("p")
axes[0].set_aspect("equal")

# --- Coherent state + Kerr gate (non-Gaussian) ---
prog_kerr = sf.Program(1)
eng_kerr = sf.Engine("fock", backend_options={"cutoff_dim": 25})
with prog_kerr.context as q:
    ops.Dgate(1.5, 0.0) | q[0]
    ops.Kgate(0.5) | q[0]  # Kerr gate with kappa=0.5
result_kerr = eng_kerr.run(prog_kerr)
W_kerr = result_kerr.state.wigner(0, xvec, pvec)
axes[1].contourf(xvec, pvec, W_kerr, levels=60, cmap="RdBu_r")
axes[1].set_title("After Kerr gate (non-Gaussian)")
axes[1].set_xlabel("x")
axes[1].set_ylabel("p")
axes[1].set_aspect("equal")

plt.tight_layout()
plt.savefig("kerr_wigner.png", dpi=150)
plt.show()

The coherent state Wigner function is a perfect round Gaussian. After the Kerr gate, it develops a crescent or banana shape, which is a hallmark of non-Gaussian character. The Wigner function also dips below zero in some regions, confirming genuine quantum non-classicality that no Gaussian state can produce.

The theoretical result: a Kerr gate together with the full set of Gaussian gates forms a universal gate set for CV quantum computing. In practice, strong Kerr nonlinearities are difficult to achieve in photonic hardware, which is why most current photonic experiments focus on Gaussian circuits with non-Gaussian measurements (photon counting) instead.

CV Quantum Key Distribution: A Simplified GG02 Protocol

Continuous-variable quantum key distribution (CV-QKD) is one of the most practical applications of CV quantum information. Unlike discrete-variable QKD (e.g., BB84), CV-QKD uses standard telecom components: lasers, homodyne detectors, and optical fibers. No single-photon sources or detectors are needed.

The GG02 protocol works as follows:

1. Alice draws random values x_A from a Gaussian distribution and prepares coherent states |alpha⟩ with alpha = x_A (real displacement in the X quadrature).
2. Bob measures the X quadrature using homodyne detection, obtaining outcomes x_B.
3. Due to vacuum noise and channel loss, x_B = sqrt(T) * x_A + noise, where T is the channel transmittance.
4. Alice and Bob publicly compare a subset of their data to estimate the noise level. An eavesdropper (Eve) would introduce excess noise beyond the vacuum level, which Alice and Bob can detect.
5. From the correlated Gaussian variables (x_A, x_B), Alice and Bob extract a shared secret key using classical post-processing.

Here is a simplified simulation:

import strawberryfields as sf
from strawberryfields import ops
import numpy as np

n_signals = 500
modulation_variance = 2.0  # Alice's modulation strength

# Alice prepares random displacements
np.random.seed(42)
alice_x = np.random.normal(0, np.sqrt(modulation_variance), n_signals)

# Bob measures X quadrature on each coherent state
bob_x = []
for x_a in alice_x:
    prog = sf.Program(1)
    eng = sf.Engine("gaussian")
    with prog.context as q:
        # Alice prepares a coherent state displaced by x_a in the X quadrature
        ops.Dgate(np.abs(x_a), 0.0 if x_a >= 0 else np.pi) | q[0]
        # Bob measures X
        ops.MeasureHomodyne(0) | q[0]
    result = eng.run(prog)
    bob_x.append(result.samples[0, 0])

bob_x = np.array(bob_x)

# Compute correlation
correlation = np.corrcoef(alice_x, bob_x)[0, 1]
print(f"Correlation coefficient: {correlation:.4f}")

# Estimate mutual information I(A;B) for Gaussian variables
# I(A;B) = -0.5 * log2(1 - rho^2) where rho is the correlation coefficient
mutual_info = -0.5 * np.log2(1 - correlation**2)
print(f"Mutual information I(A;B): {mutual_info:.4f} bits per symbol")

# In a real protocol, Alice and Bob would compare a fraction of their data
# to estimate the channel noise, then use the remaining data for key generation.
# The secret key rate is: K = I(A;B) - chi(B;E)
# where chi(B;E) is Eve's Holevo information, bounded by the observed noise.
print(f"\nAlice variance: {np.var(alice_x):.4f}")
print(f"Bob variance: {np.var(bob_x):.4f}")
print(f"Conditional variance Var(B|A): {np.var(bob_x - correlation * np.std(bob_x) / np.std(alice_x) * alice_x):.4f}")

In a noiseless channel (no Eve, no loss), the correlation should be very high and the mutual information should be close to its maximum. Any excess noise degrades the correlation and reduces the extractable key rate. If the noise exceeds a threshold, Alice and Bob abort the protocol because Eve may have too much information.

Cloud Access via Borealis

import strawberryfields as sf

# Authenticate with your Xanadu Cloud account
sf.store_account("YOUR_API_KEY")

# Connect to Borealis
eng = sf.RemoteEngine("borealis")

# Check device specifications
spec = eng.device_spec
print(f"Modes: {spec.modes}")
print(f"Max squeezing: {spec.gate_parameters['Sgate']['r']['max']}")

Borealis Hardware Architecture

The Borealis processor uses time-domain multiplexing (TDM), which is a fundamentally different approach from spatial-mode photonics. Instead of having one physical waveguide per mode, all modes travel through a single fiber loop as sequential time bins of light.

The key components:

• Squeezed light source: a single source that emits squeezed pulses at a fixed repetition rate. Each pulse is one mode.
• Three nested delay loops: with delays of 1, 6, and 36 time bins. Each loop contains a programmable beamsplitter (implemented by an electro-optic modulator) that mixes adjacent time bins separated by the loop’s delay length.
• Threshold detectors: click/no-click detectors (not photon-number-resolving). Each detector reports whether zero or at least one photon arrived.

This architecture supports 216 modes with up to approximately 6 dB of squeezing. The three delay loops create a specific connectivity pattern: each mode can interfere with modes that are 1, 6, or 36 time bins away. This means not all arbitrary interferometers are implementable on Borealis. You must check the device specification for allowed gate parameters and connectivity:

# Inspect available gate parameters
print("Gate parameters:")
for gate, params in spec.gate_parameters.items():
    print(f"  {gate}: {params}")

# Check the allowed connectivity and loop structure
print(f"\nTarget: {spec.target}")
print(f"Layout: {spec.layout}")

When designing circuits for Borealis, keep in mind:

• The circuit must follow the loop-based TDM structure. You cannot apply arbitrary beamsplitters between any pair of modes.
• Measurement is threshold detection only (not photon-number-resolving), so the output is a binary string, not photon counts.
• Compile your program using spec.compile() to verify it is compatible with the hardware before submitting.

Common Mistakes and Pitfalls

1. Fock Backend Cutoff Too Small

The cutoff dimension must be large enough to capture the significant photon number components of your state. For a coherent state with amplitude |alpha|, the mean photon number is |alpha|^2 and the standard deviation is |alpha|. A safe cutoff is |alpha|^2 + 4*|alpha| + 5.

import strawberryfields as sf
from strawberryfields import ops
import numpy as np

# BAD: alpha=3.0 has mean photon number 9, but cutoff=10 misses the tail
prog_bad = sf.Program(1)
eng_bad = sf.Engine("fock", backend_options={"cutoff_dim": 10})
with prog_bad.context as q:
    ops.Dgate(3.0, 0.0) | q[0]
result_bad = eng_bad.run(prog_bad)
probs_bad = result_bad.state.all_fock_probs()
print(f"Bad cutoff: sum of probs = {sum(probs_bad):.4f}")  # significantly less than 1.0

# GOOD: cutoff=25 captures the full distribution
prog_good = sf.Program(1)
eng_good = sf.Engine("fock", backend_options={"cutoff_dim": 25})
with prog_good.context as q:
    ops.Dgate(3.0, 0.0) | q[0]
result_good = eng_good.run(prog_good)
probs_good = result_good.state.all_fock_probs()
print(f"Good cutoff: sum of probs = {sum(probs_good):.4f}")  # very close to 1.0

If the probabilities do not sum to approximately 1.0, your cutoff is too small and the results are unreliable.

2. Gaussian Backend Limitations

The Gaussian backend cannot simulate non-Gaussian operations. Specifically:

• MeasureFock() is not available (photon counting is a non-Gaussian measurement)
• Kgate() and Vgate() (cubic phase) are not supported
• Fock() state preparation is not supported

If you try to use these operations with the Gaussian backend, Strawberry Fields raises an error. Use the Fock backend for any circuit involving non-Gaussian elements.

3. Phase Conventions

The BSgate(theta, phi) convention in Strawberry Fields uses theta as the transmissivity angle: theta=0 means full transmission (no mixing), theta=pi/4 means 50/50 splitting. The phi parameter is an additional phase.

This differs from other frameworks. For example, Perceval (Quandela’s framework) parameterizes the beamsplitter differently. Always check the convention when porting circuits between frameworks.

4. Homodyne Measurement Angle

The angle parameter in MeasureHomodyne(phi) is in radians, not degrees or fractions of pi. Common values:

• MeasureHomodyne(0): measures X quadrature
• MeasureHomodyne(np.pi / 2): measures P quadrature
• MeasureHomodyne(np.pi / 4): measures (X + P) / sqrt(2)

5. Borealis Circuit Constraints

Not all circuits can run on Borealis. The time-domain multiplexing architecture imposes specific connectivity constraints. Before submitting a job:

# Always compile against the device spec first
compiled_prog = prog.compile(device=spec)

If compilation fails, your circuit requires connectivity that the hardware does not support. Redesign the circuit to respect the loop structure (modes can only interact with modes separated by 1, 6, or 36 time bins).

What to Try Next

• Compute a graph adjacency matrix Hafnian using GBS with the thewalrus library to find dense subgraphs
• Experiment with the Kerr gate at different kappa values to see how non-Gaussian features evolve
• Simulate a CV teleportation protocol using two-mode squeezing, beamsplitters, and homodyne measurement
• Implement a variational circuit that optimizes gate parameters to prepare a target state
• Explore heterodyne detection and its role in CV quantum state tomography
• Read about continuous-variable quantum error correction using the GKP (Gottesman-Kitaev-Preskill) encoding
• See the Strawberry Fields Reference for the full gate and backend API