Perceval: A Software Platform for Discrete Variable Photonic Quantum Computing

We introduce Perceval, an open-source software platform for simulating and interfacing with discrete-variable photonic quantum computers, and describe its main features and components. Its Python front-end allows photonic circuits to be composed from basic photonic building blocks like photon sources, beam splitters, phase-shifters and detectors. A variety of computational back-ends are available and optimised for different use-cases. These use state-of-the-art simulation techniques covering both weak simulation, or sampling, and strong simulation. We give examples of Perceval in action by reproducing a variety of photonic experiments and simulating photonic implementations of a range of quantum algorithms, from Grover's and Shor's to examples of quantum machine learning. Perceval is intended to be a useful toolkit for experimentalists wishing to easily model, design, simulate, or optimise a discrete-variable photonic experiment, for theoreticians wishing to design algorithms and applications for discrete-variable photonic quantum computing platforms, and for application designers wishing to evaluate algorithms on available state-of-the-art photonic quantum computers.


Introduction
Quantum computing has gained huge interest and momentum in recent decades because of its promise to deliver computational advantages and speedups compared to the classical computing paradigm. The essential idea is that information processing takes place on physical devices, and if the components of those devices behave according to the laws of quantum rather than classical physics then it opens the door to exploiting quantum effects to process information in radically different and potentially advantageous ways.
Quantum algorithms like Shor's factorisation algorithm [1], which gives an exponential speedup over its best known classical counterpart, or Grover's search algorithm [2], which gives a quadratic speedup over its classical counterparts, are often cited, and have captured the imagination and provided motivation for the rapid developments that have taken place in quantum technologies. However, for such algorithms to have practical significance will require large-scale fault-tolerant quantum computers. Yet the quantum devices that are available commercially or in research laboratories today are somewhat more limited and belong, for now at least, to the so-called noisy-intermediate scale quantum (NISQ) regime [3].
Of course NISQ devices are a necessary step on the path to large-scale fault-tolerant quantum computers, but in principle they could already enable quantum computational advantages [4] -in which they could outperform even the most powerful classical supercomputers available today -for tasks with practical relevance. Several experiments have already claimed to demonstrate quantum computational advantage in sampling tasks [5][6][7][8][9]. The practical significance of these is still being explored [10][11][12], but meanwhile many other promising proposals for quantum algorithms that may deliver practical advantages in this regime are also being pursued. Among others, these include the quantum variational eigensolver [13], universal function approximators like that of [14] (both of which we will return to later in this paper), the quantum approximate optimisation algorithm [15], as well as a host of other algorithms [16] with applications that range from chemistry [17,18] and many-body physics [19,20] to combinatorial optimisation [21,22] and machine learning [23,24].
In fact, a number of technological routes to building quantum computers are being actively pursued, in which quantum information is encoded in very different kinds of physical systems. These include matter-based approaches that rely on superconducting circuits, cold atoms, or trapped ions, but also light-based approaches in which photons are the basic information-carrying systems. Among these, photons have a privileged status, as they are the natural and indeed only viable support for communicating quantum information, which will eventually be required for networking quantum processors and devices. As such they will necessarily be a part of any longer-term developments in quantum computational hardware and infrastructure. More than this, photons provide viable routes to both NISQ [25] and large-scale fault-tolerant quantum computing through measurement-based models [26,27] in their own right.
Perceval is a complete and efficient software platform for the discrete variable (DV) model for photonic quantum computing. 1 It especially uses Fock state descriptions of photons generated by sources, evolving through linear optical networks -composed for example of beam splitters, phase shifters, waveplates or other linear optical components -and then being detected. The familiar qubit and measurement-based models for quantum computing can be encoded within the DV model (see e.g. [25,26]). However, the DV model is also of significant interest in its own right -not least because it has led to some of the first proposals [28] for quantum computational advantage to be demonstrated with NISQ devices. This concerns a computational problem known as Boson Sampling, discussed in detail in Section 4.2, which essentially consists of sampling from the output distribution of photons that interfere in a generic interferometer. We will also demonstrate direct DV photonic approaches to quantum machine learning in Section 4.6.
Perceval 's features can be useful for designing, optimising, simulating, and eventually transpiling DV linear optical circuits and executing them on cloud-based physical processors. Although Perceval provides bridges (see Section 3.2.6) with other open-source quantum computing toolkits [30] such as Qiskit [31], it further allows users to work at a level that is closer to the photonic hardware than these toolkits regular gate-based qubit quantum circuit model, which is already the focus of a number of other software platforms. It is also more flexible and supports many more functionalities than existing packages for linear optical quantum systems [32]. The finer-grained control of the physical hardware can be valuable both for the NISQ regime, where it is crucial to achieve the maximum performance from the specific hardware resources available, and for optimising the building block photonic modules in schemes for reaching the large-scale fault-tolerant regime.
Perceval is intended to be useful for experimentalists wishing to design photonic experiments, including allowing for realistic modelling of noise and imperfections, and for computer scientists and theoreticians seeking to develop algorithms and applications for photonic quantum computers. It is an open source platform that is intended for community development via GitHub, the project forum website and an updated documentation.
Perceval integrates several state-of-the-art algorithms for running simulations optimised with low-level single instruction, multiple data (SIMD) implementations, allowing users to push close to the limits of classical simulability with desktop computers. Extensions to the framework are also intended for high-performance computing (HPC) cluster deployment which can permit simulation to scale further. Since version 0.7, Perceval is able to run samples from real photonic chips through the Cloud, giving the user access in real-time to the exact hardware characteristics of the photonic source (brightness, purity, indistinguishability), the chip, and the detectors. Remote computers are also available if the user wants to perform classical simulation with greater computational power. 2 The remainder of this white paper is structured as follows. We provide some brief background on photonic quantum computing in Section 2, before outlining the structure and key features of Perceval in Section 3, and then go on to give a number of illustrative examples of Perceval in action in Section 4. The code of the examples can be found in Appendix B. This paper is based on version 0.7.3 of Perceval .

Photonic Quantum Computing
Similar to the qubit quantum circuit model, the DV photonic model can be presented as a gate-based model, in which states are prepared, transformed as they are acted upon by gates, and measured.
We consider a number m ∈ N of spatial modes. Physically these could correspond to waveguides in an integrated circuit, optical fibres, or paths in free space. Photon sources prepare initial states in these spatial modes. These are number states |n , where n is the number of photons in the mode, or superpositions of number states. We use the shorthand notation |0, 1 for a two-mode system with state |0 in the first and |1 in the second, etc. Unless otherwise specified, photons are assumed to be indistinguishable. Sometimes we will wish to keep track of the polarisation of the photons, which can be achieved by further splitting each spatial mode into two polarisation modes, e.g. horizontal (denoted by |H ) and vertical (denoted by |V ), and recording the (superpositions of) number states for each. Perceval also allows the possibility of tracking other attributes of photons.
Transformations are performed on the states by evolving them through linear optical networks. The simplest linear optical operations (gates) are: phase shifters, which act on a single spatial mode, and beam splitters, which act on pairs of spatial modes. These operations preserve photon number and are best described as unitary matrices that act on the creation operator of each mode (see [33,34] for a more detailed treatment). A creation operator is defined by its action a † |n = √ n + 1 |n + 1 . The unitary associated with a phase shifter P φ with phase parameter φ ∈ [0, 2π] is simply the scalar e iφ , while in the Rx convention the unitary matrix associated to a beam splitter B θ is given in Equation 1. We have included for clarity the basis on which the matrix acts in gray above and on the left.
where the parameter θ relates to the reflectivity and |1, 0 denotes the state in which the photon travels in the first spatial mode. Perceval introduces all theoretically equivalent beam splitter matrix conventions as shown in Table 1. The action of any linear optical circuit is thus given by the unitary matrix obtained by composing and multiplying the matrices associated with its elementary components. Interestingly, it can be shown that conversely any unitary evolution can be decomposed into a combination of beam splitters and phase shifters, e.g. in a 'triangular' [35] or 'square' [36] array. Perceval also allows for swaps or permutations of spatial modes, and can include operations like waveplates, which act on polarisation modes of a single spatial mode and are described by the following unitary [37]: Here δ is a parameter proportional to the thickness of the waveplate and ξ represents the angle of the waveplate's optical axis in the {|H , |V } plane. Especially important is the case that δ = π/2, known as a half-wave plate, which rotates linear polarisations in the {|H , |V } plane. Quarter-wave plates (δ = π/4) convert circular polarisations to linear ones (e.g. |L = (|H + i |V )/ √ 2) and vice-versa. Polarising beam splitters will convert a superposition of polarisation modes in a single spatial mode to the corresponding equal-polarisation superposition of two spatial modes, and vice versa, and so in this sense allow us to translate between polarisation and spatial modes. The unitary matrix associated to a polarising beam splitter acting on the tensor product of the spatial mode and the polarisation mode is where |H, 0 denotes the single-photon state in which the photon travels in the first spatial mode with a horizontal polarisation. Similarly then, any unitary evolution on spatial and polarisation modes can be decomposed into an array of polarising beam splitters, beam splitters and phase shifters.
Finally measurement, or readout, is made by single photon detectors (these may be partially or fully number resolving). Both photon generation and measurement are nonlinear operations. In particular, measurements can be used to induce probabilistic or postselected non-linearities in linear optical circuits. This could be used to further implement feedforward, whereby measurement events condition parameters further along in a circuit.
The basic post-selection-, feedforward-and polarisation-free (without polarisation modes or operations) fragment of the DV model described above is precisely the model considered by Aaronson and Arkhipov in [28]. For classical computers, simulating this simple model, straightforward though its description may be, is known to be #P -hard, and up to complexity-theoretic assumptions the same is true for approximate classical simulation. The reason for this essentially comes down to the fact that calculating the detection statistics requires evaluating the permanents of unitary matrices, a problem which itself is known to be #P-hard [38]. We describe this in a little more detail when we look at Boson Sampling in Section 4.2.
Similarly the polarisation-free fragment was shown by Knill, Laflamme and Milburn to be quantum computationally universal [25]. In particular, it can be used to represent qubits and qubit logic gates.
Just as the bit -any classical two level system -is generally taken as the basic informational unit in classical computer science, the qubit -any two level quantum system -is usually taken as the basic informational unit in quantum computer science. Photons have many degrees of freedom and offer a rich variety of ways to encode qubits. One of the most common approaches is to use the dual-rail path encoding of [25]. Each qubit is encoded by one photon which may be in superposition over two spatial modes, which correspond to the qubit's computational basis states: Single-qubit unitary gates are particularly straightforward to implement in this encoding, requiring only a fully parametrisable beam splitter, as in Equation 1, and a phase shifter. An example of an entangling two-qubit gate (CNOT), which uses heralding, i.e. postselection over auxiliary modes, will be presented in Section 3. Another common qubit encoding is the polarisation encoding, and indeed versions of the KLM scheme also exist for this encoding [39]. Here each qubit is encoded in the polarisation degree of freedom of a single photon. This encoding will be used in the example of Section 4.3. In practice, encodings may be chosen based on the availability or ease of implementation of different linear optical elements in laboratory settings. For instance path-encoding is at present more accessible than polarisation encoding in integrated optics, where linear optical circuits are implemented 'on-chip'. Although we have discussed ways of encoding qubit quantum circuits into the DV model, it should be reiterated that the DV model has interesting features in its own right, and the examples of Boson Sampling in Section 4.2 and quantum machine learning in Section 4.6 give some illustrations of interesting applications that bypass qubit descriptions.
We have presented the idealised DV model, but it is important to note that in realistic implementations various noise, imperfections, and errors will arise. Perceval is also intended as a tool for designing, modelling, simulating and optimising realistic DV linear optical circuits and experiments, and to incorporate realistic noise-models. We briefly demonstrate how Perceval handles imperfect photon purity at source and distinguishability due to imperfect synchronisation in Section 4.1. Other relevant factors are losses at all stages from source to detection, imperfect or imperfectly characterised linear optical components, cross-talk effects, detector dark-counts, etc. Future versions of Perceval will contain increasingly sophisticated and parametrised modelling of such factors.

Global Architecture
Perceval is a linear optical circuit development framework whose design is based on the following core ideas: • It is simple to use, both for theoretical and experimental physicists and for computer scientists.
• It does not constrain the user to any framework-specific conventions that are theoretically equivalent (for example it was noticed early on that many different conventions for beam splitters can be found in the literature).
• It provides state-of-the-art optimised algorithms -as benchmarked on specific use cases.
• It provides -when possible and appropriate -access to symbolic calculations for finding analytical solutions.
• It provides companion tools, such as a unitary matrix toolkit, and L A T E Xor HTML rendering of algorithms.
• It incorporates realistic, parameterisable error and noise modelling.
• It aims to provide a seamless transition from simulators to actual photonic processors (QPU). As such, most of the programmatic interfaces are designed for QPU control.
In particular, fixing a specific QPU automatically hides methods giving access to properties (such as probability amplitudes) and features which would be unavailable on the actual hardware.
Perceval is a modular object-oriented Python code, with optimised functions written in C making use of SIMD vectorisation. In the following section, we give an overview of the main classes available to the user. A full documentation is maintained on GitHub and available online through the project website.

States
Information in a linear optical circuit is encoded in the state of photons in certain "modes" that are defined by the circuit designer. States are implemented in Perceval by the following two classes: • BasicState is used to describe Fock states of n photons over m modes.By default, photons are indistinguishable but each photon can be annotated, controlling its distinguishability. An annotation is a way to associate additional information to each photon and can thus represent additional degrees of freedom -for instance, polarisation is represented as a photon annotation; • StateVector extends BasicState to represent superpositions of states.

Circuits
The Circuit class provides a practical way to assemble a linear optical circuit from predefined elementary components, other circuits, or unitary matrices. It can also compute the unitary matrix associated to a given circuit, or conversely decompose a unitary into a linear optical circuit consisting of user-defined components. A library of predefined elementary components is provided (see Table 1). Those components can be added on the desired mode(s) on the right of a Circuit with the Circuit.add(modes, component) method.
An example of code using this class can be found in Code 1.

Name Unitary Matrix Representation
Beam Splitter R X convention:

Back-ends
Perceval allows the user to choose between its four back-ends -CliffordClifford2017, Naive, SLOS and Stepper -each one taking a different computational approach to circuit simulation. They perform the following tasks: • Sample individual single output states -CliffordClifford2017; • Compute the probability, or probability amplitude, of obtaining a given output state from a given input state -Naive; • Describe the exact complete output state -SLOS, Stepper.
The CliffordClifford2017 Back-end. This back-end is the implementation of the algorithm introduced in [40]. The algorithm, applied to Boson Sampling, aims to "produce provably correct random samples from a particular quantum mechanical distribution". Its time and space complexity are respectively O n2 n + mn 2 and O(m). The algorithm has been implemented in C++, and uses an adapted Glynn algorithm [41] to efficiently compute n simultaneous "sub-permanents". Recently, the same authors have proposed a faster algorithm in [42] with an average time complexity of O(nρ n θ ) for a number of modes m = θn which is linear in the number of photons n, where: For example, taking θ = 2, which corresponds to dual rail path encoding without auxiliary modes, the average performance of this algorithm is O n( 5 5 8 2 3 3 ) n ≈ O(n1.8 n ) as opposed to O(n2 n ) for the original algorithm of [40]. The implementation of [42] in Perceval is in progress.
The Naive Back-end. This back-end implements direct permanent calculation and is therefore suited for single output probability computation with small memory cost. Both Ryser's [43] and Glynn's [41] algorithms have been implemented. Extra care has been taken on the implementation of these algorithms, with usage of different optimisation techniques including native multithreading and SIMD vectorisation primitives. A benchmark of these algorithms against the implementation present in the The Walrus library [44] is provided in Figure 1.
The SLOS Back-end. The Strong Linear Optical Simulation SLOS algorithm developed by a subset of the present authors is introduced in [45]. It unfolds the full computation path in memory, leading to a remarkable time complexity of O n m+n−1 n for computing the full distribution. The current implementation also allows restrictive sets of outputs, with average computing time in O(nρ n θ ) for single output computation. As discussed in [45], it is possible to use the SLOS algorithm in a hybrid manner that can combine both weak and strong simulation, though it has not yet been implemented in the current version of Perceval . The tradeoff in the SLOS algorithm is a huge memory usage of O n m+n−1 n that limits usage to circuits with ≈ 20 photons on personal computers and with ≈ 24 photons on super-computers. The processor is a 32 core, 3.1GHz Intel Haswell. For The Walrus, version 0.19 is used and installed from pypi. The Ryser implementation is run on 4 or 32 threads. The Glynn implementation is run on a single thread. What is interesting to note is that all implementations have convergence to the theoretical performance but the factor between optimised and less optimised implementation still makes a perceptible time difference for the end-user. Based on different behaviour between Ryser and Glynn with n and potential multi-threading, Perceval has some built-in logic to switch between the two algorithms.

The
Stepper Back-end. This back-end takes a completely different approach. Without computing the circuit's overall unitary matrix first, it applies the unitary matrix associated with the components (see Table 1) in each layer of the circuit one-by-one, simulating the evolution of the state vector. The complexity of this back-end is therefore proportional to the number of components. It has the nice features that: • it can support more complex components like Time Delay; • it is very flexible with simulating noise in the circuit, like photon loss, or more generally with simulating any non-linear operation the user would wish to implement; • it simplifies the debugging of circuits by exposing intermediate states.
The Stepper back-end is really meant for circuit simulation with losses or non linear components. In all other cases, it is more efficient to directly consider the unitary matrix of the whole circuit and compute the output state with SLOS, instead of performing a computation for each component using the Stepper back-end.
Theoretical performances and specific features of the different back-ends are summarised in Table 2.

Processors
The Processor class allows the user to emulate a real photonic quantum processor, taking into account the single photon source, the circuit and the detectors. The Perceval source model allows tuning of the transmittance, multiple photon emission probability, and indistinguishability. By default, the source is perfect. The circuit can be composed of unitary or non-unitary components. Detector imperfections can be emulated -one can switch between number resolving detectors and threshold detectors. Post-selection features can   [42] is in progress; p is a polynomial function. (**) An implementation of SLOS Sampling and single output efficiency is in progress. (***) Practical limits are subjective and corresponding to a memory usage < 16Gb, and a usage time for a given function of less than a few seconds. For Stepper, it is hard to evaluate exactly the complexity since it is really proportional to the number of "components" (N c ) and the size of the output space that are circuit-specific.
also be added to a Processor via heralded modes and/or a final post-selection function. See the Appendix A.4 for more details.

Algorithms
The algorithm library provides a set of tasks which can be performed on a Processor. These tasks can be as simple as obtaining a sample result (see Appendix A.4), or slightly more complex Analyzer will output a probability table of expected input-output correlations, as illustrated in Figure 4.

Bridges to Other Quantum Computing Toolkits
In order to facilitate work across multiple platforms, Perceval provides both a catalog of predefined components and specific converters: • The predefined component catalog is available in Perceval.components.catalog, where each element is an instance of CatalogItem providing a ready-to-use circuit (as_circuit() method), a processor (as_processor() method), and a complete documentation (doc property).
• Converters to Perceval from other frameworks are available in perceval.converters. Their task is to provide a photonic Perceval equivalent of circuits and algorithms defined in other open source frameworks. Currently a qiskit 4 converter is available and converters for tket and QLM are in development.
Examples for both tools are provided in a notebook of the Perceval documentation.

Step-by-Step Example
The following code implements a simple linear optical circuit corresponding to a pathencoded CNOT gate (after post-selection in the coincidence basis) [46]. Code 1: Implementation of the CNOT of [46] Let us explain how it works. We remind that the CNOT gate is a two-qubit gate that acts in the following way: i.e. it's a controlled-not gate where the first qubit acts as control and the second qubit is the target. We use dual-rail encoding (see Equation 4) and thus need four spatial modes to represent the target and the control. We additionally need two auxiliary empty modes to implement this post-selected linear optical CNOT. The circuit is made of three central beam splitters of reflectivity 2/3, whereas the two to the left and right of the central ones have a reflectivity of 1/2. The corresponding circuit displayed by Perceval is depicted in Figure 2. The first and last spatial modes correspond to the auxiliary empty modes, the second and third to the control qubit, and the fourth and fifth to the target qubit. The associated unitary matrix acting on the six spatial modes computed with Perceval is displayed in Figure 3. We then need to post-select on measuring two (dual-rail encoded) qubits in the modes corresponding to the target and the control, i.e. on measuring exactly one photon in the second or third spatial modes and exactly one photon in the fourth or fifth spatial modes. One can check that this measurement event happens with probability 1/9 -we call this number the 'performance' of a post-selected gate. We say that the error rate is 0 if the implementation of the gate is perfect (after post-selection). This computation can be done with Perceval : we define a processor, use the SLOS back-end, and perform a full output distribution and performance analysis, as illustrated in Figure 4. analyzer . compute ( expected = { " 00 " : " 00 " , " 01 " : " 01 " , " 10 " : " 11 " , " 11 " : " 10 " } ) pcvl . pdisplay ( analyzer , output_format = pcvl . Format . LATEX ) print ( " = > performance = %s , fidelity = % . 1f % % " % ( pcvl . simple_float ( analyzer . performance ) [ 1 ] , analyzer . fidelity * 100 ) ) Code 2: Analysis of the input-output map of the CNOT from Code 1

Perceval in Action
In this section we provide examples of the Perceval software in use. We reproduce several photonic experiments that implement important quantum algorithms and then demonstrate a photonic quantum machine learning algorithm on a simulated photonic quantum processor. The Hong-Ou-Mandel (HOM) effect [47] is an interference effect between pairs of indistinguishable photons, which, when incident on a balanced beam splitter via the respective input modes, will bunch and emerge together in either of the output modes with equal probability. In 2002, Santori et al. demonstrated the HOM effect with a unique source of single photons [48], evidencing the indistinguishability of consecutive photons and paving the way to many developments in experimental quantum optics.
In this experiment, "two pulses separated by 2 ns and containing 0 or 1 photons, arrive through a single-mode fibre. The pulses are interfered with each other using a Michelsontype interferometer with a (2 ns+∆t) path-length difference. [. . . ] The interferometer outputs are collected by photon counters, and the resulting electronic signals are correlated using a time-to-amplitude converter followed by a multi-channel analyser card, which generates a histogram of the relative delay time τ = t 2 − t 1 between a photon detection at one counter (t 1 ) and the other (t 2 )" [48]. An equivalent experiment was carried out in [49] with similar results, presented in Figure 5.  [49]. The red line shows the theoretical histogram for a perfect two-photon interference.

Perceval Implementation
The implementation on Perceval uses a simple circuit composed of a beam splitter (corresponding to the main beam splitter of [48, Figure 3(a)]), the lower output mode of which then passes through a 1-period time-delay, and then a second beam splitter on the modes where the interference happens.
The circuit and result are presented in Figure 6. The resolution is less fine-grained than the original experiments due to the absence of modelling of photon length in the current implementation of Perceval . However, we clearly recognise the same distinctive peaks in the relative amplitude. The very small coincidence rate at τ = 0 is the signature of photon interference. The implementation of the algorithm is provided in Appendix.

Introduction
Boson Sampling 5 is a sampling problem originally proposed by Aaronson and Arkhipov [28]. We give a detailed description of the problem below. It essentially consists of sampling from the probability distribution of outcome detection coincidences when single photons are introduced into a random linear optical circuit. Let m, n ∈ N * be positive integers with m ≥ n. Let S m,n be the set of all possible tuples Let U be an m × m Haar-random unitary matrix. From U we construct an n × n matrix U T,S as follows. For a given S = (s 1 , . . . , s m ) ∈ S m,n , first construct an n × m matrix U S by copying s i times the i th row of U . Then, for a given (fixed) T = (t 1 , . . . , t m ) ∈ S m,n , construct U T,S from U S by taking t j copies of the j th column of U S . Let Perm(U T,S ) be the permanent [51, Chapter 7] of U T,S , and let It can be shown that P (S) ∈ [0, 1], and that the set D U = {P (S) | S ∈ S m,n } is a probability distribution over outcomes S [28]. Let ε ∈ [0, 1] be a given precision. Approximate Boson Sampling can then be defined as the problem of sampling outcomes S from a probability distributionD such that where . denotes the total variational distance between two probability distributions. Exact Boson Sampling corresponds to the case in which ε = 0. For m n 2 , it was shown in [28] that: • No classical algorithm can solve Exact Boson Sampling in time poly(n), • No classical algorithm can solve Approximate Boson Sampling in time poly(n, 1 ε ), unless some widely accepted complexity theoretic conjectures turn out to be false.
These hardness results are closely related to the hardness of computing the permanent of matrices [28,38].
On the other hand, Boson Sampling can be solved efficiently in the exact case on a photonic quantum device which is noiseless [28], while the approximate version only requires sufficiently low noise levels [52][53][54]. This is done by passing n identical single photons through a lossless m-mode universal linear optical circuit [35,36] configured in such a way that it implements the desired Haar-random unitary transformation U . Universal here is meant in the sense of photonic unitaries and refers to the ability of the circuit to generate any arbitrary unitary evolution on the creation operators of the modes. A universal linear optical circuit can be configured to implement any m × m unitary chosen from the Haar measure, for example using the recipe of [55]. The input configuration of single photons corresponds to a tuple T . We then measure the output modes of the circuit using perfect single photon detectors, and will obtain the output configuration tuple S according to the distribution D U [28]. The detectors should be photon-number resolving in general, but not necessarily when working in the no-collision regime (n 2 << m) [28].

Perceval Implementation 6
Boson Sampling with single photons in a 60 mode linear optical circuit with up to 14photon coincidences at the outputs was reported in [56]. In this Section, we report the results of a noiseless Boson Sampling simulation performed with Perceval for n = 14 single photons and m = 60 modes. For this Boson Sampling, the total number of possible output states is M = m+n−1 n = 73 14 ≈ 10 14 [28]. For an input state of 14 photons in 14 arbitrarily chosen modes k 1 , ..., k 14 ∈ {0, ..., 59}, we generated 300 Haar random 60 × 60 unitaries, and performed for each of these unitaries 5 × 10 6 runs of Boson Sampling, where each run consists of sampling a single output. In total, we collected 1.5 × 10 9 samples.
The classical algorithm for Boson Sampling integrated into Perceval and which was used for this simulation is that of Clifford and Clifford [40]. Note that faster versions of this algorithm have been found by the same authors in the case where m is proportional to n [42]. Integration of the algorithm of [42] into Perceval is a subject of on-going work. The simulations were performed on a 32-core 3.1GHz Intel Haswell processor, at the rate of 8547 runs (samples) per second. It took roughly 2 days to collect 1.5 × 10 9 samples.
Brute-force certification that our simulation has correctly implemented Boson Sampling would require, for each of the 300 Haar random unitaries, calculation of approximately 10 14 permanents (to get the ideal probability distribution of the Boson Sampler), then performing ≈ 10 14 runs of Boson Sampling using our simulation, in order to get the distribution of the simulated Boson Sampler to within the desired accuracy. Clearly, this task is computationally intractable. In order to get around this issue, while still having some level of confidence that our Boson Sampling simulation has indeed been implemented correctly, we have appealed to computing partial certificates [50,57]. These certificates are usually efficiently computable, but nevertheless can be used to rule out some common adversarial strategies designed to spoof Boson Sampling [58,59], although they cannot provably rule out every possible adversarial spoofing strategy [50,[57][58][59].
The partial certificate we used here is the probability P (K) that all n input photons are measured in the first K output modes of the Boson Sampler [57], where K ≤ m. We believe this choice of partial certificate is natural for benchmarking our simulation mainly for the following reasons. First, it can be computed straightforwardly from the output probabilities of Boson Sampling. Second, it relies on computing high order marginals, which are most likely difficult to compute efficiently classically [28]. Finally, for some values of n, m, and K, P (K) can be computed to very good accuracy by using only a polynomial number of samples from the device [57].
where, for a given m × m unitary U , P (S) is computed as in Equation 7. The average P (K) of P (K) over the Haar measure of m × m unitaries U has been computed analytically in [57]: We computed an estimate P (K) of P (K) by computing an estimateP (K) of P (K) using 5×10 6 samples for each of the 300 Haar random unitaries, then performing a uniform average over all these 300 unitaries. Our calculations for various values of K, together with the standard deviation of the distribution ofP (K) are presented in Table 3. Std. Dev. 0.00005% 0.0014% 0.015% 0.035% 0.045% Table 3: Analytical values ( P (K) , and values computed with Perceval ( P (K) ) of the probability that all 14 photons gather in the first K out of 60 output modes, for various values of K.
We observe a good agreement between the analytical values and the values computed with Perceval , in particular when the value of K is close to 60. This is to be expected, since these values of K are closest to the regime in which n(m−K) m, where it has been shown [57] that P (K) is polynomially close to one in n and m, and therefore a poly(n, m) number of samples from the Boson Sampler is enough to computeP (K) to good precision so as to allow a reliable certification.
As a final remark, note that Perceval can also simulate imperfect Boson Sampling, where the imperfections integrated so far include photon loss and the possibility of multiphoton emissions. For further details we refer the reader to the documentation.

Introduction
Grover's algorithm [2] is an O √ N -runtime quantum algorithm for searching an unstructured database of size N . It is remarkable for providing a polynomial speedup over the best known classical algorithm for this task, which runs in time O(N ). We now provide an intuitive explanation of how this algorithm works, similar to that which can be found in the Qiskit documentation or in [60].
We work in the n-qubit Hilbert space (C 2 ) ⊗n , which allows us to treat a database of size N = 2 n . In its standard formulation, the goal of Grover's algorithm is to look for an unknown computational basis state |T , which we will call the target state, with Let |R be the (normalised) uniform superposition of all the N − 1 computational basis states other than |T , which is orthogonal to |T . The uniform superposition |S of all N computational basis states can then be written as Describing the problem setting in this way allows us to describe Grover's algorithm in a geometric picture. Indeed, we can now think of the states |T , |R , and |S as lying in a 2dimensional plane with a basis {|R , |T }. For example, |S forms an angle θ = arcsin( 1 √ N ) with respect to |R (or π 2 − θ with respect to |T ). Grover's algorithm relies on the application of two unitaries noted U O and U d . The former is an oracle unitary 7 defined on computational basis states |x as while the latter is the diffusion unitary U d is defined as where 1 n the n-qubit identity matrix. On the uniform superposition, the unitary U O acts as Geometrically, U O reflects the state |S through the vector |R . Let's call the resulting reflected state |S 0 . Then, U d performs a reflection of |S 0 through the vector |S . Let the resulting reflected state be denoted |S O,d . One can see that these transformations keep the states in the same 2D-plane, and that |S O,d forms an angle of 2θ with respect to |S (or π 2 − 3θ with respect to |T ). The state after applying U d U O is thus closer to |T than the initial state |S was.
The reflections that we described correspond to one iteration of Grover's algorithm. The full algorithm first creates the state |S by applying n Hadamard gates H ⊗n to an input state |0 ⊗n . It then transforms the state |S into |T by applying k times the unitary U d U O to |S . Following the geometric reasoning above, the required number of iterations k is given by: where · is the ceiling function. As a closing remark, let us mention that the original version of Grover's algorithm presented here is non-deterministic, albeit with a probability of success approaching one exponentially quickly with increasing N [2]. Fully deterministic versions of Grover's algorithm can be obtained, for example by using different oracles and diffusion unitaries than those used here [61], or by using two different types of diffusion unitaries while keeping the same U O [60].

Perceval Implementation 8
Here we reproduce the photonic demonstration of Grover's algorithm of [62]. It uses the spatial and polarisation degrees of freedom to implement a mode realisation of the algorithm with two spatial modes. We spell out the correspondence between the marked database element and the associated quantum state in Table 4.  Quantum Circuit. Figure 7 represents the quantum circuit that implements Grover's algorithm on two qubits, encoded over two spatial modes and their polarisation modes as described in Table 4. It features a state initialisation stage creating a uniform quantum superposition over all states, an oracle stage, and a diffusion stage. Replacing N = 4 in Equation 15 gives k ≈ π 4. π 6 − 1 2 = 1. Thus, one application of the oracle and diffusion gates suffices. Linear Optical Circuit. The linear optical circuit we will use in our simulation is shown in Figure 8. This circuit was experimentally realised by Kwiat et al. [62] using bulk (i.e. free-space) optics and compiled to limit the number of optical elements introduced in the experimental setup, by moving and combining operations like phase shifts. Simulation. Here, we simulate in Perceval the linear optical circuit of [62] (see Figure  8), a photonic realisation of the quantum circuit in Figure 7, which implements Grover's algorithm for marked elements 00, 01, 10 and 11. The results are displayed in Figure 9, along with the experimental results of Kwiat et al. [62]. The simulated results are as expected and thus match the results obtained in the experiment (up to experimental error).

Introduction
The problem of factoring an integer is thought to be in NP-intermediate and the best known classical algorithms only achieve sub-exponential running times. Its classical complexity is well studied since it has been used as the basis of the most widely adopted encryption scheme, Rivest-Shamir-Adleman (RSA) [63], where the secret key consists in two large primes p and q, while their product N = pq is the corresponding public key. In this context, Shor's algorithm [1] greatly boosted interest in quantum algorithms by showing that such composite numbers can in fact be factored in polynomial time on a quantum computer.
From Period-Finding to Factoring (Classically). We start by assuming that there exists a period-finding algorithm for functions over integers given as a black-box. This algorithm is used to find the order of integers a in the ring Z N , i.e. the smallest value r such that a r ≡ 1 (mod N ). In particular, the values of a are sampled at random until the associated r is even. Once such an a has been found, 9 we have that a r − 1 = (a r 2 − 1)(a r 2 + 1) is a multiple of N but not a r 2 − 1 (otherwise the period would be r/2). We test whether a r 2 + 1 is a multiple of N and if so restart the procedure. 10 Otherwise, we have found a multiple of one of the prime factors p, q since a r 2 + 1 divides a multiple of N . Taking the greatest common divisor (GCD) of a r 2 + 1 and N easily yields this prime factor. Indeed, the GCD can be found efficiently classically for example by using Euclid's algorithm.
In summary, all the steps described are basic arithmetic operations, simple to implement on a classical machine, and most of the complexity is hidden in the period-finding algorithm.

Shor's Quantum Period-Finding Algorithm.
The key contribution of Shor's work was to show that the period-finding algorithm can be done in polynomial time on a quantum computer.
Indeed, finding the order for a value a can be reduced to a problem of phase estimation for the unitary U a implementing the Modular Exponentiation Function (MEF) x → a x (mod N ) on computational basis states. It can be shown that all eigenvalues of these operators are of the form e 2iπk r for an integer k and the value of interest r. The phase estimation algorithm uses as basis the Quantum Fourier Transform and its inverse, along with controlled versions of unitaries of the type U a 2 p up to 2 p ≈ N 2 . 11 This part of the circuit is specific to the value of a and N and can therefore be optimised once they have been chosen. Although the most costly part of the algorithm in terms of gates, the unitaries C-U a 2 p (controlled U a 2 p gates) can still be implemented in polynomial time on a quantum computer, making the overall procedure efficient as well.

Perceval Implementation 12
Here we reproduce the photonic realisation of Shor's algorithm from [64].
Quantum Circuit. The quantum circuit shown in Figure 10 for factoring N = 15 using parameter a = 2, whose order is r = 4. It acts on 5 qubits labelled x 0 , x 1 , x 2 (for the top three) and f 1 , f 2 (for the bottom two). The CNOT gates apply a version of the MEF which has been optimised for this specific value of a to the qubits x i in superposition, storing the outcome in qubits f i . The outcome is given by measuring the qubits x 1 , x 2 , followed by classical post-processing.

Linear Optical Circuit.
Since qubit x 0 remains unentangled from the other qubits it can be removed from the optical implementation of this circuit. Furthermore, the CNOT gates on x i , f i can be realised as where the indices denote the qubits to which the gate is applied (in the case of CNOT, the first qubit is the control). Finally, the inverse Quantum Fourier Transform can be performed via classical post-processing, and so does not need to be implemented as a quantum gate in the circuit. The circuit after these simplifications is given in Figure 11. The expected output state of the circuit above is We work with path encoded qubits, as in [64]. With path encoding, each H gate in the quantum circuit is implemented with a beam splitter with reflectivity R = 1/2 between the two paths corresponding to the qubit. In our implementation in Perceval , phase shifters are added to properly tune the phase between each path. CZ gates are implemented with 3 beam splitters with reflectivity R = 2/3 acting on 6 modes: one inner beam splitter creates interference between the two qubits, and two outer beam splitter balance detection probability using auxiliary modes. This optical implementation succesfully yields the output state produced by a CZ gate with probability 1/9; otherwise it creates a dummy state, which can be removed by post-selection.
The circuit implemented in Perceval is illustrated in Figure 12.
The matrix associated to the optical circuit, giving the probability amplitude of each combination of input and output modes, is given in Figure 13.
Simulation. After entering the initial Fock state associated to the input qubit state |0 x 1 |0 x 2 |0 f 1 |1 f 2 , we plot the probability amplitudes of the output state, post-selected on Fock states corresponding to a qubit state with path encoding. The results are given in Figure 14.
After re-normalisation we find that the output amplitudes computed with Perceval match the expected output state described in Equation 17 up to numerical precision.
When decomposing the expected output state in the qubit basis, the qubit states with non-zero amplitude are |0001 , |0100 , |1011 and |1110 for qubits x 1 , x 2 , f 1 , f 2 . We plot the output distribution of the circuit, post-selected on these states, without renormalisation; the result is presented in Table 5.  The distribution obtained with Perceval is uniform over each outcome, which matches the expected distribution in [64].

Outcome Distribution Interpretation.
For each outcome, the values of qubits x 2 , x 1 , x 0 (with x 0 = 0) represent a binary number between 0 and 7, here corresponding to 0, 4, 2, 6 in the order of Table 5. After sampling the circuit, obtaining outcomes 2 or 6 allows to successfully compute the order r = 4 [64]. Obtaining outcome 0 is an expected failure of the quantum circuit, inherent to Shor's algorithm. Outcome 4 is an expected failure as well, as it only allows to compute the trivial factors 1 and 15.
Since the distribution from Figure 5 is uniform the circuit successfully yields a successful outcome with probability 1/2. This probability can be amplified exponentially close to 1 by sampling the circuit multiple times [64].

Introduction
The Variational Quantum Eigensolver (VQE) introduced by [13] is an algorithm for finding eigenvalues of an operator. Applications range from finding ground state energies and properties of atoms and molecules to various combinatorial optimisation problems.
Since the eigenvector |ψ * associated to the smallest eigenvalue of H minimises the Rayleigh-Ritz quotient ψ * |H|ψ * ψ * |ψ * , the eigenvalue problem can be rephrased as a variational problem on this quantity. The VQE algorithm uses as a sub-routine the Quantum Expectation Estimation (QEE) algorithm, developed in the same paper, whose task is precisely to compute the expectation value H := ψ|H|ψ of Hamiltonian H for an input state |ψ (here assumed to be normalised). Given an ansatz represented by tunable experimental parameters {θ i }, set to some initial values {θ init i }, the VQE algorithm works by iterating a loop which first applies the QEE algorithm on a state generated using the initial parameters {θ init i }. Then, based on the outcome, it tunes these parameters using a classical minimisation algorithm. This process is repeated until a termination condition specified by the classical minimisation algorithm is satisfied.
Let n be the number of qubits of our system, σ x , σ y , and σ z the single qubit Pauli X, Y, and Z matrices, and 1 the single qubit identity matrix.

Perceval Implementation 13
Here we use Perceval to reproduce the original photonic implementation of the VQE from [13]. For small enough instances of the problem, Perceval is able to compute explicitly the complete state vector |ψ . This allows us to skip the QEE subroutine and directly compute the mean value H .
Linear Optical Circuit. We use the linear optical circuit of the original paper [13] which was first introduced in [67]. The circuit has 6 optical modes, consists of 13 beam splitters, 8 tunable phase shifters and 4 single-photon detectors, and is shown in Figure  15. This circuit is essentially a 2-qubit circuit where qubits 1 and 2 are path encoded respectively in mode pairs (1, 2) and (3, 4) (here mode numbering is from 0 to 5, from top to bottom). Modes 0 and 5 are auxiliary modes. The circuit consists of two parts. The first part being the 3 central beam splitters in Figure 15, together with the 2 beam splitters to the left and right of these central beam splitters, and which act on modes 3 and 4. The 3 central beam splitters have reflectivity 1/3, whereas the two to the left and right of the central ones have reflectivity 1/2. These 5 beam splitters are used, together with modes 0 and 5, to apply a CNOT gate to qubits 1 and 2 [46]. This CNOT is successful with probability 1/9 under post-selection. The second part of the circuit consists of the remaining 8 beam splitters and phase shifters, and is used to implement arbitrary single qubit rotations on qubits 1 and 2. All 8 phase shifters have tunable phase shifts, and all 8 beam splitters have reflectivity 1/2. Simulation. Our strategy for the simulation of the VQE is as follows. We first compute the output state vector |ψ of the linear optical circuit, which depends on the phase shifters parameters (φ i ) i∈{1,...,8} [13], for some random initial configuration of these parameters. This enables the evaluation of the Rayleigh-Ritz quotient in Equation 18. Then, we proceed to minimise this quantity by using the Nelder-Mead minimisation algorithm [68].
Using these techniques, we computed the ground-state molecular energies for the Hamiltonians of three different experiments [13] [69] [70]. In all these experiments, we used the circuit of Figure 15 to compute the output states |ψ . In line with the VQE algorithm, at each iteration we produce a different |ψ by tuning the angles of the phase shifters of the circuit in Figure 15. To test the validity of our techniques, we computed the theoretical ground-state molecular energies of the Hamiltonians of each of the 3 experiments. This was done with the NumPy [71] linear algebra package. A plot of these energies computed both theoretically and with Perceval is shown in Figure 16.    The circuit of Figure 15 was originally used in [13] to compute the ground-state energies of the Hamiltonian in Figure 16a. We observe a very good overlap between theoretical and computed energies, indicating that the Perceval simulation was succesful. In [69,70], different circuits than that of Figure 15 were used for computing the ground-state energies. Nevertheless, and in order to explore the expressivity [65] of the circuit of Figure 15, we computed and plotted in 16b,16c the results of ground-state energy computations of the Hamitonians in [69,70] using a VQE with the circuit of Figure 15. Our results in Figures  16b,16c indicate in general a good overlap between theoretical and computed energies, with the deviations between theoretical and computed energies probably due to the fact that the circuit of Figure 15 is not expressive enough to give better accuracies.

Introduction
Linear optics has proven to be a fascinating playground for exploring the advantages offered by quantum devices over their classical counterparts. For example, as discussed in Section 4.2, Boson Sampling [28], shows how a quantum device composed only of single photons, linear optical circuits, and single-photon detectors, can perform a computational task which quickly becomes unfeasible for the most powerful classical computers.
In [14], the authors show how linear optical circuits, similar to those used in Boson Sampling, can be used to solve other problems of more practical use than sampling. The basic idea is, following the work of [72,73] for qubit-based quantum circuits, to encode data points x ∈ R onto the angles of phase shifters of a universal linear optical circuit. Effectively, this allows a non-linear manipulation of these data points. Indeed, this encoding of the data points x allows to express the expectation value of some observable, computed using the linear optical circuit, as a Fourier series of the data points x [14]. The Fourier series, being a well known universal (periodic) function approximator, can be used for a variety of tasks, including approximating the solution of differential equations, which we will study here. Interestingly, Ω = {−n, . . . , n} where n is the number of photons inputted into the linear optical circuit [14]; meaning that the expressivity of the Fourier series -how well it can approximate a given function -depends (among other things) on the number of input photons of the linear optical circuit.
In the coming sections, we give an example of a quantum machine learning algorithm, using the above encoding, which solves a differential equation.

Expression of Photonic Quantum Circuit Expectation Values as Fourier Series
We focus on constructing a universal function approximator of a one-dimensional function f (x) with a photonic quantum device of n photons and m modes. This device consists of the following components: n sources of single photons, m-mode universal linear optical circuits, and m number-resolving single-photon detectors [74]. Linear optical circuits can be configured to implement m × m unitary matrices U . A linear optical circuit is called universal if it can implement any such unitary U [35]. The n single photon sources produce input Fock states of the form |n 1 , ..., n m , where n i is the number of photons in mode i, and i=1,..,m n i = n. The Fock space of n photons in m modes is isomorphic to the Hilbert space C M , with M = m+n−1 n [28]. This isomorphism, together with a homomorphism from U(m) to U(M ) (the m and M -dimensional unitary groups) detailed in [28], allows to understand the action of the m × m unitary U implemented by the linear optical circuit as an M × M unitary U acting on the input Fock state. In the rest of this section, as in [14], we will work with these M × M unitary matrices when computing the universal function approximator.
The circuit architecture of [14] implements the M × M unitary transformation The phase shift operator S(x) incorporates the x dependency of the function we wish to approximate. It is sandwiched between two universal linear optical circuits W (1) (θ 1 ) and W (2) (θ 2 ). The parameters (angles of beam splitters and phase shifters of the linear optical circuits) θ 1 and θ 2 are tunable to enable training of the circuit, θ := {θ 1 , θ 2 }. Let |n (i) = |n where the sum is taken over all M possible Fock states |n (f ) of n photons in m modes, and {λ |n (f ) } some tunable set of parameters. The expectation value of M(λ) with respect to the output state U(x, θ)|n (i) of the linear optical circuit can be computed by measuring the output modes of this circuit using number-resolving detectors. f (n) (x, θ, λ) can be rewritten as the following Fourier series [14] where Ω n = −n, n is the frequency spectrum one can reach with n incoming photons and {c ω (θ, λ)} are the Fourier coefficients. Hence this specific architecture can be used as a universal function approximator.

Application to Differential Equation Solving
The most general form of a differential equation verified by a function f (x) is being an operator acting on f (x), its derivatives and x. Given such an expression, we wish to optimise the parameters θ such that f (n) (x, θ, λ) is a good approximation to a solution f (x). More precisely, in this Quantum Machine Learning task, we aim at minimising a loss function L(θ) whose value is related to the closeness of our approximator to a solution. This is done via a classical optimisation of the quantum free parameters θ, yielding ideally in the end For the solving of differential equations, the loss function described in [75] consists of two terms The first term L corresponds to the differential equation which has been discretised over a set of M points {x i }: where L(a, b) is associated to the initial conditions of our desired solution. It is defined as where η is the weight granted to the boundary condition and f 0 is given by f (x 0 ) = f 0 , for some initial data point x 0 . These functions will be applied to the iterated approximations f (n) (x, θ, λ).

Perceval Implementation 14
Our aim is to reproduce some of the results of [75], where the authors use so-called differentiable quantum circuits, together with the classical optimisation BFGS method of SciPy [76], to provide an approximation to the solution of the nonlinear differential equation with λ, κ ∈ R, and boundary condition f (0) = f 0 ∈ R. The analytical solution of this differential equation is Here, we solve this differential equation using the linear optical circuit architectures of Section 4.6.2 simulated in Perceval , together with classical optimisation. Note that the authors of [75] use Chebyshev polynomials as universal function approximators, whereas here we will use the Fourier series. In our Perceval implementation, η is chosen empirically as η = 5, granting sufficient weight to the boundary condition. Concerning the λ parameters, each one of them is sampled uniformly randomly in the interval {−200, . . . , 200}, which is also empirically chosen. One could tune these λ parameters as well for greater accuracy, at the cost of a considerably slower minimisation procedure. The number of modes m is taken to be equal to the number of photons n. Differentiation is numerically conducted, df dx ∆f ∆x , ∆x taken as 10 −3 . The discretised version of the loss function is taken on a grid of 50 points, uniformly spaced between 0 and 1.
Results are shown for various photon numbers in Figure 17a, demonstrating the increase in expressivity of the quantum circuit with growing n. The converged solution for 6 input photons provides a convincing approximation to the analytical solution. Increasing the number of photons further should yield even more accurate results, allowing for the realisation of more complex Quantum Machine Learning tasks.
A common indicator of the performance of a machine learning algorithm is the convergence of its loss function [77]. Confirming the results from Figure 17a  To summarise, we have used Perceval to simulate linear optical circuits providing universal function approximators [14], and shown that these can be used together with techniques from Quantum Machine Learning to accurately compute the solutions of differential equations. The accuracy of these function approximators depends, among other things, on the number of input photons of the linear optical circuits. An interesting future direction we aim to pursue is using Perceval and our developed techniques to solve other types of differential equations of significant practical interest [78,79].

Conclusion
Perceval is a unique framework dedicated to linear optics and photonic quantum computing. This white paper has aimed to provide an overview of the platform, the motivations for its development, its structure and main features, and to give a variety of examples of Perceval in action. These examples are intended to be illustrative of some of the immediate uses of the platform.
Perceval 's simulation back-ends are optimised to run on local desktop devices, with extensions for HPC clusters. They can be used to run computational experiments to finetune algorithms, compare with experimental data from actual experiments and photonic quantum computing platforms, and can reproduce published articles in few lines of code.
Perceval is intended to be accessible to physicists, both experimental and theoretical, and computer scientists alike, with a goal of providing a bridge between these communities. With the intention of keeping a strong connection between software and hardware for photonic quantum computing, a major focus of future development will be on the continued development of realistic noise-models, that can describe with increasing accuracy the functioning of specific hardware components.
Perceval allows users to design algorithms and linear optical circuits through a large collection of predefined components. The collection of algorithms described here are available and presented as tutorials in the documentation. This is an open source project, with a modular architecture, and is welcoming of contributions from the community.
It is intended that future versions of Perceval will include more optimised simula-tors, noisy simulators, features for working with density matrices and cluster states, more advanced features and options for detectors to cover both threshold and photon-number resolving detectors, as well as features for treating circuits with feedforward.