Scalable on-chip quantum state tomography

TL;DR

This paper introduces a scalable quantum state tomography method that significantly reduces measurement complexity, demonstrated on a photonic chip with high fidelity for two- and three-photon states.

Contribution

The paper presents a novel measurement technique that scales linearly with qubit number, enabling efficient characterization of large quantum systems.

Findings

Measurement complexity scales linearly with qubits

Achieved 99.67% fidelity in state reconstruction

Demonstrated on integrated photonic chip for 2- and 3-photon states

Abstract

Quantum information systems are on a path to vastly exceed the complexity of any classical device. The number of entangled qubits in quantum devices is rapidly increasing and the information required to fully describe these systems scales exponentially with qubit number. This scaling is the key benefit of quantum systems, however it also presents a severe challenge. To characterize such systems typically requires an exponentially long sequence of different measurements, becoming highly resource demanding for large numbers of qubits. Here we propose a novel and scalable method to characterize quantum systems, where the complexity of the measurement process only scales linearly with the number of qubits. We experimentally demonstrate an integrated photonic chip capable of measuring two- and three-photon quantum states with reconstruction fidelity of 99.67%.

Equations

$$U_{1,j}{E}_1^{(in)}+U_{2,j}{E}_2^{(in)}=E_j^{(out)}$$

$${\left||U_{1,j}||E_1^{(in)}|+|U_{2,j}||E_2^{(in)}|e^{i({\varphi }_j+\Delta d/\lambda )}\right|}^2=|E_j^{(out)}{|}^2$$

$$|E_j^{(out)}{|}^2=\frac{a}{2}[1+\cos ({\varphi }_j+\Delta d/\lambda )]+b$$

$$|U_{1,j}|E_1^{(in)}+|U_{2,j}|e^{i({\varphi }_j+\theta )}{E}_2^{(in)}=E_j^{(out)}$$

$$|\Psi ⟩=\sum _{n,m}\int d{\omega }_{1}d{\omega }_2{\Psi }_{n,m}({\omega }_1,{\omega }_2){\hat{a}}_n^{†}({\omega }_1){\hat{a}}_m^{†}({\omega }_2)|0⟩$$

$${\Psi }_{n,m}({\omega }_1,{\omega }_2)={\varphi }^{(a)}({\omega }_1){\varphi }^{(b)}({\omega }_2){\delta }_{n,1}{\delta }_{m,2}$$

$$\begin{gathered} |{\mathbf{D}}_{n,m}⟩⟨{\mathbf{D}}_{n^{\prime },m^{\prime }}|= \\ \int d{\omega }_{3}d{\omega }_4{D}^{\ast }({\omega }_3)F^{\ast }({\omega }_4){\hat{a}}_n^{†}({\omega }_3){\hat{a}}_m^{†}({\omega }_4)|0⟩⟨0|D({\omega }_3)F({\omega }_4){\hat{a}}_{n^{\prime }}({\omega }_3){\hat{a}}_{m^{\prime }}({\omega }_4), \end{gathered}$$

$$\begin{gathered} {\hat{\rho }}_{n,m,n^{\prime },m^{\prime }}=⟨\Psi ||{\mathbf{D}}_{n,m}⟩⟨{\mathbf{D}}_{n^{\prime },m^{\prime }}||\Psi ⟩= \\ \frac{1}{2}⟨0|\int d{\omega }_{1}d{\omega }_2{\varphi }^{(a)\ast }({\omega }_1){\varphi }^{(b)\ast }({\omega }_2){\hat{a}}_1({\omega }_1){\hat{a}}_2({\omega }_2)\times \\ \int d{\omega }_{3}d{\omega }_4{D}^{\ast }({\omega }_3)F^{\ast }({\omega }_4){\hat{a}}_n^{†}({\omega }_3){\hat{a}}_m^{†}({\omega }_4)|0⟩⟨0|D({\omega }_3)F({\omega }_4){\hat{a}}_{n^{\prime }}({\omega }_3){\hat{a}}_{m^{\prime }}({\omega }_4)\times \\ \int d{\omega }_{5}d{\omega }_6{\varphi }^{(a)}({\omega }_5){\varphi }^{(b)}({\omega }_6){\hat{a}}_1^{†}({\omega }_5){\hat{a}}_2^{†}({\omega }_6)|0⟩. \end{gathered}$$

$${\hat{\rho }}_{1,2,1,2}=\int d{\omega }_{3}d{\omega }_4{\varphi }^{(a)\ast }({\omega }_3){\varphi }^{(b)\ast }({\omega }_4){\varphi }^{(a)}({\omega }_3){\varphi }^{(b)}({\omega }_4)=1/2,$$

$${\hat{\rho }}_{2,1,2,1}=\int d{\omega }_{3}d{\omega }_4{\varphi }^{(a)\ast }({\omega }_4){\varphi }^{(b)\ast }({\omega }_3){\varphi }^{(a)}({\omega }_4){\varphi }^{(b)}({\omega }_3)=1/2,$$

$$\begin{gathered} {\hat{\rho }}_{1,2,2,1}={\hat{\rho }}_{2,1,1,2}^{\ast }=\int d{\omega }_{3}d{\omega }_4{\varphi }^{(a)\ast }({\omega }_3){\varphi }^{(b)\ast }({\omega }_4){\varphi }^{(a)}({\omega }_4){\varphi }^{(b)}({\omega }_3)/2= \\ \int d{\omega }_3{\varphi }^{(a)\ast }({\omega }_3){\varphi }^{(b)}({\omega }_3)\int d{\omega }_4{\varphi }^{(b)\ast }({\omega }_4){\varphi }^{(a)}({\omega }_4)/2. \end{gathered}$$

$$\hat{\rho }=\frac{1}{2}\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 1& I_{a,b}{I}_{b,a}& 0\\ 0& I_{a,b}{I}_{b,a}& 1& 0\\ 0& 0& 0& 0\end{array}\right],$$

Full text

The standard way to characterize a quantum system is known as quantum state tomography Vogel:1989-2847:PRA ; Altepeter:2005-105:RAR . It involves measuring expectation values of a complete set of observables and using these to reconstruct the system’s density matrix James:2001-52312:PRA ; Haffner:2005-643:NAT ; Leibfried:2005-639:NAT ; Lvovsky:2009-299:RMP ; Bayraktar:2016-20105:PRA ; Shadbolt:2012-45:NPHOT . To characterize an $N$-qubit state, $2^{2N}$ different observables are measured James:2001-52312:PRA , thus the measurement apparatus must be reconfigured exponentially many times, which is impractical for large states. Furthermore many of the expectation values measured will be vanishingly small, and thus contribute little useful information. Finally, even if all the measurements can be completed, the task of reconstructing the density matrix from measurement data becomes computationally challenging for high qubit-number states Banaszek:2013-125020:NJP .

New approaches to quantum state tomography are being developed in an effort to increase its practicality and efficiency. Some approaches seek to avoid unnecessary measurements by assuming that the system is in particular low rank states, such as sparse states Gross:2010-150401:PRL ; Oren:2016-226:OPT or low dimensional matrix product states Cramer:2010-149:NCOM . Alternatively tomography can be ‘self-guided’, where real-time feedback of measurement results guides the next choice of the measurement basis Ferrie:2014-190404:PRL ; Chapman:2016-40402:PRL , helping to avoid taking measurements that have limited utility for analyzing the state. It has been shown that tomography procedures involving some global quantum measurements have increased error robustness relative to using only local qubit measurements, and thus can be completed in less time Miranowicz:2014-62123:PRA . The computational burden of inverting large data sets to find the density matrix is reduced with simple real-time optimization algorithms in self guided tomography, or can be completely avoided with systems for direct projection of density matrix parameters Lundeen:2011-188:NAT ; Bolduc:2016-10439:NCOM ; Thekkadath:2016-120401:PRL . However all these approaches rely on a common measurement paradigm, whereby different characteristics of a system are measured sequentially, thus they become exponentially complex to implement as the number of parameters in the density matrix scales exponentially with qubit number.

Here we present a quantum tomography method with complexity that scales linearly with qubit number. This is achieved by leveraging quantum systems’ greatest strength, the simultaneous occupation of exponentially many states, in the measurement process. Instead of preforming a sequence of different measurements on the state, we design a single static measurement system [Fig. 7(a)] that preforms one, many-outcome measurement. Since the state is spread coherently across all the outputs, the number of different measurement outcomes follows a similar scaling to the number of parameters in the density matrix. Thus the exponential scaling of quantum sates can be balanced by a similar scaling in the amount of information extracted from measurement of the state. This leads to the striking benefit that the required physical complexity of multi-outcome measurement only scales linearly with the number of qubits in the state being measured, in contrast to the usual exponential scaling.

In the context of photonic quantum states our approach also removes the need to build complex reconfigurable measurement systems James:2001-52312:PRA , instead allowing full quantum tomography with just static linear optical circuits, which can easily be implemented on photonic chips. This avoids the problem faced by conventional approaches, where some measurements provide little useful information about the underlying state. This is because in our approach, the full complement of measurements are performed simultaneously, and thus the most important correlation detections for reconstructing a particular state naturally have the highest count rates. Our approach is based on interfering all the photons through a special unitary transformation, thus can be optimized to incorporate nonlocal measurement, allowing the error robustness to be increased. Furthermore, we show that our approach is compatible with computationally scalable reconstruction, avoiding resource intensive direct inversion.

We explain our method for the case of $N$-photons in an arbitrary pure or mixed state featuring spatial quantum entanglement between two input ports, although a larger number of ports can be considered as well Titchener:2016-4079:OL . The $N$-photon quantum state is described by a density matrix, $\hat{\rho}_{in}$, which has $2^{2N}$ real parameters. We apply a single (fixed) linear transformation, $\hat{U}$, to map the quantum state from two inputs to a larger number of output waveguides, $M$, as illustrated in Fig. 7(a). Then, we obtain information about the quantum state by measuring coincidences, $\Gamma$, in the arrival time of photons to different combinations of the single-photon detectors. There are $M!/(M-N)!$ different $N$-photon correlations which can be measured, assuming the detectors cannot resolve photon numbers, but can distinguish photons (i.e. by arrival time for time-bin encoded states). The reconstruction of the density matrix is possible if the number of different correlations is larger than the number of unknown density matrix elements. We establish that the minimum number of output ports required scales linearly with the number of photons and approaches $M=N+3$ for large photon numbers, both in cases of distinguishable (red shaded area) and indistinguishable (blue line) photons as shown in Fig. 7(b); mathematical details are provided in the Sec. 2 of the supplementary information.

We experimentally demonstrate our approach by performing tomography of spatially entangled mixed or pure states of two indistinguishable photons, using a specially designed on-chip laser-written waveguide circuit Meany:2015-363:LPR . The action of the circuit on the input single photon state $\ket{\psi}=\ket{0}+i\ket{1}$ is shown in Fig. 8(a). The circuit allows full reconstruction of the input density matrix just by measuring the output two-photon coincidences with non-photon-number resolving single-photon detectors.

The circuit was optimized to make the tomographic reconstruction highly robust to measurement errors. Each output waveguide carries information from a different vector on the input Bloch sphere. The Bloch vectors were determined experimentally using a classical characterization method Heilmann:2015-96:RAR , and are shown graphically in Fig. 8(b). The equal spacing of vectors around the Bloch sphere gives the device maximum robustness to errors in the tomography procedure Foreman:2015-263901:PRL , by essentially realizing non-local measurements. This can be confirmed by calculating the condition number of the transfer function of the chip, lower values of which correspond to higher robustness of the state reconstruction to measurement errors. The experimentally realized device has a condition number of $\simeq 5$, which is better than the condition number of $\simeq 9$ for typical to tomography Miranowicz:2014-62123:PRA .

To test the performance of the device we prepared a range of different two-photon entangled quantum states, coupled them into the chip, and measured the output correlations as schematically shown in Fig. 8(c). We first analyze an anti-bunched state, which in an ideal form is described by the pure wavefunction $\ket{\psi}=\left(\ket{01}+\ket{10}\right)/\sqrt{2}$ and the corresponding density matrix $\hat{\rho}=\ket{\psi}\bra{\psi}$. We present in Fig. 9(a) the experimentally measured probabilities of detecting the photon-pair in a given pair of output waveguides, and the reconstructed real and imaginary parts of the density matrix, as indicated by labels. We confirm that this is indeed an anti-bunched state, with the fidelity of $95.0\%$. Furthermore, this measurement permits us to get information about the spectral overlap of the pair of photons, since the observed correlations exhibit a generalized form of Hong-Ou-Mandel interference Hong:1987-2044:PRL ; Spring:2017-90:OPT ; Mahrlein:2015-15833:OE , see Sec. 1 of the supplementary material for details.

We also prepared N00N states, with wavefunctions given by $\ket{\psi}=\left(\ket{00}+e^{i\phi}\ket{11}\right)/\sqrt{2}$, where $\phi$ is a phase shift. The phase shift is determined by the photon propagation before the chip, and because the shift is double the value that would accumulate classically Ou:1990-2957:PRA it is highly sensitive to the environment. Experimentally, we explored this quantum-enhanced sensitivity by propagating a two-photon N00N state through one meter long optical fibers before the chip. The accumulated phase was very sensitive to fiber stress, varying by up to $2\pi$ on the scale of a few minutes. With our approach we observed experimentally the temporal variation of the phase in the density matrix using an integration time of $20$ seconds with a photon pair detection rate of $30$Hz. We show typical two-photon correlations at different times in Figs. 8(b) and (c) and the corresponding reconstructed N00N states, with phases determined to be $\phi=0$ in Fig. 9(b) and $\phi=\pi/2$ in Fig. 9(c). The fidelity of both states exceeds $94\%$. Thus, we can observe with high precision the density matrix of a quantum state that is varying over time (a video of the time evolution of the density matrix is included in the supplementary material).

The reconstruction of the density matrices from the measured correlations was carried out using a computationally scalable algorithm. This is important, since although the number of detectors and waveguides in the circuit scales linearly with the number of qubits to be measured, the number of photon detections required for reconstruction still necessarily scales exponentially, and processing of measurement data can be extremely resource demanding. We employ an optimization technique known as the simultaneous perturbation stochastic approximation Spall:1992-332:RAR , similar to the algorithm formulated for self-guided quantum state tomography Ferrie:2014-190404:PRL . It minimizes a distance measure between the true state and the algorithm’s current guess. We use the least squares distance between the measured correlations and the correlations that would be produced by the current guess of the density matrix, $|\Gamma_{\text{meas.}}-\Gamma_{\text{guess}}|^{2}$. The reconstruction fidelity of the algorithm is $99.67\%$ after 1500 iterations (see Methods).

We demonstrate recovery of the mixed state $\hat{\rho}_{\text{mix}}=\left(\hat{\rho}_{\text{anti-bunched}}+\hat{\rho}_{\text{N00N}}\right)/2$ using the algorithm. The inset of Fig. 6(a) shows the experimentally measured correlations corresponding to $\hat{\rho}_{\text{mix}}$, obtained by numerically combining separately recorded raw coincidence data from an anti-bunched state and a N00N state. The main plot in Fig. 6(a) shows 1000 realizations of the algorithm as red lines, and density matrices at iteration numbers 1, 30 and 300 are shown in Fig. 6(b) for the realization highlighted in blue in Fig. 6(a). Importantly, our method offers the same computational advantage as self-guided tomography, but without a need for complex reconfigurable measurements which so far restricted this approach to only pure states Chapman:2016-40402:PRL .

The chip presented in Fig. 8(a) is also capable of tomography of degenerate three-photon states. We demonstrate this using the experimentally determined transfer function, $\hat{U}$, and simulate the propagation of a three-photon Greenberger-Horne-Zeilinger (GHZ) state [Fig. 5(a)] through the chip. The simulated output three-photon correlations are shown in Fig. 5(b), where Gaussian noise with standard deviation $5\%$ of the peak correlation value has been added to each element. Reconstruction of the input density matrix gives highly accurate results despite this noise. The real part of the recovered density matrix is shown in Fig. 5(d), which closely matches the three-photon GHZ state in Fig. 5(c). This provides a significantly simpler and more stable platform for three-photon tomography compared to previous realizations Hamel:2014-801:NPHOT .

In our work, we have demonstrated that a fixed linear optical transformation can be devised to allow complete quantum tomography of $N$-photon states. Importantly, the complexity of the transformation only scales linearly with the number of photons in the state, in contrast to the exponential scaling in the number of measurements required in usual tomography White:1999-3103:PRL . Due to its simplicity and lack of tunable elements, this approach to quantum measurement is uniquely suited to integration with on-chip single photon detectors Najafi:2015-5873:NCOM for a fully on-chip tomography scheme. This provides a promising way to facilitate the characterization and development of increasingly complex quantum communication and computation systems.

I Methods

Direct waveguide writing. We write our waveguides into transparent fused silica wafers (Corning 7980 ArF Grade), using ultrashort laser pulses ($\tau<150$ fs, $\lambda=800$ nm) that are focused 250 $\mu$m below the sample surface using a 20$\times$ microscope objective (NA $\approx$ 0.35). The actual writing speed, achieved with a high-precision positioning system (Aerotech) is 100 mm/min at a pulse energy of 200 nJ and a repetition rate of 100 kHz (Coherent Mira/Reg A). Such waveguides exhibit low propagation losses ( $<$ $0.3$ dB/cm) and almost vanishing birefringence ($\Delta n_{H,V}\approx 10^{-6}$). With a supported mode field diameter of $12$ $\mu$m $\times$ $15$ $\mu$m coupling losses of $3$ dB are obtained with standard single mode fibers.

Photon-pair generation and measurement. We generate photon-pairs at $\lambda=815$ nm using a standard type-I spontaneous parametric down-conversion source with a visibility of 93%. A BiBO crystal is pumped by a $100$ mW, $407.5$ nm laser diode producing horizontally polarized photon-pairs, which are collected by polarization maintaining fibers. Commercial V-groove fiber arrays were used to couple the photons into the chip as well as collecting them at the output facet from the individual waveguides. We used high-NA multi-mode fibers in order to feed the photons to the respective avalanche photo-diodes, ensuring low coupling losses at the output side of the chip. From the data of the photo-diodes, the photon probability distribution at the output, as well as the inter-channel correlations, were computed using a correlation device (Becker-Hickl) and standard computer programs (LabView for the data acquisition and MatLab for the data processing). The photonic chip was fabricated such that the waveguide spacing at the input and output facets matched the standard fiber array spacing of $127\,\mu$m. The typical two-photon coincidence rates after propagating through the chip were approximately $30$Hz, and an integration time of around $20$ seconds was used for coincidence measurements.

Algorithmic reconstruction. The algorithm for reconstructing the density matrix from a measured set of correlations uses the simultaneous perturbation stochastic approximation (SPSA) Spall:1992-332:RAR . We follow an optimization process very similar to Ferrie:2014-190404:PRL , except we adapt it to work with mixed quantum states. This is achieved by defining the density matrix as, $\hat{\rho}(t)=\hat{T}^{\dagger}(t)\hat{T}(t)/\text{Tr}\{\hat{T}^{\dagger}(t)\hat{T}(t)\}$, as in Eq. (4.5) of Ref. James:2001-52312:PRA . Then the optimization algorithm proceeds optimizing $t$, rather than $\psi$ as in previous approaches for pure quantum states. We determine the reconstruction accuracy using the experimentally measured device transfer function. The recovery of 500 random density matrices is simulated, each for 1500 iterations of the algorithm. The reconstruction fidelity, the average fidelity between the true density matrix and the algorithm’s final guess, is found to be $99.67\%$.

II Acknowledgments

We acknowledge support by the Australian Research Council (ARC) (DP130100135, DP160100619); Erasmus Mundus (NANOPHI 2013 5659/002-001); Alexander von Humboldt-Stiftung; Australia-Germany Joint Research Co-operation Scheme of Universities Australia; German Academic Exchange Service.

III Author contributions

J.G.T, A.S.S, and A.A.S developed the theoretical concept. J.G.T, M.G, R.H, and A.S designed and fabricated the photonic chip. J.G.T, M.G, and R.H preformed the experiments. All authors co-wrote the manuscript.

IV Supplementary material

IV.1 Classical characterization of the device transfer function

Using the photonic chip for tomography requires that the action of the chip on an input photon can be accurately predicted. The propagation of a photon wavefunction through the chip can be described by $\mathbf{\hat{U}}\ket{\psi_{in}}=\ket{\psi_{out}}$, here $\mathbf{\hat{U}}$ is the transfer function of the chip, describing the mapping from the input waveguides to the outputs. Thus in order to use the chip for tomography this transfer function needs to be determined. This could be achieved using quantum process tomography, however we used a more straightforward method Heilmann:2015-96:RAR requiring just classical light.

We use a classical laser beam coupled into the chip in order to find its transfer function. The absolute values of each element in the transfer function, $|U_{i,j}|^{2}$ are found simply by coupling the laser into one of the inputs and measuring the intensity at each of the 6 outputs. For example coupling the laser into input waveguide $1$ and measuring the intensity at output waveguide 6 will give the transfer element $|U_{1,6}|^{2}$. Doing this for both the two input ports of the chip gives the absolute values of each of the $12$ elements in the transfer function.

In a similar way the phases of the transfer function can be easily determined, up to a constant overall phase factor between the two columns of the transfer function. To do this we pass the laser through a beamsplitter, splitting it into two beams, then couple each beam into one of the two input waveguides. The two beams interfere within the chip. The output field in the $j^{\text{th}}$ waveguide is be given by,

[TABLE]

Here $E_{1}^{(in)}$ is the complex field amplitude at input waveguide number $1$ and $E_{2}^{(in)}$ is the complex field amplitude at input $2$. Since the absolute values of the transfer matrix elements are already known, it makes sense to separate the transfer matrix into its absolute values and phase components. Thus we rewrite Eq. (1) as,

[TABLE]

Here $\phi_{j}$ is the phase difference between $U_{1,j}$ and $U_{2,j}$, which we need to determine. The global phase $\Delta d/\lambda$ denotes a variable phase shift we will intentionally introduce by increasing the path length of the beam coupled into input 2. Varying the phase of one input relative to the other is ultimately what will allow the unknown transfer matrix phases, $\phi_{j}$, to be determined. The changing phase produces oscillations in the intensity of each output waveguide. These oscillations measured in our chip are shown in Fig. 6(a). The output intensity, $|E_{j}^{(out)}|^{2}$, will oscillate with changing $\Delta d$ according to

[TABLE]

Here $a$ and $b$ are constants related to absolute values of the transfer matrix and input beams. Since $\Delta d/\lambda$ at a given time is the same for all output ports, the phase shifts between transfer matrix elements, $\phi_{j}$, can to be determined, up to an unknown phase shift, $\theta$, between the two columns of the transfer function. This is achieved by comparing the phase shifts between the oscillations of the different outputs in Fig. 6(a). To determine the phase shifts efficiently we plot the output intensity of waveguide $j$ against the intensity of waveguide $k$, which yields elliptical Lissajous curves. The shape and orientation of the ellipses’ determine the phase difference between the oscillations of the intensity in waveguides $j$ and $k$. For example the Lissajous curve for waveguides $2$ and $3$ is plotted in Fig. 6(b). An ellipse is fitted to the data points and the parameters of the ellipse determine the phase difference between the oscillations in outputs $2$ and $3$, thus give the quantity $\phi_{2}-\phi_{3}$. This is repeated for all combinations of output ports to find all the phases of the transfer function, up to an unknown constant phase shift $e^{i\theta}$ between the two columns of the transfer function, as in Heilmann:2015-96:RAR .

After extracting the phases between all the elements in one column of the transfer function we can now write the transfer function up to the constant overall phase factor $e^{i\theta}$,

[TABLE]

This constant phase factor, $\theta$, is unknown because the absolute value of the phase shift between the two inputs $E_{1}^{(in)}$ and $E_{2}^{(in)}$ is never directly measured, all that is known is that the shift is changing with the displacement $\Delta d$. This unknown phase factor could have been determined using an interferometric setup, but we chose not to do this since we can already observe relative phase differences between different density matrices without needing to know the absolute phase.

IV.2 The spatial density matrix of a broadband two-photon wavefunction

In the main work, ‘Scalable on-chip quantum state tomography’, we characterized the spatial density matrix of a number of different two-photon quantum states. Of course these states also have spectral and temporal properties, and these properties have an impact on the spatial density matrices we recover. In this section we show a link between the spectral/temporal properties of the photon-pair and the resulting spatial density matrices. We do this starting from a general frequency dependent two-photon wavefunction of the type that would be produced by our SPDC source, then we integrate over the photon frequencies in order to derive the reduced density matrix of the spatial state of the two-photons. Finally we link this analytic result to real measurements of two-photon states with different temporal delays, demonstrating that our approach to tomography can recover information about spectral/temporal distinguishably of photons.

First we write the frequency dependent two-photon wavefunction as

[TABLE]

Here $\psi_{n,m}(\omega_{1},\omega_{2})$ is the joint spatial and spectral distribution of photons $a$ and $b$. The spatial mode of photon $a$ is denoted $n$ and the mode of photon $b$ is denoted $m$.

For comparison with experiment we now specify the form of the wavefunction to be an anti-bunched state, of the type that would be produced by our SPDC source. Specifically we set

[TABLE]

Giving photons $a$ and $b$ spectra $\phi^{(a)}$ and $\phi^{(b)}$, and putting them in different spatial modes to make an anti-bunched state.

Elements in the reduced spatial density matrix $\hat{\rho}_{n,m,n^{\prime},m^{\prime}}$ can be calculated by projecting the two-photon wavefunction into the operator basis for the spatial density matrix, effectively integrating over frequency to leave a purely spatial description of the density matrix. So we will take inner products with the following operators,

[TABLE]

Where $D$ and $F$ represent some spectral or temporal windows used to distinguish the two photons. In the case where the measurement apparatus is not designed to distinguish different photons we set $D^{\ast}(\omega_{3})D(\omega_{3})=F^{\ast}(\omega_{3})F(\omega_{3})=1/\sqrt{2}$, to preserve normalization since each photon will be detected by both detection windows. Otherwise $D$ and $F$ can be set to be step functions in frequency or time, depending of whether photons are distinguished by their spectrum or arrival times.

So the elements of the anti-bunched state’s density matrix are given by the expectation values,

[TABLE]

This can be evaluated to give,

[TABLE]

[TABLE]

[TABLE]

All other density matrix elements are zero. Thus general form of the anti-bunched state density matrix will be,

[TABLE]

where $I_{a,b}=\int d\omega\,\phi^{(a)\ast}(\omega)\phi^{(b)}(\omega)$ is the spectral overlap between the two photons. Clearly this spectral overlap is an important factor in determining the exact form of the two photon density matrix. It can be equal to zero in cases where the photons have non-overlapping spectrum or if there is a temporal delay between the photons. In the case where the photons are indistinguishable the overlap will be unity, resulting in a different spatial density matrix compared to the non overlapping (distinguishable case) case.

The reduced density matrix of the anti-bunched state [Eq. 12] can be recovered using our static photonic chip based approach to tomography. This can be achieved for quantum states of both distinguishable and indistinguishable photons, despite the fact that we do not use single-photon detectors to directly determine which photon was detected. Ultimately this is possible because of generalized Hong-ou-Mandel (HOM) Hong:1987-2044:PRL interference between the photons in the case where their spectra overlap.

To demonstrate this we use a two-photon source that can be tuned between producing an anti-bunched state of distinguishable or indistinguishable photon pairs. The distinguishably is achieved by introducing variable time delay between the two photons. The quality of this source is checked using traditional HOM interference, showing a visibility of $93\%$ [Fig. 7]. With this photon source we can therefore create two different anti-bunched two-photon states, one state where the two photons are distinguishable because of a time delay between them, and one state where there is no time delay so the photons can’t be distinguished. These correspond to non-degenerate and degenerate two-photon density matrices, respectively, both of which can be described by Eq. 12.

We then demonstrate recovery of both the distinguishable and indistinguishable anti-bunched states. We couple these two different states into our photonic chip, and observe the two distinct sets of two-photon correlations at the six outputs, much like in traditional HOM. For these two cases the two-photon correlations that are observed across the six output waveguides of the device are shown in Fig. 8(a) and 8(b) for the distinguishable and indistinguishable cases respectively. There is a clear difference between these two sets of measurements, especially in the elements $(2,3)$, $(3,4)$, and $(4,5)$, which are much larger in the distinguishable case than the indistinguishable case. This is due to quantum interference as some parts of the two-photon wavefunction are non-orthogonal in the indistinguishable case and thus interfere with one another. Thus we can observe a generalized form of the HOM effect in our device.

As can be seen in Eq. 12 the form of the density matrix we are recovering depends on the overlap between the (complex) spectrum of the photons being measured, so it contains information about how indistinguishable the photons are. Performing reconstruction of the input density matrices of the distinguishable and indistinguishable coincidence measurements from Fig. 8(a) and (b) allows us to see the difference between the density matrices. The real and imaginary parts of the reconstructed density matrix for the distinguishable case is shown in Figs. 9(a) and (b) and the density matrix for the indistinguishable case is shown in Figs. 9(c) and (d). Both reflect the expected form of the reduced density matrix that was derived in Eq. 12. The density matrix of the distinguishable case has small off diagonal elements, suggesting that it is a ‘classical’ state. In contrast the density matrix of the indistinguishable state has two large off diagonal elements, these show that there is quantum coherence between different states of the two particles, and thus the state will exhibit some non-classical statistics.

Therefore, even without photon distinguishing detection schemes, our approach to quantum tomography is capable of recovering distinguishable and indistinguishable multi-photon quantum states using generalized HOM interference.