Quantum state tomography of orbital angular momentum photonics qubits via a projection-based technique

While measuring the orbital angular momentum state of bright light beams can be performed using imaging techniques, a full characterization at the single-photon level is challenging. For applications to quantum optics and quantum information science, such characterization is an essential capability. Here, we present a setup to perform the quantum state tomography of photonic qubits encoded in this degree of freedom. The method is based on a projective technique using spatial mode projection via fork holograms and single-mode fibers inserted into an interferometer. The alignment and calibration of the device is detailed as well as the measurement sequence to reconstruct the associated density matrix. Possible extensions to higher-dimensional spaces are discussed.


I. INTRODUCTION
In the past decades, light beams carrying orbital angular momentum (OAM) have raised a considerable interest [1,2]. Laguerre-Gaussian beams for instance have been used in a variety of fields for applied and fundamental investigations, including hyperdense data transfer [3], particule-motion control [4,5], enhanced-sensitivity measurements [6,7], fundamental tests of quantum mechanics [8,9] or quantum information protocols [10][11][12][13][14]. Quantum information science and technology has indeed recognized single photons in OAM state superpositions as promising information carriers as the high dimensionality of the Hilbert space they live in could enable enhanced security and multiplexing [15]. In this framework, optical and quantum memories for OAM states have recently seen tremendous developments [16][17][18][19][20][21], opening the path to quantum networks and scalable communication architectures based on this degree of freedom.
In these quantum information related applications, it is crucial to measure the orbital angular momentum state of single photons, i.e. to perform the full quantum state tomography. While OAM measurement for bright beams of light can be done via imaging, it is much more challenging in the single-photon regime. The various solutions that have been proposed and tested heretofore can be classified into two categories. The earlier methods perform single-photon detection after the signal has passed through a projector on an OAM eigenstate. They typically have low efficiencies but good selectivity. The more recently-developed mode-sorting techniques consist in changing the signal path depending on its OAM value. These methods usually lead to higher efficiencies but their implementation is also often more challenging.
In this paper, we describe an interferometer-based detection setup relying on the projective method, similar to the one that Miyamoto and coworkers used in order to overcome the limitations of the "hologram shifting method" [22], but with an improved efficiency and versatility. Depending on the interferometer phase, our device allows to perform state projection into multiple bases and consequently a full state tomography. This technique has been used in our recent demonstration of quantum memory for OAM encoded qubits [21]. We provide here a detailed and quantitative study of its implementation, and investigate extensions to higher-dimensional spaces.
The paper is organized as follows. Section II first briefly reminds the fundamentals on quantum bits and on orbital angular momentum in order to understand our implementation and the required measurements. The full description of our tomography setup is then given in section III. In section IV, a detailed study of the calibration procedure is presented and benchmarks are given. In section V, we finally show examples of qubit state reconstruction. Perspectives are discussed in section VI where we propose two possible extensions of the setup to perform on higher-dimensional OAM Hilbert spaces.

II. QUBITS ENCODED IN THE ORBITAL ANGULAR MOMENTUM DEGREE OF FREEDOM
In this section, we first summarize the representation of qubit states and the required measurements to reconstruct the associated density matrix. We then discuss specifically the implementation of OAM qubits with superpositions of Laguerre-Gaussian modes.

A. Qubit representation and quantum state tomography
A qubit, i.e. a two-dimensional quantum system, evolves in a Hilbert space spanned by two basis vectors usually denoted |0 and |1 in analogy with classical information. For pure state, it can be represented by: with |α| 2 + |β| 2 = 1.
The tomography of such superposition states requires measurements performed in three different bases [23]. In addition to the logical basis {|0 ,|1 }, two additional mutually unbiased bases can be defined as the superpositions: {|± = |0 ± |1 } and {| ± i = |0 ± i|1 }. Measuring the qubit in the logical basis will yield the values of |α| 2 and |β| 2 and measurements in the two superposition bases will provide the relative phase between the coefficients α and β. This sequence of measurement actually allows to perform the reconstruction of the density matrix of any mixed state, which in the general case can be expressed as : where theσ i are the Pauli matrices and the S i = Tr(ρσ i ) coefficients are usually called Stokes parameters when dealing with polarization states. The S i coefficients indicate the relative weights of either basis state in the different bases. Indeed, they can written as where p |Ψ is the probability to measure the qubit in the state |Ψ . If the qubit is in a pure state as expressed in equation 1, then and the Stokes parameters can be easily related to the parameters α and β: In this work, we focus on qubits encoded in the orbital angular momentum degree of freedom. We now briefly review the essential features of the orbital angular momentum of light, focusing on the Laguerre-Gaussian modes that offer a convenient basis to describe it.

B. Properties of the Laguerre-Gaussian modes
The set of Laguerre-Gaussian modes is a complete orthonormal basis of solutions to the paraxial propagation equation. Their transverse shape is propagation invariant, as can be seen from the following expression of their electric field amplitude: where E 0 is the electric field amplitude, K lp = 2 π p! (l+p)! is a normalization constant, and L α n (x) are the generalized Laguerre polynomials. Here, we have chosen a wave traveling along the z axis and polar coordinates (r, θ) to parametrize the transverse plane. The parameters w(z) = w 0 1 + ( z /z R ) 2 , z R = π w 2 0/λ, R(z) = z 1 + ( z R /z) 2 , and ζ(z) = arctg( z /z R ) are the radius, Rayleigh length, radius of curvature and Gouy phase for a beam of waist w 0 at a wavelength λ.
In contrast to the standard TEM 00 = LG l=0 p=0 mode, the higher-order Laguerre-Gaussian modes have a rotating phase profile e ilθ with a singularity at the origin. Due to this rotating phase, the local Poynting vector has a non-vanishing component along the orthoradial direction and the beam thus exhibits an orbital angular momentum around this axis. The orbital angular momentum carried by each photon in such a mode is equal to the index l ∈ Z (in units), which is also equal to the circulation of the phase around the axis divided by 2π. As the LG modes are eigenfunctions of the propagation equation, their transverse shape, and hence the OAM of the photons, is preserved. This makes the OAM number l a relevant quantum number for information encoding. As required for the smoothness of the electric field amplitude, the phase singularity is associated with an intensity nulling in the middle as the 2l-th power of the radial coordinate ( r /w(z)) 2|l| . This feature gives the LG modes their characteristic doughnut-shaped intensity profiles.
The other index, p ∈ N, describes the radial shape of the beam. With each additional p unit, the amplitude has an additional sign change along the radius and the intensity an additional zero-value ring. This number has a less straightforward interpretation than the orbital index l and has therefore been subject to less investigation hitherto [24,25].

C. Superpositions of Laguerre Gaussian modes and OAM qubit encoding
One promise of the OAM degree of freedom for information encoding is the large (potentially infinite) dimension of the available Hilbert space. We first restrict ourselves to a two-dimensional space spanned by the modes LG l=±1 p=0 . Extensions of the tomography scheme to higherdimensional OAM spaces are discussed in section VI.
For a single photon living in this qubit space, we define the corresponding logical basis vectors as |0 = |l = +1 = |R and |1 = |l = −1 = |L . R and L re- spectively refer to the right and left handedness of the helical wavefront. The OAM difference of 2 between the modes defining the logical basis ensures a very good distinguishability as will be explained later. As the index p does not play a significant role in the present study, we drop it and we consider only modes with radial index p = 0.
Equally-weighted superpositions |Ψ = |R + e iφ |L , which span the equatorial plane of the Bloch sphere as shown in Figure 2, correspond to rotated Hermite-Gaussian (HG) modes of the type TEM 01 . These modes consist of two bright spots, separated by a dark line. In the plane orthogonal to the propagation axis, the angle α d of the dark axis with respect to the horizontal axis is related to the relative phase φ by: Accordingly to their spatial shapes (see Fig. 2) and in analogy to the case of a polarization basis, we denote the specific following modes as vertical, diagonal, horizontal and anti-diagonal: The bases {|R , |L }, {|H , |V } and {|D , |A } will constitute the 3 mutually unbiased bases for the tomography.

D. Detection of single photons in Laguerre-Gaussian modes
The determination of the OAM state of a bright beam of light can be done by standard imaging and wavefront measurements (either using interferometry or a microlens array, or any equivalent available technique). This intrinsically requires many photons, so that characterizing an OAM state at the single-photon level calls for other methods. The different techniques developed so far can be classified into two categories: • Projective-based techniques. In these methods, the photons to be measured impinge on a device that performs a projection onto an OAM eigenstate before being measured. The mode projectors are typically made of a hologram that converts an input mode with non-zero l value into a TEM 00 mode followed by a spatial filter (pinhole or single-mode fiber). The holograms can be either fixed [8] or dynamically programmed with a spatial light modulator [9,16], they can be either intensity [26] or phase holograms [8], they can diffract the light to all directions or be optimized for a single output direction. In any case, the projector only selects one mode and photons in other modes are lost.
• Mode-sorting techniques. Here, the propagation direction of the signal is changed depending on its OAM value. This feature overcomes the problem of losses inherent to mode projection. However, these methods are often more challenging than the previous ones. Among them, one can cite a Mach-Zehnder interferometer in the arms of which Dove prisms have been inserted [27]. Another method that has seen significant developments in the few past years relies on a log-polar coordinate interpolation realized with two phase-modulating elements. The radial and polar coordinates r, θ in one plane are mapped one by one onto the cartesian x and y coordinates in a subsequent plane. This approach was first implemented with two spatial light modulators [28] then with fixed refractive optics [29][30][31]. It is also possible to take benefit from the fact that any unitary manipulation of transverse modes (and hence mode sorting) can be achieved by multiple phase modulation steps separated by optical Fourier transforms [32].
As underlined before, quantum state tomography requires projection into different bases in order to access the coherence terms of state superpositions. In this context, building on the projective techniques in order to benefit from their simplicity and high distinction ratio, we developed an OAM tomography setup, as described in the next section.

III. SETUP FOR QUANTUM STATE TOMOGRAPHY OF OAM QUBITS
The state tomography requires to project the state to be characterized into three bases. In this section, we present an interferometric setup based on fork holograms enabling such projections. We explain how the interferometric phase defines the projection basis and we detail how to experimentally access this phase via an imaging technique. where each path includes a mode projector based on a blazed fork computer-generated hologram (phase pattern shown) and a single-mode fiber. Holograms labelled l = +1 and l = −1 are respectively placed in the so-called R and L paths. By OAM subtraction or addition, as defined by the orientation of the centered dislocation, they transform modes with orbital quantum number l = +1 (resp. l = −1) into TEM00 modes, which are coupled into single-mode fibers performing spatial filtering. The fiber paths are then recombined and the two outputs are directed towards output X and Y where single-photon counting modules are located (SPCM-AQR-14-FC). A phase reference beam (green arrows) is injected backwards and detected by a digital camera at the input beam-splitter in order to measure the phase ϕ of the interferometer. The value of ϕ defines the projection basis.

A. Experimental setup: Interferometer and mode projectors
The apparatus is schematized in Figure 3. The incoming state is first split using a non-polarizing beam splitter. Each of the subsequent paths includes a mode projector onto an OAM eigenstate. These mode projectors are based on the combination of a hologram and a singlemode fiber [8]. A blazed fork phase-hologram diffracts the light and performs OAM addition or subtraction depending on its orientation. Thus, on one path, the mode LG l=+1 is converted into a mode LG l=0 = TEM 00 , which is then efficiently coupled to the single-mode fiber, while any other mode is converted into a Laguerre-Gaussian beam with a non-vanishing l value and hence not coupled to the subsequent fiber. There are two such paths, denoted R and L, that are arranged to project the incoming state onto the |R and |L states respectively. The diffraction efficiency of the holograms is 80% and the coupling efficiency to the single-mode fiber is also around 80%, leading to an overall transmission around 65% for the corresponding mode. The rejection of the other mode is optimized using the alignment described later, with typical value around 25 dB.
As shown in Figure 3, the two paths are then recombined via a fiber beam-splitter with two outputs labelled X and Y . The difference in propagation length along the two arms of the interferometer causes a phase shift denoted ϕ. If the input is in a Laguerre-Gaussian mode |R or |L , then light will only be coupled into one of the interferometer arms, and this light will be equally distributed over X and Y regardless of ϕ. In contrast, if the input is in a superposition state, then there is a non-zero amplitude in both arms and these amplitudes will interfere. The probabilities to detect light at either outputs X or Y will vary sinusoidally with ϕ.
More specifically, let us take the example of a pure state |Ψ = a|R + be iφ |L (with a, b ∈ R) entering the device. At the output, it will be transformed into the state: where we have assumed perfect transmission and mode rejection. The events detected at the output X for instance will thus be proportional to They correspond thus to the projection of the incoming state on the state |R + e iϕ |L . By choosing ϕ, any projection basis in the equatorial plane of the Bloch sphere can thus be chosen; this is the key feature of this setup. In summary, the interferometer directs light in a given Hermite-Gaussian state towards output X and the orthogonal mode towards output Y . These two states, i.e. the basis in which this separation occurs, depend on the value of the interferometer phase ϕ. For ϕ = 0, photons in the |H state will be directed towards output X while L path R path Full interferometer blocked blocked ϕ = 0 ϕ = π/2 ϕ = π ϕ = 3π/2 photons in state |V will be directed towards output Y . In the same way, for ϕ = π/2, the incoming photons are measured with respect to the basis |A and |D . Finally, if one of the paths is blocked, then the device will act as a projector onto the OAM eigenstate |R or |L . In this case, the detectors at X and Y receive the same signal. Therefore, when the device is properly calibrated (as explained in section IV), it can be regarded as a black box performing state projection and yielding photon counts in the outputs X and Y as summarized in Table I.

B. Scanning and measurement of the interferometer phase
In order to change the projection basis, one mirror inside the interferometer is mounted on a piezoelectric transducer, allowing to vary ϕ. In the following, we explain how to access this phase using backwardspropagating reference light.

Phase reference beam.
Thermal and mechanical drifts continuously change the interferometer phase ϕ on a scale of a few degrees in a few seconds. To access this phase, a phasereference beam is injected backwards into the interferometer. While it propagates backwards, the reference beam crosses the two holograms. The TEM 00 modes emerging out of the fibers are thus converted into an LG l=−1 mode in the R path, and into an LG l=+1 mode in the L path. These modes are superimposed at the input beam-splitter with a phase difference that is precisely equal to the interferometer phase ϕ. As a result, we get the superposition LG l=+1 + e iϕ LG l=−1 . Similarly to equations 4 and 5, the equal weight superposition of LG l=+1 and LG l=−1 with a relative phase ϕ results in a rotated Hermite-Gaussian TEM 01 mode, consisting of two bright spots with a dark line between them. The dark line makes an angle with the horizontal axis. Measuring this angle enables us to access the interferometer phase ϕ. We now detail how to analyze the images taken by the digital camera to efficiently extract this information.
In a first step, all images are enhanced by applying a median filter (to reduce high-frequency noise and dead pixels) and a midtone stretching filter (to increase the contrast in the middle intensity region and to reduce variations in the high-and low-intensity regions).
When ϕ varies, the dark line in the intensity pattern of the reference beam rotates around the beam axis, according to equation 8. First, the center around which the dark line rotates has to be determined. For this purpose, many images for different values of ϕ are required, covering roughly uniformly the whole range of ϕ ∈ [0, 2π[ (corresponding to the range of [0, π[ for the angle α d of the dark line). Averaging them results in a ring-shaped image, as shown in Table II, to which a doughnut-like distribution is fitted. The fit provides two parameters: the position of the center and the radius of interest. As long as the alignment of the reference path is not changed, these values remain valid for all images. Consequently, this initial procedure has to be performed only once for a measurement series: either during the calibration, if a real-time analysis of the phase is desired, or using a part of the stored images, if post-processing is performed.
The following analysis of the individual images is illustrated in Fig. 4. The center is used as the origin of polar coordinates, while the radius of interest defines the area that will be analyzed.
• This circular area is first divided into N angular bins ("pie slices"), where N has to be divisible by 8 so that there will be angle bins corresponding to each of the four |H , |D , |V and |A modes [42].
In our experiments, a typical value of N was 120.
Image Contrast without enhancements Fitted center enhancement & projections TABLE II: Illustration of the steps involved in determining the center of symmetry of the images. The image displayed in the first column is the average of a few hundred phasereference pictures. After treating it with the median filtering and midtone stretching operations explained in the text, the second image is obtained. Their difference shows the importance of the image processing; without it, the background would be strong enough to shift the center out of its actual position. In the second image, the projections of the intensity onto x and y axes are shown. They are used to obtain the starting values for the fit: the first-order momenta give an estimated center position (red dot), the second-order momenta a starting value for the width of the ring. The last column shows the fit output: the center is depicted in red and the circle with the radius of interest in green. The intensity values are processed as explained in the text to obtain the folded and averaged intensity valuesĨ. These are plotted as blue dots. In this example, the smallest value ofĨ was found in the 48 • angle bin, which is marked by a blue box and corresponds to the angle α d of the dark axis. In panel (d), the pixels belonging to this angular bin are colored in light blue to mark the dark axis, and the area of interest is indicated by the green circle. The original image size was 330 × 330 pixels.
• For each angular bin k, the average intensity I(α k ) is calculated. The first half of this data (0 − 180 • ) is plotted in Fig. 4c as open red circles. We could now fit a sinusoidal function to I(α k ) to determine α d as its phase. However, since we have to process many images and since fits are computationally expensive, we choose the following straight-forward calculation instead.
• Since we are only interested in an axis and not a direction, the angle bins are folded: The intensities of each two opposing bins, i.e. at 180 • from each other, is added. This leaves us with only half the number of bins.
• The dark line should be along the axis of least intensity, but also orthogonal to the axis with most intensity. To account for both conditions at once, we subtract from the intensity of each bin the intensity of the (unique) bin at 90 • from it. These last two steps (folding and subtracting) also rectify slight image asymmetries that are always present (see explanations below and Fig. 5).
• To reduce the influence of remaining bin-to-bin noise, the data is smoothed: the processed inten-sityĨ(α k ) of bin k is calculated as the average over a 45 • wide sector centered around that bin. The last bin N/2 − 1 is considered being next to bin 0, i.e. the angle is taken modulo 180 • . The dataĨ is shown in Fig. 4c as blue dots.
• The angle of the bin with the smallestĨ is the angle α d of the dark axis, providing the interferometer phase ϕ via Eq. 8.
Since the alignment of the beams from the L and R paths is subjected to experimental imperfections, the image analysis routine has been tested against computergenerated images presenting various simulated defects. For this, we numerically generated superpositions of Laguerre-Gaussian beams with either positional or angular misalignment. Examples of such test images are given in Figure 5. The simulated beams have a more pronounced deformation than the experimental phase reference beam. Yet, the image analysis routine yields the correct result for each of them, i.e., the algorithm finds the bin corresponding to the angle φ that was used to calculate the respective mode superposition.

Timing and noise issues.
The reference light is running backwards through the interferometer, but a tiny fraction of it is scattered towards the single-photon detectors at the outputs. To avoid any strong illumination of the detectors and to reduce potential noise in a first step, a very low power for which our image analysis routine found ϕ = 48 • . On the right of the figure, we show a selection of simulated deformed beams. The beams were computed as the sum of two ideal LG beams with a relative phase φ. In the top line, we simulated the effect of a small mismatch in the position of the center of the beams. In the bottom line, we simulated the effect of a small angle between the beams. Both defects yield a deformed HG mode in which one lobe is brighter than the other. This asymmetry is also present in the experimental reference beam on the left, but it is less pronounced than in the displayed simulations. When the simulated images were fed to the image analysis routine, the output was the correct φ value indicated above the pictures.
for the reference beam is required. The power chosen (∼ 2 nW) and the camera exposure time (∼ 100 ms) finally resulted from a compromise between the reduction of this noise source and the recording of an image within a time shorter than the typical phase drift. Already at this low power level, the detector count rate due to the scattered reference light (10 6 Hz) is already close to saturation. During the signal measurements, we therefore interrupt the reference beam using acoustooptic modulators. This interruption is short compared to the duration of the experimental cycle, which is in turn much shorter than the exposure time of the reference camera, so this has no influence on the reference image acquisition. To illustrate this with numbers from our measurement presented later: the camera is taking a new reference image every 125 ms, with an exposure time of 100 ms. It is not synchronized to our experimental cycle, which interrupts the reference light every 15 ms for 3 ms to perform the APD measurements, leading effectively to 20% less average intensity on the camera.
While the reference light was on, we switched off (gated) the single-photon detectors. Even in this configuration, we still noted an influence of the reference light on the background counts within the measurement window: when the reference light was completely blocked, the dark count rate was at about 80 Hz. With the reference light switched on as described above, the dark noise increased to 200 − 250 Hz within the measurement window. We excluded light leaking through the AOMs as the cause of the increased background count rates by additionally switching the reference light with mechanical shutters in a test. The dark count rate of the detector decreased over tens of milliseconds after turning off the relatively strong light exposure. Phenomenologically, the decay might be described by a stretched-exponential function [33]. This behavior has already been observed in other experiments, e.g. [34], and it might be explained by delayed afterpulses of the avalanche photodiodes [35].

Polarization and wavelength of the reference beam.
We show now that it is highly desirable to use the same polarization and the same wavelength for the reference light as for the signal.
First, the polarization-maintaining fibers and fibered beam-splitter, but also dielectric mirrors are birefringent, so the interferometer phase will be different for signal and reference if they have different polarizations. Furthermore, the birefringence changes with the mechanical stress of the fibers and with the temperature drifts, thus this difference will not stay constant.
Second, if the two optical paths differ geometrically by ∆L, two beams at different wavelengths will accumulate different interferometer phases. The variation of the measured interferometer phase around a certain wavelength λ = c/ν can be easily calculated as: where ν is the light frequency. So with a path difference ∆L on the order of a few centimeters, a difference by a few hundreds of MHz in the light providing the phase reference beam is already enough to change the inferred value of ϕ by several tens of degrees. Finally, dispersion can also play a significant role since a part of the interferometer is fibered. In a single-mode fiber, the change of the effective refractive index is dominated by the dispersion of the material [36]. We can thus estimate the dispersion dn/dν of our silica fibers to be on the order of 10 −3 /100 nm for our fibers [37]. In a perfectly symmetrical situation, the first-order contribution of dispersion vanishes. However, even if the optical path lengths are precisely equal, there might be a difference in their composition in terms of free-space and fibered lengths. Let us call this difference in fiber length ∆L fib . With this, we find: This effect will be smaller than the previous one for typical configurations, but can still play a role if signal and reference are separated by several GHz. If the frequencies of signal and reference are different, but stay constant, these two contributions lead, first, to a constant offset that could be determined, and second, to a different change of the phase when varying the path length. The latter difference is proportional to the relative wavelength difference and can thus in many cases be neglected for close wavelengths. As soon as the frequencies vary however, especially with respect to each other, the correlation between signal and reference interferometer phase will be lost. We therefore avoided these problems by using light from the same source as the one used for the signal state to be measured.
In this section, we have presented a detailed description of our OAM detection setup, and discussed the relevant issues relative to the phase-reference light. We now turn to the calibration procedure and provide the main figures of merit that have been measured thoroughly.

IV. ALIGNMENT PROCEDURE AND CALIBRATION
The detection setup presented here is very sensitive to the incoming beam position on the hologram dislocation and to its direction. The required fine tuning allows to calibrate and assess the performance of the setup, as shown in the following.

A. Optimizing couplings and cross-talks
The goal of the alignment procedure is twofold: an optimal coupling of the respective LG mode into the fiber of the respective branch and a strong rejection of all the other modes, as demanded by the mode selection requirement. The alignment is performed using classical bright light fields aligned with the signal to be analyzed later, their path being defined by the same single-mode fiber. A spatial light modulator is inserted on the way and enables to send various spatial modes into the setup.
In a first step, the position of the holograms is set by monitoring the intensity distribution with a camera placed behind the hologram. This distribution strongly depends on the relative position of the incoming beam and hologram center. Sending in a LG l=+1 mode, we can thus optimize the observed intensity distribution to be close to a TEM 00 profile in the R path (Table IV A). In the same way, the position of the L path hologram is optimized by sending in a LG l=−1 mode. Even small deviation of the hologram by a few micrometers become clearly visible. Then, the coupling into the fibers is optimized. Using two mirrors behind each of the holograms, we are able to adapt the mode exiting the hologram by maximizing the coupling efficiency up to 80%.
The next stage consists in sending a TEM 00 mode and using the same mirrors to now minimize its coupling. The rejection pattern of the setup is more pronounced than the acceptance and thus allows a better approach to the optimal point. Finally, a random search in the region around the found optimum allows to do some fine tuning. Here, all 6 degrees of freedom (2 transversal positions of the hologram and 4 directions for the 2 mirrors) are slightly varied while switching between the coupled mode and the unwanted modes (such as TEM 00 ). This way, the ratio η others /η LG is minimized (where η Λ is the coupling efficiency of a mode Λ into the single-mode fiber), while the coupling η LG of the targeted mode is kept at or close to the maximum.
Each path is finally characterized by measuring the couplings for the different LG modes sent into the setup. In the R path for example, we record the transmission of the mode corresponding to that path (LG l=+1 ), the two neighbors in l number (LG l=+2 and LG l=0 = TEM 00 ), and the mode corresponding to the other path (LG l=−1 ). An average rejection of 17 dB for the next neighbors (2% transmission) is typically obtained. For the suppression for the opposite mode (at ∆l = 2), a rejection of 23 dB was obtained in the worst case, and up to 37 dB in the best case. The typical value is 25 ± 3 dB [43]. An example of detailed coupling figure of merit is given in Table IV A.

B. Calibration of the fringes
An additional characterization is performed by sending classical beams carrying the four Hermite-Gaussian modes H, V, D and A, as shown in Figure 6. Theoretically, all four modes should lead to the same power of coupled light into both fibers. However, the power balance is not strictly mode-independent. Most of this imbalance can be explained by the imperfect mode filtering given in Table IV A. In the R path for instance, 82% transmission for mode LG l=+1 and 0.5% for mode LG l=−1 make

Mode at input
LG l=+1 TEM00 LG LG l=−1 ) is converted into a TEM00 mode in the far field of the R path (resp. L path), while other modes are converted to higher l-valued modes with doughnut shapes. The goal of this alignment is to maximize the subsequent coupling of the desired mode in the optical fiber and to minimize the coupling of all other modes. (b) Coupling efficiency in each path for various Laguerre-Gaussian modes at the input.
up for a ±6% difference in transmission between different HG modes. These mode-selective losses decrease the count rate and can additionally lead to a reduction of the fringe visibility, leading in turn to a decrease in the observed fidelity (see section IV D).
Finally, the phase reference circuit is started and the interferometer phase slowly scanned. Sending the modes H, V, D and A, the respective fringes are measured and correlated to the interferometer phase obtained from the reference. These fringes allow to check first for a good fringe visibility (typical values above 93% are obtained, mainly limited by the residual imbalance in the coupling of the HG modes and by the background noise) and second for the correct phase relation between signal and reference. Figure 6 illustrates this process. The power detected at the output X should exhibit a sinusoidal dependence in ϕ = 2α d + π (eq. 4). The condition for minimum power, i.e. ϕ = θ + π ↔ 2α d = θ resulting from equation 6, allows to deduce the value of θ from the position of the fringes. We checked that this was indeed the case and found good agreement within ±3 • .
The setup enables thus to accurately project an input state on various target modes. We now discuss the overall detection efficiency of the setup.  6: (a) Fringes measured during the calibration procedure for the various HG modes at the input. Experimental points (black dots) corresponding to the power received at the output X are given together with a sinusoidal fit (solid red line). The horizontal axis corresponds to twice the angle α d of the dark line in the phase reference images. It is related to ϕ by the relation ϕ = 2α d + π. As expected, the fringes shift as the various modes are sent. The position of the minima (indicated by arrows) gives the position where θ + ϕ ≡ π (mod 2π). (b) From this measurement, we can locate the four HG modes on the equator of the Bloch sphere (seen from top). The fringe positions measured in (a) are shown as red lines. The blue cross is obtained by averaging the deviations of the measured angles (red lines) from the respective theoretical angle (gray lines). In other words, it is the cross whose arms are the closest possible to the red lines while enforcing 90 • between its arms. It shows an excellent agreement with the theoretical θ values (0, π/2, π and 3π/2 respectively for modes H, D, V and A). In this particular example, the fringe visibility reaches 93 ± 2%, and the inferred θ values were respectively -0.1 • , 86.5 • , 176.1 • and 268.6 • for the modes H, D, V and A, leading to an average deviation of 2.2 ± 1.6 • . Let us note that, during the calibration measurement, the contamination with reference light is stronger than during the APD measurements for technical reasons. Therefore, the visibility is slightly reduced here in comparison with the data shown later.

C. Detection efficiency
The overall loss of the device comes from various contributions: • Input-beam splitter and filtering: ∼ 50%. The combination of the 50/50 beam splitter at the input and the subsequent mode selection causes a 50% fraction of the signal to go "into the wrong path" whatever the input mode. This fraction is filtered out and lost. It does not account for the quality of the coupling into the mode selecting fibers.
• Beam splitter for back-propagating the phase reference beam: ∼ 25% losses. Larger splitting ratios, such as 90:10 or even 99:1 would reduce these losses, which are thus not intrinsic to our scheme.
In our implementation, the detection efficiency for the HG modes was slightly lower than for the LG modes. This was due not to the detection setup but to the imperfection in the mode preparation. Indeed, the probe modes were generated by a reflection on a phase-only spatial light modulator. Yet, this method is limited by the intrinsic mode overlap between an ideal (target) mode and the sharply phase-modulated beam with gaussian envelope produced by the SLM. This overlap is lower for HG modes (64%) than for LG modes (78%) [32]. The rest of the power goes to higher-order p > 0 modes that are not detected due to spatial filtering. In practice, this generates some additional mode-selective losses but cannot degrade the quality of the tomography (visibility and then fidelity), as they are the same for each mode in every pair of orthogonal modes.

D. Impact of setup imperfections on the measurable fidelity
Experimental imperfections prevent from measuring unit fidelities even for ideal qubit states at the input. In this section, we discuss the upper limits on the fidelity values that can be measured.
The fidelity for HG states is especially sensitive to a reduced visibility V , setting a limit of F max = 1 /2(1 + V ). The reduced visibility can originate from a contamination with background noise, an imperfect mode filtering or a mode-dependent fiber-coupling. An imbalance of ±∆η in the fiber coupling of orthogonal HG modes will decrease the visibility of the fringe by approximately (∆η) 2 /2. A 99% fringe visibility, as achieved in some of the reported results without noise subtraction, leads to a maximum fidelity of 99.5%.
For LG modes at the entrance, no fringe should be present in the ideal case when scanning ϕ. Therefore, their fidelity is insensitive to the previous visibility reduction. However, a small relative leakage of mode L in path R (and vice versa) leads to a spurious fringe with a visibility 2 √ . This leakage translates into a maximum fidelity equal to 1− . With a leakage of less than −25 dB, the resulting error on the fidelity is limited to a fraction of a percent.
Finally, a shift in the interferometer calibration by an angle ∆θ will also decrease the maximum fidelity by an amount (∆θ) 2 . Even for a ∆θ on the order of 5 • , the fidelity decreases thus by less than a percent. As can be seen from Figure 6, this phase shift was controlled to better than 3 • in our implementation. In summary, this section has shown the reliability of our setup for OAM qubit measurements. Before moving on to some example of measurements performed with (approximate) single photons, let us recall three key parameters: • Cross-talk suppression (in LG basis): ∼ 25 dB, • Precision in phase measurement: better than 3 • , • Detection efficiency (excluding detectors): 24% ± 3%.

V. EXAMPLE OF QUANTUM STATE TOMOGRAPHY
We will now give some examples of density matrix reconstructions performed with our detection setup on weak coherent pulses in the single-photon regime. In this section, the setup is regarded as a black-box, yielding photon counts at the different outputs according to the interferometer configuration as described in Table I. Only data from output X is used for simplicity.
The data presented here results from measurements on weak coherent states with a mean photon-number of 0.6 photons per pulse. The spatial modes were imprinted by reflection on a spatial light modulator and were tuned to all of the six modes |R , |L , |H , |D , |V and |A . With the device in the HG bases configuration (i.e. measuring a fringe), three million measurements were made for each of the HG input modes, while one million measurements were performed for the LG input modes. In the configurations for measuring in the LG basis (i.e. blocked interferometer arm), one million measurements were performed for each of the two basis states, irrespective of the input state.
From these measurements, the density matrix of each input state was reconstructed using the formulae given in section II : S 1 = p |0 − p |1 , S 2 = p |+ − p |− and S 3 = p |+i − p |−i , where p |Ψ is the probability to measure the qubit in the state |Ψ . In Table IV, we give the density matrix parameters resulting from these calculations, and the fidelity to the ideal targeted state. In Figure 7, we give a more graphical view of these data for the three modes |R and |D and |H . Fidelities from 96% to 99% are measured, mainly limited by the imperfect state preparation and the background noise, which has not been corrected here.  TABLE IV: Stokes parameters and state fidelities extracted from the analysis of coherent pulses with a mean photonnumber of 0.6. In this particular example, the reconstruction of mode |A suffers from the coupling imbalance between paths R and L and from our imperfect mode preparation. The fidelity of the other modes is mostly limited by the background noise, which has not been subtracted here. The error bar on the fidelities is on the order of ±1%. For the sake of completeness, we finally give in the following some detailed values for this set of data: • The measured average signal raten sig (i.e., the rate at the average height of the fringes) was 0.03 clicks/measurement.
• The background noise originated from various sources. The intrinsic source is the APD dark noise of about 200 Hz, as discussed in Sect. III B 3. The extrinsic sources are related to the environment in which our detection setup was used (a quantum memory experiment with at least one additional light beam during the measurements [21]). With about 1 kHz, they constitute the major background contribution. An integration time of 800 ns per measurement gives thus a background rate of about n bg = 10 −3 background events per measurement (n bg /n sig = 3%).
• While the total number of measurements is in the millions, the number of events in a single bin is much smaller, thus a significant Poissonian fluctuation occurs. The angular bins at different ϕ values were only roughly equally covered by measurements, showing a standard deviation of up to 35%. The average numbers of measurements per bin (statistical errors) were 560 (4%) for the LG input states and 1400 (3%) for the HG input states.
The interferometric scheme and the related calibrations enable thus to reconstruct the density matrix of qubits encoded in OAM states with fidelities close to unity. In the following, we discuss how to extend this approach to higher-dimensional space.  9: A possible alternative implementation for extending the device to higher-dimensional Hilbert spaces, using mode sorting refractive optics R1 and R2 as developed in [29,31]. The input beam-splitters and subsequent holograms can be replaced by the combination of a circulator (to separate the signal and the phase reference beam) and an OAM mode sorter made of two refractive elements R1 and R2. The reformatter R1 performs a log-polar to Cartesian coordinate mapping and the phase corrector R2 corrects for the different propagation lengths from R1 to R2. These two elements transform co-propagating LG beams with different l values into approximate TEM00 beams with a propagation direction tilted by an angle proportional to l. The fan-out separation enhancer is a specific periodic phase-only hologram that increases the angular separation between the modes by decreasing the residual mode-overlap. After being separated and losing their specific spatial shapes, the different OAM components can be brought to interfere in an array of beamsplitters forming a mesh of interferometers with controllable phases, then directed towards single-photon counting modules. A good timing of the phase reference should be ensured with appropriated shutters inserted in the interferometer array in order to record all the relevant phase differences. The superposition |l = +6 + |l = −6 is shown at the input for illustrative purpose.

VI. PERSPECTIVES AND POSSIBLE EXTENSIONS
To reach higher-dimensional space, one can introduce additional beam-splitters after the first input one. Figure 8 shows the example of an extended setup to perform quantum state tomography in a four-dimensional Hilbert space spanned by the modes LG l=−3 , LG l=−1 , LG l=+1 and LG l=+3 . In each of the subsequent paths, a mode projector on a different OAM value is inserted. Keeping an OAM difference ∆l = 2 between different modes ensures a better mode filtering. Fibers at the end of the mode projectors are connected to a series of cascading fiber beam-splitters, creating an array of nested Mach-Zehnder interferometers. Alternatively, these nested interferometers could be engraved in photonic circuitry [38][39][40], which would also provide greater simplicity and better phase stability. Similarly to the two-dimensional setup, a phase reference beam is sent backwards. The various unused output ports allow the imaging of the phase reference beams, and the determination of all the relevant phases. The previously described image-analysis routine (III B 2) can be directly used to compute the relevant phases between pairs of modes, given the fact that the phase reference is timed in order to image only two mode superpositions.  V: Expected performances of two possible extensions of the current device ("Current") to higher dimensional Hilbert spaces. The current device enters into the category "Beamsplitter cascade", although the cascade is only one step deep. The device "3BS" refers to the setup shown in Fig. 8 and extends the current device to a two-step cascade of beamsplitters. "OAM sorter" refers to the device shown in Fig. 9 exploiting a mode sorter scheme. Detection efficiencies does not take into account the intrinsic efficiency of the singlephoton counters.
One more series of beam splitters would allow quantum tomography in an eight-dimensional space, but it would become challenging. Indeed, the addition of more beam-splitters will result in degrading the count rates exponentially in the number of beam splitters (linearly in the number of detected modes). Also the phase measurement and/or stabilization can become a more serious issue if the dimension of the Hilbert space increases.
Another way to reach higher-dimensional spaces would be to take advantage of the recently developed OAM mode-sorting techniques [29][30][31]41] where OAM states are converted into transverse momentum states. As shown in Figure 9, the input beam-splitter can indeed be replaced by an OAM mode sorter made out of refractive elements as described in [29] with a separation enhancer [31]. This combination would perform both the separation and mode conversion in the same time, thus largely improving the detection efficiency and versatility. Indeed, the detection efficiency would remain (almost) constant as the number of detected modes increases. The reduction in the required number of optical elements may also provide a better phase stability. Table V summarizes the performances one can foresee using some realistic parameter estimations extracted from [41].

VII. CONCLUSION
In conclusion, we have presented and quantitatively characterized a setup enabling the quantum state tomography of a photonic qubit encoded in the orbital angular momentum degree of freedom. This interferometric setup performs the projection of an input state over different mutually unbiased bases and enables the reconstruction of the full density matrix. Various noise and imperfection sources have been discussed and the benchmarks of the experimental implementation have been given. We have finally proposed some extensions of our scheme to higher-dimensional Hilbert spaces. The combination of the presented technique with the recent and elegant mode-sorting methods can indeed be applied to protocols requiring OAM qudit tomography.