The Complexity of Linear Optics

Version: 1 Revised: 10/29/11
© 2011 Justin Dove <dove@csail.mit.edu>
Obtain permission before redistributing or editing.

This notebook implements the tools to carry out numerical simulations relating to the work of Aaronson and Arkhipov under the same title. Cells that define computational tools are marked as initializations cells and will automatically evaluate. All other cells such as examples, demonstrations, and numerical calculations must be evaluated manually.

Reck et. al. Algorithm Implementation

This section implements the notation and algorithm devloped by Reck. et. al. for implementing an arbitrary finite dimensional unitary operator via linear optics. First, we define some notation and functions that will be used for the algorithm.

Dimension=4; (*"Dimension" is a variable that holds the value of the number of modes in the optical network, or equivalently the dimenion of the space the unitary operator acts on. Here it is given an arbitraty value which can and will be changed as needed.*)

(Subscript[T, p_,q_])[ω_,ϕ_,n_Integer]:= ReplacePart[IdentityMatrix[n],{{q,q}-> (E^(I*ϕ))*Sin[ω],{q,p}-> (E^(I*ϕ))*Cos[ω],{p,q}-> Cos[ω],{p,p}-> -1*Sin[ω]}]
(*This defines a matrix describing the action of a generalized optical element via the notation developed by Reck et. al. Subscripts may be entered by typing [CTRL]+[-]. "p" and "q" denote the two modes that are being operated on where p > q. The input "ω" denotes the beam splitter angle. The input "ϕ" denotes the phase shift on the "q" output mode. The input "n" denotes the dimension of the space. The matrix generated by this function will act as an identity on all modes except p and q.*)

(Subscript[T, p_,q_])[ω_,ϕ_]:= (Subscript[T, p,q])[ω,ϕ,Dimension]
(*This defines what to do when "n" is not explicitly input. In such an event, the Dimension parameter is used.*)

Phases[phases_List]:=DiagonalMatrix[(E^(I*phases))]
(*This function takes in a list of phases in the form {(Subscript[ϕ, 1]),(Subscript[ϕ, 2]),(Subscript[ϕ, 3]),(Subscript[ϕ, 4]),...} and generates a matrix that implements the appropriate phase shifts on each mode. The first element in the list is the shift on mode 1 and so on.*)

Having defined the basic functions, we now demonstrate them via the example unitary given at the end of Reck et. al. by reconstructing it from the decomposition they give.

Dimension=4; (*It is a 4\[Cross]4 unitary.*)
ex=((Subscript[T, 4,3])[0,0].(Subscript[T, 4,2])[π/4,π/2].(Subscript[T, 3,1])[π/4,π/2].Phases[{0,0,π,π}])\[ConjugateTranspose];
(*Generates the unitary from the decomposition via Equation (7) of Reck et. al. and stores it in "ex" to use later. That tiny superscripted symbol at the end is a dagger (conjugate transpose), written by [ESC]ct[ESC]. Note, (Subscript[T, p,q])(0,0) implements a swapping of modes p and q.*)

ex//MatrixForm (*Displays the generated matrix*)

Indeed, the correct unitary is generated. Next, we define a function “Optics” that will take in an arbitrary n×n unitary and generate a set of parameters describing the beam splitters and phase shifters that implement it as an n mode linear optical network. This is the implementation of the algorithm discussed in Reck et. al. However, there is one convention to note. The way the algorithm is described in their paper, they right multiply the desired unitary by appropriate (Subscript[T, p,q])(ω,ϕ) until they get a diagonal matrix. Then the right multiply it by an appropriate phase shifting matrix to get the identity. By doing it this way, the decomposition they generate actually implements the inverse of the unitary. Thus, they tell you that you should take the conjugate transpose of this to get the actual unitary. This means that to implement it optically, you have to switch the inputs and outputs (the equivalent of reodering the product backwards) and conjugate transpose all the elements (which amounts to negating all phase shifts). Instead of putting this extra work on the experimentalists, we make a simple change of convention. We first take the conjugate transpose of the inputted unitary, and then run the Reck et. al. algorithm on that. This way, the elements that are generated will implement the original unitary with no need for any extra work. If you instead wish to generate the same exact results that Reck et. al. would have gotten, then simply input the conjugate transpose of your unitary. Otherwise, input your unitary, get out the optics.

Optics[unitary_]:=Module[{Utemp,n,optics,p,q,sol,ω,ϕ,i},(*Localizes some variables*)
Utemp=unitary\[ConjugateTranspose]; (*Inverts the input unitary and stores it in a temporary variable*)
n=Dimensions[Utemp][[1]]; 
(*Sets "n" equal to the number of rows in the unitary. "Dimensions[Utemp]" outputs a list of the two dimensions (row and column) of Utemp. "[[1]]" picks out the first element in that list*)
optics=Array[0&,{n,n}]; (*Initializes an n\[Cross]n matrix of all 0's to hole the optics parameters.*)
For[p=n,p>1,p--,
For[q=p-1,q>0,q--,
(*Runs through all values of p,q from n,n-1 down to 2,1 which denote all the intersections of modes.*)

Clear[ω,ϕ]; (*Clears any previous values of ω and ϕ for safety*)

(*Next we find the beamsplitter and phaseshifter angles needed to zero off (Subscript[Utemp, p,q]). Since this zeroing off is accomplished by right multiplying Utemp by (Subscript[T, p,q])(ω,ϕ), the problem amounts to finding ω and ϕ that solve the equation (e^(iϕ))sin(ω)(Subscript[Utemp, p,q])+cos(ω)(Subscript[Utemp, p,p])=0. To solve this equation, we need to pick a ϕ that will rotate the complex number (Subscript[Utemp, p,q]) into the real span of (Subscript[Utemp, p,p]). Then we have to find ω that will apprpriately balance out the magnitudes of the two complex numbers so that they cancel out. This is equivalent to picking ω = (Tan^(-1))(-((Subscript[Utemp, p,p])/((Subscript[Utemp, p,q]) (e^(iϕ))))) by rearranging the equation.*)

If[Utemp[[p,q]]==0,ω=-(π/2);ϕ=π;,(*If (Subscript[Utemp, p,q]) is already 0, then we don't need any element at the intersection of modes p and q. We then set ω=-π/2 and ϕ=π becuase (Subscript[T, p,q])(-π/2,π) is just the identity matrix.*)
(*Otherwise...*)
ϕ=Arg[Utemp[[p,p]]]-Arg[Utemp[[p,q]]];(*Pick ϕ to be the angle of separation between (Subscript[Utemp, p,q]) and (Subscript[Utemp, p,p]) in the complex plane. That is to say the difference in their arguments. This will then assure that (e^(iϕ)) rotates (Subscript[Utemp, p,q]) onto (Subscript[Utemp, p,p]) as necessary. There is however some choice in this angle. First of all, we can arbitrary add any multiple of 2π to it as that will not change the rotation. Less trivially, we can add any multiple of π because that will have the effect of negating the term (rotating it 180°) which will still keep it in the real span of (Subscript[Utemp, p,p]), and will still allow the two terms to be cancelled out by an appropriate ω.*)

Which[(*Here we use the available freedom in our choice of ϕ to keep it between 0 and π. This code could be written much more succinctly, but we sacrifice succinctness for the sake of transparency and ease of editing.*)
-2 π≤ϕ<-π, ϕ=ϕ+2 π;ω=ArcTan[-((Abs[Utemp[[p,p]]])/(Abs[Utemp[[p,q]] (E^(I ϕ))]))];,(*Add 2π to ϕ between -2π and -π*)

-π≤ϕ<0,ϕ=ϕ+π;ω=ArcTan[((Abs[Utemp[[p,p]]])/(Abs[Utemp[[p,q]] (E^(I ϕ))]))];,(*Add π to ϕ between -π and 0*)
0≤ϕ<π,ω=ArcTan[-((Abs[Utemp[[p,p]]])/(Abs[Utemp[[p,q]] (E^(I ϕ))]))];,(*Leave ϕ between 0 and π alone.*)
π≤ϕ≤2 π,ϕ=ϕ-π;ω=ArcTan[((Abs[Utemp[[p,p]]])/(Abs[Utemp[[p,q]] (E^(I ϕ))]))](*Subtract π from ϕ between π and 2π*)

(*After having chosen ϕ, we could pick ω = (Tan^(-1))(-((Subscript[Utemp, p,p])/((Subscript[Utemp, p,q]) (e^(iϕ))))) regardless of which case we had. However, general precision issues could cause ω to have a very small imaginary part. If we chose to do this, we would have to use the Chop[] function to keep ω real. However, explicitly using rounding techniqies in our function is an inelegant solution to circumventing rounding errors in Mathematica, so instead we look at the ArcTan of the magnitude of the appropriate terms, that way ω is assured to be completely real. In doing this, we also have to multiply by the ArcTan argument by -1 for the cases where the two terms are 180° apart in the complex plane because the Abs[] strips this away. Thus, we choose our ω appropriately in the case statement that chose ϕ. In principle taking the magnitude and multiplying by -1 when necessary is fine because the ArcTan[] arguement should be completely real due to the choice of ϕ. However, precision issues can make this only approximately true, so using Abs[] does amount to rounding in some sense. That said, the rounding is implicity and is the fault of Mathematica or the machine, as opposed to the explicit rounding that Chop[] would implement.*)
];
];

Utemp=Utemp.(Subscript[T, p,q])[ω,ϕ,n]; (*This redefines "Utemp" by right multiplying it by the appropriate element. This ammounts to doing the zeroing off of the appropriate off-diagonal elements. We use the longer form of (Subscript[T, p,q]) that specififies n so that the "Optics" function remains completely self-contained and doesn't reference or change the global "Dimension" varaible.*)
optics[[p,q]]={ω,ϕ}
(*This stores the parameters that were found in "optics", indexed by p and q.*)

](*The end of an instance of the inner for loop amounts to having reduced "Utemp" to a smaller unitary by zeroing off an entire row and column (with exception to the diagonal element). The pth row has been explicitly zeroed off at this point (with exception to the diagonal element) and the pth column is automatically guaranteed to have been taken care of by unitarity.*)
];(*The end of the outer for loop amounts to having zeroed off all off-diagonal elements.*)

(*Next we find the appropriate "final" (they are actually on the input of each mode if you use our convention) phase shifts. This corresponds to the "D" at the end of the product in Equations (6) and (7) of Reck et. al.*)

For[i=1,i≤n,i++,
(*Loop through each mode*)
optics[[i,i]]=Arg[Diagonal[Utemp\[ConjugateTranspose]]][[i]]
(*The proper phase for mode i is the argument (phase angle) of the ith diagonal element of the conjugate transpose of the "Utemp" left over from all the zeroing. We store this phase angle in spot [[i,i]] of the "optics" array*)
];
optics
(*We return the parameters in the form of the array "optics". The diagonal element (i,i) corresponds to the phase shift on the mode. The element (p,q) for p>q is a pair of parameters describing the element at the intersection of mode p and q. The first parameter is the beam splitter angle, the second parameter is the phase shift on the output mode q. The elements (p,q) where p<q are meaningless and are left to be 0.*)
]

To restate the final code comment, the form of the output of Optics[] is an array/matrix as follows:

The (i,i) diagonal element corresponds to the phase shift on the input of mode i (in our convention). The (p,q) element for p>q is a pair of parameters describing the element at the intersection of mode p and q. The first parameter is the beam splitter angle, the second parameter is the phase shift on the output mode q. The (p,q) elements where p<q are meaningless and are left to be 0. Next we demonstrate this funtion by testing it on the example matrix from Reck et. al that we previously stored in “ex”.

Optics[ex\[ConjugateTranspose]]//N//MatrixForm (*As we noted before, to generate the parameters via the Reck et. al. make-extra-work-for-experimentalists-convention we have to conugate transpose our input.*)

Indeed this agrees with the parameters that are given in the paper. Note {-π/2, π} is an identity, i.e. no action on the intersection, and {0, 0} is a swap/exchange of the two modes. Also note that the first two zeroes on the diagonal represent 0 phase shifts on those modes, while the zeroes in the upper right triangle are just meaningless place keepers. You can use this information to implement the unitary with linear optics as follows:

Examine Figure 3 of Reck et. al. Take note that mode k (for any 1≤ k≤ N) is the mode that starts at input k, goes through every splitter, and ends at exit k’. To use our convention, construct this lattice using Reck et. al.’s labeled outputs as the inputs and vice versa.  So mode 1’s input is the upper right end of the figure, and it exits at the bottom of the figure. If you like the inputs on the left (so “time flows left to right”), just vertically reflect the figure. Now, at the intersection of mode p and q, place the beam splitter with angle specified by the first parameter of the (p,q) element of “Optics” and for the phase shifter just next to it use the phase specified by the second parameter of the (p,q) element of “Optics”.  There is actually an error here with the figure. If you note, the phase shifter after each beam splitter is acting on the larger number mode, whereas in the rest of the paper and in the algorithm this phase acts on the smaller number mode (see Figure 1 [Extra note on Figure 1: if you actually go through the math, you have to use 2ω for the sandwiched phase shifter of the Mach Zehnder if you want it to implement a beam splitter angle of ω.]). To correct for this, move those phase shifters to the line below, making them act on the outgoing smaller number mode (this is actually not an error in their convention if you don’t negate these phases in the inversion process). So looking at the figure, the phase shifters to the upper left of each beam splitter should really be to the lower left. For the phase shifter on the input of mode k (near the k’ label), use the phase specified by the (k,k) element of “Optics”.

Figure 3 of Reck et. al. is not explicitly wrong in their convention. You will be fine if you follow their convention (or run our algorithm by inputing the conjugate transpose of your unitary), and follow Figure 3 using the beam splitter angles given by the algorithm, the phase shift angles give by the algorithm interpreting each phase shifter as belonging to the intersection towards its bottom right, and negate the overall phase shifts before applying them to the outputs of each mode. However, on top of keeping track of the sign of the phases, you must also be careful in trying to directly correlate each optical element with a matrix from the decomposition since the decompsition is inverted which will cause the time ordering to reverse and each element to be conjugate transposed. It is easy to get misled or confused here. It’s likely that just reading about these conventions and figures has been confusing enough already. So for the sake of simplicity we now include our own figure that depicts how the optical parameters directly outputted by out algorithm correspond to a linear optical lattice.

Note that although for the sake of the graphic we use the likeness of a symmetric beam splitter, by default our algorithm gives the parameters for hermitian beam splitters.

Now we develop the inverse function of "Optics" which we call "Unitary". This function will take in a set of optical parameters of the same form that is outputted by "Optics" and output the unitary that those optical elements are implementing.

Unitary[parameters_]:= Module[{n,unitary,p,q},(*Localizes some variables*)

n=Dimensions[parameters][[1]]; (*The dimension of the space is the number of rows in "parameters"*)
unitary = IdentityMatrix[Dimension]; (*Initializes an identity matrix*)

For[p=n,p>1,p--,
For[q=p-1,q>0,q--,
(*Runs through all values of p,q from n,n-1 down to 2,1 which denote all the intersections of modes.*)

unitary = unitary.(Subscript[T, p,q])[ parameters[[p,q,1]],parameters[[p,q,2]] ,n]
(*Right multiplies "unitary" by the matrix corresponding to the optical element described by the {p,q} element of "parameters". As in "Optics", we use the longer form of (Subscript[T, p,q]) that specififies n so that the "Unitary" function remains completely self-contained and doesn't reference or change the global "Dimension" varaible. *)

]
]

unitary=unitary.Phases[Diagonal[parameters]]; 
(*Right multiplies "unitary" by the phases specified on the diagonal of "parameters".*)
unitary
(*Returns the unitary that is implemented by the optical elements described by "parameters".*)

]

We now demonstrate “Unitary” by showing that it is indeed the inverse of “Optics”. That is to say, if you use “Optics” to find the optical parameters for any unitary, and then use those parameters as the input to “Unitary”, you will just get your unitary back out. We show this for “ex” and the Pauli matrices as examples.

Unitary[Optics[ex]]//MatrixForm
Unitary[Optics[PauliMatrix[1]]]//MatrixForm
Unitary[Optics[PauliMatrix[2]]]//MatrixForm
Unitary[Optics[PauliMatrix[3]]]//MatrixForm

This concludes our coding and demonstration of Reck et. al.’s method for converting between arbitrary finite dimensional unitary operators and triangular linear optical lattices. These functions can now be used to serve greater purposes.

Multiphoton Calculations

The unitary operators discussed by Reck et. al. and the previous section describe the action of an optical network on a single photon in a number of modes. Aaronson and Arkhipov supply a homomorphism φ(U) that describes the action of such networks on multiple photon states in terms of the single photon unitary. Here we implement that action. To do so, we define a function that takes in a single photon unitary “U” and two lists of integers “S” and “T” that correspond to multiphoton basis states and outputs the value of ⟨S|φ(U)|T⟩, or equivalently, the corresponding matrix element of the larger unitary. The length of each list corresponds to the number of modes, and the elements correspond to the number of photons in each mode. Before defining this function, we first define the matrix (Subscript[U, S,T]) described by Aaronson and Arkhipov on pages 26-27.

(Subscript[U_List, S_List,T_List]):= Module[{i,j,k,UT,UST},

UT={};(*Prepares an empty list that will become (Subscript[U, T]).*)
UST={};(*Prepares an empty list that will become (Subscript[U, S,T]).*)

For[j=1,j≤Length[T],j++,(*Loops through the slots of T*)
For[k=1,k≤T[[j]],k++,(*Loops (Subscript[T, j]) times.*)
AppendTo[UT,Transpose[U][[j]]]
](*Here we append (Subscript[T, j]) copies of the jth column of U to (Subscript[U, T]). To do this, we append the jth row of the tranpose of U to "UT" (Subscript[T, j]) times via the for loop indexed by k. This will give us the transpose of (Subscript[U, T]) in "UT".*)
];

UT=Transpose[UT];(*We now transpose "UT" so that it now stores (Subscript[U, T]).*)


For[i=1,i≤Length[S],i++,(*Loops through the slots of S.*)
For[k=1,k≤S[[i]],k++,(*Loops (Subscript[S, i]) times.*)
AppendTo[UST,UT[[i]] ];
](*Here we append (Subscript[S, i]) copies of the ith row of (Subscript[U, T]) to (Subscript[U, S,T]).*)
];
UST(*Output (Subscript[U, S,T])*)
]

We now demonstrate this definition by implementing the example given by Aaronson and Arkhipov on page 27.

Module[{U,S,T},
U=(({{0, 1, 0}, {1, 0, 0}, {0, 0, -1}}));
S={0,1,2};
T={0,1,2};
(Subscript[U, S,T])//MatrixForm
]

Next we implement the Ryser Formula for calculating the permanent of a n×n matrix. http://mathworld.wolfram.com/RyserFormula.html. We do not include the Gray code optimization. The Combinatorica package contains a “GrayCodeSubsets” function, but implementing the Gray code optimization is not quite as simple as replacing “Subsets” in our code with this function.

Permanent[A_List]:=Module[{n,S,i,j,k},
n=Dimensions[A][[1]];(*n is the size of the matrix*)
S=Delete[Subsets[Range[n]],1];(*S holds the set of all non-empty subsets of the integers 1 through n*)
(((-1))^n) (∑_(k=1)^((2^n)-1))(((-1))^(Length[S[[k]]])) (∏_(i=1)^n)(∑_(j=1)^(Length[S[[k]]]))A[[i,S[[k,j]]]] (*The Ryser Formula*)
(*The first sum sums over all of the (2^n)-1 non-empty subsets S[[k]] (the kth element of the S list). Each term is multplied by -1 raised to the power of the cardinality of the corresponding subset S[[k]]. Then there is a product over all n of a sum over all of the elements of S[[k]]. These elements take the form S[[k,j]] indexed by j. Finally, inside this final sum we have the element of A in the ith row and S[[k,j]]th column.*)
]

We demonstrate our implementation of the permanent for general 2x2 and 3x3 matrices.

Clear[a,b,c,d,e,f,g,h,i]
Simplify[Permanent[{{a,b},{c,d}}]]
Simplify[Permanent[{{a,b,c},{d,e,f},{g,h,i}}]]

With (Subscript[U, S,T]) and the permanent defined, we can now calculate the matrix elements of the multiphoton unitary via Theorem 23 of Aaronson and Arkhipov.

⟨S_List\[VerticalSeparator]φ[U_List]\[VerticalSeparator]T_List⟩:=If[Total[S]==Total[T],((Permanent[(Subscript[U, S,T])])/(Sqrt[Times@@(S!) Times@@(T!)])),0]
(*Times @@ (S!) takes the cumulative product of the list composed of the factorial of each element in S. The matrix element is 0 whenever the two states don't have the same number of photons.*)

This concludes our implementation of the tools for working with multiple photons.

The Boson Sampling Problem

For the case where the input state T is a single photon per mode (and any remaining modes empty) then the calculation of (Subscript[U, S,T]) can be simplified by skipping straight to just taking (Subscript[S, i])copies of the ith row of U and then taking the first n columns of the resulting matrix, where n is the number of photons. We reimplement this code as a function (Subscript[U, S]) (or (Subscript[A, S]) equivalently). This should not be confused with the (Subscript[U, T]) mentioned in the (Subscript[U, S,T]) calculation, which corresponds to taking (Subscript[T, j]) copied of the jth column of U.

(Subscript[A_List, S_List]):= Module[{i,j,k,AS},

AS={};(*Prepares an empty list that will become (Subscript[A, S]).*)

For[i=1,i≤Length[S],i++,(*Loops through the slots of S.*)
For[k=1,k≤S[[i]],k++,(*Loops (Subscript[S, i]) times.*)
AppendTo[AS,A[[i]] ];
](*Here we append (Subscript[S, i]) copies of the ith row of A to (Subscript[A, S]).*)
];
Take[AS,All,Total[S]](*Takes the first n columns and outputs (Subscript[A, S])*)
]

Note that this defintion implicitly assumes that the number of input photons equals the number of output photons. With this optimized definition for the single photon per mode input case, we now implement the BosonSampling problem. This function takes in a matrix A (corresponding to the single photon unitary action of an optical arrangement, or just the first n columns of such a matrix) and a list S, and outputs the probability of detecting the S basis state given a single photon in each of the first n input modes. The code is equivalent to taking |⟨S|φ[A]|{1,1,1,1,...,0,0,0,...}⟩|^2 with the code we have already provided, but using this new function will optimize the running time by eliminating the many extraneous cases of multiplying by 1. The script D can be inserted by typing [ESC]scD[ESC]. This notation is used only to conform to that of Aaronson and Arkhipov. For convenience we define a mode coding friendly notation as well.

(Underscript[Pr, (Subscript[𝒟, A_List])])[S_List]:=(((Abs[Permanent[(Subscript[A, S])]])^2)/(Times@@(S!)))
Pr[A_List,S_List]:=(Underscript[Pr, (Subscript[𝒟, A])])[S]

We now demonstrate this function by showing the Hong Ou Mandel dip.

Module[{optics,U,S},
optics=(({{0, 0}, {{π/4,0}, 0}}));
U=Unitary[optics];
{(Underscript[Pr, (Subscript[𝒟, U])])[{2,0}],(Underscript[Pr, (Subscript[𝒟, U])])[{1,1}],(Underscript[Pr, (Subscript[𝒟, U])])[{0,2}]}
]

The core of the BosonSampling problem is sampling this probability distribution over Haar random unitary matrices. Accordingly, we now develop functions for generating random optical arrangements and Haar unitaries. First we define a function called “RandomOptics” that takes in an integer “n” specifying a number of modes and outputs random optical parameters for an n mode lattice. Each parameter is an angle, uniformly picked from [0, 2π]. One should note that the unitaries implemented by the optics generated in this matter are not Haar random. This is analogous to the fact that uniformly picking 2 random angles as spherical coordinates will not give a uniform distribution of points on the surface of a sphere.

RandomOptics[n_Integer]:=Module[{optics,p,q,i},(*Localizes some variables*)


optics=Array[0&,{n,n}]; (*Initializes an nxn matrix of all 0's to hole the optics parameters.*)
For[p=n,p>1,p--,
For[q=p-1,q>0,q--,
(*Runs through all values of p,q from n,n-1 down to 2,1 which denote all the intersections of modes.*)

optics[[p,q]]={RandomReal[2π],RandomReal[2π]}
(*Generates and stores random beamsplitter and phaseshifter angles for each intersection of modes.*)
]
];


For[i=1,i≤n,i++,
(*Loop through each mode*)
optics[[i,i]]=RandomReal[2π]
(*Generates random phaseshifts for each mode.*)
];
optics (*Outputs the generated parameters.*)
]

Now we demonstrate our “RandomOptics” function by generating the parameters for a random 5 mode lattice, then showing the transformation it implements, and verifying that the transformation is indeed unitary (up to rounding error) by multiplying it by its conjugate transpose.

Module[{optics,U},

optics=RandomOptics[5];
U=Unitary[optics];

Print[optics//MatrixForm];
Print[U //MatrixForm];
Print[Chop[U.U\[ConjugateTranspose]]//MatrixForm];
]

Next we define a function called “HaarUnitary” that takes in an integer “n” specifying a dimension and outputs an n×n Haar random unitary.

HaarUnitary[n_Integer]:=Module[{unitary={},i,j,v},

For[i=1,i≤n,i++,(*Iterates once for each row*)

v=RandomVariate[NormalDistribution[],n]+I*RandomVariate[NormalDistribution[],n];
(*Creates row vector of gaussian complex numbers. This vector will later be normalized resulting in a complex vector uniformly chosen from the surface of the complex n-ball, due to the fact that a product of like gaussians is radially gaussian.*)

If[i>1,
For[j=1,j≤i-1,j++,
v=v-Conjugate[unitary[[j]]].v*unitary[[j]](*This projects the gaussian vector onto the subspace orthogonal to all of the previously generated rows of the unitary. These previous rows are already normalized. It is theoretically possible that the gaussian vector will have no projection in the orthogonal subspace, in which case this algorithm fails and at least one row of the output unitary will be the zero vector. However, this is an exceedingly unlikely event, which can be verified by generating many random matrices, multiplying them by their conjugate transpose, and observing that you always (in practice) get the identity.*)
]
];

AppendTo[unitary,Normalize[v]](*Append the normalized version of the final vector to the unitary.*)
];
unitary(*Output the unitary*)
]

We now demonstrate “HaarUnitary” by showing a randomly generated 7×7 matrix, and verifying that it is unitary (up to round error) by multiplying it by its conjugate transpose.

Module[{U},
U=HaarUnitary[7];
Print[U//MatrixForm];
Print[Chop[U.U\[ConjugateTranspose]]//MatrixForm];
]

This concludes our implementation of the BosonSampling problem.

Distinguishable Photons

In an ideal world, the BosonSampling problem would be experimentally implemented with n completely indistinguishable photons. In this case, the probabilities predicted by quantum mechanics are given by the functions defined above. However, in the real world, a source of large numbers of completely indistinguishable photons is unrealistic and beyond current technology. As a result, real multiphoton sources can allow certain pairs of photons to be partiallty distinguished by ancillary modes such as polartization. Accordingly, we now develop functions that calculate the probabilties of detecting basis states under the assumption that the photons are completely distinguishable. Although this doesn’t address the question of partial distinguishability directly, it may be insightful to compare this case to the indistinguishable case and observe how the two differ.

In the case of completely distinguishable photons, no interference occurs. Thus, we can treat this problem by pretending that each photon enters and leaves the apparatus one by one, at separate times. In which case, the calculations of concern only depend on the element-wise absolute square of the unitary describing the linear optical lattice, which describes the probability of a single photon in a given input mode to end up in a given output mode. Noting that |(Subscript[U, p,q])((|)^2) is the probabilty for a photon staring in slot q to end up in slot p, the desired results can be obtained entirely from simple traditional counting arguments. Now we define a function “DistPr” that takes in a unitary “U”, an input basis state “T”, and an output basis state “S”, and returns the probability of detecting “S” given “T” as the input to the lattice described by “U”.

DistPr[U_List,T_List,S_List]:=If[Total[S]==Total[T],((Permanent[((Abs[(Subscript[U, S,T])])^2)])/(Times@@(S!))),0]
(*Both the Abs and traditional power functions are defined element-wise on matrices in Mathematica.*)

It is left for the reader to convince themself that this formula implements the correct counting argument. A rough sketch of the argument is as follows: First note that Perm(|U((|)^2)) describes the probability to find one photon in each output mode given that there was one photon in each input mode. The permanent reflects the fact that although we are treating the photons to be fundamentally distinguishable, the detectors still can’t distinguish between them, so we need to sum over every way that the individual distinguishable input photons could be arranged amongst the output modes. Now, note that the transformation of U into (Subscript[U, S,T]) discards the uneeded empty input and output modes, and adds new “pseudo-modes” or “slots” to account for each indivdual photon. Columns are duplicated to account for multiple photons in each mode as necessary, and rows are subsequently duplicated to treat multiple photons in output modes as if they were single photons in a number of separate “slots” that had the same probabilities of ending up in. This is done entirely to simplify the calculation. In effect, we have defined separate input and output “slots” for each photon, that way Perm(|(Subscript[U, S,T])((|)^2)) corresponds to the probabilty of finding one photon in each output slot given one photon in each input slot (do not confuse these abstract slots with the physical modes). Since in fact, we looked at (Subscript[s, i]) output slots that all correspond to the same physical output mode i, we have summed over certain probabilties too many times. Fixing for this amounts to dividing by (Subscript[s, i])! since that is the number of permutations of the fake slots we introduced. By doing this correction for each output mode, we arrive at the correct formula for the probability of concern. Again, this is only a rough sketch of the argument. We encourage the reader to carry through the calculation of a simple case to see for themselves why this gives the correct value.

Now we define another form of DistPr that only takes in a unitary and an ouput state. This function assumes that the input consisted of a single photon in the first n modes, where n is the number of photons. Note that this implicitly assumes that the number of input photons is equal to the number of output photons.

DistPr[U_List,S_List]:=((Permanent[((Abs[(Subscript[U, S])])^2)])/(Times@@(S!)))

It is interesting to note that this function is almost the same as the indistinguishable case, where the only difference is the order in which the permanent and absolute-square are taken. Perhaps unsurprisingly, this summarizes the “quantumness” of the BosonSampling problem, since approximating the permanent of a matrix with positive real entries can be done in polynomial time on a classical computer. This concludes out implementation of distinguishable photon probabilities.

Experimental Tools

In this section we develop tools to assist in the process of designing an apparatus and analyzing the results of verification experiments. We begin by defining a function “PrList” that takes in a unitary matrix and an integer n specifying the number of photons, and outputs a table that lists the absolute and relative probabilties of detecting all of the possible non-collision output Fock states (assuming that the input state consists of one photon in each of the first n modes and no photons in the remaining modes).

PrList[unitary_,n_]:=Module[
{m=Length[unitary],states,probs,relprobs,table}, (*Initializes working variables. m represents the number of modes.*)

states=Permutations[Join[Array[1&,n],Array[0&,m-n]]]; (*Creates a list of all the possible non-collision states.*)
probs=Pr[unitary,#]&/@states; (*Creates a new list that contains the probability of detecting each of the states in the "states" list.*)
relprobs=probs/Total[probs]; (*Divides the "probs" list by it's sum. This represents the probability of detecting the various non-collison states conditional on the output being non-collision.*)
table=Join[{#}&/@states,Transpose[{probs}],Transpose[{relprobs}],2]; (*Creates the output table. Some transposes and brackets are used to get the desired formatting.*)

TableForm[table,TableDepth->2,TableHeadings->{None,{"State","Probability","Relative Probability"}}](*Outputs the table. The "TableDepth -> 2" is needed to get the states to display correctly since they take the form of lists.*)

]

We now demostrate the PrList function using a Haar random 4 mode unitary and a 3 photon specification.

PrList[HaarUnitary[4],3]

Next we define the function RelativePermanentSignature which takes in a unitary and a number of photons, and displays a bar chart of the relative probabilities of all the non-collision output states (given the standard input state).

RelativePermanentSignature[unitary_,n_]:= Module[
{m=Length[unitary],states,probs},

states=Permutations[Join[Array[1&,n],Array[0&,m-n]]]; (*Generates a list of the non-collision output states.*)
probs=Pr[unitary,#]&/@states;(*Finds the probabilities of the states.*)

BarChart[{probs/Total[probs]}] (*Plots the list of probabilities, normalized by their sum.*)
]

We demonstrate this function for the case of a 4 mode Haar unitary with 3 photons.

RelativePermanentSignature[HaarUnitary[4],3]

Now we define another function called CompareSignatures which takes in the same input as RelativePermanentSignatures. CompareSignatures displays the same bar graph as RelativePermanentSignatures on the left, and on the right displays the relative probabilities of the same output states assuming that each photon were completey distinguishable.

CompareSignatures[unitary_,n_]:= Module[
{m=Length[unitary],states,probs,distprobs},

states=Permutations[Join[Array[1&,n],Array[0&,m-n]]];(*Generates a list of the non-collision output states.*)
probs=Pr[unitary,#]&/@states;(*Finds the ideal probabilities of the states.*)
distprobs=DistPr[unitary,#]&/@states;(*Finds the distinguishable probabilities of the states.*)

BarChart[{probs/Total[probs],distprobs/Total[distprobs]}] (*Plots both lists of probabilities, each normalized by their respective sum.*)
]

Again we demonstrate the case of a Haar random 4 mode unitary with 3 photons.

CompareSignatures[HaarUnitary[4],3]

Now we define a very useful interactive function called TweakOptics. This function takes in an optics matrix (of the same form of the output of the various optics functions previously defined) and a number of photons and displays a UI showing CompareSignatures for that setup. It also includes sliders for every optics component in the setup (beamsplitter and phaseshifter angles, not including initial or final phaseshifts which have no effect on the standard input) which can be varied, and the effect of which can be observed in real time in the CompareSignatures display. The optics matrix is also displayed for easy copy and paste once it has been tweaked to the user’s desire. Note however, that even though the number of sliders only grows as O((m^2)) with the number of modes, that’s still a lot to display!!! So this UI will get unwieldy quickly. Also, the “real-time” CompareSignatures will slow down considerably with large numbers of modes as it has to compute numerous permanents. Also note that the size of the CompareSignatures bar chart can be manually resized. There doesn’t seem to be an easy solution for displaying both the bar chart and optics matrix without Mathermatica reducing the default display size of the bar chart.

TweakOptics[opinit_,n_]:=With[{

op=
Table[

If[p>q,{(Subscript[ω, p,q]),If[p==Length[opinit],0,(Subscript[ϕ, p,q])]},0]
,{p,1,Length[opinit]},{q,1,Length[opinit]}](*Generates an optics matrix filled with place holder variables and 0's where appropriate.*)

,
controls=Sequence@@Flatten[Table[If[p==Length[opinit],{{{(Subscript[ω, p,q]),opinit[[p,q,1]]},-π,π},Delimiter},{{{(Subscript[ω, p,q]),opinit[[p,q,1]]},-π,π},{{(Subscript[ϕ, p,q]),opinit[[p,q,2]]},0,π},Delimiter}],{p,2,Length[opinit]},{q,1,p-1}],2]},(*Generates the sliders and sepearates each pair by a delimeter. The If statement prevents the addition of sliders for controlling the inital and final phaseshifters that have no effect.*)

Manipulate[Column[{CompareSignatures[Unitary[op],n],op//MatrixForm},Center],controls]]; (*Creates the interactive display of the CompareSignatures bar chart and optics matrix using the previously defined sliders as controls.*)

We now demonstrate TweakOptics for the 4 mode Haar unitary 3 photon case. Try it out!

TweakOptics[Optics[HaarUnitary[4]],3]

Work Space

This space is reserved as working room for the user’s use of the above developed tools. The alternative option would be to make these tools a Mathematica package which could be imported into the users’ notebooks, but that would decrease user friendliness and remove the possibility of the interactive demonstrations within the development.