High-dimension experimental tomography of a path-encoded photon quantum-state

TL;DR

This paper introduces a resource-efficient method for high-dimensional quantum state tomography of path-encoded photons using a rotating optical Fourier Transform, achieving high fidelity in experimental demonstrations.

Contribution

The authors develop and experimentally validate a simple, scalable tomography technique for high-dimensional path-encoded photon states, reducing complexity compared to traditional methods.

Findings

Achieved an average fidelity of 0.9852 in reconstructing six-dimensional states.

Demonstrated the method's effectiveness across various six-dimensional density matrices.

Performance exceeds or matches existing tomographic techniques for similar systems.

Abstract

Quantum information protocols often rely on tomographic techniques to determine the state of the system. A popular method of encoding information is on the different paths a photon may take, for example, parallel waveguides in integrated optics. However, reconstruction of states encoded onto a large number of paths is often prohibitively resource intensive and requires complicated experimental setups. Addressing this, we present a simple method for determining the state of a photon in a superposition of d paths using a rotating one-dimensional optical Fourier Transform. We establish the theory and experimentally demonstrate the technique by measuring a wide variety of six-dimensional density matrices. The average fidelity of these with the expected state is as high as 0.9852 +/- 0.0008. This performance is comparable or exceeds established tomographic methods for other types of systems.

Equations

$$P(k_x)=|\tilde{\psi }(k_x){|}^2({\rho }_{11}+{\rho }_{22}+2|{\rho }_{12}|\cos (L_{12}^x{k}_x+{\varphi }_{12})),$$

$$\begin{array}{c}{F}_{k_x}\{P(k_x)\}(\bar{x})=({\rho }_{00}+{\rho }_{11})\delta (\bar{x})+{\rho }_{12}\delta (\bar{x}-L_{12}^x)+{\rho }_{21}\delta (\bar{x}+L_{12}^x).\end{array}$$

$$\begin{array}{c}{\mathbf{\rho }}_x(k_y)=⟨k_y|\mathbf{\rho }|k_y⟩=\sum _{ij}^d{\rho }_{ij}\text{innerproduct}{k}_y{{\psi }_i}_y\text{innerproduct}{\psi }_j{k}_y{|{\psi }_i⟩}_x⟨{\psi }_j|.\end{array}$$

$$\begin{array}{c}{\mathbf{\rho }}_x(k_y)=|\tilde{\psi }(k_y){|}^2\sum _{ij}^d{\rho }_{ij}{e}^{i{L}_{ij}^y{k}_y}{|{\psi }_i⟩}_x⟨{\psi }_j|.\end{array}$$

$${\mathbf{\rho }}_x(k_y)=\sum _i^d{\rho }_{ii}{|{\psi }_i⟩}_x⟨{\psi }_i|+\sum _{i\ne j}^d{\rho }_{ij}{e}^{i{L}_{ij}^y{k}_y}{|{\psi }_i⟩}_x⟨{\psi }_j|.$$

$$\begin{array}{c}P(x_m,k_y)=⟨x_m|{\mathbf{\rho }}_x(k_y)|x_m⟩=|\psi (0){|}^2\left(\sum _i^d{\rho }_{ii}{\delta }_{x_m{x}_i}+\sum _{i\ne j}^d{\rho }_{ij}{e}^{i{L}_{ij}^y{k}_y}{\delta }_{x_m{x}_i}{\delta }_{x_m{x}_j}\right).\end{array}$$

$$\begin{array}{c}{F}_{k_y}\{P(x_m,k_y)\}(\bar{y})=\int \mathrm{d}{k}_{y}P(x_m,k_y)e^{i{k}_y\bar{y}}=\\ \sum _i^d{\rho }_{ii}{\delta }_{x_m{x}_i}\delta (\bar{y})+\sum _{i\ne j}^d{\rho }_{ij}\delta (\bar{y}-L_{ij}^y){\delta }_{x_m{x}_i}{\delta }_{x_m{x}_j}.\end{array}$$

Full text

I Introduction

Determining the quantum state of a physical system is a key task in quantum physics, in particular in quantum metrology, quantum information and quantum cryptography. Since the quantum state $\boldsymbol{\rho}$ determines all the measurable properties of system, reconstructing it can aid in building and bench-marking devices in the aforementioned areas. Many techniques exist for reconstructing a variety of types of physical quantum systems, ranging from electron spin to atomic position Paris and Sacchi (2003); Lundeen et al. (2009); James et al. (2005); Rippe et al. (2008). Considerable effort has been put towards the development of these experimental techniques, and improving their scalability. These techniques have even found applications outside their original purposes, such as in classical image processing D’Ariano et al. (1996).

Early schemes for quantum information processing with photons proposed encoding information in a set of $d$ paths that a photon could be in Milburn (1989); Reck et al. (1994). These paths are sometimes called rails and usually follow a waveguide or the route a classical beam would take. The quantum state is the superposition of paths that the photon is in. If two paths are used, the state functions as qubit. Crucially, unlike a commonly used photon qubit, polarization, higher dimensional states can be formed by adding more paths. This allows more information to be encoded on a single photon, which is then termed a ‘qudit’ Bechmann-Pasquinucci and Tittel (2000). Higher-dimensional encodings can also simplify the design of quantum logic circuits Lanyon et al. (2009). Path encoding is an especially popular alternative to polarization encoding for on-chip photonic devices due to the difficulty of fabricating devices that maintain polarization Matthews et al. (2009). A seminal contribution to manipulating path states was a scheme to implement any chosen discrete $d\times d$ unitary using an interferometric array of beamsplitters and phaseshifters acting on $d$ paths Reck et al. (1994); Clements et al. (2016). On-chip implementations of the above-mentioned universal discrete unitaries are an especially active research direction Metcalf et al. (2013); Peruzzo et al. (2011) due their ability to implement a wide set of quantum information algorithms and networks with a single device. In free-space however, paths are not a commonly used method of state encoding. This is in part due to the difficulty of reconstructing the quantum state encoded in a large number of paths.

Despite the importance of the path-state encoding, relatively little work has been conducted on methods to reconstruct path states with more than two dimensions. One possibility is to use a universal unitary to rotate between a complete set of incompatible bases, making projective measurements in each. This is quantum state tomography. If one universal unitary implements a quantum information algorithm another would be required for the tomography, thus doubling the complexity of the on-chip device. Moreover, in free-space, the interferometer required for the universal unitary has never been constructed due to its complexity, phase instability and alignment issues. Another possibility is to interfere each and every possible pair of paths. This has the drawback of requiring $\mathcal{O}(d^{2})$ switchable elements such as beamsplitters and mirrors. As far as we know, neither strategy for $d$-path quantum state reconstruction has been implemented above $d=2$.

On the other hand, schemes have been devised and implemented to reconstruct the full spatial quantum state of a photon Bamber and Lundeen (2014); McAlister et al. (1995); Stoklasa et al. (2014); Smith et al. (2005). The $d$-path state is embedded in this larger and continuous Hilbert space, that of the photon’s transverse position. Consequently, one might consider using these techniques to reconstruct the $d$-path state. However, these methods provide much more information than is required, for example, the modal distribution of each path. For a fixed measurement time, one would expect these methods to have a higher noise per point than a method that begins with $a$ $priori$ information that limits the size of the state-space. A method that only determines the parameters of the $d$-path density matrix $\boldsymbol{\rho}$ would outperform these full spatial state reconstructions and may be less experimentally complicated.

We propose such a method for reconstructing a $d$-dimensional path encoded state, and implement it using both a classical and non-classical source of light. The method relies on the fact that the discrete set of states is embedded in the aforementioned continuous Hilbert space, that of the photon’s transverse position. This opens up the possibility of having many more measurement outcomes than $d$. Hence, projections in many incompatible bases can be performed with one measurement apparatus configuration. This is accomplished by either measuring position $x$, with a camera, or momentum $k$, by adding a cylindrical lens that performs a one-dimensional optical Fourier transform (OFT). When the paths are distributed in both the $x$ and $y$ transverse dimensions the lens rotates to a finite set of angles in the $xy$ plane. A discrete computational Fourier transform (DFT) of the measured momentum probability distribution then allows us to easily reconstruct high-dimensional path-encoded states. In classical optics, this is equivalent to measuring the optical coherency matrix Kagalwala et al. (2015); Horne et al. (1989).

II Theory

Before introducing our method, we begin by rigorously defining what is meant by the path-state encoding. A single photon is travelling along $z$ and has transverse position $x$ (for clarity we restrict ourselves to one transverse dimension for now). A set of $d$ states $\{\ket{\psi_{i}}_{x}\}$ are defined in this $x$-space. Each has a narrow position distribution $|\psi_{i}(x)|^{2}=|\innerproduct{x}{\psi_{i}}_{x}|^{2}$ and is displaced in $x$ from the other states. The displacement is much larger than the width of the any state distribution, which ensures that they have minimal overlap, $\innerproduct{\psi_{i}}{\psi_{j}}_{x}\approx\delta_{ij}$ (the Kronecker delta). Crucially, this means they form an approximate orthonormal set. In turn, this set is a basis that defines a discrete subspace embedded in the continuous $x$ Hilbert space. Each basis state $\ket{\psi_{i}}_{x}$ represents the photon being in a particular single path $i$. Our method additionally requires that the path-states have identical spatial distributions, other than the displacement. In the momentum basis, a displacement only affects the phase of the distribution. Hence, this identicality ensures that $|\innerproduct{k_{x}}{\psi_{i}}_{x}|^{2}=|\widetilde{\psi}(k_{x})|^{2}$ for all $i$, i.e. the states have identical probability distributions in momentum $k_{x}$. Our goal is to reconstruct states in this $d$-path subspace. Specifically, we will determine the density matrix $\boldsymbol{\rho}=\sum_{ij}^{d}\rho_{ij}\ket{\psi_{i}}\bra{\psi_{j}}_{x}.$

We now describe the basic concept behind our reconstruction method. Each diagonal element of the density $\rho_{ii}$ is the probability for the photon to be in the corresponding path $i$. These diagonals are straight-forward to determine by directly projecting onto the path states with $d$ detectors or by measuring the position $x$ distribution with a camera. The off-diagonal elements are more challenging. We call these elements ‘coherences’ since a given element $\rho_{ij}=|\rho_{ij}|\exp(\phi_{ij})$ describes the coherence between paths $i$ and $j$. These two paths play the same role as the slits in a double-slit interferometer. There, the coherence between the two slit-paths sets the resulting visibility (i.e. contrast) of the resulting interference pattern seen on a distant screen. Effectively, this pattern is the transverse momentum $k_{x}$ probability distribution of the photon just after the slits. Distinguishability, entanglement with other systems, amplitude imbalance, and technical or environmental noise can all decrease the visibility and, hence, the magnitude $|\rho_{ij}|$. The phase difference between a photon passing through slit one and two, equivalent to $\phi_{ij}$, sets the transverse offset of the interference pattern. In summary, if one can measure the visibility and phase of the interference between the paths $i$ and $j$, one can determine $\rho_{ij}$ for the path-state.

The path states can be interfered by looking in transverse momentum (i.e. Fourier) space. Rather than using a distant screen, this can be easily accomplished with a lens of focal length $f$. The transverse position $x_{f}$ of the photon $f$ distance $after$ the lens is proportional to the transverse momentum $k_{x}$ at $f$ distance $before$ the lens. The exact relationship is $x_{f}=f\lambda k_{x}/2\pi$, where $\lambda$ is the wavelength of the photon. In this way, a lens performs an optical Fourier transform. A general state $\rho(x,x^{\prime})=\bra{x}\boldsymbol{\rho}\ket{x^{\prime}}$ in the continuous position basis corresponds to $\rho(k_{x},k_{x}^{\prime})=\bra{k_{x}}\boldsymbol{\rho}\ket{k_{x}^{\prime}}$ in the momentum basis. A camera at $f$ after the lens will thus record a signal proportional to the momentum probability distribution $P(k_{x})=\rho(k_{x},k_{x})$. In the recorded camera distribution (i.e. image), two or more interfering paths will form an interference pattern composed of fringes.

With an illustrative two-path example, we now explicitly show how we extract $\rho_{ij}$ from the measured momentum distribution. Our two paths, $\ket{\psi_{1}}_{x}$ and $\ket{\psi_{2}}_{x}$ are centered at $x=x_{1}$ and $x=x_{2}$, respectively. Utilizing the Fourier shift theorem, we get an explicit form for the interference pattern in terms of the path-state density matrix $P(k_{x})=|\widetilde{\psi}(k_{x})|^{2}\sum_{i,j}^{d}{\rho}_{ij}\exp(iL^{x}_{ij}k_{x})$, where $L^{x}_{ij}=|x_{i}-x_{j}|$ is the pairwise distance between the paths $i$ and $j$ in the $x$ direction. Using the Hermiticity of $\boldsymbol{\rho}$, this can be simplified to

[TABLE]

where $\phi_{12}$ is the phase of $\rho_{12}$. This is the interference pattern; the fringes are due to the oscillating cosine term, and the envelope is given by Fourier Transform of a single path’s spatial distribution. The visibility of an interference pattern is defined as the ratio of amplitude of the oscillation to its average value. Because $\mathrm{Tr}\{\boldsymbol{\rho}\}=1$, the visibility can be found to be $2|\rho_{ij}|/(\rho_{11}+\rho_{22})=2|\rho_{ij}|$. As we motivated earlier, the visibility is directly proportional to the magnitude of $\rho_{ij}$ and the phase $\phi_{ij}$ of the cosine is the phase of $\rho_{ij}$.

We can extract the magnitude and phase information from the interference pattern in Eq. 1 by using a discrete Fourier transform, performed with a computer. For clarity, here we present the method in terms of the continuous Fourier transform $F_{k_{x}}$, where the subscript denotes the transform is along the $k_{x}$-axis. The Fourier transform of the measured momentum distribution gives

[TABLE]

Because this is the Fourier transform of momentum, the resulting function is back in terms of position. However, whereas a Fourier transform of the amplitude distribution would be in terms of $x$, this is transform of the momentum probability distribution. We distinguish the resulting position variable by using a bar, i.e. $\bar{x}$. With this calculated distribution, it becomes simple to find $\rho_{12}$. Namely, the complex value of the Fourier transform at a position equal to spacing $L^{x}_{12}$ is $\rho_{12}$. That is, $\rho_{12}=F_{k_{x}}\{P(k_{x})\}(\bar{x}=L^{x}_{12}).$ In short, a Fourier transform of the measured interference pattern directly gives the phase and magnitude of off-diagonal’s such as $\rho_{ij}$.

While interfering paths using a lens works well for two paths, extending to higher dimensions by adding more paths introduces the possibility of overlapping peaks in $F_{k_{x}}\{P(k_{x})\}(\bar{x})$. If more than one pair of paths have the same spacing they would contribute to the value of the Fourier transform at the same $\bar{x}$. This would make it impossible to distinguish the contributions from their corresponding density matrix elements. For instance, they may cancel completely if completely out of phase. We now discuss this problem and present our solution.

II.1 Extension to more paths

One way to add more path states would be to simply distribute them along $x$. In many cases however, it may be convenient to instead distribute them in both transverse directions $x$ and $y$. We begin by defining our path-states in two dimensions. As before, all the path states have identical position wavefunctions other than a displacement. Specifically, the $i$-th path-state $\ket{\psi_{i}}_{xy}$ is centered at position $\mathcal{P}_{i}=[x_{i},y_{i}]$. We call the set of all positions $\mathcal{P}_{i}$ the ‘geometry’. For simplicity, we assume that each state has a position wavefunction that is a two-dimensional Gaussian with equal width $\sigma$ in both the $x$ and $y$ directions. This ensures that the distribution is symmetric under rotations in the $xy$ plane and, also, that it is a product state: $\ket{\psi_{i}}_{xy}=\ket{\psi_{i}}_{x}\ket{\psi_{i}}_{y}$. The subscripts $x$ and $y$ distinguish the two Hilbert spaces. The position wavefunctions are, thus, $\psi_{i}(x,y)=\psi(x-x_{i})\psi(y-y_{i})$. As before, the paths must be separated by much more than $\sigma$ so that they have minimal overlap, $\innerproduct{\psi_{i}}{\psi_{j}}_{x}\approx\delta_{ij}$ and $\innerproduct{\psi_{i}}{\psi_{j}}_{y}\approx\delta_{ij}$. Again, this means the states form an approximate orthonormal basis for the embedded discrete Hilbert space. Our goal is to reconstruct the density matrix of the full $d$-path quantum state, $\boldsymbol{\rho}=\sum_{ij}^{d}\rho_{ij}\ket{\psi_{i}}_{xy}\langle\psi_{j}|$.

We must now carefully label and then distinguish every possible pair of paths since each is associated with a coherence $\rho_{ij}$. The positions of two paths, $i$ and $j$, are a connected by the line segment $\mathcal{L}_{ij}$ with end points $(\mathcal{P}_{i},\mathcal{P}_{j})$. The spacing of the path-pair is given by the length of the line segment, $L_{ij}=\sqrt{\left(L^{x}_{ij}\right)^{2}+\left(L^{y}_{ij}\right)^{2}}$, where $L^{s}_{ij}=$ $|s_{i}-s_{j}|$ is the path spacing along the $s=x,y$ axis. Just as in the simple example in the last section, we obtain the coherences $\rho_{ij}$ from the recorded signal in momentum space, $P(k_{x},k_{y})$. However, a $d$-dimensional context allows for the possibility of multi-path interference. If each line segment $\mathcal{L}_{ij}$ has a unique length and/or orientation, then each path-pair will create a signal at a distinct oscillation frequency and/or orientation in the interference pattern. In this case, all the coherences $\rho_{ij}$ can just be read out from the DFT of the two-dimensional interference pattern in $P(k_{x},k_{y})$. The latter could be measured by using a spherical lens to perform the optical Fourier Transform in both the $x$ and $y$ directions. The required path geometry for this, called a Golomb Rectangle, is discussed in the Appendix. However, consider if in our geometry there are two parallel line segments $\mathcal{L}_{ij}$ and $\mathcal{L}_{i^{\prime}j^{\prime}}$ with the same length, $L_{ij}=L_{i^{\prime}j^{\prime}}$. In this case, the contribution to the DFT from $\rho_{ij}$ and $\rho_{i^{\prime}j^{\prime}}$ can not be distinguished since they will appear at same peak position $[\bar{x}=L^{x}_{ij},\bar{y}=L^{y}_{ij}]$.

A cylindrical lens solves this problem so long as $\mathcal{L}_{ij}$ and $\mathcal{L}_{i^{\prime}j^{\prime}}$ are not both segments of the same line. As an example, in Fig. 2a, $\mathcal{L}_{16},\mathcal{L}_{25},$ and $\mathcal{L}_{34}$ are parallel and of equal length. As explained above, using a spherical lens will fail in this case. Consider instead the action of cylindrical lens with its OFT-axis parallel to these example line segments $\mathcal{L}_{ij}$. The resulting probability distribution will be of momentum in the direction of $\mathcal{L}_{ij}$ and position in the perpendicular direction (Fig. 2b). Critically, the line segments, $\mathcal{L}_{16},\mathcal{L}_{25},$ and $\mathcal{L}_{34}$, have distinct positions along the axis perpendicular to $\mathcal{L}_{ij}$. The cylindrical lens will leave these positions unchanged. Consequently, the interference oscillation corresponding to each line segment will appear at a distinct position and, hence, will be distinguishable. In our example, we will therefore be able to measure $\rho_{16}$, $\rho_{25}$, and $\rho_{34}$ by taking a one-dimensional DFT of the recorded camera image along cylindrical lens axis. This concept can be extended to other coherences of $\boldsymbol{\rho}$ by rotating the cylindrical lens and, hence, OFT-axis to all the unique angles $\theta_{ij}$ in our full set of line segments, where $\theta_{ij}=\tan^{-1}((y_{i}-y_{j})/(x_{i}-x_{j}))$. An example of the recorded image for one such angle is shown in Fig. 2c-d. So long as no two $\mathcal{L}_{ij}$ and $\mathcal{L}_{i^{\prime}j^{\prime}}$ are segments of the same line, this procedure can be used to determine the density matrix.

For clarity in our mathematical description of this procedure, instead of rotating the OFT we rotate the full quantum state $\rho(x,x^{\prime})$ in position space. Because our path wavefunctions $\psi_{i}(x,y)$ are rotationally invariant, this task is reduced to finding the rotated path positions $\mathcal{P}_{i}^{\prime}$. If $R(\theta)$ is the standard rotation matrix then $\mathcal{P}_{i}^{\prime}=R(\theta)\mathcal{P}_{i}$. After this rotation, the wavefunctions are $\psi^{\prime}_{i}(x,y)=\psi(x-x_{i}^{\prime})\psi(y-y_{i}^{\prime}).$ In the following discussion, we drop the primed notation for clarity, assuming that the paths have been rotated to the appropriate angle $\theta$. Whatever the angle, the OFT is always along the $y$-axis and the $x$-axis is left untransformed. An example of this is shown in Fig. 2b for $\theta_{34}=90^{\circ}$. The camera now records a signal proportional to probability distribution $P(x,k_{y})$.

We consider first the action of the OFT along the $y$ direction. The result is $\boldsymbol{\rho}_{x}(k_{y})$, which is still an operator in $x$-space, as indicated by the subscript $x$:

[TABLE]

Using the Fourier shift theorem, $\innerproduct{k_{y}}{\psi_{i}}_{y}=F\{\psi(y-y_{i})\}=\widetilde{\psi}(k_{y})\exp(iy_{i}k_{y})$, we simplify this to the following expression:

[TABLE]

We will drop the envelope $|\widetilde{\psi}(k_{y})|^{2}$ for the sake of brevity in the math below, but in practice it sets the width and overall height scale of the peaks in the DFT.

To analyze the sum in Eqn. 4, first we separate the terms where $i=j$,

[TABLE]

We now find the position and momentum probability distribution $P(x,k_{y})$. We evaluate it at $x=x_{m}$, where $x_{m}$ is drawn from the set of positions $\mathcal{P}_{i}$ of the paths in the rotated geometry. If $\ket{\psi_{i}}_{x}$ is not centered on $x_{m}$, then the wavefunction is nearly zero at $x=x_{m}$. In other words, $\innerproduct{x_{m}}{\psi_{i}}_{x}=\psi(x_{m}-x_{i})=\psi(0)\delta_{x_{m}x_{i}}$. Note that even though the values of $x_{m}$ may not be integers, they are from a discrete set and as such we use the Kronecker delta. With this, the probability distribution is

[TABLE]

For notational convenience, we drop the $|\psi(0)|^{2}$ factor. In practice, this factor is accounted for when the recorded camera images are normalized.

As in our two-path example in the last section, we now take a Fourier transform of $P(x_{m},k_{y})$. We do so along only the $k_{y}$ direction. In practice, this is implemented by performing a DFT of the camera image. The result is

[TABLE]

There are two cases to consider: Case 1. For a chosen $x_{m}$, one or more pairs of paths have $x_{i}=x_{j}=x_{m}$. If every pair of these has a unique $L^{y}_{ij}$ then the peak at $F_{k_{y}}\{P(x_{m},k_{y})\}(\bar{y}=L^{y}_{ij})$ has value $\rho_{ij}$. Case 2. For a chosen $x_{m}$, only a single path in the rotated frame has $x_{i}=x_{m}$. In this case, the second term is zero and the first sum gives $F_{k_{y}}\{P(x_{m},k_{y})\}(\bar{y}=0)=\rho_{ii}$, a diagonal element of the density matrix. Together, these two cases allow us to determine both the diagonal and off-diagonal elements. Repeating this for the full set of cylindrical lens angles $\theta_{ij}$, the full $d$-path density matrix can be determined. This method works as long as no two parallel line segments $\mathcal{L}_{ij}$ between paths with equal spacing are along the same line, which is equivalent to the condition in Case 2. The above theoretical description constitutes the method to determine $\boldsymbol{\rho}$ that we will later demonstrate experimentally.

III Experiment Setup

We demonstrate our technique by reconstructing density matrices produced by the setup shown in Fig. 3. Using polarized light from a diode laser at $\lambda=808$ nm, we use a series of waveplates and calcite beam displacers, with the crystal axis cut at 45*∘*, to generate $d=6$ paths that contain state $\boldsymbol{\rho}$. The beam displacers utilize birefringence to split the ordinary and extraordinary polarizations of light, producing the paths, while offering robust phase stability. A 4$f$ system ensures that all optics are within the Rayleigh length of the beam. The beam waist is $340$ $\upmu$m. The first beam displacer produces a transverse shift of $\delta_{x}=$ 2.7 mm. The crystal axis is orientated such that horizontally polarized light is shifted horizontally. The second crystal shifts vertically polarized light vertically by $\delta_{y}=$ 2.7 mm as well. Lastly, the third crystal shifts horizontally polarized light by $\delta_{X}=$ 4 mm horizontally. In total, the beam displacers produce eight paths. However, this geometry is incompatible with the tomography method. For example, $\mathcal{L}_{13}$ and $\mathcal{L}_{57}$ are both segments of the same line and of equal length. Therefore, we block two of the eight paths to obtain a two by three path geometry which meets all the conditions necessary for the reconstruction method to work. The resulting six-dimensional density matrix is a function of the half-waveplate (HWP) and quarter-waveplate (QWP) angles $\phi,\zeta,$ and $\Omega$, preceding each crystal, shown in Fig. 3a. As such we write that $\boldsymbol{\rho}=\boldsymbol{\rho}(\phi,\zeta,\Omega)$. Note that in this situation a 2D OFT could not be used to reconstruct the state as the geometry due to repeating spacings of parallel line segments, such as, $\mathcal{L}_{12},\mathcal{L}_{34},$ and $\mathcal{L}_{56}$ (among others), as shown in Fig. 3a.

The spatial distributions are imaged onto a camera using a 4$f$ system where $f_{1}=1000$ mm and $f_{2}=400$ mm, and background subtracted. The OFT is performed using a cylindrical lens ($f=250$ mm), placed one focal length before the camera. The camera has a resolution of $3088\times 2076$ and each pixel is square and has side-length $\gamma=2.40$ $\upmu$m. The lens rotates in an automated mount to each of the right angles $\theta_{ij}$. The one dimensional OFT produces sets of fringes, one for each pair of paths for which the angle between them is $\theta_{ij}$. To obtain the coherence, first we rotate the image back by $-\theta_{ij}$ so as to orient the interference patterns along $y$. We then take a single line of pixels parallel to $y$ and through the center of each interference pattern. We perform a discrete Fourier transform of the recorded intensity distribution along these pixels. Alternatively, one could average over a wider section of each interference pattern. However, averaging could potentially lower the visibility if the fringes are skewed, which could occur due to spherical aberration or misalignment of the lens.

The complex amplitude $S_{ij}$ of a peak in the DFT at spacing $\bar{y}=L^{y}_{ij}$ is proportional to $\rho_{ij}$. We define a scaling factor so that $|\rho_{ij}|=|S_{ij}/S|$, where $S=\sum_{i}^{d}S_{ii}$, is the total recorded intensity. We can obtain $S$ directly from the images used to calculate $\rho_{ij}$, without having to change the apparatus. Specifically, we note that in Eq. 7 the following holds: $\sum_{x_{m}}\sum_{i}^{d}\rho_{ii}\delta_{x_{m}x_{i}}=1$. In terms of the recorded signal, we get $\sum_{x_{m}}\sum_{i}^{d}S_{ii}\delta_{x_{m}x_{i}}=S$. That is, the sum of the $\bar{y}=0$ peaks for each $x_{m}$ gives the proper normalization constant $S$. This is convenient as no extra data is required to obtain this scaling.

The complex peaks also determine the phase of coherence $\rho_{ij}$ up to an overall convention and a constant offset. The offset could also be found by determining the location of $y=0$ in the camera image. However, this can be experimentally difficult. Instead, we input a known quantum state, reconstruct it, and use it as a reference. A second reference quantum state sets our convention for which direction along the $y$ axis is positive. This is equivalent to phase-conjugation or, equivalently, which peak height $\rho_{ij}$ is proportional to in the DFT, $\bar{y}=L^{y}_{ij}$ or $\bar{y}=-L^{y}_{ij}$ . Any two states whose coherences are different and non-zero may be used as a reference.

IV Results

The reconstruction method is demonstrated with $d=6$ dimensional states of the form $\boldsymbol{\rho}(\phi=22.5^{\circ},\zeta,\Omega=22.5^{\circ})=\boldsymbol{\rho}(\zeta)$ (i.e., all HWP are held fixed), whose coherences are a function of the QWP angle $\zeta$. Fig. 4a shows the real and imaginary parts of the fifteen off-diagonal elements $\rho_{ij}$ (not including the conjugates, $\rho_{ij}$). Note that, several elements of $\rho_{ij}$ will be equal since the set of states the experimental setup can produce is constrained. We see good agreement between the experimental values (dots) and the theoretical values (dotted lines). Fig. 4b shows the experimental and theoretical density matrix for $\zeta=30^{\circ}$ in a side by side comparison.

Next, we calculate the state fidelity, $F(\boldsymbol{\rho},\boldsymbol{\sigma})=\mathrm{Tr}\{\sqrt{\sqrt{\boldsymbol{\rho}}\boldsymbol{\sigma}\sqrt{\boldsymbol{\rho}}}\}$, which is a measure of how close two density matrices $\boldsymbol{\rho}$, and $\boldsymbol{\sigma}$ are. The fidelity is bound between zero and unity, and $F(\boldsymbol{\rho},\boldsymbol{\sigma})=1$ only if $\boldsymbol{\rho}=\boldsymbol{\sigma}$. We calculate the fidelity between the experimentally reconstructed state $\boldsymbol{\rho}(\phi,\zeta,\Omega)$ and the theoretically predicted one $\boldsymbol{\rho}_{\mathrm{th}}$. This is done for sets of states where one of the waveplate angles ${\phi}$, ${\Omega}$, or ${\zeta}$ are varied. Specifically, if a HWP is not being varied, then it is fixed at 22.5*∘*. If either HWPs ${\phi}$ or ${\Omega}$ are varied then we fix $\zeta=45^{\circ}$. Fig. 5a displays the results. The average fidelity when ${\phi}$ is varied is $0.987\pm 0.001$, $0.9893\pm 0.0008$ when ${\Omega}$ is varied, and $0.979\pm 0.001$ when ${\zeta}$ is varied. Average over all trials the fidelity is $0.9852\pm 0.0008$. Given the fidelity is consistently close to one, we conclude that the reconstructed state faithfully reconstructs the theoretical state.

The density matrix generalizes the concept of the state by accommodating statistical mixtures of states. A measure of this is the state purity Tr$\{\boldsymbol{\rho}^{2}\}$. Pure states are those for which the purity is one. Statistical mixtures decrease the purity to a minimum of $1/d$, where $d$ is the dimension of the Hilbert space. We reconstruct $\boldsymbol{\rho}$ and calculate the purity. To mix the state we introduce a series of HWP-QWP-HWPs-QWP, shown in Fig. 4a. Here the QWP are 90*∘* apart. The subscript $s$ on HWPs indicates that it is rapidly spinning in a mount. HWPs spins faster than the collection time of the camera. This mixes the polarizations, which in turn mixes the path state density matrix. The first HWP controls the degree of mixing. The waveplate angle $\tau=0^{\circ}$ produces a pure state and $\tau=22.5^{\circ}$ produces a totally mixed polarization state. The theoretical purity of the resulting path state $\boldsymbol{\rho}(\tau)$ is Tr$\{\boldsymbol{\rho}^{2}(\tau)\}=1/9(5+4\cos^{2}{(4\tau)})$. Fig. 5b shows that the measured purity agrees strongly with the theory.

To demonstrate that the technique is applicable to non-classical sources of light we switch to a source of down converted photons. Along with the significantly lowered average photon count number, single photon sources have the added challenge that the coherence length is also substantially reduced. To accommodate for this fact, we introduce glass slides in the paths of the photons to compensate for the path length difference imparted by the displacement crystals. We pump a 15 mm ppKTP crystal with a 404 nm laser to produce degenerate pairs at 808 nm. The measured autocorrelation $g^{(2)}$ of the source is 0.1979 $\pm\;0.0005$, which is much less than one thereby confirming the light is antibunched and non-classical. We generate a 4$\times$4 path encoded density matrix in a square geometry with each side length being 2.7 mm. We replace the QWP labeled by $\zeta$ with a HWP, and set all HWP at 22.5*∘*. The purity of this matrix also depends on $\tau$. The experiment is otherwise the same, except that the down-converted photons are imaged onto an electron multiplying CCD (EMCCD). The EMCCD has a resolution of $512\times 512$ with a pixel length of $\gamma=16$ $\upmu$m. The calculated purity is shown in Fig. 5b. Because the path compensation is not perfect, some elements have decohered, leading to reduced measured purity. Nevertheless, the results follow the theoretical trend well. In Fig. 5c we plot the density matrix elements for $\tau=22.5^{\circ}$, and list the experimental values. For reference, every theoretical elemlent is either 0.25 or 0. The calculated fidelity is 1.00 $\pm 0.03$, which confirms that the method works well for quantum light sources such as single photons.

The reconstruction method satisfies the unity trace and hermiticity conditions, but, like in most linear reconstruction methods, the state may not be positive semi-definite, as is required by a physically possible density matrix. Positive semi-definitivity requires that $|\rho_{ij}|\leq\sqrt{\rho_{ii}\rho_{jj}}.$ Indeed, experimental imperfections cause calculated fidelity of the single photon density matrix, 1.00 $\pm 0.03$, attains values greater than one when the uncertainty is added. Imposing positive semi-definitivity would require maximum-likelihood algorithms to fit the most likely physical state Baumgratz et al. (2013). Nevertheless, we see that even without these likelihood algorithms, metrics such as fidelity and purity, are still very high. In the Appendix, we discuss the experimental resource requirements of the present tomographic method, as well explicit algorithms for constructing compatible geometries of path states.

V Summary

The above results show the reconstruction method faithfully reconstructs several different metrics for a six dimensional density matrix. The reconstruction method exploits the fact that discrete paths are embedded in a larger, continuous space. This allows us to use Fourier transforms in a continuous space to obtain data for a discrete density matrix. To avoid having multiple signals with the same spacing overlap, we utilize a cylindrical lens to separate the interference patterns spatially along the orthogonal direction. In this sense, the method uses both position space and $k$-space simultaneously to reconstruct the state. We emphasize that the above results reflect a proof of principle; the tomography method presented above can, in theory, accommodate much larger dimensions, given they conform to the physical constraints (e.g. pixel size, camera size, aberration etc.) discussed above. We demonstrate the method experimentally by calculating several metrics for six dimensional matrices. In summary, we present an experimentally simple method using one dimensional Fourier transforms, one optical and the other digital, to reconstruct large dimensional density matrices encoded in the paths of photons.

Acknowledgements

This work was supported by the Canada Research Chairs (CRC) Program, the Canada First Research Excellence Fund (CFREF), and the Natural Sciences and Engineering Research Council (NSERC).

Appendix

Counting the resources required for the method

We briefly consider the resources needed for the method. This allows one to compare it to other quantum-state reconstruction methods such as tomography. We begin with the required number of experimental configurations $\eta$, which are sometimes called ”measurement settings”. In our method, the only experimental change from one measurement to the next is the OFT-axis angle. The number $\eta$ of angles $\theta_{ij}$ required to determine all the density matrix off-diagonal elements $\rho_{ij}$ depends on the specific geometry. In the worse case, each lens angle only retrieves one element. Thus the maximum number of angles would be $d(d-1)/2$, the number of independent off-diagonal complex elements in $\boldsymbol{\rho}$. More typically, a number of paths will lay along the same line $\mathcal{L}$. In this case, multiple off-diagonals can be found with a single lens angle aligned along $\mathcal{L}$. Consequently, judicious choice of geometry can reduce the total number of measurement settings significantly.

Consider a one-dimensional geometry. That is, the path-states are arranged along, say, the $y$-axis. In this case, only one experimental setting is needed to determine the all the density matrix off-diagonals. The lens axis is set parallel to the line of path-states along $y$. The diagonals can be found with one additional setting, the lens axis along $x$. Valid one-dimensional geometries are given by a ”Golomb ruler”, a set of points with no repeated spacings. An optimal Golomb ruler is one that minimizes the maximum spacing for a number of points, $d$. In the limit of large $d$, the Erös-Turan construction for an asymptotically optimal set is $x_{i+1}=L_{\mathrm{min}}(2di+(i^{2}\,\bmod\,d)),$ for $i=0,2,...d-1$ for any value of $d$ that is and odd prime Erdös and Turán (1941). Here, $L_{\mathrm{min}}$ is the minimum possible spacing between paths. The maximum spacing, $L_{\mathrm{max}}=x_{d}-x_{1}$, is then given by $L_{\mathrm{max}}=L_{\mathrm{min}}(2d(d-1)+1)$. Space is a resource for both the path-state encoding and our reconstruction method. The space available (e.g. the size of the camera screen) constrains the maximum possible number of path-states. Considering a one-dimensional arrangement, this maximum is $d_{\mathrm{max}}=L_{\mathrm{max}}/L_{\mathrm{min}}$. In contrast, a Golomb ruler uses space less efficiently: $d=\sqrt{L_{\mathrm{max}}/(2L_{\mathrm{min}})}=\sqrt{d_{\mathrm{max}}/2}$ for large $d$. A one dimensional geometry is optimal for the number of measurement settings required to reconstruct $\boldsymbol{\rho}$ (i.e., $\eta=2$), but sub-optimal in its use of space.

We now examine whether a two-dimensional geometry uses space more efficiently. Robinson gave a generalization of the Golomb ruler called the Golomb rectangle Robinson (1979). In the area of astronomical inteferometers, these are known as ”non-redundant configurations” Kopilovich (2013). In it, every line segment $\mathcal{L}_{ij}$ is associated with a unique vector, i.e. orientation and length. Given a non-redudant configuration, a spherical lens is sufficient to reconstruct $\boldsymbol{\rho}$. As with the one-dimensional geometry, two measurement settings are sufficient: lens present and lens absent. Given $L_{\mathrm{max}}$ and $L_{\mathrm{min}}$, explicit computational searches for optimal Golomb rectangles suggest that maximum number of path-states is $d=\mathcal{O}((L_{\mathrm{max}}/L_{\mathrm{min}}))$ Kopilovich (2013). In two dimensions, the maximum number of closely packed path-states is $d_{\mathrm{max}}=\mathcal{O}((L_{\mathrm{max}}/L_{\mathrm{min}})^{2})$. Consequently, the optimum Golomb rectangles achieve $d=\mathcal{O}(\sqrt{d_{\mathrm{max}}})$. In terms of the total dimension of the discrete path-state Hilbert space, the benefit of going from one to two dimensions is limited to the an overall constant factor.

A number of experimental factors will constrain the largest feasible $L_{\mathrm{max}}$. The most obvious is the range of positions $D$ over which OFT lens system correctly performs a Fourier Transform. This will be set by lens aberrations and/or the lens diameter. Less obvious is the pixel pitch of the camera $\gamma$, the distance between pixel centers. The finest interference pattern produced by the paths will have a fringe spacing of $2\pi L_{\mathrm{max}}/\lambda f$ at the camera, where $f$ is the focal length of the OFT lens. In order to sample the fringes without aliasing, the Nyquist condition then constrains $L_{\mathrm{max}}<\lambda f/\pi\gamma$. On the other hand, $L_{\mathrm{min}}$ is constrained primarily by the width of the paths.

Another resource to consider is the number of measurement outcomes $N$ per measurement setting. In our case, this is the number of camera pixels. As described above, the finest pixel pitch required is $\gamma=\lambda f/\pi L_{\mathrm{max}}$ . However, the camera pixel-array width $W$ must be large enough to record the largest fringe spacing, which in turn, results from the smallest spacing $L_{\mathrm{min}}$. This sets $W=\lambda f/\pi L_{\mathrm{min}}$. Together, these requirements determine the number of pixels along one camera dimension, must satisfy $N>W/\gamma=L_{\mathrm{max}}/L_{\mathrm{min}}=d_{\mathrm{max}}$. As discussed above $L_{\mathrm{max}}<D$. In addition, $L_{\mathrm{min}}$ must be greater than the path-state position-width, $L_{\mathrm{min}}>\sigma$. Taken together, we have the requirement $N>D/\sigma$. Or, using the results from the optimal Golomb ruler above, $N>2d^{2}$ for large $d$. In typical quantum state tomography, $N=d$. These arguments can be straightforwardly extended to two dimensions. This shows that, not only is our method sub-optimal in its use of space, but also in terms of the number of measurement outcomes per measurement setting. Of course, to its credit, it is, so far, the only experimentally feasible method for path-state reconstruction and it only requires two settings.

Perhaps these resources could be traded for one another by using a cylindrical lens. The use of a cylindrical lens permits a denser usage of space since, as long as they are not collinear, two or more line segments may have identical vectors. It follows that the ratio of $d$ to $N$ and $d_{\mathrm{max}}$ will improve. The trade-off is that multiple lens angles and, thus, measurement settings $\eta$ , will be necessary. As far as we know, the required path-state geometry has not been studied and, so, optimal solutions do not exist. Consequently, we are unable to evaluate this resource trade-off in more detail.