Experimental minimum-error quantum-state discrimination in high dimensions

TL;DR

This paper demonstrates the first experimental implementation of minimum-error quantum-state discrimination in high dimensions, achieving high accuracy for nonorthogonal state identification using optical transverse modes.

Contribution

It provides the first experimental realization of minimum-error measurement for high-dimensional nonorthogonal states, advancing quantum communication techniques.

Findings

Achieved high fidelity in discriminating states up to 21 dimensions

Tested nearly 14,000 states with deviations below 3.6%

Validated the scheme's suitability for quantum communication protocols

Abstract

Quantum mechanics forbids perfect discrimination among nonorthogonal states through a single shot measurement. To optimize this task, many strategies were devised that later became fundamental tools for quantum information processing. Here, we address the pioneering minimum-error (ME) measurement and give the first experimental demonstration of its application for discriminating nonorthogonal states in high dimensions. Our scheme is designed to distinguish symmetric pure states encoded in the transverse spatial modes of an optical field; the optimal measurement is performed by a projection onto the Fourier transform basis of these modes. For dimensions ranging from D = 2 to D = 21 and nearly 14000 states tested, the deviations of the experimental results from the theoretical values range from 0.3% to 3.6% (getting below 2% for the vast majority), thus showing the excellent performance of our scheme. This ME measurement is a building block for high-dimensional implementations of many quantum communication protocols, including probabilistic state discrimination, dense coding with nonmaximal entanglement, and cryptographic schemes.

Figures

Tables

Table 1. Table 1: Summary of the sets of states tested.

$D$	$N$	Sets/$N$	No. of sets	No. of states
2	2, 3, 7	15	45	180
3	3, 5	225	450	1800
4	4, 6	60	120	600
5	5, 7	80	160	960
6	6, 7	100	200	1300
7	7, 11	120	240	2160
8	8, 12	140	280	2800
9	9, 13	160	320	3520
10 to 21	$N\ge D$	1	3	$\sim 600$
Total			1851	$\sim 13920$

Equations

$$|{\psi }_j⟩=\sum _{n=0}^{D-1}{c}_n{\omega }^{jn}|n⟩,$$

$${\hat{\Pi }}_k^{\mathrm{ME}}={\hat{\mathcal{F}}}_N|k⟩⟨k|{\hat{\mathcal{F}}}_N^{-1}\equiv |{\mu }_k⟩⟨{\mu }_k|,$$

$$P_{\mathrm{corr}}=\frac{1}{N}\sum _{j=0}^{N-1}P({\delta }_j|{\psi }_j)=\frac{1}{N}{\left(\sum _{j=0}^{D-1}{c}_j\right)}^2.$$

$$\hat{\Pi }(x_k,f)\equiv {\hat{\Pi }}_k^{\mathrm{ME}},\mathrm{for}{x}_k=-\lambda f{m}_k/dN.$$

$$t(x,y)=\mathrm{rect}\left(\frac{y}{L}\right)\sum _{\ell =0}^{D-1}{t}_{\ell }\mathrm{rect}\left(\frac{x-{\xi }_{\ell }d}{2a}\right),$$

$$t_{\ell }=$$

$$c_{\max }=$$

$${\hat{\rho }}_{\mathrm{tomo}}=\sum _{\mu =0}^2\sum _{j=0}^1{p}_{\mu j}|b_{\mu j}⟩⟨b_{\mu j}|-\hat{I},$$

$$[{\mathcal{I}}_{\mathrm{est}}$$

$$={\mathcal{I}}_{\max }{\left|\sum _{j=0}^{D-1}\sqrt{I_j}{e}^{i{\varphi }_j}{e}^{ikjd{x}_{\kappa }/f}\mathrm{sinc}\left(\frac{ka{x}_{\kappa }}{f}\right)\right|}^2,$$

$$\Delta (\{{\varphi }_k^{\star }\},{\mathcal{I}}_{\max }^{\star })=$$

$$\Delta (\{{\varphi }_k\},{\mathcal{I}}_{\max })=$$

$$|{\psi }_{\mathrm{est}}⟩\propto \sum _{j=0}^{D-1}\sqrt{I_j}{e}^{i{\varphi }_j^{\star }}|j⟩.$$

$$|{\psi }_j⟩=\sum _{k=0}^{D-1}{c}_k{\omega }^{jk}|k⟩,j\in \left\{0,\cdots ,N-1\right\},$$

$$\displaystyle t_{j}=\begin{cases}\displaystyle\sqrt{1-\sqrt[D]{\frac{j-j_{0}+1}{D-j_{0}}\alpha}}~{}~{}~{}~{}&~{}{\rm if}~{}~{}j\geq j_{0}\\ {\,\,}&}{\textstyle\\ 1~{}~{}~{}~{}&~{}{\rm if}~{}~{}j<j_{0}\\ \end{cases},$$

$$c_j=\frac{t_j}{{\mathcal{N}}_{\alpha ,j_0}},{\mathcal{N}}_{\alpha ,j_0}=\sqrt{\sum _{k=0}^{D-1}{t}_k^2},$$

$$m_k=\{\begin{array}{ll}k& \mathrm{if}0\le k\le \frac{N}{2}\\ k-N& \mathrm{if}\frac{N}{2}<k\le N-1\end{array},$$

$$P({\delta }_k|{\psi }_j)=\frac{I_{kj}}{{\sum }_{l=0}^{N-1}{I}_{lj}},$$

$$P_{\mathrm{corr}}=\frac{1}{N}\sum _{j=0}^{N-1}P({\delta }_j|{\psi }_j),$$

Full text

Quantum mechanics establishes fundamental bounds to our capability of distinguishing among states with nonvanishing overlap: if one is given at random one of two or more nonorthogonal states and asked to identify it from a single shot measurement, it will be impossible to accomplish the task deterministically and with full confidence. This constraint has deep implications both foundational, underlying the debate about the epistemic and ontic nature of quantum states Spekkens07 ; Leifer14 ; Branciard14 ; Ringbauer15 , and practical, warranting secrecy in quantum key distribution Gisin02 ; Grosshans03 . Beyond that, the problem of discriminating nonorthogonal quantum states plays an important role in quantum information and quantum communications Reviews .

A wide variety of measurement strategies have been devised in order to optimize the state discrimination process according to a predefined figure of merit Reviews . The pioneering one was the minimum-error (ME) measurement Holevo73 ; Yuen75 ; HelstromBook where each outcome identifies one of the possible states and the overall error probability is minimized. Other fundamental strategies conceived later Ivanovic87 ; Chefles98 ; Croke06 ; Jimenez11 ; Sugimoto12 ; Herzog12 ; Bagan12 employ the ME discrimination in the step next to a transformation taking the input states to more distinguishable ones Nakahira12 , which enables us to identify them with any desired confidence level (within the allowed bounds) and a maximum success probability. Nowadays, the ME measurement is central to a range of applications, including quantum imaging Tan08 , quantum reading Pirandola11 , image discrimination Nair11 , error correcting codes Lloyd11 , and quantum repeaters Loock06 , thus stressing its importance.

Closed-form solutions for ME measurements are known only for a few sets of states. One of these is the set of symmetric pure states (defined below) prepared with equal prior probabilities Ban97 . Discriminating among them with minimum error sets the bounds on the eavesdropping in some quantum cryptographic schemes Phoenix00 and is crucial for optimal deterministic and probabilistic realizations of protocols like quantum teleportation Neves12 , entanglement swapping Prosser14 , and dense coding Pati05 , when the quantum channels are nonmaximally entangled.

Experimental demonstrations of ME discrimination have been provided in continuous variables for two coherent states Cook07 , and in two-dimensional Hilbert spaces for sets of two Barnett97 , three, and four states Clarke01 encoded in the light polarization, and two states encoded in the 14N nuclear spin Waldherr12 . It is of key importance to extend this to high-dimensional quantum systems (qudits) due to the many advantages they offer over qubit-based applications. Qudits provide an increase in the channel capacity for quantum communication Fujiwara03 , and a higher error rate tolerance and improved security in quantum key distribution Cerf02 ; Huber13 ; Sheridan10 . Moreover, they are a necessary resource for fundamental tests of quantum mechanics, like contextuality tests Kochen67 , and for lowering the detection efficiencies required for nonlocality tests Vertesi10 . Therefore, the ability to perform ME measurements for qudit states will enhance their potential of use in many of these practical applications. It will also be essential for exploring novel protocols of quantum-state discrimination Jimenez11 ; Sugimoto12 ; Herzog12 , bringing positive impacts for quantum communications.

In this Letter, we report the first experimental demonstration of ME discrimination among nonorthogonal states of a single qudit. This is done for equally likely symmetric states in dimensions $D$ ranging from $D=2$ to $D=21$. Using states encoded in $D$ transverse spatial modes of an optical field, known as spatial qudits Neves05 , we carried out the experiment for every dimension in that range in a total of 1851 sets of states. Up to small deviations caused by unavoidable experimental imperfections, our scheme is shown to be optimal, achieving the minimum error probability allowed by quantum mechanics.

To outline the problem, consider a set of $N$ quantum states $\{|\psi_{j}\rangle\}_{j=0}^{N-1}$ spanning a $D$-dimensional Hilbert space $\mathcal{H}_{D}$, with $D\leq N$, defined by

[TABLE]

where $\omega=\exp(2\pi i/N)$, the $c_{n}$’s are real and nonnegative ($\sum_{n}c_{n}^{2}=1$), and $\{|n\rangle\}_{n=0}^{D-1}$ is an orthonormal basis in $\mathcal{H}_{D}$. They are symmetric under $\hat{U}=\sum_{l=0}^{D-1}\omega^{l}|l\rangle\langle l|$ since $|\psi_{j}\rangle=\hat{U}|\psi_{j-1}\rangle=\hat{U}^{j}|\psi_{0}\rangle$ and $|\psi_{0}\rangle=\hat{U}|\psi_{N-1}\rangle$. If each of these states is prepared with the same a priori probability $1/N$, they can be identified with minimum error through the measurement Ban97

[TABLE]

where $k\!=0,\ldots,N-1$, and $\hat{\mathcal{F}}_{N}\!=\!\frac{1}{\sqrt{N}}\sum_{m,n=0}^{N-1}\omega^{mn}|m\rangle\langle n|$ is the discrete Fourier transform acting on an $N$-dimensional Hilbert space $\mathcal{H}_{N}$. For $N=D$ (linearly independent states), this is a projective measurement on $\mathcal{H}_{D}$. For $N>D$ (linearly dependent states), this is a projective measurement on $\mathcal{H}_{N}$ that realizes the positive operator valued measure (POVM) for ME discrimination on $\mathcal{H}_{D}$. This procedure is based on Neumark’s theorem HelstromBook for implementing POVMs by embedding the system space $\mathcal{H}_{D}$ into a larger Hilbert space $\mathcal{H}_{N}$ given by the direct sum $\mathcal{H}_{D}\oplus\mathcal{H}_{N-D}$, where $\mathcal{H}_{N-D}$ represents the $N-D$ unused extra dimensions of the original system Chen07 . From Eqs. (1) and (2), the probability of obtaining an outcome $\delta_{k}$ (associated with $\hat{\Pi}_{k}^{\rm ME}$) if the prepared state was $|\psi_{j}\rangle$ is given by $P(\delta_{k}|\psi_{j})=|\langle\mu_{k}|\psi_{j}\rangle|^{2}$. Thus, the maximum overall probability of correctly identifying the state will be

[TABLE]

Equivalently, $P_{\rm err}=1-P_{\rm corr}$ will be the ME probability.

The experimental setup used to optically demonstrate the ME discrimination among symmetric states is shown and described in Fig. 1(a). We use a bright beam from a single-mode laser as our light source, which is usual in most optical implementations of quantum-state discrimination Barnett97 ; Clarke01 ; Mohseni04 ; Mosley06 . The transverse spatial profile of this beam is proportional to the transverse probability amplitude of a single-photon multimode field. With the procedure described next, it will mimic the quantum state (1) of a single-photon qudit. The arrangement within the light shaded region enables the preparation of arbitrary pure states of spatial qudits, as recently demonstrated Solis-Prosser13 . This is possible by modulating the transverse spatial profile of the incoming field as follows: an array of $D$ rectangular slits with a blazed diffraction grating inside each one is displayed on the screen of the phase-only spatial light modulator. [Figure 1(b) shows a typical phase mask for $D=3$.] In the back focal plane of a converging lens, the beam portions impinging within the slit zones are diffracted into different orders ($0,\pm 1,\ldots$); the portions impinging outside the slit zones go to the zeroth diffraction order. If we choose one of the high orders, say $+1$, to prepare the states, the amplitude of their complex coefficients is obtained by controlling the amount of light diffracted by each slit to that order, which is a function of the phase modulation depth of the gratings. The relative phases are set either by adding a constant phase value to the gratings Solis-Prosser13 or by lateral displacements of them Varga14 (here we interchanged between these two procedures Supplement ). Finally, a slit diaphragm filters out the $+1$ diffraction order from the others and the emerging light will be a coherent superposition of the which-slit modes $\{|n\rangle\}_{n=0}^{D-1}$ modulated by proper complex coefficients, thus representing the desired qudit state of Eq. (1). A schematic of the preparation stage is shown in Fig. 1(b). The inset presents an intensity measurement at the image plane for a state $|\psi_{j}\rangle$ prepared with the phase mask shown there.

The ME measurement is performed through spatial postselection. Let $\hat{\Pi}(x,z)=|\mu(x,z)\rangle\langle\mu(x,z)|$ denote the measurement operator associated with a pointlike detector in the transverse position $x$ and a longitudinal distance $z\in[f,2f]$ from a converging lens of focal length $f$. In an $N$-dimensional space this detector postselects the state $|\mu(x,z)\rangle\propto\sum_{\ell=0}^{N-1}\varphi_{\ell}(x,z)|\ell\rangle$, where the complex coefficients $\varphi_{\ell}(x,z)$ are given by the Fresnel diffraction integral of slit $\ell$ calculated in $(x,z)$ Solis-Prosser11 . If the detector is placed at the focal plane of the lens ($z=f$), it will postselect the state $|\mu(x,f)\rangle=\frac{1}{\sqrt{N}}\sum_{\ell=0}^{N-1}\omega^{xdN\ell/\lambda f}|\ell\rangle$, where $d$ is the slit separation and $\lambda$ is the light wavelength. Now, consider an array of $N$ detectors at $z=f$ distributed along the transverse positions $x_{k}=-\lambda fm_{k}/dN$, where $k=0,\ldots,N-1$, and $m_{k}=k$ if $0\leq k\leq N/2$ or $m_{k}=k-N$ if $N/2<k\leq N-1$ Supplement . It is easy to see that each detector in this array postselects the state $|\mu(x_{k},f)\rangle=\hat{\mathcal{F}}_{N}|k\rangle$. From Eq. (2), we have

[TABLE]

Therefore, with such a detection scheme we can implement the ME discrimination of symmetric states. From the discussion following Eq. (2), for $N>D$, the $D$-slit array can be viewed as an array of $N$ slits where $D$ input modes are used to encode the states and the remaining $N-D$ modes are in their respective vacuum states. The propagation through the lens system provides the unitary coupling between these two subspaces and the projective measurement (4) accomplishes the optimal POVM.

In our setup, the ME measurement is shown in the dashed box of Fig. 1(a). Figure 1(b) shows a schematic of this stage with an example of $D=3$ and $N=5$. Note that we used the same lens to assist the preparation and measurement stages described above. This was done to simplify the setup, as our main goal was to demonstrate the optimal discrimination process. The same results would be achieved if the detection were performed after letting the spatial qudit state propagate through an arrangement of lenses taking the Fourier transform of the which-slit modes (possibly having to rescale the $x_{k}$’s).

Our implementation comprised both linearly independent and dependent states in every dimension from $D=2$ to $D=21$. Their coefficient amplitudes, $\{c_{n}\}_{n=0}^{D-1}$ in Eq. (1), were parametrized as follows. For $D=2$, $c_{0}=\cos\frac{\theta}{2}$ with $\theta\in[0,\pi]$; for $D=3$, $c_{0}=\sin\frac{\theta_{1}}{2}\cos\frac{\theta_{2}}{2}$ and $c_{1}=\sin\frac{\theta_{1}}{2}\sin\frac{\theta_{2}}{2}$ with $\theta_{1(2)}\in[0,\pi]$. These hyperspherical coordinates enabled the experiment to be performed for arbitrary symmetric states in those dimensions. However, for $D\geq 4$, this approach would prevent us to represent the results in a single plot. Thus, for $D=4$ to 9 we defined $c_{n}^{2}\propto 1-\sqrt[D]{\frac{n-j_{0}+1}{D-j_{0}}\alpha}$ if $n\geq j_{0}$, and $c_{n}^{2}\propto 1$ if $n<j_{0}$, where $j_{0}=1,\dots,D-1$ and $\alpha\in[0,1]$. Using only two parameters $(j_{0},\alpha)$ we were able to test a large diversity of sets of symmetric states in those dimensions and, at the same time, to obtain surfaces for the probabilities in Eq. (3) with good contrast Supplement . This latter aspect ensured that the tested sets were well distinguishable from one another, since $P_{\rm corr}$ in the ME discrimination is also a measure of distinguishability among sets of states Chefles02 . As $D$ gets larger, to experimentally build up those surfaces became very time consuming. Thus, for $D=10$ to 21 we found it sufficient to generate a few sets of states (three for each $D$) to perform the discrimination. In order to avoid any bias in our choice, the coefficients were randomly selected Comment .

With the preparation and measurement stages outlined above, the experiment is carried out in the following way: we first define the set of symmetric states by specifying $D$, $N$, and the coefficient amplitudes $\{c_{n}\}_{n=0}^{D-1}$ [Eq. (1)]. Afterwards, we prepare one state of this set and perform the discrimination process on it by measuring, with a CMOS camera, the light intensity at the $N$ transverse positions defined in Eq. (4). This is exemplified in Fig. 1(b). The measured intensities are integrated over $y$ [this direction, also shown in Fig. 1(b), is not relevant for us], the background noise is subtracted, and a small compensation for the detection efficiency due to diffraction is applied Supplement . These procedures are repeated for each state in the set. Denoting $I_{kj}$ as the resulting intensity at detector $k$ when the input state is $|\psi_{j}\rangle$, the probability of identifying it (correctly or incorrectly) is obtained as $[P(\delta_{k}|\psi_{j})]_{\rm expt}=I_{kj}/\sum_{\ell=0}^{N-1}I_{\ell j}$, from which we estimate the overall probability of correct identification, $P_{\rm corr}$, using the expression in the middle of Eq. (3).

In total, the experiment was performed for 1851 sets of symmetric states, comprising nearly 14000 different states Supplement . Figures 2 and 3 show a collection of our results for the values of $N$ and $D$ specified there. In Fig. 2 we plot $P_{\rm corr}$ as a function of the parameter(s) of the state coefficients defined earlier. The experimental results are given by the squares and points while the theoretical predictions—computed from the rightmost side of Eq. (3)—by the solid curves and surfaces. The error bars in these graphs were smaller than the size of the data points. In Fig. 3 we plot $P(\delta_{k}|\psi_{j})$, defined above Eq. (3). Each graph shows these probabilities for a single set of states settled by the coefficient amplitudes shown in the inset. The red stems are the experimental results and the empty bars the expected theoretical values given by $[P(\delta_{k}|\psi_{j})]_{\rm theo}=|\langle\mu_{k}|\psi_{j}\rangle|^{2}$. We observe in both figures the close agreement between theory and experiment.

For each of the tested sets we calculated the root-mean-square deviation (RMSD) of the results from the optimal theoretical values. With the residuals defined by $R_{kj}\equiv[P(\delta_{k}|\psi_{j})]_{\rm theo}-[P(\delta_{k}|\psi_{j})]_{\rm expt}$, we have ${\rm RMSD}=\sqrt{\sum_{j,k=0}^{N-1}R_{kj}^{2}/N^{2}}$. Figure 4(a) shows a histogram of the obtained deviations, which ranged from $0.3\%$ to 3.6%. For 80% of the sets the RMSDs were below 2%. The errors in our experiment are caused mainly by aberrations of the optical elements and the imperfect preparation of the input states, which is, possibly, responsible for the largest deviations. For instance, examining the results of Fig. 2(i), it can be seen that the worse ones occurred for $j_{0}=1,2$ and some values of $\alpha\neq 0$. For this particular realization, we plot in Figs. 4(b) and 4(c) the average fidelities of preparation for the states in each set and the RMSD$\times(-1)$, respectively. It is clearly observed that the largest deviations from the theoretical predictions correspond to the sets of states with lower fidelities. Despite all this, the deviations in our experiment, as shown in Fig. 4(a), were considerably small taking into account that experimental implementations are always subjected to errors, the wide variety of states tested, and the high dimensions we have explored. Therefore, it is safe to say that our scheme is indeed optimal.

In conclusion, we have presented the first experimental demonstration of minimum-error discrimination among nonorthogonal states of a $D$-dimensional qudit. Our measurement scheme was designed to distinguish symmetric pure states, an important resource in quantum information Chefles98 ; Croke06 ; Jimenez11 ; Sugimoto12 ; Herzog12 ; Nakahira12 ; Ban97 ; Phoenix00 ; Neves12 ; Prosser14 ; Pati05 . Envisaging its application to quantum protocols, a single photon implementation of this experiment would require a source with a sufficiently large transverse spatial coherence, in order to prepare high-quality states. The preparation and measurement stages would be exactly the same described earlier. The detector array, however, will be required to have single-photon counting capability. Some candidates for this include charged-coupled device cameras (CCDs), either an intensified CCD Fickler13 or an electron multiplying CCD Zhang09 , and also single-photon avalanche diode arrays based on CMOS technology Scarcella13 ; Unter16 .

The ME measurement is central not only to many applications in quantum information, but also to put forward the implementation of “complete” probabilistic discrimination protocols in high dimensions: by iterating the discrimination process in case of failed attempts Neves12 , one can significantly increase the information gain about the input states. Our demonstration constitutes a building block for future realizations of these protocols and, accordingly, will benefit a variety of tasks in high-dimensional quantum information processing Phoenix00 ; Neves12 ; Prosser14 ; Pati05 .

Appendix A Preparation and characterization of symmetric states

A.1 State preparation and phase encoding

Let us consider a pure qudit state $|\psi\rangle=\sum_{\ell=0}^{D-1}c_{\ell}|\ell\rangle$, where $\sum_{\ell=0}^{D-1}|c_{\ell}|^{2}=1$. This state can be encoded onto the one-dimensional transverse spatial modes of an optical field by modulating its transverse profile with a transmission function $t(x,y)$ given by

[TABLE]

where $\xi_{\ell}=\ell-(D-1)/2$, $d$ is the distance between the centers of two neighboring slits, $2a$ is the width of each slit, $L$ is the slit length, ${\rm rect}(\bullet)$ are rectangle functions, and $t_{\ell}$ is the transmission coefficient of each slit, defined as

[TABLE]

The procedures indicated in Refs. Solis-Prosser13 ; Varga14 , based on the use of diffraction gratings generated by a phase-only spatial light modulator (SLM) and the selection of its first diffraction order, were followed in order to generate the states to be discriminated. Briefly, the modulus of the transmission function $|t(x,y)|$ is encoded by a diffraction grating with variable phase modulation depth. Afterwards, the complex phase of the transmission function, ${\rm Arg}\bm{(}t(x,y)\bm{)}$, is encoded either by adding the phase directly to the grating Solis-Prosser13 or by displacing the grating by a certain amount of pixels for each slit Varga14 . These methods are known as phase addition (PA) and grating displacement (GD), respectively. The aforementioned procedures were complemented with a convenient choice of the grating period and orientation Marquez07 and an appropriate configuration of the SLM—i.e. phase modulation response to the addressed gray levels Moreno04 and digital electric addressing sequences Lizana10 —in order to improve the efficiency and temporal stability of the modulated phase.

Figure 5 shows an example of $|t(x,y)|$, ${\rm Arg}\bm{(}t(x,y)\bm{)}$ and the corresponding phase masks addressed at the SLM (via PA or GD method) used to prepare a given spatial qudit state with $D=7$. States ranging from $D=2$ to $D=9$ were generated using the PA method, whereas GD method was used for $D=10$ to $D=21$. This choice was motivated by the quality of the states being generated, whose estimation is explained in Sections A.2 and A.3.

A.2 State characterization for qubits

For a spatial qubit, a complete tomography can be performed since an informationally complete set of measurements can be taken from four transverse positions of the detectors in the Fourier plane and two ones located at the image plane Solis-Prosser11 . For this reason, the optical setup included a CMOS camera situated at the image plane of the SLM in order to perform the state estimation, as Fig. 6 depicts. These six measurements complete a set of three orthogonal bases.

The first basis, corresponding to the states $|b_{0j}\rangle=(|0\rangle+e^{ij\pi}|1\rangle)/\sqrt{2}$ with $j=0,1$, can be measured at positions $x_{0j}=j\pi f/kd$ in the focal plane. The second basis contains the states $|b_{1j}\rangle=(|0\rangle+ie^{ij\pi}|1\rangle)/\sqrt{2}$ and can be measured at $x_{1j}=(j-\tfrac{1}{2})\pi f/2kd$, also in the focal plane. The third basis is the computational basis $|b_{2j}\rangle=|j\rangle$, which is measured at positions $x_{2j}=-(j-\tfrac{1}{2})d$ in the image plane (the minus sign at the beginning is consequence of the inversion of the image created by any lens system).

Once the intensities were recorded, the probabilities are obtained by dividing each intensity by the total intensity within each set, i.e., $p_{\mu j}=I_{\mu j}/(I_{\mu 0}+I_{\mu 1})$. The first estimation of the state is given by Ivanovic81

[TABLE]

where $\hat{I}$ is the identity operator on a two-dimensional Hilbert space. This method for state reconstruction relies on measurements performed on three mutually unbiased bases (MUBs) and, hence, it is commonly referred to as MUB tomography.

The operator $\hat{\rho}_{\rm tomo}$ defined in Eq. (6) is, by construction, an hermitian matrix. However, due to experimental imperfections, this matrix is not always positive definite and, hence, it might not represent a physical state. The same measurements taken can be used for a second procedure known as maximum likelihood estimation (MLE). This procedure is an optimization problem in which $\hat{\rho}_{\rm tomo}$ can be used as a starting point for the numerical optimization in order to obtain the appropriate density operator which most likely describes the experimental data James01 . This second estimated state, defined as $\hat{\rho}_{\rm MLE}$, is used afterwards to compute its fidelity with respect to the target state $|\psi\rangle$, given by $F(|\psi\rangle,\hat{\rho}_{\rm MLE})=\langle\psi|\hat{\rho}_{\rm MLE}|\psi\rangle$, in order to characterize the quality of the state preparation.

A.3 State characterization for qudits

Unlike the case of spatial qubits, an arbitrary $D$-dimensional (${D>2}$) qudit state cannot be properly reconstructed only from measurements at the Fourier and image planes of a single imaging system and such a task requires the assistance of an additional SLM Lima11 . Nevertheless, the states generated in our setup are nearly pure due to the high transverse coherence length of the laser beam, which is attested by the high visibility of the interference patterns. Thus, an estimation can be performed by using least squares fitting as follows. We first write the state to be estimated as $|\psi_{\rm est}\rangle=\sum_{j=0}^{D-1}a_{j}|j\rangle$, where $|j\rangle$ is the state of the transverse optical mode associated with $j$th slit. The intensities $I_{j}$ of each slit at the image plane allows us to obtain the complex modulus of each coefficient $a_{j}$, according to $|a_{j}|^{2}=I_{j}/(\sum_{k=0}^{D-1}I_{k})$. This allows to rewrite the estimated state as $|\psi_{\rm est}\rangle\propto\sum_{j=0}^{D-1}\sqrt{I_{j}}e^{i\phi_{j}}|j\rangle$. The expected far field intensity distribution for this state is given by

[TABLE]

for every recorded transverse position $x_{\kappa}$ in the focal plane, where $\mathcal{I}_{\max}$ is a global proportionality constant for each experimental intensity distribution, and ${\rm sinc}(z)=\sin(z)/z$. Then, the phase retrieval is performed from the experimentally obtained far field intensity distribution $[\mathcal{I}_{\rm exp}]_{\kappa}$ by least squares according to

[TABLE]

where

[TABLE]

and $\{\phi_{k}^{\star}\}$ and $\mathcal{I}_{\max}^{\star}$ are the parameters that allow to minimize $\Delta\bm{(}\{\phi_{k}\},\mathcal{I}_{\max}\bm{)}$. Once the minimization has been carried out, the estimated state is completed by the phases that solve the problem of Eq. (8), i.e.,

[TABLE]

This state is then used to characterize the quality of the preparation through the fidelity $F(|\psi_{\rm est}\rangle,|\psi\rangle)=|\langle\psi_{\rm est}|\psi\rangle|^{2}$ between the estimated state and the target state $|\psi\rangle$ intended to be prepared.

A.4 State parametrization for $D=4$ to $D=9$

In general, we studied the case of discrimination among $N$ states of the form

[TABLE]

where $c_{k}$ are normalized non-negative coefficients and $\omega=\exp({2\pi i/N})$. Coefficients $c_{k}$ can be easily parametrized by using hyperspherical coordinates. In dimension 2, the Bloch parametrization is recovered, as $|\psi_{j}\rangle={\cos(\theta/2)|0\rangle}{+\omega^{j}\sin(\theta/2)|1\rangle}$, where $\theta\in[0,\pi]$. The three dimensional case is easily parametrized as well: $|\psi_{j}\rangle=\sin(\theta_{1}/2)\cos(\theta_{2}/2)|0\rangle+\omega^{j}\sin(\theta_{1}/2)\sin(\theta_{2}/2)|1\rangle+\omega^{2j}\cos(\theta_{1}/2)|2\rangle$, with $\theta_{1,2}\in[0,\pi]$. And so on for larger dimensions. For each dimension $D$, we require $D-1$ angles $\theta_{m}\in[0,\pi]$ in order to be able of preparing arbitrary states of the form given by Eq. (11).

However, for $D\geq 4$, a practical problem arises: It is not possible to represent the results in a single plot as more than 2 angles are required. Additionally, the involved times increase exponentially since, considering $M$ as the number of values taken of each angle $\theta_{m}$, a complete discrimination of $N$ states for all fiducial states needs a number of measurements equal to $N\times M^{D-1}$, becoming extremely time-consuming for large dimensions.

In order to overcome this drawback, we devised a strategy in which these $D$-dimensional fiducial states are parametrized by only two parameters, making the graphic representations of the results feasible. The challenge lies in the way the parametrization is chosen, as we imposed that (i) it must provide a reasonable diversity of fiducial states, and (ii) the surface for the plot of the theoretical minimum-error probability must have a good contrast, i.e., a reasonable difference between the minimum and maximum expected values. This last requirement, beyond aesthetic reasons, is motivated by the fact that the ME discrimination probability can also be used as a distinguishability measure within a set of quantum states CheflesDist : A surface exhibiting few changes and no slopes could have been interpreted as having measured sets with very similar distinguishabilities. Hence, the experiments would have considered a very particular subset of states from the Hilbert space rather than being a proof-of-principle of the general applicability of this method for sets constructed from arbitrary fiducial states. We defined the following biparametric representation of the transmission coefficients

[TABLE]

with $\alpha\in[0,1]$ and $j_{0}=1,\dots,D-1$, where the state coefficients computed as

[TABLE]

accomplish the aforementioned requirements. The behavior of the transmission coefficients $t_{j}$ for each slit, where $t_{j}\propto|c_{j}|$, is depicted in Fig. 7 together with the expected optimal minimum error probability. As it is shown, the transmission of the first $j_{0}$ slits ($0,\dots,j_{0}-1$) will remain unaltered as $\alpha$ increases. On the contrary, the remaining slits $(j_{0},\dots,D-1)$ will have their transmission coefficients changed as a function of $\alpha$. The measurements started, for each value of $j_{0}$, at $\alpha=0$, which corresponds to $N$ maximally distinguishable states in a $D$-dimensional space. It ended at $\alpha=1$, where linearly dependent states are being discriminated. For $0<\alpha<1$, the linear dependence will only rely on the relationship between $N$ and $D$. As $\alpha$ increases, the optimal minimum error probability decreases and, therefore, distinguishability decreases as well.

Appendix B Detectors and data analysis

Detectors were CMOS cameras Thorlabs DCC1545M, whose pixels are $5.2\,{\rm\mu m}$ wide squares. The camera gamma function was programmed to obtain, in terms of gray levels, a linear response with respect to the incoming optical power. The laser power and camera gain were adjusted to avoid saturation in any pixel. For every generated state, three images from the Fourier plane were captured and summed in order to obtain less noisy results. The exposure time of each image was set the closest possible to 1/60 s, which corresponds to the flicker period of the SLM. Once the data has been taken, an integration over the $y$ axis is performed in order to trace over this non relevant direction since the intensities of interest were distributed along the $x$ axis. In this way, only one-dimensional intensity distributions are recorded. The background noise was subtracted from each pattern taking the minimum of the diffraction envelope as a reference.

Once the one-dimensional patterns were processed, state discrimination is performed. This is done by registering the intensity at $N$ transverse positions given by $x_{k}=-\lambda fm_{k}/dN$, where $f$ is the focal length of the lens, $\lambda$ is the optical field wavelength, $N$ is the number of states being discriminated in each set, and

[TABLE]

is a number that controls the distribution of the positions $x_{k}$ around $x=0$. We shall define the intensities $I_{k}^{(j)}$ as the intensity measured at $x_{k}$ when the input state is $|\psi_{j}\rangle$. Thus, the experimentally estimated probability of having a photodetection at the $k$th detector is given by

[TABLE]

where $I_{kj}=I_{k}^{(j)}/\eta_{k}$ is the intensity registered at the $k$th transverse position after being compensated by the optical diffraction efficiency, given by $\eta_{k}={\rm sinc}^{2}(kax_{k}/f)$. We can note that the optical efficiency compensation depends solely on the detector position, independently of the state being discriminated.

For each input state, the probability of having a correct retrodiction of the state is given by $P(\delta_{j}|\psi_{j})$. The ME strategy aims to maximize the overall probability of correct retrodiction for each set, i.e.,

[TABLE]

where the factor $1/N$ corresponds to the a priori probability of each state $|\psi_{j}\rangle$ being selected, which are assumed here to be equally likely. Figures 8 to 12 show some of the results for the probabilities of correct retrodiction $P(\delta_{j}|\psi_{j})$ compared with their respective optimal expected value. The average $P_{\rm corr}$ for these cases has been shown in Fig. 2 of the main text.

B.1 Number of sets and states tested

Table 1 summarizes the dimensions ($D$), number of states per set ($N$), number of sets per $N$, and the number of sets and states per $D$ tested in our experiment. A total number of 1851 sets were tested where more than 13900 states were generated and discriminated. The process of generation, characterization and discrimination of these states was completely automated through a single portable computer controlling the CMOS cameras and the SLM simultaneously. Before the experiment, a calibration of the CMOS was made in order to set the correct positions for state discrimination and state estimation. Afterwards, the full discrimination procedure was performed: For each state, a 1080$\times$1920 grayscale image containing a phase mask, like the examples $\Phi_{\rm PA}(x,y)$ and $\Phi_{\rm GD}(x,y)$ shown in Fig. 5, is addressed to the SLM in order to prepare the state. Then, both cameras registered the intensity distribution from Fourier and image planes. State estimation is performed according to the methods described in Secs. A.2 and A.3. Discrimination procedure is described in the main text and in Sec. B.