New Monte Carlo Integration Models for Underwater Wireless Optical Communication

Underwater wireless optical communications are knowing a high interest in recent years, with many military and civil applications. The design of efficient systems relies on powerful simulation tools. The most used simulation algorithms are based on a paper from 1989, written by Prahl and designed for light propagation through human tissue. This method relies on Monte Carlo Simulation, following the path of a photon in the participating media. As such, it is difficult to propose evolution or optimization. In addition, a large amount of photons are needed to obtain good results with low error, especially for turbid water. With Monte Carlo Integration method, there exist a lot of optimization techniques to reduce the variance. This paper proposes a mathematical formalization of the propagation of light in water or any other participating medium as an integral problem. Therefore, it opens the way for a large number of future optimizations. A very straightforward and simple variance reduction technique is proposed as an example. Our simulation results show that this new technique has a lower sample variance as expected, and thus better convergence rates.


I. INTRODUCTION
Underwater wireless communications are present in many different applications in both military and civil domains, like submarine communications, fishing industry, ecology, etc. [3]. There are three main technologies for this purpose: acoustic wireless communications for long distances but with high latency and low throughput; wireless radio communications, but with low data rates due to high signal absorption by water; and Underwater Wireless Optical Communications (UWOC), a promising technology with good properties such as high data rates, security due to small propagation range, and large bandwidth. Over the past few decades, much literature has been produced on UWOC, as evidenced by the many The associate editor coordinating the review of this manuscript and approving it for publication was Giambattista Gruosso . recent surveys dealing with applications of UWOC [1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12].
The development and characterization of UWOC systems require accurate and fast simulation tools, to test the developed communication protocols on different scenarios. We need to well understand the propagation mechanisms involved in UWOC channels to provide such tools. This kind of channel is characterized by different physical phenomena, inherent to ''participating media'': light absorption, due to particles that typically absorb and convert the light flux into heat; light scattering, due to particles that deflect the luminous flux in other directions, causing the dispersion of the light. The latter is particularly true for participating media like fog, clouds, and liquids such as good old Cognac or turbid harbor water.
Absorption and scattering imply that the UWOC channel can become very difficult to simulate. Previous works VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ have proposed mainly two kinds of UWOC channel models: models based on Monte Carlo Simulations, where photons are tracked in the scene and counted when they reach the receivers; and statistical models fitting either measurements or Monte Carlo simulations. This article considers only the first category.
To the best of our knowledge, all previous solutions based on Monte Carlo historically came from the neutron transport problem, and thus from the photon transport. Using the formalization of radiative transfer, their foundation is the differential form of the Radiative Transfer Equation (RTE) [16]. Previous works have considered the propagation of a single photon or a group of photons, usually using a Monte Carlo Simulation performed with this differential form of the RTE. Even more recent works that attempt to formalize propagation of light flux by an integral equation continue to use this approach. But there is a dual form of the RTE, given in integral form and which expresses the radiance received at a given point by considering that emitted at another point [16]. In this article, we propose to start from scratch with this integral form of the RTE. This will allow to express the UWOC channel in many different ways with different integral equations, onto which it is possible to apply various variance reduction techniques.
To summarize, the objectives of this paper are: • To provide a robust mathematical approach to the UWOC channel modeling as an integration problem.
• To reveal the link with Monte Carlo algorithms existing in previous works, mainly the Prahl's one.
• To propose a new algorithm by removing useless importance sampling in previous approaches.
• To show the variance reduction obtained with this new algorithm through some simulation results. The rest of this article is organized as follows. Section II recalls the previous works on UWOC channel simulation, the main algorithm proposed by Prahl and its properties. Section III derives a mathematical formalization of the UWOC channel based on the integral form of the RTE, retrieving the Prahl's algorithm. Section IV introduces a new algorithm thanks to our new mathematical formalization, with simulation results in comparison to previous methods. At last, Section VI concludes this article, with some future works.
All the notations used in this article are summarized in Table 8 in Appendix B.

II. PREVIOUS WORKS A. ON THE RADIATIVE TRANSPORT EQUATIONS
Propagation of light in in-homogeneous or homogeneous media is an old topic. Schuster investigates the absorption and scattering of light in such media, considering half-forward and backward scattering only [13]. Later, Preisendorfer proposes mathematical foundations to the light scattering process [14]. In the sixties, Chandrasekhar and Preisendorfer produced two books of reference about radiative transfer [15], [16], that are considered as the basis of all modern works albeit the first one is rather difficult to read, focusing on calculation procedures. The second recaps all the main physical quantities allowing to express the propagation of light in participating media, from radiant flux to radiance L(x, − → ω ) at point x and in unit direction − → ω and their properties. It also defines the ''transmittance'' as the ratio of outgoing to incoming radiance for an infinitesimal volume of length l: It states the transmittance's properties: multiplication property T (l 1 + l 2 ) = T (l 1 )T (l 2 ), contraction property T (l) ≤ 1 and identity property T (0) = 1. It also discusses the ''volume attenuation function'', as the light attenuation effect linked to T (l): the loss of radiance being 1 − T (l), the attenuation is defined as: Today, σ t is called the extinction coefficient. It can be defined as the sum of the absorption and the scattering coefficients σ t = σ a + σ s . Transmittance being dimensionless, the dimension paired with these coefficients is that of inverse length. The transmittance is defined as follows for homogeneous participating medium: The scattering of light by small particles in the participating medium is introduced thanks to the Volume Scattering Function (VSF) β. Still in [16], Preisendorfer introduced in Section 21 the two forms of the RTE. Ignoring the emitted radiance by participating media (which is always the case for UWOC) and with today's notation, its expression in differential form is as follows: where − → ω is a unit direction, x is the position where radiance is expressed, and L i is the incoming radiance expressed as follows: where is the space of directions in R 3 , corresponding to unit sphere points. The simplified integral form of the RTE is: where x s = x + s − → ω is the point at distance s in direction − → ω , and x l = x + l − → ω are points on the half-line starting at x and with direction − → ω .

B. FORMER WORKS ON LIGHT PROPAGATION SIMULATION IN GENERAL PARTICIPATING MEDIUM
Solving directly these equations is quite impossible in most practical situations. Twomey et al. used the matrix methods for the computation of the reflection and diffuse transmission of radiation through multiply-scattering parallel layers [17], [18]. To our knowledge, Collins and Wells were the first to solve the RTE using Monte Carlo Simulation [19], [20] using previous code developed for neutron transport problems. They tried to calculate the propagation of light after a nuclear explosion, taking into account the different layers of the atmosphere, the sea surface and clouds. Later, Plass and Kattawar applied Monte Carlo calculations to light propagation in the atmosphere, including clouds, and then from atmosphere to the sea including reflection on the sea bed, considering uncollimated light source and laser, visible and infra-red light [21], [22], [23], [24], [25], [26], [27], [28]. It is remarkable that all these previous works tried to connect every scattering point with each receiver considered as point (surface receiver being handled in another way). This technique allows a significant reduction of the variance. It is generally called ''Next Event Estimation'' (NEE). Moreover, they made use of the differential form of the RTE, and so naturally they used importance sampling in the choice of next scattering points considering the transmittance as Probability Density Function (PDF). They also used importance sampling for the choice of the new scattering direction considering the phase function as PDF, albeit they had employed tabulated form of the VSF and so the cumulative distribution to do the sampling. All the roots of modern calculations where already set in these first works on Monte Carlo Simulation of light. Gordon et al. used a similar Monte Carlo algorithm to evaluate the light reflected by the sea [29], taking into account the sea bed [30] and two layers of ocean [33]. But contrary to previous solutions, they used in [31] a phase function obtained from measurement in the ''Sargasso sea'', and ''Tongue of the ocean'' (Bahamas island) from Petzold's measurement [32].
Poole et al. have proposed later to use a ''semianalytic'' Monte Carlo method to calculate the propagation of light in ocean, and the flux received by an atmospheric lidar system [34]. They applied NEE to reduce the variance due to the low probability of a photon crossing the receiver surface, by explicitly connecting to it each scattering point. This corresponds in fact to the Collins' method, and the Plass and Kattawar for point receivers, but here for surface receivers. Moreover, they explicitly used truncated distributions to ensure the photon stays in the ocean, and does not exit to the atmosphere.
Kirk used a Monte Carlo procedure similar to the Plass and Kattawar one in [35] and [36] with a phase function corresponding to very turbid water obtained from data measured in San Diego's harbor by Petzold [32]. From his simulation results, he deduced an analytical simple law for light propagation in this kind of sea water. The same Monte Carlo method was applied by Lermer and Summers using the Henyey-Greenstein phase function [37] for the first time to our knowledge. Furthermore, they also consider the angle and time of arrival of photons at the receiver to draw some statistics, showing that previous simple analytical propagation models were no more valid for such conditions.
The modern algorithm used in UWOC channel estimation was first proposed by Prahl in [38], albeit it was for collimated light propagation in layered tissue. They consider that Monte Carlo samples are groups of photons, and not individual photons. Notice that they used Russian roulette to stop the propagation of the photon packet, and not a threshold as in former works. They also provided their algorithm through an organigram that may explains why this work was used later as the starting point, albeit it resembles a lot to all the previous solutions starting from Collins's one.

C. PREVIOUS UWOC METHODS
Recent years have seen a lot of works for UWOC channel estimation. Jaruwatanadilok has reformulated the differential form of the RTE first, and then solved the incoming radiance L i thanks to a Fourier transformation, obtaining a new differential integral expression [39]. This new expression was then solved using Monte Carlo Integration process. However, this method first involves a discretization of the space of incoming directions to evaluate L i , which then requires balancing accuracy and computation time. Notice that the used phase function was obtained from Mie scattering.
More recently, Ghazy et al. investigate the misalignment problem in UWOC system using MIMO link [66], with simulations made by using a single scattering assumption. But more than one scattering is needed for realistic water types, as shown in Section V. For instance, clear water needs 22 scatterings for a 12 meters link using Prahl's algorithm and 1E −6 as threshold value. In the same scenario but with coastal or harbor waters, 40 and 94 scatterings are needed respectively.
Some works try to solve the RTE using finite difference scheme [67], [68]. Based on previous time dependent RTE (TD-RTE) solution [69], Illi et al. simulate the light propagation in two dimensions, for rectangular cross-section finite propagation domain. The height of their studied space is 0.2 meters only, for distance from transmitter to receiver from 1 to 22 meters. They show that in such a configuration, their finite difference scheme becomes faster than a Monte Carlo simulation algorithm. Nevertheless, the studied VOLUME 10, 2022 domain is too small for many water types and in ocean, since many significant scattered contributions come from outside. Moreover their finite difference scheme needs to subdivide the propagation domain in both geometrical space, direction space and time. For an acceptable level of accuracy and large 3-dimensional propagation domains including obstacles like in [58], the complexity of the TD-RTE method will significantly increase, and its memory requirements will become prohibitive. On the contrary, the Monte Carlo variance is related to the square root of the number of samples, regardless of the problem's dimension. Thus, in realistic configurations, Monte Carlo simulations will become the best solutions for the same level of accuracy.
Recent works explore the possibility to express the propagation of photons as a probabilistic mechanism [68], [70], [71], [72], [73]. In [72] and [73] the single scattering effect was investigated first, before to generalize to multiple scattering and to apply Monte Carlo Integration techniques. In [69], [71], and [74], multiple scattering are considered, in 2D in [68] or in 3D for others. The logic behind these articles should be noted: starting from the differential form of the RTE, the application of Monte Carlo Simulation leads to probabilities for the step size of photon's movement and the angles of deflection in a photon's trajectory, and then a combination of these probabilities leads to a sort of integral form of the RTE describing the photon arrival probability after n scattering events. However, the integral form of the RTE was known for a long time, and it is far simpler to start from it instead of through these circonvolutions. That's exactly our purpose: to restart from the integral form of the RTE and to express the light propagation with 0 to n scatterings, then to retrieve the classical Prahl's algorithm and to apply some conventional variance reduction techniques resulting in a solution close to the semianalytical old method, albeit without assuming a locally constant VSF. Notice also that our method allows using any kind of VSF.

D. MONTE CARLO SIMULATION AND INTEGRATION
Most of the previous works relies on Monte Carlo Simulation, which consists in simulating a physical process by using some parts of its analytical description as a sampling function. Here, sampling function means using a PDF and a sampling mechanism based on such a density. Such a PDF p(x) is a positive function with unit integral defined over its definition domain D: The main two physical events in UWOC channel are the scattering and the extinction, respectively modeled by the VSF β and the transmittance T . These two events are directly used into previous simulation algorithms. Unfortunately the Monte Carlo Simulation does not provide a complete mathematical expression of the radiance at the receiver after n scattering events. We think this explains why, since Prahl's Algorithm 1 Prahl's Algorithm [38] 1: x 0 = samplePoint(T x ); 2: − → ω 0 = sampleDirection(T x ); 3: ray = createFromStartAndDir(x 0 , − → ω 0 ); 4: factor = 1 / p(x 0 )p( − → ω 0 ) ; 5: length = 0; 6: for n = 1 to maxScatteringPerLink do 7: // Sample the length 8: x n = setEndOfRayWithTransmittance(ray, σ t ); 9: if isCrossingAPlane(ray, endSurfaceSensor) then 10: crosses R x ? 11: tryReceiveLineSegment(ray, factor, length); 12: return null; 13: end if 14: if useVisibility then 15: hit = sceneIntersection(ray); 16: if hit = null then 17: return hit; 18: end if 19: end if 20: length = length + getDistance(ray); 21: factor = factor * σ s / σ t ; 22: // Sample the direction 23: newDir = sampleDirection(ray, σ s ); 24: setDirectionAndStart(ray, newDir, ray.to); 25: end for 26: return null; 1989 algorithm, there have been no significant changes in the way participating environments are simulated.
On the contrary, with Monte Carlo Integration method, we start from a mathematical expression of a given problem. This allows to use many different techniques to reduce the variance of the estimator, and so the number of samples for a given error threshold [76]. One of these techniques is importance sampling, which is inherent to Monte Carlo Simulation. With Monte Carlo Integration method the best variance reduction technique can be chosen, without being limited to the only importance sampling method.
Before deriving such a mathematical formalization in Section III, the following section proposes a general algorithm based on previous works, mainly those of Dong et al. [47] and Prahl [38].

E. A DETAILED IMPLEMENTATION
This section presents a detailed algorithm for simulating the UWOC channel, based on previous works (see algorithm 1). More precisely it is based on [47], itself based on [38]. Hence, this algorithm is named as Prahl's algorithm in this article.
Lines 1 to 5 consist in some initialization, for a given path. Then, the loop from lines 6 to 25 builds a path. Actually this algorithm evaluates a single path, thus it should be called N times for a full simulation.
This for-loop starts by setting the end of the ray at line 8, thanks to a PDF p(l) = σ t T (l). Then, lines 9 to 13 manage the ray that passes behind the receiver. In such a case, it is considered that the path cannot cross again the receiver, and so the loop exits (line 12). Before to exit, the intersection between the segment of line [x n−1 , x n ] and the receiver R x is handled at line 11: the function tryReceiveLineSegment does this test and if necessary, adds the path's weight factor to a global accumulator at specific time corresponding to the path's length.
On lines 14 to 19, the intersection of the path is tested if necessary (line 14) with the obstacles in the scene (line 15). If such an intersection occurs, the intersected object is returned (line 17). Notice that this allows to handle some objects in the participating media, like some authors are doing [58].
Then the next segment of path is built from lines 20 to 24. Line 20 extends the path length to the last point x n . Line 21 updates the path's weight to x n . Line 23 samples a new direction thanks to the VSF β (i.e. the PDF is β). Line 24 resets the ray so that it starts from the previous endpoint x n and has the direction newDir sampled in line 23.

III. MATHEMATICAL FORMALIZATION
This section aims at providing a mathematical formalization of the light propagation in a participating medium as an integral problem. It will help to rewrite existing UWOC algorithms based on Monte Carlo Simulation with a Monte Carlo Integration scheme. More specifically, the goal is to provide a formalization of the Prahl's algorithm given in Section II-E. The new formalization is a key point to think about new algorithms with lower variances and thus lower computation time.

A. VOLUME LIGHT PROPAGATION
Considering only scattering events in participating media, the light propagation channel can be modeled by its impulse response as follows: where n is the number of scattering events [75] and P t is the total radiated flux from T x . Notice that the impulse response's unit is s −1 . The formalization of P n (T x , R x , t), the received flux at R x coming from T x after n scattering events at time t, can be written thanks to the integral form of the RTE (see chapter 15.1 in [74]) as follows: where A T x is the transmitter's surface, V is the volume of the participating media, A R x is the receiver's surface, y = {x 0 , x 1 , . . . , x n , x r } is a propagation path of light starting at x 0 on T x and reaching x r on R x after n scattering events at locations x i for i ∈ [1 . . . n], and f n 0 (y) is the contribution of the path y defined as: where: is the visibility function defined as: • − → n R x is the normal vector to receiver's surface.
• τ is the propagation delay of path y.
• δ(t − τ ) is the Dirac delta function. This integral can be calculated numerically using Monte Carlo Integration [75]. The corresponding estimator is: where y k is the random variable corresponding to a sampled path, generated according to the PDF p(y k ), N is the number of samples y k .
Here, p(y k ) depends on the sampling method. Considering that each point x i of y k is sampled independently, p is expressed as: This estimator can be evaluated by uniformly choosing points on their respective domains, i.e. the transmitter surface for x 0 , the receiver surface for x r , and the participating medium volume for other points x i . Unfortunately, the variance Var( P n ) of such an estimator is very large due to the high dynamics of transmittance and VSF. This high variance leads to a very poor estimator with high error: by the central limit theorem, the convergence of this method is in Var( P n )/N .
A conventional way to reduce the estimator variance and so the error for a given number of MC samples, consists in using importance sampling [75]. It uses parts of the integrated function (here f n 0 ) as PDF, such that these parts vanish in Equation (12). It is the basic idea involved in Monte Carlo Simulation algorithms, and so in previous UWOC simulation tools based onto Prahl's algorithm. It uses importance sampling to select each scattering points x i for i ∈ [1 . . . n] involved in sample y k , by considering PDF based on β and T . To provide a formalization of Prahl's algorithm, each volume integral (over the domain V ) must be transformed in a pair of integrals over the directions' space (for β) and the distances' space L (for T ). VOLUME 10, 2022 B. VOLUME TO DIRECTIONAL FORMALISM Each volume integration appearing in Equation (9) must be replaced by integration over two other domains: The directions' space , and the distances' space L. To this end, it can be noted that the choice of any volume point x i when constructing the sample path y k is done by knowing the previous path point x i−1 , and this from x 0 on the transmitter. Therefore, it is possible from a given point x i for i ≥ 1, to construct a bijection between the space of points V and that of directions and distances L. From Equation (44) in appendix A it follows: From this equality, Equation (9) can be rewritten as follows: where f n 1 (y) = f n 0 (y) Notice the product l 2 0 . . . l 2 n−1 that comes from transforming the integral domain V to × L. Hence, this equation can be simplified as follows: where each point x i is obtained from previous one as for i ∈ [1 . . . n].

C. A FURTHER STEP TOWARDS PRAHL's ALGORITHM
The last integration domain on the receiver's surface in Equation (15) can be replaced by an integral over directions R x . For this purpose, we use the following transformation coming from Equation (43) in Appendix A: It leads to the following expression of the flux received at R x coming from T x after n scattering events: (20) using: The function f n 2 can be written as the following: where and is the intersection point between the receiver surface A R x and the half-line starting at x n with direction − → ω n . The last integration domain on the receiver's solid angle R x can be replaced by full direction domain . Therefore, to ensure the exactness of such a new formulation, all the extra directions d − → ω n not being into the domain R x should be multiplied by 0. For this purpose, we introduce the function H x n , − → ω n defined as follows: To summarize our transformations, the following relations apply: Then, the flux received at R x coming from T x after n scattering events becomes: where

D. PRAHL's ALGORITHM FORMALIZATION
In his work Prahl uses Monte Carlo Simulation that relies on importance sampling. Actually he used importance sampling both for all phase and transmittance functions appearing in Equation (10), even for the last segment of path y from x n to x n+1 . This corresponds to n + 1 integrals on the doubledomain ×L. However, the Equation (25) contains n+1 integrals on but only n integrals on L. To obtain the integral corresponding to Prahl's algorithm, we need to introduce a last integration domain on L.
For this last segment, Prahl used importance sampling according to the transmittance to sample a length. This allows to find a point x n+1 onto the half-line starting at x n with direction − → ω n . As a consequence, Prahl has to check that the segment of line [x n x n+1 ] intersects the receiver. So we have to introduce a new integral in our formalization that corresponds to this behavior, while producing the same result for the estimation of the flux received at R x after n scattering events. Let us write it as follows: where the function H checks the intersection of the segment of line [x n , x n+1 ] with the receiver's surface: It should be noticed that H (x n , x n+1 ) = 1 implies H (x n , − → ω n ) = 1 (but not necessarily the reverse). Moreover, g can be any function that respects the Equation (27).
Putting all together, we obtain the following Prahl's algorithm formalization: where Then, we deduce the following Monte Carlo estimator: In his work, Prahl uses the following probability p(y k ): where: • p (x 0 ) depends on the emitter area, • p − → ω 0 depends on the radiation pattern of the emitter, where i ≥ 0 for the scattering events, • p (l i ) = σ t T (l i ) for extinction events. To end the formalization of Prahl's previous work, g(l n ) should be chosen to use importance sampling over the last integration domain L, according to p (l n ). It leads to: where c is a normalization constant introduced to respect the Equation (27). It can be computed as follows: Hence, the expression of g(l n ) is: The product g × G divided by the p(l n ) can be simplified as follows: Moreover, it should be observed that using importance sampling into Equation (29) leads to: At last, the MC estimator corresponding to Prahl's algorithm is simplified to the following expression: There exist different methods to reduce the variance in Monte Carlo methods [76], like stratified sampling, importance sampling and control variates. Since Prahl's algorithm is based on the importance sampling intrinsic to Monte Carlo Simulation, we focus in this paper only on importance sampling. Considering our first formalization of the UWOC propagation channel through the function f n 0 , there exist very few ways to apply importance sampling. Indeed, the integration domains are A T x , V and A R x . Then, we may apply importance for: • choosing starting point x 0 , • choosing scattering point x i in the volume, • choosing the reception point x r . In fact, none of the points on the path y can be chosen directly with importance from f n 0 . We wrote a first algorithm based onto this formalization with uniform sampling, and its results are very poor as expected. Actually, the more important parts of f n 0 are the VSF β, the transmittance T and the visibility function V is . The later is quite difficult to predict analytically, and since generally the UWOC channel does not include many obstacles if any, we mainly focus here on β and T . Hence, it is more interesting to choose the next point x i+1 from the point x i using β and T . The sampling spaces are for the former and L for the latter. So, we need at least the second formalization using f n 1 . But, we would also like to apply importance sampling to the last scattering point, de facto excluding the use of f n 1 . Moreover, applying importance sampling according to β on f n 2 would require to normalize β on R x , which seems to be very difficult in practice. Finally, Prahl having implicitly used the form f n 4 , we propose to define an algorithm based on f n 3 .

B. NEW UWOC ALGORITHM
There are very little differences between expressions f n 3 and f n 4 : the later just adds a final integration domain on L. So it needs an extra term g(l n ) that respects the Equation (27), plus a new version H of the hit function. From the algorithm point of view, this mainly implies that the contribution is added after sampling the length l n and when H = 1.
On the contrary, by using the expression f n 3 , the contribution should be added when the direction − → ω n is known, so one step before as depicted in the Algorithm 2, line 8. It is remarkable that the two algorithms differs only a few. Obviously, instead of checking the intersection with function H , we have to check with function H : The difference is that the former uses a segment of line [x n , x n+1 ] while the second uses a half-line [x n , − → ω n ).

Algorithm 2 New Algorithm
1: x 0 = samplePoint(T x ); 2: ω 0 = sampleDirection(T x ); 3: ray = createFromStartAndDir(x 0 , ω 0 ); 4: factor = 1 / (p(x 0 )p(ω 0 )); 5: length = 0; 6: for n = 1 to maxScatteringPerLink do 7: // Try to connect 8: tryReceiveHalfLine(ray, factor, length); 9: // Sample the length 10: setEndOfRayWithTransmittance(ray, σ t ); 11: if isCrossingAPlane(ray, endSurfaceSensor) then 12: return null; // Early stop 13: end if 14: if useVisibility then 15: hit = sceneIntersection(ray); 16: if hit = null then 17: return hit; 18: end if 19: end if 20: length = length + getDistance(ray); 21: factor = factor * σ s / σ t ; 22: // Sample the direction 23: newDir = sampleDirection(ray, σ s ); 24: setDirectionAndStart(ray, newDir, ray.to); 25: end for 26: return null; Then, the probability to add a contribution is higher with f n 3 than that with f n 4 . An interesting side-effect of this property is that while the Prahl's algorithm adds at most one contribution per path, our new algorithm potentially add up to n + 1 contributions for a single path. One may argue that each of these contributions is not sampled using the transmittance for the very last segment, and so intuitively the variance would increase as compared to the Prahl's one. But since even for clear water many scattering events are needed to reach the receiver, the higher number of contributions should compensate this aspect. This will be confirmed by the results shown in Section V.
It should be noticed that this new algorithm proposes an early exit from the loop (line 11 to 13), as does Prahl's algorithm. It corresponds to checking whether the new scattering point x n is behind the receiver plane. In such a case, we admit that the probability of returning in front of the receiver to recross it is very low.

V. PERFORMANCE EVALUATION AND DISCUSSIONS
This section explores simulation results obtained from our implementation of two UWOC impulse response estimation algorithms. The first one corresponds to the Prahl's algorithm described in Algorithm 1. The second one is our new version, described in Algorithm 2 and based on the Equation (25).
We start with a well-known configuration coming from the literature [47] to compare our new solution to Prahl's algorithm. But a single scene is not sufficient to draw general conclusions. Hence, we generate different configurations from the initial one by modifying one parameter at a time: R x 's area and FOV, distance between T x and R x , R x 's orientation and position relatively to T x .

A. METHODOLOGY
The initial configuration is based onto the literature [47] and is defined as follows.
• There is no obstacle, but only water. So, the volume of water is considered as infinite.
• The receiver R x is a disc with a default Field Of View (FOV) set to 90 • , with diameter of 20 cm, located at (12, 0, 0) and oriented to (−1, 0, 0).
• The distance between T x and R x is 12 meters by default.
• The default water is Harbor, with parameters σ a and σ s being equal to 0.295 and 1.875 respectively.
• The phase function β is Henyey-Greenstein with coefficient g equals to 0.9199.
• The maximum number of scatterings is deduced from the threshold set to 1E −6 . For Harbor water, it corresponds to 94 scattering events at most. Let us recall that in previous works, the path is terminated when its weight becomes lower than this threshold. This weight equals to σ s σ t n , where n is the number of scatterings. So, the maximum number of scatterings is the integer part of log T h log σ s −log σ t where T h denotes the threshold. Each simulation is made using 1 billion of Monte Carlo samples (paths), on a PC running with a bi-Xeon Silver 4210R at 2.40 GHz and 64 Gb of RAM, using 40 threads (the two algorithms being parallelized). These simulations produce discrete impulse responses with a time step t = 0.1 nano second. From these impulse responses, we can calculate the channel's gain C g defined as follows: The two algorithms are discussed comparing the computation time, the unbiased sample variance and the sample error. The unbiased sample variance is computed from the samples of the population, so from the paths used during the simulation. Let C k g be the contribution of the path y k to the channel's gain C g (T x , R x ), and let C g be the average of these contributions. The unbiased sample variance of C g (T x , R x ) is then as follows: From this sample variance, the estimation error of C g (T x , R x ) is computed as follows: The following sections present simulations for different configurations, first to compare impulse responses, computation time and unbiased sample variances considering the initial configuration, then by varying different parameters: R x 's area and FOV, distance between T x and R x , R x 's orientation and position relatively to T x .

B. PERFORMANCE COMPARISON WITH INITIAL CONFIGURATION
This section directly compares the two algorithms using the initial configuration presented above. Figure 1 shows the impulse responses obtained from the two algorithms with 1 billion samples.  Figure 2 shows the local variance C g (T x , R x , i) of the two simulations, corresponding to the unbiased sample variance of the channel's gain computed into each temporal bin of the impulse response, and given as follows.
Clearly, our new algorithm reduces the sample variance. It confirms what the Figure 1 seems to indicate: the new algorithm has a better convergence. Thus, we observe a reduction of both the estimator's variance and error for a small increase of computation time. This point is very important: if we want to reduce the error of the first algorithm to the one of the second, we need to launch more paths and so to increase the computation time accordingly. By denoting Var(C 1 g ) the sample variance and N 1 the number of paths for the first algorithm, and Var(C 2 g ) and N 2 the ones of the second algorithm, we can estimate this required number of paths as: Hence, we need 1.2479 billion paths for the first algorithm to obtain the same error as with the second algorithm. In such a case the computation time is 17 minutes and 3 seconds, to be compared to the 14 minutes and 6 seconds for the second algorithm with the same error. This shows that the new algorithm is faster than the previous one for a given target quality.
To enhance this result, Figure 3 shows the expected computation time for the two algorithms varying the Monte Carlo target error into the range [1E −9 . . . 9E −9 ]. Considering the minimum target error, the computation time is 20 hours and 24 minutes for the first algorithm and 17 hours and 25 minutes for the second. The corresponding impulse responses are shown on Figure 4. Notice that for the same low target error, the first algorithm needs 89 billion paths while the second algorithm needs 72 billion paths only. Considering the maximum target error used into Figure 3, the computation time is 15 minutes and 24 seconds for the first algorithm and 13 minutes and 32 seconds for the second. The corresponding impulse responses are shown on Figure 5. For high target   error, the first algorithm needs 1.1 billion paths while the second algorithm needs 884 million paths only.

C. CHANGING RECEIVER's AREA AND FOV
The effectiveness of the new algorithm may depend on different configuration's parameters. Among these parameters,   one is the area of the receiver. Indeed, the probability that a path y crosses the receiver depends on its size. Table 1 shows the impact of the receiver's area A R x on the unbiased sample variance. It can be observed that our new algorithm is more efficient with very small receiver. This probably comes from the fact that the probability to cross the receiver grows with the apparent size of the receiver from the transmitter and/or from any scattering point. Hence, the impact of our new algorithm is reduced with large receiver. Another parameter that may impact the variance of our simulation algorithm is the FOV. Table 2 shows the impact of the FOV on the variance of the Monte Carlo estimator with the Prahl's and our new algorithms. It can be observed a behavior similar to the area's impact until a 90 degrees FOV, and after a stabilization of this effect. In fact, the VSF β is so directive in the harbor water that only a small percentage of the radiated power is received beyond a 90 degrees FOV.

D. CHANGING DISTANCE BETWEEN TRANSMITTER AND RECEIVER
The distance between the transmitter and the receiver may affect the simulation algorithms' variance. Table 3 shows the results obtained by varying the distance from the transmitter to the receiver along the line starting at (0, 0, 0) and with direction (1, 0, 0). It should be noted the variation of the channel's gain done according to the distance. This explains  the variation on sample variance, which by definition is proportional to the channel's gain as stated by Equation (39). Let us recall here that these are sample variances, and so an estimation of the variance of the estimators. As such, they are noisy since evaluated thanks to a Monte Carlo process. Nevertheless, it seems that the ratio between the variance of Prahl's algorithm and the variance of our new algorithm increases almost linearly with the increase of the distance between the transmitter and the receiver.

E. CHANGING RECEIVER's ORIENTATION AND RELATIVE POSITION
Both the receiver's orientation and relative position may affect the simulation algorithms' variance. Table 4 shows the results obtained by modifying the receiver's orientation, from 0 • to 75 • . Obviously, the first effect is to reduce the received power, since the receiver's solid angle is reduced. Hence, the unbiased sample variances are reduced too. But, it should be noticed that the variance reduction of our new algorithm remains. Table 5 shows the results moving the receiver while keeping the same distance (12 meters). In other words, the receiver is rotated around the transmitter position by angle from 0 • to 75 • . Regardless of its position, the direction of the receiver always points towards the transmitter. As expected the channel's gain is reduced, since the transmitter radiates mainly in one direction (with 0 • orientation) and its beam is very small (10 • ). It is also worth noting the reduction in the ratio of the variances. Nevertheless, our new algorithm remains better in performance.

F. USING DIFFERENT WATER TYPES
The last studied parameter that may affect the variance of our new simulation algorithm is the water type. This article uses different parameters that can be found in the literature [43]. This sections studies the impact of water type, using Pure   Sea, Clear Ocean, Coastal and still Harbor. Their parameters are presented in Table 6. The right column indicates the maximum number of scatterings for the considered threshold. It can be observed that considering this parameter, these four water types are roughly linearly spaced from Pure Sea to Harbor. The gap between the Harbor and the others is quite huge, especially for the scattering coefficient, explaining the difference in channel's gain shown in Table 7. This difference is also illustrated by the Figures 6-8: most of the light flux is received in less than 1 ns for Pure Sea, Clear Ocean and Coastal waters, as compared to the 25 ns shown in Figure 1.

G. DISCUSSION
The results presented in the previous sections allow to propose some general conclusions about the effectiveness of our new algorithm. For most of the studied cases, it appears that for a given computation time we obtain better results in term of variance and so confidence to the result. The sole exception is when the receiver's area is becoming large. In this case, the probability of a ray hitting the receiver increases, reducing the difference between the two algorithms in terms of sample variance.

VI. CONCLUSION
This paper deals with the simulation of the underwater optical communication channel, which is necessary for the design of new efficient UWOC systems. The previous works are based on an algorithm stated in 1989 by Prahl, which itself is based on earlier work by Collins and Wells, Plass and Kattawar, etc. These former solutions are based on Monte Carlo simulation. The lack of a well-defined mathematical foundation using integrals, for example, did not allow for improvement of this old algorithm.
This paper therefore proposes a new mathematical formalization of the UWOC channel. This new formalism allows to retrieve the old Prahl's algorithm but using the Monte Carlo Integration, i.e. a numerical integration technique using random numbers. This new formalism can be used to apply well-known Monte Carlo variance reduction techniques, like importance sampling which is the basis of Monte Carlo simulation, but also stratified sampling, etc.
Our new formalization having given rise to different variations, we have also proposed a new algorithm to estimate the impulse response of the UWOC channel. We have performed simulations which show that this new algorithm allows for a given computation time to obtain a better result, i.e. a result with a lower sample variance and therefore a lower error. Alternatively and for a given target error, our new algorithm is faster than Prahl's in all tested situations.
This new algorithm shows that there is some space to improve the estimation of the UWOC channel, with new algorithms having even better performance. In future work, we plan to apply different Monte Carlo variance reduction techniques, such as importance sampling or next event estimation. In addition, we intend to merge these new algorithms with the classical Monte Carlo algorithm to handle both participating media with some refractions and reflections. This will allow the design of new wireless communication systems in complex environments including participating media but also surfaces, opaque or transparent objects, etc.

APPENDIX A SOLID ANGLE BETWEEN ELEMENTARY POINTS
The integrals manipulated into the optical channel modeling rely on three different kinds of domain or space: directions ∈ , surface points ∈ A, and volume points ∈ V . It is possible to move from one domain to another, using some simple transformations. The solid angle d − → ω x corresponding to an elementary surface point dx seen from a given position x is the most well known (see Fig. 9). By correcting the solid angle via a cosine to account for the tilt of the surface dx , it is expressed using the following two useful equations: The solid angle d − → ω x corresponding to an elementary volume point dx seen from position x is very similar to the surface case, except the cosine that disappears (see Fig. 9b). Indeed the elementary volume point has no surface, but just a volume. It leads to the following two useful equations: Table 8 summarizes the main notations used in this paper. It comes with units (SI) when necessary. VOLUME 10, 2022