Quantum channel correction with twisted light using compressive sensing
aa r X i v : . [ qu a n t - ph ] D ec Quantum channel correction with twisted light using compressive sensing
Chemist M. Mabena
1, 2 and Filippus S. Roux
3, 2, ∗ CSIR National Laser Centre, P.O. Box 395, Pretoria 0001, South Africa School of Physics, University of the Witwatersrand, Johannesburg 2000, South Africa National Metrology Institute of South Africa, Meiring Naud´e Road, Brummeria, Pretoria, South Africa
Compressive sensing is used to perform high-dimensional quantum channel estimation with clas-sical light. As an example, we perform a numerical simulation for the case of a three-dimensionalclassically non-separable state that is propagated through atmospheric turbulence. Using singularvalue thresholding algorithm based compressive sensing, we determine the channel matrix, whichwe subsequently use to correct for the atmospheric turbulence induced distortions. As a measure ofthe success of the procedure, we calculate the fidelity and the trace distance of the corrected densitymatrix with respect to the input state, and compare the results with those of the density matrixfor the uncorrected state. Furthermore, we quantify the amount of classical non-separability in thedensity matrix of the corrected state by calculating its negativity. The results show that compressivesensing could contribute in the development and implementation of free-space quantum and opticalcommunication systems.
I. INTRODUCTION
Quantum entanglement in higher dimensions is a prop-erty that can enhance communication security usingquantum key distribution [1–3]. The spatial degreesof freedom of single-photon fields allow infinitely manymodes, such as the orbital angular momentum (OAM)modes [4], which can be used to design and prepare suchhigher dimensional quantum systems. Optical fields ofOAM modes have an azimuthal phase dependence givenby exp( iℓφ ), with φ the azimuthal angle. They carryquantized OAM of ℓ ~ per photon. Photons that have theabove described phase dependence and carry OAM arecommonly referred to as twisted photons .Free-space optical and quantum communication sys-tems based on the OAM states of light, are adverselyaffected by turbulence in the atmosphere, which leads toa distortion of the OAM states [5–7]. This distortion fur-ther results in the decay of quantum entanglement, whichhas been studied both theoretically [8–13] and experi-mentally [14–18]. Our understanding of the behavior ofOAM-based quantum systems indicates that the poten-tial offered by the OAM states of light for free-space com-munication requires efficient methods to mitigate againstthe adverse effects of atmospheric turbulence.The standard method of characterizing a quantumchannel is called standard quantum process tomography(SQPT) [19]. The characterization and determinationof quantum processes are crucial tasks necessary for theimplementation of quantum communication and informa-tion processing systems, enabling the design of mitigationstrategies against noise and other distortions that maynegatively affect a quantum system. However, SQPT isa resource intensive process that scales dramatically withthe size of the quantum system [19, 20]. The number ofresources required for SQPT is generally O ( N ), where ∗ [email protected] N is the dimension of the Hilbert space.Scintillation, caused by atmospheric turbulence, leadsto the distortion of input OAM modes, distributing thepower into other modes [5, 21], thus increasing the di-mensionality of the photonic quantum systems when ex-panded in terms of OAM modes [22]. The increase indimensionality caused by scintillation implies that theapplication of SQPT would not be favorable for charac-terizing a turbulent atmospheric channel. The mecha-nism that produces scintillation is a pure phase modu-lation caused by the fluctuating refractive index of theatmosphere. It does not produce a loss prior to the mea-surement stage. Although the atmosphere does have ad-ditional loss mechanisms such as absorption and scatter-ing, we do not include them in our current investigation,for reasons explained below.For successful implementation of quantum communica-tion systems with twisted photons, a different method ofdetermining and characterizing the turbulent quantumchannel is needed. This is where we employ compres-sive sensing [23–25] — a data processing method thatprovides an efficient mechanism for the recovery of un-known signals from only a fraction of the required mea-surements. A generalization of this method to matrices iscalled matrix completion [26]. In traditional compressivesensing the signal is required to have certain propertiessuch as sparseness in the appropriate basis. In the matrixgeneralization case, the matrices must also obey certaincriteria. For instance, the method is more likely to besuccessful for low rank matrices. Fortunately, this is thecase for a single realization of atmospheric turbulence —a snapshot of the atmospheric conditions along the en-tire path of the channel as experienced by a single lightpulse. Scintillation is treated as a unitary process thatpreserves the purity of quantum states.Gross et al. [27] have established a method based on arandom selection of Pauli measurements for efficient re-construction of an unknown quantum state. They showedthat their method can reconstruct a rank r unknowndensity matrix with only O ( rN log N ) measurements,in contrast to the O ( N ) measurements for the standardmethod. In another study [28], a high-dimensional entan-gled state was reconstructed from a significantly smallernumber of measurements, using a related approach basedon compressive sensing. While this method for the recon-struction of signals from an underdetermined system ofequations is very popular in signal and image processingapplications, it has now attracted interest and becometopical in quantum information science related applica-tions as well [27, 29–33].The aim of this work is to show with the aid of nu-merical simulations that compressive sensing is a viabletechnique for channel estimation in the specific case ofthe transmission of twisted light through atmosphericturbulence. Moreover, inspired by a recent study onthe characterization of a quantum channel with classi-cal light [34, 35], we developed a more efficient proto-col for characterizing turbulent channels by combining itwith compressive sensing. In turn, it leads to an efficientscheme for the correction of twisted photons after passingthrough atmospheric turbulence.The proposed scheme does not suffer any detrimentaleffects due to loss, because it is based on classical light (abright coherent state). In the context of quantum states,a probability distribution of different losses would causethe quantum state to become mixed. However, the pro-posed scheme allows one to make the pulse shorter thanthe time scale given by the Greenwood frequency[36] ofthe turbulent medium. As a result, the accumulated lossthat the pulse experience does not represent a proba-bility distribution. The pulse only experiences a singleconstant accumulated loss. Hence, no mixing takes placeand the state remains pure. For this reason, we do notinclude any loss mechanisms in our current investigation.The compressive sensing part requires relatively fewrandom measurements, while the use of classical light re-moves the intrinsic limitation of quantum mechanics thatrequires repeated measurements on an ensemble. Also,since classical light is extremely bright compared to afew discrete photons, we can perform the different ran-dom measurements at the same time, thus leading toa significant speed up in time. The measurements onthe classical light upon propagating through atmosphericturbulence will be performed using sum-frequency gener-ation as previously discussed in Ref. [35]. Measurementson the output state are made by optically combining theoutput state with a tailored measurement state througha non-linear crystal, the measurement state selects outa specific component in the output state. This entireprocess leads to a photon detection which represents ameasurement. The compressive sensing model used inthis work is based on the singular value thresholding al-gorithm [37] that has been modified in a way similar tothe one used in Ref. [28].The outline of this contribution is as follows. In Sec.IIwe present the model for this work. The numerical simu-lation method is discussed briefly in Sec. III. Section IV isbased on the results of this work. Finally, the conclusions are given in Sec. V. II. MODELA. Channel matrix
The Choi-Jamiolkowski isomorphism [38] establishesa correspondence between a completely positive trace-preserving quantum map Λ and a quantum state ρ by ρ Λ = (Λ A ⊗ B )( | ψ ih ψ | ) , (1)where | ψ i = 1 √ N N X n =1 | n i A ⊗ | n i B , (2)is a maximally entangled state, is the identity operatorfor subsystem B, and N is the dimension of the Hilbertspace. This isomorphism means that the identificationand characterization of the quantum channel is tanta-mount to performing a quantum state tomography.Here, the partites of the state in Eq. (2) are representedby the spatial modes and the wavelengths of a classicaloptical field. A perfect correlation between these degreesof freedom (maximal nonseparability) gives us an exactanalogy with a maximally entangled quantum state. So,the typical input state for our consideration reads [35] | ψ i = 1 √ M + 1 M X ℓ = − M | ℓ, i A | λ ℓ i B , (3)where | ℓ, i represents a Laguerre-Gaussian (LG) modewith azimuthal index ℓ and radial index p = 0; λ ℓ is thewavelength of the corresponding LG mode; and M is aninteger representing the maximum OAM. In general, onecan have an arbitrary radial index, so that the LG modewould be | ℓ, p i . However, for the moment, we assumecorrelation between the wavelength and the azimuthalindex and therefore we set p = 0.Transmission through the atmosphere causes the OAMmodes to scatter into other modes. The atmospheric tur-bulence only affects the spatial degree of freedom andleaves the wavelength unaffected. Therefore, after propa-gating through the atmosphere, a given input OAM modebecomes | ℓ, p i → X ℓ ′ ,p ′ | ℓ ′ , p ′ i Λ ℓ ′ ,p ′ ℓ,p , (4)where | ℓ ′ , p ′ i is an LG mode with azimuthal index ℓ ′ andradial index p ′ ; Λ ℓ ′ ,p ′ ℓ,p is the tensor representation of theatmospheric turbulence Kraus operator. It is this tensorthat represents the effect of the atmosphere. We need tocharacterize it to mitigate against its effect.To simplify notation, we’ll index the input and outputmodes by single integers. Hence, the scattering processis represented by | n i → X m | m i Λ mn . (5)where the Kraus operation is now represented by an N × M matrix. The Krauss operator can thus be represented,either in terms of the tensor or in terms of a matrixΛ = X ℓ,p,ℓ ′ ,p ′ | ℓ ′ , p ′ i Λ ℓ ′ ,p ′ ℓ,p h ℓ, p | = X m,n | m i Λ mn h n | . (6)In the most general case, the matrix Λ mn is rectangular M = N . The input state is usually defined in terms of afinite number of modes. So, the dimensionality of the in-put Hilbert space is finite and determines the number ofcolumns N of the matrix. On the other hand, the numberof rows M is determined by the crosstalk induced by theatmospheric turbulence. Since there are an infinite num-ber of spatial modes and since the scintillation processcan potentially scatter the input modes into any combi-nation of output modes, one can expect that M ≫ N and that M → ∞ . However, the scattering is not uni-form — the dominant scattering tends to produce modeslying close to the input modes. One can therefore trun-cate the output space to a finite number of dimensions,depending on the accuracy that is required. B. Compressive sensing
Here we briefly review the compressing sensing proce-dure that we used for our work. For this purpose, wefollow the procedure of Tonolini et al. [28].Using the matrix completion method [26], one can re-cover a low-rank matrix when some of the elements of thematrix are unknown. Instead of recovering the densitymatrix from a sample of its elements ρ i,j , we recover thefull matrix that represents the states in terms of a sampleof the results of measurements made on the system.In the most general case, the density matrix can bedecomposed in terms of the Bloch representation ρ = N X i =1 α i τ i , (7)where α i are the elements of the Bloch vector, N repre-sents the dimensionality of the state vector and τ i denotesthe generalized Gell-Mann matrices (GGMs), includingthe identity matrix. The GGMs reduce to the Pauli ma-trices in the case where N = 2. In this formalism, theGGMs form a convenient measurement basis for the char-acterization of the state. To determine the state througha full tomography, one must perform measurements thatreveal N real parameters α i . The Bloch vector elementsare the expectation values given by the trace α i = tr { τ i ρ } . (8) In the implementation of the modified matrix comple-tion problem [27], we consider an under-sampled set ofmeasurements ( m ≪ N ) chosen at random — i.e., weconsider a situation where there is only a subset of thetotal possible measurements α i . The optimization prob-lem is thus described as follows: Minimize || ρ r || tr suchthat tr { ρ r } = 1 ,ρ r = ρ † r , tr { ρ r τ i } = α i for i = 1 ...m, (9)where ρ r is the to-be-recovered density matrix, and || · || tr is the trace norm of the matrix, given by || ρ r || tr = tr (cid:26)q ρ r ρ † r (cid:27) . (10)The compressive sensing algorithm is based on the sin-gular value thresholding (SVT) algorithm [37]. However,we modify the algorithm to take advantage of knownproperties of the state that we intend to recover [28].Essentially, in applying the SVT algorithm we performan eigenvalue decomposition of the density matrix ρ r = X j | φ j i σ j h φ j | , (11)where | φ j i are the eigenvectors of ρ r and σ j are the cor-responding eigenvalues. Given the above decomposition,we apply the thresholding operator on the eigenvalues,by selecting the eigenvalues above a certain threshold ǫ that we fixed a priori. Ensuring that the eigenvalues ofthe density matrix are real, we force the density matrixto be Hermitian. The normalization of the density ma-trix is achieved by dividing the resultant matrix by itstrace.The algorithm uses a guess matrix as a starting pointfor the optimization process. The most crucial part insetting up the guess matrix is choosing its dimensions,because the dimensions of the matrix must be represen-tative of the atmospheric turbulence conditions. The ma-trix must be big enough to accommodate most of thenon-zero elements in the OAM modal spectrum. Oncewe have determined the dimensions of the guess matrix,we can apply the algorithm.After performing the thresholding on the density ma-trix and normalizing it, we obtain a density matrix thatmay represent a real physical system, but it no longercorresponds to the density matrix whose measurementsgave the correct measurement results { α i } with respectto the GGMs τ i for i = 1 ...m [28]. To understand thereason for this issue, we use a geometrical perspective,represented in terms of hyperplanes. For this purpose,we re-express Eq. (8) in vector notation: M ¯ ρ = ¯ α, (12)where ¯ ρ is the vectorized density matrix, M is a matrixwith rows representing the vectorized GGM, and ¯ α is thevector of the measurement results. In this representation,each vectorized GGM ¯ τ i and its corresponding measure-ment result α i can be associated with a hyperplane in an N dimensional space. As such, the compressive sensingapproach attempts to solve an under-determined linearsystem of equations. The solution should be a singlepoint common to all the m hyperplanes. However, theprocedure explained above does not produce a point thatlies in the linear space of Eq. (12): M ¯ ρ r = ¯ α. (13)To resolve this issue, we project the density matrixback into the linear space, determined by τ i and theircorresponding results α i . In the process, we modify thedensity matrix ρ r to give the correct measurement results.The projection is done stepwise for each value of α i . Theprojection starts by defining a vector ¯ v i normal to thehyperplane ¯ τ i ¯ ρ = α i , and with a magnitude | ¯ v i | = α i − ˆ v · ¯ ρ i − , (14)where ˆ v is the unit vector that is normal to the hyper-plane and ¯ ρ i − is the density matric obtained in the pre-vious iteration. The new density matrix becomes¯ ρ i = ¯ ρ i − + ¯ v i . (15)The projection process is performed for all m measure-ments. At the end of this process, we start again andperform the procedure for Hermiticity and trace normal-ization. In other words, we perform the thresholding op-eration and further recompose the matrix according tothe conditions in the optimization problem statement.All these steps are performed iteratively until the normof the difference between the density matrices of two con-secutive iterations lies within a predefined tolerance. C. Channel correction
Here follows the main contribution of our work. It in-volves the generation of the Kraus operator matrix fromthe estimated density matrix, obtained from the com-pressive sense procedure described above.The full Kraus operator matrix can be obtained fromthe density matrix of the output state. This step is madepossible by the fact that, despite the randomness of themedium, propagation through a turbulent atmosphere isa unitary process. Therefore, the output quantum stateafter transmission through a single realization of atmo-spheric turbulence is always a pure state, provided thatthe input state was a pure state. We also assume that thetruncation of the output space to a finite dimension doesnot affect the unitarity of the process significantly, pro-vided that we use a large enough number for the outputdimension. As a result, we consider the reconstructedKraus operator as a unitary operator, so that its Hermi-tian adjoint represents the inverse process. After the Kraus operator Λ has been reconstructed,using the compressive sensing methods described above,we perform the correction process as follows | Ψ i in = (Λ A ⊗ B ) † | Ψ i out = X m,n | m i Λ mn h n | ⊗ B ! † | Ψ i out . (16)The validity of the unitary assumption is assessed by thequality of the correction, as discussed below. III. NUMERICAL SIMULATION ANDCOMPUTATIONS
We performed numerical simulations to test the pro-posed compressive sensing-based scheme for quantumchannel estimation and correction. For this purpose,we consider a three-dimensional bipartite input state(qutrit). The two degrees of freedom that represents thetwo partites are the spatial mode (OAM mode) and thewavelength. The input state can be expressed as | Ψ i = 1 √ (cid:2) | ℓ, i| λ ℓ i + | , i| λ i + | ¯ ℓ, i| λ ¯ ℓ i (cid:3) . (17)In the simulation, the input state is represented by three n × n sampled functions for the three spatial modes,where n = 1024. Note that Eq. (17) represents a clas-sical field that is nonseparable (as opposed to being en-tangled), expressed in terms of Dirac notation.Each of the modes are separately propagated in thesimulation process. The propagation of paraxial opti-cal fields in atmospheric turbulence is described by thestochastic parabolic equation ∂ z f ( r ) = i k ∇ ⊥ f ( r ) − ik δn ( r ) f ( r ) , (18)where ∇ ⊥ = ∂ x + ∂ y , δn ( r ) is the refractive index fluc-tuation of the atmosphere ( n = 1 + δn ), k = 2 π/λ isthe wave number, and λ is the wavelength of the opticalfield.In weak scintillation, the only effect of the atmosphericturbulence is a phase perturbation on the optical field. Itmeans that propagation through the atmosphere underweak scintillation conditions can be represented by twosteps. The first step is a random phase modulation ofthe input optical field, which represents the perturbationof the field and leads to refraction. The subsequent stepis free-space propagation (without turbulence) over thefull propagation distance.For arbitrary scintillation conditions, one can still usethese two steps to simulate propagation through turbu-lence. However, one would repeat the two steps mul-tiple times, each time propagating over a short enoughdistance to ensure weak scintillation conditions for thatstep. During each step, the optical field is modulated bya different random phase screen. Therefore, to simulatethe propagation of an optical field through a turbulent at-mosphere we use a split-step method [22], in which thesetwo steps are repeated several times.The generation of the random phase screen entailstransforming a 2D-array of random complex numbersthat have zero mean and unit variance into an array thathas the same statistics as the atmospheric turbulence.This process is also known as filtering Gaussian noiseand is given by [39, 40] θ ( R ) = (cid:0) πk ∆ z (cid:1) / F − ( χ ( K ) (cid:20) Φ n ( K )∆ k (cid:21) / ) , (19)where ∆ z is the partitioned propagation distance be-tween two consecutive phase screens, F − {·} denotesthe inverse Fourier transform, Φ n ( K ) is the Kolmogorovpower spectral density (PSD) for the refractive index,∆ k is the grid-spacing in the spatial frequency domain,and K = ( k x , k y ) is the transverse wavevector. Thenormally distributed complex random function χ ( K ) haszero mean and is δ -correlated, such that h χ ( K ) χ ∗ ( K ) i = (2 π ∆ k ) δ ( K − K ) , (20)where the angled brackets h·i denote an ensemble aver-age. The random phase function generated by Eq. (19)is complex [unless χ ∗ ( K ) = χ ( − K )]. Hence, a singleapplication of Eq. (19) produces two independent ran-dom phase screens — the real and imaginary parts of thecomplex phase function.It is important to point out that the Fourier calcula-tion in Eq. (19) does not take into account the effect oflarge eddies, which are excluded due to the discrete gridsamples in the Fourier domain. As a result, the statisticsobtained from such phase screens do not represent theKolmogorov structure function correctly. One way toimprove the accuracy is to add sub-harmonics [41] to thephase function generated by the FFT method in Eq. (19).The sub-harmonic phase function is given by θ SH ( j ∆ x, l ∆ y ) = N s X n =1 1 X p,q = − [ a ( p, q, n ) + ib ( p, q, n )] × exp (cid:20) πi (cid:18) jp n N x + lq n N y (cid:19)(cid:21) . (21)The variance of the randomly generated functions a and b is given as, h a ( p, q, n ) i = h b ( p, q, n ) i =∆ p n ∆ q n Φ θ ( p ∆ p n , q ∆ q n ) , (22)where, ∆ p n = ∆ p/ n , ∆ q n = ∆ q/ n , N x and N y are thenumbers 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 statethrough 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 iteratedseveral times to perform a multi-phase screen propaga-tion of the input state through a realization of the tur-bulent medium. The complete propagation is then donefor several realizations.The parameters that are used for the propagation pro-cess are: input beam waist radius w = 0 . z = 2 z R , where z R = πw /λ is theRayleigh range; and turbulence strength (refractive in-dex structure constant) C n = 1 × − m − / . Thewavelengths of the three modes are different, however, inpractical setups one can always make these differences tobe very small. For the purpose of the simulations, thefollowing wavelengths were used for the different modes: λ = 1 . µ m, λ = 1 . µ m, and λ = 1 . µ m.The turbulent medium causes crosstalk, transferringpower 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 ofthe output state depends on various parameters, includ-ing the strength of turbulence and the distance of prop-agation. These parameters combine to determine thestrength of the scintillation — the extent of the distortionimparted by the medium on the state.The dimension of the Hilbert space of the input stateis N in = N A · N B . For the state described in Eq. (17),we have N A = N B = 3, which implies that N in = 9.The corresponding density matrix has 81 elements andis depicted graphically in Fig. 1. A unitary process doesnot change the minimum number of modes required torepresent the state. However, the nature of these modesis unknown due to the lack of a priori knowledge of theunitary process associated with a single realization of theturbulent medium. As a result, the output state must berepresented in terms of some nominal modal basis. FIG. 1. Graphical representation of the input state densitymatrix.
Here, we use the LG modes as our chosen output modalbasis. While the input state has only three-dimensionsper degree of freedom, the dimension of the output statein terms of the LG modes is much larger. The turbu-lent medium affects the dimension of the one subsystem N A → N ′ A but leaves the that the other N B the same.The dimension of the output density matrix thus becomes N out × N out , where N out = N ′ A · N B .In our simulation, we include radial indices up to p = 7.For each value of the radial indices, we considered az-imuthal indices in the range − ≤ ℓ ≤
15. It im-plies that the output dimension in subsystem A becomes N ′ A = 210, so that the total output dimension becomes N out = 630. This number is required to represent theoutput beam profiles with adequate fidelity for the pa-rameters that we used in the simulation. The outputdensity matrix thus has almost 400 000 elements. Tospecify a density matrix obtained in this way, one needsas many measurements. F i d e li t y W = 0.5
CorrectedUncorrected F i d e li t y W = 2
CorrectedUncorrected (a)(b)
FIG. 2. Fidelity of the corrected and the uncorrected den-sity matrices for the different atmospheric turbulence realiza-tions and for (a) W = 0 . W = 2. Triangle markersrepresent the fidelities for the uncorrected, truncated densitymatrices. Diamond markers represent the fidelities for thecorrected density matrices. However, the state remains pure for each realizationof the atmospheric turbulence, because the process isunitary. As a result, the Kraus operator matrix asso-ciated with that realization has a low rank. The purityof a state implies that it is represented by a state vector | ψ out i . The density matrix of such a pure state has rankequal to unity, which allows us to use compressive sensingfor estimating the output state and thereby determiningthe Kraus operator matrix. Using compressive sensing,we can determine the state of the output density matrixto a high level of accuracy, using a significantly smallernumber of measurements.For this work, we use approximately 5% of the totalnumber of measurements (which in this case is roughly20000) to reconstruct the output density matrix reliably. Using the reconstructed density matrix, we extract theelements of the Kraus operator for the atmospheric tur-bulence. To test the reliability of the reconstruction weuse the Kraus operator to generate the correction matrixand apply it to the output density matrix. The resultis then tested for fidelity and trace distance against theinput qutrit state Eq. (17). T r a c e d i s t a n c e W = 0.5
CorrectedUncorrected T r a c e d i s t a n c e W = 2
CorrectedUncorrected (a)(b)
FIG. 3. Trace distance of the corrected and the uncorrecteddensity matrices for the different atmospheric turbulence real-izations and for (a) W = 0 . W = 2. Triangle markersrepresent the uncorrected, truncated density matrices. Dia-mond markers represent the corrected density matrices. Upon application of the compressive sensing algorithmdescribed in Sec. II B, a density matrix of the output stateis obtained. Since a single realization of atmospheric tur-bulence is a unitary process, the output state is pure.Furthermore, given that the output density matrix is theouter product of the state vector and its adjoint, it followsthat a single element of the density matrix is a product oftwo elements. Hence, by choosing one column or row ofthe output density matrix and dividing it by the relevantelement, the elements of the state vector can be obtained,apart from an overall phase constant. These elements arethen rearranged to form the Kraus operator Λ mn , which issubsequently used to correct the scintillation as describedin Sec. II C. IV. RESULTS
The numerical simulations of the propagation and cor-rection of OAM states in atmospheric turbulence wasperformed several times for different realizations of aturbulent medium. The performance of the compressivesensing-based correction scheme is assessed by calculat-ing the fidelity and the trace distance from these results. (a)(b)
FIG. 4. Graphic representations of the (a) uncorrected, trun-cated density matrix, and the (b) corrected density matrix,after passing through the same simulated turbulence repre-sented by the scintillation strength W = 1. The fidelity with respect to the initial maximally non-separable state in Eq. (17) is given by F ( ρ c , | Ψ ih Ψ | ) = tr (cid:26)q √ ρ c | Ψ ih Ψ |√ ρ c (cid:27) , (23)where | Ψ i is the input state, and ρ c is the corrected den-sity matrix from our compressive sensing algorithm. Theresults of these calculations are shown in Fig. 2 for twodifferent scintillation strengths: W = 0 . W = 2 in Fig. 2(b). Here, the scintillation strengthis represented by the dimensionless number W = w /r ,where r is the Fried parameter [42], which is given by r = 0 . (cid:18) λ C n z (cid:19) / . (24)The trace distance, which is shown in Fig. 3, for thesame two scintillation strengths, is defined as D ( ρ c , | Ψ ih Ψ | ) = 12 tr {| ρ c − | Ψ ih Ψ ||} , (25) where the magnitude | A | of a matrix is given as | A | = √ A † A . The data in Figs. 2 and 3 are plotted against thedifferent realizations of the atmospheric turbulence.In Fig. 2, it is observed that the corrected density ma-trices are close to the ideal maximally nonseparable inputstate, in contrast to the uncorrected, truncated densitymatrices. The mean fidelity for the corrected density ma-trices over the 100 turbulence realizations in the first casewith W = 0 . . ± . . ± .
02. Here, theerror is given as the standard error. For the second case, W = 2, the mean fidelity of the corrected density matriesover 100 turbulence realizations is 0 . ± .
005 and thatof the uncorrected density matrices is 0 . ± . W = 0 . . ± .
02, while that of the uncorrected density matri-ces is 0 . ± . W = 2, the mean trace distance for the corrected densitymatrix is 0 . ± .
011 and for the uncorrected densitymatrices it is 0 . ± .
30 40 50Turbulence realization00.20.40.60.81 N e g a t i v i t y W = 0.5
CorrectedUncorrected
FIG. 5. The negativity of the correct and the uncorrected den-sity matrices for the different atmospheric turbulence realiza-tions. Triangle markers represent the uncorrected, truncateddensity matrices. Diamond markers represent the correcteddensity matrices.
To investigate how the entanglement (or classical non-separability) is affected by the compressive sensing cor-rection technique, we calculate the negativity of the den-sity matrices. The negativity of a state is given by E ( ρ c ) = 12 X n ( | λ n | − λ n ) , (26)where λ n are the eigenvalues of the partially transposeddensity matrix. The partial transpose of a density ma-trix is obtained by performing a transpose on the densitymatrix of one subsystem, leaving the other the same.Figure 5 shows that the negativity of the state im-proves with application of the correction procedure thatis based on the compressive sensing technique. The meannegativity for the corrected density matrices over 50 re-alizations is 0 . ± .
04. This value can be contrastedwith that of the uncorrected, truncated density matri-ces, which is 0 . ± .
18, as evidence of the improvement.For the practical design of compressive sensing-basedchannel correction for optical and quantum communica-tion systems, we considered the mean negativity of thecorrected density matrices as a function of the outputdimension N out . In Fig. 6, we display three curves fordifferent scintillation strengths W (by changing the prop-agation distance). N _out00.20.40.60.81 N e g a t i v i t y W = 0.5W = 1W = 2
FIG. 6. The negativity of the corrected density matrix as afunction of N out for different scintillation strengths W . One can see that the mean negativity of the correcteddensity matrices increases with N out , indicating the effectof the chosen size of the output Hilbert space. Further-more, we observe that the mean negativity becomes sat-urated. The saturation level and the rate of saturationdepends on the scintillation strength. These observations suggest a possible future improve-ment to the proposed scheme. By applying methods suchas deep learning [43, 44] or other generic machine learn-ing algorithms [45, 46], one may be able to train the sys-tem to determine an optimal value for N out , given certainturbulence conditions. V. CONCLUSIONS
Using a numerical analysis, we demonstrate the per-formance of a proposed compressive sensing-based chan-nel correction method. A classically non-separable state,consisting of three OAM modes with different wave-lengths, was used as input to a numerically simulatedturbulent free-space channel. Using compressive sensing-based state tomography, we reconstructed the outputstate and used it to determine the Kraus operator matrixfor the channel. The singular value thresholding tech-nique was used for the compressive sensing algorithm.The results show that compressive sensing drastically re-duces the number of measurements required for the char-acterization of the turbulent channel. Although the chan-nel estimation process uses classical light, it determinesthe Kraus operator matrix for the channel and thus al-lows it to be used for quantum communication. Conse-quently, the proposed scheme would be useful in the de-sign of quantum communication systems that are basedon the OAM states of light.As a further study, we intend to investigate the use ofdeep learning or machine learning as a method to informthe system on the appropriate dimension of the Krausoperator matrix. It could be done by training a model todetermine an appropriate value of this dimension, basedon the parameters of the turbulent free-space channel.
ACKNOWLEDGEMENTS
CMM acknowledges support from the CSIR NationalLaser Centre. The research for this work was supportedin part by the National Research Foundation (NRF) ofSouth Africa (Grant Number: 118532). [1] G. Gibson, J. Courtial, M. J. Padgett, M. Vasnetsov,V. Pas’ko, S. M. Barnett, and S. Franke-Arnold, “Free-space information transfer using light beams carrying or-bital angular momentum,” Opt. Express , 5448 (2004).[2] J. T. Barreiro, T.-C. Wei, and P. G. Kwiat, “Beatingthe channel capacity limit for linear photonic superdensecoding,” Nature physics , 282 (2008).[3] M. Mirhosseini, O. S. Maga˜na-Loaiza, M. N. O’Sullivan,B. Rodenburg, M. Malik, M. P. J. Lavery, M. J. Pad-gett, D. J. Gauthier, and R. W. Boyd, “High-dimensionalquantum cryptography with twisted light,” New J. Phys. , 033033 (2015). [4] L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, andJ. P. Woerdman, “Orbital angular momentum of lightand the transformation of laguerre-gaussian laser mode,”Phys. Rev. A , 8185 (1992).[5] C. Paterson, “Atmospheric turbulence and orbital angu-lar momentum of single photons for optical communica-tion,” Phys. Rev. Lett. , 153901 (2005).[6] L. C. Andrews and R. L. Phillips, Laser Beam Propaga-tion Through Random Media (SPIE, Washington, 1998).[7] C. M. Mabena and F. S. Roux, “Optical orbital angularmomentum under strong scintillation,” Phys. Rev. A ,013828 (2019). [8] B. J. Smith and M. G. Raymer, “Two-photon wave me-chanics,” Phys. Rev. A , 062104 (2006).[9] C. Gopaul and R. Andrews, “The effect of atmo-spheric turbulence on entangled orbital angular momen-tum states,” New J. Phys. , 94 (2007).[10] A. K. Jha, G. A. Tyler, and R. W. Boyd, “Effects ofatmospheric turbulence on the entanglement of spatialtwo-qubit states,” Phys. Rev. A , 053832 (2010).[11] F. S. Roux, “Infinitesimal-propagation equation fordecoherence of an orbital-angular-momentum-entangledbiphoton state in atmospheric turbulence,” Phys. Rev. A , 053822 (2011).[12] N. D. Leonhard, V. N. Shatokhin, and A. Buchleit-ner, “Universal entanglement decay of photonic-orbital-angular-momentum qubit states in atmospheric turbu-lence,” Phys. Rev. A , 012345 (2015).[13] T. Br¨unner and F. S. Roux, “Robust entangled qutritstates in atmospheric turbulence,” New J. Phys. ,063005 (2013).[14] B.-J. Pors, C. H. Monken, E. R. Eliel, and J. P. Woerd-man, “Transport of orbital-angular-momentum entangle-ment through a turbulent atmosphere,” Opt. Express ,6671 (2011).[15] M. Malik, M. O’Sullivan, B. Rodenburg, M. Mirhosseini,J. Leach, M. P. J. Lavery, M. J. Padgett, and R. W. Boyd,“Influence of atmospheric turbulence on optical commu-nications using orbital angular momentum for encoding,”Opt. Express , 13195 (2012).[16] A. Hamadou Ibrahim, F. S. Roux, M. McLaren, T. Kon-rad, and A. Forbes, “Orbital-angular-momentum entan-glement in turbulence,” Phys. Rev. A , 012312 (2013).[17] S. K. Goyal, A. Hamadou Ibrahim, F. S. Roux, T. Kon-rad, and A. Forbes, “The effect of turbulence onentanglement-based free-space quantum key distributionwith photonic orbital angular momentum,” J. Opt. ,064002 (2016).[18] Y. Zhang, S. Prabhakar, A. Hamadou Ibrahim, F. S.Roux, A. Forbes, and T. Konrad, “Experimentally ob-served decay of high-dimensional entanglement throughturbulence,” Phys. Rev. A , 032310 (2016).[19] M. A. Nielsen and I. L. Chuang, Quantum Computationand Quantum Information (Cambridge University Press,Cambridge, England, 2000).[20] M. Mohseni, A. T. Rezakhani, and D. A. Lidar,“Quantum-process tomography: Resource analysis of dif-ferent strategies,” Phys. Rev. A , 032322 (2008).[21] G. A. Tyler and R. W. Boyd, “Influence of atmosphericturbulence on the propagation of quantum states of lightcarrying orbital angular momentum,” Opt. Lett. , 142(2009).[22] J. A. Anguita, M. A. Neifeld, and B. V. Vasic,“Turbulence-induced channel crosstalk in an orbital an-gular momentum-multiplexed free-space optical link,”Appl. Opt. , 2414 (2008).[23] D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf.Theory , 1289 (2006).[24] E. J. Cand`es, J. K. Romberg, and T. Tao, “Stable signalrecovery from incomplete and inaccurate measurements,”Commun. Pure Appl. Math. , 1207 (2006).[25] E. J. Cand`es, in Proceedings of the InternationalCongress of Mathematicians (2006), vol. 3, pp. 1433–1452.[26] E. J. Cand`es and B. Recht, “Exact matrix completionvia convex optimization,” Found. Comput. Math. , 717 (2009).[27] D. Gross, Y.-K. Liu, S. T. Flammia, S. Becker, and J. Eis-ert, “Quantum state tomography via compressed sens-ing,” Phys. Rev. Lett. , 150401 (2010).[28] F. Tonolini, S. Chan, M. Agnew, A. Lindsay, andJ. Leach, “Reconstructing high-dimensional two-photonentangled states via compressive sensing,” Sci. Rep. ,6542 (2014).[29] G. A. Howland, P. B. Dixon, and J. C. Howell, “Photon-counting compressive sensing4laser radar for 3d imag-ing,” Appl. Opt. , 5917 (2011).[30] G. A. Howland and J. C. Howell, “Efficient high-dimensional entanglement imaging with a compressive-sensing double-pixel camera,” Phys. Rev. X , 011013(2013).[31] G. A. Howland, D. J. Lum, M. R. Ware, and J. C. How-ell, “Photon counting compressive depth mapping,” Opt.Express , 23822 (2013).[32] R. L. Kosut, A. Shabani, and D. A. Lidar, “Robust quan-tum error correction via convex optimization,” Phys.Rev. Lett. , 020502 (2008).[33] A. Shabani, R. L. Kosut, M. Mohseni, H. Rabitz, M. A.Broome, M. P. Almeida, A. Fedrizzi, and A. G. White,“Efficient measurement of quantum dynamics via com-pressive sensing,” Phys. Rev. Lett. , 100401 (2011).[34] B. Ndagano, B. Perez-Garcia, F. S. Roux, M. McLaren,C. Rosales-Guzman, Y. Zhang, O. Mouane, R. I.Hernandez-Aranda, T. Konrad, and A. Forbes, “Char-acterizing quantum channels with non-separable statesof classical light,” Nature Physics , 397 (2017).[35] C. M. Mabena and F. S. Roux, “High-dimensional quan-tum channel estimation using classical light,” Phys. Rev.A , 053860 (2017).[36] D. P. Greenwood, “Bandwidth specification for adaptiveoptics systems,” J. Opt. Soc. Am. A , 390 (1977).[37] J.-F. Cai, E. J. Cand`es, and Z. Shen, “A singular valuethresholding algorithm for matrix completion,” SIAM J.Optim. , 1956 (2010).[38] M. Jiang, S. Luo, and S. Fu, “Channel-state duality,”Phys. Rev. A , 022310 (2013).[39] J. M. Martin and S. M. Flatt´e, “Intensity images andstatistics from numerical simulation of wave propagationin 3-d random media,” Appl. Opt. , 2111 (1988).[40] D. L. Knepp, “Multiple phase-screen calculation of thetemporal behavior of stochastic waves,” Proc. IEEE ,722 (1983).[41] R. G. Lane, A. Glindemann, and J. C. Dainty, “Simula-tion of a kolmogorov phase screen,” Waves in RandomMedia , 209 (1992).[42] D. L. Fried, “Optical resolution through a randomly in-homogeneous medium for very long and very short expo-sures,” J. Opt. Soc. Am. , 1372 (1966).[43] I. Goodfellow, Y. Bengio, and A. Courville, Deep learning (MIT press, 2016).[44] M. A. Nielsen,
Neural networks and deep learning , vol. 25(Determination press San Francisco, CA, USA:, 2015).[45] E. Alpaydin,
Introduction to Machine Learning (MITPress, 2004).[46] Y. Ismail, I. Sinayskiy, and F. Petruccione, “Integratingmachine learning techniques in quantum communicationto characterize the quantum channel,” J. Opt. Soc. Am.B36