Strawberry Fiels

Strawberry Fields:
A Software Platform for Photonic Quantum Computing

Nathan Killoran, Josh Izaac, Nicolás Quesada, Ville Bergholm, Matthew Amy, and
Christian Weedbrook
Xanadu, 372 Richmond St W, Toronto, M5V 1X6, Canada
We introduce Strawberry Fields, an open-source quantum programming architecture for light-based quantum computers,
and detail its key features. Built in Python, Strawberry Fields is a full-stack library for design, simulation, optimization, and
arXiv:1804.03159v2 [quant-ph] 4 Mar 2019
quantum machine learning of continuous-variable circuits. The platform consists of three main components: (i) an API for
quantum programming based on an easy-to-use language named Blackbird; (ii) a suite of three virtual quantum computer
backends, built in NumPy and TensorFlow, each targeting specialized uses; and (iii) an engine which can compile Black-
bird programs on various backends, including the three built-in simulators, and – in the near future – photonic quantum
information processors. The library also contains examples of several paradigmatic algorithms, including teleportation,
(Gaussian) boson sampling, instantaneous quantum polynomial, Hamiltonian simulation, and variational quantum circuit
optimization.
ing theme for these software efforts is to target qubit-

based quantum devices. In reality, there are several com-
Introduction peting models of quantum computing which are equiva-
lent in computational terms, but which are conceptually
The decades-long worldwide quest to build practical quite distinct. One prominent approach is the continuous
quantum computers is currently undergoing a critical pe- variable (CV) model of quantum computing [48–50]. In
riod. During the next few years, a number of different the CV model, the basic information-processing unit is an
quantum devices will become available to the public. While infinite-dimensional bosonic mode, making it particularly
fault-tolerant quantum computers will one day provide sig- well-suited for implementations and applications based on
nificant computational speedups for problems like factor- light. The CV model retains the computational power of
ing [1], search [2], or linear algebra [3], near-term quan- the qubit model (cf. Chap. 4 of Ref. [51]), while offering a
tum devices will be noisy, approximate, and sub-universal number of unique features. For instance, the CV model is a
[4]. Nevertheless, these emerging quantum processors are natural fit for simulating bosonic systems (electromagnetic
expected to be strong enough to show a computational ad- fields, trapped atoms, harmonic oscillators, Bose-Einstein
vantage over classical computers for certain problems, an condensates, phonons, or optomechanical resonators) and
achievement known as quantum computational supremacy. for settings where continuous quantum operators – such as
As we approach this milestone, work is already under- position & momentum – are present. Even in classical com-
way exploring how such quantum processors might best puting, recent advances from deep learning have demon-
be leveraged. Popular techniques include variational quan- strated the power and flexibility of a continuous represen-
tum eigensolvers [5, 6], quantum approximate optimiza- tation of computation [52, 53] in comparison to the discrete
tion algorithms [7, 8], sampling from computationally hard computational model which has historically dominated.
probability distributions [9–17], and quantum annealing Here we introduce Strawberry Fields1 , an open-source
[18, 19]. Notably, these methodologies can be applied to software architecture for photonic quantum computing.
tackle important practical problems in chemistry [20–23], Strawberry Fields is a full-stack quantum software plat-
finance [24], optimization [25, 26], and machine learning form, implemented in Python, specifically targeted to the
[27–34]. These known applications are very promising, yet CV model. Its main element is a new quantum program-
it is perhaps the unknown future applications of quantum ming language named Blackbird. To lay the groundwork
computers that are most tantalizing. We may not know the for future photonic quantum computing hardware, Straw-
best applications of quantum computers until these devices
become available more widely to researchers, students, en-
trepreneurs, and programmers worldwide.
1
To this end, a nascent quantum software ecosystem has This document refers to Strawberry Fields version 0.9. Full documenta-
recently begun to develop [35–47]. However, a prevail- tion is available online at strawberryfields.readthedocs.io.
2
berry Fields also includes a suite of three CV quantum simu-

lator backends implemented using NumPy [54] and Tensor-
Flow [55]. Strawberry Fields comes with a built-in engine
to convert Blackbird programs to run on any of the simula-
Quantum Computation with Continuous
tor backends or, when they are available, on photonic quan- Variables
tum computers. To accompany the library, an online service
for interactive exploration and simulation of CV circuits is
available at strawberryfields.ai. Many physical systems in nature are intrinsically contin-
Aside from being the first quantum software frame- uous, with light being the prototypical example. Such sys-
work to support photonic quantum computation with tems reside in an infinite-dimensional Hilbert space, offer-
continuous-variables, Strawberry Fields provides additional ing a paradigm for quantum computation which is distinct
computational features not presently available in the quan- from the discrete qubit model. This continuous-variable
tum software ecosystem: model takes its name from the fact that the quantum op-
erators underlying the model have continuous spectra. It
1. We provide two numeric simulators; a Gaussian back- is possible to embed qubit-based computations into the CV
end, and a Fock-basis backend. These two formula- picture [57], so the CV model is as powerful as its qubit
tions are unique to the CV model of quantum compu- counterparts.
tation due to the use of an infinite Hilbert space, and
came with their own technical challenges.
From Qubits to Qumodes
(a) The Gaussian backend provides state-of-the-art
methods and functions for calculating the fi-
delity and Fock state probabilities, involving cal- A high-level comparison of CV quantum computation
culations of the classically intractable hafnian with the qubit model is depicted in Table I. In the remainder
[56]. of this section, we will provide a basic presentation of the
key elements of the CV model. A more detailed technical
(b) The Fock backend allows operations such as overview can be found in Appendix A. Readers experienced
squeezing and beamsplitters to be performed with CV quantum computing can safely skip to the next sec-
in the Fock-basis, a computationally intensive tion.
calculation that has been highly vectorized and
benchmarked for performance.
2. We provide a suite of important circuit decomposi- CV Qubit

tions appearing in quantum photonics – such as the
Williamson, Bloch-Messiah, and Clements decompo- Basic element Qumodes Qubits
sitions.
Relevant Quadratures x̂, p̂ Pauli operators
operators Mode operators â, â† σ̂ x , σ̂ y , σ̂z
3. The Fock-basis backend written using the TensorFlow
machine learning library allows for symbolic calcu- Coherent states |α〉
Pauli eigenstates
lations, automatic differentiation, and backpropaga- Common states Squeezed states |z〉
|0/1〉 , |±〉 , |±i〉
Number states |n〉
tion through CV quantum simulations. As far as we
are aware, this is the first quantum simulation library Rotation,
written using a high-level machine learning library, Displacement, Phase shift,
with support for dataflow programming and auto- Common gates Squeezing, Hadamard, CNOT,
matic differentiation. Beamsplitter, Cubic T-Gate
Phase
The remainder of this white paper is structured as fol- Homodyne |x φ 〉〈x φ |, Pauli eigenstates
Common
lows. Before presenting Strawberry Fields, we first provide Heterodyne π1 |α〉〈α|, |0/1〉〈0/1|, |±〉〈±|,
measurements
a brief overview of the key ingredients for CV quantum com- Photon-counting |n〉〈n| | ± i〉〈±i|
putation, specifically the most important states, gates, and
measurements. We then introduce the Strawberry Fields ar- Table I: Basic comparison of the CV and qubit settings.
chitecture in full, presenting the Blackbird quantum assem-
bly language, outlining how to use the library for numerical
simulation, optimization, and quantum machine learning.
Finally, we discuss the three built-in simulators and the in- The most elementary CV system is the bosonic harmonic
ternal representations that they employ. In the Appendices, oscillator, defined via the canonical mode operators â and
we give further mathematical and software details and pro- â† . These satisfy the well-known commutation relation
vide full example code for a number of important CV quan- [â, â† ] = I. It is also common to work with the quadra-
tum computing tasks. ture operators (also called the position & momentum oper-
3
ators)2 ,
v
tħh
x̂ := (â + â† ), (1)
2
v
tħh
p̂ := −i (â − â† ), (2)
2
where [x̂, p̂] = iħ

hÎ. We can picture a fixed harmonic oscil-
lator mode (say, within an optical fibre or a waveguide on a
photonic chip) as a single ‘wire’ in a quantum circuit. These
qumodes are the fundamental information-carrying units of
CV quantum computers. By combining multiple qumodes –
each with corresponding operators âi and âi† – and interact-
ing them via sequences of suitable quantum gates, we can
FIG. 1: Schematic representation of a Gaussian state for a
implement a general CV quantum computation.
single mode. The shape and orientation are parameterized
by the displacement α and squeezing z = r exp(iφ).
CV States
The dichotomy between qubit and CV systems is perhaps
most evident in the basis expansions of quantum states: the centre of the distribution, while the squeezing deter-
mines the variance and rotation of the distribution (see Fig.
Qubit |φ〉 = φ0 |0〉 + φ1 |1〉 , (3) 1). Multimode Gaussian states, on the other hand, are pa-
Z
rameterized by a vector of displacements r̄ and a covariance
Qumode |ψ〉 = d x ψ(x) |x〉 . (4) matrix V. Many important pure states in the CV model are
special cases of the pure Gaussian states; see Table II for a
For qubits, we use a discrete set of coefficients; for CV sys- summary.
tems, we can have a continuum. The states |x〉 are the eigen-
states of the x̂ quadrature, x̂ |x〉 = x |x〉, with x ∈ R. These
State family Displacement Squeezing
quadrature states are special cases of a more general family
of CV states, the Gaussian states, which we now introduce. Vacuum state |0〉 α=0 z=0
Gaussian states Coherent states |α〉 α∈C z=0
Our starting point is the vacuum state |0〉. Other states Squeezed states |z〉 α=0 z∈C
can be created by evolving the vacuum state according to Displaced squeezed
α∈C z∈C
states |α, z〉
|ψ〉 = exp(−i t H) |0〉 , (5)
α ∈ C,q
x̂ eigenstates |x〉 φ = 0, r → ∞
where H is a bosonic Hamiltonian (i.e., a function of the x = 2 ħ2h Re(α)
operators âi and âi† ) and t is the evolution time. States α ∈ C,q
where the Hamiltonian H is at most quadratic in the oper- p̂ eigenstates |p〉 φ = π, r → ∞
p = 2 ħ2h Im(α)
ators âi and âi† (equivalently, in x̂ i and p̂i ) are called Gaus-
sian. For a single qumode, Gaussian states are parameter- Fock states |n〉 N.A. N.A.
ized by two continuous complex variables: a displacement
parameter α ∈ C and a squeezing parameter z ∈ C (of- Table II: Common single-mode pure states and their relation
ten expressed as z = r exp(iφ), with r ≥ 0). Gaussian to the displacement and squeezing parameters. All listed
states are so-named because we can identify each Gaussian families are Gaussian, except for the Fock states. The n = 0
state with a corresponding Gaussian distribution. For single Fock state is also the vacuum state.
qumodes, the identification proceeds through its displace-
ment and squeezing parameters. The displacement gives
Fock states
Complementary to the continuous Gaussian states are the
2
discrete Fock states (or number states) |n〉, where n are non-
It is common to picture ħ h as a (dimensionless) scaling parameter for
negative integers. These are the eigenstates of the number
the x̂ and p̂ operators rather than a physical constant. However, there
are several conventions for the scaling value in common use [58]. operator n̂ = â† â. The Fock states form a discrete (count-
These self-adjoint operators are proportional to the Hermitian and anti- able) basis for qumode systems. Thus, each of the Gaussian
Hermitian parts of the operator â. Strawberry Fields allows the user to states considered in the previous section can be expanded
h = 2.
specify this value, with the default ħ in the Fock-state basis. For example, coherent states have
4
the form universality. A number of other useful CV gates are listed in

∞
Appendix B.

|α|2
X αn
|α〉 = exp − 2 p |n〉 , (6)
n=0 n! Gate Unitary Symbol
while (undisplaced) squeezed states only have even number Displacement Di (α) = exp (αâi† − α∗ âi ) D
states in their expansion:
Rotation R i (φ) = exp (iφ n̂i ) R
∞ p
1 X (2n)!
|z〉 = p [−e iφ tanh(r)]n |2n〉 . (7) Squeezing Si (z) = exp ( 21 (z ∗ âi2 − z âi†2 )) S
cosh r n=0 2n n!
BSi j (θ , φ) =
Mixed states Beamsplitter BS
exp (θ (e iφ âi â†j − e−iφ âi† â j ))
Mixed Gaussian states are also important in the CV pic- γ

Cubic phase Vi (γ) = exp i 3ħh x̂ i3 V
ture, for instance, the thermal state
∞ Table III: Some important CV model gates. All listed gates
X nn except the cubic phase gate are Gaussian.
ρ(n) := |n〉〈n|, (8)
n=0
(1 + n)n+1
which is parameterized via the mean photon number n :=

CV Measurements
Tr(ρ(n)n̂). Starting from this state, we can consider a
mixed-state-creation process similar to Eq. (5), namely As with CV states and gates, we can distinguish between
Gaussian and non-Gaussian measurements. The Gaus-
ρ = exp(−i t H)ρ(n) exp(i t H). (9) sian class consists of two (continuous) types: homodyne
and heterodyne measurements, while the key non-Gaussian
Analogously to pure states, by applying Hamiltonians of measurement is photon counting. These are summarized in
second-order (or lower) to thermal states, we generate the Table IV.
family of Gaussian mixed states.
Homodyne measurements
CV Gates
Ideal homodyne detection is a projective measurement
Unitary operations can be associated with a generating onto the eigenstates of the quadrature operator x̂. These
Hamiltonian H via the recipe (cf. Eqs. (5) & (9)) states form a continuum, so homodyne measurements are
inherently continuous, returning values x ∈ R. More gen-
U := exp (−i t H). (10) erally, we can
consider projective measurement onto the
eigenstates x φ of the Hermitian operator
For convenience, we classify unitaries by the degree of their
generator. A CV quantum computer is said to be universal x̂ φ := cos φ x̂ + sin φ p̂. (11)
if it can implement, to arbitrary precision and with a finite
number of steps, any unitary which is polynomial in the This is equivalent to rotating the state clockwise by φ and
mode operators [48]. We can build a multimode unitary by performing an x̂-homodyne measurement. If we have a
applying a sequence of gates from a universal gate set, each multimode Gaussian state and we perform homodyne mea-
of which acts only on one or two modes. We focus on a surement on one of the modes, the conditional state of the
universal set made from the following two subsets: unmeasured modes remains Gaussian.
Gaussian gates: Single and two-mode gates which are at
Heterodyne measurements
most quadratic in the mode operators, e.g., Displace-
ment, Rotation, Squeezing, and Beamsplitter gates. Whereas homodyne detection is a measurement of x̂,
heterodyne detection can be seen as a simultaneous mea-
Non-Gaussian gate: A single-mode gate which is degree 3 surement of both x̂ and p̂. Because these operators do not
or higher, e.g., the Cubic phase gate. commute, they cannot be simultaneously measured without
A number of fundamental CV gates are presented in Ta- some degree of uncertainty. Equivalently, we can picture
ble III. Many of the Gaussian states from the previous sec- heterodyne measurement as projection onto the coherent
tion are connected to a corresponding Gaussian gate. Any states, with measurement operators π1 |α〉〈α|. Because the
multimode Gaussian gate can be implemented through a coherent states are not orthogonal, there is a corresponding
suitable combination of Displacement, Rotation, Squeezing, lack of sharpness in the measurements. If we perform het-
and Beamsplitter Gates [50], making these gates sufficient erodyne measurement on one mode of a multimode state,
for constructing all quadratic unitaries. The cubic phase the conditional state on the remaining modes stays Gaus-
gate is presented as an exemplary non-Gaussian gate, but sian.
any other non-Gaussian gate could also be used to achieve
5
Measurement Measurement These elements provide access points for users to design
Measurement quantum circuits. These circuits are then linked to a back-
Operators values
end via a quantum compiler engine. For a backend, the
Homodyne |x φ 〉〈x φ | x ∈R engine can currently target one of three included quantum
computer simulators. When CV quantum processors be-
1 α∈C
Heterodyne π |α〉〈α| come available in the near future, the engine will also build
and run circuits on those devices. Further, high-level quan-
Photon counting |n〉〈n| n∈N
tum computing applications can be built by leverging the
Strawberry Fields frontend API. Existing examples include
Table IV: Key measurement types for the CV model. The
the Strawberry Fields Interactive website, the Quantum Ma-
‘-dyne’ measurements are Gaussian, while photon-counting is
chine Learning Toolbox (for streamlining the training of
non-Gaussian.
variational quantum circuits), and SFOpenBoson (an inter-
face between the electronic structure library OpenFermion
[45] and Strawberry Fields).
Photon Counting In the remainder of this section, the key elements of
Photon counting (also known as as photon-number re- Strawberry Fields will be presented in more detail. Proceed-
solving measurement), is a complementary measurement ing through a series of examples, we show how CV quantum
method to the ‘-dyne’ measurements, revealing the particle- computations can be defined using the Blackbird language,
like, rather than the wave-like, nature of qumodes. This then compiled and run on a quantum computer backend.
measurement projects onto the number eigenstates |n〉, re- We also outline how to use Strawberry Fields for optimiza-
turning non-negative integer values n ∈ N. Except for the tion and machine learning on quantum circuits. Finally, we
outcome n = 0, a photon-counting measurement on a sin- discuss the suite of quantum computer simulators included
gle mode of a multimode Gaussian state will cause the re- within Strawberry Fields.
maining modes to become non-Gaussian. Thus, photon-
counting can be used as an ingredient for implementing Blackbird: A Quantum Programming Language
non-Gaussian gates. A related process is photodetection, As classical computers have become progressively more
where a detector only resolves the vacuum state from non- powerful, the languages used to program them have also
vacuum states. This process has only two measurement op- undergone considerable paradigmatic changes. Machine
erators, namely |0〉〈0| and I − |0〉〈0|. code gave way to human-readable assembly languages, fol-
lowed by higher-level procedural and object-oriented lan-
guages. With each generation, the trend has been towards
higher levels of abstraction, separating the programmer
The Strawberry Fields Software Platform more and more from details of the actual computer hard-
ware. Quantum computers are still at an early stage of
The Strawberry Fields library has been designed with sev- development, so while we can imagine what higher-level
eral key goals in mind. Foremost, it is a standard-bearer for quantum programming might look like, in the near term we
the CV model, laying the groundwork for future photonic first need to build languages which are conceptually closer
quantum computers. As well, Strawberry Fields is designed to the quantum hardware.
to be simple to use, giving entry points for as many users Blackbird is a standalone domain specific language (DSL)
as possible. Finally, since the potential applications of near- for continuous-variable quantum computation. With a well-
term quantum computers are still being worked out, it is defined grammar in extended Backus-Naur form, and both
important that Strawberry Fields provides powerful tools Python and C++ parsers available, Blackbird provides op-
to easily explore many different use-cases and applications. erations that match the basic CV states, gates, and mea-
Strawberry Fields has been implemented in Python, a surements, and maps directly to low-level hardware instruc-
modern language with a gentle learning curve which is al- tions. The abstract syntax keeps a close connection between
ready familiar to many programmers and scientific practi- code and the quantum operations that they implement; this
tioners. The accompanying quantum simulator backends syntax is modeled after that of ProjectQ [41], but special-
are built upon the widely used Python packages NumPy ized to the CV setting. Blackbird can be used as part of
and TensorFlow. All Strawberry Fields code is open source. the Strawberry Fields stack, but also directly with photonic
Strawberry Fields can be accessed programmatically as a quantum computing hardware systems.
Python package, or via a browser-based interface for de- Within the Strawberry Fields framework, we have built
signing quantum circuits. an implementation of Blackbird using Python 3 as the em-
A pictorial outline of Strawberry Fields’ key elements and bedding language — an ‘embedded’ DSL. This ‘Python-
their interdependencies is presented in Fig. 2. Conceptu- enhanced’ Blackbird language provides the same core oper-
ally, the software stack is separated into two main pieces: a ations and follows the same grammar and syntactical rules
user-facing frontend layer and a lower-level backends com- as the standalone DSL, but, by nature, may also contain
ponent. The frontend encompasses the Strawberry Fields valid Python constructs. Furthermore, Strawberry Fields’
Python API and the Blackbird quantum assembly language. ‘Python-enhanced’ Blackbird provides users with additional
6
FRONT-ENDS BACK-ENDS
Strawberry Fields
Quantum Processor
Interactive Server
Simulators
Blackbird
Strawberry Quantum COMPILER Fock Representation
Fields API Programming ENGINE NumPy Tensorflow
Language
Gaussian Representation
NumPy
Applications
FIG. 2: Outline of the Strawberry Fields software stack. The Strawberry Fields Interactive server is available online at
strawberryfields.ai.
quantum operations that are decomposed into lower-level Conceptually, the vertical bar symbol ‘|’ separates Opera-
Blackbird assembly commands. We will introduce the ele- tions – like state preparation – from the registers that they
ments of Blackbird through a series of basic examples, dis- act upon. Notice that we can use Operations inline, or con-
cussing more technical aspects as they arise. struct them separately and reuse them several times.
After creating states, we will want to transform these us-
Operations ing quantum gates.
Quantum computations consist of four main ingredients:
1 # Apply the Displacement gate to qumode 0
state preparations, application of gates, performing mea-
2 alpha = 2.0 + 1j
surements, and adding/removing subsystems. In Blackbird, 3 Dgate(alpha) | q[0]
these are all considered as Operations, and share the same 4
basic syntax. In the following code examples, we use the 5 # Apply the Rotation gate
variable q for a set of qumodes (more specifically, a quan- 6 phi = 1.157
tum register), the details of which are deferred until the next 7 Rgate(phi) | q[0]
section. 8
Our first considered Operation is state preparation. By 9 # Apply the Squeezing gate
default, qumodes are initialized in the vacuum state. Var- 10 Sgate(2.0, 0.17) | q[0]
11
ious other important CV states can be created with simple
12 # Apply the Beamsplitter gate to qumodes 0 & 1
Blackbird commands. 13 BSgate(0.314, 0.223) | (q[0], q[1])
14
1 # Create the vacuum state in qumode 0 15 # Apply Cubic phase gate (VGate) to qumode 0
2 Vac | q[0] 16 gamma = 0.1
3 17 Vgate(gamma) | q[0]
4 # Create a coherent state in qumode 1 18
5 alpha = 2.0 + 1j 19 # Apply Hermitian conjugate of a gate

6 Coherent(alpha) | q[1] 20 V = Vgate(gamma)
7 21 V.H | q[0]
8 # Create squeezed states in qumodes 0 & 1
9 S = Squeezed(2.0) Codeblock 2: Blackbird code for applying various CV gates.
10 S | q[0]
11 S | q[1] Blackbird supports all of the gates listed in the previous sec-
12 tion as well as a number of composite gates, each of which
13 # Create a Fock state in qumode 1 can be decomposed using the universal gates. The sup-
14 Fock(4) | q[1] ported composite gates are: controlled X (CXgate), con-
trolled Z (CZgate), quadratic phase (Pgate), and two-
Codeblock 1: Blackbird code for creating various CV mode squeezing (S2gate). A full list of gates currently
quantum states. supported in Blackbird can be found in Appendix B.
Finally, we can specify measurement Operations using
Blackbird state preparations such as those used in Code-
Blackbird.
block 1 implicitly reset the existing state of the qumodes.
7
1 # Homodyne measurement at angle phi 24 eng.run(backend="gaussian")

2 phi = 0.785 25 vals = [reg.val for reg in q]
3 MeasureHomodyne(phi) | q[0]
4 Codeblock 4: Code for declaring and running Blackbird
5 # Special homodyne measurements programs using the Strawberry Fields library.
6 MeasureX | q[0]
7 MeasureP | q[1] The above code example is runnable and carries out a com-
8 plete quantum computation, namely the preparation and
9 # Heterodyne measurement measurement of an EPR entangled state. Note that Oper-
10 MeasureHeterodyne() | q[0] ations can be declared outside of the Engine, but their ac-
11 MeasureHD | q[1] # shorthand tion on the quantum registers must come within the En-
12
gine context. Also notice that our register has a length
13 # Number state measurements of various qumodes
14 MeasureFock() | q[0] of 2, so any single-mode Operations must act on specific
15 MeasureFock() | (q[1], q[2]) # multiple modes elements, i.e., q[i], while two-mode Operations can act
16 Measure | q[3] # shorthand on q directly. Finally, the user must specify a backend –
as well as any backend-dependent settings – when calling
Codeblock 3: Blackbird code for carrying out CV eng.run(). We will discuss the Strawberry Fields back-
measurements. ends further in a later section.
Measurements have several effects. For one, the numeri- Quantum and Classical Registers
cal result of the measurement is placed in a classical regis-
ter. As well, the state of all remaining qumodes is projected When a Strawberry Fields Engine is constructed, the
to the (normalized) conditional state for that measurement user must specify the number of qumode subsystems to be-
value. Finally, the state of the measured qumode is reset gin with. This number is required for the initialization of
to the vacuum state. This is the typical behaviour of pho- the Engine, but may change within a computation (e.g.,
tonic hardware, where measurements absorb all the energy when temporary ancilla modes are used). Qumodes can
of the measured qumode. be added/deleted by using the New and Del Operations.
Running Blackbird Programs in Python 1 # A Blackbird circuit where gates

2 # are added and deleted
Within Python, Blackbird programs are managed by an 3 alice = q[0]
Engine. The function Engine in the Strawberry Fields API 4 with eng:
will instantiate an Engine, returning both the Engine and 5 Sgate(1) | alice
its corresponding quantum register. The Engine is used as 6 bob, charlie = New(2)
a Python context manager, providing a convenient way to 7 BSgate(0.5) | (alice, bob)
encapsulate Blackbird programs. 8 CXgate(1) | (alice, charlie)
9 Del | alice
1 # Create Engine and quantum register 10 S2gate(0.4) | (charlie, bob)
2 import strawberryfields as sf 11
3 eng, q = sf.Engine(num_subsystems=2) 12 # Attempting to act on registers which have

4
13 # been removed will raise an IndexError
5 # The register is also available via 14 try:
,→ eng.register 15 with eng:
6 assert q == eng.register 16 Dgate(0.1) | alice
7
17 except Exception as e:
8 # Put Blackbird Operations in namespace 18 assert isinstance(e, IndexError)
9 from strawberryfields.ops import *
10 Codeblock 5: Adding and deleting qumode subsystems.
11 # Declare a Blackbird program
12 from math import pi An Engine maintains a unique numeric indexing for the
13 z = 4. quantum registers based on the order they were added.
14 S = Sgate(z) When a subsystem is deleted from the circuit, no further
15 B = BSgate(pi / 4, 0) gates can act on that register.
16 with eng: As stated earlier, measurement Operations produce clas-
17 S | q[0] sical information (the measurement result) and manipu-
18 S.H | q[1] late the corresponding register. Compared to previously
19 B | q released quantum programming frameworks, such as Pro-
20 MeasureP | q[0] jectQ, PyQuil, and Qiskit, Strawberry Fields has been de-
21 MeasureP | q[1]
signed so that the notation q[i], while principally denoting
22
23 # Execute Blackbird program and extract values the quantum register, may also encapsulate classical infor-
mation or a classical register. This behaviour is contextual,
8
and unique to Strawberry Fields. For example, prior to mea- explore the conditional state created by a specific value, de-
surement, q[i] simply references a quantum register. Once termine the measurement-dependent corrections we need
a measurement operation is performed, q[i] continues to to make in a teleportation circuit, or even design an algo-
represent a quantum register — now reset to the vacuum rithm which inherently contains post-selection. This func-
state — as well as storing the numerical value of the mea- tionality is supported in Strawberry Fields through the op-
surement, accessible via the attribute q[i].val. Note that tional keyword argument select, which can be supplied
this numerical value is only available if a computation has for any measurement Operation. The measurement out-
been run up to the point of measurement. We may also use come will then return exactly this value, while the remain-
a classical measurement result symbolically as a parameter ing modes will be projected into the conditional state cor-
in later gates without first running the computation. To do responding to this value3 .
this, we simply pass the measured register (e.g., q[i]) ex-
plicitly as an argument to the required gate. As before, the 1 with eng:
Strawberry Fields quantum register object is contextual — 2 Fock(3) | q[0]
when passed as a gate argument, Strawberry Fields implic- 3 Fock(2) | q[1]
itly accesses the encapsulated classical register. 4 BSgate() | (q[0], q[1])
5 MeasureFock(select=4) | q[0]
6 MeasureFock() | q[1]
1 # Numerical evaluation of a measurement 7 eng.run("fock", cutoff_dim=6)
2 # result using eng.run() 8 assert q[0].val == 4
3 with eng:
4 MeasureX | q[0] Codeblock 8: Selecting a specific desired measurement
5 eng.run("gaussian")
outcome.
6 val = q[0].val
7
8 # Use a measured register symbolically Decompositions

9 # in another gate
10 with eng: In addition to the core CV operations discussed above,
11 MeasureX | q[0] Strawberry Fields also provides support for some impor-
12 Dgate(q[0]) | q[1] tant decompositions frequently used in quantum optics.
13 eng.run("gaussian") These include the (a) Williamson decomposition [59],
for decomposing arbitrary Gaussian states to a symplectic
Codeblock 6: Evaluating measurement results numerically transformation acting on a thermals state, (b) the Bloch-
and using them symbolically. Messiah decomposition [60–62], for decomposing the ac-
In quantum algorithms, it is common to process a measure- tion of symplectic transformations to interferometers and
ment result classically and use the post-processed value as single-mode squeezing, and (c) the Clements decomposi-
a parameter for further operations in a circuit. Strawberry tion [63], for decomposing multi-mode linear interferom-
Fields provides the convert decorator to transform a user- eters into arrays of beamsplitters and rotations of fixed
specified numerical function into one which acts on regis- depth. In all cases, the resulting decomposition into the
ters. universal CV gate set may be viewed via the engine method
eng.print_applied(). Strawberry Fields thus provides
a natural environment for embedding graphs and matrices
1 @sf.convert
2 def neg(x): in quantum optical circuits, and viewing the resulting phys-
3 return -x ical components.
4
5 # A Blackbird computation using classical 1 U = np.array([[1-1j, np.sqrt(2)],

6 # data processing 2 [-np.sqrt(2), 1+1j]])/2
7 with eng: 3
8 MeasureX | q[0] 4 eng, q = sf.Engine(2)

9 Dgate(neg(q[0])) | q[1] 5
10 eng.run("gaussian") 6 with eng:

7 Squeezed(0.43) | q[0]
Codeblock 7: Symbolically processing a measured value 8 Interferometer(U) | (q[0], q[1])
before using it. 9
10 eng.run("gaussian")
11 eng.print_applied()
Post-selection
The measurement Operations in Strawberry Fields are
stochastic in nature, with outcomes determined by some
underlying quantum probability distribution. Often it is 3
Users should be careful to avoid post-selection on measurement values
convenient to select specific values for these measurements which have no probability of occuring given the current circuit state. In
rather than sampling them. For instance, we might want to this case, the expected behaviour of a backend is not defined.
9
When we do this, any registers measured in the circuit will

12 # >> Squeezed(0.43, 0) | (q[0])
13 # >> Rgate(2.356) | (q[0]) be populated with unevaluated Tensor objects rather than
14 # >> BSgate(0.7854, 0) | (q[0], q[1]) numerical values (without eval=False, the TensorFlow
15 # >> Rgate(-3.142) | (q[0]) backend returns purely numerical results, similar to the
16 # >> Rgate(0.7854) | (q[1]) other simulators). These Tensors can still be evaluated nu-
merically by running them in a TensorFlow Session. In
Codeblock 9: Using the in-built Clements decomposition to this case, measurement results will be determined stochas-
decompose a 2 × 2 Interferometer into beamsplitters tically on each evaluation.
and phase rotations.
18 # Evaluate Tensors numerically
19 sess = tf.Session()
Optimization and Quantum Machine Learning 20 sess.run(tf.global_variables_initializer())
21 print(sess.run(q[0].val))
Strawberry Fields can perform quantum circuit simula-
22 # >> a numerical measurement value
tions using both numerical and symbolic representations. 23 print(sess.run(q[0].val))
Numerical computation is the default operating mode and 24 # >> a (different) numerical measurement value
is supported by all three supplied backends. Symbolic com-
putation is enabled only for the TensorFlow backend. In Codeblock 11: Numerically evaluating Tensors.
this section, we outline the main TensorFlow functionali-
When specifying a circuit in Blackbird, we can make
ties accessible through the Strawberry Fields frontend in-
use of various special symbolic TensorFlow classes, such as
terface. More details about the corresponding TensorFlow
backend are discussed in the next section. TensorFlow [55]
Variable, Tensor, placeholder, or constant.
models calculations abstractly using a computational graph,
26 # Declare circuits using Tensorflow objects
where individual operations are represented as nodes and 27 alpha = tf.Variable(0.1)
their dependencies by directed edges. This viewpoint sep- 28 phi = tf.constant(0.5)
arates the symbolic representation of a computation from 29 theta_bs = tf.placeholder(tf.float64)
its numerical evaluation, and makes optimization and ma- 30 phi_bs = tf.nn.relu(phi) # a Tensor object
chine learning more amenable. On top of this, TensorFlow 31
provides a number of advanced functionalities, including 32 eng, q = sf.Engine(num_subsystems=2)

automatic gradient computation, GPU utilization, built-in 33 with eng:
optimization algorithms, and various other machine learn- 34 Coherent(alpha) | q[0]
ing tools. Note that the term ‘quantum machine learning’ 35 Measure | q[0]
36 BSgate(theta_bs, phi_bs) | (q[0], q[1])
will be used here in a hybrid sense, i.e., applying conven-
37 MeasureHomodyne(phi) | q[1]
tional machine learning methods to quantum systems. 38 eng.run("tf",
To build a TensorFlow computational graph using Straw- 39 cutoff_dim=5,
berry Fields, we instantiate an Engine, declare a circuit in 40 eval=False)
Blackbird code, then execute eng.run on the TensorFlow 41
("tf") backend. To keep the underlying simulation fully 42 sess = tf.Session()

symbolic, the extra argument eval=False must be given. 43 sess.run(tf.global_variables_initializer())
44 feed_dict = {theta_bs: 0.5, q[0].val: 2}
1 # Create Engine and run symbolic computation 45 print(sess.run(q[1].val, feed_dict=feed_dict))
2 import strawberryfields as sf 46 # >> a numerical measurement value
Codeblock 12: Using abstract TensorFlow classes as circuit
4 import tensorflow as tf
5
parameters.
6 eng, q = sf.Engine(num_subsystems=1) In the above example, we supplied an additional
7 with eng: feed_dict argument when evaluating. This is a Python
8 Dgate(0.5) | q[0]
dictionary which specifies the numerical values (typically,
9 MeasureX | q[0]
10 eng.run("tf", coming from a dataset) for every placeholder that
11 cutoff_dim=5, appears in a circuit. However, as can be seen from the
12 eval=False) example, it is also possible to substitute desired values for
13 other nodes in the computation, including the values stored
14 # Registers contain symbolic measurements in quantum registers. This allows us to easily post-select
15 print(q[0].val) measurement values and explore the resulting conditional
16 # >> Tensor("Measure_homodyne/Meas_result:0", states.
,→ shape=(), dtype=float64) We can also perform post-processing of measurement re-
sults when working with TensorFlow objects. In this case,
Codeblock 10: Creating a TensorFlow computational graph the functions decorated by convert should be written us-
for a quantum circuit. ing TensorFlow Tensors, Variables, and operations.
10
By taking advantage of these additional functionalities of

48 # Use a simple neural network as
49 # a processing function the TensorFlow backend, we can straightforwardly perform
50 @sf.convert optimization and machine learning on quantum circuits in
51 def NN(x): Strawberry Fields [64, 65]. A complete code example for
52 weight = tf.Variable(0.5, dtype=tf.float64) optimization of a quantum circuit is located in Appendix C.
53 bias = tf.Variable(0.1, dtype=tf.float64)
54 return tf.sigmoid(weight * x + bias)
55
56 eng, q = sf.Engine(num_subsystems=2) Strawberry Fields’ Quantum Simulators

57 with eng:
58 MeasureP | q[0] The ultimate goal is for Blackbird programs to be carried
59 Dgate(NN(q[0])) | q[1]
out on photonic quantum computers. To lay the ground-
60 MeasureX | q[1]
61 eng.run("tf", work for these emerging devices, Strawberry Fields comes
62 cutoff_dim=5, with a suite of three CV quantum computer simulators spe-
63 eval=False) cially designed for the CV model. These simulators tar-
64 print([r.val for r in q]) get different use-cases and support different functionality.
65 # >> [<tf.Tensor For example, many important algorithms in the CV formal-
,→ 'Measure_homodyne_2/Meas_result:0' ism involve only Gaussian states, operations, and measure-
,→ shape=() dtype=float64>, <tf.Tensor ments. We can take advantage of this structure to more
,→ 'Measure_homodyne_3/Meas_result:0' efficiently simulate such computations. Other circuits have
,→ shape=() dtype=float64>] inherently non-Gaussian elements to them; for these, the
66
Fock basis provides the standard description. These repre-
Codeblock 13: Processing a measurement result using a sentations are available in Strawberry Fields in the Gaussian
neural network. For compactness, the example uses a width backend and the Fock backend, respectively. The third built-
1 perceptron, but any continuous processing function in backend is the TensorFlow backend. Also using the Fock
supported by TensorFlow can be used. representation, this backend is geared primarily towards
optimization and machine learning applications.
The TensorFlow backend additionally supports the use Most Blackbird operations are supported across all three
of batched processing, allowing for many evaluations of backends. A small subset, however, are not supported uni-
a quantum circuit to potentially be computed in paral- formly due to mathematical incompatibility. For example,
lel. Scalars are automatically broadcast to the specified the Gaussian backend does not generally support non-
batch size. Finally, we can easily run circuit simulations on Gaussian gates or the preparation/measurement of Fock
special-purpose hardware like GPUs or TPUs. states. Sometimes it is also useful to work with a backend
directly. To allow this, Strawberry Fields provides a back-
68 batch_size = 3 end API, giving access to additional methods and properties
69 eng, q = sf.Engine(num_subsystems=1) which are not part of the streamlined frontend API. A stan-
70
dalone backend can be created in Strawberry Fields using
71 alpha = tf.Variable([0.1] * batch_size)
72 # scalars are automatically cast to batch size
strawberryfields.backend.load_backend(name)
73 beta = tf.Variable(0.5) where name is one of "gaussian", "fock", or "tf" (for
74
comparison, the backend associated to an Engine eng is
75 sess = tf.Session() available via eng.backend).
76 sess.run(tf.global_variables_initializer()) Three important backend methods, common to all simu-
77 lators, are begin_circuit, reset, and state. The first
78 with tf.device("/gpu:0"): command instantiates the circuit simulation, the second re-
79 with eng: sets the simulation back to an initial vacuum state, clearing
80 Dgate(alpha) | q[0] all previous operations, and the third returns a class which
81 Dgate(beta) | q[0] encapsulates the current quantum state of the simulator.
82 MeasureX | q[0]
In addition to containing the numerical (or symbolic) state
83 eng.run("tf",
84 cutoff_dim=10, data, state classes also contain a number of useful meth-
85 batch_size=batch_size, ods and attributes for further exploring the quantum state,
86 eval=False) such as fidelity, mean_photon, or wigner. As a conve-
87 print(q[0].val) nience for the user, all simulations carried out via eng.run
88 # >> will return a state class representing the final circuit state
,→ Tensor("Measure_homodyne_4/Meas_result:0", (see Appendix C for examples).
,→ shape=(3,), dtype=float64,
,→ device=/device:GPU:0) Gaussian Backend
Codeblock 14: Running a batched computation and This backend, written in NumPy, uses the symplectic for-
explicitly placing the computation on a GPU. malism to represent CV systems. At a high level, this repre-
11
sentation tracks the quantum state of an N -mode quantum low that all states fit within the specified cutoff, then sim-
system using two Gaussian components: a 2N -dimensional ulations with the Fock and Gaussian backends will be in
displacement vector r̄ and a 2N × 2N -dimensional covari- numerical agreement.
ance matrix V (a deeper technical overview is located Like the Gaussian backend, the Fock backend has a
in Appendix A). After we have created a Gaussian state state method which encapsulates the numerical state,
(either via state = eng.run(backend="gaussian") while also providing a number of methods and attributes
or by directly calling the state method of a Gaussian specific to the Fock representation (such as ket, trace,
backend), we can access r̄ and V via state.means and and all_fock_probs.). Unlike the Gaussian representa-
state.cov, respectively. Other useful Gaussian state tion, mixed state simulations take up more resources than
methods are displacement and squeezing, which re- pure states. Pure states are represented in the Fock back-
turn the Gaussian parameters associated to the underlying end by an D N -dimensional complex vector and mixed states
state. by a D2N -dimensional density matrix. Because of this extra
The scaling of the symplectic representation with the overhead, by default the Fock backend will internally rep-
number of modes is O (N 2 ). On one hand, this is quite pow- resent a quantum circuit as long as possible as a pure state,
erful. It allows us to to efficiently simulate any computa- switching to the mixed state representation only when it
tions which are fully Gaussian. On the other, the formalism becomes necessary. Most importantly, for N > 2 qumodes,
is not expressive enough to simulate more general quan- all state preparation Operations (Vacuum, Squeezed, Fock,
tum computations. Only a small number of non-Gaussian etc.) cause the representation to become mixed. This is be-
methods are available for this backend. These are auxiliary cause the mode where the state is prepared could be entan-
methods where we extract some non-Gaussian information gled with other modes. To keep physically consistent, the
from a Gaussian state, but do not update the state of the Fock backend will first trace out the relevant mode, neces-
circuit. One such method is fock_prob, which is imple- sitating a mixed state representation. When possible, it is
mented using an optimized – yet still exponentially scaling – recommended to apply gates to the (default) vacuum state in
algorithm. This method enables simulation of the Gaussian order to efficiently prepare pure states. If desired, a mixed
boson sampling algorithm [10] using the Gaussian backend; state simulation can be enforced by passing the argument
see Appendix C for a complete code example. pure=False when calling begin_circuit.
Fock Backend TensorFlow Backend

This backend, also written in NumPy, uses a fundamen- The other built-in backend for Strawberry Fields is coded
tally different description for qumodes than the Gaussian using TensorFlow. As a simulator, it uses the same internal
representation. As discussed in the introductory sections, representation as the Fock backend (Fock basis, finite cut-
the Fock representation encodes quantum computation in off, pure vs. mixed state representations, etc.) and has the
a countably infinite-dimensional Hilbert space. This repre- same methods. It can operate as a numerical simulator sim-
sentation is faithful, allowing a precise description of CV ilar to the other backends. Its main purpose, however, is to
systems, in particular non-Gaussian circuits. Yet simulating leverage the many powerful tools provided by TensorFlow
infinite-dimensional systems leads to some computational to enable optimization and machine learning on quantum
tradeoffs which are not present for qubit simulators. Most circuits. Much of this functionality was presented in the
importantly, we impose a cutoff dimension D for simula- previous section, so we will not repeat it here.
tions (chosen by the user), so the Fock backend only covers Like the other simulators, users can query the TensorFlow
a restricted set of number states {|0〉 , ..., |D − 1〉} for each backend’s state method to access the internal representa-
mode. The size of simulated quantum systems thus depends tion of a circuit’s quantum state. This functions similarly
on both the number of subsystems N and the cutoff, being to the Fock backend’s state method, except that the state
O (D N ). We contrast this with qubit systems, where the base returned can be an unevaluated Tensor object when the
is fixed, i.e., O (2N ). While these scalings are both expo- keyword argument eval is set to False. This state Ten-
nential, in practice simulating qumode systems for D > 2 is sor can be combined with any supported TensorFlow op-
more computationally demanding than qubits. This is be- erations (norm, self_adjoint_eig, inv, etc.) to enable
cause the (truncated) qumode subsystems have a higher di- optimization and machine learning on various properties of
mension and thus encode more information than their qubit quantum circuits and quantum states.
counterparts.
Physically, imposing a cutoff is a reasonable strategy since
higher photon-number states must have higher energy and,
in practice, quantum-optical systems will have bounded en- Conclusions
ergy (e.g., limited by the power of a laser). On the other
hand, there are certainly states which can be easily pre- We have introduced Strawberry Fields, a multi-faceted
pared in the lab, yet would not fit accurately on the sim- software platform for continuous-variable quantum com-
ulator. Thus, some care must be taken to trade off between puting. The main components of this library – a custom
the numerical cutoff value, the number of modes, and the quantum programming language (Blackbird), a compiler
energy scale of the circuit. If the energy scale is sufficiently engine, and a suite of quantum simulators targeting distinct
12
applications – have been presented in detail. Further infor- The stage is now set for the broader community to use
mation is available in both the Appendices and the Straw- Strawberry Fields for exploration, research, and develop-
berry Fields online documentation. ment of new quantum algorithms, specialized circuits, and
machine learning models. We anticipate the creation of
further software applications and backend modules for the
Strawberry Fields platform (developed both internally and
externally), providing advanced functionality and applica-
tions for quantum computing and quantum machine learn-
ing.
Acknowledgements
We thank our colleagues at Xanadu for testing Strawberry
Fields, reviewing this white paper, and providing helpful
feedback. In particular, we thank Patrick Rebentrost for
valuable discussions and suggestions.
13
bic phase and rotation gates in Table III, while the two-
mode generator corresponds to the beamsplitter gate. Fi-
Appendix A: The CV model nally, the x̂ 2 generator corresponds to a quadratic phase
gate, P̂(s) = exp(is x̂ 2 /(2ħ h)). This gate can be written in
In this Appendix, we provide more technical and mathe- terms of the single-mode squeezing gate and the rotation
matical details about the CV model, the quantum comput- gate as follows: P̂(s) = R̂(θ )Ŝ(r e iφ ), where cosh(r) =
ing paradigm underlying Strawberry Fields.
p
1 + (s/2)2 , tan(θ ) = s/2, φ = −θ − sign(s)π/2. Equiv-
alently, we could have included the squeezing generator
Universal Gates for CV Quantum Computing x̂ p̂ + p̂ x̂ in place of the quadratic phase and still had a uni-
In discrete qubit systems, the notion of a universal gate set versal set.
has the following meaning: given a set of universal gates, We have just outlined an efficient method to construct
we can approximate an arbitrary unitary by composition of any gate of the form exp (−iH t), where the generator H is
said gates. In the CV setting, we have a similar situation: we a polynomial of the quadrature (or mode) operators. How
can approximate a broad set of generators – i.e., the Hamil- can this be used for quantum computations? As shown in
tonians appearing in Eq. (10) – by combining elements of Eq. (4), the eigenstates of the x̂ quadrature form an (or-
a CV universal gate set. However, unlike the qubit case, we thogonal) basis for representing qumode states. Thus, these
do not try to approximate all conceivable unitaries. Rather, states are often taken as the computational basis of the CV
we seek to create all generators that are a polynomial func- model (though other choices are also available). By apply-
tion of the quadrature (or mode) operators of the system ing gates as constructed above and performing measure-
[48, 49]. Remember that generators of second-degree or ments in this computational basis (i.e., homodyne measure-
lower belong to the class of Gaussian operations, while all ments), we can carry out a CV quantum computation. One
higher degrees are non-Gaussian. primitive example [49] is to compute the product of two
We can create a higher-order generator out of lower- numbers – whose values are stored in two qumode regis-
order generators Â and B̂ by using the following two con- ters – and store the result in the state of a third qumode.
catenation identities [48]: Consider the generator x̂ 1 x̂ 2 p̂3 /ħ
h, which will cause the “po-
sition” operators to evolve according to
e−i Âδt e−i B̂δt e i Âδt e i B̂δt = e[Â,B̂]δt + O(δt 3 ), (A1a)
2
x̂ 1 → x̂ 1 , x̂ 2 → x̂ 2 , x̂ 3 → x̂ 3 + x̂ 1 x̂ 2 t. (A3)
e i Âδt/2 e i B̂δt/2 e i B̂δt/2 e i Âδt/2 = e i(Â+B̂)δt + O(δt 2 ). (A1b)
If we have two second-degree generators, such as û = A measurement in the computational basis of mode 3 will
x̂ 2 + p̂2 (the generator for the rotation gate) and ŝ = reveal the value of the product x 1 x 2 . Note that these encod-
x̂ p̂ + p̂ x̂ (the generator for the squeezing gate), and a ings of classical continuous degress of freedom into quan-
third-degree (or higher) generator, such as x̂ 3 (the gen- tum registers allows for the generalization of neural net-
erator for the cubic phase gate), we can easily construct works into the quantum regime [66]. In the following sec-
generators of all higher-degree polynomials, e.g., x̂ 4 = tions we show how more complicated algorithms are con-
h2 ). This reasoning can be extended by
−[x̂ 3 , [x̂ 3 , û]]/(18ħ structed.
induction to any finite-degree polynomial in x̂ and p̂ (equiv-
alently, in â and â† ) [48, 49]. Multiport Gate Decompositions
In the above argument, it is important that at least one One important set of gates for which it is critical to de-
of the generators is third-degree or higher. Indeed, com- rive a decomposition in terms of universal gates is the set of
mutators of second-degree polynomials of x̂ and p̂ are also multiport interferometers. A multiport interferometer, rep-
second-degree polynomials and thus their composition us-
resented by a unitary operator Û acting on N modes, will
ing Eq. (A1) cannot generate higher-order generators. The
map (in the Heisenberg picture) the annihilation operator
claim can be easily extended to N -mode systems and mul-
of mode i into a linear combination of all other modes
tivariate polynomials of the operators
X
âi → Û † âi Û = (Ui j â j ). (A4)
r̂ = (x̂ 1 , . . . , x̂ N , p̂1 , . . . , p̂N ). (A2)
j
Combining single-mode universal gates (including at least In order to preserve the commutation relations of different
one of third-degree or higher) with some multimode in- modes, the matrix U must also be unitary U U † = U † U = IN .
teraction, e.g., the beamsplitter interaction generated by Note that every annihilation operator is mapped to a linear
b̂i, j = p̂i x̂ j − x̂ i p̂ j , we can construct arbitrary-degree poly- combination of annihilation operators and thus annihilation
nomials of the quadrature operators of an N -mode system. and creation operators are not mixed. Because of this, mul-
With the above discussion in mind, we can combine tiport interferometers generate all passive (in the sense that
the set of single-qumode gates generated by {x̂ i , x̂ i2 , x̂ i3 , ûi } no photons are created or destroyed) linear optics transfor-
and the two-mode gate generated by b̂i, j for all pairs of mations. In a pioneering work by Reck et al. [67], it was
modes into a universal gate set. The first, third and fourth shown that any multiport interferometer can be realized us-
single-mode generators correspond to the displacement, cu- ing N (N − 1)/2 beamsplitter gates distributed over 2N − 3
14
layers. Recently this result was improved upon by Clements of a quantum N mode state ρ:
et al. [63], who showed that an equivalent decomposition
can be achieved with the same number of beamsplitters but D̂(ξ) = exp (ir̂Ωξ) , χ(ξ) = 〈 D̂(ξ)〉ρ , (A6)
using only N + 1 layers.
where ξ ∈ R2N . We can now consider the Fourier transform
Gaussian Operations of the characteristic function to obtain the Wigner function
of the state ρ̂
As mentioned in the previous subsections, generators
d 2N ξ
Z
that are at most quadratic remain closed under composi-
tion of their associated unitaries. In the Heisenberg picture W (r) = exp(−irΩξ)χ(ξ). (A7)
R2N
(2π)2N
these quadratic generators will produce all possible linear
transformations between the quadrature (or mode) opera- The 2N real arguments r of the Wigner function are the
tors, eigenvalues of the quadrature operators from r̂.
The above recipe maps an N -mode quantum state liv-
X
âi → Ui j â j + Vi j â†j . (A5)
ing in a Hilbert space to the real symplectic space K :=
j
(R2N , Ω), which is called phase space. The Wigner function
These operations are known as Gaussian operations. All is an example of a quasiprobability distribution. Like a prob-
gates in Table III are Gaussian operations except for the cu- ability distribution over this phase space, the Wigner func-
bic phase gate. Pure Gaussian states are the set of states tion is normalized to one; however, unlike a probability dis-
that can be obtained from the (multimode) vacuum state tribution, it may take negative values. Gaussian states have
by Gaussian operations [50, 68]. Mixed Gaussian states the special property that their characteristic function (and
are obtained by applying Gaussian operations to thermal hence their Wigner function) is a Gaussian function of the
states or marginalizing pure Gaussian states. Analogous to variables r. In this case, the Wigner function takes the form
Gaussian states, we can also define Gaussian measurements
exp − 21 (r − r̄)V−1 (r − r̄)

as the set of measurements whose positive-operator valued
W (r) = p (A8)
measure (POVM) elements can be obtained from vacuum (2π)N detV
via Gaussian transformations. Homodyne and heterodyne
measurements are prominent examples of Gaussian mea- where r̄ = 〈r̂〉ρ = Tr (r̂ρ̂) is the displacement or mean vector
surements, whereas photon counting and photodetection and Vi j = 12 〈∆ri ∆r j + ∆ri ∆r j 〉ρ with ∆r̂ = r̂ − r̄. Note that
are prominent examples of non-Gaussian measurements. the only pure states that have non-negative Wigner func-
An important result for the CV formalism is that Gaussian tions are the pure Gaussian states [58].
quantum computation, i.e., computation that occurs with Each type of Gaussian state has a specific form of co-
Gaussian states, operations and measurements, can be effi- variance matrix V and mean vector r̄. For the single-mode
ciently simulated on a classical computer (this is the foun-
vacuum state, we have V = ħ2h I2 and r̄ = (0, 0) T . A ther-
dation for the Gaussian backend in Strawberry Fields). This
mal state (Eq. (8)) has the same (zero) displacement but
result is the CV version [69] of the Gottesman-Knill theorem
a covariance matrix V = (2n̄ + 1) ħ2h I2 , where n̄ is the mean
of discrete-variable quantum information [49]. Hence we
photon number. A coherent state (Eq. (6)), obtained by
need non-Gaussian operations in order to achieve quantum
displacing vacuum, has the same V asq the vacuum state
supremacy in the CV model. Interestingly, even in the re-
stricted case where all states and gates are Gaussian, with but a nonzero displacement vector r̄ = 2 ħ2h (Re(α), Im(α)).
only the final measurements being non-Gaussian, there is Lastly, a squeezed state (Eq. (7)) has zero displacement
strong evidence that such a circuit cannot be efficiently sim- and covariance matrix V = ħ2h diag(e−2r , e2r ). In the limit
ulated classically [10, 70]. More discussion and example r → ∞, the squeezed state’s variance in the x̂ quadrature
code for this situation (known as Gaussian boson sampling) becomes zero and the state becomes proportional to the x̂-
is provided in Appendix C. eigenstate |x〉 with eigenvalue 0. Consistent with the un-
certainty principle, the squeezed state’s variance in p̂ blows
Symplectic Formalism up. We can also consider the case r → −∞, where we find
a state proportional to the eigenstate |p〉 of the p̂ quadrature
In this section we review the symplectic formalism which with eigenvalue 0. In the laboratory and in numerical simu-
lies at the heart of the Gaussian backend of Strawberry lation we must approximate every quadrature eigenstate us-
Fields. The symplectic formalism is an elegant and compact ing a finitely squeezed state (being careful that the variance
description of Gaussian states in terms of covariance matri- of the relevant quadrature is much smaller than any other
ces and mean vectors [50, 68]. To begin, the commutation uncertainty relevant to the system). Any other quadrature
relations of the 2N position and momentum operators of eigenstate can be obtained from the x = 0 eigenstate by
Eq. (A2) can be easily summarized as [r̂i , r̂ j ] = ħ
hiΩi j where applying suitable displacement and rotation operators. Fi-
0 I

Ω = −IN 0N is the symplectic matrix. Using the symplectic nally, note that Gaussian operations will transform the vec-
matrix, we can define the Weyl operator D(ξ) (a multimode tor of means via affine transformations and the covariance
displacement operator) and the characteristic function χ(ξ) matrix via congruence transformations; for a detailed dis-
15
cussion of these transformations, see Sec. 2 of [50].

Given a 2N × 2N real symmetric matrix, how can we
check that it is a valid covariance matrix? If it is valid,
which operations (displacement, squeezing, multiport in-
terferometers) should be performed to prepare the corre-
sponding Gaussian state? To answer the first question: a
2N × 2N real symmetric matrix Ṽ corresponds to a Gaus-
sian quantum state if and only if Ṽ + i ħ2h Ω ≥ 0 (the matrix
inequality is understood in the sense that the eigenvalues
of the quantity Ṽ + i ħ2h Ω are nonnegative). The answer to
the second question is provided by the Bloch-Messiah reduc-
tion [60–62]. This reduction shows that any N -mode Gaus-
sian state (equivalently any covariance matrix and vector
of means) can beN constructed by starting with a product of
N thermal states i ρi (n̄i ) (with potentially different mean
photon numbers), then applying a multiport interferometer
Û , followed by single-mode squeezing operations i Si (zi ),
N
followed by another multiport V̂ , followed by single-mode

displacement operations i Di (αi ). Explicitly,
N

O
ρGaussian = Ŵ ρi (n̄i ) Ŵ † , (A9)
i

O O
Ŵ = Di (αi ) V̂ Si (zi ) Û . (A10)
i i
Note that if the Gaussian state is pure (which happens

if and only if det(V) = (ħ h/2)2N ), the occupation number
of the thermal states in the Bloch-Messiah decomposition
are all zero and the first interferometer will turn the vac-
uum to vacuum again. Thus for pure Gaussian states we
need only generate N single-mode squeezed states and send
them through a single multiport interferometer V̂ before
displacing. For a recent discussion of this decomposition
see Ref. [71, 72]. More generally, the occupation numbers
of the different thermal states in Eq. (A9) ni = (νi − 1)/2
can be obtained by calculating the symplectic eigenvalues
νi of the covariance matrix V . The symplectic eigenvalues
come in pairs and are just the standard eigenvalues of the
matrix |i(2/ħh)ΩV| where the modulus is understood in the
operator sense (see Sec. II.C.1. of Ref. [50]).
16
Appendix B: Strawberry Fields Operations

In this Appendix, we present a complete list of the CV states, gates, and measurements available in Strawberry Fields.
Operation Name Definition
Vacuum() Vacuum state The vacuum state |0〉, representing zero photons
Coherent(a) Coherent state A displaced vacuum state, |α〉 = D(α) |0〉, α ∈ C
A squeezed vacuum state, |z〉 = S(z) |0〉,

Squeezed(r,phi) Squeezed state where z = r e iφ , r, φ ∈ R, r ≥ 0, φ ∈ [0, 2π)
A squeezed then displaced vacuum state,

DisplacedSqueezed(a,r,phi) Displaced squeezed state |α, z〉 = D(α)S(z) |0〉
P∞
ρ(n̄) = n=0 [nn /(1 + n)n+1 ]|n〉〈n|,
Thermal(n) Thermal state where n̄ ∈ R+ is the mean photon number
Fock(n)* Fock state or number state |n〉, where n ∈ N0 represents the photon number
p1 (|α〉 + exp(iπp) |−α〉), where p = 0, 1 gives an

N
Catstate(a,p)* Cat state even/odd cat state and N is the normalization
Prepare an arbitrary multi-mode pure state, repre-

Ket(x)* Arbitrary Fock-basis ket sented by array x, in the Fock basis.
Prepare an arbitrary multi-mode mixed state, repre-

DensityMatrix(x)* Arbitrary Fock basis state
sented by a density matrix array x in the Fock basis.
Table V: State preparations available in Strawberry Fields. Those indicated with an asterisk (*) are non-Gaussian.

D(α) = exp(αâ† − α∗ â)
Dgate(a) Displacement gate Æ
= exp −i(Re(α)p̂ − Im(α)x̂) 2/ħ h , α∈C
p
Xgate(x) Position displacement gate X (x) = D(x/ 2ħ h) = exp (−i x p̂/ħ
h), x ∈ R
p
Zgate(p) Momentum displacement gate Z(p) = D(ip/ 2ħ h) = exp (ip x̂/ħ
h), p ∈ R
2

S(z) = exp (z ∗ â2 − z â† )/2 , where
Sgate(r,phi) Squeezing gate z = r exp (iφ), r, φ ∈ R, r ≥ 0, φ ∈ [0, 2π)

Rgate(theta) Rotation gate R(θ ) = exp iθ â† â , where θ ∈ [0, 2π)

Fouriergate() Fourier gate F = R(π/2) = exp i(π/2)â† â

Pgate(s) Quadratic phase gate P(s) = exp is x̂ 2 /(2ħh) , where s ∈ R

Vgate(g)* Cubic phase gate V (γ) = exp iγx̂ 3 /(3ħ h) , where γ ∈ R
2
Kgate(k)* Kerr interaction gate K(κ) = exp iκ â† â , where κ ∈ R
Table VI: Single mode gate operations available in Strawberry Fields. Those indicated with an asterisk (*) are non-Gaussian
and can only be used with a backend that uses the Fock representation.
17

B(θ , φ) = exp θ (e iφ â1 â2† − e−iφ â1† â2 ) ,

BSgate(theta,phi) Beamsplitter where the transmissivity and reflectivity

amplitudes are t = cos θ , r = e iφ sin θ
S2 (z) = exp z ∗ â1 â2 − z â1† â2† ,

S2gate(r,p) Two-mode squeezing gate
where z = r exp (iφ)
CXgate(s) Controlled-X or addition gate C X (s) = exp (−is x̂ 1 p̂2 /ħ
h), s ∈ R
CZgate(s) Controlled phase shift gate C Z(s) = exp (is x̂ 1 x̂ 2 /ħ
h), s ∈ R
C K(κ) = exp iκâ1 â1 â2† â2 , κ ∈ R
†

CKgate(k)* Controlled Kerr interaction gate
Table VII: Two-mode gate operations available in Strawberry Fields.

Prepares an arbitrary N -mode Gaussian state de-
Gaussian(cov,mu) Gaussian state preparation fined by covariance matrix V ∈ R2N ×2N and means
vector µ ∈ R2N using the Williamson decomposition
Applies an N -mode Gaussian transformation de-
GaussianTransform(S) Gaussian transformation fined by symplectic matrix S ∈ R2N ×2N using the
Bloch-Messiah decomposition
Applies an N -mode interferometer defined by
Interferometer(U) Multi-mode linear interferometer unitary matrix U ∈ CN ×N using the Clements
decomposition
Table VIII: Multi-mode Gaussian decompositions available in Strawberry Fields.

Projects the state onto |x φ 〉〈x φ |
MeasureHomodyne(phi) Homodyne measurement
where x̂ φ = cos φ x̂ + sin φp̂
Projects the state onto the coherent states, sam-
MeasureHeterodyne() Heterodyne measurement
pling from the joint Husimi distribution π1 〈α|ρ|α〉
MeasureFock()* Photon counting Projects the state onto |n〉〈n|
Table IX: Measurement operations available in Strawberry Fields. Those indicated with an asterisk (*) are non-Gaussian and
can only be used with a backend that uses the Fock representation.
p
where cosh(r) = 1 + (s/2)2 , tan(θ ) = s/2, and φ =
−sign(s) π2 − θ .
Gate Decompositions
In addition, the Strawberry Fields frontend can be used to Two-mode squeeze gate
provide decompositions of certain compound gates. The The two-mode squeeze gate S2gate(z) is decomposed
following gate decompositions are currently supported. into a combination of beamsplitters and single-mode
squeezers
Quadratic phase gate S2 (z) = B † (π/4, 0) [S(z) ⊗ S(−z)] B(π/4, 0)
The quadratic phase shift gate Pgate(s) is decomposed

into a squeezing and a rotation,
P(s) = R(θ )S(r e iφ ),

18
Controlled addition gate

The controlled addition or controlled-X gate CXgate(s) is
decomposed into a combination of beamsplitters and single-
mode squeezers
C X (s) = B(φ, 0) [S(r) ⊗ S(−r)] B(π/2 + φ, 0),
where sin(2φ) = −1/cosh(r), cos(2φ) = − tanh(r), and

sinh(r) = −s/2.
Controlled phase gate

The controlled phase shift gate CZgate(s) is decomposed
into a controlled addition gate, with two rotation gates act-
ing on the target mode,

C Z(s) = [I ⊗ R(π/2)] C X (s) I ⊗ R(π/2)†
19
14 # 50-50 beamsplitter
Appendix C: Quantum Algorithms 15 BS = BSgate(pi/4, 0)
16
17 # maximally entangled states

In this Appendix, we present full example code for sev-
18 Squeezed(-2) | alice
eral important algorithms, subroutines, and use-cases for 19 Squeezed(2) | bob
Strawberry Fields. These examples are presented in more 20 BS | (alice, bob)
detail in the online documentation located at strawberry- 21
fields.readthedocs.io. 22 # Alice performs the joint measurement
23 # in the maximally entangled basis
Quantum Teleportation 24 BS | (psi, alice)
25 MeasureX | psi
Quantum teleportation — sometimes referred to as state
26 MeasureP | alice
teleportation to avoid confusion with gate teleportation — 27
is the reliable transfer of an unknown quantum state across 28 # Bob conditionally displaces his mode
spatially separated qubits or qumodes, through the use of a 29 # based on Alice's measurement result
classical transmission channel and quantum entanglement 30 Xgate(scale(psi, sqrt(2))) | bob
[73]. Considered a fundamental quantum information pro- 31 Zgate(scale(alice, sqrt(2))) | bob
tocol, it has applications ranging from quantum commu- 32 # end circuit
nication to enabling distributed information processing in 33
quantum computation [74]. 34 state = eng.run("gaussian")

In general, all quantum teleportation circuits work on 35 # view Bob's output state and fidelity
the same basic principle. Two distant operators, Alice and 36 print(q[0].val,q[1].val)
37 print(state.displacement([2]))
Bob, share a maximally entangled quantum state (in dis-
38 print(state.fidelity_coherent([0,0,1+0.5j]))
crete variables, any one of the four Bell states, and in CVs,
a maximally entangled state for a fixed energy), and have Codeblock 15: Quantum teleportation of a coherent state in
access to a classical communication channel. Alice, in pos- Strawberry Fields. Note: while the CV quantum teleportation
session of an unknown state which she wishes to transfer algorithm relies on infinitely squeezed resource states, in
to Bob, makes a joint measurement of the unknown state practice a squeezing magnitude of -2 (∼ 18 dB) is sufficient.
and her half of the entangled state, by projecting onto the
Bell or quadrature basis. By transmitting the results of her
measurement to Bob, Bob is then able to transform his half Gate Teleportation
of the entangled state to an accurate replica of the origi-
In the quantum state teleportation algorithm discussed in
nal unknown state, by performing a conditional phase flip
the previous section, the quantum state is transferred from
(for qubits) or displacement (for qumodes) [75]. The CV
the sender to the receiver. However, quantum teleportation
teleportation circuit is shown in Fig. 3.
can be used in a much more powerful manner, by simulta-
neously transforming the teleported state — this is known
q0 : |ψ〉 • |x〉 〈x| = m1
as gate teleportation.
BS
• |p〉 〈p| = m2 In gate teleportation, rather than applying a quantum
q1 : |0〉 P
unitary directly to the state prior to teleportation, the
BS
q2 : |0〉 P |ψ〉 unitary is applied indirectly, via the projective measure-
X Z
ment of the first subsystem onto a particular basis. This
FIG. 3: State teleportation of state |ψ〉 from mode q0 to measurement-based approach provides significant advan-
mode q2 . tages over applying unitary gates directly, for example by
reducing resources, and in the application of experimen-
tally hard-to-implement gates [74]. In fact, gate teleporta-
tion forms a universal quantum computing primitive, and is
1 import strawberryfields as sf
2 from strawberryfields.ops import * a precursor to cluster state models of quantum computation
3 from strawberryfields.utils import scale [76].
4 from numpy import pi, sqrt First described by Gottesman and Chuang [77] in the case
5 of qubits, gate teleportation was generalised for the CV case
6 eng,q = sf.Engine(3) by Bartlett and Munro [78] (see Fig. 4). In an analogous
7 process to the discrete-variable case, we begin with the al-
8 with eng: gorithm for local state teleportation. Unlike the spatially-
9 psi, alice, bob = q[0], q[1], q[2] separated quantum state teleportation we considered in the
10
previous section, local teleportation can transport the state
11 # state to be teleported:
using only two qumodes; the state we are teleporting is en-
12 Coherent(1+0.5j) | psi
13
tangled directly with the squeezed vacuum state in the mo-
mentum space through the use of a controlled phase gate.
20
q0 : |n1 〉 R n̂
BS BS BS
q1 : |n2 〉 R n̂
BS BS
q2 : |n3 〉 R n̂
BS BS BS
q3 : n4 R n̂
FIG. 4: Gate teleportation of a quadratic phase gate P onto FIG. 5: 4-mode boson sampling.
input state |ψ〉.
surements, all ‘stacked vertically’ (i.e., coupled to each con-

To recover the teleported state exactly, we must perform secutive qumode via a conditional phase gate). It is from
Weyl-Heisenberg corrections to the second mode; here, that this primitive that the model of cluster state quantum com-
would be F † X (m)† , where m is the correction based on the putation can be derived [76].
Homodyne measurement. However, for convenience and
simplicity, it is common to simply write the circuit without Boson Sampling
the corrections applied explicitly. Introduced by Aaronson and Arkhipov [9], boson sam-
pling presented a slight deviation from the general ap-
1 import strawberryfields as sf proach in quantum computation. Rather than presenting a
2 from strawberryfields.ops import * theoretical model of universal quantum computation (i.e.,
3
a framework that enables quantum simulation of any arbi-
4 eng, q = sf.Engine(4)
trary Hamiltonian [51]), boson sampling-based devices are
5
6 with eng: instead an example of an intermediate quantum computer,

7 # create initial states designed to experimentally implement a computation that
8 Squeezed(0.1) | q[0] is intractable classically [79–82].
9 Squeezed(-2) | q[1] Boson sampling proposes the following quantum linear
10 Squeezed(-2) | q[2] optics scheme. An array of single photon sources is set up,
11 designed to simultaneously emit single photon states into a
12 # apply the gate to be teleported multimode linear interferometer; the results are then gen-
13 Pgate(0.5) | q[1] erated by sampling from the probability of single photon
14
measurements from the output of the linear interferometer.
# conditional phase entanglement
Fock states, |ψ〉 =
15
For example, consider N single photon
16 CZgate(1) | (q[0], q[1])
|m1 , m2 , . . . , mN 〉, composed of b = i mi photons, incident
P
17 CZgate(1) | (q[1], q[2])
18
on an N -mode linear interferometer described by the uni-
19 # projective measurement onto tary U acting on the input mode creation and annihilation
20 # the position quadrature operators. It was shown that the probability of detecting n j
21 Fourier.H | q[0] photons at the jth output mode is given by [9]
22 MeasureX | q[0]
23 Fourier.H | q[1] |Per(Ust )|2
24 MeasureX | q[1] Pr(n1 , . . . , nN ) = ; (C1)
m1 ! · · · mN !n1 ! · · · nN !
25 # compare against the expected output
26 # X(q1).F.P(0.5).X(q0).F.|z> i.e., the sampled single photon probability distribution is
27 # not including the corrections
proportional to the permanent of Ust , a submatrix of the
28 Squeezed(0.1) | q[3]
29 Fourier | q[3] interferometer unitary, dependent upon the input and out-
30 Xgate(q[0]) | q[3] put Fock states. Whilst the determinant can be calculated
31 Pgate(0.5) | q[3] efficiently on classical computers, calculation of the per-
32 Fourier | q[3] manent belongs to the computational complexity class #P-
33 Xgate(q[1]) | q[3] Hard problems [83], which are strongly believed to be clas-
34 # end circuit sically hard to calculate. This implies that simulating bo-
35 son sampling is an intractable task for classical comput-
36 state = eng.run("gaussian") ers, providing an avenue for the demonstration of quantum
37 print(state.reduced_gaussian([2])) supremacy.
38 print(state.reduced_gaussian([3])) Continuous-variable quantum computation is ideally
Codeblock 16: Gate teleportation of a quadratic phase gate suited to the simulation and demonstration of the boson
in Strawberry Fields. sampling scheme, due to its grounding in quantum optics.
In quantum linear optics, the multimode linear interferom-
Note that additional gates can be added to the gate tele- eter is commonly decomposed into two-mode beamsplitters
portation scheme described above simply by introducing (BSgate) and single-mode phase shifters (Rgate) [67], al-
additional qumodes with the appropriate projective mea- lowing for a straightforward translation into a CV quantum
21
q0 : |0〉 S R n̂ scalability, due to its dependence on an array of simultane-

BS BS BS ously emitting single photon sources. Currently, most phys-
q1 : |0〉 S R n̂
BS BS ical implementations of boson sampling make use of a pro-
q2 : |0〉 S R n̂ cess known as Spontaneous Parametric Down-Conversion
BS BS BS (SPDC) to generate the single-photon source inputs. How-
q3 : |0〉 S R n̂
ever, this method is non-deterministic – as the number of
FIG. 6: 4-mode Gaussian boson sampling. modes in the apparatus increases, the average time required
until every photon source emits a simultaneous photon in-
creases exponentially.
circuit. In order to allow for arbitrary linear unitaries on m In order to simulate a deterministic single photon source
imput modes, we must have a minimum of m + 1 columns array, several variations on boson sampling have been pro-
in the beamsplitter array [63]. posed; the most well-known being scattershot boson sam-
pling [70]. However, a recent boson sampling variation by
1 import numpy as np [10] negates the need for single photon Fock states alto-
2 import strawberryfields as sf gether, by showing that incident Gaussian states – in this
3 from strawberryfields.ops import * case, single mode squeezed states – can produce problems
4 in the same computational complexity class as boson sam-
5 # initialise the engine and register pling. Even more significantly, this mitigates the scalabil-
6 eng, q = sf.Engine(4) ity problem with single photon sources, as single mode
7
squeezed states can be simultaneously generated experi-
8 with eng: mentally.
9 # prepare the input fock states
10 Fock(1) | q[0] With an input ensemble of N single-mode squeezed states
11 Fock(1) | q[1] with squeezing parameter z = r ∈ R, incident on a linear-
12 Vac | q[2] interferometer described by unitary U, it can be shown
13 Fock(1) | q[3] that the probability of detecting an output photon pattern
14 (n1 , . . . , nN ), where nk ∈ {0, 1}, is given by [10]
15 # rotation gates
Rgate(0.5719) Haf[(U L tanh(r )U T ) ]2

16
i i S
17 Rgate(-1.9782) Pr(n1 , . . . , nN ) = , (C2)
i cosh(ri )
Q
18 Rgate(2.0603)
19 Rgate(0.0644)
20 where S denotes the subset of modes where a photon was
21 # beamsplitter array detected and Haf[·] is the Hafnian [84, 85]. That is, the
22 BSgate(0.7804, 0.8578) | (q[0], q[1]) sampled single photon probability distribution
L is propor-
23 BSgate(0.06406, 0.5165) | (q[2], q[3]) tional to the hafnian of a submatrix of U i tanh(ri )U T .
24 BSgate(0.473, 0.1176) | (q[1], q[2]) The hafnian is known to be a generalization of the per-
25 BSgate(0.563, 0.1517) | (q[0], q[1]) manent, and can be used to count the number of perfect
26 BSgate(0.1323, 0.9946) | (q[2], q[3]) matchings of an arbitrary graph [86]. The formula above
27 BSgate(0.311, 0.3231) | (q[1], q[2])
can be generalized to pure Gaussian states with finite dis-
28 BSgate(0.4348, 0.0798) | (q[0], q[1])
29 BSgate(0.4368, 0.6157) | (q[2], q[3]) placements by using the loop hafnian function which counts
30 # end circuit the number of perfect matchings of graphs with loops[56].
31 Since any algorithm that could calculate the hafnian could
32 # run the engine also calculate the permanent, it follows that calculating the
33 state = eng.run("fock", cutoff_dim=7) hafnian remains a classically hard problem; indeed, the best
34 known classical algorithm for the calculation of a hafnian of
35 # extract the joint Fock probabilities an arbitrary symmetric complex matrix of size N scales like
36 probs = state.all_fock_probs() O(N 3 2N /2 ) [87]. The hardness of approximate GBS, under
37
imperfections such as loss, is a subject of current research
38 # print the joint Fock state probabilities [88, 89].
39 print(probs[1,1,0,1])
40 print(probs[2,0,0,1])
1 import numpy as np
Codeblock 17: 4-mode boson sampling example in 2 import strawberryfields as sf
Strawberry Fields. Parameters are chosen arbitrarily. 3 from strawberryfields.ops import *
4
5 # initialise the engine and register

Gaussian Boson Sampling 6 eng, q = sf.Engine(4)
7
While boson sampling allows the experimental imple-
8 with eng:
mentation of a sampling problem that is hard to simulate 9 # prepare the input squeezed states
classically, one of the setbacks in experimental setups is
22
|0〉 S • Z • V p̂
10 S = Sgate(1)
11 S | q[0] |0〉 S • • V p̂
12 S | q[1]
13 S | q[2] |0〉 S V Z • • Z p̂
14 S | q[3]
15 |0〉 S Z Z Z p̂
16 # rotation gates
17 Rgate(0.5719) | q[0] FIG. 7: 4-mode instantaneous quantum polynomial (IQP)
18 Rgate(-1.9782) | q[1] CV circuit example.
19 Rgate(2.0603) | q[2]
20 Rgate(0.0644) | q[3]
21
22 # beamsplitter array however, the IQP protocol was not constructed with quan-
23 BSgate(0.7804, 0.8578) | (q[0], q[1]) tum linear optics in mind. Nevertheless, the IQP model
24 BSgate(0.06406, 0.5165) | (q[2], q[3]) was recently extended to the CV formulation of quantum
25 BSgate(0.473, 0.1176) | (q[1], q[2]) computation by Douce et al. [17], taking advantage of
26 BSgate(0.563, 0.1517) | (q[0], q[1]) the ability to efficiently prepare large resource states, and
27 BSgate(0.1323, 0.9946) | (q[2], q[3]) the higher efficiencies afforded by homodyne detection.
28 BSgate(0.311, 0.3231) | (q[1], q[2]) Furthermore, the computational complexity results of the
29 BSgate(0.4348, 0.0798) | (q[0], q[1]) discrete-variable case apply equally to the so-called CV-IQP
30 BSgate(0.4368, 0.6157) | (q[2], q[3])
model, assuming a specific input squeezing parameter de-
31 # end circuit
32
pendent on the circuit size. The CV-IQP model is defined as
33 # run the engine follows:
34 state = eng.run("gaussian")
35
1. inputs are squeezed momentum states |z〉, where z =
36 # Fock states to measure at output r ∈ R and r < 0;
37 measure_states = [[0,0,0,0], [1,1,0,0],
,→ [0,1,0,1], [1,1,1,1], [2,0,0,0]] 2. the unitary transformation is diagonal in the x̂
38
quadrature basis, by randomly selecting from the set
39 # extract the probabilities of calculating of gates U = {Z(p), C Z(s), V (γ)};
40 # several different Fock states at the output
41 # and print them to the terminal 3. and finally, homodyne measurements are performed
42 for i in measure_states: on all modes in the p̂ quadrature.
43 prob = state.fock_prob(i)
44 print("|{}>: {}".format("".join(str(j) for
,→ j in i), prob)) 1 import strawberryfields as sf
3
Codeblock 18: 4-mode Gaussian boson sampling example
4 # initialise the engine and register
in Strawberry Fields. Parameters are chosen arbitrarily.
5 eng, q = sf.Engine(4)
6
Instantaneous Quantum Polynomial (IQP) 7 with eng:

8 # prepare the input squeezed states
Like boson sampling and Gaussian boson sampling before 9 S = Sgate(-1)
it, the instantaneous quantum polynomial (IQP) protocol is 10 S | q[0]
a restricted, non-universal model of quantum computation, 11 S | q[1]
designed to implement a computation that is classically in- 12 S | q[2]
tractable and verifiable by a remote adjudicator. First intro- 13 S | q[3]
duced by Shepherd and Bremner [90], IQP circuits are de-
14
15 # gates diagonal in the x quadrature

fined by the following conditions: (i) there are N inputs in 16 CZgate(0.5719) | (q[0], q[1])
state |0〉 acted on by Hadamard gates, (ii) the resulting com- 17 Vgate(-1.9782) | q[2]
putation is diagonal in the computational
p basis by randomly 18 Zgate(2.0603) | q[3]
selecting from the set U = {R(π/4), C Z} (hence the term 19 Zgate(-0.7804) | q[0]
‘instantaneous’; since all gates commute, they can be ap- 20 Zgate(0.8578) | q[2]
plied in any temporal order), and (iii) the output qubits are 21 CZgate(1.321) | (q[1], q[2])
measured in the Pauli-X basis. Efficient classical sampling 22 Zgate(0.473) | q[3]
of the resulting probability distribution H ⊗N U H ⊗N |0〉⊗N — 23 CZgate(0.9946) | (q[0], q[2])
even approximately [13] or in the presence of noise [91] — 24 Zgate(0.1223) | q[3]
has been shown to be #P-hard, and would result in the col- 25 Vgate(0.6157) | q[0]
Vgate(0.3110) | q[1]
lapse of the polynomial hierarchy to the third level [12, 92].
26
27 Zgate(-1.243) | q[2]
Unlike boson sampling and Gaussian boson sampling,
23
the Lie product formula, we find that

28 # end circuit
29
Jt † † Ut 2
30 # run the engine e = exp −i (â1 â2 + â2 â1 ) exp −i
iH t
n̂
31 eng.run("fock", cutoff_dim=5) k 2k 1
k
Ut 2 Ut Ut
exp −i n̂ exp i n̂1 exp i n̂2
2k 2 2k 2k
Codeblock 19: 4-mode instantaneous quantum polynomial
(IQP) example in Strawberry Fields. + O t 2 /k , (C4)

Moreover, the resulting probability distributions have where O t 2 /k is the order of the error term, derived from
been shown to be given by integrals of oscillating functions the Lie product formula. Comparing this to the form of var-
in large dimensions, which are believed to be intractable ious gates in the CV circuit model, we can write this as the
to compute by classical computers. This leads to the result product of beamsplitters, Kerr gates, and rotation gates:
that even in the case of finite squeezing and reduced mea-
e iH t = {BS (θ , φ) [K(r)R(−r) ⊗ K(r)R(−r)]}k + O t 2 /k

surement precision, approximate sampling from the output
CV-IQP model remains classically intractable [17, 93]. (C5)
where θ = −J t/k, φ = π/2, and r = −U t/2k. Using J = 1,

Hamiltonian Simulation U = 1.5, k = 20, and a timestep of t = 1.086, this can be
easily implemented in Strawberry Fields using only 20×3 =
The simulation of atoms, molecules and other biochemi- 60 gates.
cal systems is another application uniquely suited to quan-
tum computation. For example, the ground state energy 1 import strawberryfields as sf
of large systems, the dynamical behaviour of an ensemble 2 from strawberryfields.ops import *
of molecules, or complex molecular behaviour such as pro- 3 from numpy import pi
tein folding, are often computationally hard or downright 4
impossible to determine via classical computation or exper- 5 # initialise the engine and register
imentation [94, 95]. 6 eng, q = sf.Engine(2)
7
In the discrete-variable qubit model, efficient methods of 8 # set the Hamiltonian parameters
Hamiltonian simulation have been discussed at length, pro- 9 J = 1 # hopping transition
viding several implementations depending on properties of 10 U = 1.5 # on-site interaction
the Hamiltonian, and resulting in a linear simulation time 11 k = 20 # Lie product decomposition
[96, 97]. Efficient implementations of Hamiltonian simula- ,→ terms
tion also exist in the CV formulation [98], with specific ap- 12 t = 1.086 # timestep
plication to Bose-Hubbard Hamiltonians (describing a sys- 13 theta = -J*t/k
tem of interacting bosonic particles on a lattice of orthog- 14 r = -U*t/(2*k)
onal position states [99]). As such, this method is ideally 15
16 with eng:
suited to photonic quantum computation.
17 # prepare the initial state
18 Fock(2) | q[0]
19
q0 : |2〉 K R K R 20 # Two node tight-binding
BS BS 21 # Hamiltonian simulation
q1 : |0〉 K R K R
22
FIG. 8: Bose-Hubbard Hamiltonian simulation for a 2-node

23 for i in range(k):
24 BSgate(theta, pi/2) | (q[0], q[1])
lattice, with k = 2. The error is of the order O (t 2 /k). 25 Kgate(r) | q[0]
26 Rgate(-r) | q[0]
27 Kgate(r) | q[1]
Considering a lattice composed of two adjacent nodes 28 Rgate(-r) | q[1]
characterised by Adjacency matrix A, the Bose-Hubbard 29 # end circuit
Hamiltonian is given by 30
31 # run the engine

XX 1 X 32 state = eng.run("fock", cutoff_dim=5)
H=J Ai j âi† â j + U n̂i (n̂i − 1) 33 # the output state probabilities
2
i j i 34 print(state.fock_prob([0,2]))
1 35 print(state.fock_prob([1,1]))
= J(â1† â2 + â2† â1 ) + U(n̂21 − n̂1 + n̂22 − n̂2 ), (C3) 36 print(state.fock_prob([2,0]))
2
where J represents the transition of the boson between Codeblock 20: Tight-binding Hamiltonian simulation for a
nodes, and U is the on-site interaction potential. Applying 2-node lattice in Strawberry Fields.
24
For more complex Hamiltonian CV decompositions, in-

cluding those with nearest-neighbour, see Kalajdzievski
et al. [98]. This decomposition is also implemented in
1 import strawberryfields as sf
the SFOpenBoson plugin [45], which provides an interface 2 from strawberryfields.ops import *
between OpenFermion, the quantum electronic structure 3 import tensorflow as tf
package, and Strawberry Fields. This allows arbitrary Bose- 4
Hubbard Hamiltonians, generated in OpenFermion, to be 5 eng, q = sf.Engine(1)
simulated using Strawberry Fields. 6
7 alpha = tf.Variable(0.1)
Optimization of Quantum Circuits 8 with eng:
9 Dgate(alpha) | q[0]
One of the unique features of Strawberry Fields is that it 10 state = eng.run("tf", cutoff_dim=7,
has been designed from the start to support modern com- ,→ eval=False)
putational methods like automatic gradients, optimization, 11
and machine learning by leveraging a TensorFlow backend. 12 # loss is probability for the Fock state n=1
We have already given an overview in the main text of how 13 prob = state.fock_prob([1])
these features are accessed in Strawberry Fields. We present 14 loss = -prob # negative sign to maximize prob
here a complete example for the optimization of a quantum 15
circuit. Our goal in this circuit is to find the Dgate param- 16 # Set up optimization
eter which leads to the highest probability of a n = 1 Fock 17 optimizer = tf.train.GradientDescentOptimizer(
,→ learning_rate=0.1)
state output. This simple baseline example can be straight-
18 minimize_op = optimizer.minimize(loss)
forwardly extended to much more complex circuits. As op-
timization is a key ingredient of machine learning, this ex-
19
ample can also serve as a springboard for more advanced
20 # Create TF Session and initialize variables
data-driven modelling tasks.
21 sess = tf.Session()
We note that optimization here is in a variational sense, 22 sess.run(tf.global_variables_initializer())
i.e., we choose an ansatz for the optimization by fixing the 23
discrete structure of our circuits (the selection and arrange- 24 # Carry out optimization
ment of specific states, gates, and measurements). The op- 25 for step in range(50):
timization then proceeds over the parameters of these op- 26 prob_val, _ = sess.run([prob, minimize_op])
erations. Finally, we emphasize that for certain circuits and 27 print("Value at step {}: {}".format(step,
objective functions, the optimization might be non-convex ,→ prob_val))
in general. Thus we should not assume (without proof) that
the solution obtained via a Strawberry Fields optimization Codeblock 21: Example of optimizing quantum circuit
is always the global optimum. However, it is often the case parameters in Strawberry Fields. This can be extended to a
in machine learning that local optima can still provide ef- machine learning setting by incorporating data and a more
fective solutions, so this may not be an issue, depending on sophisticated loss function.
the application.
[1] P. W. Shor. Polynomial-time algorithms for prime fac- algorithms. New Journal of Physics, 18(2):023023, 2016.
torization and discrete logarithms on a quantum com- doi:10.1088/1367-2630/18/2/023023.
puter. SIAM review, 41(2):303–332, 1999. doi: [7] E. Farhi, J. Goldstone, and S. Gutmann. A quan-
10.1137/S0036144598347011. tum approximate optimization algorithm. arXiv preprint
[2] L. K. Grover. A fast quantum mechanical algorithm for arXiv:1411.4028, 2014.
database search. In Proceedings of the twenty-eighth annual [8] E. Farhi and A. W. Harrow. Quantum supremacy through
ACM Symposium on Theory of Computing, pages 212–219. the quantum approximate optimization algorithm. arXiv
ACM, 1996. preprint arXiv:1602.07674, 2016.
[3] A. W. Harrow, A. Hassidim, and S. Lloyd. Quantum algorithm [9] S. Aaronson and A. Arkhipov. The computational complexity
for linear systems of equations. Physical Review Letters, 103 of linear optics. In Proceedings of the forty-third annual ACM
(15):150502, 2009. doi:10.1103/PhysRevLett.103.150502. symposium on Theory of computing, pages 333–342. ACM,
[4] J. Preskill. Quantum Computing in the NISQ era and beyond. 2011. doi:10.1145/1993636.1993682.
Quantum, 2:79, 2018. doi:10.22331/q-2018-08-06-79. [10] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen,
[5] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, C. Silberhorn, and I. Jex. Gaussian boson sam-
P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien. A variational pling. Physical Review Letters, 119:170501, 2017. doi:
eigenvalue solver on a photonic quantum processor. Nature 10.1103/PhysRevLett.119.170501.
Communications, 5, 2014. doi:10.1038/ncomms5213. [11] L. Chakhmakhchyan, R. Garcia-Patron, and N. J. Cerf. Boson
[6] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru- sampling with Gaussian measurements. Physical Review A,
Guzik. The theory of variational hybrid quantum-classical 96:032326, 2017. doi:10.1103/PhysRevA.96.032326.
25
[12] M. J. Bremner, R. Jozsa, and D. J. Shepherd. Classical simu- in Physics, 2:5, 2014. doi:10.3389/fphy.2014.00005.
lation of commuting quantum computations implies collapse [26] F. Neukart, D. Von Dollen, G. Compostella, C. Seidel,
of the polynomial hierarchy. In Proceedings of the Royal Soci- S. Yarkoni, and B. Parney. Traffic flow optimization using
ety of London A: Mathematical, Physical and Engineering Sci- a quantum annealer. Frontiers in ICT, 4:29, 2017. doi:
ences, page rspa20100301. The Royal Society, 2010. doi: 10.3389/fict.2017.00029.
10.1098/rspa.2010.0301. [27] H. Neven, V. S. Denchev, G. Rose, and W. G. Macready. Train-
[13] M. J. Bremner, A. Montanaro, and D. J. Shepherd. Average- ing a large scale classifier with the quantum adiabatic algo-
case complexity versus approximate simulation of commut- rithm. arXiv preprint arXiv:0912.0779, 2009.
ing quantum computations. Physical Review Letters, 117(8): [28] K. L. Pudenz and D. A. Lidar. Quantum adiabatic machine
080501, 2016. doi:10.1103/PhysRevLett.117.080501. learning. Quantum Information Processing, 12(5):2027–
[14] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, 2070, 2013. doi:10.1007/s11128-012-0506-4.
Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven. Char- [29] D. Crawford, A. Levit, N. Ghadermarzy, J. S. Oberoi, and
acterizing quantum supremacy in near-term devices. Nature P. Ronagh. Reinforcement learning using quantum Boltz-
Physics, 14(6):595, 2018. doi:10.1038/s41567-018-0124-x. mann machines. Quantum Information & Computation, 18
[15] S. Aaronson and L. Chen. Complexity-Theoretic Foundations (1-2):0051–0074, 2018. doi:10.26421/QIC18.1-2.
of Quantum Supremacy Experiments. In R. O’Donnell, edi- [30] M. H. Amin, E. Andriyash, J. Rolfe, B. Kulchytskyy, and
tor, 32nd Computational Complexity Conference (CCC 2017), R. Melko. Quantum boltzmann machine. Phys. Rev. X, 8:
volume 79 of Leibniz International Proceedings in Informat- 021050, May 2018. doi:10.1103/PhysRevX.8.021050.
ics (LIPIcs), pages 22:1–22:67. Schloss Dagstuhl–Leibniz- [31] D. Ristè, M. P. Da Silva, C. A. Ryan, A. W. Cross, A. D. Cór-
Zentrum fuer Informatik, 2017. ISBN 978-3-95977-040-8. coles, J. A. Smolin, J. M. Gambetta, J. M. Chow, and B. R.
doi:10.4230/LIPIcs.CCC.2017.22. Johnson. Demonstration of quantum advantage in machine
[16] C. Neill, P. Roushan, K. Kechedzhi, S. Boixo, S. V. Isakov, learning. npj Quantum Information, 3(1):16, 2017. doi:
V. Smelyanskiy, A. Megrant, B. Chiaro, A. Dunsworth, 10.1038/s41534-017-0017-3.
K. Arya, et al. A blueprint for demonstrating quan- [32] G. Verdon, M. Broughton, and J. Biamonte. A quantum al-
tum supremacy with superconducting qubits. Science, 360 gorithm to train neural networks using low-depth circuits.
(6385):195–199, 2018. doi:10.1126/science.aao4309. arXiv preprint arXiv:1712.05304, 2017.
[17] T. Douce, D. Markham, E. Kashefi, E. Diamanti, T. Coudreau, [33] J. Otterbach, R. Manenti, N. Alidoust, A. Bestwick, M. Block,
P. Milman, P. van Loock, and G. Ferrini. Continuous- B. Bloom, S. Caldwell, N. Didier, E. S. Fried, S. Hong, et al.
variable instantaneous quantum computing is hard to Unsupervised machine learning on a hybrid quantum com-
sample. Physical Review Letters, 118(7), 2017. doi: puter. arXiv preprint arXiv:1712.05771, 2017.
10.1103/PhysRevLett.118.070503. [34] M. Schuld and N. Killoran. Quantum machine learning in
[18] A. Finnila, M. Gomez, C. Sebenik, C. Stenson, and J. Doll. feature hilbert spaces. Physical Review Letters, 122:040504,
Quantum annealing: a new method for minimizing multi- Feb 2019. doi:10.1103/PhysRevLett.122.040504.
dimensional functions. Chemical Physics Letters, 219(5-6): [35] A. S. Green, P. L. Lumsdaine, N. J. Ross, P. Selinger, and B. Val-
343–348, 1994. doi:10.1016/0009-2614(94)00117-0. iron. Quipper: a scalable quantum programming language.
[19] M. W. Johnson, M. H. Amin, S. Gildert, T. Lanting, F. Hamze, In ACM SIGPLAN Notices, volume 48, pages 333–342. ACM,
N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, 2013. doi:10.1145/2491956.2462177.
et al. Quantum annealing with manufactured spins. Nature, [36] D. Wecker and K. M. Svore. Liqui|>: A software design ar-
473(7346):194–198, 2011. doi:10.1038/nature10012. chitecture and domain-specific language for quantum com-
[20] M.-H. Yung, J. Casanova, A. Mezzacapo, J. McClean, puting. arXiv preprint arXiv:1402.4467, 2014.
L. Lamata, A. Aspuru-Guzik, and E. Solano. From transistor [37] A. JavadiAbhari, S. Patil, D. Kudrow, J. Heckey, A. Lvov, F. T.
to trapped-ion computers for quantum chemistry. Scientific Chong, and M. Martonosi. ScaffCC: Scalable compilation
Reports, 4:3589, 2014. doi:10.1038/srep03589. and analysis of quantum programs. Parallel Computing, 45:
[21] P. O’Malley, R. Babbush, I. Kivlichan, J. Romero, J. Mc- 2–17, 2015. doi:10.1016/j.parco.2014.12.001.
Clean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, [38] M. Smelyanskiy, N. P. Sawaya, and A. Aspuru-Guzik. qHiP-
et al. Scalable quantum simulation of molecular en- STER: the quantum high performance software testing envi-
ergies. Physical Review X, 6(3):031007, 2016. doi: ronment. arXiv preprint arXiv:1601.07195, 2016.
10.1103/PhysRevX.6.031007. [39] S. Pakin. A quantum macro assembler. In High Performance
[22] Y. Shen, X. Zhang, S. Zhang, J.-N. Zhang, M.-H. Yung, Extreme Computing Conference (HPEC), 2016 IEEE, pages 1–
and K. Kim. Quantum implementation of the unitary 8. IEEE, 2016. doi:10.1109/HPEC.2016.7761637.
coupled cluster for simulating molecular electronic struc- [40] R. S. Smith, M. J. Curtis, and W. J. Zeng. A practi-
ture. Physical Review A, 95(2):020501, 2017. doi: cal quantum instruction set architecture. arXiv preprint
10.1103/PhysRevA.95.020501. arXiv:1608.03355, 2016.
[23] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, [41] D. S. Steiger, T. Häner, and M. Troyer. ProjectQ: an open
J. M. Chow, and J. M. Gambetta. Hardware-efficient varia- source software framework for quantum computing. Quan-
tional quantum eigensolver for small molecules and quantum, 2:49, 2018. doi:10.22331/q-2018-01-31-49.
tum magnets. Nature, 549(7671):242–246, 2017. doi: [42] A. W. Cross, L. S. Bishop, J. A. Smolin, and J. M. Gam-
10.1038/nature23879. betta. Open quantum assembly language. arXiv preprint
[24] G. Rosenberg, P. Haghnegahdar, P. Goddard, P. Carr, K. Wu, arXiv:1707.03429, 2017.
and M. L. de Prado. Solving the optimal trading trajectory [43] E. S. Fried, N. P. Sawaya, Y. Cao, I. D. Kivlichan, J. Romero,
problem using a quantum annealer. IEEE Journal of Selected and A. Aspuru-Guzik. qtorch: The quantum tensor con-
Topics in Signal Processing, 10(6):1053–1060, 2016. doi: traction handler. PloS one, 13(12):e0208510, 2018. doi:
10.1109/JSTSP.2016.2574703. 10.1371/journal.pone.0208510.
[25] A. Lucas. Ising formulations of many NP problems. Frontiers [44] A. J. McCaskey, E. F. Dumitrescu, D. Liakh, M. Chen, W.-c.
26
Feng, and T. S. Humble. Extreme-scale programming model trix for multimode systems: U(n) invariance, squeezing, and
for quantum acceleration within high performance comput- normal forms. Physical Review A, 49(3):1567, 1994. doi:
ing. arXiv preprint arXiv:1710.01794, 2017. 10.1103/PhysRevA.49.1567.
[45] J. R. McClean, I. D. Kivlichan, D. S. Steiger, Y. Cao, E. S. [63] W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S.
Fried, C. Gidney, T. Häner, V. Havlíček, Z. Jiang, M. Neeley, Kolthammer, and I. A. Walsmley. Optimal design for uni-
et al. OpenFermion: The electronic structure package for versal multiport interferometers. Optica, 3(12):1460–1465,
quantum computers. arXiv preprint arXiv:1710.07629, 2017. 2016. doi:10.1364/OPTICA.3.001460.
[46] S. Liu, X. Wang, L. Zhou, J. Guan, Y. Li, Y. He, R. Duan, [64] J. M. Arrazola, T. R. Bromley, J. Izaac, C. R. Myers, K. Bradler,
and M. Ying. Q|S I〉: A quantum programming environment. and N. Killoran. Machine learning method for state prepa-
In Symposium on Real-Time and Hybrid Systems, pages 133– ration and gate synthesis on photonic quantum computers.
164. Springer, 2018. doi:10.1007/978-3-030-01461-2_8. Quantum Science and Technology, 2018. doi:10.1088/2058-
[47] K. Svore, A. Geller, M. Troyer, J. Azariah, C. Granade, 9565/aaf59e.
B. Heim, V. Kliuchnikov, M. Mykhailova, A. Paz, and M. Roet- [65] N. Quesada and A. M. Brańczyk. Gaussian functions
teler. Q#: Enabling scalable quantum computing and de- are optimal for waveguided nonlinear-quantum-optical pro-
velopment with a high-level DSL. In Proceedings of the Real cesses. Physical Review A, 98:043813, 2018. doi:
World Domain Specific Languages Workshop 2018, page 7. 10.1103/PhysRevA.98.043813.
ACM, 2018. doi:10.1145/3183895.3183901. [66] N. Killoran, T. R. Bromley, J. M. Arrazola, M. Schuld, N. Que-
[48] S. Lloyd and S. L. Braunstein. Quantum computation over sada, and S. Lloyd. Continuous-variable quantum neural net-
continuous variables. Physical Review Letters, 82(8):1784, works. arXiv preprint arXiv:1806.06871, 2018.
1999. doi:10.1103/PhysRevLett.82.1784. [67] M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani.
[49] S. L. Braunstein and P. van Loock. Quantum information with Experimental realization of any discrete unitary opera-
continuous variables. Reviews of Modern Physics, 77:513– tor. Physical Review Letters, 73:58–61, 1994. doi:
577, 2005. doi:10.1103/RevModPhys.77.513. 10.1103/PhysRevLett.73.58.
[50] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. [68] A. Serafini. Quantum Continuous Variables: A Primer of The-
Ralph, J. H. Shapiro, and S. Lloyd. Gaussian quantum infor- oretical Methods. CRC Press, 2017.
mation. Reviews of Modern Physics, 84(2):621, 2012. doi: [69] S. D. Bartlett, B. C. Sanders, S. L. Braunstein, and K. Nemoto.
10.1103/RevModPhys.84.621. Efficient classical simulation of continuous variable quantum
[51] M. A. Nielsen and I. Chuang. Quantum computation and information processes. Physical Review Letters, 88:097904,
quantum information. Cambridge University Press, 2002. 2002. doi:10.1103/PhysRevLett.88.097904.
[52] A. Graves, G. Wayne, and I. Danihelka. Neural Turing ma- [70] A. Lund, A. Laing, S. Rahimi-Keshari, T. Rudolph, J. L.
chines. arXiv preprint arXiv:1410.5401, 2014. O’Brien, and T. Ralph. Boson sampling from a Gaussian
[53] A. Graves, G. Wayne, M. Reynolds, T. Harley, I. Danihelka, state. Physical Review Letters, 113(10):100502, 2014. doi:
A. Grabska-Barwińska, S. G. Colmenarejo, E. Grefenstette, 10.1103/PhysRevLett.113.100502.
T. Ramalho, J. Agapiou, et al. Hybrid computing using a [71] G. Cariolaro and G. Pierobon. Bloch-Messiah reduction of
neural network with dynamic external memory. Nature, 538 Gaussian unitaries by Takagi factorization. Physical Review
(7626):471–476, 2016. doi:10.1038/nature20101. A, 94(6):062109, 2016. doi:10.1103/PhysRevA.94.062109.
[54] S. van der Walt, S. C. Colbert, and G. Varoquaux. The NumPy [72] G. Cariolaro and G. Pierobon. Reexamination of Bloch-
array: a structure for efficient numerical computation. Com- Messiah reduction. Physical Review A, 93(6):062115, 2016.
puting in Science & Engineering, 13(2):22–30, 2011. doi: doi:10.1103/PhysRevA.93.062115.
10.1109/MCSE.2011.37. [73] C. H. Bennett, G. Brassard, C. Crépeau, R. Jozsa, A. Peres,
[55] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, and W. K. Wootters. Teleporting an unknown quantum
C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, et al. state via dual classical and Einstein-Podolsky-Rosen chan-
Tensorflow: Large-scale machine learning on heterogeneous nels. Physical Review Letters, 70:1895–1899, 1993. doi:
distributed systems. arXiv preprint arXiv:1603.04467, 2016. 10.1103/PhysRevLett.70.1895.
[56] N. Quesada. A faster calculation of franck-condon factors [74] A. Furusawa and P. van Loock. Quantum Teleportation and
and fock matrix elements of gaussian unitaries using loop Entanglement: A Hybrid Approach to Optical Quantum Infor-
hafnians. arXiv preprint arXiv:1811.09597, 2018. mation Processing. Wiley, 2011.
[57] D. Gottesman, A. Kitaev, and J. Preskill. Encoding a qubit in [75] W. Steeb and Y. Hardy. Problems and Solutions in Quan-
an oscillator. Physical Review A, 64(1):012310, 2001. doi: tum Computing and Quantum Information. World Scientific,
10.1103/PhysRevA.64.012310. 2006.
[58] A. Ferraro, S. Olivares, and M. G. Paris. Gaussian states [76] M. Gu, C. Weedbrook, N. C. Menicucci, T. C. Ralph, and
in continuous variable quantum information. arXiv preprint P. van Loock. Quantum computing with continuous-variable
quant-ph/0503237, 2005. clusters. Physical Review A, 79:062318, 2009. doi:
[59] J. Williamson. On the algebraic problem concerning the nor- 10.1103/PhysRevA.79.062318.
mal forms of linear dynamical systems. American Journal of [77] D. Gottesman and I. L. Chuang. Demonstrating the viabil-
Mathematics, 58(1):141, 1936. doi:10.2307/2371062. ity of universal quantum computation using teleportation
[60] C. Bloch and A. Messiah. The canonical form of an anti- and single-qubit operations. Nature (London), 402:390–393,
symmetric tensor and its application to the theory of su- 1999. doi:10.1038/46503.
perconductivity. Nuclear Physics, 39:95–106, 1962. doi: [78] S. D. Bartlett and W. J. Munro. Quantum teleportation of
10.1016/0029-5582(62)90377-2. optical quantum gates. Physical Review Letters, 90:117901,
[61] S. L. Braunstein. Squeezing as an irreducible re- 2003. doi:10.1103/PhysRevLett.90.117901.
source. Physical Review A, 71(5):055801, 2005. doi: [79] M. Tillmann, B. Dakić, R. Heilmann, S. Nolte, A. Szameit,
10.1103/PhysRevA.71.055801. and P. Walther. Experimental boson sampling. Nature Photon-
[62] R. Simon, N. Mukunda, and B. Dutta. Quantum-noise ma- ics, 7(7):540–544, 2013. doi:10.1038/nphoton.2013.102.
27
[80] N. Spagnolo, C. Vitelli, M. Bentivegna, D. J. Brod, A. Crespi, [97] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders. Effi-
F. Flamini, S. Giacomini, G. Milani, R. Ramponi, P. Mataloni, cient quantum algorithms for simulating sparse Hamiltoni-
R. Osellame, E. F. Galvão, and F. Sciarrino. Experimental ans. Communications in Mathematical Physics, 270(2):359–
validation of photonic boson sampling. Nature Photonics, 8 371, 2006. doi:10.1007/s00220-006-0150-x.
(8):615–620, 2014. doi:10.1038/nphoton.2014.135. [98] T. Kalajdzievski, C. Weedbrook, and P. Rebentrost.
[81] A. Crespi, R. Osellame, R. Ramponi, D. J. Brod, E. F. Galvão, Continuous-variable gate decomposition for the bose-
N. Spagnolo, C. Vitelli, E. Maiorino, P. Mataloni, and F. Scia- hubbard model. Physical Review A, 97:062311, Jun 2018.
rrino. Integrated multimode interferometers with arbitrary doi:10.1103/PhysRevA.97.062311.
designs for photonic boson sampling. Nature Photonics, 7(7): [99] T. Sowiński, O. Dutta, P. Hauke, L. Tagliacozzo, and
545–549, 2013. doi:10.1038/nphoton.2013.112. M. Lewenstein. Dipolar molecules in optical lat-
[82] J. B. Spring, B. J. Metcalf, P. C. Humphreys, W. S. Koltham- tices. Physical Review Letters, 108:115301, 2012. doi:
mer, X.-M. Jin, M. Barbieri, A. Datta, N. Thomas-Peter, N. K. 10.1103/PhysRevLett.108.115301.
Langford, D. Kundys, J. C. Gates, B. J. Smith, P. G. R.
Smith, and I. A. Walmsley. Boson sampling on a pho-
tonic chip. Science, 339(6121):798–801, 2012. doi:
10.1126/science.1231692.
[83] L. Valiant. The complexity of computing the permanent.
Theoretical Computer Science, 8(2):189–201, 1979. doi:
10.1016/0304-3975(79)90044-6.
[84] E. R. Caianiello. Combinatorics and renormalization in quan-
tum field theory, volume 38. Benjamin, 1973.
[85] A. Barvinok. Combinatorics and complexity of partition func-
tions, volume 274. Springer, 2016.
[86] K. Brádler, P.-L. Dallaire-Demers, P. Rebentrost, D. Su, and
C. Weedbrook. Gaussian boson sampling for perfect match-
ings of arbitrary graphs. Physical Review A, 98:032310, Sep
2018. doi:10.1103/PhysRevA.98.032310.
[87] A. Björklund, B. Gupt, and N. Quesada. A faster hafnian
formula for complex matrices and its benchmarking on the
titan supercomputer. arXiv preprint arXiv:1805.12498, 2018.
[88] N. Quesada, J. M. Arrazola, and N. Killoran. Gaussian boson
sampling using threshold detectors. Physical Review A, 98:
062322, 2018. doi:10.1103/PhysRevA.98.062322.
[89] B. Gupt, J. M. Arrazola, N. Quesada, and T. R. Bromley. Clas-
sical benchmarking of gaussian boson sampling on the titan
supercomputer. arXiv preprint arXiv:1810.00900, 2018.
[90] D. Shepherd and M. J. Bremner. Temporally unstruc-
tured quantum computation. Proceedings of the Royal
Society of London A: Mathematical, Physical and En-
gineering Sciences, 465(2105):1413–1439, 2009. doi:
10.1098/rspa.2008.0443.
[91] M. J. Bremner, A. Montanaro, and D. J. Shepherd. Achieving
quantum supremacy with sparse and noisy commuting quan-
tum computations. Quantum, 1:8, 2017. doi:10.22331/q-
2017-04-25-8.
[92] A. P. Lund, M. J. Bremner, and T. C. Ralph. Quantum sam-
pling problems, BosonSampling and quantum supremacy.
npj Quantum Information, 3(1), 2017. doi:10.1038/s41534-
017-0018-2.
[93] J. M. Arrazola, P. Rebentrost, and C. Weedbrook. Quantum
supremacy and high-dimensional integration. arXiv preprint
arXiv:1712.07288, 2017.
[94] A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-
Gordon. Simulated quantum computation of molecular
energies. Science, 309(5741):1704–1707, 2005. doi:
10.1126/science.1113479.
[95] J. D. Whitfield, J. Biamonte, and A. Aspuru-Guzik. Simu-
lation of electronic structure Hamiltonians using quantum
computers. Molecular Physics, 109(5):735–750, 2011. doi:
10.1080/00268976.2011.552441.
[96] A. M. Childs and N. Wiebe. Hamiltonian simulation using
linear combinations of unitary operations. Quantum Infor-
mation and Computation, 12(11-12):901–924, 2012. doi:
10.26421/QIC12.11-12.

Strawberry Fiels

Uploaded by

Copyright:

Available Formats

Strawberry Fiels

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Strawberry Fiels

Uploaded by

Copyright:

Available Formats

Strawberry Fields:

A Software Platform for Photonic Quantum Computing

ing theme for these software efforts is to target qubit-

berry Fields also includes a suite of three CV quantum simu-

2. We provide a suite of important circuit decomposi- CV Qubit

where [x̂, p̂] = iħ

Gaussian states Coherent states |α〉 α∈C z=0

the form universality. A number of other useful CV gates are listed in

which is parameterized via the mean photon number n :=

5 alpha = 2.0 + 1j 19 # Apply Hermitian conjugate of a gate

1 # Homodyne measurement at angle phi 24 eng.run(backend="gaussian")

Running Blackbird Programs in Python 1 # A Blackbird circuit where gates

3 eng, q = sf.Engine(num_subsystems=2) 12 # Attempting to act on registers which have

8 # Use a measured register symbolically Decompositions

5 # A Blackbird computation using classical 1 U = np.array([[1-1j, np.sqrt(2)],

8 MeasureX | q[0] 4 eng, q = sf.Engine(2)

10 eng.run("gaussian") 6 with eng:

When we do this, any registers measured in the circuit will

provides a number of advanced functionalities, including 32 eng, q = sf.Engine(num_subsystems=2)

("tf") backend. To keep the underlying simulation fully 42 sess = tf.Session()

By taking advantage of these additional functionalities of

56 eng, q = sf.Engine(num_subsystems=2) Strawberry Fields’ Quantum Simulators

Fock Backend TensorFlow Backend

cussion of these transformations, see Sec. 2 of [50].

followed by another multiport V̂ , followed by single-mode

Note that if the Gaussian state is pure (which happens

Appendix B: Strawberry Fields Operations

Operation Name Definition

Coherent(a) Coherent state A displaced vacuum state, |α〉 = D(α) |0〉, α ∈ C

A squeezed vacuum state, |z〉 = S(z) |0〉,

A squeezed then displaced vacuum state,

p1 (|α〉 + exp(iπp) |−α〉), where p = 0, 1 gives an

Prepare an arbitrary multi-mode pure state, repre-

Prepare an arbitrary multi-mode mixed state, repre-

Operation Name Definition

Operation Name Definition

BSgate(theta,phi) Beamsplitter where the transmissivity and reflectivity

Table VII: Two-mode gate operations available in Strawberry Fields.

Operation Name Definition

Table VIII: Multi-mode Gaussian decompositions available in Strawberry Fields.

Operation Name Definition

Quadratic phase gate S2 (z) = B † (π/4, 0) [S(z) ⊗ S(−z)] B(π/4, 0)

The quadratic phase shift gate Pgate(s) is decomposed

P(s) = R(θ )S(r e iφ ),

Controlled addition gate

C X (s) = B(φ, 0) [S(r) ⊗ S(−r)] B(π/2 + φ, 0),

where sin(2φ) = −1/cosh(r), cos(2φ) = − tanh(r), and

Controlled phase gate

17 # maximally entangled states

quantum computation [74]. 34 state = eng.run("gaussian")

surements, all ‘stacked vertically’ (i.e., coupled to each con-

6 with eng: instead an example of an intermediate quantum computer,

q0 : |0〉 S R n̂ scalability, due to its dependence on an array of simultane-

5 # initialise the engine and register

Instantaneous Quantum Polynomial (IQP) 7 with eng:

15 # gates diagonal in the x quadrature

the Lie product formula, we find that

where θ = −J t/k, φ = π/2, and r = −U t/2k. Using J = 1,