SOQCS: A Stochastic Optical Quantum Circuit Simulator

We present Stochastic Optical Quantum Circuit Simulator (SOQCS) C++/Python library for the simulation of quantum optical circuits, and we provide its implementation details. SOQCS offers a framework to define, simulate and study quantum linear optical circuits in the presence of various imperfections. These come from partial distinguishability of photons, lossy propagation media, unbalanced beamsplitters and non-ideal emitters and detectors for example. SOQCS is developed as a series of different modules which provide quantum circuits, different simulator cores and tools to analyze the output. Quantum circuits can be defined from basic components, including emitters, linear optical elements, delays and detectors. Post-selection can be configured straightforwardly as part of detector definitions. An important attribute of SOQCS is its modularity which allows for its further development in the future.


Introduction
Quantum computing promises a considerable speed up of certain computational problems [1,2].It has been shown that universal quantum computation can be attained using optical circuits with linear elements and post-selection [3,4,5].A state of a quantum bit or qubit can be easily encoded using photonic quantum states.The initial photonic state is created by means of emitters and transformed by the optical circuit into the final state that is measured using detectors.A quantum optical circuit capable of complex operations is defined in this context as a network of simpler linear optical elements, such as phase shifters and beamsplitters.The schematic of a general optical circuit is illustrated in fig. 1.
There has been a considerable development of software platforms for quantum computing in the last few years.For example, IBM introduced the open source environment Qiskit [6].Also, tools such as Openquasm [7], Cirq [8], Quil [9] and QP [10] concern the simulation of the logical operations of quantum computers.Additionally, there are platforms like D-Wave development kit [11] that are more oriented to interface with specific hardware.On the other hand, software focused on simulation and interface with optical platforms is primarily centered on continuous variable models [12,13].
Figure 1: Schematic of a quantum optical circuit.In the context of this article an optical circuit is defined as a set of interconnected linear optical elements, complemented by additional operations such as measurement and delays.
In this paper we present the implementation details of the C++ and Python library SOQCS [14] for the simulation of quantum optical circuits where photonic states are represented using a Fock state description.Different methods or cores are implemented in SOQCS to calculate the probability amplitudes and other characteristics at the circuit output.
In the time span of SOQCS development other similarly focused software packages have arisen reflecting the interest in these kinds of simulations.For example, these include Perceval [15] and PhotoniQLAB [16] libraries.Compared with those platforms SOQCS has a strong focus in simulating the effect of various circuit imperfections while still being friendly to users with limited programming expertise or physical knowledge.Previous libraries either do not support imperfections or they do it with less depth than SOQCS in a manner that requires an additional coding effort.
The paper is organized as follows.In section 1 we will introduce some basic concepts and the overall architecture of the SOQCS library.In sections 2 and 3, we will explain the basic data structures and how they are used to perform a simulation.In section 4, different simulation methods available in SOQCS are reviewed.In section 5, we explain some of the main physical imperfections implemented in SOQCS, such as partial photon distinguishability, presence of mixed photonic states in the emission process and photon losses.Finally, in section 6 we will discuss how quantum measurement is treated in SOQCS.This includes the physical effects in the detection process and provides also details on how post-selection is defined and what data structures are needed to analyze the output.The last section concludes the paper with a brief summary and outlook.SOQCS is built in layers on top of Eigen3 [17] library for linear algebra, with the level of abstraction increasing upwards.Each layer contains C++ classes that use other more basic classes defined on lower layers.

SOQCS Python
SOQCS library is built with the objective to simulate general quantum optical circuits where imperfections are present.We define a model for an optical circuit from its interconnected components characterized by a set of suitable parameters.These can be configured to fit particular simulation objectives that may include imperfections.We are particularly interested in a modular approach on how circuits are built and configured by a library user.Moreover we are also interested in the modularity of the library architecture itself, so simulations can be easily adapted to different objectives.
Upon configuration, a user may choose from a catalog of devices and imperfections.The library is built in a way that allows these catalogs to grow and be integrated with the rest of the library in the easiest way possible.This modular philosophy is behind many of the design decisions of SOQCS and will be revisited throughout the paper.
For this reason, SOQCS has been conceived as a stack of layers built on top of the Eigen3 C++ template library for linear algebra and matrix manipulation [17] (see fig. 2).Each layer contains definitions of C++ classes that use other classes defined on lower layers while at the same time they are used by classes on higher layers.At the bottom of the library stack the custom made util (utilities) sub-library resides which contains mathematical and dictionary management tools to perform multiple basic operations all across the software stack.
In the next layer, above util, the C++ classes state and qocircuit are found.These classes represent a photonic quantum state and a quantum optical circuit respectively.In general, a state will be represented as a list of photon occupations by mode.Those modes represent the channel, polarization and wavepacket that a photon may take in the circuit.On the other hand, the circuit class contains the matrix definition of the linear optical circuit and a dictionary that relates mode numbers with their meaning.Additionally, the circuit object contains methods to define its individual elements.Linear optical elements are represented internally by matrices while non-linear elements like delays or detectors are defined by virtual circuit elements that inform and drive the simulator actions (more on this in the next sections).
Furthermore, qocircuit also handles the definition of wavepackets (equivalent to single particle wavefunctions).In SOQCS those wavepackets are limited to Gaussian or Exponential shapes parametrized by real numbers.There is an infinite number of those wavepackets in principle.However, only a discrete set is used in each simulation to make the problem manageable by a computer with finite memory.
The next layer contains the C++ class qodev (or "quantum device").The different C++ classes that describe states, circuits and packets are integrated into a single class to allow for a more transparent and easy to use configuration of a simulation.Sometimes we will refer to this layer as the abstraction layer because it creates an interface to define circuits, packets and detectors in a more human friendly form than some of the classes it encompasses.
On top of qodev are defined the classes dmatrix and p_bin to handle the representation of mixed states and measurement outcomes in the form of density matrices and probability distributions respectively.Finally, the class simulator defines the simulator that constitutes the core of the library.This class uses all the different data structures and methods defined below to transform an input state into an output one and to handle this result in an automated manner according to a set configuration.
At the top of the library we find the layer with the class mthread that allows the execution of various simulations in parallel to speed up for example the calculations of density matrices representing the output of a circuit.Finally, at the highest level, an interface that encapsulates the C++ library is located allowing a Python interpreter to call SOQCS as a Python library.This functionality has the double advantage of making SO-QCS accessible to both C++ and Python users and also use a more efficient C++ implementation in the more widespread and flexible Python environment.

Data structures
The simulator transforms photonic input states of a circuit | ⃗ Ψ in ⟩ into output states | ⃗ Ψ out ⟩.We restrict ourselves to optical circuits made of linear optical elements that can be represented by matrices [3,4,5].In the next sub-sections we discuss how states and linear elements are represented in SOQCS data structures.Additionally, we will show how linear elements are integrated into optical circuits.

States
The photonic states are defined as a linear superposition of kets where each ket in the basis describes a different configuration of the photon occupations by mode.That is, where n i, j ≥ 0 is the occupation of mode j in the ket i and α i is the relevant complex coefficient.A mode j labels a combination of quantum numbers for the channel, polarization and other degrees of freedom.A state can be represented in a computer as a list of kets where each ket consists of a 2-tuple made of a complex amplitude and a list of photon occupations by mode, The simulator C++ class transforms the input state into the output state according to a circuit definition, where Note that the implementation of this transformation in the simulator class is oblivious of the meaning of each mode.We need to maintain a dictionary to keep track of the relationship between the labelling of the different modes and their meaning (channel, polarization, etc).This dictionary is stored as part of the circuit class, so it can be used by all states related to that circuit.The simulator itself only requires the input state and the matrix defining the circuit transformation to perform the calculation.The circuit matrix is also stored as part of the circuit data structure.

Circuits 3.2.1. Basic concepts
In SOQCS a quantum optical circuit is a C++ class.It contains both the dictionary that relates each mode index with its meaning (channel, polarization and wavepacket) and the transformation matrix [3,4,5] that relates the input modes to the output ones.
For example, if we consider a circuit made of a single beamsplitter parametrized by two angles (θ and ϕ), then the corresponding transformation matrix is, This matrix is encoding the relationship between input and output single photon creation operators, This transformation can be used to calculate an output state of this simple circuit from any given input state, For two photons inputs the corresponding outputs of this beamsplitter are, To obtain an output state from an input state given any general circuit it is necessary to: • Describe the circuit providing a list of individual elements and their interconnections.
• Build the circuit transformation matrix from that description.
• Apply a method to perform the transformation operation using the circuit matrix (and the input state).
The two first points will be covered in the rest of this section while the last point is the purpose of the simulator class and it will be treated in the next section.Naturally, a transformation matrix for each individual optical element has to be available in the library.The transformation matrices of many different linear optical elements can be found in the literature [3,4,5].The list of these elements implemented in SOQCS is provided in appendix A. Note that some of them, like detectors or delays perform non-linear operations that are treated in a different way than the linear elements explained here.We will discuss those below.

Circuit matrix construction
The circuit transformation requires a matrix U that maps input into output creation operators for a particular optical circuit.The operation of building the matrix U from an enumeration of individual optical elements is divided in two steps: embedding and composition.
• Embedding: An individual optical element only connects a subset of the available modes in the full circuit.Here we must identify the relationship between the local definition of an optical element and the degrees of freedom of the entire circuit.For example, we may recall that the matrix representation of a beamsplitter is given as where each row and column refers to the channels A and B. If we assume a four channel circuit (without additional considerations about polarization or other degrees of freedom) then we have to specify to which channels the beamsplitter is connected.If A = 2 and B = 3 then, Embedding consists of replacing the appropriate submatrix of the identity matrix by the matrix representation of the optical element.
• Composition: The total circuit matrix U is the product ..U n are the matrices of each of the individual optical elements.Note that, here we are neglecting effects of backscattering.Additionally, we assume short photon packets in relation with the distance between optical elements thus also neglecting resonance effects due standing waves.
As an example, the circuit depicted in fig 3 can be described in SOQCS by the set of instructions, circuit.beasmplitter(0,1,theta1,phi1);circuit.beasmplitter(2,3,theta2,phi2);circuit.beasmplitter(1,2,theta3,phi3); that will translate into the matrix U = U 12 U 23 U 01 where U i j is a beamsplitter matrix for channels i and j.In this set of instructions we list the components of the circuit, their connections but also the parameters of each of the beamsplitters.Overall, SOQCS may represent any optical linear circuit just by listing its elements (with their corresponding parameters) and the connections between them.The construction of the entire circuit matrix is performed in an automated manner.
To perform the calculation it is needed to create and run a simulator C++ object: simulator=soqcs.simulator()outputst=simulator.run_st(inputst,circuit) We will cover the inner workings of the simulator class in the next section.Finally, the output is post selected by the condition |1, 0⟩ in channels one and two.This condition is specified above as a part of the detectors definition in the circuit,

outputst=circuit.apply_condition(outputst)
The final state is printed using the command, outputst.prnt_state(column=1)which provides the intended output state with the phase flip operation applied to the term with two photons, | 0 >: 0.49999999 + 0.00000000 j | 1 >: 0.50000001 -0.00000000 j | 2 >: -0.50000000 + 0.00000000 j Note that the state is intentionally not normalized to unity in order to provide information about the probability of success of post-selection which is P = 0.25.This is not the only way to perform post-selection in SOQCS.Any projector can be defined in a manner similar to definition of states in order to perform the post-selection operation.
In appendices B and C, examples for the CZ [5] and CNOT [18] gates are provided.These examples present different ways to perform simulations and also show how SOQCS can be used to encode qubits into photonic states.

Simulator core
As explained above, the simulator C++ class is used to transform input states into output states.The eq. 3 that describes this transformation can be written in algorithmic form, ψ j } end for end for where ⃗ Ψ in , ⃗ Ψ out , and ⃗ Ψ aux are the objects representing the input, output and auxiliary states as list of kets in the memory of a computer.Note that zero amplitude kets do not need to be stored.
The transformation of the input state is made ket by ket and each ket transformation T ( ⃗ ϕ i , U) depends on the circuit matrix U. Note that, the sum of two states "add" is equivalent to the concatenation of two lists of kets provided that no ket is repeated.In the case of a repetition the amplitudes of those duplicated are summed.We store a dictionary as part of the state data structure to check for repetitions in a straightforward manner.
The function T ( ⃗ ϕ i , U) defines the core operation of the simulator and can be implemented using different methods.SOQCS implements several of them and provides an interface that allows for new ones to be added in the future.The algorithms implementing those methods have a common structure, where a probability amplitude β j is calculated for each possible output ket ⃗ ψ j (given an input one ⃗ ϕ i and the circuit matrix U).The transformation f is different to each chosen method.The output kets are taken from an ensemble B made of an orthonormal subset of kets relevant for the current simulation.At this moment in SOQCS there are implemented two methods to perform this operation and these methods may be used in three different ways depending on B.
• B is made of kets with the same number of photons as the input ket.In this case, the full output distribution is calculated.
• B also assumes photon number conservation but is restricted to include only kets with a maximum of one photon per mode.The overall computation time with respect of the previous option is improved because we are restricting ourselves to a smaller Hilbert subspace.This is useful for circuits implementing qubit operations because the usual encodings of qubits into photonic states do not involve more than one photon by channel.
• B is chosen specifically by the user.Only outcomes that are relevant for a particular simulation are calculated.Those that are considered unimportant are neglected.For example, in the calculation of the fidelity of a particular gate with respect to some parameters only those kets involved in the ideal solution are of interest.Restricting the calculation to even a smaller subset of the Hilbert space makes the overall algorithm faster compared to the previous options.
Next, we review the algorithms which implement T ( ⃗ ϕ i , U) and which are part of SOQCS.

Direct method
The most straightforward way to perform the single ket transformation T ( ⃗ ϕ i , U) is to imitate how the analytical calculation is performed.First, it is convenient to re-interpret the circuit matrix transformation as just a list of rules of the form, where a † i is the creation operator of a photon in the mode i and U j,i is the probability amplitude of where l is the number of modes, ⃗ p is and index list ⃗ p = {p 1 , p 2 , ..., p n } and n = n 1 +n 2 +...+n l is the total number of photons.Thus, the simulation consists of generating all the possible sequences of indexes ⃗ p from 1, 1, ..1 to l, l, ...l that belong to B. Note that different sequences of indexes may lead to the same photon occupation by mode.This straightforward method scales quite poorly.To calculate a single probability amplitude of the output it is needed to consider n!/n 1 !n 2 !..n n !sequences of indexes.Each sequence implies the multiplication of n coefficients to obtain the probability amplitude.This calculation can be optimized if zeros are present in the multiplication.In a worst case scenario, the computational cost scales as O(nn!) .In spite of its unfavorable scaling properties, it offers the best cost for problems with a small number of photons due its simplicity.

Permanent method
A known approach to obtain the probability amplitude of a photonic circuit outcome is based on the calculation of the permanent [19] of a square matrix.This matrix is built from the same coefficients as the circuit matrix but with a different number of rows and columns.If the amplitude of an output ket is calculated given an input one, columns that correspond to the transformation of a determined mode are repeated as many times as that mode is occupied in the input ket.On the other hand, rows are repeated as many times as their corresponding mode is occupied in the output ket.Various algorithms are available in the literature to perform a permanent calculation [20,21].In SOQCS we calculate permanents using Ryser [22] perm(A) = (−1) n S ⊆{1,...,n} and the Balasubramanian-Bax-Franklin-Glynn [21] formulas implemented in gray code, where for each amplitude the computational cost is O(n2 n ) where n is the number of photons.

Benchmarks
Here we present benchmark results of SOQCS for single thread simulations performed on an Intel I7-10750H@2.60GHzwith 16GB of RAM memory.These benchmarks have been carried out using circuits represented by matrices of random coefficients.These circuits are initialized with states that do not contain more than one photon by channel.The first benchmark consists in measuring the time that takes to calculate an output state.
The results are presented in table 1 where in the Direct F and Glynn F columns are shown the computation times for the straightforward and permanent methods respectively.These methods calculate the amplitude distribution of the output state considering a base of kets with the same number of photons as the input.On the other hand, Direct R and Glynn R are their "restricted" counterparts where only outputs with a maximum of one photon per mode are considered.The smaller the Hilbert space of possible solutions the faster the result is obtained.Both the direct and the permanent method perform similarly up to four photons.Permanent methods scale better for a large number of photons.
Times needed to obtain a selected single probability amplitude are shown in table 2. These results are obtained using similar random circuits and inputs as in the previous benchmarks.
It is feasible to calculate in a reasonable time a single probability amplitude for circuits up to 30 photons.This limit can be improved using parallelism in the permanent calculation.However, the calculation of an output state implies the calculation of the full distribution of probability amplitudes.This imposes a more realistic limit of 7 to 9 photons independently of the use of parallelism.The overhead of performing a parallel calculation is usually unfavorable for this range of photons.Therefore we prefer to parallelize the calculation of different output states.This is useful when calculating output density matrices.

Sampling
Sampling means generating instances from the probability distribution of the output of the circuit.In SO-QCS, it is implemented using the Clifford A algorithm described in ref. [23] and a metropolis algorithm presented in ref. [24].

Clifford algorithm
The Clifford A algorithm [23] relies on the fact that the probability to measure n photons in a particular set of modes depends on a chain of conditional probabilities, p(r 1 , r 2 , ..., r n ) = p(r 1 )p(r 2 |r 1 )...p(r n |r n−1 , ...r 2 , r 1 ) , where r i represents one of the l available modes.
The algorithm to obtain a single sample is structured as a loop.In each iteration of the loop the probability of a photon i to be in a given mode is calculated assuming knowledge of the modes occupied by the previous i − 1 photons.The mode corresponding to photon i is chosen randomly from the probability distribution of that photon to be at each mode.These probabilities are also calculated using permanents.We refer to ref. [23] for further details regarding this algorithm.The loop is finished when all photons have been allocated.This algorithm has a computational cost of O(nl3 n ).

Metropolis algorithm
In the metropolis algorithm proposed in ref. [24] successive samples s with probability P(s) are proposed with probability Pc(s).The transition between these samples is accepted with probability, where s is the previous accepted sample, s ′ the new sample proposed and Pc(s) is taken as the classical distribution where photons are distinguishable particles.The result is a Markov chain that in its stationary regime produces samples that approximate the desired distribution P(s).This algorithm has a computational cost of O(n2 n ) (like Clifford B) and it needs to calculate only one permanent by sample.Note that a burn period and a thinning procedure may be needed to make the chain independent of the initial state and avoid sample correlations respectively.

Example: Sampling comparison with the exact result
In the following example we perform a comparison between the probability distribution obtained by the direct method and the Clifford A sampler.The circuit is initialized in a way that provides probability distributions automatically in both cases.

Imperfections
Simulation methods are used to obtain output states or samples of ideal circuits.However, different kinds of imperfections are relevant to photon emission, photon transfer across a circuit or to photon detection.Some circuit imperfections can be handled without considering any additional mechanism.For example, unbalanced beamsplitters will imply θ 45 • .However, some other imperfections require additional support from the SOQCS library.
In this section we will explain how SOQCS deals with three of the most common imperfections in optical circuits.That is, partial distinguishability of photons, imperfect sources and losses.Additional imperfections related to the detection process will be treated in section 6.

Basic mechanism
Partial distinguishability between photons is related to the partial overlap between their wavefunctions.If the overlap is total they are indistinguishable from each other.Partial distinguishability may happen for many reasons like small relative delays between photons or differences in their frequencies.Using again as an example the beamsplitter, we can extend the definition of a ket as |n ch=0,d=0 , n ch=0,d=1 , n ch=1,d=0 , n ch=1,d=1 ⟩ where ch is the channel degree of freedom and d is a label given to a particular wavepacket that characterizes the photons in that mode.If all photons are indistinguishable only modes that share the same d label will be populated.This results in the following transformation: where photons perfectly bunch in one of the output modes as explained by the Hong-Ou-Mandel (HOM) effect [25].On the other hand, if photons are distinguishable, inputs with different d labels are populated.This leads to an output, where all possible outcomes have the same probability.
Note that if photons are distinguishable the output probability distribution of a circuit can be calculated by classical means.
In the general case of partial distinguishability the wavepacket states are not orthogonal ⟨P i |P j ⟩ 0. They can be transformed into an orthonormal basis using the Gram-Schmidt orthonormalization procedure as suggested in ref. [26]: A software implementation of partial distinguishability of photons is quite involved.We provide an extended discussion of this topic in ref. [27].The key points are summarized below.
If photon wavepackets are included in the mode definition then the simulation methods presented above for ideal circuits can also be used for partially distinguishable photons.Photon wavepackets are treated on the same footing as the channel or polarization degrees of freedom.As a consequence, the implementation of the simulation methods is independent of the distinguishability model.This is an advantage because this independence allows for both calculations to be developed separately as different reusable modules.
In SOQCS packets may be configured to be of Gaussian or Exponential shape.Gaussian packets, are characterized by a central time t i a central frequency ω i and width ∆ω i where t is time.On the other hand, exponential packets, are characterized by a characteristic time t i and a decay time τ x i .Physically, there is an infinite set of wavepackets that represent the spatial degrees of freedom of the photons.To implement partial indistinguishability in SO-QCS photons are restricted to a discrete finite subset.This discretization allows to define both an emitter and a delay optical gates in a similar manner as a linear optical element.This means that they are represented by their respective matrices.
The emitter matrix contains the transformation to orthogonal packets derived from eq. 18, On the other hand, time may be split into different periods where the same set of packets is defined in each period but displaced in time.Therefore, the orthonormalization of packets is also the same at those different periods.We assume that the time between periods is large in comparison with packets width.A delay of one photon by one period of time implies the substitution of its packet label for the equivalent label in the next period, TD = The delay matrix reflects the relabelling of packets at different discrete periods.
Although those matrices are non-unitary these elements can be handled in the same way as any other linear optical element.The reason is that in this kind of simulation optical circuits are interpreted as one way operations from the input to the output.For further details see ref. [27].

Photon as an abstraction
To calculate partial distinguishability effects in SO-QCS the Gram-Schmidt coefficients have to be known.However, the user would likely find cumbersome to define those coefficients directly or to create a code to do so.Instead, SOQCS only requires the user to define how many photons exist in each of the channels and which are the parameters given to those photons.To allow for this, we have implemented a stack of abstraction levels with the photon abstraction at the top as shown in fig. 5 and the emitter and the delay definitions given by the Gram-Schmidt coefficients at the bottom.Each of the levels is responsible in organizing the information and performing the calculations needed to fill the gap between the two equivalent sets of information.The photon description is understandable for a user while the Gram-Schmidt coefficients are the information needed by the computer to perform the actual calculation.These levels are (in a brief description), • Level 3: Photon abstraction level.This is the upper abstraction level where photons are described by their paramenters (channel, time, frequency, etc).This information is splitted into the definition of the initial state of the circuit and a list of packets with their parameters.
• Level 2: Packet abstraction level.In this level packets are treated as virtual circuit elements used to create a configuration table with the information of all packets.
• Level 1: Shape model.In this level the information stored in the configuration table is used to calculate the overlap matrix between packets.
• Level 0: Finally, in the bottom level the overlap coefficients are used to perform the Gram-Schmidt orthonormalization procedure.As a result, the emitter matrix to perform the transformation from the non-orthonormal basis to the orthonormal one is obtained.
The photon abstraction level is implemented in SO-QCS as part of the device abstraction.The class qodev is defined as an ensemble of an input state, a quantum optical circuit and a list of packets.In this context, photons can be defined as virtual elements of the quatum circuit in a natural way without the need of defining separately the initial state and the packets related to the circuit.These elements are declared in the same way as the linear optical elements with the difference that they serve the purpose of configuring the simulation and have no matrix representation attached.The class qocircuit implements the circuit itself at a more fundamental level.In this class packets may be added to the circuit definition directly (level 2) also as virtual elements.
In general, the function of the qodev class is to allow for a more direct definition of photons.However, it also automates some aspects of the measurement process configuration.We will expand this matter further in section 6.

Example: Partial distinguishability of photons
Next, we present the code for a circuit simulation where two groups of five photons are set up in each of the input channels of a beamsplitter.These two groups of photons arrive to the beamsplitter with a relative delay ∆t.The wavefunctions of both sets of photons are modeled as Gaussians of frequency ω = 1 and uncertainty ∆ω = 1 in arbitrary adimensional units.The code is written in a function called HOMP that returns the probability of the outcome specified by the parameter args for a delay ∆t given by the parameter dt.Detectors are configured as counters that measure the total number of photons independently of their detection time.Outcome probabilities are calculated automatically from the amplitudes of the output state.We will expand on how detectors are implemented in SOQCS in section 6.

Sources
Physical photon emitters have imperfections.Therefore, emitted photons have to be represented by a density matrix instead of a pure state.In particular, SO-QCS implements a parametrized quantum dot able to simulate the creation of pairs of entangled photons of features compatible with those emitted from a physical bi-exciton XX exciton X cascade.

Quantum dot emission model
From [28,29,30] the density matrix of photons emitted from a quantum dot cascade is, This matrix is written as the weighted sum of pure states density matrices.Each of these pure state density matrices are related with a physical process in the quantum dot.Here, ρ en is the density matrix of the entangled photons emitted in the ideal bi-exciton and exciton cascade.ρ pd is the density matrix of the photons where some cross-dephasing event played a role in their emission and finally, ρ noise is the density matrix created by random photons from external sources of noise or photons emitted through any non-coherent process.The probability of each outcome is a function of k that is the fraction of photons emitted by the quantum dot source.On the other hand, (1 − k) is the fraction caused by noise.p s = e − t tss is the probability of the dot to emit coherently, t ss being the characteristic time of coherence.Finally, 1 − p d = 1 − e − t t HV is the probability of a dephasing event where t HV is the cross dephasing characteristic time.
Here is where the stochastic nature of the simulation plays an important role.To calculate the output density matrix of a circuit, the possible input pure states are sampled according to their corresponding probabilities and simulated individually.The output states are then used to calculate the resulting density matrix.The main advantage of this scheme is that it can be implemented using computation methods developed for pure states and ideal circuits while achieving our objective of modularity.
The density matrix of the entangled photons ρ en is given by, in the basis {|HH⟩, |HV⟩, |V H⟩, |VV⟩} where H and V are the "horizontal" and "vertical" polarizations respectively.The density matrix above represents the pure state, where S is the fine structure splitting energy between the states in the horizotally HH and vertically VV polarized photon emission cascades and ∆ t is generated randomly with an exponential distribution.
On the other hand, the density matrix ρ pd does not contain linear superpositions between differently polarized states but photons are still emitted correlated in polarization, The density matrix is therefore equivalent to a mixture of the states, each one emitted with 1/2 probability.Finally, the noise density matrix ρ noise , models the random emission of photons in any of the four possible pure basis states with equal probability 1/4.These are, The process of generating a list of states compatible with the quantum dot emission density matrix consists in choosing the different paths of the next branching tree with the stated probabilities, The process is repeated each time a new pair of photons is emitted.The output density matrix of a circuit obtained from stochastically sampling an imperfect source is calculated as the sum of all the individual output states, Note that in the case of experiments with fluctuating emission parameters, it is needed to structure the output to avoid the infinite growth in the number of packets.That is, if in each simulation the wavepackets are different, then after a large number of samples that number becomes also large.Fortunately, there are many ways to handle this situation.If the detectors are counters, the particular packets, modeling the time of arrival of photons to those detectors, are irrelevant.On the other hand, detection may be performed or treated in a discrete manner.This discretization into intervals of detection reduces the number of possible outcomes in relation to time (or any other quantity) to a finite manageable size.
A density matrix is represented in SOQCS by the C++ object dmatrix.Pure states can be stored/added to dmatrix like in a list but internally a matrix representation is maintained.The matrix can be printed at any moment, updated with more states and printed again with improved statistics.Entries in the matrix (rows and columns indexes) are created dynamically to represent only entries with non-zero values.Internally a dictionary is maintained to recognize if there is already a preexisting entry to a particular coefficient that has to be updated or if a new entry (row and column) has to be created.
Note too that density matrices may also be needed in cases where the input to the circuit is a pure state.This may happen when post-selection conditions are multivalued.Once again, we will expand upon this in section 6.

Entanglement swapping protocol
In the entanglement swapping protocol [30] two pairs of photons are emitted by a quantum dot in entangled Bell states at two different times.Each pair is emitted from a cascade where one of the photons is created from a bi-exciton XX and the other from an exciton X decay.
The exciton X photon of each pair is sent to a beamsplitter, each one into a different channel.The first exciton photon is delayed using a Mach-Zehnder interferometer.This can be modeled using a delay or, in this case, assuming for simplicity a delayed time as the emission time of the first exciton photon.If one photon is detected in each output channel of the beamsplitter then the remaining two photons become entangled.Ideally, this condition is only met when the two photons that arrive to the beamsplitter are in a Ψ − state.However, the quantum dot does not emit perfect Bell states and some fluctuations in the photons emission times also applies.Therefore, the resulting density matrix of the remaining bi-exciton XX photons will not correspond to a pure state of a perfectly entangled pair of photons.
The simulation is repeated as many times as necessary with different initial states and the results are stored in the density matrix to reconstruct the statistics.The protocol is successful if two photons are measured at channels 1 and 2 independently of their detection time.
On the other hand, the delay between photons in channel 0 is large enough to label them while disregarding small variations in their particular detection times among different simulations.
The density matrix below is obtained from repeated runs of the presented code.The result is consistent with analytical calculations using equivalent random and cross-dephasing noises and fine structure splitting values of the quantum dot while disregarding pure dephasing effects, | H(0)0, H(2)0 > 0.2251 0.0000 0.0000 0.0000 | H(0)0, V(2)0 > 0.0000 0.2801 -0.0225 0.0000 | H(2)0, V(0)0 > 0.0000 -0.0225 0.2703 0.0000 | V(0)0, V(2)0 > 0.0000 0.0000 0.0000 0.2245 where the first column indicates the basis vectors corresponding to each row of the density matrix.Figure 7: Schematic of a circuit implementing the entanglement swapping protocol as discussed in ref. [30].XX are bi-exciton and X exciton photons emitted by a quantum dot modeled as an imperfect source of photons.Detectors in channels one and two impose a post-selection condition of one photon detected in each channel. .

Model
Photon losses imply that the total probability of obtaining an outcome with the same number of photons as the input is smaller than one.If the simulator calculation only considers outputs with the same number of photons as the input then output states are obtained normalized to values smaller than one.Therefore, a circuit with losses is described by a non-unitary matrix.An example can be provided by a lossy beamsplitter which is discussed below.
SOQCS implements a beamsplitter model of losses in which an extra virtual loss mode is added for each physical mode present in the circuit.The non-unitary matrix of the circuit is thus extended to include those extra degrees of freedom.Even if losses are defined locally (for each optical element) the loss model has to consider the entire circuit matrix.The reason is that we can not attribute a photon loss to any particular circuit element with certainty but all the possibilities have to be considered simultaneously and in linear superposition.Therefore, the matrix representing the entire circuit extended to include the virtual loss modes has to be unitary.If a photon is not transmitted then it is lost and has to appear in one of the virtual loss modes thus conserving probability.This way, we can use the already mentioned simulation methods for ideal circuits in which the number of photons at the output is the same as in the input.As a result, the implementation of the calculation of losses is independent of the chosen simulation method allowing us to achieve our goal of modularity.
It is always possible to construct a 2n × 2n unitary matrix U from a smaller n × n non-unitary matrix M if the singular values of M are smaller than one.The singular value decomposition can be expressed as M = RDV where D is a diagonal matrix and R and V are unitary matrices.The 2n × 2n matrix is built as, Finally, after a simulation, virtual loss channels can be traced out from the output state to obtain the circuit density matrix.This method can be applied independently of the circuit structure and parameters and it is used to obtain the amplitude of probability of m photon outcomes given n > m photon inputs.

Example: Lossy beasmplitter
We can test the library for the particular case of a small thin dielectric film using the following equations [31], âout = t âin + r bin + Fa bout = t bin + r âin + Fb , where we re-interpret the loss operators Fa and Fb as virtual channel destruction operators (or at least their normalized counterparts Fa = 1 − |r| 2 − |t| 2 fa ).Fa and Fb do not commute and their self commutators are not one, This thin dielectric act as an imperfect beamsplitter with losses.The dielectric is modeled as the non-unitary gate, that is extended following the previously established prescription, where . For this particular case, the dielectric matrix is normal therefore R = V † .D is the eigenvalue matrix and R contains the eigenmodes of the circuit.Hence, in this case we can interpret √ 1 − D 2 as the coupling coefficients between each of those eigenmodes and a single loss mode.
The analytical result presented in ref. [31] for different transmission amplitudes can be reproduced numerically in SOQCS using this matrix as shown in fig.8.

Detectors
SOQCS handles detectors as virtual circuit elements that can be defined in the same manner as the linear optical elements.However, the process of detection is a non-unitary irreversible operation.Their virtual definitions are used to specify the parameters needed to configure the post-selection conditions of the circuit and the different models of imperfections in the detection process.

Post-selection
In general, the probability of measurement of a determined outcome is, where ρ is the density operator and Π i = |Φ i ⟩⟨Φ i | is the projector operator.Post-selection is a partial measurement where only some of the degrees of freedom that describe a state are measured.If a certain outcome is fulfilled then the output obtained in the remaining modes is accepted (otherwise it is rejected).Thus, the result of post-selection is a reduced density matrix for the remaining degrees of freedom, If the output of a circuit is calculated repeatedly in presence of noise or with different inputs (ex.non ideal source) then we will obtain a different output each time.
In this case, the original unreduced density matrix can be re-interpreted as the average of density matrices of each of the output pure states ρ = 1 N N j=1 ρ j .As a consequence, the reduced density matrix can be rewritten as the average of each pure state matrix reduced by postselection, This means that we can calculate the overall reduced density matrix as the sum of already post-selected individual state projectors, Therefore to perform post-selection in SOQCS, • First, the output state of an optical circuit is obtained • Second, post-selection is performed applying a projector to the output state ( defined for a subset of degrees of freedom of that state).
• Third, the resulting reduced state is added to a density matrix in the same manner as it is done in circuits with no post-selection.
We can interpret the post-selection operation as a function that relates a state with another state defined in a reduced set of degrees of freedom.For example, if we perform a post-selection operation over the state, projecting the 2nd and 3rd degrees of freedom (usually channels) to the projector |Ψ − ⟩ 12 the resulting reduced state is | Ψ⟩ = |Ψ − ⟩ 03 .However, post-selection is not always defined for a particular projector.It is usual to demand certain conditions (for example occupation values) to be fulfilled in certain channels independently of their polarization, frequency or other characteristics.In this case more than one projector will fulfill the output condition.For example, in the entanglement swapping protocol explained above the post-selection condition is the detection of one photon in both channels one and two independently of photon polarization.Therefore, the projectors considered can be |HH⟩, |HV⟩, |V H⟩, |VV⟩ where photons may be "vertically" V or "horizontally" H polarized or, alternatively, |ϕ + ⟩, |ϕ − ⟩, |ψ + ⟩, |ψ − ⟩.The reduced density matrix is calculated as the sum of reduced density matrices obtained from applying each projector, SOQCS determines automatically the projectors needed, given the output condition, and performs the calculation of the density matrix automatically.

Detector abstraction
Measurement may depend on the kind of detector used.For example, if the detector is only a counter of photons to be read at the end of the experiment then that detector is blind to the particular characteristics of the photon wavepackets (time, frequency, etc).Therefore, the degrees of freedom of the resulting density matrix have to be reduced to obtain a new state description without the packet degree of freedom.This does not mean that the packet degree of freedom is unimportant in the circuit outcome.For example, in fig.6 we can see the probability of having five photons in two different channels of a circuit where detectors are configured as counters.The counting in different channels is different depending on whether they are synchronized due to HOM.A relabeling of the output to remove undesired degrees of freedom is always possible for the diagonal elements of the density matrix while off diagonal elements are more involved and depend on the particular detector definition.
In general, there are three kinds of detectors defined in SOQCS, • Counters.All photons in a mode are summed independently of their packet degrees of freedom.
• Timed detection.The packet index is dropped in favor of a time label.
• Full detection.We do not perform any further process to the output and the result is expressed in terms of the orthonormalized packets.
Configuring a quantum device to use the proper kind of detector will provide automatically the desired probability distribution (and in some cases also the density matrix).

Imperfect physical detectors
Detectors are physical devices with their own imperfections.Detector imperfections affect the outcome probabilities of a circuit.For example, the detector efficiency is related to losses within the detector.In SOQCS, this can be calculated automatically by adding lossy medium to the circuit when the detector is defined.
With the exception of losses the rest of the detection imperfections can be calculated with classical error models.These errors alter the distribution of probability of the circuit outcomes.(Distribution probabilities are the diagonal coefficients of the density matrix.Errors are undefined for the off diagonal coefficients because they have no classical meaning).In SOQCS there is a second C++ class that handles distributions of probabilities in a equivalent way to the diagonal elements of the density matrix.The probability bins class p_bin is more efficient and faster if we are uninterested in the coherences for a particular problem.It is in this class where physical detector effects are implemented.
The classical error models implemented in SOQCS are, • Dead time.After a photon measurement the detector is off for a certain time.We model this effect by rejecting a detector measurement randomly with a certain probability.We define this probability as the fraction of time that a detector is off with respect to the simulation time.This model provides a very simple approximation.We are assuming that the main cause of dead time is noise with no correlation between detectors.This is normally true if the simulation time is large enough in comparison with the dead time.
• Dark counts.A detector may observe photons emitted by itself.The dark count rate is the average rate of registered counts without any incident light.These false detections are mainly of thermal origin.Therefore, photons are added to the resulting outcome with a probability determined by a Poissonian distribution, where n is the number of photons (dark counts) to be added and λ is the rate (average value) of dark counts in the simulation time.
• Noise.Finally, the output probability may be altered with a random Gaussian noise to account for a wide range of sources of noise present in experimental setups.
The overall detection procedure computed in SOQCS is implemented as, 1. Calculate the probability distribution.Note that the computation of losses adds some auxiliary extra loss modes.These modes can be eliminated from the probability distribution in step 4th performing a post-selection procedure for all the possible outcomes in those channels.

Conclusions
In this paper we report the implementation of the optical circuit simulator SOQCS that uses an approach akin to stochastic quantum trajectories.The simulator has been built in a modular way and allows for multiple simulation methods.These methods provide as results output states, probability distributions and density matrices depending on user's configuration.Post-selection conditions and measurement configurations are implemented using detectors as virtual circuit elements.Furthermore, partial distinguishability, delays, losses, circuit and detector imperfections are parametrized and calculated automatically.Simulations can be carried out using SOQCS as a C++ library or by means of its Python port that also contains a graphical circuit plotter.In conclusion, SOQCS is a framerwork for the simulation of quantum optical circuits that has been built with flexibility in mind and with a modular code that can be further expanded to include additional simulation methods or new imperfection models.csign10.detector(2)csign10.detector(3) In this particular case, we are interested in a phase flip, and hence we instruct SOQCS to provide a state as output.Other options like probability distributions or density matrices are available using different instructions.The resulting output is shown next.Note that in this case we have instructed SOQCS to provide a normalized output, | 0, 0 >: 0.50000000 + 0.00000000 j | 0, 1 >: 0.50000000 + 0.00000000 j | 1, 0 >: 0.50000000 + 0.00000000 j | 1, 1 >: -0.50000000 + 0.00000000 j

Figure 2 :
Figure 2: Diagram of the software structure of the SOQCS library.SOQCS is built in layers on top of Eigen3[17] library for linear algebra, with the level of abstraction increasing upwards.Each layer contains C++ classes that use other more basic classes defined on lower layers.

Figure 3 :
Figure 3: Diagram of an arbitrary optical circuit made of three beamsplitters.

|𝜓Figure 4 :
Figure 4: NSX circuit diagram as discussed in ref. [3].Except for the ket |Ψ⟩ symbol this figure has been created by SOQCS graphical module.Emitters are plotted as triangles while detectors are implemented as square like shapes with one round side.The optical elements are represented as boxes.PS stands for phase shifter while BS for beamsplitter.

Figure 5 :
Figure 5: Stack of abstraction levels used to implement partial distinguishability of photons in SOQCS library.

Figure 8 :
Figure8: Outcome probability of a lossy beamsplitter a) in case that two identical photons are set-up in the same input channel |2, 0⟩ and t = ir.b) in case that two identical photons are set up, one in each input channel, |1, 1⟩ and t = r.Both calculations carried out in SOQCS reproduce the results for the same device presented in ref.[31].

Figure B. 9 :Figure
Figure B.9: The CZ circuit [5].Plot is generated automatically by SOQCS.The modes 1 and 3 are not acted upon by the circuit and hence are omitted from the plot to avoid a clutter.

Table 1 :
Calculation time of each method implemented in SOQCS to obtain the output state of a random circuit.

Table 2 :
Calculation time of each method implemented in SOQCS to obtain the probability amplitude of a single outcome.
[26,27]nt outcome probabilities for each delay time ∆t are calculated by means of repeated calls to HOMP and plotted in fig.6.This example is similar, but with more photons, to examples discussed in[26,27].