Quantum channel correction with twisted light using compressive sensing

TL;DR

This paper demonstrates that compressive sensing can effectively estimate and correct atmospheric turbulence effects on high-dimensional quantum channels using classical light, enhancing free-space optical communication.

Contribution

It introduces a novel application of compressive sensing for quantum channel correction with twisted light, including a numerical simulation demonstrating its effectiveness.

Findings

High fidelity in corrected states

Reduced atmospheric turbulence effects

Quantified non-separability post-correction

Abstract

Compressive sensing is used to perform high-dimensional quantum channel estimation with classical light. As an example, we perform a numerical simulation for the case of a three-dimensional classically non-separable state that is propagated through atmospheric turbulence. Using singular value thresholding algorithm based compressive sensing, we determine the channel matrix, which we subsequently use to correct for the atmospheric turbulence induced distortions. As a measure of the success of the procedure, we calculate the fidelity and the trace distance of the corrected density matrix with respect to the input state, and compare the results with those of the density matrix for the uncorrected state. Furthermore, we quantify the amount of classical non-separability in the density matrix of the corrected state by calculating its negativity. The results show that compressive sensing could contribute in the development and implementation of free-space quantum and optical communication systems.

Equations

$${\rho }_{\Lambda }=({\Lambda }_A\otimes \text{mathds}{1}_B)(|\psi ⟩⟨\psi |),$$

$$|\psi ⟩=\frac{1}{\sqrt{N}}\sum _{n=1}^N|n{⟩}_A\otimes |n{⟩}_B,$$

$$|\psi ⟩=\frac{1}{\sqrt{2M+1}}\sum _{\ell =-M}^M|\ell ,0{⟩}_A|{\lambda }_{\ell }{⟩}_B,$$

$$|{\ell,p}\rangle\rightarrow\sum_{{\ell}^{\prime},p^{\prime}}|{\ell^{\prime},p^{\prime}}\rangle\Lambda_{{\ell},p}^{{\ell}^{\prime},p^{\prime}},\$$

$$|{n}\rangle\rightarrow\sum_{m}|{m}\rangle\Lambda_{n}^{m}.\$$

$$\Lambda =$$

$$=$$

$$\rho =\sum _{i=1}^{N^2}{\alpha }_i{\tau }_i,$$

$${\alpha }_i=\text{tr}\left\{{\tau }_i\rho \right\}.$$

$$\displaystyle\begin{split}\text{tr}\{\rho_{\text{r}}\}&=1,\\ \rho_{\text{r}}&=\rho_{\text{r}}^{\dagger},\\ \text{tr}\{\rho_{\text{r}}\tau_{i}\}&=\alpha_{i}~{}~{}\text{for}~{}~{}i=1...m,\end{split}\$$

$$||{\rho }_{\text{r}}|{|}_{\text{tr}}=\text{tr}\left\{\sqrt{{\rho }_{\text{r}}{\rho }_{\text{r}}^{†}}\right\}.$$

$${\rho }_{\text{r}}=\sum _j|{\varphi }_j⟩{\sigma }_j⟨{\varphi }_j|,$$

$$M\bar{\rho }=\bar{\alpha },$$

$$M{\bar{\rho }}_{\text{r}}\ne \bar{\alpha }.$$

$$|{\bar{v}}_i|={\alpha }_i-\hat{v}\cdot {\bar{\rho }}_{i-1},$$

$${\bar{\rho }}_i={\bar{\rho }}_{i-1}+{\bar{v}}_i.$$

$$|\Psi {⟩}_{\text{in}}=$$

$$=$$

$$|\Psi ⟩=\frac{1}{\sqrt{3}}\left[|\ell ,0⟩|{\lambda }_{\ell }⟩+|0,0⟩|{\lambda }_0⟩+|\bar{\ell },0⟩|{\lambda }_{\bar{\ell }}⟩\right].$$

$${\partial }_{z}f(\mathbf{r})=\frac{i}{2k_0}{\nabla }_{\perp }f(\mathbf{r})-i{k}_0\delta n(\mathbf{r})f\left(\mathbf{r}\right),$$

$$\theta (\mathbf{R})={\left(2\pi k_0^2\Delta z\right)}^{1/2}{\mathcal{F}}^{-1}\left\{\chi (\mathbf{K}){\left[\frac{{\Phi }_n(\mathbf{K})}{{\Delta }_k^2}\right]}^{1/2}\right\},$$

$$⟨\chi ({\mathbf{K}}_1){\chi }^{\ast }({\mathbf{K}}_2)⟩={\left(2\pi {\Delta }_k\right)}^2\delta ({\mathbf{K}}_1-{\mathbf{K}}_2),$$

$${\theta }_{SH}(j\Delta x,l\Delta y)=$$

$$\times \exp \left[2\pi i\left(\frac{jp}{3^n{N}_x}+\frac{lq}{3^n{N}_y}\right)\right].$$

$$⟨a^2(p,q,n)⟩=$$

$$=$$

$$F({\rho }_{\text{c}},|\Psi ⟩⟨\Psi |)=\text{tr}\left\{\sqrt{\sqrt{{\rho }_{\text{c}}}|\Psi ⟩⟨\Psi |\sqrt{{\rho }_{\text{c}}}}\right\},$$

$$r_0=0.185{\left(\frac{{\lambda }^2}{C_n^{2}z}\right)}^{3/5}.$$

$$D({\rho }_{\text{c}},|\Psi ⟩⟨\Psi |)=\frac{1}{2}\text{tr}\left\{\left|{\rho }_{\text{c}}-|\Psi ⟩⟨\Psi |\right|\right\},$$

$$\mathcal{E}({\rho }_{\text{c}})=\frac{1}{2}\sum _n(|{\lambda }_n|-{\lambda }_n),$$

Full text

I Introduction

Quantum entanglement in higher dimensions is a property that can enhance communication security using quantum key distribution Gibson et al. (2004); Barreiro et al. (2008); Mirhosseini et al. (2015). The spatial degrees of freedom of single-photon fields allow infinitely many modes, such as the orbital angular momentum (OAM) modes Allen et al. (1992), which can be used to design and prepare such higher dimensional quantum systems. Optical fields of OAM modes have an azimuthal phase dependence given by $\exp(i\ell\phi)$, with $\phi$ the azimuthal angle. They carry quantized OAM of $\ell\hbar$ per photon. Photons that have the above described phase dependence and carry OAM are commonly referred to as twisted photons.

Free-space optical and quantum communication systems based on the OAM states of light, are adversely affected by turbulence in the atmosphere, which leads to a distortion of the OAM states Paterson (2005); Andrews and Phillips (1998); Mabena and Roux (2019). This distortion further results in the decay of quantum entanglement, which has been studied both theoretically Smith and Raymer (2006); Gopaul and Andrews (2007); Jha et al. (2010); Roux (2011); Leonhard et al. (2015); Brünner and Roux (2013) and experimentally Pors et al. (2011); Malik et al. (2012); Hamadou Ibrahim et al. (2013); Goyal et al. (2016); Zhang et al. (2016). Our understanding of the behavior of OAM-based quantum systems indicates that the potential offered by the OAM states of light for free-space communication requires efficient methods to mitigate against the adverse effects of atmospheric turbulence.

The standard method of characterizing a quantum channel is called standard quantum process tomography (SQPT) Nielsen and Chuang (2000). The characterization and determination of quantum processes are crucial tasks necessary for the implementation of quantum communication and information processing systems, enabling the design of mitigation strategies against noise and other distortions that may negatively affect a quantum system. However, SQPT is a resource intensive process that scales dramatically with the size of the quantum system Nielsen and Chuang (2000); Mohseni et al. (2008). The number of resources required for SQPT is generally $O(N^{4})$, where $N$ is the dimension of the Hilbert space.

Scintillation, caused by atmospheric turbulence, leads to the distortion of input OAM modes, distributing the power into other modes Paterson (2005); Tyler and Boyd (2009), thus increasing the dimensionality of the photonic quantum systems when expanded in terms of OAM modes Anguita et al. (2008). The increase in dimensionality caused by scintillation implies that the application of SQPT would not be favorable for characterizing a turbulent atmospheric channel. The mechanism that produces scintillation is a pure phase modulation caused by the fluctuating refractive index of the atmosphere. It does not produce a loss prior to the measurement stage. Although the atmosphere does have additional loss mechanisms such as absorption and scattering, we do not include them in our current investigation, for reasons explained below.

For successful implementation of quantum communication systems with twisted photons, a different method of determining and characterizing the turbulent quantum channel is needed. This is where we employ compressive sensing Donoho (2006); Candès et al. (2006); Candès (2006) — a data processing method that provides an efficient mechanism for the recovery of unknown signals from only a fraction of the required measurements. A generalization of this method to matrices is called matrix completion Candès and Recht (2009). In traditional compressive sensing the signal is required to have certain properties such as sparseness in the appropriate basis. In the matrix generalization case, the matrices must also obey certain criteria. For instance, the method is more likely to be successful for low rank matrices. Fortunately, this is the case for a single realization of atmospheric turbulence — a snapshot of the atmospheric conditions along the entire path of the channel as experienced by a single light pulse. Scintillation is treated as a unitary process that preserves the purity of quantum states.

Gross et al. Gross et al. (2010) have established a method based on a random selection of Pauli measurements for efficient reconstruction of an unknown quantum state. They showed that their method can reconstruct a rank $r$ unknown density matrix with only $O(rN\log^{2}N)$ measurements, in contrast to the $O(N^{2})$ measurements for the standard method. In another study Tonolini et al. (2014), a high-dimensional entangled state was reconstructed from a significantly smaller number of measurements, using a related approach based on compressive sensing. While this method for the reconstruction of signals from an underdetermined system of equations is very popular in signal and image processing applications, it has now attracted interest and become topical in quantum information science related applications as well Howland et al. (2011); Howland and Howell (2013); Howland et al. (2013); Kosut et al. (2008); Shabani et al. (2011); Gross et al. (2010).

The aim of this work is to show with the aid of numerical simulations that compressive sensing is a viable technique for channel estimation in the specific case of the transmission of twisted light through atmospheric turbulence. Moreover, inspired by a recent study on the characterization of a quantum channel with classical light Ndagano et al. (2017); Mabena and Roux (2017), we developed a more efficient protocol for characterizing turbulent channels by combining it with compressive sensing. In turn, it leads to an efficient scheme for the correction of twisted photons after passing through atmospheric turbulence.

The proposed scheme does not suffer any detrimental effects due to loss, because it is based on classical light (a bright coherent state). In the context of quantum states, a probability distribution of different losses would cause the quantum state to become mixed. However, the proposed scheme allows one to make the pulse shorter than the time scale given by the Greenwood frequencyGreenwood (1977) of the turbulent medium. As a result, the accumulated loss that the pulse experience does not represent a probability distribution. The pulse only experiences a single constant accumulated loss. Hence, no mixing takes place and the state remains pure. For this reason, we do not include any loss mechanisms in our current investigation.

The compressive sensing part requires relatively few random measurements, while the use of classical light removes the intrinsic limitation of quantum mechanics that requires repeated measurements on an ensemble. Also, since classical light is extremely bright compared to a few discrete photons, we can perform the different random measurements at the same time, thus leading to a significant speed up in time. The measurements on the classical light upon propagating through atmospheric turbulence will be performed using sum-frequency generation as previously discussed in Ref. Mabena and Roux (2017). Measurements on the output state are made by optically combining the output state with a tailored measurement state through a non-linear crystal, the measurement state selects out a specific component in the output state. This entire process leads to a photon detection which represents a measurement. The compressive sensing model used in this work is based on the singular value thresholding algorithm Cai et al. (2010) that has been modified in a way similar to the one used in Ref. Tonolini et al. (2014).

The outline of this contribution is as follows. In Sec.II we present the model for this work. The numerical simulation method is discussed briefly in Sec. III. Section IV is based on the results of this work. Finally, the conclusions are given in Sec. V.

II Model

II.1 Channel matrix

The Choi-Jamiolkowski isomorphism Jiang et al. (2013) establishes a correspondence between a completely positive trace-preserving quantum map $\Lambda$ and a quantum state $\rho$ by

[TABLE]

where

[TABLE]

is a maximally entangled state, $\mathds{1}$ is the identity operator for subsystem B, and $N$ is the dimension of the Hilbert space. This isomorphism means that the identification and characterization of the quantum channel is tantamount to performing a quantum state tomography.

Here, the partites of the state in Eq. (2) are represented by the spatial modes and the wavelengths of a classical optical field. A perfect correlation between these degrees of freedom (maximal nonseparability) gives us an exact analogy with a maximally entangled quantum state. So, the typical input state for our consideration reads Mabena and Roux (2017)

[TABLE]

where $|{\ell,0}\rangle$ represents a Laguerre-Gaussian (LG) mode with azimuthal index $\ell$ and radial index $p=0$; $\lambda_{\ell}$ is the wavelength of the corresponding LG mode; and $M$ is an integer representing the maximum OAM. In general, one can have an arbitrary radial index, so that the LG mode would be $|{\ell,p}\rangle$. However, for the moment, we assume correlation between the wavelength and the azimuthal index and therefore we set $p=0$.

Transmission through the atmosphere causes the OAM modes to scatter into other modes. The atmospheric turbulence only affects the spatial degree of freedom and leaves the wavelength unaffected. Therefore, after propagating through the atmosphere, a given input OAM mode becomes

[TABLE]

where $|{\ell^{\prime},p^{\prime}}\rangle$ is an LG mode with azimuthal index $\ell^{\prime}$ and radial index $p^{\prime}$; $\Lambda_{{\ell},p}^{{\ell}^{\prime},p^{\prime}}$ is the tensor representation of the atmospheric turbulence Kraus operator. It is this tensor that represents the effect of the atmosphere. We need to characterize it to mitigate against its effect.

To simplify notation, we’ll index the input and output modes by single integers. Hence, the scattering process is represented by

[TABLE]

where the Kraus operation is now represented by an $N\times M$ matrix. The Krauss operator can thus be represented, either in terms of the tensor or in terms of a matrix

[TABLE]

In the most general case, the matrix $\Lambda_{n}^{m}$ is rectangular $M\neq N$. The input state is usually defined in terms of a finite number of modes. So, the dimensionality of the input Hilbert space is finite and determines the number of columns $N$ of the matrix. On the other hand, the number of rows $M$ is determined by the crosstalk induced by the atmospheric turbulence. Since there are an infinite number of spatial modes and since the scintillation process can potentially scatter the input modes into any combination of output modes, one can expect that $M\gg N$ and that $M\rightarrow\infty$. However, the scattering is not uniform — the dominant scattering tends to produce modes lying close to the input modes. One can therefore truncate the output space to a finite number of dimensions, depending on the accuracy that is required.

II.2 Compressive sensing

Here we briefly review the compressing sensing procedure that we used for our work. For this purpose, we follow the procedure of Tonolini et al. Tonolini et al. (2014).

Using the matrix completion method Candès and Recht (2009), one can recover a low-rank matrix when some of the elements of the matrix are unknown. Instead of recovering the density matrix from a sample of its elements $\rho_{i,j}$, we recover the full matrix that represents the states in terms of a sample of the results of measurements made on the system.

In the most general case, the density matrix can be decomposed in terms of the Bloch representation

[TABLE]

where $\alpha_{i}$ are the elements of the Bloch vector, $N$ represents the dimensionality of the state vector and $\tau_{i}$ denotes the generalized Gell-Mann matrices (GGMs), including the identity matrix. The GGMs reduce to the Pauli matrices in the case where $N=2$. In this formalism, the GGMs form a convenient measurement basis for the characterization of the state. To determine the state through a full tomography, one must perform measurements that reveal $N^{2}$ real parameters $\alpha_{i}$. The Bloch vector elements are the expectation values given by the trace

[TABLE]

In the implementation of the modified matrix completion problem Gross et al. (2010), we consider an under-sampled set of measurements ($m\ll N^{2}$) chosen at random — i.e., we consider a situation where there is only a subset of the total possible measurements $\alpha_{i}$. The optimization problem is thus described as follows: Minimize $||\rho_{\text{r}}||_{\text{tr}}$ such that

[TABLE]

where $\rho_{\text{r}}$ is the to-be-recovered density matrix, and $||\cdot||_{\text{tr}}$ is the trace norm of the matrix, given by

[TABLE]

The compressive sensing algorithm is based on the singular value thresholding (SVT) algorithm Cai et al. (2010). However, we modify the algorithm to take advantage of known properties of the state that we intend to recover Tonolini et al. (2014). Essentially, in applying the SVT algorithm we perform an eigenvalue decomposition of the density matrix

[TABLE]

where $|{\phi_{j}}\rangle$ are the eigenvectors of $\rho_{\text{r}}$ and $\sigma_{j}$ are the corresponding eigenvalues. Given the above decomposition, we apply the thresholding operator on the eigenvalues, by selecting the eigenvalues above a certain threshold $\epsilon_{0}$ that we fixed a priori. Ensuring that the eigenvalues of the density matrix are real, we force the density matrix to be Hermitian. The normalization of the density matrix is achieved by dividing the resultant matrix by its trace.

The algorithm uses a guess matrix as a starting point for the optimization process. The most crucial part in setting up the guess matrix is choosing its dimensions, because the dimensions of the matrix must be representative of the atmospheric turbulence conditions. The matrix must be big enough to accommodate most of the non-zero elements in the OAM modal spectrum. Once we have determined the dimensions of the guess matrix, we can apply the algorithm.

After performing the thresholding on the density matrix and normalizing it, we obtain a density matrix that may represent a real physical system, but it no longer corresponds to the density matrix whose measurements gave the correct measurement results $\{\alpha_{i}\}$ with respect to the GGMs $\tau_{i}$ for $i=1...m$ Tonolini et al. (2014). To understand the reason for this issue, we use a geometrical perspective, represented in terms of hyperplanes. For this purpose, we re-express Eq. (8) in vector notation:

[TABLE]

where $\bar{\rho}$ is the vectorized density matrix, $M$ is a matrix with rows representing the vectorized GGM, and $\bar{\alpha}$ is the vector of the measurement results. In this representation, each vectorized GGM $\bar{\tau}_{i}$ and its corresponding measurement result ${\alpha}_{i}$ can be associated with a hyperplane in an $N^{2}$ dimensional space. As such, the compressive sensing approach attempts to solve an under-determined linear system of equations. The solution should be a single point common to all the $m$ hyperplanes. However, the procedure explained above does not produce a point that lies in the linear space of Eq. (12):

[TABLE]

To resolve this issue, we project the density matrix back into the linear space, determined by ${\tau}_{i}$ and their corresponding results $\alpha_{i}$. In the process, we modify the density matrix $\rho_{\text{r}}$ to give the correct measurement results. The projection is done stepwise for each value of $\alpha_{i}$. The projection starts by defining a vector $\bar{v}_{i}$ normal to the hyperplane $\bar{\tau}_{i}\bar{\rho}=\alpha_{i}$, and with a magnitude

[TABLE]

where $\hat{v}$ is the unit vector that is normal to the hyperplane and $\bar{\rho}_{i-1}$ is the density matric obtained in the previous iteration. The new density matrix becomes

[TABLE]

The projection process is performed for all $m$ measurements. At the end of this process, we start again and perform the procedure for Hermiticity and trace normalization. In other words, we perform the thresholding operation and further recompose the matrix according to the conditions in the optimization problem statement. All these steps are performed iteratively until the norm of the difference between the density matrices of two consecutive iterations lies within a predefined tolerance.

II.3 Channel correction

Here follows the main contribution of our work. It involves the generation of the Kraus operator matrix from the estimated density matrix, obtained from the compressive sense procedure described above.

The full Kraus operator matrix can be obtained from the density matrix of the output state. This step is made possible by the fact that, despite the randomness of the medium, propagation through a turbulent atmosphere is a unitary process. Therefore, the output quantum state after transmission through a single realization of atmospheric turbulence is always a pure state, provided that the input state was a pure state. We also assume that the truncation of the output space to a finite dimension does not affect the unitarity of the process significantly, provided that we use a large enough number for the output dimension. As a result, we consider the reconstructed Kraus operator as a unitary operator, so that its Hermitian adjoint represents the inverse process.

After the Kraus operator $\Lambda$ has been reconstructed, using the compressive sensing methods described above, we perform the correction process as follows

[TABLE]

The validity of the unitary assumption is assessed by the quality of the correction, as discussed below.

III Numerical simulation and computations

We performed numerical simulations to test the proposed compressive sensing-based scheme for quantum channel estimation and correction. For this purpose, we consider a three-dimensional bipartite input state (qutrit). The two degrees of freedom that represents the two partites are the spatial mode (OAM mode) and the wavelength. The input state can be expressed as

[TABLE]

In the simulation, the input state is represented by three $n\times n$ sampled functions for the three spatial modes, where $n=1024$. Note that Eq. (17) represents a classical field that is nonseparable (as opposed to being entangled), expressed in terms of Dirac notation.

Each of the modes are separately propagated in the simulation process. The propagation of paraxial optical fields in atmospheric turbulence is described by the stochastic parabolic equation

[TABLE]

where $\nabla_{\perp}=\partial^{2}_{x}+\partial^{2}_{y}$, $\delta n(\mathbf{r})$ is the refractive index fluctuation of the atmosphere ($n=1+\delta n$), $k_{0}=2\pi/\lambda$ is the wave number, and $\lambda$ is the wavelength of the optical field.

In weak scintillation, the only effect of the atmospheric turbulence is a phase perturbation on the optical field. It means that propagation through the atmosphere under weak scintillation conditions can be represented by two steps. The first step is a random phase modulation of the input optical field, which represents the perturbation of the field and leads to refraction. The subsequent step is free-space propagation (without turbulence) over the full propagation distance.

For arbitrary scintillation conditions, one can still use these two steps to simulate propagation through turbulence. However, one would repeat the two steps multiple times, each time propagating over a short enough distance to ensure weak scintillation conditions for that step. During each step, the optical field is modulated by a different random phase screen. Therefore, to simulate the propagation of an optical field through a turbulent atmosphere we use a split-step method Anguita et al. (2008), in which these two steps are repeated several times.

The generation of the random phase screen entails transforming a 2D-array of random complex numbers that have zero mean and unit variance into an array that has the same statistics as the atmospheric turbulence. This process is also known as filtering Gaussian noise and is given by Martin and Flatté (1988); Knepp (1983)

[TABLE]

where $\Delta z$ is the partitioned propagation distance between two consecutive phase screens, $\mathcal{F}^{-1}\{\cdot\}$ denotes the inverse Fourier transform, $\Phi_{n}(\mathbf{K})$ is the Kolmogorov power spectral density (PSD) for the refractive index, $\Delta_{k}$ is the grid-spacing in the spatial frequency domain, and $\mathbf{K}=(k_{x},k_{y})$ is the transverse wavevector. The normally distributed complex random function $\chi(\mathbf{K})$ has zero mean and is $\delta$-correlated, such that

[TABLE]

where the angled brackets $\langle\cdot\rangle$ denote an ensemble average. The random phase function generated by Eq. (19) is complex [unless $\chi^{*}(\mathbf{K})=\chi(-\mathbf{K})$]. Hence, a single application of Eq. (19) produces two independent random phase screens — the real and imaginary parts of the complex phase function.

It is important to point out that the Fourier calculation in Eq. (19) does not take into account the effect of large eddies, which are excluded due to the discrete grid samples in the Fourier domain. As a result, the statistics obtained from such phase screens do not represent the Kolmogorov structure function correctly. One way to improve the accuracy is to add sub-harmonics Lane et al. (1992) to the phase function generated by the FFT method in Eq. (19). The sub-harmonic phase function is given by

[TABLE]

The variance of the randomly generated functions $a$ and $b$ is given as,

[TABLE]

where, $\Delta{p_{n}}=\Delta p/3^{n}$, $\Delta{q_{n}}=\Delta q/3^{n}$, $N_{x}$ and $N_{y}$ are the numbers of points in the $x$ and $y$ directions, respectively, and $N_{s}$ is the number of sub-harmonics.

In the simulation, we propagate the input state through the turbulent medium by performing the split-step process on each of the three input modes, with their associated wavelengths. The two-step process is iterated several times to perform a multi-phase screen propagation of the input state through a realization of the turbulent medium. The complete propagation is then done for several realizations.

The parameters that are used for the propagation process are: input beam waist radius $w_{0}=0.1$ m; propagation distance $z=2z_{R}$, where $z_{R}=\pi w_{0}^{2}/\lambda$ is the Rayleigh range; and turbulence strength (refractive index structure constant) $C_{n}^{2}=1\times 10^{-16}~{}\text{m}^{-2/3}$. The wavelengths of the three modes are different, however, in practical setups one can always make these differences to be very small. For the purpose of the simulations, the following wavelengths were used for the different modes: $\lambda_{1}=1.0~{}\mu\text{m}$, $\lambda_{2}=1.020~{}\mu\text{m}$, and $\lambda_{3}=1.040~{}\mu\text{m}$.

The turbulent medium causes crosstalk, transferring power to numerous higher order LG modes. Therefore, the dimensions of the Kraus operator could be very large. The strength of the crosstalk and the dimensionality of the output state depends on various parameters, including the strength of turbulence and the distance of propagation. These parameters combine to determine the strength of the scintillation — the extent of the distortion imparted by the medium on the state.

The dimension of the Hilbert space of the input state is $N_{\text{in}}=N_{A}\cdot N_{B}$. For the state described in Eq. (17), we have $N_{A}=N_{B}=3$, which implies that $N_{\text{in}}=9$. The corresponding density matrix has $81$ elements and is depicted graphically in Fig. 1. A unitary process does not change the minimum number of modes required to represent the state. However, the nature of these modes is unknown due to the lack of a priori knowledge of the unitary process associated with a single realization of the turbulent medium. As a result, the output state must be represented in terms of some nominal modal basis.

Here, we use the LG modes as our chosen output modal basis. While the input state has only three-dimensions per degree of freedom, the dimension of the output state in terms of the LG modes is much larger. The turbulent medium affects the dimension of the one subsystem $N_{A}\rightarrow N_{A}^{\prime}$ but leaves the that the other $N_{B}$ the same. The dimension of the output density matrix thus becomes $N_{\text{out}}\times N_{\text{out}}$, where $N_{\text{out}}=N_{A}^{\prime}\cdot N_{B}$.

In our simulation, we include radial indices up to $p=7$. For each value of the radial indices, we considered azimuthal indices in the range $-14\leq\ell\leq 15$. It implies that the output dimension in subsystem $A$ becomes $N_{A}^{\prime}=210$, so that the total output dimension becomes $N_{\text{out}}=630$. This number is required to represent the output beam profiles with adequate fidelity for the parameters that we used in the simulation. The output density matrix thus has almost $400~{}000$ elements. To specify a density matrix obtained in this way, one needs as many measurements.

However, the state remains pure for each realization of the atmospheric turbulence, because the process is unitary. As a result, the Kraus operator matrix associated with that realization has a low rank. The purity of a state implies that it is represented by a state vector $|{\psi_{\text{out}}}\rangle$. The density matrix of such a pure state has rank equal to unity, which allows us to use compressive sensing for estimating the output state and thereby determining the Kraus operator matrix. Using compressive sensing, we can determine the state of the output density matrix to a high level of accuracy, using a significantly smaller number of measurements.

For this work, we use approximately $5\%$ of the total number of measurements (which in this case is roughly 20000) to reconstruct the output density matrix reliably. Using the reconstructed density matrix, we extract the elements of the Kraus operator for the atmospheric turbulence. To test the reliability of the reconstruction we use the Kraus operator to generate the correction matrix and apply it to the output density matrix. The result is then tested for fidelity and trace distance against the input qutrit state Eq. (17).

Upon application of the compressive sensing algorithm described in Sec. II.2, a density matrix of the output state is obtained. Since a single realization of atmospheric turbulence is a unitary process, the output state is pure. Furthermore, given that the output density matrix is the outer product of the state vector and its adjoint, it follows that a single element of the density matrix is a product of two elements. Hence, by choosing one column or row of the output density matrix and dividing it by the relevant element, the elements of the state vector can be obtained, apart from an overall phase constant. These elements are then rearranged to form the Kraus operator $\Lambda^{m}_{n}$, which is subsequently used to correct the scintillation as described in Sec. II.3.

IV Results

The numerical simulations of the propagation and correction of OAM states in atmospheric turbulence was performed several times for different realizations of a turbulent medium. The performance of the compressive sensing-based correction scheme is assessed by calculating the fidelity and the trace distance from these results.

The fidelity with respect to the initial maximally nonseparable state in Eq. (17) is given by

[TABLE]

where $|{\Psi}\rangle$ is the input state, and $\rho_{\text{c}}$ is the corrected density matrix from our compressive sensing algorithm. The results of these calculations are shown in Fig. 2 for two different scintillation strengths: $W=0.5$ in Fig. 2(a) and $W=2$ in Fig. 2(b). Here, the scintillation strength is represented by the dimensionless number $W=w_{0}/r_{0}$, where $r_{0}$ is the Fried parameter Fried (1966), which is given by

[TABLE]

The trace distance, which is shown in Fig. 3, for the same two scintillation strengths, is defined as

[TABLE]

where the magnitude $|A|$ of a matrix is given as $|A|=\sqrt{A^{\dagger}A}$. The data in Figs. 2 and 3 are plotted against the different realizations of the atmospheric turbulence.

In Fig. 2, it is observed that the corrected density matrices are close to the ideal maximally nonseparable input state, in contrast to the uncorrected, truncated density matrices. The mean fidelity for the corrected density matrices over the 100 turbulence realizations in the first case with $W=0.5$ is $0.942\pm 0.003$, while the uncorrected density matrices give a mean fidelity of $0.38\pm 0.02$. Here, the error is given as the standard error. For the second case, $W=2$, the mean fidelity of the corrected density matries over 100 turbulence realizations is $0.925\pm 0.005$ and that of the uncorrected density matrices is $0.265\pm 0.012$.

A similar conclusion about the success of the correction method can be reached by looking at Fig. 3, which shows the trace distance. The mean trace distance for the corrected density matrices for $W=0.5$ is $0.31\pm 0.02$, while that of the uncorrected density matrices is $0.892\pm 0.011$. For the second scintillation strength, $W=2$, the mean trace distance for the corrected density matrix is $0.359\pm 0.011$ and for the uncorrected density matrices it is $0.955\pm 0.004$.

An example of the real parts of the elements of an uncorrected and corrected density matrix is shown in Fig. 4 for an arbitrary turbulence realization. Comparison of the unperturbed density matrix in Fig. 1 with the density matrices in Fig. 4 shows the potential advantage of the correction process.

To investigate how the entanglement (or classical non-separability) is affected by the compressive sensing correction technique, we calculate the negativity of the density matrices. The negativity of a state is given by

[TABLE]

where $\lambda_{n}$ are the eigenvalues of the partially transposed density matrix. The partial transpose of a density matrix is obtained by performing a transpose on the density matrix of one subsystem, leaving the other the same.

Figure 5 shows that the negativity of the state improves with application of the correction procedure that is based on the compressive sensing technique. The mean negativity for the corrected density matrices over 50 realizations is $0.83\pm 0.04$. This value can be contrasted with that of the uncorrected, truncated density matrices, which is $0.56\pm 0.18$, as evidence of the improvement.

For the practical design of compressive sensing-based channel correction for optical and quantum communication systems, we considered the mean negativity of the corrected density matrices as a function of the output dimension $N_{\text{out}}$. In Fig. 6, we display three curves for different scintillation strengths $W$ (by changing the propagation distance).

One can see that the mean negativity of the corrected density matrices increases with $N_{\text{out}}$, indicating the effect of the chosen size of the output Hilbert space. Furthermore, we observe that the mean negativity becomes saturated. The saturation level and the rate of saturation depends on the scintillation strength.

These observations suggest a possible future improvement to the proposed scheme. By applying methods such as deep learning Goodfellow et al. (2016); Nielsen (2015) or other generic machine learning algorithms Alpaydin (2004); Ismail et al. (2019), one may be able to train the system to determine an optimal value for $N_{\text{out}}$, given certain turbulence conditions.

V Conclusions

Using a numerical analysis, we demonstrate the performance of a proposed compressive sensing-based channel correction method. A classically non-separable state, consisting of three OAM modes with different wavelengths, was used as input to a numerically simulated turbulent free-space channel. Using compressive sensing-based state tomography, we reconstructed the output state and used it to determine the Kraus operator matrix for the channel. The singular value thresholding technique was used for the compressive sensing algorithm. The results show that compressive sensing drastically reduces the number of measurements required for the characterization of the turbulent channel. Although the channel estimation process uses classical light, it determines the Kraus operator matrix for the channel and thus allows it to be used for quantum communication. Consequently, the proposed scheme would be useful in the design of quantum communication systems that are based on the OAM states of light.

As a further study, we intend to investigate the use of deep learning or machine learning as a method to inform the system on the appropriate dimension of the Kraus operator matrix. It could be done by training a model to determine an appropriate value of this dimension, based on the parameters of the turbulent free-space channel.

Acknowledgements

CMM acknowledges support from the CSIR National Laser Centre. The research for this work was supported in part by the National Research Foundation (NRF) of South Africa (Grant Number: 118532).