Experimental study of quantum thermodynamics using optical vortices

TL;DR

This paper presents an experimental platform using optical vortices to study non-equilibrium quantum thermodynamics, verifying fluctuation relations and demonstrating Maxwell's demon operation in high-dimensional quantum systems.

Contribution

It introduces a novel optical setup employing Laguerre-Gaussian beams to emulate quantum harmonic oscillators and experimentally investigates quantum thermodynamic properties.

Findings

Verified the Jarzynski fluctuation relation experimentally.

Reconstructed work distribution using two-point measurement scheme.

Demonstrated Maxwell's demon operation with optical vortices.

Abstract

Non-equilibrium thermodynamics and quantum information theory are interrelated research fields witnessing an increasing interest, both theoretical and experimental. This is manly due to the broadness of these theories, which found applications in many different fields of science, ranging from biology to the foundations of physics. Here, by employing the orbital angular momentum of light, we propose a new platform for studying non-equilibrium properties of high dimensional quantum systems. Specifically, we use Laguerre-Gaussian beams to emulate the energy eigenstates of a two-dimension quantum harmonic oscillator having angular momentum. These light beams are subjected to a process realized by a spatial light modulator and the corresponding work distribution is experimentally reconstructed employing a two-point measurement scheme. The Jarzynski fluctuation relation is then verified. We also demonstrate the operation of the system as a Maxwell's demon.

Equations

$$⟨e^{-\beta W}⟩\equiv \int dWP(W)e^{-\beta W}=e^{-\beta \Delta F},$$

$$H$$

$$L_z$$

$$\text{Energy: }{\epsilon }_{\ell p}=(|\ell |+2p+1)\hslash \omega$$

$$\text{Angular momentum: }{\lambda }_{\ell }=\hslash \ell$$

$${\epsilon }_{\ell }=(|\ell |+1)\hslash \omega$$

$$P(W)=\sum _{\ell ,{\ell }^{\prime }}{p}_{\ell {\ell }^{\prime }}\delta \left(W-W_{\ell {\ell }^{\prime }}\right),$$

$${\rho }_{\beta }=\sum _{\ell =-\infty }^{+\infty }{p}_{\ell }|\ell ⟩⟨\ell |,\text{with}{p}_{\ell }=\frac{\text{mbox}{e}^{-\beta {\epsilon }_{\ell }}}{Z}\text{and}Z=e^{\beta \hslash \omega }\tanh \frac{\beta \hslash \omega }{2}.$$

$$p(|\ell |)=N{\text{e}}^{-\beta {\epsilon }_{|\ell |}}=N{\text{e}}^{-\beta (|\ell |+1)\hslash \omega },$$

$${\mathcal{L}}_{\pm 5}|\ell ⟩=|\ell \pm 5⟩.$$

$$⟨e^{-\sigma -I}⟩=1,$$

$$\mathbf{Y}=\mathbf{AX}.$$

$$\underset{\mathbf{A}}{\min }||\mathbf{Y}-\mathbf{AX}|{|}^2$$

Full text

I Introduction

The orbital angular momentum (OAM) of light is a property of the topology of the optical modes, and are characterized by discrete numbers associated to the amount of orbital angular momentum per photon in the mode Padgett . The natural family of optical modes with orbital angular momentum are the Laguerre-Gaussian (LG) modes, a set of solutions of the paraxial wave equation salehteich that are described by their radial number $p$ and the azimuthal number $\ell$. The study and application of these modes is relatively recent and has increased considerably in the last two decades Padgett ; molina ; Babiker .

Single photons populating modes with OAM are physical realizations of high-dimensional quantum states Vaziri ; Agnew ; Dada ; lofller ; Giovannini ; Krenn1 ; Krenn , leading to the possibility of encoding more than one bit of information per photon. Such photonic qudits can be explored in order to improve quantum communication schemes and quantum information processing Molina ; Bourennane ; Huber ; Groblacher ; Walborn ; Mafu ; Mirhosseini ; Mair . Moreover, the transverse amplitude profiles of LG light modes are formally identical to the energy eigenstates of the two-dimension quantum harmonic oscillator. Therefore, they stand as a platform for the emulation of these quantum systems in a variety of interesting problems. In the present work, we employ these light modes to experimentally study some thermodymical aspects of a high dimensional quantum system. A similar approach has been successfully used to study the quantum limits of a chaotic harmonic oscillator GBL . Non-equilibrium thermodynamics is fundamentally concerned to the characterization of the response of a system under external perturbations. The theory of the linear-response regime was developed in Refs. Callen ; Green ; Kubo , based on earlier works such as Refs. Einstein ; Johnson ; Nyquist . The information about the complete nonlinear response is contained in the so called fluctuation theorems, which have been proved for classical Bochkov ; Bochkov1 ; Jarzynski and for quantum systems Tasaki ; Hanggi .

Fluctuation relations can be understood as a quantification of the probability of observing a violation of the second law of thermodynamics for small systems (when fluctuations come into play) and short time-scales. Considering the new trend in miniaturization, such fluctuations and time-scales are becoming more important for the development of new technological devices Hanggi_rev . Therefore, the theoretical and experimental study of quantum fluctuation relations are of primary interest, both for fundamental issues and for understanding the limitations of implementing quantum information processing and communication devices.

The quantum versions of the classical fluctuation theorems are possible only due to the two-point measurement approach for defining work. Work performed on (or by) the system is defined as the difference between two energy measurements, one before and one after the considered process takes place. To be specific, let us consider an externally driven system $\mathcal{S}$, whose time-dependent Hamiltonian is denoted by $H_{\mathcal{S}}(t)$, initially in the thermal state $\rho_{\beta}$, with $\beta=1/k_{B}T$, where $T$ is the temperature of the system and $k_{B}$ is the Boltzmann constant. The scenario considered here can be divided into three steps: $i)$ projective measurement on the initial Hamiltonian, $H_{\mathcal{S}}(0)$, eigenbasis; $ii)$ unitary (driven) evolution for a time interval $\tau$; $iii)$ projective measurement on the final Hamiltonian, $H_{\mathcal{S}}(\tau)$, eigenbasis. Defining the two-point measurement variable $W_{mn}=\varepsilon_{m}-\varepsilon_{n}$, where $\varepsilon_{m}$ and $\varepsilon_{n}$ are the eigenvalues of $H_{\mathcal{S}}(\tau)$ and $H_{\mathcal{S}}(0)$, respectively, it is not difficult to show that this stochastic variable must obey the general fluctuation relation known as the Jarzynski equality Esposito_RMP_81_1665 ; Campisi_RMP_83_771

[TABLE]

where $P(W)$ is the probability density distribution associated with the random variable $W$ and $\Delta F=F_{\tau}-F_{0}$ is the variation of the free energy over the time interval $\tau$ during which the system is subjected to the process. It has been shown that Eq. (1) is also valid for unital processes Zyczkowsky , i.e. quantum maps that do not change the identity.

Note that the final state is not necessarily a thermal state, since it is generated by a projective measurement followed by an evolution. However, what appears in the right-hand side of Eq. (1) is the equilibrium quantity $F_{\tau}$, which concerns the state the system ends up in if it is allowed to thermalize with a reservoir at the same temperature as the initial one. By defining the entropy production as $\sigma=\beta\left(W-\Delta F\right)$ we can rewrite the fluctuation relation as $\langle\operatorname{e}^{-\sigma}\rangle=1$.

The experimental investigation of such relation is new, specially in quantum systems. Regarding classical systems we can mention the experiments reported in Refs. Hummer ; Liphardt ; Collin ; Blickle ; Harris ; Saira . In the quantum regime, experiments tend to get trickier or more complicated due to the difficulty in performing energy projective measurements on arbitrary systems. Only recently, based on an alternative scheme that avoid such measurements Dorner ; Mazzola , an experimental reconstruction of the work distribution associated with a process performed on a spin-$1/2$ system was reported Batalhao . Considering the projective measurements, the only experiment to date, as far as we know, was reported in Ref. An , where the authors employed trapped ions in order to investigate the work statistics associated with a harmonic oscillator. Here we contribute to this line by employing the projective measurement scheme to reconstruct the probability distribution associated with a process performed on the two-dimensional harmonic oscillator with angular momentum using an optical setup, thus providing a new experimental platform for the investigation of thermodynamic processes in the quantum regime.

II Simulating a quantum system with classical light

It is possible to simulate a class of quantum systems using classical light and the analogy between the paraxial wave equation and the two-dimension Schrödinger equation. This analogy has been explored experimentally to investigate the quantum limit of a chaotic harmonic oscillator GBL and to propose a study, similar to the one done here, in which the characteristic function of the work distribution could be measured Talarico . OAM optical modes emulate the energy eigenstates of the harmonic oscillator in the sense that the transverse distribution of the electric field of LG beams has the same form as the energy eigenfunctions of the 2-D quantum harmonic oscillator. Moreover, under appropriate conditions, the propagation of the light beams is equivalent to the Hamiltonian evolution of the harmonic oscillator GBL ; Talarico ; marcuse .

This type of simulation accounts for all oscillatory aspects of quantum systems, such as state superposition, coherence and decoherence. The intrinsic quantum properties of light itself do not come into play in this scenario, since we are exclusively interested in light’s modal structure, rather than in its photonic content, which is usually explored by using detectors like avalanche photodiodes (for single photons) or low-reverse-bias photodiodes (for the continuous variables regime).

In the scheme we present here, we use the OAM modes to represent the wave functions of the two-dimension quantum harmonic oscillator, for which the Hamiltonian and angular momentum operators $H$ and $L_{z}$ form a Complete Set of Commuting Observables. They are written in terms of the number operators for right ($N_{r}$) and left ($N_{l}$) circular quanta as

[TABLE]

and their eigenvalues are

[TABLE]

where $\ell$ and $p$ are the azimuthal and radial quantum numbers, respectively, which are analogous to the azimuthal and radial indices used to identify the elements of the LG basis of modes. Note that, for the subset of states having the quantum number $p=0$, i.e. whenever either $N_{r}$ or $N_{l}$ has eigenvalue zero, the state energy

[TABLE]

depends only on the azimuthal quantum number $\ell$. Thus, in this case, projections in the OAM basis are equivalent to projections in the energy eigenbasis.

If we restrict ourselves to processes acting on the harmonic oscillator that only change $\ell$ and we project the system’s final state onto an eigenstate of the angular momentum, the work done in one experiment run depends solely on the change in $|\ell|$. Indeed, when the system goes from an initial $\ell$ to a final $\ell^{\prime}$, the work done is $W_{\ell\ell^{\prime}}=(|\ell^{\prime}|-|\ell|)\hbar\omega$. The work probability distribution is then

[TABLE]

where $p_{\ell\ell^{\prime}}=p_{\ell}\ p_{\ell^{\prime}|\ell}$ is the probability of observing the transition $\ell\rightarrow\ell^{\prime}$, with $p_{\ell}$ being the probability of having $\ell$ at the input and $p_{\ell^{\prime}|\ell}$ the probability of observing $\ell^{\prime}$ at the output given that the input is $\ell$.

In the context of Jarzynski equality, $p_{\ell}$ is found in the expression of the initial thermal state $\rho_{\beta}=\operatorname{e}^{-\beta H}/Z$, where $Z$ is the partition function. This state may be explicitly written as:

[TABLE]

Note that, since the states $|\ell\rangle$ and $|-\ell\rangle$ have the same energy, their probabilities are the same: $p_{\ell}=p_{-\ell}$. In fact, every energy level has degeneracy 2, except for the ground state $|\ell=0\rangle$.

III Experimental setup

The sketch of the experimental setup is shown in Fig. 1. The light from a He-Ne laser is sent through a beam expander consisting of two lenses in a confocal arrangement, with focal lengths $f_{1}=50$ mm and $f_{2}=300$ mm, resulting in an expansion factor of 6. The expanded beam is sent to the first Spatial Light Modulator (SLM1), where an OAM mode is prepared with the usual approach with a forked hologram Padgett .

SLM1 generates modes with OAM $\hbar\ell$ per photon and they are sent to a second spatial light modulator, SLM2, where another phase mask realizes some operation on them depending on the protocol. The resulting light beam is sent to a device called mode sorter modesorter , which litterally sorts the different OAM components of the beam along the horizontal axis of a CCD camera. The measurement scheme is calibrated sending OAM modes one by one and measuring the intensity distribution at the output of the mode sorter with no phase modulation on SLM2 (flat phase mask). Typical calibration curves are shown in gray in Fig. 2, where each curve represents the intensity distribution at the output for a given value of $\ell$, from -15 to +15 in this case.

Step i of the two-point measurement protocol consists in preparing the thermal state described by Eq. (8) and performing a projective measurement in the initial Hamiltonian eigenbasis or, equivalently, in the OAM basis.

The preparation of the thermal state is made by sending a Gaussian mode with $\ell=0$ to the spatial light modulator SLM1, which applies masks that generate OAM states with $\ell$ ranging from $-7\leq\ell\leq 7$. Each mask is turned on for 3 s, according to a random choice of $\ell$ with weight $p_{\ell}$. The resulting light beam is sent to SLM2, that acts just as mirror in this case, and then to the mode sorter, which analyses the OAM components. The light intensity at the output of the mode sorter is measured with a CCD camera, and the images are analyzed as explained in detail in Appendix A. A typical result is shown in Fig. 3, where the final distribution is obtained from 300 runs of the experiment. The distribution obtained for the absolute value of OAM is normalized and fitted to the function

[TABLE]

which represents the Boltzmann Distribution, leaving $N$ (the normalization factor) and $\beta\hbar\omega$ as free parameters. For the results shown in Fig. 3, we obtained an excellent agreement with the fitted function, with $\beta\hbar\omega$ = 0.67 $\pm$ 0.01.

The parameter $\beta\hbar\omega=\hbar\omega/k_{B}T$ can be interpreted as the ratio between the ground state energy and the typical scale of thermal energy at temperature $T$. So, for a given system, the greater $\beta\hbar\omega$, the lower the temperature.

We show how to prepare a thermal state and perform the projective measurements in the energy eigenbasis in order to illustrate this procedure. However, in the second step of the protocol, we prepare each energy eigenstate separately and submit it to the process in order to obtain the transition probabilities. This is strictly equivalent to using the states resulting from the first measurement and means a considerable simplification in the set up.

Step ii consists in sending input modes prepared with SLM1 and having OAM ranging from $\ell=-7$ to $\ell=+7$ to SLM2, where a phase mask is applied, realizing the process whose work distribution will be measured. This mask couples the input mode to other OAM modes, thus inducing OAM transitions. The energy spectrum of the system is discrete and infinite. However, the thermal weight of states corresponding to higher energies can be made negligible by choosing sufficiently low temperatures, so that we can truncate the initial distribution of states, as we did.

Using the mode sorter, we are able to measure the final distribution of OAM modes and their corresponding weights. This device implements a projective measurement in the orbital angular momentum basis, which in our case is equivalent to the energy under the assumption that the radial number $p=0$. Observing the mode sorter output, we can compute the transition probabilities and, consequently, reconstruct the work distribution. Nonetheless, as the calibration figure shows, there is a considerable overlap between adjacent curves (adjacent orbital angular momenta). This appreciably reduces the resolution of the OAM sorter. Newer generations of mode sorters, as well as other strategies resolve , minimise this technical inconvenient.

In the present proof of principle experiment, we overcome this issue using a process in SLM2 which generates superpositions of OAM modes that can be easily resolved by the mode sorter. Specifically, the process in our experiment implements the linear operation $(\mathcal{L}_{+5}+\mathcal{L}_{-5})/\sqrt{2}$, where we define

[TABLE]

In this way, the overlap between the two components at the output becomes negligible. Typical measurement results are shown in Figs. 2(a) and 2(b).

For each measured output, we performed a linear least squares regression in order to obtain the values of the orbital angular momenta and their respective weights (see Appendix A for details). The normalised set of all OAM weights for all outputs is exhibited in the matrix of Fig. 4(a). This is a density plot where the index of the input (output) modes are the labels of the vertical (horizontal) axis. In other words, these are the transition probabilities for a typical run of the experiment. Fig. 4(b) shows the matrix for the ideal process.

IV Results and Discussion

In our experiment, we measured the conditional transition probabilities $p_{\ell^{\prime}|\ell}$ as shown in Fig. 4. As explained earlier, we truncated the OAM space and limited the input modes to $|\ell|\leq 7$. That is to say we operate in a regime of low temperatures where the Boltzmann weights for $|\ell|>7$ can be neglected, i.e. $\beta\hbar\omega\gtrsim 1$, i.e. $\hbar\omega\gtrsim k_{B}T$.

Fig. 5 shows a plot of the quantity $\left\langle\operatorname{e}^{-\beta W}\right\rangle$ as a function of $\beta\hbar\omega$. The inset displays the probability distributions of work computed from the measurement results for $\beta\hbar\omega=2$. The probabilities in the vertical axis are obtained summing up all values of $p_{\ell\ell^{\prime}}$ (given in Eq.7) for which the corresponding transition results in a given value of work $W$. Regarding the Jarzynski equality, Eq. (1), the considered process gives $\Delta F=0$ as it does not change the Hamiltonian of the system. In other words, transitions are induced, but the energy levels do not change. Therefore, Eq. (1) becomes simply $\left\langle\operatorname{e}^{-\beta W}\right\rangle=1$.

The curve named theory represents $\left\langle\operatorname{e}^{-\beta W}\right\rangle$ computed for an ideal (but truncated) process. The corresponding ideal transition probabilities are illustrated in Fig. 4(b). The curve named experiment represents the same quantity computed from the measured transition probabilities shown in Fig. 4(a). Within the range $\beta\hbar\omega\lesssim 1$ (gray area), we obtain $\left\langle\operatorname{e}^{-\beta W}\right\rangle<1$ even for the theory curve, due to the truncation of the input states in $|\ell|\leq 7$. However, for larger values of $\beta\hbar\omega$, we see that the theory curve is essentially constant and equal to 1, while the experiment curve is always below unity. We interpret this result as a consequence of entropy increase due to natural technical limitations present in a real world measurement (which includes classical fluctuations coming from laser pointing instability, mechanical vibrations on the set-up and camera thermal noise) and experimental imperfections (such as misalignment, limited pixel resolution on the SLM and on the camera and limited optical resolution on the mode sorter). The uncertainty band shown in Fig. 5 accounts only for the fluctuations detected upon several subsequent identical measurements (see Appendix B for details). For instance, for $\beta\hbar\omega=2$, we have found $\operatorname{e}^{-\beta W}=0.910\pm 0.046$ (95%-confidence interval), clearly different from 1. We atribute the remaining difference to the experimental imperfections listed above that are not captured by our error estimation procedure, but that is captured by the Jarzynski’s fluctuation relation.

V Maxwell’s demon

While Eq. (1) is valid for any unital process Zyczkowsky , in a general context including measurements and feedback, i.e., when Maxwell’s demon come into play, a new equality holds Ueda :

[TABLE]

where $I$ is the information acquired by the demon due to the measurement process. This implies a modification of the second law of thermodynamics as $\sigma\geq-I$, highlighting the demon’s main feature, which is the use of information to reduce entropy production.

We propose an experimental scheme for realizing Maxwell’s demon using OAM modes. Fig. 7 shows the sketch of the suggested scheme. A laser beam is sent to spatial light modulator SLM1, which produces a thermal state of OAM modes as discussed in Sec. III. The light beam prepared in the thermal state is then sent through mode sorter MS1. The modes with positive angular momentum will be deflected to the right, while those with negative angular momentum will drift to the left. These two groups of beams are separated and directed to mode sorters MS2 and MS3 working in reverse reverse1 ; reverse2 , converting them back to OAM modes. With this approach, it is possible to separate OAM modes according to sign of $\ell$. SLM2 is used to apply the operation $\mathcal{L}_{+5}$ to the modes with negative $\ell$ and SLM3 is used to apply $\mathcal{L}_{-5}$ to the modes with positive $\ell$. They are finally sent to mode sorters MS4 and MS5 that perform the final measurement of the OAM using CCD cameras CCD1 and CCD2. The measurement and feedback, in this case, increases the probability of lowering the absolute value of orbital angular momentum ($|\ell|$) of the system, thus extracting work from the initial thermal state. As an example, let us start with $\ell=+3$: the operation $(\mathcal{L}_{+5}+\mathcal{L}_{-5})/\sqrt{2}$, and subsequent measurement, results in either $\ell=-2$ or $\ell=+8$ with equal probabilities, leading to an average work of $2\hbar\omega$. Now, if Maxwell’s demon takes action, $\ell$ goes invariably from +3 to -2 and $W=-\hbar\omega<0$, which means that work is extracted from the system.

As the measurement of OAM sign provides Maxwell’s demon with one bit of information, $I=\ln 2$, and Eq. (11) leads to $\left\langle\operatorname{e}^{-\sigma}\right\rangle=2$ for all $\beta\hbar\omega$, which is the Jarzynski’s fluctuation relation with demon’s action for the ideal case. This scenario has been computed and plotted in Fig. 6 (curve named theory). Note that the values of $\left\langle\operatorname{e}^{-\sigma}\right\rangle$ below 2 at high temperatures (shaded gray area on the left) are simply due to truncation of the OAM space dimension.

We would also like to have some insight on the effect of noise in the Maxwell’s demon scenario. In order to do that, we have added the same noise that appears in Fig. 4 to the probability transitions of the Maxwell’s demon scheme. The result is the curve named noise in Fig. 6, that carries an uncertainty band calculated in the same way as for the experimental curve in Fig. 5. Notice that the random noise decreases the average value of $\operatorname{e}^{-\sigma}$, meaning that the decrease in entropy production caused by Maxwell’s demon would be affected by experimental noise.

The inset shows the probability distributions for the possible values of work in the case of $\beta\hbar\omega=$ 2. Even though the value of $\left\langle\operatorname{e}^{-\sigma}\right\rangle$ changes dramatically from 1 to 2 when the demon takes action, the work distribution and the average work do not change too much: $\langle W\rangle=5.0\,\hbar\omega$ without Maxwell’s demon and $\langle W\rangle=4.8\,\hbar\omega$ with it. This is due to the fact that, in this sytem, for low enough temperatures such as $\beta\hbar\omega=2$, there is a large contribution from the state component with $\ell=0$ in the input thermal state. For an input with $\ell=0$, both $\mathcal{L}_{+5}$ and $\mathcal{L}_{-5}$ contribute for positive work, which dominates the work distribution.

VI Conclusion

In conclusion, we have experimentally investigated the quantum version of thermodynamic work and the Jarzynski’s fluctuation relation using the orbital angular momentum of light, a discrete degree of freedom with infinite dimension usually employed in the single-photon regime to realize a qudit, with applications in quantum communication and quantum information processing. Here, by exploring the analogy between the Paraxial wave equation and Schrödinger equation, we have used OAM of light to simulate the eigenstates of the two-dimensional quantum harmonic oscillator and their evolution through a given process. We measured the work distribution associated with this process and obtained the experimental Jarzynski’s fluctuation relation. We have also proposed a scheme for implementing a Maxwell’s demon measurement and feedback action. These results illustrate the usefulness of Laguerre-Gaussian beams as a practical platform to investigate aspects of the growing field of quantum thermodynamics in high dimensional Hilbert spaces. Given the versatility of this platform, one can consider using it, for example, in the study of the role of multipartite entanglement in thermodynamic processes, as well as the role of the environment, i.e., non-unitary and non-unital processes.

VII Acknowledgments

We acknowledge financial support from the Brazilian funding agencies CNPq, CAPES, and the National Institute for Quantum Information (INCT-IQ). M.P.J.L. acknowledges support from EPSRC awards EP/N032853/1. M.J.P. acknowledges support from EPSRC QuantIC (EP/M01326X/1) and ERC TWISTS (Grant No. 192382).

VIII Appendix

VIII.1 Data fitting

When an Laguerre-Gaussian mode passes through a mode sorter, its intensity profile (initially presented in a donut-like shape) becomes an elongated spot along the vertical direction on the camera. The integration of such an image along the vertical axis gives a curve (among the 31 calibration curves on Figs. 2(a) and 2(b)) showing the marginal distribution of intensity along the horizontal direction of the sorted LG mode. The right side of Fig. 1 outlines this idea.

The horizontal size of an image on the camera ($\sim$ 80 pixels) sets the extension of the horizontal axis in Fig. 2, which shows collections of lists of 80 elements. Let us call the lists corresponding to the LG inputs $\mathbf{x}_{j}=\{x_{j,k}\}_{0\leq k<80}$ and the lists at the output $\mathbf{y}_{i}=\{y_{i,k}\}_{k}$, where $-7\leq i\leq 7$ and $-15\leq j\leq 15$ may be associated with $\ell$ and $\ell^{\prime}$, respectively. We model the process implemented by SLM2 as a linear operation represented by a $15\times 31$ matrix $\mathbf{A}$ acting on the set of possible input modes $\mathbf{X}=\{\mathbf{x}_{i}\}_{i}$ and leading to the set of its outcomes $\mathbf{Y}=\{\mathbf{y}_{j}\}_{j}$, i.e.

[TABLE]

In order to find the matrix that best fits our experimental data, we apply a linear least squares approach, numerically solving the minimization problem

[TABLE]

with the additional constraints that all elements of matrix $\mathbf{A}$ must be non-negative and that the sum of the elements in each line must equal 1. The result is a matrix similar to the ones shown in Fig. 3.

The non-negativity is equivalent to the assumption that the overall phases of each OAM component of the output are the same or, at least, that the OAM components are far enough from each other so that they do not interfere and their relative phases do not play a significant role. The sum over each line equaling 1 stands for the unitarity of the process (no optical loss), which can be assumed without loss of generality.

VIII.2 Measurement uncertainty

A measurement consists in taking a picture, integrating it and obtaining its marginal intensity distribution $\mathbf{y}_{i}=\{y_{i,k}\}_{k}$. The goal here is to assess the uncertainty on each $y_{i,k}$ measured and, from this, to calculate the uncertainty:

• •

on each element of the matrix of conditional transition probabilities of Fig. 3(a) and

• •

on the mean value $\langle\operatorname{e}^{-\beta W}\rangle$ as a function of $\beta\hbar\omega$ (uncertainty band on Fig. 4).

Let us call $\sigma_{i,k}$ the uncertainty on $y_{i,k}$. In order to assess these $\sigma$’s, we performed a series of ten identical measurements on the same transformed mode over a time window of a few minutes and obtained a set of intensity distributions fluctuating, for each $k$, around a mean value $\mu_{k}$ with standard a deviation $\sigma_{k}$. These standard deviations turn out to be dependent on $\mu_{k}$ and on $k$ itself, but mostly on $\mu_{k}$. We noticed, for instance, that the relative standard deviation ($\sigma/\mu$) is always smaller than 10%, for any $k$. These measurements were used to model the typical error associated to a $\mathbf{y}_{i}$ measurement. This procedure led to a set of numbers $\sigma_{i,k}$ used as input for our model, in which we assume each $\mathbf{y}_{i}$ measured is a realization of a random variable following the multivariate normal distribution with estimated mean values $y_{i,k}$ and standard deviations $\sigma_{i,k}$. This procedure above allows us to simulate sets of measurements, realizing Monte Carlo experiments.

Ten different experimental matrices $\mathbf{Y}$ were randomly generated in this manner, from each of which we numerically solved the minimization problem above to find a different probability matrix $\mathbf{A}$. We could see from the set os matrices $\mathbf{A}$ that the relative uncertainty on each matrix element was never bigger than 2%.

Similarly (and finally), we performed 1000 Monte Carlo experiments in order to estimate the uncertainty on $\langle\operatorname{e}^{-\beta W}\rangle$ for each $\beta\hbar\omega$ ranging from 0.05 to 5. We observed that the random variable $\operatorname{e}^{-\beta W}$ nearly follows a normal distribution for all values of $\beta\hbar\omega$. For instance, for $\beta\hbar\omega=2$, we have found $\langle\operatorname{e}^{-\beta W}\rangle=0.910$ with a standard deviation $\sigma=0.022$. From the $1.96\sigma$ rule, we established our 95%-confidence interval for $\operatorname{e}^{-2W/\hbar\omega}$ to be $0.910\pm 0.046$. By doing that same for all values of $\beta\hbar\omega$, we were able to plot the uncertainty band shown in Fig. 4.