Entanglement transmission through dense scattering medium

Scattering effects are ubiquitous in practical wireless optical links. Here a transmission model with complete consideration of scattered light and beam wandering effects for underwater link is developed, with the aim to completely characterize the received quantum state of light through dense scattering medium. Based on this model, we show the influence of scattered photons on the improvement of the entanglement after transmission through turbid water may vary for different copropagation scenarios, i.e., the contribution of scattered light on entanglement transmission may be turned from positive to negative, with increase of the strength of underwater beam wandering. And the attenuation coefficient and aperture size are found to be the dominant factors affecting the entanglement through underwater link. While for the counterpropagation scenario, the scattered photons will severely deteriorate the entanglement transmission especially for the high-loss scattering links. These findings may shed light on quantum entanglement transmission and help to develop its applications through dense scattering medium.


Introduction
Entanglement is one of the most important nonclassical properties of quantum light, which has been of great interest due to the key fundamentals of quantum physics and for various practical applications, such as quantum communication [1,2], quantum metrology [3] and quantum sensing [4]. In particular, the studies of propagations of quantum light through scattering medium are very necessary for satellite-based and submarine-based quantum communications and other applications, since the scattering effects commonly exist in atmospheric and underwater wireless optical links.
Literally, except for the irreversible light absorption, the dense scattering medium in atmosphere or water will lead to remarkable photon scattering, where the well-collimated laser beam may become largely diffused when it reaches the receiver. Despite the scattering effects, another main factor that affects the propagation of light through atmosphere and water is turbulence [5][6][7], which will lead to beam wandering [8][9][10] and even signal interruption. In particular, the turbulence is caused by variations in the refractive index due to inhomogeneities in temperature and pressure in the atmosphere and due to large scale variations in temperature and salinity in the water [5]. In atmosphere, long range (>1 km) wireless optical links are affected mainly by turbulence [7]. While for most cases of underwater scenarios, particulate scattering dominates except for clear water and short ranges [5,6]. However, realistic unstable adjustment between the radiation source and receiver in both atmosphere and water will also lead to beam wandering [8]. It should be noted that the absorption and scattering on particulates are not changed randomly. For instance, scattering of light on water drops in the atmosphere has a random character, which can be described by the corresponding spectrum of the correlation function [11][12][13][14].
So far, the nonclassical properties of light propagating through atmosphere have been successively studied in references [8,[15][16][17][18][19][20], where some interesting phenomena have been found that the permanently fluctuations of channel losses caused by beam wandering may benefit the preservation of nonclassical properties, including the quantum entanglement [8,[16][17][18][19]. However, the scattering effects are not treated thoroughly in these theoretical models, since in some scenarios the inconspicuous influences of scattered effects on light transmissions are neglected. While for the underwater links, the studies mainly focus on the classical effects of the water on propagating optical signals [21][22][23][24], the related experimental demonstration of the feasibility of entanglement transmission [25], the performance analysis and related experimental demonstrations of quantum key distribution [26][27][28][29][30] and quantum state transmission [26,31,32]. The knowledges of the influences of nontrivial scattering effects on the nonclassical properties, especially the entanglement, for atmospheric and underwater environments are almost blank.
In this paper, the analytical solution of irradiance distribution model for underwater wireless optical channel is developed by considering the contribution of scattered light with finite-size aperture and the effects of beam wandering. And the effects of underwater scattering and beam wandering on violations of Bell inequalities are investigated with numerical simulation and curve fitting methods. We show in the scenario of copropagation, there is a transition for the influence of scattered photons on entanglement transmission, i.e., the contribution of scattered light on transmission may be changed from slightly improving to deteriorating the entanglement of the quantum state transmitted through turbid water, with increase of the mean link losses or fundamentally the strength of beam wandering effects. And the entanglement degradation is mainly attributed to the attenuation coefficient of the scattering medium and the aperture size of receiver. While for the counterpropagation scenario, the contribution of scattered light will mainly deteriorate the entanglement transmission, and this effect will be more severe for high-loss links. Moreover, the double-click events, where each receiver registers for the reasons of dark counts, stray light, and multiphoton entangled pairs, will generally lead to negative effects on the preservation of entanglement transmission for both of copropagation and counterpropagation scenarios, and the effects are more obvious for the latter scenario. These results provide a more complete insight into entanglement transmission and its applications through realistic dense scattering medium.
The paper is organized as follows. In section 2 we first introduce the transmission model for underwater wireless optical channel with the scattering and beam wandering effects. Then the output quantum states after transmission are given. The transmission efficiencies and violations of Bell inequalities after propagation through the underwater wireless optical links are investigated and analyzed under realistic conditions in section 3. Summary and conclusions are given in section 4.

Transmission model
In order to study the scattering effects on the nonclassical properties for the wireless optical links, the initial challenges lie in understanding and modeling the issues due to particulate scattering in wireless optical transmission environments. Fortunately, by using the solution of analytical computation of the irradiance distribution with assumptions of small angle approximation (SAA) [33], a theoretical prediction model for spatial spreading effects of scattering on the underwater optical communication signals is derived [21], which lays an analytical foundation for further investigation of scattering effects on the nonclassical properties.
The geometry for the underwater wireless optical communication links is shown in figure 1. Here, we aim to compute the analytical irradiance distribution of this source transmitting through water on the plane perpendicular to the beam axis at z = z rec . It should be mentioned that SAA assumptions contain three aspects [21]. The first one is that the phase function is highly peaked at small angles and there are no contributions from scattering events at large angles, which has been validated for realistic environment in reference [34]. The second one is that the irradiance distribution is symmetric with the beam axis. The last one is there are no contributions from light scattered for the positions z > z rec . Based on the assumptions of SAA, the irradiance distribution at a particular wavelength can be obtained as where J 0 is the Bessel function, S(ν, z rec ) is the medium optical transfer function (OTF), and I 0 (ν, z rec ) = P 0 exp[−V 0 (z rec )ν 2 /2] is the Gaussian irradiance distribution of laser source in free space at position z = z rec with initial power P 0 and variance of the Gaussian source V 0 (z rec ). The irradiance distribution I(r, z rec ) can be decomposed into nonscattered I N (r, z rec ) and scattered I S (r, z rec ) components I(r, z rec ) = I N (r, z rec ) + I S (r, z rec ).
In particular, the irradiance distribution for the nonscattered light can be expressed as Figure 1. Geometry for the underwater wireless optical links with beam wandering and scattering effects. (a) The geometry for irradiance distribution with scattering effects; (b) the receiving plane for finite-size aperture with beam wandering and scattering effects. r denotes the deflection distance from beam center to the aperture center; a is the radius of aperture; z is the propagation distance.
where I 0 (r, z rec ) is the Hankel transform of I 0 (ν, z rec ) with the expression I 0 (ν, z rec ) = P 0 2πV 0 (z rec ) exp −r 2 /2V 0 (z rec ) , and τ = z src z rec c(z)dz reflects the effects of nonhomogeneity. Here c(λ) = a(λ) + b(λ) is the beam attenuation coefficient with the absorption coefficient a(λ) and the scattering coefficient b(λ) in the units of inverse meter. Similarly, the scattered component is where p(ν) is the Hankel transform of the scattering phase function p(ν) = 1 2 π 0 p(θ)J 0 (νθ)θdθ. It governs the amount of the scattered light which reaches a receiver at a position r meters away from the beam axis due to multiple small angle forward scattering.
The density of transmission efficiency with respect to the position r through water can be calculated as The transmission efficiency can be obtained by integration of equation (5) in the aperture area as where w is the radius of aperture of the receiver detector. While ignoring the contribution of scattering light, the density of transmission efficiency T 2 (1) in equation (7) should be degraded to I N (r, z)/P 0 .

Output quantum state
Without loss of generality, we assume that the propagation of quantum light through the scattering medium is a linear attenuating process. The relation for the annihilation operators of the input and output fields, â in and â out can be expressed asâ out = Tâ in + √ 1 − T 2ĉ [20] in linear quantum optics, whereĉ is the modes of environment. The transmission losses are caused by light absorption and scattering as well as by the diffusing due to particulate scattering. Thus we can obtain the output quantum state based on Glauber-Sudarshan P representation [35], where the quantum-state input-output relation resembles the corresponding relation in classical optics with the expression where P(T) is the probability distribution of the transmission coefficient (PDTC), which reflects the fluctuations of transmission losses, and Φ out (β) and Φ in (Tβ) are the characteristic function of output and attenuated input fields, respectively. The characteristic function Φ(β) is related to the corresponding density operatorρ as For further considerations we assume the absence of dephasing, and the transmission coefficient is real and positive with T ∈ [0, 1], which reflects that the commutation relations preserve for the underwater transmission links. When beam deflection happens due to the imperfect adjustment of the radiation source or by weak turbulence, the transmission coefficient will fluctuate. Here we assume beam incidence is normal to the aperture plane, and the position of the beam-center (x 0 , y 0 ) randomly fluctuates according to the two-dimensional Gaussian distribution where d = (d x , d y ) is the mean position of beam center, and σ is the standard variance of beam deflection.
The beam deflection distance r = x 2 0 + y 2 0 will fluctuate according to the Rice distribution where d = d 2 x + d 2 y is the distance between aperture and fluctuation centers, I 0 is the modified Bessel function.
The random variables r and T 2 are related to each other via equation (7). If we denote the inverse function of T 2 as r(T), and the derivative of r(T) as r (T). This will yield that the PDTC reads as Thus, the output quantum state through underwater scattering links can be obtained through the equations (8) and (11). Moreover, we can get the mean transmission coefficient and transmission efficiency as It should be noted here we consider the beam deflections are mainly originated from the integrated effects of water turbulence and the unstable adjustments between the radiation source and receiver, i.e., both variances of these beam deflections should be added up. In particular, the later factor is one of the major contributions for submarine-based communication, and also the critical temperature difference in water will lead to the transition of Rayleigh-Bénard (RB) convection [36] to remarkable turbulent, noted as RB turbulence [10]. For the beam deflection caused by imperfect adjustment of the radiation source, the variance is σ 2 ≈ σ 2 θ z 2 with σ 2 θ the variance of the source-deflection angle [8]. While for the RB [10,37], where C 2 n is the index of refraction structure constant which indicates the strength of optical turbulence, μ 0 is the beam waist radius, and κ 0 = 2π/z 0 with z 0 the outer scale of RB turbulence. For more generalized considerations of the proposed transmission model in different application scenarios and water turbulence models, the strength of beam wandering effects are directly described with different deviations of beam deflections as analyzed in reference [8] in the followings.

Underwater transmission efficiency
In order to compute the transmission coefficient and PDTC for underwater optical wireless links. We should firstly obtain the analytical expression of scattering phase function of water. The two-parameter approximation for the marine phase function derived by Fournier and Forand (FF) [38] was recognized as one of the most accurate representation of real oceanic scattering functions. In turbid water, where the contribution of pure seawater is generally negligible, the FF function can usually be fitted to the scattering data for the whole seawater [39]. The recent form of FF phase function [39,40] can be expressed as 3(n−1) 2 , u = 2 sin θ 2 , and 1 < n < 1.35 is the refractive index of the particles relative to seawater and 3.5 < m < 5 is the slope of the power-law size distribution f(D) = kD −m . The minimum error for the second-order fitting of the FF function to the turbid seawater particle scattering function is located at m = 3.59, n = 1.073 [39]. In the followings, we use these parameters in simulations.
Here the wavelength of the transmitted light is 532 nm, the beam waist radius μ 0 is set as 2 mm, and propagation distance is z = 5 m. The variance of the Gaussian source can be obtained as where μ(z) is the variation of the beam waist μ(z) = μ 0 1 + λz . We can firstly calculate the function of f zp (ν) by using equations (6) and (13). Considering the computation complexity of integration, here we obtain an approximate analytical expression of f zp (ν) by using numerical simulation and high-precise curve fitting as (see appendix A for details) In order to explore the scattering effects on the transmission efficiency, we consider here the cases of two different scattering coefficients, i.e., b(λ, m, n) = 0.45 and b(λ, m, n) = 0.65 for the same attenuation coefficient c(λ, m, n) = 0.8 with and without the contribution of scattering photons. And the aperture radiuses are set as 10 mm and 30 mm with 180 • field of view (FOV), respectively. We further assume that the water is homogenate. Based on the equations (7) and (14), the transmission efficiencies T 2 for different values of beam-deflection distance r are obtained numerically as the discrete circle points shown in figure 2. In order to take further analysis of the entanglement after transmission through the underwater optical channel, we also perform high-precise curve fitting as shown the curves in figure 2 and obtain an approximate analytical expressions for transmission efficiencies (see appendix B for details). Apparently, the transmission efficiency will decrease dramatically with a certain small value of beam-deflection distance r when the contribution of scattered light is neglected. Actually, the realistically scattered photons can slightly improve the transmission efficiency, which is more obviously for larger scattering coefficient b. And it shows apparently in figure 2 that the existed scattering effects can help to build the communication links when the beam-deflections become severe. But the main factor which affects the transmission efficiency is still the attenuation coefficient c. Moreover, increase of the aperture radius is shown to be positive to improve the transmission efficiency.

Simulations of entanglement transmission 3.2.1. Copropagation scenario
Based on the input-output relation in equation (8) and PDTC in equation (11), we can consistently characterize the quantum states of the transmitted light through the scattering channel in water. In the following we firstly investigate the disturbances of the entanglement by scattering and beam wandering effects for the copropagation scenario, i.e., the entangled light copropagating in the same direction through the completely uncorrelated underwater channels, and without the incorporation of so-called double-click events [41][42][43]. It should be mentioned that the double-click events will lead to an indefinite value of the qubit for the corresponding measurements [41]. As demonstrated in references [1,8,16,17], the violations of Bell inequalities for realistic parametric down-conversion (PDC) source [1,2] can be well improved after transmitting through the turbulent atmosphere, which imply that the fluctuations of channel losses may preserve entanglement properties of light much better than standard loss channels [8,16].
The calculations of the violation of Bell inequalities here is similar to reference [8,16,17] but based on underwater channel (see appendix C for details). We can find the Bell parameter for the realistic PDC source is directly related to the averaging functions with respect to T 2 and T 4 . Here the formula for Bell parameter in Clauser-Horne-Shimony-Holt (CHSH) Bell-type inequality [44] is given by B ) with i, j = 1, 2 are the correlation coefficients for two sets of polarization angles. And the correlation coefficient E(θ A , θ B ) can be defined as where P same (θ A , θ B ) and P different (θ A , θ B ) are the probabilities to get clicks on detectors in the same and different channels of polarization analyzers, respectively. Thus, the violations of Bell inequalities can be easily evaluated with integrals of these functions. In particular, we consider in this context fluctuating losses are only caused by beam wandering. And all of the influences to the statistics of channel loss are modeled by PDTC in equation (11). We suppose here the total efficiency of the receiver module is η = 0.125, the mean number of stray-light and dark counts is N = 10 −5 .   When the mean losses of the underwater link are specified as T 2 = 6.3 × 10 −4 (32 dB), the violations of Bell inequalities for different cases with the same mean loss are depicted in figure 3. We can firstly verify that the fluctuation of transmission losses may also preserve entanglement properties of light even better than standard loss for underwater channels. And for specific aperture radius w and attenuation coefficient c in figure 3, it clearly shows that the contribution of scattered quantum states may deteriorate the entanglement after underwater transmission. While this relative deterioration phenomenon is more significant for larger aperture radius. Moreover, it can be seen from figure 3 that for the same aperture radius, the Bell parameters are almost the same for the cases b = 0.45, c = 0.8 and b = 0.65, c = 0.8, and the increase of scattering coefficient may lead to slightly more deterioration of entanglement for the same attenuation coefficient c. These results imply that the attenuation coefficient may be the dominant factor which affects the entanglement transmission through the scattering medium.
For further illustration, we consider a more realistic case for comparison, i.e., the coastal ocean seawater case [21], where the scattering coefficient is b = 0.219, attenuation coefficient is c = 0.398 in the FF function. The comparison of Bell parameters are shown in figure 4 with the same mean loss of 32 dB. The results further confirm that the entanglement properties of light can be preserved better for the scenarios of smaller attenuation coefficient c, whenever considering the contribution of scattered light. Moreover, we can find from figure 3 when the underwater channel and transmission distance are specified, increase the aperture radius may be naturally positive to preserve the entanglement property, whenever considering the scattering effects. From above simulations and analysis, we can find that the scattering effects may be positive and negative for the preservation of entanglement transmission through underwater wireless optical links, and this impact is dependent on the mean losses of the underwater links, which is quite related to the strength of underwater beam wandering. Also, the characteristics of scattering medium and aperture size will also affect the mean losses, and the exact bound for the mean loss needs further investigation. These findings are of vital importance for discrete-variable underwater quantum communication and other applications based on entanglement.
When considering the double-click events, the calculations of the Bell parameters should be renewed (see appendix C for details). The comparisons of the violations of Bell inequalities for the scattering channels (b = 0.45, c = 0.8, w = 0.01m) and constant-loss channels with the same mean loss of 20 dB are shown in figure 6. We can find that the scattered photons have relatively more positive effects on the preservation of entanglement transmission for the scenario of incorporation of double-click events. And when considering the constant-loss transmission, the Bell parameters will be almost the same for the cases with and without the incorporation of double-click events. However, the double-click events will generally lead to negative effects on the preservation of entanglement transmission, which in general accords with the atmospheric channel simulated in reference [17].

Counterpropagation scenario
Now we investigate the disturbances of the entanglement by scattering and beam wandering effects for the counterpropagation scenario, i.e., the entangled light propagating in different directions through the completely uncorrelated underwater channels. We first consider the scenario without the incorporation of double-click events. Without loss of generality, here we consider the case that only one photon of the entangled pair is sent through the underwater channel as the experiment in reference [2]. And for simplicity, one of the transmission efficiencies is assumed to be 1, the other is T 2 .
The calculations of the violation of Bell inequalities here is similar to reference [17] but based on underwater channel with and without the contribution of scattered light (see appendix D for details). We can find the Bell parameter for the realistic PDC source is also directly related to the averaging functions with respect to T 2 . Thus, the violations of Bell inequalities can be also evaluated with integrals of these functions. In particular, we also consider here fluctuating losses are only caused by beam wandering. And all of the influences to the statistics of channel loss are modeled by PDTC in equation (11). The total efficiency of the receiver module is also η = 0.125, aperture radius w = 0.01 m, and the mean number of stray-light and dark counts is N = 10 −5 .
When the mean losses of the underwater link are specified to 32 dB and 20 dB, the violations of Bell inequalities for the same mean losses are depicted in figure 7. For the case of high mean loss (32 dB), we can find that, in contrast to the copropagation scenario, the contribution of scattering photons here may lead to more serious deterioration of entanglement transmission (the Bell parameters are almost zero), and the increase of scattering coefficient b may lead to slightly more deterioration of entanglement for the same attenuation coefficient c. Generally, the entanglement properties of light will be preserved better for the copropagation scenario whenever considering the contribution of the scattered photons, even only one mode of the entangled light is considered to be transmitted here through the underwater channel for the counterpropagation. While for the case of lower mean loss (20 dB), the scattered photons have fewer effects on the preservation of entanglement transmission, but there is not a transition here for the influence of scattered photons on the entanglement transmission (the dotted and dashed curves overlap in figure 7(b)). Now we consider the double-click events, the calculations of the Bell parameters should be renewed (see appendix D for details). The comparisons of the violations of Bell inequalities for the scattering channels (b = 0.45, c = 0.8, w = 0.01m) and constant-loss channels with the same mean loss of 20 dB are shown in figure 8. We can find that the scattered photons also have obvious positive effects on the preservation of entanglement transmission as the copropagation case when considering the double-click events, such that the Bell parameters will be even larger than the case without considering the contribution of scattered photons. While discarding the double-click events, we can see from figure 7(b) that the Bell parameters are slightly smaller when considering the contribution of scattered states. This can be explained that the effect of double-click events may lead to more deterioration of entanglement transmission for the case without considering the contribution of scattered photons. And for the case of constant-loss transmission, the effect of double-click events will obviously decrease the Bell parameters. Generally, the negative influence of the double-click events on the entanglement transmission also in general accords with the atmospheric case in reference [17].

Conclusion
The transmission model for underwater scattering wireless optical channel is firstly developed based on the analytical solution of irradiance distribution with finite-size aperture and the beam wandering effects. Especially, the PDTC expression for underwater wireless optical communication link has been obtained based on the transmission model when assuming normal distribution of the beam positioning. Based on the consistent description of quantum-state input-output relations and assumption of absence of dephasing for the transmitted and scattered states entered into the detectors, we study the effects of realistic scattering and beam wandering on entanglement transmission through the underwater scattering links with numerical simulation and curve fitting methods.
We show that in the scenario of copropagation, the contribution of scattered light may deteriorate the entanglement of the quantum state when it is transmitted through high-loss scattering underwater channel, and the entanglement degradation should be mainly attributed to the attenuation coefficient of water. However, the scattered light may be also positive to preserve the entanglement when the mean losses of the underwater channel are relatively lower. Moreover, we can increase the aperture size to preserve the transmitted entanglement through the underwater scattering links. While for the counterpropagation scenario, the influence of the scattered light is mainly negative to the preservation of entanglement transmission, and this effect is more obvious for high-loss scattering underwater channel. Furthermore, we find that the effects of double-click events will generally deteriorate the transmitted entanglement for both of copropagation and counterpropagation scenarios, and the effects are more obvious for counterpropagation scenario. It should be mention that we are aware of that the analysis based on the method of numerical simulation and curve fitting may lead to inaccuracy of evaluations of Bell parameters. However, the incurred extremely small error will not affect the tendency of the entanglement transmission.

Appendix A. The calculation of f zp (ν)
We can calculate the function f zp (ν) by using equation (6). The related parameters are set as m = 3.59, n = 1.073, transmission distance is z = 5 m, beam waist radius μ 0 = 2 mm. We first calculate the values of f zp (ν) for several discrete points of ν, which are divided into two parts ν > 1 and 0 < ν 1. Here we use the curve fitting APP of MATLAB R2018a to perform curve fitting.
In particular, for ν > 1, a general model Gauss2 can well fit the relationship between ln(ν) and the corresponding calculated results f zp (ν) with the goodness of fit: SSE 1.66 × 10 −5 , R-square 1, adjusted R-square 1, and RMSE 0.000 9346. And for < ν 1, a general model Gauss3 can well fit the relationship between ν and the corresponding calculated results f zp (ν) with the goodness of fit: SSE 2.017 × 10 −7 , R-square 1, adjusted R-square 1, and RMSE 0.000 1246. The obtained piecewise function for f zp (ν) can be expressed as equation (14). The discrete points and the fitted curves are plotted in figure A1. We can find when ν become large, f zp (ν) will be close to zero.

Appendix B. The calculation of transmission efficiency T 2 (r)
Based on the obtained analytical piecewise function of f zp (ν), we can first compute the transmission efficiency by using the equations (5) and (7) in the paper for several discrete points of r in several intervals. Then we can also use the curve fitting APP of MATLAB R2018a to perform curve fitting for two parts, i.e., the nonscattering part and the scattering contributions.
We can then obtain the total transmission efficiency as T 2 N1 + T 2 S1 .
• For w = 0.03 m, z = 5 m, c = 0.8, b = 0.65. Then we consider the case without the contribution of scattering light. Since the total attenuation is same, the transmission efficiency T 2 N2 here is equal to T 2 N1 . Consider the case for the contribution of scattering light. For r ∈ [0, 0.035], a general model Exp2 can well fit the relationship between r and transmission efficiency with the goodness of fit: SSE 1.224 × 10 −8 , R-square 0.9918, adjusted R-square 0.9902, and RMSE 2.857 × 10 −5 . For r ∈ [0.035, ∞], a general model Power2 can well fit the relationship between r and the natural logarithm of transmission efficiency with the goodness of fit: SSE 0.006 091, R-square 0.9999, adjusted R-square 0.9999, and RMSE 0.0179. The fitted piecewise function for T 2 S2 can be expressed as  can be expressed as We can then obtain the total transmission efficiency as T 2 N5 + T 2 S5 .

Appendix C. The calculation of violation of Bell inequalities for copropagating
The Bell test experimental configuration for the copropagation scenario is originated from reference [1], where the quantum states are transmitted through the atmospheric channel. Here the transmission channel is changed to underwater one. As shown in figure A2, a PDC source generates entangled photon pairs, which are separated by a small time interval, and they are then sent to the receiver through the underwater channel. We assume that the time interval is much less than the time for which the water is changed. As a result, the two channels can be considered to be completely correlated, where the Figure A2. The Bell test scheme derived from reference [1], where the atmospheric channel is here replaced as an underwater one. HWP denotes half-wave plate, which changes the polarization angles A and B; PBS denotes polarizing beam splitter; D T A , D T B are the on/off detectors for the transmitted light, and D R A , D R B are the ones for the reflected light. Figure A3. The equivalent Bell test schemes derived from reference [8,16], where the atmospheric channel is replaced with underwater one. The principle of these tests are equivalent to the original one proposed in reference [1]. The two scenarios of the Bell test (a) and (b) correspond to two kinds of the events in the original scheme.
transmission coefficients are both T. In the receiver side B, the photons are collected by the telescope and split by a beam splitter and finally sent to the corresponding polarization analyzers.
Since the scheme in [1] only considers the events that the photons left different output ports of the beam splitter. So the time-separated modes in reference [1] can be replaced by spatial ones [8,16]. In particular, the resolved events in the original scheme correspond to the following events in figure A3: (a) clicks at the reflection port in side A and the transmission port in side B, (b) clicks at the transmission port in the side A and the reflection port in the side B. Considering the additional 50% detection losses caused by the beam-splitters, the final expression of Bell parameter is same as the one in [1]. The technical details of the calculations of the Bell parameter for the light generated by a PDC source is from the equivalent schemes presented in [8,16,17], which is depicted in figure A3. It should be noted that the transmission coefficients T A and T B for the two propagation directions should be completely correlated for the copropagation cases, and thus the PDTC in this case can be expressed as where we can then take T A = T B = T. In the following, we will first introduce the computation of Bell parameter without consideration of double-click events based on the equivalent schemes in [8,16,17]. The Bell parameter with the expression cannot exceed the value of 2 in any local realistic theory. But quantum light may violate this inequality. The correlation coefficients E(θ A , θ B ) in equation (C2) are defined as where P same (θ A , θ B ) = P T A ,T B (θ A , θ B ) + P R A ,R B (θ A , θ B ) is the probability to get clicks on detectors in the same channels of polarization analyzers, and P different (θ A , θ B ) = P T A ,R B (θ A , θ B ) + P R A ,T B (θ A , θ B ) is the probability to get clicks on detectors in different channels. Here P i A ,i B (θ A , θ B ) is the probability of registering clicks at the detectors i A = {T A , R A } and i B = {T B , R B } for the polarization angles θ A and θ B . According to the photodetection theory, these probabilities are given by where i A = j A , i B = j B ,ρ is the density operator, and is the positive operator-valued measure for the detector i A(B) , N is the mean number of stray-light and dark counts, η is the detection efficiency, : · · · : means normal ordering, are the photon-number operators at the outputs of the polarizing beam-splitters. The operator input-output relations of the polarization analyzers is given bŷ Finally, we can use this input-output realtion for the underwater channel, and the operatorâ in at the transmitter,â whereĉ i A(B) is the operator of environment modes, and the fluctuating transmission coefficient T is assumed to be equal for all four modes. Consider the parametric down conversion source which irradiates the states where χ is the squeezing parameter, and When n = 1, is the Bell state, where (θ (1) A , θ (1) B , θ (2) A , θ (2) B ) = 0, π 8 , π 4 , 3π 8 maximally violates the Bell inequality. By using Glauber-Sudarshan representation, the corresponding characteristic function for the state in equation (C9) can be expressed as By using the equations (C6) and (C7) and the quantum-state input-output relations Φ out (β) = 1 0 P(T)Φ in (Tβ)dT, we can get the characteristic function at the outputs of the polarization analyzers as Based on the characteristic function of equation (C12), the equation (C4) can be expressed as where d 8 β = d 2 β T A d 2 β R A d 2 β T B d 2 β R B , and K 0 (β) = 1 π exp − |β| 2 η − N , K C (β) = β(β) − K 0 (β). Based on the equations (C12) and (B13), the probabilities for simultaneous clicks at both sides and be obtained as where · denotes averaging with respect to T, C 0 = η 2 T 4 tanh 2 χ − 1 + (ηT 2 − 1)tanh 2 χ where p(r; d, σ) is the Rice distribution. The above equations (C14) and (C15) can be used in equations (C2) and (C3) to compute the Bell parameters. One can get the plots for the Bell parameter as a function of the squeezing parameter.
When considering the double-click events, the equation (C4) should be renewed (see reference [17] for more details) and the probabilities for simultaneous clicks at both sides will be changed to where i = {same, different} denotes the same and different transmitted, T A(B) , and reflected, R A(B) modes for the inputs of the polarization-analyzer detectors, respectively, and i = j. Here C same = C T A ,T B = C R A ,R B , and C different = C T A ,R B = C R A ,T B are given in equation (C15).

Appendix D. The calculation of violation of Bell inequalities for counterpropagating
The Bell test scheme for the scenario of counterpropagation is shown in figure A4, where the two links with transmission coefficients T A and T B are completely uncorrected. So the PDTC here can be written as where P A (T A ) and P B (T B ) are the single-mode PDTC for the corresponding underwater channels. Without loss of generality, here we consider a particular case similar as the experimental scheme implemented in reference [2], where only one mode A is sent through the scattering channel while the other mode B is operated in the source station. Thus we can get where T 0 is the deterministic transmission coefficient of channel B. For simplicity, we set T 0 = 1 in this paper.
When discarding the double-click events, the Bell parameter for the counterpropagating can be also calculated by the equation (C2), but the probabilities for simultaneous clicks at both sides in equation (C14) should be renewed as  . The Bell test scheme for the counterpropagation scenario derived from reference [17], where the atmospheric channel is here replaced as an underwater one. HWP denotes half-wave plate, which changes the polarization angles A and B; PBS denotes Polarizing beam splitter; D T A , D T B are the on/off detectors for the transmitted light, and D R A , D R B are the ones for the reflected light.