BosonSampling.jl: A Julia package for quantum multi-photon interferometry

We present a free open source package for high performance simulation and numerical investigation of boson samplers and, more generally, multi-photon interferometry. Our package is written in Julia, allowing C-like performance with easy notations and fast, high-level coding. Underlying building blocks can easily be modified without complicated low-level language modifications. We present a great variety of routines for tasks related to boson sampling, such as statistical tools, optimization methods, classical samplers and validation tools.


Introduction
Quantum computing holds strong hopes for a profound change of computational paradigm and power in the coming years [1,2,3].A large number of experimental and theoretical efforts took place in the last decades, leading to the fast growth in experiments' complexity.While a general purpose, universal quantum computer seems so far out of reach, non-universal models already claim a quantum advantage for specialized tasks on a quantum system that cannot be simulated in a reasonable amount of time with the most advanced classical super-computers [4,5].One of these restricted models is that of boson sampling.In their seminal paper [6], Aaronson and Arkhipov describe a quantum version of Galton's board where n bosons -generally photons -are sent through a m-modes linear interferometer implementing a unitary transformation U .The task consists in sampling output mode patterns of particles, such as finding 2 photons in the first mode, none in the second, etc. Depending on the importance of the noise sources from the experimental components, the output distribution D can be hard or easy to sample from 1 .By hard we understand that it takes a super-polynomial number of steps to generate a sample from D on a classical computer.Indeed, for modest experimental noise, AA prove that this task remains hard, under widely believed conjectures in complexity theory.However, when sufficiently strong noise is present, e.g.due to partial distinguishability or particle loss, then classical algorithms can sample from D efficiently [7,8,9,10,11,12,13].
Boson sampling sparked a great interest both from theoreticians and experimentalists.Various alternative schemes were constructed, such as Scattershot and Gaussian boson sampling (GBS) [14,15] to alleviate the technical difficulties in generating single photons.These efforts culminated in multiple claims of quantum advantage in GBS [16,17,18,19], therefore questioning the validity of the extended Church-Turing thesis [20].Standard boson sampling saw experimental implementations with n = 20 photons in m = 60 modes [21].Other experimental platforms than photonics were also considered [22,23].
These factors imply a great need for optimized, high performance numerical simulations of multiphoton interferometry.As the fight for quantum advantage is played between quantum devices and classical computers, scalable, versatile and fast implementations of the boson sampling algorithms become even more important.In particular, an adaptive framework for the boson sampling community is of great interest given the speed of growth of this field of research.Indeed, many new boson sampling paradigms were introduced theoretically, and a great variety of experimental techniques were refined, enhancing the demand for evolutive code bases.
In this paper we expose our package, Boson-Sampling.jl,written in the Julia programming language to tackle classical tasks related to boson sampling.It intends to solve the challenges discussed above by providing an adaptable and extendable framework for multi-photon interferometry designed to study new boson sampling paradigms.The choice of Julia solves the twolanguage problem: instead of using a high-level wrapper (e.g.Python) with time-critical functions written in a low-level languages (e.g.C, Fortran), Julia is fast out-of-the-box while being easy to write.As a byproduct, all functions are fast, and not only time-critical ones.This is especially relevant for time-crunched researchers who are not professional-programmers but still require a fast execution time for this type of computationally intensive research.The package is evolutive and obeys the Open Source Software standards.It welcomes contributions from the community through its GitHub portal [24].A complete documentation as well as a pedagogical tutorial are provided [25].
The paper is structured as follows.We first explain how the specific choice of the Julia programming language is, in the eyes of the authors, optimal for both theoreticians and experimental-ists and how it affects the structure of our package and benefits users.In the second section of this paper, we explore the capabilities of our codebase and expose how our package is structured, including more specific descriptions and examples in the third section.It contains worked-through examples of well-known physical phenomena such as the Hong-Ou-Mandel effect.We also show how one can make use of the package architecture to take into account new sources of noise for instance.Finally, we compare our package to other existing software and discuss possible future features.

Julia
Julia is a high-level and dynamic programming language originally designed for technical computing, in particular, the simulation of physical systems.It is catching more and more attention from the scientific community and has been adopted in several high class projects such as working with dynamical systems [26], satellites simulations [27], the Celeste project (1.54 petaFLOPS) [28], the Climate Modelling Alliance [29], as well as NASA who saw a 15.000 times speed improvement over their previous Simulink/MATLAB code [30].The field of quantum information and computation is already well covered, with several packages such as Yao.jl and Quantum Optics.jlrivalling in scope and efficiency with the leaders of the market [31,32,33,34].
The main strength of Julia, in particular for scientists, is its goal to solve the "two-languages problem".It is very common that code is first prototyped in a high-level, easy to write, language such as Python and then made fast to execute in low-level, and hard to write languages such as C++.Julia allows the coding to be convenient and efficient, while being as fast as lower level languages (such as C and Fortran) out of the box.This performance comes from its just-intime (JIT) compilation and specific type architecture, more akin to functional programming than Object Oriented Programming (OOP) found in most languages (C++, Java, Python,...).This bottom-up approach is especially suited for scientists, whose code evolve organically as results are found, without a set roadmap ("proofs then theorems") as would be preferred for OOP programming [35].
Another strength of Julia's type structure is multiple dispatch: functions are written without strong restrictions to the type of their arguments, and can be reused with custom types with no modifications [36] (provided that the function makes sense for these new types).For instance, a function computing the square root would be restricted to specific types (say, Float64) in strongly typed languages.A new type, say BigFloat, would require a rewriting of the function, while in Julia the same square root algorithm will work automatically.This design philosophy, is, once again, a blessing for scientists who may not know in advance what the future of their code looks like.
2 Contents and architecture

Package overview
Our package, BosonSampling.jl together with its paired package Permanents.jl,provides tools to classically simulate multi-photon interferometry experiments.This includes the wellknown cases of standard boson sampling, scattershot and Gaussian boson sampling.
Our platform is aimed at high-performance computations for scientists.This is reflected in the choice of the language, Julia, and the modular architecture.Our goal is to make it both easy and fast to implement new models while maintaining all the provided methods and tools adapted to the new features with fast execution time.
The package is evolutive, and welcomes outside contributions from the community.It is intended to be a toolkit for both experimentalists and theoreticians, allowing to easily bridge between new theoretical models and their realistic implementation.
We provide state-of-the-art tools addressing relevant boson sampling problems.Let us now give an overview of these tools made available.
1. Boson-samplers, which allow to get samples as if performing an experiment.These samplers are based on the best known algorithms from the literature.Those samplers, include experimental noise such as photons' partial distinguishability and loss.This noise is critical as it will affect the classical hardness of simulating a given experiment.
2. Event probability computation routines, allow to observe specific input/output patterns.Special interest is given to noise, and in particular to a fast implementation in the case of partial distinguishability.
3. Bunching tools, aimed at investigating the tendency of bosons to bunch and of fermions to anti-bunch, as exemplified by the Hong-Ou-Mandel effect.We provide functions to investigate the link between bunching and matrix permanent conjectures [37,38].

4.
Validation tools for boson-samplers.These allow, given experimental as well as blackbox samples, to discuss the degree of multiphoton interference observed in a boson sampler device, and to estimate the level of noise.This gives indications as to whether the device outputs samples that are classically hard to obtain, i.e. if quantum advantage can be claimed.Standard tools from the literature are unified in a recent validation protocol [39] along with the common method of correlators [40,41], suppression laws [42,43,44,45], and full bunching [38,46] and also allowing to recover marginal probabilities [47].
5. Optical circuits built from optical elements.While in the standard boson sampling problem, networks, described by a unitary matrix U , are chosen from the Haar measure, specific interferometers are also implemented such as the Fourier and Hadamard transformations.In addition, a suite of linear-optical elements is provided, so that networks can be constructed as in an experimental setting.
6. Optimization routines to maximize cost functions over unitary matrices.Maximizing certain properties of interferometers, such as photonic bunching, can be convenient for the design of experiments.We include algorithms from Ref [48] that can optimize costfunctions on the Haar measure.
7. Time-bin interferometry types and circuits to simulate the experimental boson sampling proposal of [49], where photons are separated in time bins and sent through a timedependant beam-splitter and delay lines.

Code architecture
The two main procedures achievable with this package are the simulation of photonic experiments and validation for multi-photon quantum experiments.We will first describe the simulation architecture and then the validation procedure.

Simulation
The present package takes advantage of Julia's type system through a hierarchy used to classify the building blocks of a quantum interferometry experiment, namely: a photonic Input state sent through an Interferometer and measured according to an OutputMeasurement.This last category also includes (classical) sampling.Those three elements are represented by abstract types, that is, they cannot be instantiated and form the backbone of the package structure.They can be seen as containers of sets of related concrete types, which can be instantiated.Concrete types will represent the specific building blocks of a given interferometric setup.For instance, an interferometer defined from a Haar distributed unitary, RandHaar, is a concrete type of the abstract type Interferometer as displayed in Fig. 2. The latter structure is at the core of the code efficiency and tunability.Indeed, the generic programming allowed by Julia combined to the inherent type hierarchy makes the package easy to adapt depending on the user's needs.As a new model can simply be embedded as a concrete type of an existing abstract type, implementing new elements can be done easily without modifying the overall package structure, or any function that acts in the same conceptual fashion on the new type.For instance, a function adding noise could work similarly on any type of Interferometer, be it RandHaar or Hadamard.We will now describe the three main abtract building blocks used in simulating an experiment.

Input
When considering photons in Fock states at the input, BosonSampling.jl offers three families of input types depending on the partial distinguishability of the particles.We refer to indistinguishable photons (identical arrival time, polarization) using Bosonic.When they are partially distinguishable, they fall under PartDist.Finally, we use Distinguishable for completely distinguishable input photons.The type PartDist allows for a general description of the partial distinguishability.Common models used in the literature are also implemented as subtypes [58].
The input state can also be Gaussian.In order to take into account partial distinguishability, we adopt the model introduced in [59].More precisely, any Gaussian input will be model of fully indistinguishable states while, partial distinguishability can be introduced at the level of the interferometer.The input state is decomposed into interfering and non-interfering modes which evolve independently through the interferometer.

Interferometers
The second building block of a quantum multi-photon experiment is the interferometer.A common practice when simulating boson sampling is sample a unitary matrix according to the Haar random measure.Specific interferometers, such as the discrete Fourier and the Hadamard transforms are built-in, together with elementary linear optical circuit elements such as beamsplitters or phase shifts.This allows the user to build networks from scratch.Generic interferometers are supported.Loss can be accounted for in generality.

Measurement/Sample
This object represents the output of the experiment, which we classify between a measurement that can be the probability of a specific output pattern, property or a randomly generated sample.For Fock states, a possibility is to observe a given output mode occupation s = (s 1 , s 2 , ...), where 0 ≤ s i ≤ n is the number of photons in mode i.Likewise, we can observe a specific partition photon count k = (k 1 , k 2 , ...), where 0 ≤ k i ≤ n is the number photons in a bin K i gathering multiple output modes.
Different samplers are available depending on the level of distinguishability of the input particles or the use of a lossy interferometer.For exact sampling, we provide Clifford-Clifford's algorithm [60].For imperfect sampling, we rely on the works of Moylett et al. [61].Their detailed implementations are discussed in the package documentation.

Validation
Experimental data can be analysed within this package to validate boson sampling devices.The procedure is straightforward for the user: data is held in form of Events, comprising all known experimental parameters, such as the detector outputs.Hypotheses are provided, such as that the input is made of indistinguishable particles.An alternative hypothesis could be that they are distinguishable.Another task that can be addressed with this package is to estimate the amount of noise, e.g due to partial distinguishability.
Common validation protocols are implemented, with a comprehensive statistical analysis, providing p-values, etc.A detailed de-scription of these protocols can be found in [62,39,13,63].

Examples
This section presents some basic experiments that can be simulated by BosonSampling.jl.Further details about the package, as well as a complete, step-by-step tutorial, can be found in the documentation, hosted on the GitHub page of the package.As for any registered Julia package, BosonSampling.jl can be installed by simply executing the following commands in the Julia REPL: using Pkg Pkg.add("BosonSampling") Code sample 1: Installing BosonSampling.jl.

Hong-Ou-Mandel
It is well known since the experimental work of Hong, Ou and Mandel [64] that two indistinguishable photons impinging on the two input modes of a balanced beamsplitter will bunch in one of the two output modes, leaving a zero probability to observe one photon at each output mode (the coincidence probability).As an introductory example, we review the effect of partial distinguishability in the HOM effect.In this context, it is a common practice to use the time delay ∆τ between the two incoming beams as a source of partial distinguishability, while the distinguishability parameter itself is ∆ω∆τ with ∆ω the uncertainty of the frequency distribution.In order to make the parallel with the built-in interpolating model, we substitute its linear parameter x by e −∆ω 2 ∆τ 2 .In this way, a distinguishable input is recovered for ∆ω∆τ → ∞ and a bosonic input for ∆τ = 0.
We reproduce in the code sample 2 the HOM experiment as described in Fig. 3(a), with a variable overlap between the two initial wave functions (Fig. 3(b)) and compute the coincidence probability.

Sampling
As described in Sec.2.1, several classical algorithms for boson sampling are implemented, such as the Clifford and Clifford's algorithm [60] for exact sampling and algorithms taking into account typical sources of noise when considering realistic schemes [61].To simulate a boson sampling experiment, one defines an input, the interferometer and a measurement.As in the previous example, an Event type is used as a container to hold all interferometric parameters.Given n photons and m modes, the package proposes methods to compute the entire photoncounting distribution exactly or approximately when different sources of noise are included [12] through the noisy_distribution function.It returns by default the exact distribution, an approximated distribution in which the computation of the probabilities are truncated and a distribution reconstructed from a Metropolis Independent Sampler (MIS).Below we compute those three distributions for n = 3 photons in m = 5 modes with a distinguishability parameter x = 0.74 and a loss l = 0.63 (see Fig. 4).

Photon counting in partitions
One of the novel theoretical tools implemented in this package is the ability to efficiently compute photon counting probabilities in partitions of the output modes of the interferometer.A detailed theoretical and numerical investigation can be found in Ref. [39].
Let us now show how we can obtain the probability of finding k = 0, ..., n photons in a subset consisting of the first half of the modes of a Haarrandomly chosen an interferometer with n = 25 photons and m = 400 modes.With these parameters, a brute force calculation would be pratically impossible, while the new methods used in this package allow for a fast computation.

Validation
We now show how experimental data can be validated.In this example, we compare the null hypothesis H 0 = (bosonic input) to the alternative H a = (distinguishable input).Multiple validation protocols are implemented in this package, as shown for instance in Fig. 2. We choose to use the new formalism developed by some of the authors of photon counting in binned output modes.
In Fig. 1, we display an example of validation protocol.Data is assessed using a Bayesian procedure on the event probabilities (instead of partition photon counting in the code sample).As explained in detail in [63,39], each sample updates the confidence level ξ.To be satisfied that the null hypothesis is true, we could require ξ = 95% for instance.We plot this curve averaged over n trials = 500 with n samples = 300 events for each run, in a randomly chosen interferometer with m = 30 modes and n = 10 photons.The input photons are set to be indistinguishable.The heatmap shows the density of the ξ-curves, on a logarithmic scale.

Incorporating new detectors: dark counts
We now show one of the main advantages of this package and Julia: the ability to create new models while keeping most of the structure untouched.
Consider a simple model of threshold detector.With a probability p the detector clicks even though no photon was present, so-called dark count.
We start by defining a new type of output measurement that we label DarkCountFockSample, which is a sub-type of OutputMeasurement.It takes as arguments the probability p to observe a dark count.A variable holds the sampled ModeOccupation, but is first initialised empty at creation of DarkCountFockSample.It will be filled later.This is an interesting feature of Julia: the absence of an argument is integrated efficiently by the compiler while a typical, variablesize array would not.The second step is to define a new method sample! that implements a sampling algorithm for this specific type of detector.
As mentioned before, Julia enables multiple dispatch, which means that the new method will not overwrite the already existing ones but will extend the method to the new detectors.That is, when calling sample!, Julia will automatically apply the relevant method for a given type (FockSample, DarkCountFockSample,...).The new algorithm first extracts a sample without dark counts, then adds their influence on the photon counts:  Code sample 8: Usage of the newly added detectors.

Julia is fast!
When simulating a boson sampling experiment, the most time consuming part comes from the computation of the transition probabilities, since they are related to the evaluation of the permanent.The best known exact algorithm for its computation is due to Ryser [65] and runs in O 2 n−1 n time for a matrix of size n (using Gray ordering).Having an intensive usage of Ryser's algorithm, we compare the running time of its implementation in Julia with two other interpreted languages: Matlab and Python.Fig. 5 shows how Julia outperforms the two other implementations by two orders of magnitude in computation time.

Comparison to existing softwares
The present package is designed for multi-photon interference and boson sampling experiments, taking into account sources of noise and providing a wide variety of validation tools while remaining highly flexible to new paradigms.Other packages for the numerical evaluations of matrix functions and simulation bosonic systems were published in the last few years.Among the most popular softwares are The Walrus [66] and Perceval [67], both written in Python with a C++ backend.
We illustrate how the Julia programming language is efficient and how, indeed, it can reach Clike performance while keeping the coding as easy, fast and convenient as in Python.We benchmark common functionalities in the three softwares.We first compute permanents of Haar distributed unitary matrices2 (Fig. 6.a) showing that our implementation of Ryser's algorithm remains faster even on single-core.Our simulations of perfect boson sampling through Clifford and Clifford's algorithm achieves comparable performance to what is possible on Perceval, and is more performant only via multithreading (Fig.  In the final stages of writing this paper, SO-QCS [68] was released, a C++ library centered around performance at cost of handiness and ease of modification.

Development path
This package is a continuously expanding, and many more features are expected to be implemented in the coming years.Current development areas are Gaussian boson sampling related functionalities, in particular simulating GBS with partial distinguishability in a lossy Many of the functions used in this package will be made faster, such as computing permanents and hafnians via GPU computing, and HPC-friendly deployment through multithreading and distributed computing.
Potential users are welcome to contact the authors for more details or specific development paths, either through the emails listed above or directly through the Github interface.
Outside contributions are welcome, and multiple improvement paths are outlined on the Github discussion page.

Conclusion
In this paper we introduced a new, free open source Julia package for the high-performance simulation of multi-photon interference.We motivated how the choice of language is particularly relevant for scientists, and allows for a package structure that is easy to modify while keeping execution time as fast as in low-level languages.This package is therefore aimed at experimentalists, who need to implement specifics of their hardware and at theoreticians, who want to develop new models and boson sampling paradigms.We showed how the performance positively compares to similar packages written in other languages.

Figure 1 :
Figure 1: Validating boson sampling using Bayesian methods (see Sec. 3.4).Is displayed the confidence ξ that the input is made of 10 indistinguishable photons, averaged over n trials = 500, as more samples n s are provided to the validation protocol.The color scale represents the density of ξ curves.White curves are examples of validation runs.

Figure 2 :
Figure 2: Architecture of BosonSampling.jl(a) Simulation framework: the abstract types representing the experimental setup are displayed as the colored boxes, containing they related concrete types.(b) Validation framework.An extensive description of every type and their usage can be found in the documentation and the associated tutorials.

Figure 3 :
Figure 3: The Hong-Ou-Mandel effect: (a) experimental setup (b) time delay between the incoming wave packets (c) HOM dip translating the bosonic interference appears when the overlap between the beams is perfect.

#
Set the model of partial distinguishability T = OneParameterInterpolation # Define an input of 3 photons placed # among 5 modes x = 0.74 i = Input{T}(first_modes(3,5), x) # Interferometer l = 0.63 U = RandHaar(i.m)# Compute the full output statistics p_exact, p_truncated, p_sampled = noisy_distribution(input=i,loss=l,interf=U) Code sample 4: Computing the output mode occupation statistics with partial distinguishability and loss for 3 photons at the input of a (5,5) Haar distributed unitary.Both computing p_truncated and p_sampled have by default a probability of failure and a maximal error equal to 10 −4 .It is still possible to refine their computations by setting the keywords error and failure_probability when calling noisy_distribution.

Figure 4 :
Figure 4: Photon-counting statistics at the output of a 5-mode Haar distributed unitary matrix.(a) Comparison between the ideal distribution and the approximated one.(b) Comparison between the theoretical distribution and a sampled distribution from 10 5 samples.

Figure 5 :
Figure 5: Benchmarks for the single-core running time of Ryser's algorithm in Julia, Matlab and Python on a 3.2GHz Apple M1 processor. 6.b).

Figure 6 :
Figure 6: (a) Benchmarks measuring the average time elapsed when computing the permanent of matrices with increasing size n.Are displayed the running times of (orange curve) The Walrus, version 0.19, on a single-core, (blue curve) Ryser implementation of Permanent.jl,version 0.1.1,single-core, (green dots) Glynn implementation of Perceval, version 0.6.2(Quandelibc version 0.5.3), on a single-core, (purple dots) Ryser implementation of Perceval, 4 threads.(b) Average time elapsed when simulating perfectly indistinguishable bosons via Clifford and Clifford's algorithm in Julia on 1, 4 and 8 threads compared to Perceval's implementation.Both benchmarks are realized with the same configuration as in Fig. 5.