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

D. Curic; L. Giner; J. S. Lundeen

doi:10.1364/PRJ.7.000A27

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, e.g., 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 \pm 0.0008

. This performance is comparable to or exceeds established tomographic methods for other types of systems.

1. 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 $ρ$ determines all the measurable properties of a system, reconstructing it can aid in building and benchmarking devices in the aforementioned areas. Many techniques exist for state reconstruction in a variety of types of physical quantum systems, ranging from electron spin to atomic position ^[1–4]. Considerable effort has been put toward the development of these experimental techniques and improving their scalability. These techniques have found applications even outside their original purposes, such as in classical image processing ^[5].

Early schemes for quantum information processing with photons proposed encoding information in a set of $d$ paths that a photon could be in ^[6,7]. 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 a 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” ^[8]. Higher-dimensional encodings can also simplify the design of quantum logic circuits ^[9]. 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 ^[10]. A seminal contribution to manipulating path states was a scheme to implement any chosen discrete $d \times d$ unitary using an interferometric array of beam splitters and phase shifters acting on $d$ paths ^[7,11]. On-chip implementations of the above-mentioned universal discrete unitaries are an especially active research direction ^[12,13] due to 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 $O (d^{2})$ switchable elements such as beam splitters and mirrors. As far as we know, neither strategy for $d$ -path quantum state reconstruction has been implemented above $d = 2$ .

Sign up for Photonics Research TOC. Get the latest issue of Photonics Research delivered right to you！Sign up now

On the other hand, schemes have been devised and implemented to reconstruct the full spatial quantum state of a photon ^[14–17]. 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 are ill suited to reconstructing path-encoded quantum states. Either they do not provide enough information (e.g., the gradient rather than the value of the phase ^[16]), or they provide too much information (e.g., the spatial distribution of each path). For a fixed measurement time, one would expect these latter 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 determines only the parameters of the $d$ -path density matrix $ρ$ 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 much more measurement outcomes than $d$ . Hence, projections in many incompatible bases can be performed with one measurement apparatus configuration. This is accomplished by measuring either position $x$ , with a camera, or momentum $k$ , by adding a cylindrical lens that performs a one-dimensional optical Fourier transform (OFT). A conceptual diagram of the technique is shown in Fig. 1. When the paths are distributed in both $x$ and $y$ transverse dimensions, the lens rotates to a finite set of angles in the $x y$ 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 ^[18,19].

Figure 1.The working principle of path-encoded quantum state reconstruction. (a) Geometry: the spatial arrangement of paths. (b), (c) Paths passing through a cylindrical lens to an image sensor. Along one direction, the paths are interfered by the lens. Along the other direction, the paths are unaltered. The off-diagonal elements of the density matrix $ρ$ are found by performing a discrete Fourier transform of the recorded interference pattern. The cylindrical lens allows for only chosen sets of paths to be interfered at a time. This allows the method to accommodate duplicate path spacings in the geometry.

2. THEORY

Before introducing our method, we begin by rigorously defining what is meant by path-state encoding. A single photon is traveling along $z$ and has transverse position $x$ (for clarity, we restrict ourselves to one transverse dimension for now). A set of $d$ states ${{| ψ_{i} ⟩}_{x}}$ is defined in this $x$ space. Each has a narrow position distribution ${| ψ_{i} (x) |}^{2} = | ⟨ x {| ψ_{i} ⟩}_{x} |^{2}$ and is displaced in $x$ from the other states. The displacement is much larger than the width of any state distribution, which ensures that they have minimal overlap, $⟨ ψ_{i} {| ψ_{j} ⟩}_{x} \approx δ_{i j}$ (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 ${| ψ_{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 affects only the phase of the distribution. Hence, this identicality ensures that $| ⟨ k_{x} {| ψ_{i} ⟩}_{x} |^{2} = {| \tilde{ψ} (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 $ρ = \sum_{i j}^{d} ρ_{i j} | ψ_{i} ⟩ ⟨ ψ_{j} |_{x}$ .

We now describe the basic concept behind our reconstruction method. Each diagonal element of the density $ρ_{i i}$ is the probability for the photon to be in the corresponding path $i$ . These diagonals are straightforward 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 $ρ_{i j} = | ρ_{i j} | \exp (ϕ_{i j})$ 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 $| ρ_{i j} |$ . The phase difference between a photon passing through slits one and two, equivalent to $ϕ_{i j}$ , sets the transverse offset of the interference pattern. In summary, if one can measure the visibility and phase of the interference between paths $i$ and $j$ , one can determine $ρ_{i j}$ 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 λ k_{x} / 2 π$ , where $λ$ is the wavelength of the photon. In this way, a lens performs an OFT. A general state $ρ (x, x^{'}) = ⟨ x | ρ | x^{'} ⟩$ in the continuous position basis corresponds to $ρ (k_{x}, k_{x}^{'}) = ⟨ k_{x} | ρ | k_{x}^{'} ⟩$ 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}) = ρ (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 $ρ_{i j}$ from the measured momentum distribution. Our two paths, ${| ψ_{1} ⟩}_{x}$ and ${| ψ_{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}) = {| \tilde{ψ} (k_{x}) |}^{2} \sum_{i, j}^{d} ρ_{i j} \exp (i L_{i j}^{x} k_{x})$ , where $L_{i j}^{x} = | x_{i} - x_{j} |$ is the pairwise distance between the paths $i$ and $j$ in the $x$ direction. Using the Hermiticity of $ρ$ , this can be simplified to where $ϕ_{12}$ is the phase of $ρ_{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 $Tr {ρ} = 1$ , the visibility can be found to be $2 | ρ_{i j} | / (ρ_{11} + ρ_{22}) = 2 | ρ_{i j} |$ . As we noted earlier, the visibility is directly proportional to the magnitude of $ρ_{i j}$ , and the phase $ϕ_{i j}$ of the cosine is the phase of $ρ_{i j}$ .

We can extract the magnitude and phase information from the interference pattern in Eq. (1) by using a DFT, 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 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 $ρ_{12}$ . Namely, the complex value of the Fourier transform at a position equal to spacing $L_{12}^{x}$ is $ρ_{12}$ , i.e., $ρ_{12} = F_{k_{x}} {P (k_{x})} (\bar{x} = L_{12}^{x})$ . In short, a Fourier transform of the measured interference pattern directly gives the phase and magnitude of the off-diagonal such as $ρ_{i j}$ .

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.

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 ${| ψ_{i} ⟩}_{x y}$ is centered at position $P_{i} = (x_{i}, y_{i})$ . We call the set of all positions $P_{i}$ the “geometry.” For simplicity, we assume that each state has a position wavefunction that is a two-dimensional Gaussian with equal width $σ$ in both $x$ and $y$ directions. This ensures that the distribution is symmetric under rotations in the $x y$ plane and, also, that it is a product state: ${| ψ_{i} ⟩}_{x y} = {| ψ_{i} ⟩}_{x} {| ψ_{i} ⟩}_{y}$ . The subscripts $x$ and $y$ distinguish the two Hilbert spaces. The position wavefunctions are, thus, $ψ_{i} (x, y) = ψ (x - x_{i}) ψ (y - y_{i})$ . As before, the paths must be separated by much more than $σ$ so that they have minimal overlap, $⟨ ψ_{i} {| ψ_{j} ⟩}_{x} \approx δ_{i j}$ and $⟨ ψ_{i} {| ψ_{j} ⟩}_{y} \approx δ_{i j}$ . 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, $ρ = \sum_{i j}^{d} ρ_{i j} {| ψ_{i} ⟩}_{x y} ⟨ ψ_{j} |$ .

We must now carefully label and then distinguish every possible pair of paths, since each is associated with coherence $ρ_{i j}$ . The positions of two paths, $i$ and $j$ , are connected by the line segment $L_{i j}$ with end points $(P_{i}, P_{j})$ . The spacing of the path pair is given by the length of the line segment, $L_{i j} = \sqrt{{(L_{i j}^{x})}^{2} + {(L_{i j}^{y})}^{2}}$ , where $L_{i j}^{s} = | 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 coherences $ρ_{i j}$ 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 $L_{i j}$ 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 coherences $ρ_{i j}$ 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 OFT in both $x$ and $y$ directions. The required path geometry for this, called a Golomb rectangle, is discussed in Appendix A. However, consider whether in our geometry, there are two parallel line segments $L_{i j}$ and $L_{i^{'} j^{'}}$ with the same length, $L_{i j} = L_{i^{'} j^{'}}$ . In this case, the contribution to the DFT from $ρ_{i j}$ and $ρ_{i^{'} j^{'}}$ cannot be distinguished, since they will appear at the same peak position $(\bar{x} = L_{i j}^{x}, \bar{y} = L_{i j}^{y})$ .

A cylindrical lens solves this problem as long as $L_{i j}$ and $L_{i^{'} j^{'}}$ are not both segments of the same line. As an example, in Fig. 2(a), $L_{16}$ , $L_{25}$ , and $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 a cylindrical lens with its OFT axis parallel to these example line segments $L_{i j}$ . The resulting probability distribution will be of momentum in the direction of $L_{i j}$ and position in the perpendicular direction [Fig. 2(b)]. Critically, the line segments, $L_{16}$ , $L_{25}$ , and $L_{34}$ , have distinct positions along the axis perpendicular to $L_{i j}$ . 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 $ρ_{16}$ , $ρ_{25}$ , and $ρ_{34}$ by taking a one-dimensional DFT of the recorded camera image along the cylindrical lens axis. This concept can be extended to other coherences of $ρ$ by rotating the cylindrical lens and, hence, OFT axis to all the unique angles $θ_{i j}$ in our full set of line segments, where $θ_{i j} = \arctan [(y_{i} - y_{j}) / (x_{i} - x_{j})]$ . In Appendix A, we discuss the experimental resource requirements of the present tomographic method, as well explicit algorithms for constructing compatible geometries of path states. An example of the recorded image for one such angle is shown in Figs. 2(c) and 2(d). As long as no two $L_{i j}$ and $L_{i^{'} j^{'}}$ are segments of the same line, this procedure can be used to determine the density matrix.

Figure 2.State reconstruction method. (a) Six paths are shown in the figure and encode the state $ρ$ . The optical Fourier transform (OFT) axis (blue solid line) rotates to particular angles $θ_{i j}$ , of which a few are shown, to interfere each pair of paths at a time (angles are with respect to a horizontal axis along the bottom-most paths). We assign to each pair of points a line segment $L_{i j}$ . (b) Corresponding OFT for the eight angles required to reconstruct this particular density matrix. As only paths with angle $θ_{i j}$ between them interfere, the diagonal elements can be recovered from the remaining paths. The $k_{y}$ axis is always in the direction of the interference, and the $x$ axis is perpendicular to it (example shown for $θ_{15}$ ). (c) Each pattern is recorded and analyzed one at a time via discrete Fourier transform (DFT) by taking a one pixel wide “slice” through the interference pattern. This process is repeated for every interference pattern present in a given image. (d) Fourier transform of the interference pattern (for illustrative purposes, we plot the magnitude). The magnitude $ρ_{i j}$ is recovered from the height of the DFT at the position $\bar{y} = L_{i j}$ . The normalization is obtained by summing the zero frequency peaks of each interference pattern present in the $θ_{65}$ subpanel (in this example) in panel (b). All panels contain real data.

For clarity in our mathematical description of this procedure, instead of rotating the OFT, we rotate the full quantum state $ρ (x, x^{'})$ in position space. Because our path wavefunctions $ψ_{i} (x, y)$ are rotationally invariant, this task is reduced to find the rotated path positions $P_{i}^{'}$ . If $R (θ)$ is the standard rotation matrix, then $P_{i}^{'} = R (θ) P_{i}$ . After this rotation, the wavefunctions are $ψ_{i}^{'} (x, y) = ψ (x - x_{i}^{'}) ψ (y - y_{i}^{'})$ . In the following discussion, we drop the primed notation for clarity, assuming that the paths have been rotated to the appropriate angle $θ$ . 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. 2(b) for $θ_{34} = 90 °$ . 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 $ρ_{x} (k_{y})$ , which is still an operator in $x$ space, as indicated by the subscript $x$ : Using the Fourier shift theorem, $⟨ k_{y} {| ψ_{i} ⟩}_{y} = F {ψ (y - y_{i})} = \tilde{ψ} (k_{y}) \exp (i y_{i} k_{y})$ , we simplify this to the following expression: We will drop the envelope ${| \tilde{ψ} (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 Eq. (4), first we separate the terms where $i = j$ : 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 $P_{i}$ of the paths in the rotated geometry. If ${| ψ_{i} ⟩}_{x}$ is not centered on $x_{m}$ , then the wavefunction is nearly zero at $x = x_{m}$ . In other words, $⟨ x_{m} {| ψ_{i} ⟩}_{x} = ψ (x_{m} - x_{i}) = ψ (0) δ_{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 For notational convenience, we drop the ${| ψ (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 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_{i j}^{y}$ , then the peak at $F_{k_{y}} {P (x_{m}, k_{y})} (\bar{y} = L_{i j}^{y})$ has value $ρ_{i j}$ . 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) = ρ_{i i}$ , 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 $θ_{i j}$ , the full $d$ -path density matrix can be determined. This method works as long as no two parallel line segments $L_{i j}$ 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 $ρ$ that we will later demonstrate experimentally.

3. 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 $λ = 808 nm$ , we use a series of wave plates and calcite beam displacers, with the crystal axis cut at 45°, to generate $d = 6$ paths that contain state $ρ$ . 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 μm. The first beam displacer produces a transverse shift of $δ_{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 $δ_{y} = 2.7 mm$ as well. Last, the third crystal shifts horizontally polarized light by $δ_{x} = 4 mm$ horizontally. In total, the beam displacers produce eight paths. However, this geometry is incompatible with the tomography method. For example, $L_{13}$ and $L_{57}$ are both segments of the same line and of equal length. Therefore, we block two of the eight paths to obtain a $2 \times 3$ path geometry that meets all the conditions necessary for the reconstruction method to work. The resulting six-dimensional density matrix is a function of the half-wave plate (HWP) and quarter-wave plate (QWP) angles $ϕ$ , $ζ$ , and $Ω$ , preceding each crystal, shown in Fig. 3(a). As such, we write that $ρ = ρ (ϕ, ζ, Ω)$ . 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, $L_{12}$ , $L_{34}$ , and $L_{56}$ (among others), as shown in Fig. 3(a).

Figure 3.Experiment demonstrating the state reconstruction method. (a) State preparation in blue box: the Rayleigh length of an 808 nm diode laser is set by a beam expander. A series of displacement crystals (xtal) and half- and quarter-wave plates (labeled by the angles $ϕ$ , $ζ$ , and $Ω$ ) generate the state $ρ$ . The resulting eight-path geometry is not compatible with the tomography method, and so two paths are blocked to produce a compatible six-dimensional state. A set of HWP and QWP may be inserted to form a mixed state by rapidly spinning ${HWP}_{s}$ . The purity is a function of the wave-plate angle $τ$ . We can also produce photon pairs via SPDC at 808 nm using a diode at $λ = 404 nm$ to pump a 15 mm ppKTP crystal. The measured $g^{(2)} (0)$ of the source is $0.1979 \pm 0.0005$ . (b) Analysis is presented in the purple box: lenses $f_{1} = 1000 mm$ and $f_{2} = 400 mm$ image the six paths onto a camera (an electron multiplying CCD in the case of down-converted photons). A rotating cylindrical lens ( $f = 250 mm$ ) performs the optical Fourier transform (OFT) along the OFT axis. A one pixel wide slice of each interference pattern is analyzed with DFT on a computer. Wider slices can be used and averaged over; however, this may reduce the visibility if there are imperfections in the interference pattern. This would include tilting of the dark fringes, or, as can be seen in the figure, if the intensity in each bright fringe is not evenly distributed. The coherences are obtained by the heights of the Fourier transform peaks, normalized by the total intensity. No filters are applied to the raw data.

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 $γ = 2.40 μm$ . The lens rotates in an automated mount to each of the right angles $θ_{i j}$ . The one-dimensional OFT produces sets of fringes, one for each pair of paths for which the angle between them is $θ_{i j}$ . To obtain the coherence, first we rotate the image back by $- θ_{i j}$ 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 DFT 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_{i j}$ of a peak in the DFT at spacing $\bar{y} = L_{i j}^{y}$ is proportional to $ρ_{i j}$ . We define a scaling factor so that $| ρ_{i j} | = | S_{i j} / S |$ , where $S = \sum_{i}^{d} S_{i i}$ is the total recorded intensity. We can obtain $S$ directly from the images used to calculate $ρ_{i j}$ , without having to change the apparatus. Specifically, we note that in Eq. (7), the following holds: $\sum_{x_{m}} \sum_{i}^{d} ρ_{i i} δ_{x_{m} x_{i}} = 1$ . In terms of the recorded signal, we get $\sum_{x_{m}} \sum_{i}^{d} S_{i i} δ_{x_{m} x_{i}} = S$ , i.e., 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 are required to obtain this scaling.

The complex peaks also determine the phase of coherence $ρ_{i j}$ 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 $ρ_{i j}$ is proportional to in the DFT, $\bar{y} = L_{i j}^{y}$ or $\bar{y} = - L_{i j}^{y}$ . Any two states whose coherences are different and non-zero may be used as a reference.

4. RESULTS

The reconstruction method is demonstrated with $d = 6$ dimensional states of the form $ρ (ϕ = 22.5 °, ζ, Ω = 22.5 °) = ρ (ζ)$ (i.e., all HWPs are held fixed), whose coherences are a function of the QWP angle $ζ$ . Figure 4(a) shows the real and imaginary parts of the 15 off-diagonal elements $ρ_{i j}$ (not including the conjugates, $ρ_{i j}$ ). Note that several elements of $ρ_{i j}$ 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 theoretical values (dotted lines). Figure 4(b) shows the experimental and theoretical density matrix for $ζ = 30 °$ in a side-by-side comparison.

Figure 4.Experimental results. (a) Experimental (dots) and theoretical (curves) coherences $ρ_{i j}$ of the density matrix $ρ$ . These are produced by varying the QWP angle $ζ$ in Fig. 3. As the coherences are constrained by the experimental setup, only a few unique values appear in any given matrix. As such, data points for multiple coherences overlap. Note that error bars, obtained by averaging over multiple pictures, were omitted for clarity but range from $10^{- 3}$ to $10^{- 2}$ . (b) Experimentally reconstructed six-dimensional state $ρ (ζ = 30 °)$ . Each diagram represents a 6×6 matrix, with theoretical elements to the right of each experimental element. The fidelity with the nominal input state is 0.9911 (fidelity is one if the states are identical).

Next, we calculate the state fidelity, $F (ρ, σ) = Tr {\sqrt{\sqrt{ρ} σ \sqrt{ρ}}}$ , which is a measure of how close two density matrices $ρ$ and $σ$ are. The fidelity is bound between zero and unity, and $F (ρ, σ) = 1$ only if $ρ = σ$ . We calculate the fidelity between the experimentally reconstructed state $ρ (ϕ, ζ, Ω)$ and the theoretically predicted one $ρ_{th}$ . This is done for sets of states where one of the wave-plate angles $ϕ$ , $Ω$ , or $ζ$ is varied. Specifically, if an HWP is not being varied, then it is fixed at 22.5°. If either HWP $ϕ$ or $Ω$ is varied, then we fix $ζ = 45 °$ . Figure 5(a) displays the results. The average fidelity when $ϕ$ is varied is $0.987 \pm 0.001$ , $0.9893 \pm 0.0008$ when $Ω$ is varied, and $0.979 \pm 0.001$ when $ζ$ is varied. Averaged 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.

Figure 5.Experimental results. (a) Fidelity as a function of the wave-plate angles $ϕ$ , $ζ$ , and $Ω$ , shown in Fig. 3. The fidelity is close to unity, meaning $ρ$ and $ρ_{th}$ are nearly equal. Averaged over all points, the fidelity is $0.9852 \pm 0.0008$ (dashed line). (b) Purity Tr ${ρ^{2}}$ as a function of the HWP angle $τ$ for a classical source and single-photon source. The single-photon source deviates due to the much shorter coherence length. (c) Reconstructed density matrix of single photons in four paths. The experimental values are labeled. The corresponding theoretical values are either 0.25 or 0. The calculated fidelity is $1.00 \pm 0.03$ .

The density matrix generalizes the concept of the state by accommodating statistical mixtures of states. A measure of this is the state purity Tr ${ρ^{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 $ρ$ and calculate the purity. To mix the state, we introduce a series of HWP-QWP- ${HWP}_{s}$ -QWP, shown in Fig. 4(a). Here, the QWPs are 90° apart. The subscript $s$ on ${HWP}_{s}$ indicates that it is rapidly spinning in a mount. ${HWP}_{s}$ 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 wave-plate angle $τ = 0 °$ produces a pure state, and $τ = 22.5 °$ produces a totally mixed polarization state. The theoretical purity of the resulting path state $ρ (τ)$ is $T r {ρ^{2} (τ)} = 1 / 9 [5 + 4 \cos^{2} (4 τ)]$ . Figure 5(b) 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 $ζ$ with an HWP, and set all HWP at 22.5°. The purity of this matrix also depends on $τ$ . 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 $γ = 16 μm$ . The calculated purity is shown in Fig. 5(b). 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. 5(c), we plot the density matrix elements for $τ = 22.5 °$ and list the experimental values. For reference, every theoretical element 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 self-adjoint conditions required of a physically possible density matrix. However, it will not necessarily satisfy the positive semi-definite condition. In order to estimate the implications of this, we alter the measured density matrix to be positive semi-definite by setting every negative eigenvalue to zero and renormalizing. Doing so leads to fidelity of 0.9967, which is within the uncertainty of our raw fidelity. This suggests that the fidelity would not be significantly changed even if more sophisticated methods were used to enforce this constraint. These methods, such as maximum likelihood, could be readily used with the raw camera data from our experimental method but can also introduce unwanted bias ^[20–22] and would obscure the physics behind the method. Consequently, we do not describe them further.

5. SUMMARY

The above results show the reconstruction method faithfully reconstructs several different metrics for a six-dimensional density matrix. As reconstruction of spatially encoded states has previously proven to be difficult to do, it has scarcely been done in the past. For example, the work demonstrated in Ref. ^[18] provides the closest analogue. However, the results shown are lower in dimension, and it is not clear how the technique would be extended to higher dimensions, making a true comparison difficult to do. The technique shown in this work presents a feasible and intuitive method to reconstruct high-dimensional states. 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) 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.

Appendix A: 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 $η$ , 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 $η$ of angles $θ_{i j}$ required to determine all the density matrix off-diagonal elements $ρ_{i j}$ depends on the specific geometry. In the worse case, each lens angle retrieves only one element. Thus the maximum number of angles would be $d (d 1) / 2$ , the number of independent off-diagonal complex elements in $ρ$ . More typically, a number of paths will lie along the same line $L$ . In this case, multiple off-diagonals can be found with a single lens angle aligned along $L$ . Consequently, judicious choice of geometry can reduce the total number of measurement settings significantly.

Consider a one-dimensional geometry, i.e., the path states are arranged along, say, the $y$ axis. In this case, only one experimental setting is needed to determine 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 Erd s–Turan construction for an asymptotically optimal set is $x_{i + 1} = L_{\min} [2 d i + (i^{2} \mod d)]$ , for $i = 0, 2, \dots, d 1$ for any value of $d$ that is an odd prime [23]. Here, $L_{\min}$ is the minimum possible spacing between paths. The maximum spacing, $L_{\max} = x_{d} x_{1}$ , is then given by $L_{\max} = L_{\min} [2 d (d 1) + 1]$ . Space is a resource for both 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_{\max} = L_{\max} / L_{\min}$ . In contrast, a Golomb ruler uses space less efficiently: $d = \sqrt{L_{\max} / (2 L_{\min})} = \sqrt{d_{\max} / 2}$ for large $d$ . A one-dimensional geometry is optimal for the number of measurement settings required to reconstruct $ρ$ (i.e., $η = 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 [24]. In the area of astronomical interferometers, these are known as “non-redundant configurations” [25]. In it, every line segment $L_{i j}$ is associated with a unique vector, i.e., orientation and length. Given a non-redundant configuration, a spherical lens is sufficient to reconstruct $ρ$ . As with the one-dimensional geometry, two measurement settings are sufficient: lens present and lens absent. Given $L_{\max}$ and $L_{\min}$ , explicit computational searches for optimal Golomb rectangles suggest that maximum number of path states is $d = O (L_{\max} / L_{\min})$ [25]. In two dimensions, the maximum number of closely packed path states is $d_{\max} = O [{(L_{\max} / L_{\min})}^{2}]$ . Consequently, the optimum Golomb rectangles achieve $d = O (\sqrt{d_{\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 overall constant factor.

A number of experimental factors will constrain the largest feasible $L_{\max}$ . The most obvious is the range of positions $D$ over which the 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 $γ$ , the distance between pixel centers. The finest interference pattern produced by the paths will have a fringe spacing of $2 π L_{\max} / λ 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_{\max} < λ f / π γ$ . On the other hand, $L_{\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 $γ = λ f / π L_{\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_{\min}$ . This sets $W = λ f / π L_{\min}$ . Together, these requirements determine the number of pixels along one camera dimension to satisfy $N > W / γ = L_{\max} / L_{\min} = d_{\max}$ . As discussed above, $L_{\max} < D$ . In addition, $L_{\min}$ must be greater than the path-state position width, $L_{\min} > σ$ . Taking them together, we have the requirement $N > D / σ$ . Or, using the results from the optimal Golomb ruler above, $N > 2 d^{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 requires only 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_{\max}$ will improve. The tradeoff is that multiple lens angles and, thus, measurement settings $η$ , 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 tradeoff in more detail.

References

[1] M. G. Paris, M. F. Sacchi. Quantum tomagraphy. Advanced Imaging and Electron Physics, 128, 206-309(2003).

[2] J. Lundeen, A. Feito, H. Coldenstrodt-Ronge, K. Pregnell, C. Silberhorn, T. Ralph, J. Eisert, M. Plenio, I. Walmsley. Tomography of quantum detectors. Nat. Phys., 5, 27-30(2009).

[3] D. F. James, P. G. Kwiat, W. J. Munro, A. G. White. On the measurement of qubits. Asymptotic Theory of Quantum Statistical Inference: Selected Papers, 509-538(2005).

[4] L. Rippe, B. Julsgaard, A. Walther, Y. Ying, S. Kröll. Experimental quantum-state tomography of a solid-state qubit. Phys. Rev. A, 77, 022307(2008).

[5] G. M. D’Ariano, C. Macchiavello, M. G. Paris. A fictitious photons method for tomographic imaging. Opt. Commun., 129, 6-12(1996).

[6] G. J. Milburn. Quantum optical Fredkin gate. Phys. Rev. Lett., 62, 2124-2127(1989).

[7] M. Reck, A. Zeilinger, H. J. Bernstein, P. Bertani. Experimental realization of any discrete unitary operator. Phys. Rev. Lett., 73, 58-61(1994).

[8] H. Bechmann-Pasquinucci, W. Tittel. Quantum cryptography using larger alphabets. Phys. Rev. A, 61, 062308(2000).

[9] B. P. Lanyon, M. Barbieri, M. P. Almeida, T. Jennewein, T. C. Ralph, K. J. Resch, G. J. Pryde, J. L. O’brien, A. Gilchrist, A. G. White. Simplifying quantum logic using higher-dimensional Hilbert spaces. Nat. Phys., 5, 134-140(2009).

[10] J. C. Matthews, A. Politi, A. Stefanov, J. L. O’brien. Manipulation of multiphoton entanglement in waveguide quantum circuits. Nat. Photonics, 3, 346-350(2009).

[11] W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer, I. A. Walmsley. Optimal design for universal multiport interferometers. Optica, 3, 1460-1465(2016).

[12] B. J. Metcalf, N. Thomas-Peter, J. B. Spring, D. Kundys, M. A. Broome, P. C. Humphreys, X.-M. Jin, M. Barbieri, W. S. Kolthammer, J. C. Gates, B. J. Smith, N. K. Langford, P. G. R. Smith, I. A. Walmsley. Multiphoton quantum interference in a multiport integrated photonic device. Nat. Commun., 4, 1356(2013).

[13] A. Peruzzo, A. Laing, A. Politi, T. Rudolph, J. L. O’brien. Multimode quantum interference of photons in multiport integrated devices. Nat. Commun., 2, 224(2011).

[14] C. Bamber, J. S. Lundeen. Observing Dirac’s classical phase space analog to the quantum state. Phys. Rev. Lett., 112, 070405(2014).

[15] D. McAlister, M. Beck, L. Clarke, A. Mayer, M. Raymer. Optical phase retrieval by phase-space tomography and fractional-order Fourier transforms. Opt. Lett., 20, 1181-1183(1995).