Monte Carlo-Based Performance Analysis for Underwater Continuous-Variable Quantum Key Distribution

: There is a growing interest in the security of underwater communication with the increasing demand for undersea exploration. In view of the complex composition and special optical properties of seawater, this paper deals with a performance analysis for continuous-variable quantum key distribution (CVQKD) over an underwater link. In particular, we focus on analyzing the channel transmittance and detection efﬁciency based on Monte Carlo simulation for different water types, link distances and transceiver parameters. A comparison between the transmittance obtained by simple Beer’s law and Monte Carlo simulation reveals that the transmittance of underwater link may be severely underestimated in the previous underwater CVQKD research. The effect of the receiver aperture and ﬁeld of view (FOV) on detection efﬁciency under different water types is further evaluated based on Monte Carlo. Simulation results show that the transmission distance of the underwater CVQKD system obtained by Monte Carlo simulation in pure sea water, clear ocean water and coastal ocean water is larger than that obtained by Beer’s law, while the key rate of the system in all types of water is smaller than that obtained by Beer’s law because the size and FOV of the receiver aperture are taken into account. By considering the practical system parameters, this paper establishes a comprehensive model for evaluating the security of underwater CVQKD systems with different system conﬁgurations.


Introduction
Underwater communication [1] has attracted much attention in recent years driven by the ever-increasing demand for ocean exploration and investigation. Compared to other communication methods, underwater wireless optical communication has been demonstrated to be of highest transmission data rate and lowest link delay [2,3]. Although most of the underwater optical links are implemented in line-of-sight configuration and more difficult to be eavesdropped, the potential security vulnerabilities still exists [4]. Fortunately, quantum communication provides a reliable way to achieve unconditional secure communication [5,6], and the feasibility of quantum communication over seawater has been demonstrated both theoretically and experimentally [7][8][9]. A famous application of quantum communication is quantum key distribution (QKD), in which any eavesdropping attempts by a eavesdropper Eve will unavoidably bring disturbance and be captured by the legitimated parties, Alice and Bob [10].
According to different implementation methods, QKD can be divided into discrete-variable (DV) QKD [11,12] and continuous-variable (CV) QKD [13][14][15][16]. Compared to DVQKD, CVQKD is easier to integrate with the standard off-the-shelf optical communication components and has higher secret key rates [17][18][19]. Recently, some analyses on underwater free space CVQKD have presented by researchers, such as channel parameter estimation for satellite-to-submarine links [20], security simulation for different kinds of seawater [21], and performance improvement via photon substraction [22]. All these studies used the well-known Beer's law to model the attenuation effects of the seawater due to its simplicity. However, it has two implicit assumptions that the transceiver setup is perfectly aligned and all scattered photons are lost, which may severely underestimate the received optical power since some of the scattered photons can still reach the receiver after being scattered multiple times and contribute to the final received power [2,23].
In this paper, we propose an underwater CVQKD performance evaluation method by using Monte Carlo simulation, which simulates the trajectories of the emitted photons in water. Specifically, we consider two main areas: the transmittance of the channel in different water types and the detection efficiency determined by different receiver aperture and FOV. In addition, the light source parameters, such as beam width and maximum beam divergence, are also taken into consideration to evaluate the secret key rate of the underwater CVQKD system. Simulation results show that the secret key rate of the system is obviously different from that calculated by simple Beer's law. By comprehensively considering the practical system parameters, this Monte Carlo-based model provides a feasible way to estimate the performance of various CVQKD systems with different system configurations. Figure 1 shows the schematic of a Gaussian modulated coherent state (GMCS) CVQKD system over an underwater link, which was proposed by Grosshans and Grangier in 2002 [13] and has been proven to be secure against general collective attacks and coherent attacks [24][25][26]. In this system, the sender Alice separates coherent light pulses generated by a laser diode into weak signal pulses and strong local oscillator (LO) pulses. Then the signal pulses are loaded with key information by using an amplitude modulator and a phase modulator, so that the X and P quadratures of the coherent state respectively follow a centred Gaussian distribution with a variance of V A . The LO pulses are delayed by a delay line so as to be transmitted together with the signal pulses through time multiplexing. At the receiving end, Bob first separate the LO and signal pulses by using a beam splitter and then randomly measures one of the quadratures of the signal states by homodyne detection, with the help of the LO as a phase reference. The phase modulator in Bob is used to randomly select the quadrature by imposing a π/2 phase shift to the LO. During this process, the pulses propagating in seawater are affected by two main factors: absorption and scattering. Absorption depends on the water's index of refraction and leads to energy loss of the light, while scattering changes the direction of photons and leads to broadening of a pulse [27]. Both absorption and scattering are concerned with the turbidity of the water and the type of particles in water. The total effects of these two process can be described by beam extinction coefficient c, which is defined as c(λ) = a(λ) + b(λ), where λ is the wavelength of the light, a(λ) is the absorption coefficient and b(λ) is the scattering coefficient. As a common practice, the chlorophyll concentration C (in mg.m −3 ) is employed as the free parameter to compute the values of a(λ) and b(λ), based on bio-optical models proposed in [28,29].

Underwater Channel Modeling Based on Monte Carlo
The behaviour of light propagation in a medium is described by the radiative transfer equation (RTE), which can be written as [27] dL(z, θ, φ, λ) where L(z, θ, φ, λ) denotes the light radiance at a point in units of Wm −2 sr −1 nm −1 , z is the distance of the point from the transmitter and r = z/cosθ represents the radial distance to the source, θ and φ are the polar and azimuthal angles, respectively. The first term in Equation (1) accounts for the annihilation of photons from the beam, the second and third terms denote the elastic scattering and inelastic scattering process, respectively. The effect of inelastic scattering is generally neglected because of its relatively low contribution to the solution of the RTE. Many previous works also neglect the contribution of elastic scattering and consider a straight-line propagation, describing the RTE by the simple Beer's law, as L(z) = L(0)exp(−cz).
Nevertheless, the received optical power may be severely underestimated in this way because it assumes that all the photons undergoing scattering are lost. In practice, the scattered photons may still be received, especially when the water has a high scattering coefficient. Therefore, in this paper we take the scattered photons into account and use a Monte Carlo method to solve the RTE.
Monte Carlo simulation is a rigorous and flexible method for modeling the underwater light propagation and has been adopted in many underwater optical communication research [30][31][32][33]. The main function of Monte Carlo is to model a complex probability density function based on its known individual components, which provides a feasible way for studying the characteristics of the channel under different water types, link distances, link geometries and transceiver parameters. By tracing the fate of a large number of photons according to the probability of absorption and scattering, this method can obtain a solution of the RTE. In 2008 [34], the robustness of Monte Carlo for underwater channel modeling was validated in experiment, as the results of Monte Carlo shown reasonable agreement with the water-tank experimental results. In 2014 [35], W. Cox and J. Muth also validated the accuracy of Monte Carlo simulation by comparing the simulation results with the experimental data from [36,37]. The major process of Monte Carlo simulation for photon tracking is detailed in Appendix A. In our simulation, we consider a standard underwater point-to-point link which is performed in an ideal isotropic and homogeneous seawater without flowing and turbulence, as many underwater optical communication analyses [31,35,38]. The practical underwater optical channel could be more complex if we consider water turbulence caused by the variations of water refractive index, temperature, salinity and pressure. However, in deep sea the level of sanity and temperature variations are very small and the channel fading caused by water turbulence is negligible [39]. Similarly, there is no blockage caused by bubbles, fish, or large suspended particles in deep sea. The wavelength of the light is defined as λ = 520 nm because that the seawater has a relatively low attenuation property for light with a wavelength of 450 nm to 550 nm [40]. Monte Carlo simulation is essentially a statistical method and the results rely on calculating of a large number of photons, so we generate 10 6 photons in each experiment and repeat each experiment at least 10 3 times to obtain reliable results.

Transmittance Analysis
As defined in Figure 1, the coherent states are sent from Alice to Bob through a quantum channel, which is characterized by its transmission efficiency T and its excess noise ε. Therefore, the total channel added noise referred to the channel input can be expressed as a combination of the noise caused by loss and excess noise, given by χ line = (1 − T)/T + ε, where ε is related to Eve's attack and T is related to the extinction coefficient c. As introduced in Section 2.1, the value of c is determined by the chlorophyll concentration C, which is a parameter that can describe the type of the water medium. In this section, we investigate the transmittance T of the underwater links in four typical types of water, including pure sea water, clear ocean water, coastal ocean water, and turbid harbor water.
The coefficient values of a, b and c of these four water types are listed in Table 1 [29,41]. For each water types, the transmittance is collected in two ways: Monte Carlo simulation and Beer's law calculation, and the results are shown in Figure 2. We can find that for the case of pure sea water, the difference between the transmittance obtained by Monte Carlo simulation and Beer's law calculation is negligible, because there are few particles in the pure sea water and the probability of photons being scattered is very low. The trajectory of photons in pure sea water is approximately a straight line. However, for the cases of clear ocean water, coastal ocean water, and turbid harbor water, the results of Beer's law are much smaller than the results of Monte Carlo simulation, which is due to the scattering coefficient of these water types is higher than pure sea water, and a significant part of scattered photons may still be received by the receiver, while Beer's law considers all the scattered photons to be lost. Therefore, the more turbid seawater photons are scattered more times, and accordingly, the greater the difference between the transmittance estimated by Beer's law and Monte Carlo. Figure 1. Schematic diagram of an underwater CVQKD system. The sender Alice separates coherent light pulses generated by a laser diode into weak signal pulses and strong LO pulses. Then the signal pulses are loaded with key information by using an amplitude modulator and a phase modulator, and the LO pulses are delayed by a delay line so as to be transmitted together with the signal pulses through time multiplexing. After being transmitted through the seawater, the pulses reach Bob's telescope. Bob first separates the LO and signal pulses by using a beam splitter, and then randomly measures one of the quadratures of the signal states by homodyne detection, with the help of the LO as a phase reference. The phase modulator in Bob is used to randomly select the quadrature by imposing a π/2 phase shift to the LO. The dark spots between Alice and Bob indicate the water components such as phytoplankton, suspended sediments and detritus, and colored dissolved organic matter. AM: amplitude modulator. PM: phase modulator. Tele.: telescope. BS: beam splitter. PD: photodiode. LO: local oscillator.

Detection Efficiency Analysis
In the CVQKD system, a practical homodyne detector is characterized by an efficiency η and an electronic noise v el . Therefore, the noise introduced by the detector referred to Bob's input is given by χ hom = (1 + v el )/η − 1. However, in an underwater scenario, the photons arriving the receiver plane may not be all received by the detector because of the finite size of the receiver aperture and FOV. As depicted in Figure 3, a series of photons propagate along z-axis toward a receiver with an aperture radius r, although several scattered photons reach the receiving plane, they are still considered lost due to exceeding the size of the aperture. To take the effect of the receiver aperture size and FOV into consideration, we define the total detection efficiency of Bob as η t = η r η, where η r is defined by where w i denotes the weight of the photons arriving the receiver plane, w j denotes the weight of the photons successfully received, and N a and N r represent the number of the photons arriving and received, respectively. Figure 4a shows the value of η r calculated by Monte Carlo simulation versus different size of receiver aperture under four types of water. We can see that the value of η r is extremely affected by the type of water. For example, in the turbid harbor water, η r is very small, while in the pure sea water, the value of η r is approximately equal to 1 when the receiver aperture is larger than 2 inches. This is due to the fact that in turbid harbor water, photons are scattered many times before reaching the receiving plane, which causes a large displacement in their position on the x/y plane. Furthermore, when the size of the receiver aperture is smaller than 2 inches, the value of η r is sensitive to the size of the aperture. This is because that a larger aperture can capture a larger number of photons scattered out of the incident beam. Figure 4b shows the value of η r versus different receiver FOV under four types of water, and we notice that in pure sea water, clear ocean water and coastal ocean water, the value of η r is almost unaffected by the receiver FOV, while in turbid harbor water the value of η r varies greatly with FOV. In addition to water types, the effects of the beam width ω 0 and the maximum beam divergence θ 0,max are also considered. The intensity distributions of the light with different beam widths on the receiver after transmitting through clear ocean water are shown in Figure 5a-c. We can see that with the increase of the beam width, the center of the beam spot becomes weaker, because a large beam width makes the initial positions of the photons more dispersed. Figure 5d-f show the intensity distribution of the light with different values of θ 0,max on the receiver after transmitting through clear ocean water. We notice that with the increase of θ 0,max , the beam spot gets larger and the intensity of the light becomes weaker, which indicates that the received power of the light is affected by the maximum initial divergence angle. Therefore, the detection efficiency should be calculated by integrated considering the water types and the characteristics of the light source.
x y x y z  z  x  y  0,max / 2 FOV r Figure 3. Schematic diagram showing the geometry of the transceiver and the moving trajectory of photons in seawater. θ x , θ y , θ z are the angles between the direction vector of the photon and the x, y, and z-axis, respectively. θ 0,max represents the maximum initial divergence angle of the beam and FOV represents the field of view of the receiver.

Performance Analysis for Underwater CVQKD
In this section, we analyze the performance of the underwater CVQKD system based on the parameters discussed above. The theoretical asymptotic secret key rate under collective attacks is obtained by [42,43] where β is the reverse reconciliation efficiency, I AB is the Shannon mutual information between Alice and Bob, and χ BE is the Holevo quantity which bounds the maximum information available to Eve. The calculations about I AB and χ BE are the same as fiber-based CVQKD system because in this paper we only consider a deep sea environment without turbulent effects. Therefore, the mutual information I AB is calculated by is the variance of the quadratures measured by Bob, and V B|A is the conditional variance. The Holevo quantity χ BE is given by where p(m B ) is the probability distribution of Bob's measurement, ρ m B E is Eve's state conditional on Bob's measurement m B , and S(ρ E ) is the Von Neumann entropy of the quantum state ρ E . In the case of Gaussian attack, the above equation can be further simplified as follows: where G(x) = (x + 1) log 2 (x + 1) − x log 2 (x). λ 1,2 are the symplectic eigenvalues given by with λ 3,4 are the symplectic eigenvalues given by with The last symplectic eigenvalue λ 5 = 1. In our simulation, the fixed parameters mentioned above are set to V A = 20, β = 0.95, ε = 0.01, η = 0.6, and v el = 0.01, corresponding to the standard realistic assumption for CVQKD implementations [44]. Figure 6a compares the secret key rate for different types of water with parameters calculated by Beer's law and Monte Carlo. Because Beer's law case only considers the influence of transmittance and does not consider other parameters. To better compare the two methods, in Monte Carlo case we analyze two situations. One considers the impact of the receiver aperture and FOV, and the other does not consider the impact of them with η r = 1. We find that when η r = 1, the key rate and transmission distance of the Monte Carlo method are higher than that of the Beer's law method, while when considering the influence of receiver aperture and FOV, the key rate and transmission distance are significant reduced. In particular, in pure sea water, clear ocean water and coastal ocean water, the secret key rate calculated by Monte Carlo method with considering receiver parameters is slightly smaller than Beer's law method, but the transmission distance is still longer than that of the Beer's law method. However, the results are very different in turbid harbor water. We can see that the key rate of harbor water calculated by our method is apparently lower than that obtained by Beer's law. This is because the value of η r in turbid coastal water is so small that almost all of the photons arriving the reception plane can not be successfully detected. Figure 6b compares the secret key rate for different beam width ω 0 and maximum beam divergence θ 0,max when the receiver aperture and FOV are set to 2 inches and 180 • , respectively. We find that the secret key rate and transmission distance are less affected by the beam width but are greatly affected by the maximum initial divergence angle, which corresponds to the results of the light intensity distribution in Figure 5. A large maximum beam divergence makes the beam reaching the reception plane more divergent, thereby resulting in a decreased value of η r .

Discussion and Conclusions
Our analysis shows that the there is a significant difference between the performance evaluated by simple Beer's law and Monte Carlo simulation, as the latter comprehensively considers the contribution of scattered photons and the system transceiver parameters. Monte Carlo simulation is a flexible and effective approach for studying the characteristics of underwater channels. Nevertheless, there is still a discrepancy between the results obtained by Monte Carlo and the actual implementation. Firstly, in this paper we only consider a perfectly aligned link performed in deep sea without transmitter-receiver misalignment, water turbulence, and sea surface reflection, which may affect the performance of the system. Secondly, despite the high flexibility and accurate solution, the Monte Carlo method also suffers from random statistical errors and low simulation efficiency. It is only a simulation which cannot be completely equal to real performance. Therefore, this discrepancy is worthy of further investigation and demonstration with experiments.
In conclusion, we proposed a realistic model to evaluate the performance of a CVQKD system over an underwater link by using Monte Carlo simulation, which took different water types, link distances, and transceiver characteristics into consideration to estimate the transmittance and detection efficiency of the system. Simulation results shown that the secret key rate and transmission distance of the underwater CVQKD system obtained by Monte Carlo simulation are distinctly different from those obtained by simple Beer's law. This research is beneficial for building a realistic underwater CVQKD system to achieve the necessary secure distance in different oceanic types, and also enables system designer to select the appropriate light source and transceiver devices for practical underwater CVQKD systems.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Major Processes in Monte Carlo Simulation
In this Appendix, we introduce the three major processes of Monte Carlo simulation for photon tracking, including initialization, photon propagation and photon reception, which are based on the description in [27,45].
(I) Initialization. As depict in Figure 3, the simulation geometry is defined by a cartesian coordinate system where a set of photons are assumed to be launched on x/y plane and propagate along the z-axis. Each photon is determined by its location (x, y, z) and direction µ x , µ y , µ z , where µ x = cos θ x , µ y = cos θ y , µ z = cos θ z . θ x , θ y , θ z are the angles between the direction vector of the photon and the x, y, and z-axes, respectively. The initial weight of each photon is defined as 1 to record the power loss of them. The initial position of each photon is generated according to U [0, ω 0 ], where ω 0 represents the beam width and U [a, b] represents the uniform distribution between a and b. The initial direction of each photon is generated according to U [−θ 0,max , θ 0,max ] for θ 0 and U [0, 2π] for φ 0 , where θ 0,max represents the maximum initial divergence angle, θ 0 denotes the initial zenith angle and φ 0 denotes the initial azimuthal angle. Therefore the starting direction of a photon can be calculated by µ x = sin θ 0 cos φ 0 , µ y = sin θ 0 sin φ 0 , µ z = cos θ 0 .
(A1) (II) Photon Propagation. The transmitted photon travels a random distance of δ before interacting with a particle in water. The value of δ for each photon is selected according to δ = −log(χ δ )/c, where χ δ is a random variable of distribution U [0, 1]. Besides, a part of photon weight is lost when interacting with the particle, and the new weight W post is calculated by W post = W pre (1 − a/c). After interacting with the particle, the propagation direction of photon is changed and the new direction is determined by regenerating the zenith angle θ and azimuthal angle φ. The regenerated angle φ is considered to be a random variable of distribution U [0, 2π], and the regenerated angle θ is selected based on Henyey-Greenstein model, which is a commonly used scattering phase function to model the trajectory of the scattered photon in seawater [27,35].
(III) Photon Reception. All the steps mentioned in II are cyclically repeated until the photon reaches the receiver plane or the photon weight is too small and negligible. The receiver is defined on the x/y plane and perpendicular to the beam axis. Please note that only the arrived photons within the receiver aperture and with incident angles less than the half angle of receiver FOV are considered to be received, otherwise they are considered to be lost.
After the above processes, we can obtain the weights of the photons that reach the receiving plane and the photons that are successfully received, as well as the Cartesian coordinates and incident angles of these photons when intersecting with the receiving plane.