Accurate and precise characterization of linear optical interferometers

Linear optics is important in quantum computation and communication. The simulation of a linear optical interferometer is computationally hard classically subject to reasonable conjectures [1]. Single-photon detectors and linear optical interferometers allow for efficient universal quantum computation via linear optical quantum computing (LOQC) [2]. Linear optics can simulate the quantum quincunx [3] and quantum random walks [4]. Linear optics coupled with laser-manipulated atomic ensembles enables long-distance quantum communication [5]. A wide class of communication protocols can be realized with coherent states and linear optics [6].

Recent advances in photonic technology including photonic circuits on silicon chips [7–11], noise-free high-efficiency photon number-resolving detectors [12–16], high-fidelity single-photon sources [17–20] have engendered the experimental implementation of multi-port linear optical interferometry. Reconfigurable linear optical interferometers that can perform arbitrary unitary transformations have been demonstrated [21–23].

The accurate and precise characterization of linear optics is important in quantum information processing tasks such as BosonSampling, LOQC and quantum walks. BosonSampling involves sampling from the output photon coincidence distribution of an interferometer on single photon inputs to each mode. Sampling from this distribution is computationally hard classically but is easy with a linear-optical interferometer. The classical hardness of the BosonSampling problem crucially depends on the error in the linear optical interferometer [24]. Similarly, the practical applications of BosonSampling, in quantum metrology and in the computation of molecular vibronic spectra, rely on the accurate implementation and characterization of linear optics [25, 26].

Accurate and precise characterization is important in LOQC because a high success probability of the employed non-deterministic linear-optical gates relies on implementing the desired gates with high fidelity [27]. Furthermore, linear interferometers used in photonic quantum walks, which display strong non-classical correlations, require accurate characterization especially if quantum walks are employed for solving classically hard problems [28–30]. In other words, that accurate and precise characterization of interferometers enables a verifiable quantum speedup of linear-optical protocols over classical computers.

Classical-light procedures [31, 32] for linear optics characterization are unsuitable for Fock-state based experiments because the interferometer parameters change during the coupling and decoupling of classical light sources and of homodyne detectors at the interferometer ports. This change could result from drift of interferometer parameters in the time required to couple sources and detectors or as the result of mechanical process of coupling itself. Characterization procedures that rely on Fock-state (rather than classical-light) inputs are thus more desirable in bosonsampling and LOQC implementations; such procedures would enable interferometer characterization without altering the experimental setup and would thus be accurate.

The Laing–O'Brien procedure [33] uses one- and two-photons for characterizing linear optical interferometers and is stable to the length scale of a photon packet. This procedure assumes perfect matching in source field and large-number statistics on the detected photons. Hence, implementations of this procedure are inaccurate due to spatiotemporal and polarization mode mismatch in the source field and imprecise due to shot noise.

We aim to devise an accurate and precise procedure that uses one- and two-photons for the characterization of linear optical interferometers and to devise a rigorous method to estimate the standard deviation in the linear optical interferometer parameters [34]. Furthermore, we aim to provide a correct alternative to the ${\chi }^{2}$-test, which has been used to estimate the confidence in the characterized interferometer parameters in current BosonSampling implementations [11, 35, 36] ⁵ .

Here we devise a procedure to characterize a linear optical interferometer accurately and precisely using one- and two-photon interference. Four strengths of our approach over the Laing–O'Brien procedure [33] are that our procedure (i) accounts for and corrects systematic error from spatial and polarization source-field mode mismatch via a calibration procedure (section 3); (ii) increases accuracy and precision by fitting experimental coincidence data to curve simulated using measured source spectra (section 3); (iii) accurately estimates the error bars on the characterized interferometer parameters via a bootstrapping procedure (section 4); and (iv) reduces the experimental time required to characterize interferometers using a scattershot procedure (section 5).

This section provides the background for our one- and two-photon characterization procedure. The action of a multi-port linear optical interferometer on single photons entering one or two input ports and vacuum entering the other ports is detailed. Specifically, we calculate the probability of detecting a photon at a given output port when a single photon is incident at a given input port. The section concludes with expressions for the probability of detecting a coincidence measurement when two controllably delayed photons are incident on the interferometer.

2.1. Action of a linear optical interferometer

In this subsection, we define linear optical interferometers by their action on single photons. We parameterize the unitary transformation effected by an interferometer and present our treatment of losses and dephasing at the interferometer ports.

Consider a single photon entering the ith mode of an m-mode interferometer. The monochromatic photonic creation and annihilation operators acting on the ith and the jth ports obey the canonical commutation relation⁶

$$\begin{eqnarray}&&[{a}_{i}({\omega }_{1}),{a}_{j}^{\dagger }({\omega }_{2})]={\delta }_{{ij}}\delta ({\omega }_{1}-{\omega }_{2}){\mathbb{1}}\end{eqnarray} \tag{ 1 }$$

for positive real frequencies ${\omega }_{1},{\omega }_{2}$. The state of a single photon entering the ith mode is

$$\begin{eqnarray}&&{| 1\rangle }_{i}={\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}\omega {f}_{i}(\omega ){a}_{i}^{\dagger }(\omega )| 0\rangle ,\end{eqnarray} \tag{ 2 }$$

where ${f}_{i}(\omega )$ is the normalized square integrable spectral function, $| 0\rangle$ is the m-mode vacuum state. The state of two photons entering modes i and $j\ne i$ of the interferometer is

$$\begin{eqnarray}&&{| 11\rangle }_{{ij}}={\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{1}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{2}\;{f}_{i}({\omega }_{1}){f}_{j}({\omega }_{2}){a}_{i}^{\dagger }({\omega }_{1}){a}_{j}^{\dagger }({\omega }_{2})| 0\rangle \end{eqnarray} \tag{ 3 }$$

with exchange symmetry holding if ${f}_{i}(\omega )={f}_{j}(\omega )$. One- and two-photon states are transformed into superpositions of one- and of two-photon states respectively under the action of the linear interferometer.

We treat linear interferometers as unitary quantum channels acting on the state of the incoming light. The interferometer transforms the photonic creation and annihilation operators according to

$$\begin{eqnarray}&&{a}_{j}^{\dagger }(\omega )\to \displaystyle \sum _{i=1}^{m}{V}_{{ij}}(\omega ){a}_{i}^{\dagger }(\omega )\end{eqnarray} \tag{ 4 }$$

and its complex conjugate, where $V(\omega )$ is the transformation matrix of the interferometer. Photon-number conservation imposes unitarity

$$\begin{eqnarray}&&{V}^{\dagger }(\omega )V(\omega )={\mathbb{1}}\end{eqnarray} \tag{ 5 }$$

of the transformation matrix $V(\omega )$ for all real ω. In general, the elements $\{{V}_{{ij}}(\omega )\}$ of the transformation matrix depend on the frequency of transmitted light. We assume that the spectral functions of the incoming light are narrow compared to frequencies over which the entries $\{{V}_{{ij}}\}$ change noticeably and thus treat V to be frequency-independent.

If only Fock states are incident at the interferometer and only photon-number-counting detection is performed on the outgoing light, then the measurement outcomes are invariant under phase shifts at each input and output port. That is, interferometer $\hat{V}={D}_{1}{{VD}}_{2}^{\dagger }$ produces the same measurement outcome as V for any diagonal unitary matrices D₁ and D₂. Mathematically, if ${D}_{1},{D}_{2}$ are diagonal unitary matrices, then

$$\begin{eqnarray}&&V\sim \hat{V}\ \Longleftrightarrow \ \hat{V}={D}_{1}{{VD}}_{2}^{\dagger }\end{eqnarray} \tag{ 6 }$$

is an equivalence relation. Members of the same equivalence class defined by this equivalence relation produce the same number-counting measurement outcomes on receiving Fock-state inputs.

Each equivalence class can be represented by a unique matrix U whose first row and first column consist of real elements. The complex matrix entries of the class representative $U\sim V$ are

$$\begin{eqnarray}&&\{{U}_{{ij}}={t}_{{ij}}{{\rm{e}}}^{{\rm{i}}{\theta }_{{ij}}}\;:\;{{\rm{t}}}_{{ij}}\in {{\mathbb{R}}}^{+},{\theta }_{{ij}}\in (-\pi ,\pi ],\\ &&\quad {\theta }_{i1}\equiv 0,{\theta }_{1j}=0\;\forall i,j\in \{1,2,...,m\}\}.\end{eqnarray} \tag{ 7 }$$

The constraints ${\theta }_{i1}\equiv 0,{\theta }_{1i}\equiv 0\;\forall i\in \{1,2,...,m\}$ on the input and output phases of the transformation matrix are obeyed in the following parameterization of U

$$\begin{eqnarray}U & = & L\times A\times M,\\ & & \overset{{\rm{def}}}{=}\left(\begin{array}{cccc}1 & 0 & \cdots & 0\\ 0 & \sqrt{{\lambda }_{2}} & \cdots & 0\ \\ \vdots & \vdots & \ddots & 0\\ 0 & 0 & \cdots & \sqrt{{\lambda }_{m}}\end{array}\right)\left(\begin{array}{ccc}1 & \cdots & 1\\ 1 & \cdots & {\alpha }_{2m}{{\rm{e}}}^{{\rm{i}}{\theta }_{2m}}\\ \vdots & \ddots & \vdots \\ 1 & \cdots & {\alpha }_{{mm}}{{\rm{e}}}^{{\rm{i}}{\theta }_{{mm}}}\end{array}\right)\\ & & \times \left(\begin{array}{cccc}\sqrt{{\mu }_{1}} & 0 & \cdots & 0\\ 0 & \sqrt{{\mu }_{2}} & \cdots & 0\ \\ \vdots & \vdots & \ddots & 0\\ 0 & 0 & \cdots & \sqrt{{\mu }_{m}}\end{array}\right).\end{eqnarray} \tag{ 8 }$$

Thus, the values $\{{\lambda }_{i}\},\{{\alpha }_{{ij}}\},\{{\theta }_{{ij}}\},\{{\mu }_{j}\}$ completely parametrize the class representative matrix U.

The input and output ports of the interferometer are amenable to time-dependent linear loss and dephasing. We model losses using parameters ${\nu }_{j}$ and ${\kappa }_{i}$, which are the respective probabilities of transmission at the input mode j and output mode i. Dephasing is modelled using parameters ${\xi }_{j}$ and ${\phi }_{i}$, which are the arbitrary multiplicative phases at the input and output ports. Hence, the actual transformation effected by the interferometer is given by the matrix U^lossy, which has matrix elements

$$\begin{eqnarray}{U}_{{ij}}^{\mathrm{lossy}} & = & {{\rm{e}}}^{{\rm{i}}{\phi }_{i}}\sqrt{{\kappa }_{i}}\;{{\rm{U}}}_{{ij}}\sqrt{{\nu }_{j}}{{\rm{e}}}^{{\rm{i}}{\xi }_{j}}\\ & = & {{\rm{e}}}^{{\rm{i}}{\phi }_{i}}\sqrt{{\kappa }_{i}}\sqrt{{\lambda }_{i}}\;{\alpha }_{{ij}}{{\rm{e}}}^{{\rm{i}}{\theta }_{{ij}}}\sqrt{{\mu }_{j}}\sqrt{{\nu }_{j}}{{\rm{e}}}^{{\rm{i}}{\xi }_{j}}.\end{eqnarray} \tag{ 9 }$$

Figure 1 depicts the relation between the representative matrix U and the actual transformation ${U}_{{ij}}^{\mathrm{lossy}}$ that is effected by the interferometer.

Zoom In Zoom Out Reset image size

This completes our parameterization of the linear optical interferometer. Our characterization procedure employs one- and two-photon inputs to estimate the values of parameters $\{{\lambda }_{i}\},\{{\alpha }_{{ij}}\},\{{\theta }_{{ij}}\},\{{\mu }_{j}\}$ of (9). In the next subsection, we recall the expectation values of measurements performed on interferometer outputs when one- and two-photon states are incident at the input ports.

2.2. One- and two-photon inputs to linear optical interferometer

Our characterization procedure employs single-photon counting to estimate the amplitudes $\{{\alpha }_{{ij}}\}$ of the representative matrix U entries. The arguments $\{{\theta }_{{ij}}\}$ of U are estimated using two-photon coincidence counts. In this subsection, we give expressions for one- and two-photon transmission probabilities, which are employed in our characterization procedure (section 3).

We first consider the case of single-photon transmission. The interferometer transforms the single-photon input state (2) to the state at the output ports according to (9). A photon is detected at the ith output port with a probability

$$\begin{eqnarray}&&{P}_{{ij}}={| {U}_{{ij}}^{\mathrm{lossy}}| }^{2}={\kappa }_{i}{\lambda }_{i}{\alpha }_{{ij}}^{2}{\mu }_{j}{\nu }_{j}.\end{eqnarray} \tag{ 10 }$$

When a single-photon is incident on the jth input port.

Whereas the values of $\{{\alpha }_{{ij}}\}$ are estimated using single photon counting, $\{{\theta }_{{ij}}\}$ values are estimated using two-photon coincidence measurement. We now present probabilities of detecting two-photon coincidence at the interferometer outputs when controllably delayed pairs of photon are incident at the input ports.

If a controllably delayed photon pair is incident at input ports j and ${j}^{\prime }$, then the probability ${C}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\tau )$ of coincidence measurement at detectors placed at output ports i and ${i}^{\prime }$ is

$$\begin{eqnarray}&&{C}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\tau )=\;{\kappa }_{i}{\kappa }_{{i}^{\prime }}{\nu }_{j}{\nu }_{{j}^{\prime }}[({t}_{{ij}}^{2}{t}_{{i}^{\prime }{j}^{\prime }}^{2}+{t}_{{{ij}}^{\prime }}^{2}{t}_{{i}^{\prime }j}^{2})\\ &&\quad \times \displaystyle \int {\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{| {f}_{j}({\omega }_{1}){f}_{{j}^{\prime }}({\omega }_{2})| }^{2}\\ &&\quad +2{t}_{{ij}}{t}_{{{ij}}^{\prime }}{t}_{{i}^{\prime }j}{t}_{{i}^{\prime }{j}^{\prime }}\displaystyle \int {\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{f}_{j}({\omega }_{1}){f}_{{j}^{\prime }}({\omega }_{2}){f}_{j}({\omega }_{2}){f}_{{j}^{\prime }}({\omega }_{1})\\ &&\quad \times \mathrm{cos}({\omega }_{2}\tau -{\omega }_{1}\tau +{\theta }_{{ij}}-{\theta }_{{{ij}}^{\prime }}-{\theta }_{{i}^{\prime }j}+{\theta }_{{i}^{\prime }{j}^{\prime }})].\end{eqnarray} \tag{ 11 }$$

On substituting according to (7), we obtain [40]

$$\begin{eqnarray}&&{C}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\tau )=\;{\kappa }_{i}{\kappa }_{{i}^{\prime }}{\lambda }_{i}{\lambda }_{{i}^{\prime }}{\mu }_{j}{\mu }_{{j}^{\prime }}{\nu }_{j}{\nu }_{{j}^{\prime }}[({\alpha }_{{ij}}^{2}{\alpha }_{{i}^{\prime }{j}^{\prime }}^{2}+{\alpha }_{{{ij}}^{\prime }}^{2}{\alpha }_{{i}^{\prime }j}^{2})\\ &&\quad \times \displaystyle \int {\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{| {f}_{j}({\omega }_{1}){f}_{{j}^{\prime }}({\omega }_{2})| }^{2}\\ &&\quad +2\gamma {\alpha }_{{ij}}{\alpha }_{{{ij}}^{\prime }}{\alpha }_{{i}^{\prime }j}{\alpha }_{{i}^{\prime }{j}^{\prime }}\displaystyle \int {\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{f}_{j}({\omega }_{1}){f}_{{j}^{\prime }}({\omega }_{2}){f}_{j}({\omega }_{2}){f}_{{j}^{\prime }}({\omega }_{1})\\ &&\quad \times \mathrm{cos}({\omega }_{2}\tau -{\omega }_{1}\tau +{\theta }_{{ij}}-{\theta }_{{{ij}}^{\prime }}-{\theta }_{{i}^{\prime }j}+{\theta }_{{i}^{\prime }{j}^{\prime }})],\end{eqnarray} \tag{ 12 }$$

where τ is the time delay between the two photons, ${f}_{j}(\omega ),{f}_{{j}^{\prime }}(\omega )$ describe the spectrum of light just before it enters the detectors and γ is the mode-matching parameter, which we described in the remainder of this section.

Two-photon coincidence probabilities (12) depend on the mode matching in the source field. Spatial and polarization mode mismatch is quantified by the mode-matching parameter γ [40]. Perfectly indistinguishable light sources, such as light from a single-mode fibre, have relative mode matching $\gamma =1$ whereas $\gamma =0$ indicates that the sources are completely distinguishable. Figure 4 depicts how imperfect mode matching, i.e., $\gamma \lt 1$, alters the observed two-photon coincidence counts. Our calibration procedure estimates and accounts for imperfect mode matching, which is assumed to be constant over the runtime of the characterization experiment.

The calculation of the expected coincidence probabilities as a function of the time delay between the photons is detailed in algorithm 1. The next section describes how single-photon transmission probabilities (10) and two-photon coincidence probabilities (12) are used for characterizing the linear optical interferometer.

In this section, we describe our procedure to characterize linear optical interferometers. The outline of this section is as follows. Section 3.1 describes the experimental data required by our characterization procedure. This experimental data are processed by various algorithms to determine the transformation matrix (8). The algorithm to determine the amplitudes $\{{\alpha }_{{ij}}\}$ of the transformation-matrix elements is presented in section 3.2. In section 3.3, we describe the calibration of the source field by determining the mode-matching parameter  γ. The estimation of $\{{\theta }_{{ij}}\}$ using two-photon interference is detailed in section 3.4. Maximum-likelihood estimation is employed to find the unitary matrix U that best fits the calculated $\{{\alpha }_{{ij}}\},\{{\theta }_{{ij}}\}$ values and serves as the representative matrix (8). We discuss the calculation of the best-fit unitary representative matrix in section 3.5.

3.1. Experimental procedure and inputs to algorithms

Our characterization procedure relies on measuring (i) the spectral function f_j of the source light, (ii) single-photon detection counts, (iii) two-photon coincidence counts from a beam splitter and (iv) two-photon coincidence counts from the interferometer. The measurement data constitute the inputs to our algorithms, which then yield the representative matrix. Before presenting the algorithms, we detail the experimental procedure and the inputs received by the algorithm in this subsection.

We characterize the spectral function $f({\omega }_{i})$ of the incoming light for a discrete set ${\rm{\Omega }}=\{{\omega }_{1},{\omega }_{2},...,{\omega }_{k}\}$ of frequencies. The integer k of frequencies at which the spectral function is characterized is commonly equal to the ratio of the bandwidth to the frequency step of the characterization device. The characterized spectral function $f({\omega }_{i})$ is used to calculate the coincidence probabilities as detailed in algorithm 1.

Algorithm 1.  Coincidence: Calculates the expected coincidence rate for two-photon interference for a given 2 × 2 submatrix of an arbitrary SU $(m)$ transformation.

Input:

• $k,{\rm{\Omega }}=\{{\omega }_{1},{\omega }_{2},...{\omega }_{k-1},{\omega }_{k}\}\in {({{\mathbb{R}}}^{+})}^{k}$ ▹ Frequencies at which ${f}_{1},{f}_{2}$ are given.

• ${f}_{1},{f}_{2}:{\rm{\Omega }}\to {{\mathbb{R}}}^{+}$            ▹ measured spectra.

• ${\ell },T=\{{\tau }_{1},{\tau }_{2},...,{\tau }_{{\ell }}\}\in {({\mathbb{R}}\cup 0)}^{{\ell }}$           ▹ Time delay values.

• $A\leftarrow \{{\alpha }_{{ij}},{\alpha }_{{{ij}}^{\prime }},{\alpha }_{i{\rm{J}}},{\alpha }_{{i}^{\prime }{j}^{\prime }}\}\in {({{\mathbb{R}}}^{+}\cup 0)}^{4}$       ▹ Amplitudes of 2 × 2 submatrix of A (8).

• ${\rm{\Theta }}\leftarrow {\theta }_{{ij}},{\theta }_{{{ij}}^{\prime }},{\theta }_{{i}^{\prime }j},{\theta }_{{i}^{\prime }{j}^{\prime }}\in (-\pi ,\pi ]$      ▹Phases of 2 × 2 submatrix of A (8).

• $\gamma \in [0,1]$   ▹ Mode-matching parameter of photon source.

Output:

• $C:T\to {{\mathbb{R}}}^{+}$ ▹ Two-photon coincidence probabilities correct up to multiplicative factor.

1: procedure ${\rm{C}}{\rm{OINCIDENCE}}$ $(k,{\rm{\Omega }},{f}_{1},{f}_{2},{\ell },T,A,{\rm{\Theta }},\gamma )$

2: for τ in T do

3:    $C(\tau )\leftarrow {\rm{I}}{\rm{NTEGRATE}}[F(A,{\rm{\Phi }},{f}_{1},{f}_{2},\gamma ,{\omega }_{i},{\omega }_{j},\tau ),\{{\omega }_{a}\in {\rm{\Omega }},{\omega }_{b}\in {\rm{\Omega }}\}].$

▹Numerically integrate rhs of (12) over ${\omega }_{i},{\omega }_{j}$ with ${\kappa }_{i}={\kappa }_{{i}^{\prime }}={\nu }_{j}={\nu }_{{j}^{\prime }}=1$.

4: end for

5: return C

6: end procedure

The amplitudes $\{{\alpha }_{{ij}}\}$ are determined by impinging single photons at the interferometer and counting single-photon detections at the outputs. Single-photon counting is repeated multiple $(B\in {{\mathbb{Z}}}^{+})$ times in order to estimate the precision of the obtained $\{{\alpha }_{{ij}}\}$ values. Specifically, the number

$$\begin{eqnarray}&&{N}_{{{ijb}}_{j}}\in {{\mathbb{Z}}}^{+}\ :\;i,j\in \{1,...,m\},{b}_{j}\in \{1,...,B\}\end{eqnarray} \tag{ 13 }$$

of single-photon detection events are counted at all m output ports $\{i\}$ for single photons impinged at the jth input ports in the b_jth repetition. The counting is then performed for each of the input ports $j\in \{1,...,m\}$ of the interferometer. Algorithm 2 uses $\{{N}_{{{ijb}}_{j}},{b}_{j}\in \{1,...,B\}\}$ values to estimate ${\alpha }_{{ij}}$ and the standard deviation of the estimate. The experimental setup for $\{{\alpha }_{{ij}}\}$ measurement is depicted in figure 2.

Zoom In Zoom Out Reset image size

Arguments $\{{\theta }_{{ij}}\}$ are calculated by fitting curves of measured coincidence counts to curves calculated using measured spectra according to (12). Appendix B elucidates the inputs and outputs of the curve-fitting procedure, such as the Levenberg–Marquardt algorithm [41, 42], employed by our algorithms. Before calculating $\{{\theta }_{{ij}}\}$, we calibrate the source field for imperfect mode matching by measuring coincidence counts on a beam splitter of known reflectivity. Controllably delayed single-photon pairs are incident at the two input ports of the beam splitter and coincidence counting is performed on the light exiting from its two output ports. Algorithm 3 details the estimation of γ using coincidence counts ${C}^{\mathrm{cal}}(\tau )$ for time delay τ between the incoming photons.

The absolute values and the signs of the arguments $\{{\theta }_{{ij}}\in (-\pi ,\pi ]\}$ are calculated separately. To estimate the absolute values $\{| {\theta }_{{ij}}| \}$ of the arguments, pairs of single photons are incident at two input ports 1 and $j\in \{2,...,m\}$ and coincidence measurement is performed at two output ports 1 and $i\in \{2,...,m\}$. The choice of the input and output ports labelled by index 1 is arbitrary. The signs

$$\begin{eqnarray}\mathrm{sgn}{\theta }_{{ij}}\overset{{\rm{def}}}{=}\left\{\begin{array}{cc}-1 & \mathrm{if}\;{\theta }_{{ij}}\lt 0,\\ 0 & \mathrm{if}\;{\theta }_{{ij}}=0,\\ 1 & \mathrm{if}\;{\theta }_{{ij}}\gt 0\end{array}\right.\end{eqnarray} \tag{ 14 }$$

of the arguments are estimated using an additional ${(m-1)}^{2}$ coincidence measurements. Algorithm 6 details the choice of input and output ports for estimating $\{\mathrm{sgn}{\theta }_{{ij}}\}$. A schematic diagram of the experimental setup for $\{{\theta }_{{ij}}\}$ estimation is presented in figure 3.

Zoom In Zoom Out Reset image size

3.2. Single-photon transmission counts to estimate $\{{\alpha }_{{ij}}\}$ (algorithm 2)

Now we present our procedure to estimate $\{{\alpha }_{{ij}}\}$ values using single-photon counting. Single-photon transmission probabilities are connected to the amplitudes $\{{\alpha }_{{ij}}\}$ according to the relation ${P}_{{ij}}={\kappa }_{i}{\lambda }_{i}{\alpha }_{{ij}}^{2}{\mu }_{j}{\nu }_{j}$ (10). Although the $\{{\alpha }_{{ij}}\}$ values can be calculated from single-photon transmission counts, the factors $\{{\lambda }_{i}\},\{{\mu }_{j}\}$ cannot. The transmission probabilities depend on the products of the factors $\{{\lambda }_{i}\},\{{\mu }_{j}\}$ and the loss terms $\{{\kappa }_{i}\},\{{\nu }_{j}\}$, so $\{{\lambda }_{i}\},\{{\mu }_{j}\}$ cannot be measured without prior knowledge of the losses. The loss terms are usually unknown and can change between experiments. Hence, we calculate the values of $\{{\alpha }_{{ij}}\}$ from single-photon measurements and choose $\{{\lambda }_{i}\}$ and $\{{\mu }_{j}\}$ such that U = LAM is unitary.

Zoom In Zoom Out Reset image size

The amplitudes $\{{\alpha }_{{ij}}\}$ are determined by estimating transmission probabilities. The probabilities ${P}_{11},{P}_{i1},{P}_{1j},{P}_{{ij}}$ of single-photon detection at output ports $1,i$ when single photons are incident at input ports $1,j$ are expresses in terms of the ${\alpha }_{{ij}}$ values according to

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{11}{P}_{{ij}}}{{P}_{1j}{P}_{i1}}=\displaystyle \frac{{| {r}_{1}{\lambda }_{1}{\alpha }_{11}{\mu }_{1}{s}_{1}| }^{2}}{{| {r}_{1}{\lambda }_{1}{\alpha }_{1j}{\mu }_{j}{s}_{j}| }^{2}}\displaystyle \frac{{| {r}_{i}{\lambda }_{i}{\alpha }_{{ij}}{\mu }_{j}{s}_{j}| }^{2}}{{| {r}_{i}{\lambda }_{i}{\alpha }_{i1}{\mu }_{1}{s}_{1}| }^{2}}={\left|\displaystyle \frac{{\alpha }_{11}{\alpha }_{{ij}}}{{\alpha }_{1j}{\alpha }_{i1}}\right|}^{2}.\end{eqnarray} \tag{ 15 }$$

The probabilities ${P}_{11},{P}_{i1},{P}_{1j},{P}_{{ij}}$ are estimated by counting transmitted photons. The definition (8) of ${\alpha }_{{ij}}$ implies that ${\alpha }_{11}={\alpha }_{i1}={\alpha }_{1j}=1$. Hence, the values of ${\alpha }_{{ij}}$ are connected to the single-photon transmission probabilities according to

$$\begin{eqnarray}&&{\alpha }_{{ij}}=\sqrt{\displaystyle \frac{{P}_{11}{P}_{{ij}}}{{P}_{1j}{P}_{i1}}},\end{eqnarray} \tag{ 16 }$$

which is independent of the losses at the input and the output ports.

The transmission probabilities P_ij are estimated by counting transmitted photons as follows. The estimated values of $\{{\alpha }_{{ij}}\}$ are random variables that are amenable to random error from under-sampling and experimental imperfections. Thus, data collection is repeated multiple times. For accurate estimation of ${\alpha }_{{ij}}$ and its standard deviation $\delta {\alpha }_{{ij}}$, the number B of repetitions is chosen such that the standard deviation of $\{{N}_{{{ijb}}_{j}}\ :{b}_{j}\in \{1,...,B\}\}$ converges in B for all $i,j\in \{1,...,m\}$. The mean and standard deviation of $\{{N}_{{{ijb}}_{j}}\ :{b}_{j}\in \{1,...,B\}\}$ converge for large enough B if the cumulants of the distribution are finite [43].

Algorithm 2.  AmplitudeEstimation: Uses single-photon detection counts to calculate the amplitudes of the complex entries of the transformation matrix. $\tilde{\bullet }$ represents our estimate of •.

Input:

• $m\in {{\mathbb{Z}}}^{+}$,       ▹ Number of modes of interferometer.

• ${N}_{{{ijb}}_{j}}\quad :\{1,...,m\}\times \{1,...,m\}\times \{1,...,B\}\to {{\mathbb{Z}}}^{+}$   ▹ Single-photon detection counts.

• $B\in {{\mathbb{Z}}}^{+}$       ▹ Number of times single-photon counting is repeated.

Output:

• $\{{\tilde{\alpha }}_{{ij}}\}\in {({{\mathbb{R}}}^{+}\cup 0)}^{{m}^{2}}$    ▹ Estimate of $\{{\alpha }_{{ij}}\}$ (8).

1: procedure ${\rm{A}}{\rm{MPLITUDE}}$ ${\rm{E}}{\rm{STIMATION}}$ $(m,{N}_{{{ijb}}_{j}},B)$

2: for $i,j\in \{1,...,m\}\times \{1,...,m\}$ do

3:   ${\tilde{\alpha }}_{{ij}}\;\leftarrow$ ${\rm{M}}{\rm{EAN}}$ $(\sqrt{{N}_{11{b}_{1}}{N}_{{{ijb}}_{j}}/{N}_{1{{jb}}_{j}}{N}_{i1{b}_{1}}}:{b}_{1},{b}_{j}\in \{1,...,B\})$

4: end for

5: return $\{{\tilde{\alpha }}_{{ij}}\}$

6: end procedure

The probabilities P_ij are estimated by counting single-photon detection events. Suppose ${N}_{{{ijb}}_{j}}$ photons are transmitted from input port j to the detector at output port i when ${N}_{{b}_{j}}$ photons are incident and ${b}_{j}\in \{1,...,B\}$ . For large enough B, the transmission probability converges according to

$$\begin{eqnarray}&&{P}_{{ij}}\leftarrow \mathrm{mean}\left\{\displaystyle \frac{{N}_{{{ijb}}_{j}}}{{N}_{{b}_{j}}}\ :\;{b}_{j}\in \{1,...,B\}\right\}.\end{eqnarray} \tag{ 17 }$$

Likewise, the amplitudes $\{{\alpha }_{{ij}}\}$ are estimated by averaging the single-photon detection counts according to

$$\begin{eqnarray}&&{\alpha }_{{ij}}=\sqrt{\displaystyle \frac{{P}_{11}{P}_{{ij}}}{{P}_{1j}{P}_{i1}}}\leftarrow \mathrm{mean}\\ &&\quad \times \left\{\sqrt{\displaystyle \frac{{N}_{11{b}_{1}}}{{N}_{{b}_{1}}}\displaystyle \frac{{N}_{{{ijb}}_{j}}}{{N}_{{b}_{j}}}\displaystyle \frac{{N}_{{b}_{j}}}{{N}_{1{{jb}}_{j}}}\displaystyle \frac{{N}_{{b}_{1}}}{{N}_{i1{b}_{1}}}}\ :\;{b}_{1},{b}_{j}\in \{1,...,B\}\right\}\\ &&\quad =\mathrm{mean}\left\{\sqrt{\displaystyle \frac{{N}_{11{b}_{1}}{N}_{{{ijb}}_{j}}}{{N}_{1{{jb}}_{j}}{N}_{i1{b}_{1}}}}\ :\;{b}_{1},{b}_{j}\in \{1,...,B\}\right\}.\end{eqnarray} \tag{ 18 }$$

The estimate of ${\alpha }_{{ij}}$ relies on single-photon counts measured by impinging photons at the first input port repeatedly (repetition index ${b}_{1}\in \{1,...,B\}$) and independently at the jth input port (with repetitions labelled by a different index ${b}_{j}\in \{1,...,B\}$).

Henceforth, we represent our estimate of any parameter • by $\tilde{\bullet }$. The estimate ${\tilde{\alpha }}_{{ij}}$ calculated using (18) is independent of ${N}_{{b}_{j}}$ and thus resistant to variations in the incident-photon number ${N}_{{b}_{j}}$ over different input modes j and different repetitions b_j. Thus, our estimates $\{{\tilde{\alpha }}_{{ij}}\}$ are accurate in the realistic case of fluctuating light-source strength and coupling efficiencies.

Finally, the standard deviations $\sigma ({\tilde{\alpha }}_{{ij}})$ of our estimates are calculated according to

$$\begin{eqnarray}&&\sigma ({\tilde{\alpha }}_{{ij}})\leftarrow \mathrm{std}.\;\mathrm{dev}.\left(\sqrt{\displaystyle \frac{{N}_{11{b}_{1}}{N}_{{{ijb}}_{j}}}{{N}_{1{{jb}}_{j}}{N}_{i1{b}_{1}}}}:\;{b}_{1},{b}_{j}\in \{1,...,B\}\right),\end{eqnarray} \tag{ 19 }$$

which converges for a large enough B. In line with standard nomenclature, we refer to these standard deviations as error bars. Algorithm 2 details the estimation of $\{{\tilde{\alpha }}_{{ij}}\}$ and error bars on the obtained estimates.

3.3. Calibration to estimate mode-matching parameter γ (algorithm 3)

In this subsection, we describe the procedure to calibrate our light sources for imperfect mode matching. The mode-matching parameter γ is estimated using one- and two-photon interference on an arbitrary beam splitter. First, the reflectivity of the beam splitter is determined using single-photon counting [33]. Next, controllably delayed photon pairs are incident at the beam splitter inputs and coincidence counting is performed on the beam splitter output. We introduce a curve-fitting procedure to estimate the value of γ such that (12) best fits the measured coincidence counts.

The beam-splitter reflectivity, which is denoted by $\mathrm{cos}\vartheta$, is estimated as follows. A beam splitter of reflectivity $\mathrm{cos}\vartheta$ effects the 2 × 2 transformation

$$\begin{eqnarray}{U}_{\mathrm{bs}} & = & \left(\begin{array}{cc}\mathrm{cos}\vartheta & {\rm{i}}\;\mathrm{sin}\vartheta \\ {\rm{i}}\;\mathrm{sin}\vartheta & \mathrm{cos}\vartheta .\end{array}\right)\\ & = & \left(\begin{array}{cc}1 & 0\\ 0 & {\rm{i}}\end{array}\right)\left(\begin{array}{cc}1 & 0\\ 0 & \mathrm{tan}\vartheta \end{array}\right)\left(\begin{array}{cc}1 & 1\\ 1 & -{\mathrm{cot}}^{2}\vartheta \end{array}\right)\\ & & \times \left(\begin{array}{cc}\mathrm{cos}\vartheta & 0\\ 0 & \mathrm{sin}\vartheta \end{array}\right)\left(\begin{array}{cc}1 & 0\\ 0 & {\rm{i}}\end{array}\right),\end{eqnarray} \tag{ 20 }$$

which is in the form of (8) with ${\alpha }_{22}\overset{{\rm{def}}}{=}{\mathrm{cot}}^{2}\vartheta$. The value of ${\alpha }_{22}$ is estimated using single-photon counting as described in algorithm 2. The estimated beam-splitter reflectivity is

$$\begin{eqnarray}&&\mathrm{cos}\tilde{\vartheta }=\sqrt{\displaystyle \frac{{\alpha }_{22}}{1-{\alpha }_{22}}}.\end{eqnarray} \tag{ 21 }$$

The error bar on $\mathrm{cos}\tilde{\vartheta }$ is estimated by repeating the photon counting along the lines of algorithm 2.

Algorithm 3.  Calibration Calculates the mode-matching parameter γ of source-field using a beam splitter of known reflectivity.

Input:

• $k,{\rm{\Omega }}=\{{\omega }_{1},{\omega }_{2},...{\omega }_{k-1},{\omega }_{k}\}\in {({{\mathbb{R}}}^{+})}^{k}$ ▹ Frequencies at which ${f}_{1},{f}_{2}$ are given.

• ${f}_{1},{f}_{2}:{\rm{\Omega }}\to {{\mathbb{R}}}^{+}$   ▹ Given spectral functions.

• ${\ell },T=\{{\tau }_{1},{\tau }_{2},...,{\tau }_{{\ell }}\}\in {({\mathbb{R}}\cup 0)}^{{\ell }}$ ▹ Time delay values coincidence is measured at.

• ${C}^{\mathrm{cal}}:T\to {{\mathbb{R}}}^{+}$    ▹ Measured coincidence curve.

• $\vartheta \in (-\pi ,\pi ]$ ▹ $\mathrm{cos}\vartheta$ is reflectivity of calibrating beam splitter.

Output:

• $\tilde{\gamma }\in [0,1]$ ▹ Estimate of mode-matching parameter of photon source.

1: procedure ${\rm{C}}{\rm{ALIBRATION}}$ $(k,{\rm{\Omega }},{f}_{1},{f}_{2},{\ell },T,{C}^{\mathrm{cal}},\vartheta )$

2:   $A\leftarrow \{\mathrm{cos}\vartheta ,\mathrm{sin}\vartheta ,\mathrm{sin}\vartheta ,\mathrm{cos}\vartheta \}$ ▹ Beamsplitter of reflectivity R (20)

3:   ${\rm{\Phi }}\leftarrow \{0,\pi /2,\pi /2,0\}$ ▹ Beamsplitter of reflectivity R (20)

4:    $C(\tau ,\gamma )\overset{{\rm{def}}}{=}{\rm{C}}{\rm{OINCIDENCE}}({\rm{\Omega }},{f}_{1},{f}_{2},T,A,{\rm{\Phi }},\gamma )$ ▹ The quantities ${\rm{\Omega }},{f}_{1},{f}_{2},{R}^{\mathrm{cal}}$ are given. γ is unknown. ${\rm{C}}{\rm{OINCIDENCE}}({\rm{\Omega }},{f}_{1},{f}_{2},T,{R}^{\mathrm{cal}},\gamma )$ depends on γ and τ

5:   return $\tilde{\gamma }\leftarrow$ ${\rm{FIT}}$ $(C(\tau ,\gamma ),{C}^{\mathrm{cal}}(\tau ),1/{C}^{\mathrm{cal}}(\tau ),\mathrm{InitGuesses})$ ▹ Least-squares curve fitting to obtain the value of γ that minimizes $\frac{{\sum }_{\tau \in T}{| {C}^{\mathrm{cal}}(\tau )-C(\tau ,\gamma )| }^{2}}{{C}^{\mathrm{cal}}(\tau )}$. The argument $1/{C}^{\mathrm{cal}}(\tau )$ is the weight function [44] that accounts for experimental noise, which is assumed to be proportional to $\sqrt{C(\tau )}$. Ignore values of τ at which $C(\tau )=1$. Appendix B details the choice of initial guesses to the algorithm.

6: end procedure

Next we estimate γ using two-photon coincidence counting. Controllably delayed pairs of photons are incident at the two input ports of the beam splitter. Coincidence measurement is performed at the output ports for different values of time delay between the two photons. A curve-fitting algorithm is employed to find the best-fit value of γ, i.e., the value $\tilde{\gamma }$ that minimizes the squared sum of residues between the measured counts and the coincidence counts expected from (12) for the beam splitter matrix (20). Algorithm 3 details the calculations of $\tilde{\gamma }$, which is used to estimate $\{{\theta }_{{ij}}\}$ values accurately.

3.4. Two-photon interference to estimate $\{{\theta }_{{ij}}\}$ (algorithms 4–6)

In this subsection, we describe our procedure to estimate the arguments $\{{\theta }_{{ij}}\}$ of the representative matrix U (8). Our procedure requires the measurement of coincidence counts for $2{(m-1)}^{2}$ different choices of input and output ports. Of these measurements, ${(m-1)}^{2}$ are used to estimate the absolute values $\{| {\theta }_{{ij}}| \}$ of the arguments and the remaining ${(m-1)}^{2}$ are used to estimate the signs $\{\mathrm{sgn}{\theta }_{{ij}}\}$.

The absolute values $\{| {\theta }_{{ij}}| \}$ are estimated as follows. Single-photon pairs are incident at input ports 1 and j and coincidence measurements are performed at output ports 1 and i for $i,j\in \{2,...,m\}$. The state (3) of a photon pair is transformed under the action of the 2 × 2 submatrix

$$\begin{eqnarray}{U}_{i1j1}=\left(\begin{array}{cc}\sqrt{{\kappa }_{1}} & 0\\ 0 & \sqrt{{\kappa }_{i}}\end{array}\right)\left(\begin{array}{cc}\sqrt{{\lambda }_{1}} & 0\\ 0 & \sqrt{{\lambda }_{i}}\end{array}\right)\left(\begin{array}{cc}1 & 1\\ 1 & {\alpha }_{{ij}}{{\rm{e}}}^{{\rm{i}}{\theta }_{{ij}}}\end{array}\right)\\ \quad \times \left(\begin{array}{cc}\sqrt{{\mu }_{1}} & 0\\ 0 & \sqrt{{\mu }_{j}}\end{array}\right)\left(\begin{array}{cc}\sqrt{{\nu }_{1}} & 0\\ 0 & \sqrt{{\nu }_{j}}\end{array}\right)\end{eqnarray} \tag{ 22 }$$

of U labelled by the rows 1 and i and columns 1 and j. The probability of detecting a coincidence at the output ports $1,i$ is

$$\begin{eqnarray}&&{C}_{i1j1}(\tau )={\kappa }_{j}{\kappa }_{1}{\lambda }_{j}{\lambda }_{1}{\nu }_{i}{\nu }_{1}{\mu }_{i}{\mu }_{1}[\{{\alpha }_{{ij}}^{2}+1\}\\ &&\quad \times \displaystyle \int {\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{| {f}_{j}({\omega }_{1}){f}_{1}({\omega }_{2})| }^{2}\\ &&\quad +2\gamma {\alpha }_{{ij}}\displaystyle \int {\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{f}_{j}({\omega }_{1}){f}_{1}({\omega }_{2}){f}_{j}({\omega }_{2}){f}_{1}({\omega }_{1})\\ &&\quad \times \mathrm{cos}({\omega }_{2}\tau -{\omega }_{1}\tau +{\theta }_{{ij}})],\end{eqnarray} \tag{ 23 }$$

which is obtained by setting ${i}^{\prime }={j}^{\prime }=1$ in (12).

The measured coincidence counts are used to estimate the value of $| {\theta }_{{ij}}|$ as follows. The shape of the coincidence-versus-τ curve (23) depends on the values of ${\alpha }_{{ij}}$ and ${\theta }_{{ij}}$ . The shape does not depend on the parameters ${\kappa }_{1},{\kappa }_{i},{\lambda }_{1},{\lambda }_{i},{\mu }_{1},{\mu }_{j},{\nu }_{1},{\nu }_{j}$, which lead to a constant multiplicative factor to the coincidence expression. Furthermore, the shape is unchanged under the transformation ${\theta }_{{ij}}\to -{\theta }_{{ij}}$ for ${\theta }_{{ij}}\in (-\pi ,\pi ]$ if the spectral functions are identical. Hence, $| {\theta }_{{ij}}|$ can be estimated using the shape of the coincidence function (23) and the values $\{{\tilde{\alpha }}_{{ij}}\}$ estimated using algorithm 2. A curve-fitting algorithm estimates the value $| {\tilde{\theta }}_{{ij}}| \in [0,\pi ]$ that best fits the measured coincidence counts. The calculation of $\{| {\tilde{\theta }}_{{ij}}| \}$ is detailed in algorithm 4.

Algorithm 4.  Argument2Port: Calculates the unknown complex argument in the entries of a 2 × 2 transformation using a two-photon coincidence curve.

Input:

• $k,{\rm{\Omega }}=\{{\omega }_{1},{\omega }_{2},...{\omega }_{k-1},{\omega }_{k}\}\in {({{\mathbb{R}}}^{+})}^{k}$ ▹ ${f}_{1},{f}_{2}$ are measured at frequencies Ω.

• ${f}_{1},{f}_{2}:{\rm{\Omega }}\to {{\mathbb{R}}}^{+}$  ▹ measured spectra.

• ${\ell },T=\{{\tau }_{1},{\tau }_{2},...,{\tau }_{{\ell }}\}\in {({\mathbb{R}}\cup 0)}^{{\ell }}$ ▹ Time delay values coincidence is measured at.

• ${C}^{\mathrm{exp}}:T\to {{\mathbb{R}}}^{+}$   ▹ Measured coincidence curve.

• $A\leftarrow \{{\alpha }_{{ij}},{\alpha }_{{{ij}}^{\prime }},{\alpha }_{{i}^{\prime }j},{\alpha }_{{i}^{\prime }{j}^{\prime }}\}$ ▹ Complex amplitudes of 2 × 2 submatrix of A (8).

• ${\rm{\Theta }}\leftarrow \{{\theta }_{{{ij}}^{\prime }},{\theta }_{{i}^{\prime }j},{\theta }_{{i}^{\prime }{j}^{\prime }}\in (-\pi ,\pi ]\}$ ▹ Three complex arguments of submatrix.

• γ    ▹ Mode-matching parameter of photon source.

Output:

• $| {\tilde{\theta }}_{{ij}}|$ ▹ Estimated magnitude of the unknown complex argument.

1: procedure ${\rm{A}}{\rm{rgument}}$2${\rm{P}}{\rm{ort}}$ $(k,{\rm{\Omega }},{f}_{1},{f}_{2},{\ell },T,{C}_{\mathrm{exp}},A,{\rm{\Theta }},\gamma )$

2:   ${\rm{\Phi }}\overset{{\rm{def}}}{=}\{{\theta }_{{ij}},{\theta }_{{{ij}}^{\prime }},{\theta }_{{i}^{\prime }j},{\theta }_{{i}^{\prime }{j}^{\prime }}\}$ ▹ Set of three known phases and one unknown phase.

3:  $C(\tau ,{\theta }_{{ij}})\overset{{\rm{def}}}{=}{\rm{C}}{\rm{OINCIDENCE}}({\rm{\Omega }},{f}_{1},{f}_{2},T,A,{\rm{\Phi }},\gamma )$

4:  return ${\tilde{\theta }}_{{ij}}\leftarrow | {\rm{LM}}(C(\tau ,{\theta }_{{ij}}),{C}^{\mathrm{exp}}(\tau ),1/{C}^{\mathrm{exp}}(\tau ))|$ ▹ Use curve fitting to compute the ${\theta }_{{ij}}$ value that minimizes $\frac{{\sum }_{\tau \in T}{| {C}_{\mathrm{exp}}(\tau )-C(\tau ,\gamma )| }^{2}}{{C}_{\mathrm{exp}}(\tau )}.$

5: end procedure

Our procedure computes the signs by using an additional ${(m-1)}^{2}$ coincidence measurements. First we arbitrarily set ${\theta }_{22}$ as positive

$$\begin{eqnarray}&&\mathrm{sgn}{\theta }_{22}=1\end{eqnarray} \tag{ 24 }$$

because of the invariance⁷ of one- and two-photon statistics under complex conjugation $U\to {U}^{*}\;$[33]. The signs of the remaining arguments $\{{\theta }_{{ij}}\}$ are set using the coincidence counts between output ports $\{i,{i}^{\prime }\}$ when photon pairs are incident at input ports $\{j,{j}^{\prime }\}$ for a suitable choice of $\{{i}^{\prime },{j}^{\prime }\}$ as we describe below. The coincidence probability at the output ports $i,{i}^{\prime }$ is

$$\begin{eqnarray}&&{C}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\tau )=\;{\kappa }_{i}{\kappa }_{{i}^{\prime }}{\lambda }_{i}{\lambda }_{{i}^{\prime }}{\mu }_{j}{\mu }_{{j}^{\prime }}{\nu }_{j}{\nu }_{{j}^{\prime }}[({\alpha }_{{ij}}^{2}{\alpha }_{{i}^{\prime }{j}^{\prime }}^{2}+{\alpha }_{{{ij}}^{\prime }}^{2}{\alpha }_{{i}^{\prime }j}^{2})\\ &&\quad \times \displaystyle \int {\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{| {f}_{j}({\omega }_{1}){f}_{{j}^{\prime }}({\omega }_{2})| }^{2}\\ &&\quad +2\gamma {\alpha }_{{ij}}{\alpha }_{{{ij}}^{\prime }}{\alpha }_{{i}^{\prime }j}{\alpha }_{{i}^{\prime }{j}^{\prime }}\displaystyle \int {\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{f}_{j}({\omega }_{1}){f}_{{j}^{\prime }}({\omega }_{2}){f}_{j}({\omega }_{2}){f}_{{j}^{\prime }}({\omega }_{1})\\ &&\quad \times \mathrm{cos}({\omega }_{2}\tau -{\omega }_{1}\tau +{\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }})],\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}&&{\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}\overset{{\rm{def}}}{=}| {\theta }_{{i}^{\prime }{j}^{\prime }}-{\theta }_{{{ij}}^{\prime }}-{\theta }_{{i}^{\prime }j}+{\theta }_{{ij}}| \in [0,\pi ].\end{eqnarray} \tag{ 26 }$$

Curve fitting is employed to estimate the value of ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}$ that best fits the measured coincidence counts.

Algorithm 5.  SignCalc: Calculates the complex-phase sign of an element of the 2 × 2 submatrix of an interferometer transformation matrix.

Input:
• $\beta \equiv | {\theta }_{{i}^{\prime }{j}^{\prime }}-{\theta }_{{{ij}}^{\prime }}-{\theta }_{{i}^{\prime }j}+{\theta }_{{ij}}|$    ▹ As defined in (27).
• ${\theta }_{{i}^{\prime }{j}^{\prime }},{\theta }_{{{ij}}^{\prime }},{\theta }_{{i}^{\prime }j},| {\theta }_{{ij}}|$     ▹ Equations (27–28).
Output:
• $\mathrm{sgn}{\theta }_{{ij}}$  ▹ Sign of $\theta \in (-\pi ,\pi ]$ is defined in (14)
1: procedure ${\rm{S}}{\rm{IGN}}{\rm{C}}{\rm{ALC}}$ $(\beta ,{\theta }_{{i}^{\prime }{j}^{\prime }},{\theta }_{{{ij}}^{\prime }},{\theta }_{{i}^{\prime }j},| {\theta }_{{ij}}| )$
2: ${\beta }^{+}\leftarrow | {\theta }_{{i}^{\prime }{j}^{\prime }}-{\theta }_{{{ij}}^{\prime }}-{\theta }_{{i}^{\prime }j}+| {\theta }_{{ij}}| |$ ▹ If ${\theta }_{{gh}}\lt 0$, then $\beta ={\beta }^{-}$.
3: ${\beta }^{-}\leftarrow | {\theta }_{{i}^{\prime }{j}^{\prime }}-{\theta }_{{{ij}}^{\prime }}-{\theta }_{{i}^{\prime }j}-| {\theta }_{{ij}}| |$ ▹ If ${\theta }_{{gh}}\gt 0$, then $\beta ={\beta }^{+}$.
4: $\mathrm{sgn}\;{\theta }_{{ij}}\leftarrow \mathrm{sgn}| \beta -{\beta }^{-}| -| \beta -{\beta }^{+}|$
5: return $\mathrm{sgn}\;{\theta }_{{ij}}$
6: end procedure

The estimated value of ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}$ is employed by algorithm 5 to ascertain the sign of ${\theta }_{{ij}}$. Algorithm 5 relies on the identity

$$\begin{eqnarray}&&\mathrm{sgn}{\theta }_{{ij}}=\mathrm{sgn}(| {\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}-{\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}^{-}| -| {\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}-{\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}^{+}| ),\end{eqnarray} \tag{ 27 }$$

and on known values of

$$\begin{eqnarray}&&{\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\pm }\overset{{\rm{def}}}{=}| {\theta }_{{i}^{\prime }{j}^{\prime }}-{\theta }_{{{ij}}^{\prime }}-{\theta }_{{i}^{\prime }j}\pm | {\theta }_{{ij}}| | ,{\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\pm }\in [0,\pi ]\end{eqnarray} \tag{ 28 }$$

to ascertain the sign of ${\theta }_{{ij}}$. If the sign of ${\theta }_{{ij}}$ is positive, then ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}={\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}^{+}$ and (27) returns a positive $\mathrm{sgn}{\theta }_{{ij}}$. Otherwise, ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}={\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}^{-}$, in which case (27) gives a negative sign.

Algorithm 6 iteratively chooses indices $i,{i}^{\prime },j,{j}^{\prime }$ such that the signs of ${\theta }_{{{ij}}^{\prime }},{\theta }_{{i}^{\prime }j},{\theta }_{{i}^{\prime }{j}^{\prime }}$ have already been ascertained before ascertaining the sign of ${\theta }_{{ij}}$. In each iteration, the values of ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\pm }$ are calculated by substituting $| {\theta }_{{ij}}| ,| {\theta }_{{{ij}}^{\prime }}| ,| {\theta }_{{i}^{\prime }j}| ,| {\theta }_{{i}^{\prime }{j}^{\prime }}| ,\mathrm{sgn}{\theta }_{{{ij}}^{\prime }},\mathrm{sgn}{\theta }_{{i}^{\prime }j},\mathrm{sgn}{\theta }_{{i}^{\prime }{j}^{\prime }}$. The algorithm estimates ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}$ by curve fitting measured coincidence counts to (25). Algorithm 5 is ascertains the sign of ${\theta }_{{ij}}$ using the estimates of ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}$ and ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\pm }$. One suitable ordering of indices ${{ii}}^{\prime }{{jj}}^{\prime }$, which we depict in figure 5, is

• set ${i}^{\prime }=2,{j}^{\prime }=1$ to determine $\mathrm{sgn}{\theta }_{i2}$ for $i\in \{3,...,m\}$ (figure 5(b)),
• set ${i}^{\prime }=1,{j}^{\prime }=2$ to determine $\mathrm{sgn}{\theta }_{2j}$ for $j\in \{3,...,m\}$ (figure 5(c)),
• set ${i}^{\prime }=2,{j}^{\prime }=2$ to determine $\mathrm{sgn}{\theta }_{{ij}}$ for $(i,j)\in \{3,...,m\}\times \{3,...,m\}$ (figure 5(d)).

In summary, $\mathrm{sgn}{\theta }_{{ij}}$ is determined using the values of ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}$, which are estimated by curve fitting, and of ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\pm }$ , which are computed using the signs and amplitudes of ${\theta }_{{{ij}}^{\prime }},{\theta }_{{i}^{\prime }j},{\theta }_{{i}^{\prime }{j}^{\prime }}$. Algorithms 4–6 detail the step-by-step procedure to determine the absolute values and the signs of $\{{\theta }_{{ij}}\}$.

Zoom In Zoom Out Reset image size

For certain interferometers U, the ordering of indices ${{ii}}^{\prime }{{jj}}^{\prime }$ depicted in figure 5 can lead to instability in the characterization procedure. Appendix A elucidates on this instability and presents strategies to counter the instability. This completes our procedure to characterize the matrix A for representative matrix U = LAM. In the next subsection, we present a procedure to estimate the matrix that is most likely for the characterized matrix A.

Algorithm 6.  ArgumentCalc: Calculate $\{{\theta }_{{ij}}\}$ using two-photon coincidences.

Input:

• $k,{\rm{\Omega }}=\{{\omega }_{1},{\omega }_{2},...{\omega }_{k-1},{\omega }_{k}\}\in {({{\mathbb{R}}}^{+})}^{k}$     ▹ ${f}_{1},{f}_{2}$ are measured at frequencies Ω.

• ${f}_{1},{f}_{2}:{\rm{\Omega }}\to {{\mathbb{R}}}^{+}$    ▹ measured spectra.

• ${\ell },T=\{{\tau }_{1},{\tau }_{2},...,{\tau }_{{\ell }}\}\in {({\mathbb{R}}\cup 0)}^{{\ell }}$   ▹ Time delay values coincidence is measured at.

• ${C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{exp}}(\tau )$ for $(i,{i}^{\prime },j,{j}^{\prime })\in \{1,2\}\times \{1,...,m\}\times \{1,2\}\times \{1,...,m\},i\ne {i}^{\prime },j\ne {j}^{\prime }$

▹ Measured coincidence at output ports ${i}^{\prime },{j}^{\prime }$ when photons that have mutual delay τ are incident at input ports $i,j$.

• $\tilde{\alpha }:\{2,...,m\}\times \{2,...,m\}\to {{\mathbb{R}}}^{+}$    ▹ Complex amplitudes (8).

• $\gamma \in [0,1]$ ▹ Mode-matching parameter estimated using Algorithm 3.

Output:

• ${\tilde{\theta }}_{{ij}}:\{1,...,m\}\times \{1,...,m\}\to (-\pi ,\pi ]$  ▹ Complex Arguments (8).

1: procedure ${\rm{A}}{\rm{RGUMENT}}$ ${\rm{C}}{\rm{ALC}}$ $(k,{\rm{\Omega }},{f}_{1},{f}_{2},{\ell },T,{C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{exp}}(\tau ),{\alpha }_{{ij}},\gamma )$

2: for i in $\{1,...,m\}$ do

3: ${\theta }_{i1},{\theta }_{1i},\mathrm{sgn}{\theta }_{i1},\mathrm{sgn}{\theta }_{1i}\leftarrow 0$    ▹ The first row, column are real valued.

4: end for

5: for (i, j) in $\{2,...,m\}\times \{2,...,m\}$ do

6: $A\leftarrow \{1,1,1,{\alpha }_{{gh}}\}$, ${\rm{\Phi }}\leftarrow \{0,0,0\}$ ▹ 2 × 2 matrix: rows $1,i$, columns $1,j$.

7:   $| {\tilde{\theta }}_{{ij}}| \leftarrow {\rm{A}}{\rm{RGUMENT}}2{\rm{P}}{\rm{ORT}}\quad ({C}_{1i1j}^{\mathrm{exp}}T,{\rm{\Omega }},{f}_{1},{f}_{2},T,A,{\rm{\Phi }},\gamma )$

8:    end for

9:   $\mathrm{sgn}{\theta }_{22}\leftarrow 1$      ▹ The sign of ${\theta }_{22}$ is positive by definition.

10:  for ($i,{i}^{\prime },j,{j}^{\prime })\in \{2\}\times \{3,...,m\}\times \{2\}\times \{3,...,m\}\cup \{1\}\times \{2\}\times \{2\}\times \{3,...,m\}\cup \{2\}\times \{3,...,m\}\times \{1\}\times \{2\}$ do

11: $A\leftarrow \{0,0,0\},{\rm{\Phi }}\leftarrow \{0,0,0\}$

12: ${\beta }_{i,{i}^{\prime },j,{j}^{\prime }}\leftarrow {\rm{A}}{\rm{RGUMENT}}2{\rm{P}}{\rm{ORT}}({C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{exp}}(\tau ),{\rm{\Omega }},{f}_{1},{f}_{2},T,A,{\rm{\Phi }},\gamma )$

13: end for

14: for i in $\{3,...,m\}$ do

15:   ${\tilde{\theta }}_{i2}\leftarrow | {\tilde{\theta }}_{i2}|$ ${\rm{S}}{\rm{IGN}}{\rm{C}}{\rm{ALC}}$ $({\beta }_{122i},0,{\theta }_{22},0,| {\theta }_{i2}| ,));$

16:   ${\tilde{\theta }}_{2i}\leftarrow | {\tilde{\theta }}_{2i}|$ ${\rm{S}}{\rm{IGN}}{\rm{C}}{\rm{ALC}}$ $({\beta }_{2i12},0,{\theta }_{22},0,| {\theta }_{2i}| ,));$

17:    end for

18:    for (i, j) in $\{3,...,m\}\times \{3,...,m\}$ do

19:   ${\tilde{\theta }}_{{ij}}\leftarrow | {\tilde{\theta }}_{{ij}}|$ ${\rm{S}}{\rm{IGN}}{\rm{C}}{\rm{ALC}}$ $({\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }},{\theta }_{22},{\theta }_{i2},{\theta }_{2j},| {\theta }_{{ij}}| )$

20:      end for

21:      return $\{{\theta }_{{ij}}\}$

22: end procedure

3.5. Maximum-likelihood estimation for finding unitary matrix

At this stage, we have estimated the matrix A (8). The diagonal matrices L and M can be uniquely determined from A as follows. The representative matrix U = LAM is unitary so we have

$$\begin{eqnarray}&&{{UU}}^{\dagger }={\mathbb{1}},\end{eqnarray} \tag{ 29 }$$

which, upon substitution U = LAM, implies that

$$\begin{eqnarray}{{LAMM}}^{*}{A}^{\dagger }{L}^{*} & = & {\mathbb{1}}\\ \Longrightarrow {{AMM}}^{*}{A}^{\dagger } & = & {L}^{-1}{L}^{*-1}.\end{eqnarray} \tag{ 30 }$$

Considering the first columns of the matrices (30) gives

$$\begin{eqnarray}&&{A}_{{ij}}{M}_{{jj}}^{*}{M}_{{jj}}{A}_{j1}^{\dagger }=\left(\begin{array}{c}1\\ 0\\ \vdots \\ 0\end{array}\right)\end{eqnarray} \tag{ 31 }$$

or

$$\begin{eqnarray}&&A\left(\begin{array}{c}{\mu }_{1}\\ {\mu }_{2}\\ \vdots \\ {\mu }_{m}\end{array}\right)=\left(\begin{array}{c}1\\ 0\\ \vdots \\ 0\end{array}\right).\end{eqnarray} \tag{ 32 }$$

Similarly, using ${U}^{\dagger }U={\mathbb{1}}$ we obtain

$$\begin{eqnarray}&&{A}^{\dagger }\left(\begin{array}{c}1\\ {\lambda }_{2}\\ \vdots \\ {\lambda }_{m}\end{array}\right)=\displaystyle \frac{1}{{\mu }_{1}}\left(\begin{array}{c}1\\ 0\\ \vdots \\ 0\end{array}\right).\end{eqnarray} \tag{ 33 }$$

Equations (32) and (33) are systems of linear equations that can be solved for L and M respectively using standard methods [45]. The solutions L and M of the linear systems and the characterized matrix A give us the representative matrix U = LAM.

Algorithm 7.  MaxLikelyUnitary: Calculates unitary matrix that has maximum likelihood of generating estimated $\{{\alpha }_{{ij}}\},\{{\theta }_{{ij}}\}$.

Input:
• $\tilde{\alpha }:\{1,...,m\}\times \{1,...,m\}\to {{\mathbb{R}}}^{+}\cup 0$   ▹ Estimated amplitudes of A (8).
• $\tilde{\theta }:\{1,...,m\}\times \{1,...,m\}\to (-\pi ,\pi ]$          ▹ Estimated arguments of A (8).
Output:
• $W\in \mathrm{SU}(n)$    ▹ Unitary matrix with maximum likelihood of generating A.
1: procedure ${\rm{M}}{\rm{ax}}$ ${\rm{L}}{\rm{ikely}}$ ${\rm{U}}{\rm{nitary}}$ $({\alpha }_{{ij}},{\theta }_{{ij}})$
2:   ${\lambda }_{1}\leftarrow 1$
3:   $\{{\tilde{\mu }}_{i}:i\in \{1,...,m\}\}\;\leftarrow$ solution of system (32) of linear equations.
4:   $\{{\tilde{\lambda }}_{i}:i\in \{2,...,m\}\}\;\leftarrow$ solution of system (33) of linear equations.
5:   ${\tilde{U}}_{{ij}}\leftarrow {\tilde{\lambda }}_{i}{\tilde{\alpha }}_{{ij}}{{\rm{e}}}^{{\rm{i}}{\tilde{\theta }}_{{ij}}}{\tilde{\mu }}_{j}$
6:   $W\leftarrow {(\tilde{U}{\tilde{U}}^{\dagger })}^{-\frac{1}{2}}\tilde{U}$ ▹ Assumption: ${U}_{{ij}}-{\tilde{U}}_{{ij}}$ is an iid Gaussian random variable with zero mean for all $i,j\in \{1,...,m\}$.
7: end procedure

The experimentally determined $\tilde{A}$ is different from the actual A because of random and systematic error in the experiment, where we denote the experimentally determined values of interferometer parameter • by $\tilde{\bullet }$. Similarly, the $\tilde{L}$ and $\tilde{M}$ matrices obtained by solving equations (32) and (33) for $\tilde{A}$ (rather than A) differ from the actual L and M respectively. The estimated $\tilde{U}=\tilde{L}\tilde{A}\tilde{M}$ is thus a non-unitary matrix and is not equal to U in general. Furthermore, $\tilde{U}$ is a random matrix, which depends on the random errors in the one- and two-photon experimental data.

We employ maximum-likelihood estimation to calculate the unitary matrix W that best fits the collected data. First, bootstrapping techniques are used to estimate the probability-density function (pdf) of the entries of the random matrix $\tilde{U}$ [46, 47]. Next, standard methods in maximum-likelihood estimation [48] are employed to find the unitary matrix W. Maximum-likelihood estimation simplifies under the assumption that the error on $\tilde{U}$ is a Gaussian random matrix ensemble, i.e, that the matrix entries $\{{\tilde{U}}_{{ij}}\}$ are complex independent and identically distributed (iid) Gaussian random variables centred at the correct matrix entries. In this case, the most likely unitary matrix W is the one that minimizes the Frobenius distance⁸ from $\tilde{U}\;$[49]. The unitary matrix

$$\begin{eqnarray}&&W={(\tilde{U}{\tilde{U}}^{\dagger })}^{-\frac{1}{2}}\tilde{U},\end{eqnarray} \tag{ 36 }$$

minimizes the Frobenius-norm distance from $\tilde{U}\;$[50]. Thus, if the random errors $\{{U}_{{ij}}-{\tilde{U}}_{{ij}}\}$ in the matrix elements are iid Gaussian random variables with mean zero, then W is the best-fit unitary matrix. Figure 6 is a depiction of the actual, the estimated and the most likely transformation matrices. Algorithm 7 computes W.

Zoom In Zoom Out Reset image size

This completes our procedure to estimate the most-likely unitary matrix W that represents the linear optical interferometer. In the next section, we present a procedure to estimate the error bars on the entries of the estimated representative matrix W accurately.

In this section, we present a procedure to estimate the error bars on the matrix entries $\{{W}_{{ij}}\}$ of the characterized representative matrix W. The entries $\{{W}_{{ij}}\}$ computed by algorithms 1–7 are random variables because of random error in experiments. Obtaining accurate error bars on these random variables is important for using characterized linear optical interferometers in quantum computation and communication. Current procedures compute error bars under the assumption that Poissonian shot noise is the only source of error in experiment [21, 23].

We choose to employ bootstrapping on the data determine error bars [46, 47, 51–53]. Monte-Carlo simulation is widely used but this technique is not applicable here because the Poissonian shot noise assumption is not reliable given the presence of other sources of error some of which are not understood. Bootstrapping is preferred because the nature of the error need not be characterized and instead relies on random sampling with replacement from the measured data. Bootstrapping can be employed toyield estimators such as bias, variance and error bars.

Algorithm 8 calculates the error bars $\sigma ({W}_{{ij}})$ using estimates of the $\{{W}_{{ij}}\}$ pdf's, which are obtained using bootstrapping as follows. The algorithm simulates N characterization experiments using the one- and two-photon data, i.e., the inputs to algorithms 1–7. In each of the N rounds, the one- and two-photon data are randomly sampled with replacement (resampled) to generate simulated data. The data thus simulated are given as inputs to algorithms 1–7, which return the simulated representative matrices

$$\begin{eqnarray}&&\{{{W}^{\prime }}^{b}\;:b\in \{1,...,N\},N\in {{\mathbb{Z}}}^{+}\}.\end{eqnarray} \tag{ 37 }$$

The pdf's of the simulated-matrix entries $\{W{{}^{\prime }}_{{ij}}^{b}\;:b\in \{1,...,N\}\}$ converge to the pdf's of the respective elements $\{{W}_{{ij}}\}$ for large enough N [54, 55].

The simulated data are obtained in each round by resampling from the one- and two-photon experimental data as follows. Single-photon detection counts are simulated by resampling from the set $\{{N}_{{{ijb}}_{j}}\;:\;{b}_{j}\in \{1,...,B\}\}$ of experimental detection counts (line 17 of algorithm 8). Two-photon coincidence counts are simulated by shuffling residuals obtained on curve-fitting experimental data. Specifically, the algorithm (line 12) resamples from the set

$$\begin{eqnarray}&&\{r(\tau )={C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{exp}}(\tau )-{C}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\tau )\ :\tau \in T\}\end{eqnarray} \tag{ 38 }$$

of residuals obtained by fitting experimentally measured coincidence counts to function ${C}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\tau )$ (12). The resampled residuals are added to the fitted curve to generate the simulated data (line 14)⁹ . Algorithms 1–7 are used to obtain the simulated elements W_ij of the representative matrix. Finally, the error bars on the $\{{W}_{{ij}}\}$ are estimated by the standard deviation of the pdf of the elements.

This completes the characterization of representative matrix W and the error bars on its elements. The next section details a procedure for the scattershot characterization of the interferometer to reduce the experimental time required for characterizing a given interferometer.

Algorithm 8.  Bootstrap: Estimate error bars on $\{{W}_{{ij}}\}$.

Input:

• $k,{\rm{\Omega }},{f}_{1},{f}_{2}\;:{\rm{\Omega }}\to {{\mathbb{R}}}^{+}$   ▹ Spectral functions: same as algorithm 3.

• ${\ell },T=\{{\tau }_{1},{\tau }_{2},...,{\tau }_{{\ell }}\}\in {({\mathbb{R}}\cup 0)}^{{\ell }}$   ▹ Time delay values.

• $m,{C}^{\mathrm{cal}}(\tau ),{C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{exp}}(\tau )$ for $\tau \in T$ and $(i,{i}^{\prime },j,{j}^{\prime })$        ▹ Same as algorithms 3 and 6.

• $B,{N}_{{{ijb}}_{j}}\quad :\{1,...,m\}\times \{1,...,m\}\times \{1,...,B\}\to {{\mathbb{Z}}}^{+}$ as in algorithm 2.

• N     ▹ Number of bootstrapping samples.

Output:

• $\sigma (\mathrm{re}{W}_{{ij}}),\sigma (\mathrm{im}{W}_{{ij}})\;:\{1,...,m\}\times \{1,...,m\}\to {{\mathbb{R}}}^{+}$     ▹ Error in W elements.

1: procedure ${\rm{B}}{\rm{ootstrap}}$ $(k,{\rm{\Omega }},{f}_{1},{f}_{2},{\ell },T,\gamma ,{C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{exp}}(\tau ),{\theta }_{{ij}},B,{N}_{{ijb}},N)$

2:  $A\leftarrow \{\mathrm{cos}\vartheta ,\mathrm{sin}\vartheta ,\mathrm{sin}\vartheta ,\mathrm{cos}\vartheta \}$, ${\rm{\Phi }}\leftarrow \{0,\pi /2,\pi /2,0\}$

3:  ${\mathrm{Residuals}}^{\mathrm{cal}}(\tau )\leftarrow {C}^{\mathrm{cal}}(\tau )-{\rm{C}}{\rm{OINCIDENCE}}(k,{\rm{\Omega }},{f}_{1},{f}_{2},{\ell },T,A,{\rm{\Theta }},\gamma )$

4:  ${\mathrm{NormalResiduals}}^{\mathrm{cal}}(\tau )\leftarrow \frac{{\mathrm{Residuals}}^{\mathrm{cal}}(\tau )}{{C}^{\mathrm{fit}}(\tau )}$  ▹ Assumption: ${\mathrm{Residuals}}^{\mathrm{cal}}(\tau )$ pdf width $\propto \;{C}^{\mathrm{fit}}(\tau )$.

5:  for $(i,{i}^{\prime },j,{j}^{\prime })\in \{1,2\}\times \{1,...,m\}\times \{1,2\}\times \{1,...,m\},i\ne {i}^{\prime },j\ne {j}^{\prime }$ do

6: $A\leftarrow \{{\alpha }_{{i}^{\prime }{j}^{\prime }},{\alpha }_{{i}^{\prime }j},{\alpha }_{{{ij}}^{\prime }},{\alpha }_{{ij}}\}$, ${\rm{\Phi }}\leftarrow \{{\theta }_{{i}^{\prime }{j}^{\prime }},{\theta }_{{i}^{\prime }j},{\theta }_{{{ij}}^{\prime }},{\theta }_{{ij}}\}$

7: ${C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{fit}}\leftarrow$ ${\rm{C}}{\rm{OINCIDENCE}}$($k,{\rm{\Omega }},{f}_{1},{f}_{2},{\ell },T,A,{\rm{\Theta }},\gamma$)

8: ${\mathrm{Residuals}}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\tau )\leftarrow {C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{exp}}(\tau )-{C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{fit}}(\tau )$

9: ${\mathrm{NormalResiduals}}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\tau )\leftarrow \frac{{\mathrm{Residuals}}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\tau )}{{C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{fit}}(\tau )}$

10:  end for

11:  for n  =  1 to N do

12: ${\mathrm{ShuffledNormalResiduals}}^{\mathrm{cal}}(\tau )\leftarrow$ Resample $| T|$ residuals from ${\mathrm{NormalResiduals}}^{\mathrm{cal}}(\tau )$

13: ${\mathrm{ShuffledResidual}}^{\mathrm{cal}}(\tau )\leftarrow {C}^{\mathrm{fit}}(\tau )\times {\mathrm{ShuffledNormalResiduals}}^{\mathrm{cal}}(\tau )$

14: ${C}^{n}(\tau )={C}^{\mathrm{fit}}(\tau )+\mathrm{ShuffledResidual}(\tau )$

15: ${\gamma }^{n}\leftarrow$ ${\rm{C}}{\rm{ALIBRATION}}$ $(k,{\rm{\Omega }},{f}_{1},{f}_{2},{\ell },T,{C}^{{\rm{b}}},\vartheta )$

16: for $(i,{i}^{\prime },j,{j}^{\prime })\in \{1,2\}\times \{1,...,m\}\times \{1,2\}\times \{1,...,m\},i\ne {i}^{\prime },j\ne {j}^{\prime }$ do

17:   ${\alpha }_{{ij}}^{n}\leftarrow$ ${\rm{MEAN}}$ $\sqrt{{N}_{11{b}_{1}}{N}_{{{ijb}}_{j}}/{N}_{1{{jb}}_{j}}{N}_{i1{b}_{1}}}$ from B values each of ${b}_{1},{b}_{j}$ drawn with replacement from $\{1,...,B\}$

18:   ${\mathrm{ShuffledNormalResiduals}}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\tau )\leftarrow | T|$ entries in ${\mathrm{NormalResiduals}}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\tau )$

19:   ${\mathrm{ShuffledResidual}}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\tau )\leftarrow {C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{fit}}(\tau )\times {\mathrm{ShuffledNormalResiduals}}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\tau )$

20:   ${C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{n}(\tau )={C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{fit}}(\tau )+{\mathrm{ShuffledResidual}}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\tau )$

21: end for

22: $\{{\theta }_{{ij}}^{n}\}$ = ${\rm{A}}{\rm{RGUMENT}}{\rm{C}}{\rm{ALC}}$($k,{\rm{\Omega }},{f}_{1},{f}_{2},{\ell },T,{C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{n}(\tau ),{\alpha }_{{ij}},\gamma$)

23: $\{{W}_{{ij}}^{\prime b}\}=$ ${\rm{M}}{\rm{AX}}{\rm{L}}{\rm{IKEY}}{\rm{U}}{\rm{NITARY}}$ $(\{{\alpha }_{{ij}}^{n}\},\{{\theta }_{{ij}}^{n}\})$

24:  end for

25:  for (i,j) in $\{1,...,m\}\times \{1,...,m\}$ do

26: $\sigma (\mathrm{re}{W}_{{ij}})=\mathrm{std}.\mathrm{dev}.(\{\mathrm{re}{W}_{{ij}}\});\sigma (\mathrm{im}{W}_{{ij}})=\mathrm{std}.\mathrm{dev}.(\{\mathrm{im}{W}_{{ij}}\})$

27:  end for

28: end procedure

In this section, we present a scattershot-based characterization approach to effect a reduction in the characterization time [58, 59]. Our scattershot approach reduces the time required to characterize an m-mode interferometer from ${\rm{O}}({m}^{4})$ to ${\rm{O}}({m}^{2})$ with constant error in the interferometer-matrix entries.

The straightforward approach of characterization involves coupling and decoupling light sources successively for each one- and two-photon measurement. In contrast, the scattershot characterization relies on coupling heralded nondeterministic single-photon sources to each of the input ports of the interferometer and detectors to each of the output ports. Controllable time delays are introduced at two input ports, which are labelled as the first and second ports. All sources and detectors are switched on and the controllable time-delay values are changed first for the first port and then for the second port.

Single-photon data are collected by selecting the events in which exactly one of the heralding detectors and exactly one of the output detectors register a photon simultaneously. Two-photon coincidence events at the outputs are counted when two heralding detectors register photons. The controllable time delays introduced at the first and second input ports ensure that each of the $2{(m-1)}^{2}$ coincidence measurements is performed. Note that our characterization procedure (algorithms 1–8) yields accurate estimates of interferometer parameters even when photon sources with different spectral functions are used. In summary, the required characterization data are collected by selectively recording one- and two-photon events. The setup for the scattershot characterization of an interferometer is depicted in figure 7.

Zoom In Zoom Out Reset image size

Now we quantify the experimental time required in the characterization of a linear optical interferometer. Our characterization procedure requires Bm² single-photon counting measurements and $2{(m-1)}^{2}$ coincidence-counting measurements to characterize an m-mode interferometer. We estimate the time required for each of these measurements such that random errors in the $\{{\alpha }_{{ij}}\},\{{\theta }_{{ij}}\}$ estimates remain unchanged with increasing m. To ensure constant error in the $\{{\alpha }_{{ij}}\},\{{\theta }_{{ij}}\}$ estimates, we require that the number of one- and two-photon detection counts remain unchanged with increasing m. The probability of photon detection at the output decreases with increasing m because of the concomitant decrease in the transmission amplitudes $\{{\alpha }_{{ij}}\}$.

The amplitudes $\{{\alpha }_{{ij}}\}$ drop as ${\rm{O}}(1/\sqrt{m})$ because of the unitarity of U [60]. Hence, one- and two-photon transmission probabilities (10) and (12) decrease as $1/m$ and $1/{m}^{2}$, respectively. More photons need to be incident at the interferometer input ports to offset this decrease in transmission probabilities. Therefore, maintaining a constant standard deviation in the $\{{\alpha }_{{ij}}\}$ and $\{{\theta }_{{ij}}\}$ measurements requires ${\rm{O}}(m)$ and ${\rm{O}}({m}^{2})$ scaling respectively in the number of incident photons, which amounts to an overall ${\rm{O}}({m}^{4})$ scaling in the experimental time requirement. Scattershot characterization allows ${(m-1)}^{2}$ different sets of the one- and two-photon data to be collected in parallel thereby reducing the time required to characterize the interferometer by a factor of ${(m-1)}^{2}$. The overall time required for the characterization decreases from ${\rm{O}}({m}^{4})$ to ${\rm{O}}({m}^{2})$ if the scattershot approach is employed.

Our analysis of scattershot characterization assumes that the coupling losses are small and that weak single-photon sources are used, i.e., that the probability of multi-photon emissions from the heralded sources is small as compared to single-photon emission probabilities. These assumptions are expected to hold for on-chip implementations of linear optics that have integrated single-photon sources and detectors.

Light sources used at each input port in our scattershot-based characterization procedure differ spectrally in generally. Our characterization procedure is accurate despite this difference because we measure source-field spectra and using these data in the curve-fitting procedure.

We have developed the scattershot approach which has advantages and disadvantages but on balance is a superior experimental approach to consecutive measurement. The advantage is that the time requirement for characterization if reduced by a factor that scales as ${\rm{O}}({m}^{2})$ . The disadvantage is the overhead of requiring one source at each input port and one detector at each output port. The disadvantage is not daunting because these requirements are commensurate with other active investigations of QIP such as LOQC and scattershot BosonSampling. In fact, state of the art implementations [59] meet our increased requirements for scattershot characterization.

In this section, we summarize our characterization procedure for a less formally oriented audience. We describe the processing of the collected experimental data by the various algorithms presented in section 3. We compare our procedure with the existing procedure for the characterization of linear optics using one- and two-photons [33]. We provide numerical evidence that our characterization procedure promises enhanced accuracy and precision even in the presence of shot noise and mode mismatch.

The experimental data required by our procedure to characterize an m-mode interferometer includes the following one- and two-photon measurements. The number ${N}_{{{ijb}}_{j}}$ (13) of single-photon detection events is counted at the jth output port when single photons are incident at the ith input port. This single-photon counting is repeated B times for each of the input ports and output ports, where B is chosen such that the cumulants of the set $\{{N}_{{{ijb}}_{j}}\ :{b}_{j}\in \{1,...,B\}\}$ converge. The single-photon counts $\{{N}_{{{ijb}}_{j}}\}$ are received by algorithm 2, which returns the $\{{\alpha }_{{ij}}\}$ (8) estimates using equation (18).

The spectral function ${f}_{j}(\omega )$ (2) of the light incident at each input port j is measured. This function is used by algorithm 1 to calculate the expected two-photon coincidence curves using equation (12). Fitting experimental data to these coincidence curves yields an accurate estimate of the mode-matching parameter during calibration and the arguments $\{{\theta }_{{ij}}\}$ in the argument-estimation procedure. Thus, the spectral function ${f}_{j}(\omega )$ serves as an input to the algorithms for the estimation of the mode-matching parameter and of the arguments $\{{\theta }_{{ij}}\}$ (algorithms 3–6).

The mode-matching parameter γ is estimated by performing coincidence measurement on a beam splitter that is separate from the interferometer but is constructed using the same material. First, we use single-photon data to estimate the reflectivity $\mathrm{cos}\vartheta$ of the beam splitter according to equation (21). Imperfect mode-matching changes the shape of the coincidence curve, and we find γ by comparing the shapes of (i) the curve expected for reflectivity $\mathrm{cos}\vartheta$ and (ii) the curve obtained experimentally. The estimated beam splitter reflectivity, the measured spectra and the coincidence counts are received as inputs by algorithm 3, which returns an estimate of γ.

Algorithm 6 uses two-photon coincidence counts to estimate the arguments $\{{\theta }_{{ij}}\}$. Coincidence counts are measured for the input ports $j,{j}^{\prime }$ and output ports $i,{i}^{\prime }$ for the $2{(m-1)}^{2}$ sets

$$\begin{eqnarray}&&(i,{i}^{\prime },j,{j}^{\prime })\in \{1,2\}\times \{1,...,m\}\\ &&\quad \times \{1,2\}\times \{1,...,m\},\quad i\ne {i}^{\prime },j\ne {j}^{\prime }\end{eqnarray} \tag{ 39 }$$

of input and output ports. In other words, coincidence counts are measured for different choices of two input ports and two output ports, such that each of the choices includes (i) either the first or the second input ports and (ii) either the first or the second output port. Algorithm 6 receives as input the measured spectra, the $\{{\alpha }_{{ij}}\}$ values estimated by algorithm 2, the γ value estimated by algorithm 3 and the two-photon coincidence data for the choice (39) of input ports. The algorithm returns the $\{{\theta }_{{ij}}\}$ estimates. The computed estimates of $\{{\alpha }_{{ij}}\}$ and of $\{{\theta }_{{ij}}\}$ yield the representative unitary matrix W (36) that has maximum likelihood of describing the characterized interferometer (algorithm 7). This completes a summary of our procedure for characterization of the interferometer.

Algorithm 8 employs bootstrapping to find the error bars on the elements $\{{W}_{{ij}}\}$ of the characterized unitary matrix. The bootstrapping procedure uses the experimental data that is received by algorithms 1–7 and repeatedly simulates experiments by resampling from the experimental data. The number N of repetitions is chosen such that the pdf's of the $\{{W}_{{ij}}\}$ elements over many rounds of simulation converge. The error bars on the $\{{W}_{{ij}}\}$ elements are computed based on the estimated pdf's of the elements. Our procedure thus enables the estimation of meaningful error bars on the characterized unitary matrix.

Bootstrapping is employed to test the goodness of fit between the experimental curve and expected curves [61]. Experiments [11, 36] can employ bootstrapping instead of the incorrect ${\chi }^{2}$-confidence measure to test if the data are consistent with quantum predictions or with the classical theory.

Finally, we recommend a scattershot approach for reducing the experimental time required to characterize interferometers. The approach involves coupling heralded non-deterministic single-photon sources at each of the input ports and single-photon detectors at each of the output ports. All the sources and the detectors are switched on in parallel. Single-photon counts are recorded selectively as two-photon coincidences between the heralding detectors and the output detectors, while two-photon events are recorded when two heralding detectors and two output detectors record photons. Controllable time delays are introduced at the first and second input ports so coincidences between each of the $2{(m-1)}^{2}$ choices (39) of input and output ports are recorded. The scattershot approach reduces the experimental time required to characterize an m-mode interferometer from ${\rm{O}}({m}^{4})$ to ${\rm{O}}({m}^{2})$.

Now we compare and contrast our procedure with the Laing–O'Brien procedure [33]. Our procedure is inspired by the Laing–O'Brien procedure in that it employs (i) a 'real-bordered' parameterization (8) of the representative matrix and modelling of linear losses at the interferometer ports, (ii) a ratio of single-photon data to estimate the complex amplitudes of the matrix elements and (iii) an iterative procedure that uses two-photon data to estimate the amplitudes of the complex arguments and to estimate the signs of the complex arguments.

Our procedure differs from the Laing–O'Brien [33] procedure in that we use averaged value (18) of the ratio of single-photon detections over many runs rather than the ratio of averaged values. This difference ensures accuracy of our procedure even under fluctuation in the number of incoming photons. Such fluctuations might arise from fluctuations in pump strength of the single photon source or in the strength of coupling between photon source and interferometer.

Another advance in our method is the curve-fitting procedure for estimating complex arguments of interferometer matrices. The Laing–O'Brien procedure requires coincidence-curve visibilities to estimate complex arguments ${\alpha }_{{ij}}$ . Whereas the Laing–O'Brien procedure recommends coincidence probabilities be measured at zero time delay and also at time delays large as compared to the temporal spread of the wave-packet, in practice, current implementations determine the visibilities by fitting experimental data to Gaussian curves [35, 63–67]. These implementations are flawed because source spectra differ from Gaussian in general. Our procedure is accurate because the data are fit to curves computed from spectral functions, rather than fitting to Gaussians. Figure 8 illustrates the distinction between fitting experimental coincidence counts to the coincidence function (12) simulated using spectra and fitting to Gaussian functions. Figure 9 demonstrates the increase in accuracy and precision of characterization by using the correct curve-fitting function.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We introduce the calibration subroutine, which relies on the estimation of the mode mismatch in the source field. Spatial and polarization mode mismatch is not an issue of major concern in waveguide-based interferometers, which typically operate in the single-photon regime. In these interferometers, the calibration step of our procedure can be neglected without decreasing accuracy. The mode-mismatch parameter γ, which is an input of the curve-fitting procedure, is set to unity.

Zoom In Zoom Out Reset image size

In the context of bulk-optics, our calibration step ensures accuracy and precision if (i) γ is identified as the maximum-possible source overlap in the spatial and polarization degrees of freedom and (ii) the experimentalist adjusts the setup to maximize coincidence visibility for the calibrating beam splitter and for each choice of interferometer inputs ports. Such an adjustment will ensure that the source overlap acquires its maximum-possible value γ in each of the coincidence-curve measurements. This maximum value is a property of the sources used and is independent of source alignment and focus so is expected to remain unchanged between different confidence measurements. Figure 10 demonstrates the increase in accuracy and precision of characterization by using the calibration procedure.

Other advances made in our characterization procedure over existing procedures include (i) a maximum likelihood estimation approach to determine the unitary matrix that best fits the data (ii) a bootstrapping based procedure to obtain meaningful estimates of precision and (iii) a scattershot-based procedure to improve the experimental requirements of characterization.

In conclusion, we devise a one- and two-photon interference procedure to characterize any linear optical interferometer accurately and precisely. Our procedure provides an algorithmic method for recording experimental data and computing the representative transformation matrix with known error.

The procedure accounts for systematic errors due to spatiotemporal mode mismatch in the source field by means of a calibration step and corrects these errors using an estimate of the mode-matching parameter. We measure the spectral function of the incoming light to achieve good fitting between the expected and measured coincidence counts, thereby achieving high precision in characterized matrix elements. We introduce a scattershot approach to effect a reduction in the experimental requirement for the characterization of interferometer. The error bars on the characterized parameters are estimated using bootstrapping statistics.

Bootstrapping computes accurate error bars even when the form of experimental error is unknown and is, thus, advantageous over the Monte Carlo method. Hence, our bootstrapping-based procedure for estimating error bars can replace the Monte Carlo method used in existing linear-optics characterization procedures. We open the possibility of applying bootstrapping statistics for the accurate estimation of error bars in photonic state and process tomography.

We thank Matthew A Broome, Jens Eisert, Dirk R Englund, Sandeep K Goyal, Hubert de Guise and Michael J Steel for useful comments. ID and BCS acknowledge AITF, China Thousand Talent Plan and NSERC and USARO for financial support. AK acknowledges CryptoWorks21 and NSERC for financial support.

In this section, we describe an instability in our characterization procedure, which can yield large error in the $\{{W}_{{ij}}\}$ output for small error in the experimental data ${C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{exp}}(\tau )$ in case of certain interferometers W. We present a strategy to circumvent this instability by means of collecting and processing additional experimental data.

The instability in the characterization procedure arises because of an instability in estimation of $\{\mathrm{sgn}{\theta }_{{ij}}\}$ (algorithm 5). Small error in the measured coincidence counts can lead to the wrong inference of $\mathrm{sgn}{\theta }_{{ij}}$,which can lead to a large error $\parallel W-U\parallel$ in the characterized matrix W. Recall that algorithm 5 uses the identity $\mathrm{sgn}{\theta }_{{ij}}=\mathrm{sgn}(| {\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}-{\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}^{+}| -| {\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}-{\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}^{-}| )$ (27) to determine the sign of the arguments,where ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\pm }\overset{{\rm{def}}}{=}| {\theta }_{{i}^{\prime }{j}^{\prime }}-{\theta }_{{{ij}}^{\prime }}-{\theta }_{{i}^{\prime }j}\pm | {\theta }_{{ij}}| |$ and the values of ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }},{\theta }_{{i}^{\prime }{j}^{\prime }},{\theta }_{{i}^{\prime }j},{\theta }_{{{ij}}^{\prime }},| {\theta }_{{ij}}|$ are estimated by curve fitting.

Random and systematic error in measured coincidence counts can lead to estimate of variables ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }},{\theta }_{{i}^{\prime }{j}^{\prime }},{\theta }_{{i}^{\prime }j},{\theta }_{{{ij}}^{\prime }},| {\theta }_{{ij}}|$ differing from their actual values. The estimation of $\mathrm{sgn}{\theta }_{{ij}}$ is unstable if the ${\theta }_{{i}^{\prime }{j}^{\prime }}-{\theta }_{{i}^{\prime }j}-{\theta }_{{{ij}}^{\prime }}$ term (27) is close to 0 or π because, in this case, a small error in the ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}$ estimate can lead to an incorrect $\mathrm{sgn}{\theta }_{{ij}}$ estimate. In other words, the sign estimates are unstable if the values of

$$\begin{eqnarray}&&{\theta }_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{ref}}=\mathrm{min}[{\theta }_{{i}^{\prime }{j}^{\prime }}-{\theta }_{{{ij}}^{\prime }}-{\theta }_{{i}^{\prime }j},\pi -({\theta }_{{i}^{\prime }{j}^{\prime }}-{\theta }_{{{ij}}^{\prime }}-{\theta }_{{i}^{\prime }j})],\end{eqnarray} \tag{ A.1 }$$

are small compared to the error in our ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }},{\theta }_{{i}^{\prime }{j}^{\prime }},{\theta }_{{i}^{\prime }j},{\theta }_{{{ij}}^{\prime }},| {\theta }_{{ij}}|$ estimates.

We mitigate the sign-inference instability by making two modifications to our characterization procedure; the first modification removes instability from the sign-inference of the second row and second column elements whereas the second modification prevents incorrect inference of the remaining signs. The inference of $\{\mathrm{sgn}{\theta }_{i2}\},\{\mathrm{sgn}{\theta }_{2j}\}$ (lines 14–17, figures 5(b) and (c)) is unstable if

$$\begin{eqnarray}&&{\theta }_{i2j2}^{\mathrm{ref}}=\mathrm{min}({\theta }_{22},\pi -{\theta }_{22})\end{eqnarray} \tag{ A.2 }$$

is small as compared to the error in the ${\beta }_{i2j2},{\theta }_{22},{\theta }_{2j},{\theta }_{i2},| {\theta }_{{ij}}|$ estimates. Hence, we relabel the interferometer ports such that ${\theta }_{22}$ is as far away from 0 and π as possible. Specifically, after the amplitudes of the phases have been estimated (line 8 of algorithm 6), we choose $i,j$ for which $| {\theta }_{{ij}}-\pi /2|$ is minimum, and we swap the labels of input ports $2,j$ and output ports $2,i$. We measure two-photon coincidence counts based on this new labeling and process it using algorithm 6. The instability in the procedure for estimation of the $\{{\theta }_{i2}\},\{{\theta }_{2j}\}$ signs is removed as a result of the relabeling.

The second modification is aimed at removing the instability in the remaining signs. The procedure estimates the remaining signs by using $\{{C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{exp}}(\tau )\}$ values for ${i}^{\prime }={j}^{\prime }=2$. The estimation of ${\theta }_{{ij}}$ is unstable if ${\theta }_{i2j2}^{\mathrm{ref}}$ is small as compared to the error in the ${\beta }_{i2j2},{\theta }_{22},{\theta }_{2j},{\theta }_{i2},| {\theta }_{{ij}}|$ estimates. We make a heuristic choice of a threshold angle ${\theta }^{{\rm{T}}}$ that accounts for the error in these variables, and we reject any $\mathrm{sgn}{\theta }_{{ij}}$ inferred using ${\theta }_{i2j2}^{\mathrm{ref}}\leqslant {\theta }^{{\rm{T}}}$. Additional two-photon coincidence counting is performed and employed to estimate these values of ${\theta }_{{ij}}$, as detailed in the following lines that can be added to the algorithm to remove the instability

17 + 1: for (i, j) in $\{3,...,m\}\times \{3,...,m\}$ do

2: if ${\theta }_{i2j2}^{r}\lt {\theta }_{T}$ then

3: Choose ${i}^{\prime }\ne 1,i$ and ${j}^{\prime }\ne 1,j$ such that $| {\theta }_{{i}^{\prime }{j}^{\prime }}-{\theta }_{{{ij}}^{\prime }}-{\theta }_{{j}^{\prime }j}|$ is closest to $\pi /2$.

4:  ${C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{exp}}(\tau )\;\leftarrow \;$Coincidence counts for input ports $j,{j}^{\prime }$and output ports $i,{i}^{\prime }$.

5: $A\leftarrow \{0,0,0\},{\rm{\Phi }}\leftarrow \{0,0,0\}$

6: ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}\leftarrow {\rm{A}}{\rm{RGUMENT}}2{\rm{P}}{\rm{ORT}}({C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{exp}}(\tau ),{\rm{\Omega }},{f}_{1},{f}_{2},T,A,{\rm{\Phi }},\gamma )$

7: ${\theta }_{{ij}}\leftarrow {\theta }_{{ij}}$ ${\rm{S}}{\rm{IGN}}{\rm{C}}{\rm{ALC}}$openup-->$({\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }},{\theta }_{{i}^{\prime }{j}^{\prime }},{\theta }_{{{ij}}^{\prime }},{\theta }_{{i}^{\prime }j},| {\theta }_{{ij}}| )$

8: end If

9: end for

Figure A1 illustrates the rejection of those $i,j$ choices for which ${\theta }_{i2j2}^{\mathrm{ref}}\leqslant {\theta }^{{\rm{T}}}$ and the use of ${C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{exp}}(\tau ),{j}^{\prime }\ne 2$ counts to obtain a correct estimate of $\mathrm{sgn}{\theta }_{{ij}}.$ We thus remove the instability in the estimation of $\{{\theta }_{{ij}}\}$ and in the estimation of the representative matrix W.

Zoom In Zoom Out Reset image size

Our characterization procedure employs curve fitting in algorithm 3 to estimate the mode-matching parameter γ and in algorithms 4–6 to estimate $\{{\theta }_{{ij}}\}$ values. The curve-fitting procedure determines those values of unknown parameters that maximize the fitting between experimental and expected coincidence data. In this section, we describe the inputs and outputs of the curve-fitting subroutines. We present heuristics to compute good initial guesses of the fitted parameters.

The curve-fitting subroutine receives as input (i) the choice of parameters to be fitted; (ii) the coincidence counts $\{{C}_{{{ii}}^{\prime }{{jj}}^{\prime }}^{\mathrm{exp}}(\tau )\};$ (iii) an objective function, which characterizes the least-square error between expected and experimental counts; and (iv) the initial guesses for each of the fitted parameters. The output of the curve-fitting subroutine is the set of parameter values that optimize the objective function.

The first input to the subroutine is the choice of the parameters to be fit. The curve-fitting subroutine fits three parameters. One of these three (namely the mode-matching parameter γ in algorithm 3 or the $| {\theta }_{{ij}}|$ or ${\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}$ value in algorithm 6) is related to the shape of the curve, whereas the other two are related to the ordinate scaling and the abscissa shift of the curve respectively. The ordinate scaling factor comprises the unknown losses $\{{\kappa }_{i},{\nu }_{j}\}$, transmission factors $\{{\lambda }_{i},{\mu }_{j}\}$ and the incident photon-pair count. The horizontal shift factor accounts for the unknown zero of the time delay between the incident photons. The algorithm returns the values of the shape parameter, the abscissa shift and the ordinate scaling that best fit the given coincidence curves.

The objective function quantifies the goodness of fit between the experimental data and the parameterized curve. We use a weighted sum

$$\begin{eqnarray}&&\displaystyle \sum _{\tau \in T}w(\tau )| {C}^{\mathrm{exp}}(\tau )-{C}^{\prime }(\tau ){| }^{2}\end{eqnarray} \tag{ B.1 }$$

of squares between the experimental data and the fitted curve as the objective function [44] for weighs $w(\tau )$. We assume that the pdfs of the residues are proportional to $\sqrt{{C}^{\mathrm{exp}}(\tau )}$ and we assign the weights

$$\begin{eqnarray}w(\tau )=\left\{\begin{array}{cc}1/{C}^{\mathrm{exp}}(\tau ) & \mathrm{if}\ {C}^{\mathrm{exp}}(\tau )\ne 0\\ \quad 1 & \mathrm{if}\ {C}^{\mathrm{exp}}(\tau )=0\end{array}\right.\end{eqnarray} \tag{ B.2 }$$

to the squared sum of residues. In case the pdf's of the residuals for different values of τ is not known, standard methods for non-parametric estimation of residual distribution can be employed to estimate the pdf's [56, 57]. Thus, the curve fitting algorithm returns those values of the fitting parameters that that minimize weighted sum of squared residues between experimental and fitted data.

The curve-fitting procedure optimizes the fitness function over the domain of the fitting parameter values. Like other optimization procedures, the convergence of curve fitting is sensitive to the initial guesses of the fitting parameters. The following heuristics give good guesses for the three fitting parameters. We guess the ordinate scaling as the ratio

$$\begin{eqnarray}&&\displaystyle \frac{{C}^{\mathrm{exp}}(\infty )}{{C}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\infty )}\end{eqnarray} \tag{ B.3 }$$

of the experimental coincidence counts

$$\begin{eqnarray}&&{C}^{\mathrm{exp}}(\infty )\overset{{\rm{def}}}{=}\displaystyle \frac{{C}^{\mathrm{exp}}({\tau }_{1})+{C}^{\mathrm{exp}}({\tau }_{{\ell }})}{2},\end{eqnarray} \tag{ B.4 }$$

to the coincidence probability ${C}_{{{ii}}^{\prime }{{jj}}^{\prime }}(\infty )$ for large (compared to the temporal length of the photon) time-delay values. The γ value is guessed for algorithm 3 as the ratio of the visibility of the experimental curve to the expected visibility in the curve. The initial guesses for $\vartheta \equiv | {\theta }_{{ij}}|$ and $\vartheta \equiv {\beta }_{{{ii}}^{\prime }{{jj}}^{\prime }}$ are based on the known estimate of γ and the visibility

$$\begin{eqnarray}&&V=\displaystyle \frac{2\gamma {\mathrm{cos}}^{2}\vartheta {\mathrm{sin}}^{2}\vartheta }{{\mathrm{cos}}^{4}\vartheta +{\mathrm{sin}}^{4}\vartheta },\end{eqnarray} \tag{ B.5 }$$

of the curve. As there are four kinds of curves (see figure B1 ) possible for different values of the shape parameter ($\gamma ,| {\theta }_{{ij}}| ,{\beta }_{{{ii}}^{\prime }{j}^{\prime }}$ ), another approach is to perform curve fitting four times, each time with a value from the set $\pi /4,3\pi /4,5\pi /4,7\pi /4$ of initial guesses and choose the fitted parameters that optimize the objective function. Finally, the initial value of the abscissa shift parameter is guessed such that the global maxima or minima (whichever is further from the mean of the coincidence-count values over τ) of the coincidence curve is at zero time delay.

Zoom In Zoom Out Reset image size

In summary, the curve fitting procedure uses the measured coincidence counts, the objective function and the initial guesses to compute the best fit parameters. This completes our description of the curve-fitting procedure and of heuristics that can be employed to computed the initial guesses for the fitted parameters.