Strawberry Fields: A Software Platform for Photonic Quantum Computing

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 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 Blackbird 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.


Introduction
The decades-long worldwide quest to build practical quantum computers is currently undergoing a critical period. During the next few years, a number of different quantum devices will become available to the public. While fault-tolerant quantum computers will one day provide significant computational speedups for problems like factoring [1], search [2], or linear algebra [3], near-term quantum devices will be noisy, approximate, and sub-universal [4]. Nevertheless, these emerging quantum processors are expected to be strong enough to show a computational advantage over classical computers for certain problems, an achievement known as quantum computational supremacy.
To this end, a nascent quantum software ecosystem has recently begun to develop [35][36][37][38][39][40][41][42][43][44][45][46][47]. However, a prevail-ing theme for these software efforts is to target qubitbased quantum devices. In reality, there are several competing models of quantum computing which are equivalent in computational terms, but which are conceptually quite distinct. One prominent approach is the continuous variable (CV) model of quantum computing [48][49][50]. In the CV model, the basic information-processing unit is an infinite-dimensional bosonic mode, making it particularly well-suited for implementations and applications based on light. The CV model retains the computational power of the qubit model (cf. Chap. 4 of Ref. [51]), while offering a number of unique features. For instance, the CV model is a natural fit for simulating bosonic systems (electromagnetic fields, trapped atoms, harmonic oscillators, Bose-Einstein condensates, phonons, or optomechanical resonators) and for settings where continuous quantum operators -such as position & momentum -are present. Even in classical computing, recent advances from deep learning have demonstrated the power and flexibility of a continuous representation of computation [52,53] in comparison to the discrete computational model which has historically dominated.
Here we introduce Strawberry Fields 1 , an open-source software architecture for photonic quantum computing. Strawberry Fields is a full-stack quantum software platform, implemented in Python, specifically targeted to the CV model. Its main element is a new quantum programming language named Blackbird. To lay the groundwork for future photonic quantum computing hardware, Straw-berry Fields also includes a suite of three CV quantum simulator 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 simulator backends or, when they are available, on photonic quantum computers. To accompany the library, an online service for interactive exploration and simulation of CV circuits is available at strawberryfields.ai.
Aside from being the first quantum software framework to support photonic quantum computation with continuous-variables, Strawberry Fields provides additional computational features not presently available in the quantum software ecosystem: 1. We provide two numeric simulators; a Gaussian backend, and a Fock-basis backend. These two formulations are unique to the CV model of quantum computation due to the use of an infinite Hilbert space, and came with their own technical challenges.
(a) The Gaussian backend provides state-of-the-art methods and functions for calculating the fidelity and Fock state probabilities, involving calculations of the classically intractable hafnian [56].
(b) The Fock backend allows operations such as squeezing and beamsplitters to be performed in the Fock-basis, a computationally intensive calculation that has been highly vectorized and benchmarked for performance.
2. We provide a suite of important circuit decompositions appearing in quantum photonics -such as the Williamson, Bloch-Messiah, and Clements decompositions.
3. The Fock-basis backend written using the TensorFlow machine learning library allows for symbolic calculations, automatic differentiation, and backpropagation through CV quantum simulations. As far as we are aware, this is the first quantum simulation library written using a high-level machine learning library, with support for dataflow programming and automatic differentiation.
The remainder of this white paper is structured as follows. Before presenting Strawberry Fields, we first provide a brief overview of the key ingredients for CV quantum computation, specifically the most important states, gates, and measurements. We then introduce the Strawberry Fields architecture in full, presenting the Blackbird quantum assembly 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 internal representations that they employ. In the Appendices, we give further mathematical and software details and provide full example code for a number of important CV quantum computing tasks.

Quantum Computation with Continuous Variables
Many physical systems in nature are intrinsically continuous, with light being the prototypical example. Such systems reside in an infinite-dimensional Hilbert space, offering a paradigm for quantum computation which is distinct from the discrete qubit model. This continuous-variable model takes its name from the fact that the quantum operators underlying the model have continuous spectra. It is possible to embed qubit-based computations into the CV picture [57], so the CV model is as powerful as its qubit counterparts.

From Qubits to Qumodes
A high-level comparison of CV quantum computation with the qubit model is depicted in Table I. In the remainder of this section, we will provide a basic presentation of the key elements of the CV model. A more detailed technical overview can be found in Appendix A. Readers experienced with CV quantum computing can safely skip to the next section. The most elementary CV system is the bosonic harmonic oscillator, defined via the canonical mode operatorsâ and a † . These satisfy the well-known commutation relation [â,â † ] = . It is also common to work with the quadrature operators (also called the position & momentum oper-ators) 2 where [x,p] = iħ hˆ . We can picture a fixed harmonic oscillator 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 qumodeseach with corresponding operatorsâ i andâ † i -and interacting them via sequences of suitable quantum gates, we can implement a general CV quantum computation.

CV States
The dichotomy between qubit and CV systems is perhaps most evident in the basis expansions of quantum states: Qumode For qubits, we use a discrete set of coefficients; for CV systems, we can have a continuum. The states |x〉 are the eigenstates of thex quadrature,x |x〉 = x |x〉, with x ∈ . These quadrature states are special cases of a more general family of CV states, the Gaussian states, which we now introduce.

Gaussian states
Our starting point is the vacuum state |0〉. Other states can be created by evolving the vacuum state according to where H is a bosonic Hamiltonian (i.e., a function of the operatorsâ i andâ † i ) and t is the evolution time. States where the Hamiltonian H is at most quadratic in the operatorsâ i andâ † i (equivalently, inx i andp i ) are called Gaussian. For a single qumode, Gaussian states are parameterized by two continuous complex variables: a displacement parameter α ∈ and a squeezing parameter z ∈ (often expressed as z = r exp(iφ), with r ≥ 0). Gaussian states are so-named because we can identify each Gaussian state with a corresponding Gaussian distribution. For single qumodes, the identification proceeds through its displacement and squeezing parameters. The displacement gives 2 It is common to picture ħ h as a (dimensionless) scaling parameter for thex andp operators rather than a physical constant. However, there are several conventions for the scaling value in common use [58]. These self-adjoint operators are proportional to the Hermitian and anti-Hermitian parts of the operatorâ. Strawberry Fields allows the user to specify this value, with the default ħ h = 2.

FIG. 1:
Schematic representation of a Gaussian state for a single mode. The shape and orientation are parameterized by the displacement α and squeezing z = r exp(iφ).
the centre of the distribution, while the squeezing determines the variance and rotation of the distribution (see Fig.  1). Multimode Gaussian states, on the other hand, are parameterized by a vector of displacementsr and a covariance matrix V. Many important pure states in the CV model are special cases of the pure Gaussian states; see Table II for a summary.

State family Displacement Squeezing
Vacuum state |0〉 while (undisplaced) squeezed states only have even number states in their expansion:

Mixed states
Mixed Gaussian states are also important in the CV picture, for instance, the thermal state which is parameterized via the mean photon number n := Tr(ρ(n)n). Starting from this state, we can consider a mixed-state-creation process similar to Eq. (5), namely Analogously to pure states, by applying Hamiltonians of second-order (or lower) to thermal states, we generate the family of Gaussian mixed states.

CV Gates
Unitary operations can be associated with a generating Hamiltonian H via the recipe (cf. Eqs.
For convenience, we classify unitaries by the degree of their generator. A CV quantum computer is said to be universal if it can implement, to arbitrary precision and with a finite number of steps, any unitary which is polynomial in the mode operators [48]. We can build a multimode unitary by applying a sequence of gates from a universal gate set, each of which acts only on one or two modes. We focus on a universal set made from the following two subsets: Gaussian gates: Single and two-mode gates which are at most quadratic in the mode operators, e.g., Displacement, Rotation, Squeezing, and Beamsplitter gates.

Non-Gaussian gate:
A single-mode gate which is degree 3 or higher, e.g., the Cubic phase gate.
A number of fundamental CV gates are presented in Table III. Many of the Gaussian states from the previous section are connected to a corresponding Gaussian gate. Any multimode Gaussian gate can be implemented through a suitable combination of Displacement, Rotation, Squeezing, and Beamsplitter Gates [50], making these gates sufficient for constructing all quadratic unitaries. The cubic phase gate is presented as an exemplary non-Gaussian gate, but any other non-Gaussian gate could also be used to achieve universality. A number of other useful CV gates are listed in Appendix B.

Gate
Unitary Symbol

CV Measurements
As with CV states and gates, we can distinguish between Gaussian and non-Gaussian measurements. The Gaussian class consists of two (continuous) types: homodyne and heterodyne measurements, while the key non-Gaussian measurement is photon counting. These are summarized in Table IV.

Homodyne measurements
Ideal homodyne detection is a projective measurement onto the eigenstates of the quadrature operatorx. These states form a continuum, so homodyne measurements are inherently continuous, returning values x ∈ . More generally, we can consider projective measurement onto the eigenstates x φ of the Hermitian operator This is equivalent to rotating the state clockwise by φ and performing anx-homodyne measurement. If we have a multimode Gaussian state and we perform homodyne measurement on one of the modes, the conditional state of the unmeasured modes remains Gaussian.

Heterodyne measurements
Whereas homodyne detection is a measurement ofx, heterodyne detection can be seen as a simultaneous measurement of bothx andp. Because these operators do not commute, they cannot be simultaneously measured without some degree of uncertainty. Equivalently, we can picture heterodyne measurement as projection onto the coherent states, with measurement operators 1 π |α〉〈α|. Because the coherent states are not orthogonal, there is a corresponding lack of sharpness in the measurements. If we perform heterodyne measurement on one mode of a multimode state, the conditional state on the remaining modes stays Gaussian.

Photon Counting
Photon counting (also known as as photon-number resolving measurement), is a complementary measurement method to the '-dyne' measurements, revealing the particlelike, rather than the wave-like, nature of qumodes. This measurement projects onto the number eigenstates |n〉, returning non-negative integer values n ∈ . Except for the outcome n = 0, a photon-counting measurement on a single mode of a multimode Gaussian state will cause the remaining modes to become non-Gaussian. Thus, photoncounting can be used as an ingredient for implementing non-Gaussian gates. A related process is photodetection, where a detector only resolves the vacuum state from nonvacuum states. This process has only two measurement operators, namely |0〉〈0| and − |0〉〈0|.

The Strawberry Fields Software Platform
The Strawberry Fields library has been designed with several key goals in mind. Foremost, it is a standard-bearer for the CV model, laying the groundwork for future photonic quantum computers. As well, Strawberry Fields is designed to be simple to use, giving entry points for as many users as possible. Finally, since the potential applications of nearterm quantum computers are still being worked out, it is important that Strawberry Fields provides powerful tools to easily explore many different use-cases and applications.
Strawberry Fields has been implemented in Python, a modern language with a gentle learning curve which is already familiar to many programmers and scientific practitioners. The accompanying quantum simulator backends are built upon the widely used Python packages NumPy and TensorFlow. All Strawberry Fields code is open source. Strawberry Fields can be accessed programmatically as a Python package, or via a browser-based interface for designing quantum circuits.
A pictorial outline of Strawberry Fields' key elements and their interdependencies is presented in Fig. 2. Conceptually, the software stack is separated into two main pieces: a user-facing frontend layer and a lower-level backends component. The frontend encompasses the Strawberry Fields Python API and the Blackbird quantum assembly language.
These elements provide access points for users to design quantum circuits. These circuits are then linked to a backend via a quantum compiler engine. For a backend, the engine can currently target one of three included quantum computer simulators. When CV quantum processors become available in the near future, the engine will also build and run circuits on those devices. Further, high-level quantum computing applications can be built by leverging the Strawberry Fields frontend API. Existing examples include the Strawberry Fields Interactive website, the Quantum Machine Learning Toolbox (for streamlining the training of variational quantum circuits), and SFOpenBoson (an interface between the electronic structure library OpenFermion [45] and Strawberry Fields).
In the remainder of this section, the key elements of Strawberry Fields will be presented in more detail. Proceeding through a series of examples, we show how CV quantum computations can be defined using the Blackbird language, then compiled and run on a quantum computer backend. We also outline how to use Strawberry Fields for optimization and machine learning on quantum circuits. Finally, we discuss the suite of quantum computer simulators included within Strawberry Fields.

Blackbird: A Quantum Programming Language
As classical computers have become progressively more powerful, the languages used to program them have also undergone considerable paradigmatic changes. Machine code gave way to human-readable assembly languages, followed by higher-level procedural and object-oriented languages. With each generation, the trend has been towards higher levels of abstraction, separating the programmer more and more from details of the actual computer hardware. Quantum computers are still at an early stage of development, so while we can imagine what higher-level quantum programming might look like, in the near term we first need to build languages which are conceptually closer to the quantum hardware.
Blackbird is a standalone domain specific language (DSL) for continuous-variable quantum computation. With a welldefined grammar in extended Backus-Naur form, and both Python and C++ parsers available, Blackbird provides operations that match the basic CV states, gates, and measurements, and maps directly to low-level hardware instructions. The abstract syntax keeps a close connection between code and the quantum operations that they implement; this syntax is modeled after that of ProjectQ [41], but specialized to the CV setting. Blackbird can be used as part of the Strawberry Fields stack, but also directly with photonic quantum computing hardware systems.
Within the Strawberry Fields framework, we have built an implementation of Blackbird using Python 3 as the embedding language -an 'embedded' DSL. This 'Pythonenhanced' Blackbird language provides the same core operations and follows the same grammar and syntactical rules as the standalone DSL, but, by nature, may also contain valid Python constructs. Furthermore, Strawberry Fields' 'Python-enhanced' Blackbird provides users with additional quantum operations that are decomposed into lower-level Blackbird assembly commands. We will introduce the elements of Blackbird through a series of basic examples, discussing more technical aspects as they arise.

Operations
Quantum computations consist of four main ingredients: state preparations, application of gates, performing measurements, and adding/removing subsystems. In Blackbird, these are all considered as Operations, and share the same basic syntax. In the following code examples, we use the variable q for a set of qumodes (more specifically, a quantum register), the details of which are deferred until the next section.
Our first considered Operation is state preparation. By default, qumodes are initialized in the vacuum state. Various other important CV states can be created with simple Blackbird commands. 1  Blackbird state preparations such as those used in Codeblock 1 implicitly reset the existing state of the qumodes.
Conceptually, the vertical bar symbol '|' separates Operations -like state preparation -from the registers that they act upon. Notice that we can use Operations inline, or construct them separately and reuse them several times.
After creating states, we will want to transform these using quantum gates. 1  Blackbird supports all of the gates listed in the previous section as well as a number of composite gates, each of which can be decomposed using the universal gates. The supported composite gates are: controlled X (CXgate), controlled Z (CZgate), quadratic phase (Pgate), and twomode squeezing (S2gate). A full list of gates currently supported in Blackbird can be found in Appendix B.
Finally, we can specify measurement Operations using Blackbird. Measurements have several effects. For one, the numerical result of the measurement is placed in a classical register. As well, the state of all remaining qumodes is projected to the (normalized) conditional state for that measurement value. Finally, the state of the measured qumode is reset to the vacuum state. This is the typical behaviour of photonic hardware, where measurements absorb all the energy of the measured qumode.

Running Blackbird Programs in Python
Within Python, Blackbird programs are managed by an Engine. The function Engine in the Strawberry Fields API will instantiate an Engine, returning both the Engine and its corresponding quantum register. The Engine is used as a Python context manager, providing a convenient way to encapsulate Blackbird programs. 1 # Create Engine and quantum register 2 import strawberryfields as sf 3 eng, q = sf.Engine(num_subsystems=2) The above code example is runnable and carries out a complete quantum computation, namely the preparation and measurement of an EPR entangled state. Note that Operations can be declared outside of the Engine, but their action on the quantum registers must come within the Engine context. Also notice that our register has a length of 2, so any single-mode Operations must act on specific elements, i.e., q[i], while two-mode Operations can act on q directly. Finally, the user must specify a backendas well as any backend-dependent settings -when calling eng.run(). We will discuss the Strawberry Fields backends further in a later section.

Quantum and Classical Registers
When a Strawberry Fields Engine is constructed, the user must specify the number of qumode subsystems to begin with. This number is required for the initialization of the Engine, but may change within a computation (e.g., when temporary ancilla modes are used). Qumodes can be added/deleted by using the New and Del Operations. An Engine maintains a unique numeric indexing for the quantum registers based on the order they were added. When a subsystem is deleted from the circuit, no further gates can act on that register. As stated earlier, measurement Operations produce classical information (the measurement result) and manipulate the corresponding register. Compared to previously released quantum programming frameworks, such as Pro-jectQ, PyQuil, and Qiskit, Strawberry Fields has been designed so that the notation q[i], while principally denoting the quantum register, may also encapsulate classical information or a classical register. This behaviour is contextual, and unique to Strawberry Fields. For example, prior to measurement, q[i] simply references a quantum register. Once a measurement operation is performed, q[i] continues to represent a quantum register -now reset to the vacuum state -as well as storing the numerical value of the measurement, accessible via the attribute q[i].val. Note that this numerical value is only available if a computation has been run up to the point of measurement. We may also use a classical measurement result symbolically as a parameter in later gates without first running the computation. To do this, we simply pass the measured register (e.g., q[i]) explicitly as an argument to the required gate. As before, the Strawberry Fields quantum register object is contextualwhen passed as a gate argument, Strawberry Fields implicitly accesses the encapsulated classical register. 1 # Numerical evaluation of a measurement 2 # result using eng.run() 3 with eng: # Use a measured register symbolically 9 # in another gate 10 with eng: 13 eng.run("gaussian") Codeblock 6: Evaluating measurement results numerically and using them symbolically.
In quantum algorithms, it is common to process a measurement result classically and use the post-processed value as a parameter for further operations in a circuit. Strawberry Fields provides the convert decorator to transform a userspecified numerical function into one which acts on registers.

Post-selection
The measurement Operations in Strawberry Fields are stochastic in nature, with outcomes determined by some underlying quantum probability distribution. Often it is convenient to select specific values for these measurements rather than sampling them. For instance, we might want to explore the conditional state created by a specific value, determine the measurement-dependent corrections we need to make in a teleportation circuit, or even design an algorithm which inherently contains post-selection. This functionality is supported in Strawberry Fields through the optional keyword argument select, which can be supplied for any measurement Operation. The measurement outcome will then return exactly this value, while the remaining modes will be projected into the conditional state corresponding to this value 3 .
Codeblock 8: Selecting a specific desired measurement outcome.

Decompositions
In addition to the core CV operations discussed above, Strawberry Fields also provides support for some important decompositions frequently used in quantum optics. These include the (a) Williamson decomposition [59], for decomposing arbitrary Gaussian states to a symplectic transformation acting on a thermals state, (b) the Bloch-Messiah decomposition [60][61][62], for decomposing the action of symplectic transformations to interferometers and single-mode squeezing, and (c) the Clements decomposition [63], for decomposing multi-mode linear interferometers into arrays of beamsplitters and rotations of fixed depth. In all cases, the resulting decomposition into the 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 in quantum optical circuits, and viewing the resulting physical components.

Optimization and Quantum Machine Learning
Strawberry Fields can perform quantum circuit simulations using both numerical and symbolic representations. Numerical computation is the default operating mode and is supported by all three supplied backends. Symbolic computation is enabled only for the TensorFlow backend. In this section, we outline the main TensorFlow functionalities accessible through the Strawberry Fields frontend interface. More details about the corresponding TensorFlow backend are discussed in the next section. TensorFlow [55] models calculations abstractly using a computational graph, where individual operations are represented as nodes and their dependencies by directed edges. This viewpoint separates the symbolic representation of a computation from its numerical evaluation, and makes optimization and machine learning more amenable. On top of this, TensorFlow provides a number of advanced functionalities, including automatic gradient computation, GPU utilization, built-in optimization algorithms, and various other machine learning tools. Note that the term 'quantum machine learning' will be used here in a hybrid sense, i.e., applying conventional machine learning methods to quantum systems.
To build a TensorFlow computational graph using Strawberry Fields, we instantiate an Engine, declare a circuit in Blackbird code, then execute eng.run on the TensorFlow ("tf") backend. To keep the underlying simulation fully symbolic, the extra argument eval=False must be given. 1 # Create Engine and run symbolic computation 2 import strawberryfields as sf 3 from strawberryfields.ops import * 4 import tensorflow as tf 5 6 eng, q = sf.Engine(num_subsystems=1) 7 with eng:  In the above example, we supplied an additional feed_dict argument when evaluating. This is a Python dictionary which specifies the numerical values (typically, coming from a dataset) for every placeholder that appears in a circuit. However, as can be seen from the example, it is also possible to substitute desired values for other nodes in the computation, including the values stored in quantum registers. This allows us to easily post-select measurement values and explore the resulting conditional states.
We can also perform post-processing of measurement results when working with TensorFlow objects. In this case, the functions decorated by convert should be written using TensorFlow Tensors, Variables, and operations. 48 # Use a simple neural network as 49 # a processing function 50 @sf.convert 51  The TensorFlow backend additionally supports the use of batched processing, allowing for many evaluations of a quantum circuit to potentially be computed in parallel. Scalars are automatically broadcast to the specified batch size. Finally, we can easily run circuit simulations on special-purpose hardware like GPUs or TPUs. 68  By taking advantage of these additional functionalities of the TensorFlow backend, we can straightforwardly perform optimization and machine learning on quantum circuits in Strawberry Fields [64,65]. A complete code example for optimization of a quantum circuit is located in Appendix C.

Strawberry Fields' Quantum Simulators
The ultimate goal is for Blackbird programs to be carried out on photonic quantum computers. To lay the groundwork for these emerging devices, Strawberry Fields comes with a suite of three CV quantum computer simulators specially designed for the CV model. These simulators target different use-cases and support different functionality. For example, many important algorithms in the CV formalism involve only Gaussian states, operations, and measurements. We can take advantage of this structure to more efficiently simulate such computations. Other circuits have inherently non-Gaussian elements to them; for these, the Fock basis provides the standard description. These representations are available in Strawberry Fields in the Gaussian backend and the Fock backend, respectively. The third builtin backend is the TensorFlow backend. Also using the Fock representation, this backend is geared primarily towards optimization and machine learning applications.
Most Blackbird operations are supported across all three backends. A small subset, however, are not supported uniformly due to mathematical incompatibility. For example, the Gaussian backend does not generally support non-Gaussian gates or the preparation/measurement of Fock states. Sometimes it is also useful to work with a backend directly. To allow this, Strawberry Fields provides a backend API, giving access to additional methods and properties which are not part of the streamlined frontend API. A standalone backend can be created in Strawberry Fields using strawberryfields.backend.load_backend(name) where name is one of "gaussian", "fock", or "tf" (for comparison, the backend associated to an Engine eng is available via eng.backend).
Three important backend methods, common to all simulators, are begin_circuit, reset, and state. The first command instantiates the circuit simulation, the second resets the simulation back to an initial vacuum state, clearing all previous operations, and the third returns a class which encapsulates the current quantum state of the simulator. In addition to containing the numerical (or symbolic) state data, state classes also contain a number of useful methods and attributes for further exploring the quantum state, such as fidelity, mean_photon, or wigner. As a convenience for the user, all simulations carried out via eng.run will return a state class representing the final circuit state (see Appendix C for examples).

Gaussian Backend
This backend, written in NumPy, uses the symplectic formalism to represent CV systems. At a high level, this repre-sentation tracks the quantum state of an N -mode quantum system using two Gaussian components: a 2N -dimensional displacement vectorr and a 2N × 2N -dimensional covariance matrix V (a deeper technical overview is located in Appendix A). After we have created a Gaussian state (either via state = eng.run(backend="gaussian") or by directly calling the state method of a Gaussian backend), we can accessr and V via state.means and state.cov, respectively. Other useful Gaussian state methods are displacement and squeezing, which return the Gaussian parameters associated to the underlying state.
The scaling of the symplectic representation with the number of modes is O(N 2 ). On one hand, this is quite powerful. It allows us to to efficiently simulate any computations which are fully Gaussian. On the other, the formalism is not expressive enough to simulate more general quantum computations. Only a small number of non-Gaussian methods are available for this backend. These are auxiliary methods where we extract some non-Gaussian information from a Gaussian state, but do not update the state of the circuit. One such method is fock_prob, which is implemented using an optimized -yet still exponentially scalingalgorithm. This method enables simulation of the Gaussian boson sampling algorithm [10] using the Gaussian backend; see Appendix C for a complete code example.

Fock Backend
This backend, also written in NumPy, uses a fundamentally different description for qumodes than the Gaussian representation. As discussed in the introductory sections, the Fock representation encodes quantum computation in a countably infinite-dimensional Hilbert space. This representation is faithful, allowing a precise description of CV systems, in particular non-Gaussian circuits. Yet simulating infinite-dimensional systems leads to some computational tradeoffs which are not present for qubit simulators. Most importantly, we impose a cutoff dimension D for simulations (chosen by the user), so the Fock backend only covers a restricted set of number states {|0〉 , ..., |D − 1〉} for each mode. The size of simulated quantum systems thus depends on both the number of subsystems N and the cutoff, being O(D N ). We contrast this with qubit systems, where the base is fixed, i.e., O(2 N ). While these scalings are both exponential, in practice simulating qumode systems for D > 2 is more computationally demanding than qubits. This is because the (truncated) qumode subsystems have a higher dimension and thus encode more information than their qubit 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 energy (e.g., limited by the power of a laser). On the other hand, there are certainly states which can be easily prepared in the lab, yet would not fit accurately on the simulator. Thus, some care must be taken to trade off between the numerical cutoff value, the number of modes, and the energy scale of the circuit. If the energy scale is sufficiently low that all states fit within the specified cutoff, then simulations with the Fock and Gaussian backends will be in numerical agreement. Like

TensorFlow Backend
The other built-in backend for Strawberry Fields is coded using TensorFlow. As a simulator, it uses the same internal representation as the Fock backend (Fock basis, finite cutoff, pure vs. mixed state representations, etc.) and has the same methods. It can operate as a numerical simulator similar to the other backends. Its main purpose, however, is to leverage the many powerful tools provided by TensorFlow to enable optimization and machine learning on quantum circuits. Much of this functionality was presented in the previous section, so we will not repeat it here. Like

Conclusions
We have introduced Strawberry Fields, a multi-faceted software platform for continuous-variable quantum computing. The main components of this library -a custom quantum programming language (Blackbird), a compiler engine, and a suite of quantum simulators targeting distinct applications -have been presented in detail. Further information is available in both the Appendices and the Strawberry Fields online documentation.
The stage is now set for the broader community to use Strawberry Fields for exploration, research, and development 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 applications for quantum computing and quantum machine learning.

Appendix A: The CV model
In this Appendix, we provide more technical and mathematical details about the CV model, the quantum computing paradigm underlying Strawberry Fields.

Universal Gates for CV Quantum Computing
In discrete qubit systems, the notion of a universal gate set has the following meaning: given a set of universal gates, we can approximate an arbitrary unitary by composition of said gates. In the CV setting, we have a similar situation: we can approximate a broad set of generators -i.e., the Hamiltonians appearing in Eq. (10) -by combining elements of a CV universal gate set. However, unlike the qubit case, we do not try to approximate all conceivable unitaries. Rather, we seek to create all generators that are a polynomial function of the quadrature (or mode) operators of the system [48,49]. Remember that generators of second-degree or lower belong to the class of Gaussian operations, while all higher degrees are non-Gaussian.
We can create a higher-order generator out of lowerorder generatorsÂ andB by using the following two concatenation identities [48]: If we have two second-degree generators, such asû = x 2 +p 2 (the generator for the rotation gate) andŝ = xp +px (the generator for the squeezing gate), and a third-degree (or higher) generator, such asx 3 (the generator for the cubic phase gate), we can easily construct generators of all higher-degree polynomials, e.g.,x 4 = −[x 3 , [x 3 ,û]]/(18ħ h 2 ). This reasoning can be extended by induction to any finite-degree polynomial inx andp (equivalently, inâ andâ † ) [48,49].
In the above argument, it is important that at least one of the generators is third-degree or higher. Indeed, commutators of second-degree polynomials ofx andp are also second-degree polynomials and thus their composition using Eq. (A1) cannot generate higher-order generators. The claim can be easily extended to N -mode systems and multivariate polynomials of the operatorŝ r = (x 1 , . . . ,x N ,p 1 , . . . ,p N ). (A2) Combining single-mode universal gates (including at least one of third-degree or higher) with some multimode interaction, e.g., the beamsplitter interaction generated bŷ b i, j =p ix j −x ip j , we can construct arbitrary-degree polynomials of the quadrature operators of an N -mode system. With the above discussion in mind, we can combine the set of single-qumode gates generated by {x i ,x 2 i ,x 3 i ,û i } and the two-mode gate generated byb i, j for all pairs of modes into a universal gate set. The first, third and fourth single-mode generators correspond to the displacement, cu-bic phase and rotation gates in Table III, while the twomode generator corresponds to the beamsplitter gate. Finally, thex 2 generator corresponds to a quadratic phase gate,P(s) = exp(isx 2 /(2ħ h)). This gate can be written in terms of the single-mode squeezing gate and the rotation gate as follows:P(s) =R(θ )Ŝ(r e iφ ), where cosh(r) = 1 + (s/2) 2 , tan(θ ) = s/2, φ = −θ − sign(s)π/2. Equivalently, we could have included the squeezing generator xp +px in place of the quadratic phase and still had a universal set.
We have just outlined an efficient method to construct any gate of the form exp (−iH t), where the generator H is a polynomial of the quadrature (or mode) operators. How can this be used for quantum computations? As shown in Eq. (4), the eigenstates of thex quadrature form an (orthogonal) basis for representing qumode states. Thus, these states are often taken as the computational basis of the CV model (though other choices are also available). By applying gates as constructed above and performing measurements in this computational basis (i.e., homodyne measurements), we can carry out a CV quantum computation. One primitive example [49] is to compute the product of two numbers -whose values are stored in two qumode registers -and store the result in the state of a third qumode. Consider the generatorx 1x2p3 /ħ h, which will cause the "position" operators to evolve according tô A measurement in the computational basis of mode 3 will reveal the value of the product x 1 x 2 . Note that these encodings of classical continuous degress of freedom into quantum registers allows for the generalization of neural networks into the quantum regime [66]. In the following sections we show how more complicated algorithms are constructed.

Multiport Gate Decompositions
One important set of gates for which it is critical to derive a decomposition in terms of universal gates is the set of multiport interferometers. A multiport interferometer, represented by a unitary operatorÛ acting on N modes, will map (in the Heisenberg picture) the annihilation operator of mode i into a linear combination of all other modeŝ In order to preserve the commutation relations of different modes, the matrix U must also be unitary U U † = U † U = N . Note that every annihilation operator is mapped to a linear combination of annihilation operators and thus annihilation and creation operators are not mixed. Because of this, multiport interferometers generate all passive (in the sense that no photons are created or destroyed) linear optics transformations. In a pioneering work by Reck et al. [67], it was shown that any multiport interferometer can be realized using N (N − 1)/2 beamsplitter gates distributed over 2N − 3 layers. Recently this result was improved upon by Clements et al. [63], who showed that an equivalent decomposition can be achieved with the same number of beamsplitters but using only N + 1 layers.

Gaussian Operations
As mentioned in the previous subsections, generators that are at most quadratic remain closed under composition of their associated unitaries. In the Heisenberg picture these quadratic generators will produce all possible linear transformations between the quadrature (or mode) operators,â These operations are known as Gaussian operations. All gates in Table III  An important result for the CV formalism is that Gaussian quantum computation, i.e., computation that occurs with Gaussian states, operations and measurements, can be efficiently simulated on a classical computer (this is the foundation for the Gaussian backend in Strawberry Fields). This result is the CV version [69] of the Gottesman-Knill theorem of discrete-variable quantum information [49]. Hence we need non-Gaussian operations in order to achieve quantum supremacy in the CV model. Interestingly, even in the restricted case where all states and gates are Gaussian, with only the final measurements being non-Gaussian, there is strong evidence that such a circuit cannot be efficiently simulated classically [10,70]. More discussion and example code for this situation (known as Gaussian boson sampling) is provided in Appendix C.

Symplectic Formalism
In this section we review the symplectic formalism which lies at the heart of the Gaussian backend of Strawberry Fields. The symplectic formalism is an elegant and compact description of Gaussian states in terms of covariance matrices and mean vectors [50,68]. To begin, the commutation relations of the 2N position and momentum operators of Eq. (A2) can be easily summarized as [r i ,r j ] = ħ hiΩ i j where is the symplectic matrix. Using the symplectic matrix, we can define the Weyl operator D(ξ) (a multimode displacement operator) and the characteristic function χ(ξ) of a quantum N mode state ρ: where ξ ∈ 2N . We can now consider the Fourier transform of the characteristic function to obtain the Wigner function of the stateρ The 2N real arguments r of the Wigner function are the eigenvalues of the quadrature operators fromr.
The above recipe maps an N -mode quantum state living in a Hilbert space to the real symplectic space K := ( 2N , Ω), which is called phase space. The Wigner function is an example of a quasiprobability distribution. Like a probability distribution over this phase space, the Wigner function is normalized to one; however, unlike a probability distribution, it may take negative values. Gaussian states have the special property that their characteristic function (and hence their Wigner function) is a Gaussian function of the variables r. In this case, the Wigner function takes the form wherer = 〈r〉 ρ = Tr (rρ) is the displacement or mean vector and V i j = 1 2 〈∆r i ∆r j + ∆r i ∆r j 〉 ρ with ∆r =r −r. Note that the only pure states that have non-negative Wigner functions are the pure Gaussian states [58].
Each type of Gaussian state has a specific form of covariance matrix V and mean vectorr. For the single-mode vacuum state, we have V = ħ h 2 2 andr = (0, 0) T . A thermal state (Eq. (8)) has the same (zero) displacement but a covariance matrix V = (2n + 1) ħ h 2 2 , wheren is the mean photon number. A coherent state (Eq. (6)), obtained by displacing vacuum, has the same V as the vacuum state but a nonzero displacement vectorr = 2 ħ h 2 (Re(α), Im(α)). Lastly, a squeezed state (Eq. (7)) has zero displacement and covariance matrix V = ħ h 2 diag(e −2r , e 2r ). In the limit r → ∞, the squeezed state's variance in thex quadrature becomes zero and the state becomes proportional to thexeigenstate |x〉 with eigenvalue 0. Consistent with the uncertainty principle, the squeezed state's variance inp blows up. We can also consider the case r → −∞, where we find a state proportional to the eigenstate |p〉 of thep quadrature with eigenvalue 0. In the laboratory and in numerical simulation we must approximate every quadrature eigenstate using a finitely squeezed state (being careful that the variance of the relevant quadrature is much smaller than any other uncertainty relevant to the system). Any other quadrature eigenstate can be obtained from the x = 0 eigenstate by applying suitable displacement and rotation operators. Finally, note that Gaussian operations will transform the vector of means via affine transformations and the covariance matrix via congruence transformations; for a detailed dis-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 interferometers) should be performed to prepare the corresponding Gaussian state? To answer the first question: a 2N × 2N real symmetric matrixṼ corresponds to a Gaussian quantum state if and only ifṼ + i ħ h 2 Ω ≥ 0 (the matrix inequality is understood in the sense that the eigenvalues of the quantityṼ + i ħ h 2 Ω are nonnegative). The answer to the second question is provided by the Bloch-Messiah reduction [60][61][62]. This reduction shows that any N -mode Gaussian state (equivalently any covariance matrix and vector of means) can be 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 S i (z i ), followed by another multiportV , followed by single-mode displacement operations i D i (α i ). Explicitly, 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 vacuum to vacuum again. Thus for pure Gaussian states we need only generate N single-mode squeezed states and send them through a single multiport interferometerV 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) n i = (ν 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]).

Appendix B: Strawberry Fields Operations
In this Appendix, we present a complete list of the CV states, gates, and measurements available in Strawberry Fields.

Arbitrary Fock basis state
Prepare an arbitrary multi-mode mixed state, represented by a density matrix array x in the Fock basis.

Gaussian state preparation
Prepares an arbitrary N -mode Gaussian state defined by covariance matrix V ∈ 2N ×2N and means vector µ ∈ 2N using the Williamson decomposition

Gaussian transformation
Applies an N -mode Gaussian transformation defined by symplectic matrix S ∈ 2N ×2N using the Bloch-Messiah decomposition

Multi-mode linear interferometer
Applies an N -mode interferometer defined by unitary matrix U ∈ N ×N using the Clements decomposition

Operation
Name Definition

Heterodyne measurement
Projects the state onto the coherent states, sampling from the joint Husimi distribution 1 π 〈α|ρ|α〉

Gate Decompositions
In addition, the Strawberry Fields frontend can be used to provide decompositions of certain compound gates. The following gate decompositions are currently supported.

Two-mode squeeze gate
The two-mode squeeze gate S2gate(z) is decomposed into a combination of beamsplitters and single-mode squeezers

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

Appendix C: Quantum Algorithms
In this Appendix, we present full example code for several important algorithms, subroutines, and use-cases for Strawberry Fields. These examples are presented in more detail in the online documentation located at strawberryfields.readthedocs.io.

Quantum Teleportation
Quantum teleportation -sometimes referred to as state teleportation to avoid confusion with gate teleportationis the reliable transfer of an unknown quantum state across spatially separated qubits or qumodes, through the use of a classical transmission channel and quantum entanglement [73]. Considered a fundamental quantum information protocol, it has applications ranging from quantum communication to enabling distributed information processing in quantum computation [74].
In general, all quantum teleportation circuits work on the same basic principle. Two distant operators, Alice and Bob, share a maximally entangled quantum state (in discrete variables, any one of the four Bell states, and in CVs, a maximally entangled state for a fixed energy), and have access to a classical communication channel. Alice, in possession of an unknown state which she wishes to transfer to Bob, makes a joint measurement of the unknown state 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 of the entangled state to an accurate replica of the original unknown state, by performing a conditional phase flip (for qubits) or displacement (for qumodes) [75]. The CV teleportation circuit is shown in Fig. 3. 1 import strawberryfields as sf 2 from strawberryfields.ops import * 3 from strawberryfields.utils import scale 4 from numpy import pi, sqrt

Gate Teleportation
In the quantum state teleportation algorithm discussed in the previous section, the quantum state is transferred from the sender to the receiver. However, quantum teleportation can be used in a much more powerful manner, by simultaneously transforming the teleported state -this is known as gate teleportation.
In gate teleportation, rather than applying a quantum unitary directly to the state prior to teleportation, the unitary is applied indirectly, via the projective measurement of the first subsystem onto a particular basis. This measurement-based approach provides significant advantages over applying unitary gates directly, for example by reducing resources, and in the application of experimentally hard-to-implement gates [74]. In fact, gate teleportation forms a universal quantum computing primitive, and is a precursor to cluster state models of quantum computation [76].
First described by Gottesman and Chuang [77] in the case of qubits, gate teleportation was generalised for the CV case by Bartlett and Munro [78] (see Fig. 4). In an analogous process to the discrete-variable case, we begin with the algorithm for local state teleportation. Unlike the spatiallyseparated quantum state teleportation we considered in the previous section, local teleportation can transport the state using only two qumodes; the state we are teleporting is entangled directly with the squeezed vacuum state in the momentum space through the use of a controlled phase gate.

FIG. 4: Gate teleportation of a quadratic phase gate P onto input state |ψ〉.
To recover the teleported state exactly, we must perform Weyl-Heisenberg corrections to the second mode; here, that would be F † X (m) † , where m is the correction based on the Homodyne measurement. However, for convenience and simplicity, it is common to simply write the circuit without the corrections applied explicitly. 1 import strawberryfields as sf 2 from strawberryfields.ops import *   Note that additional gates can be added to the gate teleportation scheme described above simply by introducing additional qumodes with the appropriate projective mea- surements, all 'stacked vertically' (i.e., coupled to each consecutive qumode via a conditional phase gate). It is from this primitive that the model of cluster state quantum computation can be derived [76].

Boson Sampling
Introduced by Aaronson and Arkhipov [9], boson sampling presented a slight deviation from the general approach in quantum computation. Rather than presenting a theoretical model of universal quantum computation (i.e., a framework that enables quantum simulation of any arbitrary Hamiltonian [51]), boson sampling-based devices are instead an example of an intermediate quantum computer, designed to experimentally implement a computation that is intractable classically [79][80][81][82].
Boson sampling proposes the following quantum linear optics scheme. An array of single photon sources is set up, designed to simultaneously emit single photon states into a multimode linear interferometer; the results are then generated by sampling from the probability of single photon measurements from the output of the linear interferometer.
For example, consider N single photon Fock states, |ψ〉 = |m 1 , m 2 , . . . , m N 〉, composed of b = i m i photons, incident on an N -mode linear interferometer described by the unitary U acting on the input mode creation and annihilation operators. It was shown that the probability of detecting n j photons at the jth output mode is given by [9] Pr(n 1 , . . . , n N ) = |Per(U st )| 2 m 1 ! · · · m N !n 1 ! · · · n N ! ; i.e., the sampled single photon probability distribution is proportional to the permanent of U st , a submatrix of the interferometer unitary, dependent upon the input and output Fock states. Whilst the determinant can be calculated efficiently on classical computers, calculation of the permanent belongs to the computational complexity class #P-Hard problems [83], which are strongly believed to be classically hard to calculate. This implies that simulating boson sampling is an intractable task for classical computers, providing an avenue for the demonstration of quantum supremacy. Continuous-variable quantum computation is ideally suited to the simulation and demonstration of the boson sampling scheme, due to its grounding in quantum optics. In quantum linear optics, the multimode linear interferometer is commonly decomposed into two-mode beamsplitters (BSgate) and single-mode phase shifters (Rgate) [67], allowing for a straightforward translation into a CV quantum circuit. In order to allow for arbitrary linear unitaries on m imput modes, we must have a minimum of m + 1 columns in the beamsplitter array [63].

Gaussian Boson Sampling
While boson sampling allows the experimental implementation of a sampling problem that is hard to simulate classically, one of the setbacks in experimental setups is scalability, due to its dependence on an array of simultaneously emitting single photon sources. Currently, most physical implementations of boson sampling make use of a process known as Spontaneous Parametric Down-Conversion (SPDC) to generate the single-photon source inputs. However, this method is non-deterministic -as the number of modes in the apparatus increases, the average time required until every photon source emits a simultaneous photon increases exponentially.
In order to simulate a deterministic single photon source array, several variations on boson sampling have been proposed; the most well-known being scattershot boson sampling [70]. However, a recent boson sampling variation by [10] negates the need for single photon Fock states altogether, by showing that incident Gaussian states -in this case, single mode squeezed states -can produce problems in the same computational complexity class as boson sampling. Even more significantly, this mitigates the scalability problem with single photon sources, as single mode squeezed states can be simultaneously generated experimentally.
With an input ensemble of N single-mode squeezed states with squeezing parameter z = r ∈ , incident on a linearinterferometer described by unitary U, it can be shown that the probability of detecting an output photon pattern (n 1 , . . . , n N ), where n k ∈ {0, 1}, is given by [10] Pr(n 1 , . . . , where S denotes the subset of modes where a photon was detected and Haf[·] is the Hafnian [84,85]. That is, the sampled single photon probability distribution is proportional to the hafnian of a submatrix of U i tanh(r i )U T . The hafnian is known to be a generalization of the permanent, and can be used to count the number of perfect matchings of an arbitrary graph [86]. The formula above can be generalized to pure Gaussian states with finite displacements by using the loop hafnian function which counts the number of perfect matchings of graphs with loops [56]. Since any algorithm that could calculate the hafnian could also calculate the permanent, it follows that calculating the hafnian remains a classically hard problem; indeed, the best known classical algorithm for the calculation of a hafnian of an arbitrary symmetric complex matrix of size N scales like

Instantaneous Quantum Polynomial (IQP)
Like boson sampling and Gaussian boson sampling before it, the instantaneous quantum polynomial (IQP) protocol is a restricted, non-universal model of quantum computation, designed to implement a computation that is classically intractable and verifiable by a remote adjudicator. First introduced by Shepherd and Bremner [90], IQP circuits are defined by the following conditions: (i) there are N inputs in state |0〉 acted on by Hadamard gates, (ii) the resulting computation is diagonal in the computational basis by randomly selecting from the set U = {R(π/4), C Z} (hence the term 'instantaneous'; since all gates commute, they can be applied in any temporal order), and (iii) the output qubits are measured in the Pauli-X basis. Efficient classical sampling of the resulting probability distribution H ⊗N U H ⊗N |0〉 ⊗Neven approximately [13] or in the presence of noise [91] has been shown to be #P-hard, and would result in the collapse of the polynomial hierarchy to the third level [12,92].
however, the IQP protocol was not constructed with quantum linear optics in mind. Nevertheless, the IQP model was recently extended to the CV formulation of quantum computation by Douce et al. [17], taking advantage of the ability to efficiently prepare large resource states, and the higher efficiencies afforded by homodyne detection. Furthermore, the computational complexity results of the discrete-variable case apply equally to the so-called CV-IQP model, assuming a specific input squeezing parameter dependent on the circuit size. The CV-IQP model is defined as follows: 1. inputs are squeezed momentum states |z〉, where z = r ∈ and r < 0; 2. the unitary transformation is diagonal in thex quadrature basis, by randomly selecting from the set of gates U = {Z(p), C Z(s), V (γ)}; 3. and finally, homodyne measurements are performed on all modes in thep quadrature.
1 import strawberryfields as sf 2 from strawberryfields.ops import * 3 4 # initialise the engine and register 5 eng, q = sf.Engine(4) Moreover, the resulting probability distributions have been shown to be given by integrals of oscillating functions in large dimensions, which are believed to be intractable to compute by classical computers. This leads to the result that even in the case of finite squeezing and reduced measurement precision, approximate sampling from the output CV-IQP model remains classically intractable [17,93].

Hamiltonian Simulation
The simulation of atoms, molecules and other biochemical systems is another application uniquely suited to quantum computation. For example, the ground state energy of large systems, the dynamical behaviour of an ensemble of molecules, or complex molecular behaviour such as protein folding, are often computationally hard or downright impossible to determine via classical computation or experimentation [94,95].
In the discrete-variable qubit model, efficient methods of Hamiltonian simulation have been discussed at length, providing several implementations depending on properties of the Hamiltonian, and resulting in a linear simulation time [96,97]. Efficient implementations of Hamiltonian simulation also exist in the CV formulation [98], with specific application to Bose-Hubbard Hamiltonians (describing a system of interacting bosonic particles on a lattice of orthogonal position states [99]). As such, this method is ideally suited to photonic quantum computation. Considering a lattice composed of two adjacent nodes characterised by Adjacency matrix A, the Bose-Hubbard Hamiltonian is given by = J(â † 1â 2 +â † 2â 1 ) + 1 2 U(n 2 1 −n 1 +n 2 2 −n 2 ), where J represents the transition of the boson between nodes, and U is the on-site interaction potential. Applying the Lie product formula, we find that where O t 2 /k is the order of the error term, derived from the Lie product formula. Comparing this to the form of various gates in the CV circuit model, we can write this as the product of beamsplitters, Kerr gates, and rotation gates: where θ = −J t/k, φ = π/2, and r = −U t/2k. Using J = 1, U = 1.5, k = 20, and a timestep of t = 1.086, this can be easily implemented in Strawberry Fields using only 20×3 = 60 gates.
1 import strawberryfields as sf 2 from strawberryfields.ops import * 3 from numpy import pi  For more complex Hamiltonian CV decompositions, including those with nearest-neighbour, see Kalajdzievski et al. [98]. This decomposition is also implemented in the SFOpenBoson plugin [45], which provides an interface between OpenFermion, the quantum electronic structure package, and Strawberry Fields. This allows arbitrary Bose-Hubbard Hamiltonians, generated in OpenFermion, to be simulated using Strawberry Fields.

Optimization of Quantum Circuits
One of the unique features of Strawberry Fields is that it has been designed from the start to support modern computational methods like automatic gradients, optimization, and machine learning by leveraging a TensorFlow backend. We have already given an overview in the main text of how these features are accessed in Strawberry Fields. We present here a complete example for the optimization of a quantum circuit. Our goal in this circuit is to find the Dgate parameter which leads to the highest probability of a n = 1 Fock state output. This simple baseline example can be straightforwardly extended to much more complex circuits. As optimization is a key ingredient of machine learning, this example can also serve as a springboard for more advanced data-driven modelling tasks.
We note that optimization here is in a variational sense, i.e., we choose an ansatz for the optimization by fixing the discrete structure of our circuits (the selection and arrangement of specific states, gates, and measurements). The optimization then proceeds over the parameters of these operations. Finally, we emphasize that for certain circuits and objective functions, the optimization might be non-convex in general. Thus we should not assume (without proof) that the solution obtained via a Strawberry Fields optimization is always the global optimum. However, it is often the case in machine learning that local optima can still provide effective solutions, so this may not be an issue, depending on the application. 1 import strawberryfields as sf 2 from strawberryfields.ops import * 3 import tensorflow as tf