Basic usage

This tutorial will introduce you to the definition of the building blocks for your boson sampling experiment. The general workflow for a simple simulation is to define an Input that enters into a Interferometer and ask what is the probability to get a defined OutputMeasurement.

They are linked together through an Event type, which holds the respective probabilities. As the computation of probabilities is often the most time consuming step, you need to explicitly ask for it through compute_probability! which updates the EventProbability data.

Input

BosonSampling.jl provides three distinct types of input depending on the distinguishability of the particles we want to make interfere: Bosonic, PartDist and Distinguishable. The type PartDist is a container for different models of partial distinguishability. Currently available models are:

• OneParameterInterpolation
• RandomGramMatrix
• UserDefinedGramMatrix
• Undef

In order to define the input, we first need to provide a ModeOccupation that describes the repartition of the particles among the modes.

julia> n = 3; # photon number

julia> m = 6; # mode number

julia> my_mode_occupation = ModeOccupation(random_occupancy(n,m))
 state = [1, 0, 0, 1, 1, 0]

In the example above, my_mode_occupation has been created thanks to random_occupancy that randomly places n particles among m modes. Here we have one particle in the first, fourth and fifth modes. Let's build an input made off indistinguishable photons by using the type Bosonic

julia> my_input = Input{Bosonic}(my_mode_occupation)
 Type:Input{Bosonic}
 r:state = [1, 1, 0, 0, 1, 0]
 n:3
 m:6
 G:GramMatrix{Bosonic}(3, ComplexF64[1.0 + 0.0im 1.0 + 0.0im 1.0 + 0.0im; 1.0 + 0.0im 1.0 + 0.0im 1.0 + 0.0im; 1.0 + 0.0im 1.0 + 0.0im 1.0 + 0.0im], nothing, nothing, OrthonormalBasis(nothing))
 distinguishability_param:nothing

where my_input holds the information defined above and an additional field, the GramMatrix:

help?> GramMatrix

 Fields:
  - n::Int: photons number
  - S::Matrix: Gram matrix
  - rank::Union{Int, Nothing}
  - distinguishability_param::Union{Real, Nothing}
  - generating_vectors::OrthonormalBasis

which contains everything about the distinguishability of the particles within my_input. The matrix itself can be accessed via the field S:

julia> my_input.G.S
 3×3 Matrix{ComplexF64}:
 1.0+0.0im  1.0+0.0im  1.0+0.0im
 1.0+0.0im  1.0+0.0im  1.0+0.0im
 1.0+0.0im  1.0+0.0im  1.0+0.0im

One can do the same for Distinguishable particles placed in the first_modes

julia> my_mode_occupation = first_modes(n,m);

julia> my_input = Input{Distinguishable}(my_mode_occupation);

julia> my_input.G.S
 3×3 Matrix{ComplexF64}:
 1.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  1.0+0.0im

We can move now to the PartDist case with a model of partially distinguishable particles defined by a RandomGramMatrix

julia> my_input = Input{RandomGramMatrix}(first_modes(n,m));

where my_input.G.S is a randomly generated Gram matrix.

Finally, one can resort to a OneParameterInterpolation model taking a linear distinguishability parameter as an additional argument in the definition of my_input:

julia> my_mode_occupation = ModeOccupation(random_occupancy(n,m));

julia> my_distinguishability_param = 0.7;

julia> my_input = Input{OneParameterInterpolation}(my_mode_occupation, my_distinguishability_param)
 Type:Input{OneParameterInterpolation}
 r:state = [1, 1, 1, 0, 0, 0]
 n:3
 m:6
 G:GramMatrix{OneParameterInterpolation}(3, [1.0 0.7 0.7; 0.7 1.0 0.7; 0.7 0.7 1.0], nothing, 0.7, OrthonormalBasis(nothing))
 distinguishability_param:0.7

Notice that the Bosonic Gram matrix is recovered for my_distinguishability_param = 1 while we find the Distinguishable case for my_distinguishability_param = 0.

Interferometer

The second building block of our boson sampler is the interferometer we want to apply on my_input. A common practice to study boson sampling is to pick up at random a Haar distributed unitary matrix that will represent the interferometer. This can be done as follow:

julia> my_random_interf = RandHaar(m);

julia> my_random_interf.U
 6×6 Matrix{ComplexF64}:
 -0.398793-0.162392im     0.500654+0.239935im   -0.0822224-0.189055im     0.361289+0.048139im   -0.0639807+0.521608im   -0.0608676-0.226158im
 0.331232+0.312918im    -0.362433-0.0585051im    0.172619+0.157846im       0.6224-0.0656408im   -0.186285+0.215539im    -0.233294-0.274952im
 -0.085377+0.00861048im  -0.427763+0.1271im      -0.140581-0.541889im   -0.0407117+0.219472im     0.499523-0.0486383im  -0.0711764-0.416309im
 -0.180636+0.294002im    -0.376353+0.0943096im   -0.489498-0.206616im   0.00169099-0.221399im    -0.260862+0.305118im     0.313454+0.373737im
 0.619915-0.0776736im     0.29879-0.323099im    -0.398401-0.371214im    -0.132369+0.0411683im   -0.242868+0.0913286im  -0.0282651-0.179273im
 0.179322+0.240927im    0.0369079+0.110875im     0.100647-0.0206654im   -0.153966+0.577663im     0.154379+0.372127im    -0.366647+0.481086im

where we have accessed to the matrix thanks to the field .U. We may also need to use a specific interferometer such as a Discrete Fourier Transform or the Hadamard transform:

julia> my_fourier_interf = Fourier(m);

julia> is_unitary(my_fourier_interf.U)
 true

julia> my_hadamard_tf = Hadamard(2^m);

julia> is_unitary(my_hadamard_tf.U)
 true

where we have checked the unitarity thanks to is_unitary.

The implemented interferometers are listed in Interferometer but it is still possible to define our own unitary by resorting to the type UserDefinedInterferometer:

julia> sigma_y = [0 -1im; 1im 0];

julia> my_interf = UserDefinedInterferometer(sigma_y)
 Type : UserDefinedInterferometer
 m : 2
 Unitary :
 ┌────────┬────────┐
 │ Col. 1 │ Col. 2 │
 ├────────┼────────┤
 │  0+0im │  0-1im │
 │  0+1im │  0+0im │
 └────────┴────────┘

OutputMeasurement

Now consider what you want to observe, in this numerical experiment. If looking at the case of a single output, we would use an OutputMeasurement type called FockDetection. Other types are currently defined such as PartitionCount, which would evaluate the probability of finding a photon count in a partition of the output modes.

Similary to the definition of the Input, it is also possible to define an output configuration from a ModeOccupation

julia> n = 3;

julia> m=n^2;

julia> my_mode_occupation = first_modes(n,m);

julia> my_input = Input{Bosonic}(my_mode_occupation)
 Type:Input{Bosonic}
 r:state = [1, 1, 1, 0, 0, 0, 0, 0, 0]
 n:3
 m:9
 G:GramMatrix{Bosonic}(3, ComplexF64[1.0 + 0.0im 1.0 + 0.0im 1.0 + 0.0im; 1.0 + 0.0im 1.0 + 0.0im 1.0 + 0.0im; 1.0 + 0.0im 1.0 + 0.0im 1.0 + 0.0im], nothing, nothing, OrthonormalBasis(nothing))
 distinguishability_param:nothing

julia> out = FockDetection(my_mode_occupation)
 FockDetection(state = [1, 1, 1, 0, 0, 0, 0, 0, 0])

using FockDetection. Additionally, we can define an Event that stores our input-interferometer-output content

julia> my_interf = Fourier(my_input.m)
 Type : Fourier
 m : 9

julia> ev = Event(my_input, out, my_interf)
 Event{Bosonic, FockDetection}(Input:

 Type:Input{Bosonic}
 r:state = [1, 1, 1, 0, 0, 0, 0, 0, 0]
 n:3
 m:9
 G:GramMatrix{Bosonic}(3, ComplexF64[1.0 + 0.0im 1.0 + 0.0im 1.0 + 0.0im; 1.0 + 0.0im 1.0 + 0.0im 1.0 + 0.0im; 1.0 + 0.0im 1.0 + 0.0im 1.0 + 0.0im], nothing, nothing, OrthonormalBasis(nothing))
 distinguishability_param:nothing, FockDetection(state = [1, 1, 1, 0, 0, 0, 0, 0, 0]), EventProbability(nothing, nothing, nothing), Interferometer :

 Type : Fourier
 m : 9)

and then one can compute the probability that this event occurs

julia> compute_probability!(ev)
 0.015964548319225575

Those steps can be repeated for different types of input. Let's say we want to compute the probability that partially distinguishable photons populating the n=3 first modes of m=9 modes end up in the n last output modes when interfering through a random interferometer:

julia> my_input = Input{RandomGramMatrix}(first_modes(n,m)); # input from a randomly generated Gram matrix

julia> out = FockDetection(ModeOccupation([0,0,0,0,0,0,1,1,1]));

julia> my_interf = RandHaar(m);

julia> ev = Event(my_input, out, my_interf);

julia> compute_probability!(ev)
 0.014458823860031098

Using the BosonSampling types

Julia allows to define functions that act on new types, such as ModeOccupation defined in this package, through a syntax that would otherwise be reserved for core-objects such as Float, Int.

This allows to intuitively act on custom types. For instance, two ModeOccupation can see their state summed by simply using +

n=3
m=4
s1 = ModeOccupation([1,2,3,4])
s2 = ModeOccupation([1,0,1,0])

julia> s1+s2
state = [2, 2, 4, 4]

Some functions of interest are zeros(mo::ModeOccupation), Base.cat(s1::ModeOccupation, s2::ModeOccupation) for instance.