Boson Sampling on a Photonic Chip

References

I Supplementary Information

II Materials and Methods

Multiphoton states generation. An $80$ MHz Ti-Sapphire oscillator outputs 100 fs pulses at 830 nm, which undergo type-I second harmonic generation in a $700$ $\mu$m $\beta$-BaB2O4 (BBO) crystal. This $415$ nm light is then used used to pump type-II collinear parametric downconversion (PDC) in $8$ mm-long AR-coated Potassium Dihydrogen Phosphate (KDP) crystals MosleyPJ2008hgu . Both modes of the resulting squeezed state are passed through spectral filters (Semrock, $\Delta\lambda{=}3$ nm), to maximize photon indistinguishability. In the three-photon experiment, two such PDC sources are coupled to four polarization-maintaining (PM) fibers Metcalf2012 . One of the photons is sent directly to a heralding detector, while the other three are coupled into the photonic circuit via a PM v-groove array (VGA) on a 6-axis alignment stage at the chip input. Programmable optical delay is provided by motor-controlled translation stages on each single-photon mode preceding the circuit, allowing the photons to arrive temporally coincident at the interferometric network in Fig 2A. Another PM VGA at the chip output couples the output modes into avalanche photodiode (APD) single photon counting modules (PerkinElmer SPCM-AQ4C), which are monitored by a home-built coincidence counting program loaded onto a commercially available FPGA development board (Xilinx SP605) operating with a 5 ns coincidence window. In the four-photon experiment, the higher order emission ($|22\rangle$) from a single PDC crystal was launched into two spatial modes of a different chip, though of identical geometry, which had VGAs glued to input and output. This reduced the coupling efficiency uncertainty for the four-photon QBSM, which reduced the size of the error bars in Fig. 3B. In both cases, to minimize undesired PDC emission terms the sources were pumped as lightly as possible while maintaining reasonable count rates, yielding $\lambda^{2}{=}0.011$ and $0.023$ in the three- and four-photon sampling experiments respectively. Photonic circuit fabrication. The boson sampling is performed on a UV-written silica-on-silicon integrated chip Metcalf2012 . The chip was fabricated by focusing a continuous wave UV laser (244 nm) onto a photosensitive layer to create a local increase in the index of refraction. The chip is then moved, via computer-controlled precision translation stages, transversely to the incident UV beam to trace out the desired waveguide network geometry. One can fabricate cross-coupling beamsplitters of a certain reflectivity by crossing waveguides at a specific angle. Such beamsplitters take up less space than traditional evanescent couplers in such circuits, allowing reduced chip size and lower loss Kundys2009 . While there are eight spatial modes in the middle of the circuit (Fig. 2A), the outer modes are not fabricated to the edge of the chip and are thus not accessible to the experimenter. They can then be accurately treated as losses on the neighboring modes, leaving a six (accessible) mode interferometric network. Boson sampling data collection. For the four-photon QBSM, the FPGA simultaneously counts all possible combinations of coincidence detection events (one detection event, two coincident detections etc) amongst the six APDs monitoring the output modes of the integrated circuit. For the three-photon QBSM, a similar set of coincidences is taken, but conditioned on detection of the herald photon by a separate APD. There is a nonzero background count rate on each detector which contributes to erroneous $N$-fold detection events when $N{-}1$ APDs detect photons and another APD erroneously fires. The $N{-}1$ coincidence rate is included in the above set of statistics collected by the FPGA and, by temporarily blocking all input modes to the circuit, one can estimate the background rate on each detector. The resulting background contribution to $N$-fold coincidences, comprising approximately $5\%$ of the total counts, are then subtracted from the raw data, with the results plotted in Fig. 3. Circuit characterization. Our photonic chip performs a non-unitary, complex-valued linear mapping of input to output modes $\Lambda_{ij}{=}\tau_{ij}e^{i\phi_{ij}}$. We follow the method outlined in Laing2012 to use a set of one- and two-photon data to reconstruct $\mathbf{\Lambda}$.

We can find $\tau_{i_{1},j_{1}}$ by coupling light (single photons) into mode $i_{1}$ and monitoring the probability of output in $j_{1}$ yielding

It is sufficient to measure the ratio of values between output ports, $\left(\tau_{i_{1},j_{1}}/\tau_{i_{1},j_{2}}\right)^{2}$, as this describes the transformation induced by our circuit up to a constant factor for each input, $x_{i}$. If $\mathbf{\Lambda}$ is the true transformation, then this procedure gives us $\mathbf{X}\mathbf{\Lambda}$ where $\mathbf{X}$ is a diagonal matrix with entries $X_{ii}{=}x_{i}$. However, the properties of the matrix permanent give $\textrm{Per}(\mathbf{X}\,\mathbf{\Lambda}){=}(\prod_{i}x_{i})\textrm{Per}(\mathbf{\Lambda})$. Therefore, Eq. 1 of the main text shows, since we always launch the same input state, every portion of our boson distribution is multiplied by a constant factor that cancels in the normalization. This data is collected by periodically pausing boson sampling data collection, blocking two of the three photon inputs, and measuring the relative power in each output mode. We thus obtain accurate values for $\tau$ as well as a variance that is important for calculating the error bars in Fig. 3.

One can then use two photon interference to find $\phi_{i,j}$ Laing2012 . If one inputs two photons into modes $i_{1},i_{2}$ and detect in modes $j_{1},j_{2}$, then we have $|S\rangle{=}\hat{a}^{\dagger}_{j_{1}}\hat{a}^{\dagger}_{j_{2}}|0\rangle$ and $|T\rangle{=}\hat{a}^{\dagger}_{i_{1}}\hat{a}^{\dagger}_{i_{2}}|0\rangle$. If single, indistinguishable photons are used, then the probability of post-selecting this output is

If the photons launched into the individual modes are distinguishable, then we get the incoherent sum of their individual statistics. We again take a matrix permanent, but as this is an incoherent process one finds

One can then perform a Hong-Ou-Mandel experiment and find that the resulting interference visibility is

We perform this two-photon interference experiment and fit a Gaussian to the resulting data to determine the visibility. Using the known $\tau_{i,j}$, one can then find $|\phi_{i_{1},j_{1}}+\phi_{i_{2},j_{2}}-\phi_{i_{1},j_{2}}-\phi_{i_{2},j_{1}}|$. Repeating this procedure for all accessible two photon dips, and applying additional constraints one can determine $\phi_{i,j}$ Laing2012 .

For our circuit geometry, we encounter an overconstrained problem as we measure more interference visibilities than unknown $\phi$. Therefore, we run a least squares minimization to find the set of $\phi$ that best fits our two photon interference data. To find the error bars in Fig. 3, we use a Monte Carlo method where the elements of $\mathbf{\Lambda}$ are selected from a normal distribution with an appropriate variance for each element. The one photon measurement was repeated periodically during the $160$ hour long boson sampling data collection, thus yielding the variance in the $\tau_{ij}$ characterization, which was determined to dominate the uncertainty in the predicted boson distribution in Fig. 3. This explains why the four-photon QBSM, which had VGAs glued to the ends and thus was much less susceptible to changes in the coupling, has significantly smaller error bars in the predicted distribution shown in Fig. 3B. The uncertainty in the Gaussian fit to the two-photon interference patterns and in $\tau_{ij}$ was then used to determine the variance in $\phi_{ij}$.

This characterization process is efficient, as a general linear transformation over $M$ modes can be described by $\mathcal{O}{(M^{2})}$ parameters, requiring $\mathcal{O}{(M^{2})}$ measurements with this technique. While photons from the PDC sources were used for the circuit characterization here, one can also use classical coherent states Laing2012 ; Rahimi-Keshari2012 . However, the single-photon based technique outlined here benefits from not requiring phase-stable path length matching and, because we use the same sources for boson sampling and characterization, the photonic degrees of freedom (polarization, spectrum etc) for the characterization match that used in the experiment.

With $\Lambda_{i,j}$ experimentally determined over the accessible modes, one can predict the post-selected boson distribution for any input/output $|S\rangle$, $|T\rangle$ by constructing $\mathbf{\Lambda}^{\mathbf{(S,T)}}$ from $\mathbf{\Lambda}$ and taking the permanent according to Eq. (1) in the main text. Thus, we effectively use the one and two photon boson distributions to characterize our non-unitary operation over the accessible modes. One can then take various matrix permanents of this non-unitary matrix to predict the boson distribution for any $N$.

III Supplementary Text

In this section we first outline how the boson distribution is given by a set of matrix permanents. We then show how the boson distribution can be accurately predicted by the permanents associated with a non-unitary matrix describing a lossy channel, $\mathbf{\Lambda}$, if one post-selects on trials where no photons are lost. Finally, the principal sources of error in this experiment, namely the photon distinguishability and higher order terms from our PDC photon sources, are discussed.

We assume $N$ bosons are injected into a network that performs a unitary transformation over $M$ modes. We consider the special case, appropriate to our experiment, where the input (and output) states contain no more than one boson per mode, though the general case is treated elsewhere Scheel2004 . Without loss of generality, let modes $1$ to $N$ contain an input boson, while modes $N{+}1$ to $M$ have vacuum inputs. The input state can then be described by

The unitary transformation allows one to evolve the operators according to

where $\hat{a}^{\dagger}_{i}$ and $\hat{b}^{\dagger}_{j}$ are creation operators on the $i$-th input and $j$-th output mode respectively. We then obtain the output state

To find the boson distribution, we project our output onto a state $|\mathbf{S}\rangle$ which, in the number state basis, we describe by an $N$ element vector $\mathbf{S}$, where $S_{j}$ gives the mode of the $j$-th boson. The probability of measuring this state is then

The term in square brackets can be expanded and includes $M^{N}$ terms, as one is selecting $N$ bosons from $M$ modes where repetitions are allowed ($>1$ boson in a mode). One can rewrite this term in square brackets to give

where $\mathbf{\widetilde{V}}$ is the set of $M^{N}$ permutations of $N$ photons amongst $M$ modes, repetitions allowed. The tilde notation will be used throughout this paper for a set of permutations. Then, $\widetilde{V}_{k}^{j}$ indicates the mode of the $k$-th boson in the $j$-th permutation. As an example, consider the case with $M=3$ modes and $N=2$ input bosons, then

Let us denote all $N!$ permutations of $\mathbf{S}$ by $\mathbf{\widetilde{S}}$ where $\widetilde{S}_{k}^{j}$ indicates the mode of the $k$-th boson in the $j$-th permutation. For example, if we project onto the state $|S\rangle{=}|011\rangle$ then $\mathbf{S}{=}$ and

It is clear we only retain terms from the summation in Eq. 10 where $\widetilde{V}^{j}\in\widetilde{S}$, otherwise at least one annihilation operator will act on vacuum and give $P_{S}{=}0$. This then leaves us with

The formula for the permanent of an $n\times n$ matrix $A$ with elements $a_{ij}$ is

where $\widetilde{\sigma}_{j}^{i}$ gives the $j$-th element of the $i$-th permutation of the numbers $1,2,....,n$. The term inside the modulus of Eq. 13 has the same form as the matrix permanent in Eq. 14. Our original unitary, $U$, can be described by an $M\times M$ matrix. However, it is obvious from Eq. 13 that, in general, we take the permanent of a subsection of $U$. Specifically, we only keep rows $1\rightarrow N$, those rows corresponding to modes with input photons. In addition, we only keep columns corresponding to the elements of $S$. Let us call this modified subsection of our original unitary $U^{(\mathbf{S,T})}$. Then, using the definition of the permanent we can rewrite Eq. 13

If one allows the possibility of more than one photon per input and output mode, then a similar analysis yields Eq. 1 in the main text Scheel2004 .

We also note that the above treatment assumes the bosonic commutation relation $[\hat{b}^{\dagger}_{i},\hat{b}^{\dagger}_{j}]=0\>\forall\>i,j$. If the system consisted of indistinguishable fermions, then the corresponding anticommutation relation $\{\hat{b}^{\dagger}_{i},\hat{b}^{\dagger}_{j}\}=0\>\forall\>i,j$ would be used, leading to alternating plus and minus signs introduced in the summation in Eq. 13, yielding the easily classically computable determinant.

III.2 Effects of loss

In any experimental implementation of boson sampling there will be losses. These losses, regardless of where they occur, can be modeled as beam splitters that link accessible to inaccessible modes Thomas-Peter2011a . When these losses are considered, it is important to ask whether the boson distribution is still given by a set of matrix permanents and, if so, what linear transformation does that matrix describe.

Let us adopt the convention that modes $1$ to $M$ are accessible modes while inaccessible loss modes are given the labels $M+1$ to $L$. There is an $L\times L$ unitary operation describing the evolution of our pure input state, though we must trace over these loss modes at the output, yielding a mixed state over the accessible modes. Let us assume that photons are input into modes $1$ to $N$ where $N\leq M$ and we post-select on cases where $N$ photons are detected, by definition, in the accessible modes $1$ to $M$. Then Eq. 10 becomes

but $S_{i}\leq M$ as we can only project on accessible modes. Therefore, even though $U$ is an $L\times L$ matrix and the elements of $\widetilde{V}$ range from $1\rightarrow L$, when we project onto $\widetilde{S}$ (all the permutations of $S$) we are left with

where $U^{\mathbf{(S,T)}}$ is again a modified version of the original unitary but only keeping rows $1\rightarrow N$ and columns in $S_{i}^{\prime}\leq M$. Since these elements always describe the accessible modes, then we can equivalently work in terms of $\Lambda$, where $\Lambda_{i,j}{=}U_{i,j}$ but $i,j\leq M$. In summary, when post-selecting on no bosons being lost, one can work in terms of $\Lambda$, a non-unitary linear transformation that is simply the subsection of $U$ over the accessible modes. The matrix permanents of such a non-unitary linear transformation lead to the theoretical predictions in Fig. 3 of the main text.

Equivalently, we can describe our system as a noisy (lossy) quantum channel in the operator sum representation Nielson04 . This formalism will be useful later when we discuss sources of error from higher order PDC terms. In this picture, the accessible modes are in the system $Q$, while all inaccessible loss modes form the environment system $E$. We can describe the transformation induced by our circuit over the full space $Q\otimes E$ by a unitary operation $U$. Let $\rho$ and $\sigma$ be the inputs to $Q$ and $E$ respectively, then the output in $Q$ after a projective measurement $P_{m}$ and tracing over the environment is described by

Let the basis for $E$ be described by $|e_{k}\rangle$ and the initial state of the environment be $\sigma=\sum_{j}q_{j}|j\rangle{\langle j|}$, then we can express Eq. 18 as

where $E_{jk}=\sqrt{q_{j}}{\langle e_{k}|}P_{m}U|j\rangle$ are the Kraus operators. We do not directly characterize $U$, as it extends over the environment which is inaccessible to the experimenter. However, with photons we can assume $\sigma=|0\rangle{\langle 0|}$, and all of our boson sampling results post-select on the case where no photons are lost to the environment. Therefore, the summation in Eq. 19 reduces to only one term with postselection

Experimental boson sampling efforts will sample a non-unitary transformation that is equivalent, when post-selecting on no bosons being lost, to the $E_{00}$ Kraus operator.

III.3 Sources of error

The computational difficulty for a classical machine to sample a boson distribution increases as the maximum error threshold is lowered. Therefore, while a QBSM need not sample the true boson distribution perfectly Aaronson , it will be easier to beat a classical machine if future QBSMs designs minimize their sampling errors. In this paper, we have benchmarked the accuracy of our QBSMs by inferring the probability distribution from our data, labeled $\rm{\mathbf{P}}^{\rm{exp}}$, and comparing it to the distribution obtained from Eq. 1, labeled $\rm{\mathbf{P}}^{\rm{th}}$. Throughout the text, we quantify the distance between two probability distributions via the $L_{1}$ distance, $d(\rm{\mathbf{P}}^{(1)},\rm{\mathbf{P}}^{(2)}){=}\frac{1}{2}\sum_{i}|\rm{P}^{(1)}_{i}-\rm{P}^{(2)}_{i}|$.

Our method of benchmarking QBSM accuracy will always yield a nonzero $d$ due to the finite number of collected samples. We perform a Monte Carlo simulation of $\mathbf{P}^{\rm{exp}}$ from a QBSM that perfectly samples $\mathbf{P}^{\rm{th}}$ as a function of the number of counts collected (Fig. 4, A and B), to show the rate at which $d$ asymptotically approaches zero as the number of samples collected increases. In the inset histograms, we show the range of outputs for this ideal QBSM for the actual number of experimental counts collected in the three and four photon cases, while the red dots indicate $d(\mathbf{P}^{\rm{th}},\mathbf{P}^{\rm{exp}})$. These two probability distributions show close agreement in Fig. 3, however a comparison of the histogram and red dots in Fig. 4, A and B indicates there is an additional source of error beyond the finite number of samples.

Due to experimental limitations, occasionally we sample distributions other than $\mathbf{P}^{\rm{th}}$. We model the effect of two such imperfections, photon distinguishability and higher order terms from our PDC sources, and form a new distribution $\mathbf{P}^{\rm{mod}}$ that accounts for these effects. We ignore the effects of photon impurity, as our post-selected data collection and use of nearly-spectrally factorable photon sources MosleyPJ2008hgu , minimizes this contribution. We find that the new $d(\mathbf{P}^{\rm{exp}},\mathbf{P}^{\rm{mod}})$, indicated with the green dot in Fig. 4, A and B, comes well within the output variance of an ideal machine. This indicates we have correctly diagnosed and modeled the principal sources of experimental error, which will be important in guiding designs of future, larger $N$, QBSMs.

The parametric downconversion sources we use to generate our photons actually generate a two-mode squeezed state that is given by

where $0\leq\lambda<1$ is the squeezing parameter whose magnitude is determined, in part, by the type of crystal and pump power used. For our experiment we wish to minimize higher order terms ($|22\rangle$ and $|33\rangle$ for the three- and four-photon QBSMs respectively) and so we deliberately lower our pump power as much as possible while maintaining a feasible count rate. For the three-photon experiment, $\lambda{=}\sqrt{0.011}$ and for the four-photon experiment $\lambda{=}\sqrt{0.023}$. However, even in this case we will sometimes inject more than $N$ photons into our circuit which, due to losses, could be observed as an $N$-fold detection at the output.

Due to the circuit characterization method employed, we have no information about what such terms will be. In the operator sum representation, our characterization only determines the $E_{00}$ Kraus operator, which describes the transformation of our input state when the environment, which includes all loss modes, starts and ends with zero photons. For example, if we instead inject five photons and lose two, then this process is described by the $E_{20}$ operator, about which we have no information.

To model the effect of using such squeezed sources, we start with a circuit with the same geometry used in the experiment. Non-uniform losses throughout the circuit can then be modeled by adding beam splitters that link the depicted accessible modes shown in Fig. 2A to inaccessible loss modes, where the beam splitter reflectivity indicates the loss in that channel Thomas-Peter2011a . Such ‘loss beamsplitters’ are added throughout the circuit. We cannot directly characterize these losses in our current circuit, however we can estimate them numerically. The characterized linear transformation $\mathbf{\Lambda}$ is a function of these losses as well as the fabricated interferometers shown in Fig. 2A. We have performed a loss-independent characterization of the interferometers Metcalf2012 , and then input this data into a genetic algorithm to find the relative losses between modes that best reproduces the characterized $\mathbf{\Lambda}$. We then apply three independent scaling factors to the relative losses at the sources, circuit and detectors. The source scaling factor is chosen to match the known source heralding efficiency, which is a measure of the loss in each source arm. The detector losses are scaled such that no detector has an efficiency greater than $50\%$, which is appropriate for APDs detecting photons at $830$ nm. Finally, we scale the relative losses in the circuit to reproduce the known overall system transmission observed experimentally.

With knowledge of these losses, we reproduce an actual unitary linear transformation $\mathbf{U}$ that extends over both the $N{=}6$ modes as well as all loss modes and can be used to accurately predict the effect of higher order PDC terms. Taking the three-photon QBSM as an example, we use Eq. 1 of the main text to find the probability distributions when the input is $|\mathbf{T}\rangle{=}|011010\rangle$ (desired single-photon input), as well as $|011020\rangle$ or $|022010\rangle$ (the first higher-order terms from our two sources), which we label $\mathbf{P}^{\rm{111}}$, $\mathbf{P}^{\rm{112}}$ and $\mathbf{P}^{\rm{221}}$ respectively. The boson distributions are then found by summing the terms where three photons appear in the desired accessible modes and zero, one, or two (respectively) photons appear in any combination of loss modes. For example, the probability of obtaining $|\mathbf{S}\rangle{=}|111000\rangle$ given input $|\mathbf{T}\rangle=|022010\rangle$ is the summation of all terms where three photons appear in the first three accessible modes and two photons appear in any combination of loss modes. The higher order probability distributions, $\mathbf{P}^{\rm{221}}$ and $\mathbf{P}^{\rm{112}}$, are weighted by $\lambda^{2}$ which is obtained via a conditional second order correlation measurement, $g^{(2)}(0)$SmithBJ2009ppg . As $\lambda^{2}$ is small, we only consider the first higher order terms from each source.

III.3.2 Photon Distinguishability

Boson sampling assumes indistinguishable bosons, while experimental implementations will always have some distinguishability. Assuming pure inputs we follow the notation of Rohde2012b to write our input state as

where $N$ is again the number of photons, $\alpha$ is a distinguishability parameter, and $A^{\dagger}_{\xi_{j},i}$ is the creation operator for photon $i$ in mode $\xi_{j}$. Each photon is in a superposition of a desired mode $\xi_{0}$ and another mode $\xi_{i}$. By analyzing the reduction in HOM dip visibility at a beamsplitter inside our circuit we find $\alpha=0.974$ on average in our experiment.

If one photon is distinguishable from the others, then the new probability distribution is given by the permanents of $N{-}1$ matrices which are incoherently summed. For example, assume an input state $|\mathbf{T}\rangle{=}|1\rangle^{\tau}|11000\rangle$ where $\tau$ labels a distinguishable photon, then for a unitary transformation $U$ the probability of obtaining an output $|\mathbf{S}\rangle$ is

where the terms in parentheses are permanents of $2\times 2$ matrices. We calculate these probability distributions when one photon is distinguishable and weight them by $|\alpha^{2}\sqrt{1-\alpha^{2}}|^{2}$ and $|\alpha^{3}\sqrt{1-\alpha^{2}}|^{2}$, the probability that one photon is distinguishable from the others for the three- and four-photon cases respectively. As $\alpha$ is large, we ignore the case when two photons are distinguishable.