Robust entangled qutrit states in atmospheric turbulence

The entangled quantum state of a photon pair propagating through atmospheric turbulence suffers decay of entanglement due to the scintillation it experiences. Here we investigate the robustness against this decay for different qutrit states. We use an infinitesimal-propagation equation to obtain the density matrix as a function of the propagation distance and we use the tangle to quantify the entanglement between a pair of qutrits. The evolution of various initial states as they propagate through turbulence is considered. Using optimization of the initial parameters, we obtain expressions for bipartite qutrit states that retain their initial entanglement longer than the initially maximally entangled states.


Introduction
Free-space quantum communication is one of the major components in the new quantum information technology revolution. A challenge that confronts free-space quantum communication is the decay of entanglement that entangled photons experience due to scintillation while propagating through a turbulent atmosphere. The approaches to overcome this challenge are either based on methods to correct the output optical field, such as using adaptive optics [1], or on choosing an input optical field that is to some extent robust against entanglement decay in atmospheric turbulence. Here we consider the latter approach. In other words, we investigate to what extent a quantum state can be optimized so that it will retain as much of its initial entanglement as possible, while propagating through turbulence. We'll refer to such an optimized quantum state as a robust state.
In this paper we consider bipartite qutrit states in the Laguerre-Gaussian (LG) modal basis, which we restrict to the three elements that have a radial modal index of zero (p = 0) and azimuthal modal indices of ℓ = 1, 0, −1, which gives a three-dimensional Hilbert space per photon H 3 . The orbital angular momentum (OAM) associated with an LG mode is proportional to ℓ. If, instead of choosing |ℓ| ≤ 1, we use larger values of ℓ (higher OAM), the entanglement is expected to last longer [2]. This follows from the observation that the coupling from one OAM mode into another is stronger the closer the two OAM values are to each other [3,4].
Here the infinitesimal-propagation equation (IPE) [5] is used to calculate the zdependence of the corresponding density matrix. We assume that the two photons propagate through different uncorrelated regions of atmospheric turbulence. The turbulence is modeled by the Kolmogorov power spectral density [6].
We quantify the entanglement between the two photons by the tangle τ [7,8] as a function of z. The tangle is equal to the square of the concurrence for pure states and gives a lower bound for the square of the concurrence of mixed quantum states. To find the most robust qutrit states, we optimize the tangle at a given propagation distance z > 0. The result is expressed in terms of the parameters for the initial pure state that will give the maximum tangle at that propagation distance.
It has been shown that if only one partie of a bipartite (or multipartite) state passes through a dissipative channel, the final amount of entanglement is proportional to the initial entanglement of the state [9,10,11,12]. For such cases the initially maximally entangled states are also the most robust states. This argument has been extended to the case where all parties of an entangled state pass through dissipative channels, provided that these channels are trace preserving [13]. However, if both parties of an entangled bipartite state pass through uncorrelated dissipative channels that are not trace preserving, only an upper bound exists for the evolution of the entanglement [9,10,11]. Here we'll show in particular that for two photons propagating through a turbulent atmosphere, corresponding to two qutrit states passing through uncorrelated dissipative channels that are not trace preserving, the initially maximally entangled state is not also the most robust state.

Initial state
The most general pure bipartite qutrit state, defined in terms of our Hilbert space where c m,n represents complex coefficients and |m and |n are the OAM eigenstates for the respective subsystems, with m, n ∈ {1, 0, −1}. The normalization of the initial state implies that m,n |c m,n | 2 = 1.
We use a parameterization that implicitly obeys this normalization condition. It also removes an overall phase factor and incorporates the symmetry associated with an interchange of the ℓ = 1 and ℓ = −1 states. This parameterization is given by Due to the normalization condition and the removal of the overall phase we have 16 real parameters: 8 angles (k a , ..., k h ) and 8 phases (q a , ..., q h ). The density matrix ρ(z) is obtained as a function of the propagation distance by solving the IPE [5] for an initial density matrix defined as ρ(0) = |ψ ψ|, where |ψ is given in (1), using the coefficients defined in (3). Although, we assume that the initial state (1) is a pure state, the scintillation process causes it to become mixed during propagation, which necessitates a density matrix approach.

The IPE
The IPE [5] is a set of coupled first order differential equations that describes the evolution of the density matrix for a biphoton in the OAM basis as a function of the propagation distance z through a turbulent medium. It represents a multiple phase screen approach, as opposed to the single phase screen approach [3]. As a result it can simulate both phase and intensity fluctuations, whereas the single phase screen approach can only simulate phase fluctuations.
In the IPE the turbulence model is specified in terms of a power spectral density. Here we'll use the Kolmogorov power spectral density [6], given by where C 2 n is the refractive index structure constant, which quantifies the strength of the turbulence. In the solution of the density matrix, obtained from the IPE the strength of the turbulence is contained in a dimensionless parameter defined by where w 0 is the radius of the Gaussian envelope at the waist of the beam, λ is the wavelength of the photons and z R is the Rayleigh range, defined by For the three-dimensional bipartite case the density matrix is a 9 × 9 matrix. Therefore, the IPE for this density matrix is given as 81 first order differential equations that are divided into several sets of coupled equations and a few uncoupled equations. To solve the sets of coupled differential equations we use a perturbative approach, exploiting the fact that the terms that couple the different equations are multiplied by constants that are invariably an order of magnitude smaller that the uncoupled terms.
After solving the differential equations to obtain the density matrix as a function of the propagation distance, one finds that its dependence on the propagation distance is governed by two functions that are defined by the following integrals and where we define a normalized propagation distance t = z/z R . For convenience we separate H(t) into its real and imaginary parts H(t) ≡ H r (t) + iH i (t). These integrals can be solved to give expressions in terms of hypergeometric functions , ·) represents Barnes's extended hypergeometric functions [14].  In figure 1 we show the exact solutions for Z(t), H r (t) and H i (t) in comparison with their approximations, indicating that if propagation is limited to t 1/3, one can use the approximate solutions. The requirement that t < 1/3 implies that the entanglement decays quickly, which calls for strong turbulence conditions. In what follows we'll consider both the cases where the entanglement decays at t < 1/3 (strong turbulence), using the approximate expressions in (9), (10) and (11), and the cases where the entanglement only decays at t > 1/3 (weak turbulence), using the full expressions in (9), (10) and (11).

Trace
The solution that is obtained from the IPE represents a truncated density matrix, because only three of the infinite number of basis elements are retained. This is experimentally analogous to post selection [15] of measurement results in the three dimensional subspace. During propagation through turbulence the energy in the initial modes is scattered into higher order modes that are not contained in the truncated density matrix. As a result the density matrix is subnormalized -its trace is smaller that 1. This, in turn, causes the trace of the truncated density matrix to decrease as a function of the propagation distance. Moreover, the function of the trace depends on the initial state. However, due to the symmetries among the chosen OAM basis elements, the trace function only contains the parameters k h and k d . The evolution of the trace as a function of t is plotted in figure 2 for three extreme cases: k h = 0, which is independent of k d ; k h = π/2, k d = 0; and k h = k d = π/2. These cases are considered for strong turbulence (σ = 25) in figure 2a and for weak turbulence (σ = 0.25) in figure 2b. The plots show that states with ℓ = 0 scatter slower than those with ℓ = ±1. In fact, the slowest decay of the trace is obtained for k h = 0, which represents a state consisting of only |0, 0 and, therefore, is completely separable.

Tangle
For pure two-dimensional bipartite systems the concurrence [16] is a suitable measure of entanglement, but for mixed states in higher dimensions it is computationally demanding to calculate the concurrence directly through the construction of a convex roof [17]. Instead, we use the tangle τ {ρ}, which is equal to the square of the concurrence for pure states and gives a lower bound for that of mixed quantum states [7]. The tangle is calculated via the purities of the (reduced) density matrices where ρ A and ρ B are the respective reduced density matrices of the two subsystems. For a maximally entangled state τ = τ max = 2(d − 1)/d, where d is the dimension of the subsystems' Hilbert spaces and for a separable state τ = τ min = 0.

Bell states
In the two-dimensional case, where diametric azimuthal indices ℓ = ±1 are used, it can be shown that the entanglement of the four maximally entangled Bell states, quantified by the concurrence, all decay equally [5]. This can be understood as a result of the fact that the scattering is symmetric with respect to ℓ = 1 and ℓ = −1. In the threedimensional case there is a third possible azimuthal index, which does not share this symmetry. In our case we chose the third index to be ℓ = 0, and we quantify the entanglement by the tangle. As a result, we find that the 12 Bell states in the twodimensional subspaces of our three-dimensional Hilbert space form three sets based on their decay curves are the Bell states of the two-dimensional subspaces spanned by r and s. The three sets of subspace Bell states give three different decay curves, as shown in figure 3. Figure 3 shows that the initial amount of entanglement does not uniquely determine the remaining entanglement after a propagation distance t, as was found for the twodimensional situation [5]. The Bell states of Set 3 maintain a non-zero entanglement for longer than the states of Sets 1 and 2. This can be understood as a consequence of the fact that there is a stronger coupling in this case between modes having azimuthal indices with different magnitudes than between those with identical magnitudes, which in turn follows from how the coupling strength between two OAM modes depends on the difference between their OAM-values.

Initially maximally entangled states
Next we investigate the decay of entanglement of initially maximally entangled qutrit states. Depending on their decay curves, one can identify two sets of initial qutrit states, given by where j, k ∈ {−1, 1}, and φ 1 , φ 2 , φ 3 , ψ 1 , ψ 2 and ψ 3 are arbitrary phases. Calculating the tangle for states with the initial conditions given as State 1 in (18), we find that the phases always appear in the combination φ = φ 1 + φ 2 − 2φ 3 in the expression of the tangle. The expression for the combined phase that maximizes the tangle at t = t 0 is given by while φ = φ opt (t 0 ) − π minimizes the tangle. Note that, although the combined phase that minimizes or maximizes the tangle depends on the propagation distance, this dependence is very small in the case of strong turbulence, because the entanglement decays quickly, for t < 1/3. As a result we can use the small t approximation, which implies that φ ≈ 0 and φ ≈ π, respectively, maximizes and minimizes the tangle. We'll denote the cases with minimum and maximum tangle as State 1a and State 1b, respectively. In the case where the initial state is given by State 2 the tangle is totally independent of the phases ψ 1 , ψ 2 and ψ 3 . The decay curves of State 1a, State 1b and State 2 are shown in figure 4, both in weak and strong turbulence. One finds that states represented by State 2 are more robust than those represented by State 1. This again follows from the stronger coupling between differing magnitudes of azimuthal indices compared to the coupling strength between modes with identical magnitudes of azimuthal indices.
Comparing the curves for strong turbulence (figure 4a) to the ones for weak turbulence (figure 4b), one can see that the entanglement of states for larger propagation distances decays comparatively faster -if we would have used the approximations for Z(t), H r (t) and H i (t), given in (9), (10) and (11), respectively, to calculate the tangle evolution for weak turbulence (i.e. when σ = 0.25), we would have found that the entanglement was maintained much longer than it does for the exact expressions. This can be seen from the fact that the tangle for weak turbulence decays nearly linearly all the way to zero entanglement, whereas for strong turbulence it decays linearly up to t ≈ 0.005, followed by a nonlinear slower decay. Thus, if one can operate in a scenario where σ and z R have larger values, then one would have robust states for comparatively larger propagation distances.

Most robust entangled state
Upon comparing the curves for initially maximally entangled qutrit states in figure 4 to those for the Bell states of the two-dimensional subsystems in figure 3, we find that the maximally entangled qutrit states are not the most robust qutrit states. The Bell states of Set 3 maintain their entanglement longer than any of the initially maximally entangled qutrit states in (18) and (19), even though the initial entanglement of the Bell states is smaller. As a result, if one wants to find the most robust entangled state for a specific t, it is not enough to consider only those states that are initially maximally entangled.  Figure 5: Optimization parameters µ 1 (t) for State 1 and µ 2 (t) for State 2 that gives the most robust entangled states |µ 1 (t) and |µ 2 (t) as a function of the normalized propagation distance t, in strong turbulence (a) and weak turbulence (b).
To find the most robust states we optimize the tangle at a particular point t = t 0 , for the general initial pure state, using the parameterization from (3). This optimization process is similarly to [8]. From the optimization results, expressed in terms of the parameters k a , ..., k h and q a , ..., q h , we found that all the states have one of the following two forms where we introduced the optimization parameters µ 1 and µ 2 in addition to the phases φ 1 , φ 2 , φ 3 , ψ 1 , ψ 2 and ψ 3 . These parameters are used to replace those given in (3) for the sake of simplicity. The optimization parameters µ 1 and µ 2 determine the degree of relative weighting between the maximally entangled states in (18) and (19) and the Bell states of the two-dimensional subspaces given in (13), (14) and (15). The curves of the optimized parameters µ 1 (t) and µ 2 (t) are shown in figure 5, for both weak and strong turbulence. The fact that the parameters are not independent of t means that there does not exist one single state that is the most robust state for all propagation distances. Instead there are different states that have the highest remaining entanglement at each specific value of t. For t = 0 the optimized states |µ 1 (0) and |µ 2 (0) are the maximally entangled states of (18) and (19), respectively.
Substituting (20) and the optimal values for µ 1 and µ 2 at t 0 = 0.03 (t 0 = 1.6) for strong (weak) turbulence into (21) and (22), we obtain states that maintain a nonzero entanglement much longer than the initially maximally entangled states of (18) and (19) or the Bell states in (13), (14) and (15 Figure 6: Comparison of the tangle evolution of initially maximally entangled states (State 1 and State 2) to those of the optimized states (|µ 1 (t) and |µ 2 (t) ) for strong turbulence optimized at t 0 = 0.03 (a) and for weak turbulence optimized at t 0 = 1.6 (b).
tangle are shown in figure 6. Note that State 1, which is less robust than State 2, can be optimized to produce state |µ 1 that is more robust for large propagation distances than the optimized state |µ 2 . It turns out that, when optimized, |µ 1 is the most robust qutrit state within our three-dimensional Hilbert space for propagation through turbulence.

Summary and conclusions
We have shown that for three-dimensional systems, the initial amount of entanglement does not uniquely define the amount of entanglement left after a non-zero propagation through a turbulent medium. This highlights the fact that the upper bound for the entanglement evolution [10,11] has to be treated with caution. The reason is that in our three-dimensional Hilbert space the coupling strength between differing magnitudes of azimuthal indices is different from the coupling strength between identical magnitudes. This follows from our choice of Hilbert space and the fact that the coupling strength decreases as the difference between the azimuthal modal indices increases.
By optimizing the tangle, starting from the most general initial pure state, we found that neither the initially maximally entangled states, nor one of the Bell states of the two-dimensional subsystems are the most robust states. Instead it is a coherent superposition between both. For each propagation distance one can find a specific initial entangled state that retains the most entanglement up to that propagation distance. Thus, there is not one single most robust state for all distances.
The effect of the optimization of the robustness of the qutrit states is still rather small (it only increases the propagation distance by a few percent). However, it is reasonable to expect that the effect of such an optimization would become more significant for higher dimensional Hilbert spaces, making it worth doing the optimization to find robust quantum states in higher dimensional Hilbert spaces and thereby improve the performance of quantum communication systems. As a result we believe that the optimization of the robustness of higher dimensional quantum states will have a significant impact in practical free-space quantum communication systems.