Propagation and Depolarization of a Short Pulse of Light in Sea Water

We present the results of a theoretical study of underwater pulse propagation. The vector radiative transfer equation (VRTE) underlies our calculations of the main characteristics of the scattered light field in the pulse. Under the assumption of highly forward scattering in seawater, three separate equations for the basic modes are derived from the exact VRTE. These three equations are further solved both within the small-angle approximation and numerically. The equation for the intensity is analyzed for a power-law parametrization of the wings of the sea water phase function. The distribution of early arrival photons in the pulse, including the peak intensity, is calculated. Simple relations are also presented for the variance of the angular distribution of radiation, the effective duration of the signal and other parameters of the pulse. For linearly and circularly polarized pulses, the temporal profile of the degree of polarization is calculated for actual data on the scattering matrix elements. The degree of polarization is shown to be described by the self-similar dependence on some combination of the transport scattering coefficient, the temporal delay and the source-receiver distance. Our results are in agreement with experimental and Monte-Carlo simulation data. The conclusions of the paper offer a theoretical groundwork for application to underwater imaging, communication and remote sensing.

The significant attenuation and distortion of optical signals in sea water due to absorption and multiple scattering severely constrain the maximum operating range of laser-based systems. This places rather stringent requirements on a possibility to predict the optical channel characteristics, depending on actual sea water parameters. Of prime interest is information on the temporal stretching, spatial and angular broadening of a pulsed beam at the source-receiver distances up to several tens of the photon mean free paths in the medium.
Currently, numerical methods [3][4][5][6]9,11] (regarding earlier results see [15][16][17][18][19]) based primarily on Monte Carlo simulation are applied to computations of the light field distribution in sea water from non-stationary (pulsed, modulated) sources. These methods enable taking into account essentially all factors inherent in actual experimental conditions (the angular profile of the sea water phase function, the mutual source-receiver orientation, limited aperture and field-of-view of the receiver etc.). However, the results obtained in such a way, as a rule, are related to a particular set of parameters and can not be transfered directly to other cases. In this connection, there remains interest in finding analytical solutions to the radiative transfer equation that would be based on adequate approximations and satisfactory reproduce the light field distribution in sea water for actual values of the optical parameters.
In propagation of a short pulse through the highly forward scattering medium (1 − ⟨cos γ⟩ ≪ 1, where ⟨cos γ⟩ is the average cosine of the single-scattering angle [2]) such as sea water, early arrival photons move along weakly curved trajectories and, as a consequence, exhibit relatively narrow angular distribution and small delay, ∆ = ct − z ≪ z (ct is the photon path length, c is the speed of light in the medium, z is the source-receiver distance). These photons are the least depolarized. Just this component of light field is most important to underwater optical communication and imaging [2,3,6,11].
For the early arrival component of light field, the small-angle approximation can be applied to the non-stationary radiative transfer equation. Within this approximation, many studies of multiple scattering of a pulsed beam were carried out (see, e.g., [2,18,[20][21][22][23][24][25][26][27][28][29][30][31][32][33]). Most of them [20,[23][24][25][26][27]29,[31][32][33] were based on the small-angle version of the Fokker-Planck approximation (or the diffusion approximation in the angular domain) for the collision integral in the radiative transfer equation. The results obtained within this exactly solvable model correlate only qualitatively with experimental data on the structure of a pulsed beam in sea water. This model based on the assumption that the mean square ⟨γ 2 ⟩ of the single-scattering angle exists within the small-angle phase function representation (i.e., that the phase function wings falls off rapidly [28][29][30]) can not describe adequately many of the pulsed beam characteristics (e.g., the distance dependence of the power attenuation, the temporal, angular and spatial broadening etc.). As analysis shows, the reason is that the phase function of sea water decreases rather slowly with increasing the scattering angle γ (roughly as 1 γ 3 [7,9,15,28]). This feature should be allowed for in solving the radiative transfer equation for sea water.
The aforesaid refers equally to the description of the polarization state of multiple-scattered light in sea water as depolarization results from scattering through relatively large angles [34][35][36]. The results derived within the small-angle diffusion approximation (see, e.g., [31,32,[37][38][39]) can not be applied directly to sea water and other highly forward scattering media, such as aqueous suspension of polystyrene particles for which the depolarization of pulses has been studied both experimentally [40][41][42] and numerically using Monte Carlo [17,41] and DISORT [19] codes. Information on the degree of light polarization in the transmitted pulse is important for applying the polarization difference technique (see, e.g., [42,43]). When measuring the intensity difference between cross-polarized components, the pulse tail (i.e., the contribution from delaying photons) is cut off, which can be used in optical communication and imaging [42,43].
In what follows the distribution of light field and its polarization state in a δ-pulse propagating through sea water are studied. The basic mode approximation in the vector radiative transfer equation (VRTE) which we have applied previously to a stationary case (see, e.g., [35,36,44,45]) is generalized to calculating the temporal profiles of the intensity and the degree of polarization in the pulse. Our approach are based on a solution of the eigenvalue problem for the basic modes. The results of both analytical, within the small-angle approximation, and numerical calculations are presented. When performing the analytical calculations, we rely on the small-angle version of the Reynolds-McCormick phase function [46] and the Rayleigh approximation for the scattering matrix. Actual data on the sea water phase function [47,48] and the Voss-Fry scattering matrix [49] are used in our numerical calculations. The power-law parametrization is shown to be valid for the eigenvalues appearing in the expressions for the Laplace transform of the basic modes with respect to time. This enable expressing the distribution of light field in the pulse and the degree of polarization in a self-similar form. Our results are in agreement with data of experiment [7] and Monte Carlo simulations [17].

Vector Radiative Transfer Equation
Consider a pulsed beam of polarized light propagating in a scattering and absorbing medium. The medium is assumed to be a statistically isotropic disordered ensemble of scatterers. The intensity and the polarization state of scattered light are described by the Stokes column vector [50,51] The Stokes parameters I, Q, U and V, and the components E ∥ and E ⊥ of the electric field appearing in Equation (1) are defined in the system of unit vectors e ∥ = ∂Ω ∂θ, e ⊥ = [Ω, e ∥ ], and Ω. The unit vector Ω=(sin θ cos ϕ, sin θ sin ϕ, cos θ) is the direction of propagation of the transverse electromagnetic wave, the vector e ∥ lies in the plane formed by the vectors Ω 0 and Ω (the unit vector Ω 0 is directed along z-axis of the reference frame), the vector e ⊥ is perpendicular to this plane. The brackets ⟨. . .⟩ denote statistical averaging.
The intrinsic properties of the medium is described by the scattering matrixF(cos γ) where cos γ = ΩΩ ′ . For a macroscopically isotropic and symmetric medium, the scattering matrixF(cos γ) has the block-diagonal form (see, e.g., [50,51]): For the forward scattering (Ω = Ω ′ , cos γ = 1), theF matrix is diagonal:F(1) = diag(a 1 , a 2 , a 3 , a 4 ) and a 2 (1) = a 3 (1) [51]. The matrix element a 1 appearing in matrix (2) is the scattering phase function normalized by the relation Under multiple scattering conditions, it is convenient to go from the linear basis (see Equation (1)) to the circular basis [52], where the electric field is defined in the system of unit vectors e ± = (e ∥ ± ie ⊥ ) √ 2 as superposition of waves with right (+) and left (−) circular polarizations. For the circular basis, column vector [51,52] is an analog of the Stokes vector (1). The vectorÎ obeys the non-stationary VRTE of the form [51,52] 1 c where c is the speed of light in the medium, σ tot = σ + σ a is the coefficient of total extinction, σ and σ a are the coefficients of scattering and absorption, respectively. The phase matrixd(Ω, Ω ′ ) entering into Equation (5) is given by [51] d(Ω, where χ ± = π − (β ± β ′ ), the angles β and β ′ (see Figure 1) are defined by formulas Functions cos 2β ′ and sin 2β ′ are obtained from Equation (7) by interchanging µ and µ ′ .

Basic Mode Approximation
In propagation of light through sea water and other media with weakly refractive large (the size d is larger than the wavelength λ) inhomogeneities single scattering through small angles dominates. In this case the off-diagonal elements of the phase matrix appearing in Equations (5) and (6) turn out to be small as compared to the diagonal ones b 1 , b 2 , a ∼ a 4 (see, e.g., [35,36,45]). In the first approximation we can neglect the off-diagonal elements. Then Equations (5) and (6) are reduced to three independent transfer equations [34][35][36]44,45,53]. These equations describe propagation of the individual modes I, W = (Q ± iU) exp (±2iϕ) and V appearing in Equation (4).
The specific intensity I is subject to the ordinary scalar transfer equation (see, e.g., [2,21]) Quantity W was named [44] by the basic mode of linear polarization. The transfer equation for W has the form The separate equation for W was also derived in a different way, without resorting to the circular representation, in [53]. For a circularly polarized light (I = V in the source), the basic mode approximation is similar to that discussed above. When neglecting the off-diagonal elements, we arrive at the individual transfer equation for the basic mode V of circular polarization, Equations (11)-(13) are valid for arbitrary initial polarization of light. In particular, for the elliptically polarized light, the initial values of the functions W and V entering into Equations (12) and (13) are equal to W = I ⋅ 1 − P 2 C , V = I ⋅ P C , where I and P C are the intensity and the degree of circular polarization of the source.
For a δ-pulsed unidirectional source emitting light along z-axis, the initial condition to Equation (11) can be written in the form (see, e.g., [2,11,54]) where unit total energy of radiation is assumed in the pulse. The solution to Equation (11) with the initial condition (14) determines the Green function of the problem. The light field from a source with given time and angular characteristics is expressed in terms of convolution of the Green function with the initial distribution [2,21].
Applying the Laplace transform with respect to path length ct, we can be represent a solution to the scalar transfer equation (11) as follows where µ = ΩΩ 0 , ∆ = ct − z is the "excess" path length (i.e., the difference between the photon path length ct and the distance z), ε For typical scattering and absorption coefficients of sea water, the error in Equation (16) does not exceed a few percents (see Appendix B in [36]). The corresponding solutions to the transfer equations for the basic modes W and V can also be represented in the form similar to Equation (15) . The equations for W and V (see Equations (12) and (13)) differ from the scalar transfer Equation (11) by the form of the effective phase functions, a (2,3) + ⋅ exp(2i(χ + − ψ) and a 4 , respectively. The difference between these phase functions and phase function a 1 entering into Equation (11) results in the fact that the integral of either of the effective phase functions over directions is less than unity (as opposed to Equation (3)). This gives rise to nonzero effective "absorption" in Equations (12) and (13) (even in the absence of true absorption). The effective "absorption" in Equations (12) and (13) is responsible for the additional attenuation of W and V as compared to intensity I and describes the effect of depolarization of linearly and circularly polarized light [31,32,[34][35][36]44,45].
There are two different reasons (see, e.g., [36,45]) for wave depolarization in a scattering medium. These reasons were first pointed out in the context of wave propagation through a turbulent atmosphere [55][56][57].
The "geometrical" depolarization is due to the Rytov rotation [58] of the polarization plane. The depolarization of linearly polarized light results from superposition of randomly oriented polarizations of the waves propagating along different random paths. The situation is different for circularly polarized light. Circularly polarized light can be presented as a superposition of two linearly cross-polarized waves shifted in phase by π 2 . The Rytov rotation has no effect on the phase shift between them. The wave remains circularly polarized. The pure "geometrical" depolarization can be realized in the limit a 1 (cos γ) = a 2 (cos γ) = a 3 (cos γ) = a 4 (cos γ) when the effective phase function appearing in Equation (12) is expressed in terms of a 1 (cos γ) and the angle χ + [35,36,44].
The difference between the diagonal elements of the scattering matrix (2) a i (cos γ), i = 1 ÷ 4 are responsible for the "dynamical" mechanism of depolarization [35,36,44]. The "dynamical" depolarization occurs independently of the initial polarization of light (circularly polarized light depolarizes only due to the "dynamical" reason).

Model of Single Scattering of Light in Sea Water
In sea water, single scattering of light by suspended particles is directed mainly forward, which is reflected in the angular profile of the phase function. Deflection of light by turbulent perturbations occurs through extremely small angles and does not affect the distribution of light scattered by suspended particles. In the radiative transfer theory, the effect of turbulence is not taken into consideration. Analysis of experimental data and numerical calculations (see, e.g., [1,7,9,15,47,48,59]) shows that the phase functions of sea water and other media with large inhomogeneities (e.g., aqueous suspension of polystyrene particles which is used to simulate scattering properties of sea water) fall off with increasing angle γ by a power law a 1 (γ) ∼ (γ 0 γ) α , where the exponent α varies between 2.5 and 3.5, and typical values of the angle γ 0 are in the range 1 ○ -5 ○ [1, 7,9,15,47,48,59]. A number of examples of a 1 (γ) is illustrated in Figure 2. For theoretical modeling the angular profile of single scattering in sea water, the Henyey-Greenstein phase function and its modification proposed be Reynolds and McCormick [46] a 1 (γ) = K α 4π are also widely used [7,9,15,21,50] (e.g., in Monte Carlo simulations [15], the phase function of sea water was fitted by Equation (18) with α = 2.7 and g = 0.96).
The scattering matrix elements of sea water can be inferred from data of measurements [49] and calculations (see, e.g., [61] and references therein). According to [49,62], the values of these elements for sea water are somewhat different from their Rayleigh values (see Equation (10)). The angular dependence of the elements a (2,3) + and a 4 appearing in Equations (12) and (13) is shown in Figure 2. The data for sea water presented in Figure 3 can be approximated by the relations where α + = 0.72, β + = 1.38, α 4 = 3.77, β 4 = 4.14. The phase function a 1 (γ) is conveniently expanded in a series of Legendre polynomials The expansion coefficients a 1 (l) are used in solving the radiative transfer equation (see Appendix A and, e.g., [50]). The scattering matrix element a 4 (γ) appearing in Equation (13) can be expanded similarly. A somewhat different representation is valid for the scattering matrix element a For the highly forward scattering medium such as sea water, the expansion coefficients a 1 (l), a (2,3) + (l) and a 4 (l) differ from each other only at relatively small values of l. This is illustrated in Figure 4 where the results of numerical calculations of these coefficients for sea water and its imitations are presented.

Distribution of Light Field in a Pulse
Multiple scattering of light in sea water occurs predominantly through small angles. This primarily relates to early arrival photons in the pulse. In this case, the radiative transfer Equations (11)-(13) can be from the outset transformed within the small-angle approximation.
To apply the small-angle approximation to Equation (11), we expand the coefficients on the left-hand-side of Equation (11) in a series in the small quantities θ x and θ y (θ x and θ y are the projections of the unit vector Ω on the x and y axises, Ω = θ x , θ y , 1 − (θ 2 x + θ 2 y ) ) and retain only the first nonvanishing terms. As a result, we arrive at the following equation [28]: where θ = θ x , θ y andĨ = I exp(σ a ct). The condition (14) takes the form Using the Laplace transform with respect to the variable ∆ and the Fourier transform with respect to the angular variable θ, the solution of the problem (22) and (23) can be written as (see, e.g., [28,63]): where ε (n) I and Φ (n) I are the eigenvalues and normalized eigenfunctions of the equation with J 0 (x) is the Bessel function. Equations (24) and (25) can be cosidered as small-angle versions of the expression (15) and the characteristic Equation (A4) (see Appendix A). Correspondingly, Equation (26) is a small-angle limit in the well-known expression for the coefficients of the phase function expansion in a series in the Legendre polynomials (see Equation (A7)). The intensity of non-scattered light can be excluded from Equation (24) In this case, the expansion of 1 − a 1 (l) in powers of l begins with the term proportional to l α−2 where Γ(x) is the gamma-function. Under the assumption that the multiple scattering angles θ are greater than the characteristic single-scattering angle γ 0 , only the first term in the expansion (28) should be taken into account (simultaneously, we neglect from here on the contribution from the non-scattered light, exp(−(σ + σ a )z)δ(t − z c)δ(θ)). In this case, the eigenvalues ε where c (n) α is the α-dependent numerical coefficient, and σ tr = σ(1 − ⟨cos γ⟩) is the transport scattering coefficient, for the phase function (18) σ tr = 2σ((α − 2) (4 − α))(γ 0 2) α−2 . Then the intensity can be expressed in terms of dimensionless variable ∆ z(σ tr z) 2 (α−2) . Integrating Equation (22) over angles and taking into account the dependence of the intensity on the variable ∆ z(σ tr z) 2 (α−2) , we can find the mean square of the angle θ at given z and t, For early arrival photons, ∆ < z(σ tr z) 2 (α−2) , only the minimum eigenvalue ε (n=0) I and, correspondingly, the first term in the sum (24) make a contribution to the intensity I(z, θ, t). In this case, the integral over p is governed by rather great values of p, p ≫ σ tr . The inverse Laplace transform in Equation (24) can be performed using the stationary phase method (the saddle point is equal , the numerical values of c (n=0) α are given in Figure 5). As a result, the distribution of light field in the pulse can be presented in the form where ζ 1 (α) = ((α − 2) 2π) 3 The approach based on Equations (22)-(25) can also be extended to a narrow beam geometry (see Appendix B).
The distribution (31) is valid for early arrival photons, including the peak of the pulse. The tail of the pulse (∆ > z(σ tr z) 2 (α−2) ) is described by the Laplace transform of the intensity at relatively small values of p, p < σ tr (σ tr z) α (α−2) . In this case, we can take advantage of the results of calculations for relatively weak absorption, σ 1−2 α a σ 2 α tr z ≪ 1 [64]. Generalizing the results [64] to the pulse propagation problem, we find that the integral-over-angle intensity I(z, t) decreases with increasing the "excess" path ∆ as where ζ 3 (α) = 2 α 2−2 (4 − α) α. Thus, the behavior of the temporal sweep of the pulse at relatively great values of the "excess" path ∆ is directly related to the angular profile of the phase functions at γ > γ 0 . The results presented above (see Equations (29)-(32)) are applicable within the range γ 2 0 ≪ ∆ z ≪ 1. The first inequality is the multiple scattering condition under which only the first term proportional to l α−2 is kept in Equation (28). The second inequality is due to the small-angle approximation.

Degree of Polarization
For the basic modes of linear and circular polarizations, solutions of Equations (12) and (13) can also be presented in the form that is similar to Equation (24). The corresponding solutions differ from Equation (24) only by the specific eigenvalues, ε where ∆a W,V (l) is the difference between a 1 (l) and the Bessel transforms of the "effective" phase functions appearing in Equations (12) and (13). The small-angle transfer equation for the basic mode of linear polarization differs from Equation (25) only by the "perturbation" term σδa W (l) [35,36,44], The first term in Equation (34) is responsible for the "geometrical" depolarization due to the Rytov effect. The second term describes the "dynamical" depolarization. Substituting Equation (34) into Equation (33), we obtain where η W (α) = Γ(4 − α 2) 2 4−α 2 . When deriving Equation (35), we use the Gaussian ansatz for the eigenfunction Φ (p) p. For the basic mode of circular polarization, the difference δa V is only due to the "dynamical" mechanism of depolarization. The value of δa V turns out to be two times greater than δa Wdyn appearing in (34) [35,36,44]. The quantity δε V is determinate by where η V (α) = (4 − α)Γ(3 − α 2) 2 4−α 2 . The results (35) and (36) rely on the small-angle phase function (27) and the Rayleigh approximation for the scattering matrix elements (see Equation (10)). Evaluating the integral in Equation (24) and the integrals in the corresponding representations for the modes W and V by the saddle-point method, we obtain the expressions for the degree of linear and circular polarization of the pulse where ⟨θ 2 ⟩ z,t is determined by Equation (30). From Equations (37) and (38) it follows that the degree of polarization depends only on one dimensionless variable and decreases with increasing the "excess" path ∆.

Results of Numerical Calculations
The results obtained above in Sections 3 and 4 can be validated using the numerical solution of the eigenvalue problem (see Appendix A) without resort to the small-angle approximation. In what follows it is assumed that all eigenvalues and eigenfunctions for the basic modes I, W and V correspond to the number n = 0.
First we outline the numerical results for ε I (p). The p-dependence of ε I (p) at p > σ tr can be approximated by a power law, The values of c ν and ν for sea water (Kopelevich's model [47,48]) and its imitation are given in Table 1. The calculations of ε I (p) were carried out with Equation (A4) (see Appendix A). Table 1. The values of cν and ν for sea water (Kopelevich's model [47,48]) and its imitation. The eigenvalues ε I , ε W and ε V , and the differences δε W = ε W − ε I , δε V = ε V − ε I as functions of p σ tr are shown in Figure 6. For sea water, the results of measurements by Voss and Fry [49] were used for calculating the scattering matrix elements a (2,3) + and a 4 . For the Henyey-Greenstein phase function, the calculations were carried out within the Rayleigh approximation (see Equation (10)). The calculations of a The angular dependence of the eigenfunctions Φ I , Φ W and Φ V is illustrated in Figure 7 for the reference case of the Henyey-Greenstein phase function. The difference between the eigenfunctions is small in the forward hemisphere, suggesting that the angular profile is unaffected by the initial polarization of light or, in other words, the decay of the initial polarization is virtually angle independent. The contribution from the geometrical depolarization to δε W is shown in Figure 8. This contribution is the main for aqueous suspension of relatively large polystyrene spheres (d = 0.993 µm, λ = 633 nm), and for the reference case involving the combination of the Henyey-Greenstein phase function with the Rayleigh scattering matrix and is comparatively modest for sea water. The features of the Voss-Fry scattering matrix [49] that are inherent in scattering by non-spherical particles make the dynamical depolarization principal for sea water. As p increases, δε W and δε V decreases according to a power law, The numerical values of D ν 1,2 and ν 1,2 are presented in Table 2. For the latter two lines of Table 2 the values of ν 1 and ν 2 coincide to the corresponding values of 3 − α 2, correlating with Equations (35) and (36).

Discussion
The results obtained above give an insight into the distribution of light field in the pulse and enable calculating semi-analytically its main characteristics, including the polarization state, in sea water.
For the intensity of radiation in the pulse, the small-angle approximation can be applied in combination with an adequate parametrization of the sea water phase function. A peculiarity of light scattering by sea water is that the wings of the phase function at γ > γ 0 (γ 0 is of the order of a few degrees) fall off rather slowly, according to a power law. Therefore, the parametrization such as the Henyey-Greenstein function or the Reynolds-McCormick one grasps the key feature of light scattering by sea water in the forward hemisphere. This is confirmed by the data presented in Figures 2 and 4. Additionally, to validate the power-law approximation (28) in application to a realistic situation, we compare the results of our calculations with experiment [7]. Substituting Equation (28) with the same parameters as in Monte Carlo simulation [7,9] (i.e., α = 3, γ 0 = 0.07) to the small-angle solution of the transfer equation (see, e.g., [21]), we obtain the attenuation curve which agrees with the experimental data (see Figure 9). The relative contribution of multiple-scattered radiation to the received signal increases with the optical distance, resulting in the decrease of the attenuation rate as compared to Beer's law. Optical distance ( + a )z Normalized received power Figure 9. Normalized received power versus optical source-receiver distance (symbols [7]). The results of calculations within the small-angle approximation for the Henyey-Greenstein phase function (γ 0 = 1 − g = 0.07) are shown by solid curve. The attenuation of non-scattered radiation, Beer's law, is displayed by straight line.
Using the power-law approximation for the phase function, we can represent the distribution of light field in the pulse in the self-similar form as a function of dimensionless variables (see Equations (31) and (32)). These variables are certain combinations of the "excess" path ∆, the transport scattering coefficient σ tr , the source-receiver distance z and the exponent α appearing in the parametrization of the phase function. Correspondingly, the main characteristics of the light field distribution (the effective time spread, dispersion in angles and etc.) are expressed in terms of these variables.
The important result for the variance of the angular distribution of light field, ⟨θ⟩ z,t (see Equation (30)), enables controlling the validity of the small-angle approximation at given z and t.
From the results obtained it follows that the peak in the temporal sweep of the pulse is observed at small values of the "excess" path, ∆ ≪ z, or, for a narrow beam, ∆ − ρ 2 (2z) ≪ z (ρ is the transverse displacement from the beam axis, see Appendix B). This is illustrated in Figure 10 where the distribution of light field in the forward direction and the integral-over-angle distribution are shown. The latter correlates well with the Monte Carlo simulation data [11] (see Figures 3 and 4 in [11]). The calculations of light field with the use of the numerical results outlined in Section 5 are performed in the following way. First, the integral over p in Equation (15) is evaluated by the stationary phase (or saddle-point) method, and, as before, only the contribution from the minimum eigenvalue ε I = ε (n=0) I is taken into account. The saddle-point p 0 is sought from the equation The value of p 0 can be found using the power-law approximation for ε I (p) (see Equation (40)) or, directly, from the numerical values ε I (p) (see Figure 6). In the latter case, the results obtained are valid even beyond the small-angle approximation. Further the integral in Equation (15) is calculated routinely (all p-dependent factor appearing in Equation (15) are taken at the saddle-point p 0 ). Our approach to calculating the light field distribution in the pulse differs from the analytical solutions proposed for the same problem previously (see, e.g., [18,20,[22][23][24][25][26][27]29,33]) in that the specific angular profile of the sea water phase function is allowed for in our calculations.
The basic modes of linear and circular polarization, and the corresponding degree of polarization are calculated similarly to the light field distribution. The degree of polarization can be written as where δε W,V = ε W,V − ε I and p 0 is the root of Equation (42). The numerical values of δε W and δε V as functions of p are given in Figure 6. Their analytical parametrization is described by Equations (35) and (36) or (41). The applicability of the expression for the degree of polarization in the form (43) is not limited by the small-angle approximation.
The results of comparison of our calculations with data of Monte Carlo simulations for aqueous suspension of polystyrene spheres [17] are shown in Figure 11. The numerical simulation [17] was carried out for the linearly polarized pulse propagating through a plane slab of different transport optical thickness σ tr z. Depolarization ratio (1 − P L ) (1 + P L ) was presented in [17] as a function of the normalized "excess path" ∆ z. When going from ∆ z to novel variable σ tr z(∆ z) 3−α 2 (see Equation (39)) we obtain virtually universal pattern of the simulation data (see Figure 11b,d). The agreement between our calculations and the data [17] shows possibility of application of the power-law parametrization (41) to realistic cases.  As noted in Section 2.2 (see also [36,45]), two mechanisms underlie the depolarization of light in scattering media. For sea water, as follows from the numerical results presented in Figures 6  and 8, the dynamical mechanism of depolarization proves to be dominant. This is due to specific feature of the Voss-Fry scattering matrix at relatively small scattering angles γ. From the Voss-Fry measurements (see Figure 3 and also the approximation (19) and the parametrization proposed in [62]) it follows that, for small γ, both (a (2,3) + − a 1 ) a 1 and (a 4 − a 1 ) a 1 are proportional to γ 2 . In the Rayleigh approximation (see Equation (10)) and in the case of aqueous suspension of polystyrene spheres, contrastingly, these quantities are proportional to γ 4 . This distinction can be explained by the contribution of large non-spherical particles to the scattering matrix of sea water. For this reason the effect of circular polarization memory that is typical for aqueous suspension of polystyrene spheres (see, e.g., [36,45,[65][66][67]) can not be observed in propagation of light through sea water. The linearly and circularly polarized beams are depolarized equally (see Figure 6).
It should be noted that the results presented above can also be applied to a steady-state beam of light. To do this, the expressions for the basic modes should be integrated over the "excess" path ∆. This is reduced to putting p = σ a in the corresponding formulas (e.g., in Equations (15), (16) or (43)).

Conclusions
We have studied the distribution of light field and the degree of light polarization in transmission of a pulsed beam through sea water. A novel method has been proposed that enables calculating the distribution and the degree of polarization of light field in a pulsed signal, based on a numerical solution of the eigenvalue problem for the basic modes of the VRTE. The small-angle approximation underlies our analytical results that rely also on the power-law parametrization of the phase function and on the Rayleigh scattering matrix. The numerical calculations beyond this approximation have been performed for the actual data on the sea water phase function and the scattering matrix. The eigenvalues entering into the expressions for the Laplace transform of the basic modes with respect to the time variable have been shown to be written in a power-law form. This enables expressing the light field distribution and the degree of polarization in the pulse in terms of dimensionless variables in a self-similar form. The variance of the multiple-scattering angle has been shown to be proportional to the ratio of the "excess" path ∆ = ct − z to the source-receiver distance z, justifying the application of the small-angle approximation to calculating the early arrival component of light field (∆ ≪ z) in the pulse. We have also found that the depolarization properties of sea water differ from those (see, e.g., [36,45,[65][66][67]) of aqueous suspension of polystyrene spheres. Contrary to [36,45,[65][66][67], as follows from our calculations based on the actual data on the sea water scattering matrix [49], the depolarization rates for linearly and circularly polarized pulses prove to be virtually equal to each other. This can be explained by scattering from non-spherical particles suspended in sea water.
The analysis and results presented above can find application to the problems involving propagation of non-stationary laser beams through sea water, such as underwater communication, imaging and remote sensing.
where P l (µ) and P l 22 (µ) are the Legendre polynomials and the generalized spherical functions, respectively. Detailed properties of P l 22 (µ) can be found in [51]. Substituting expansions of the type of Equation (15) into Equations (11)- (13) and taking into account the expansions of the scattering matrix elements in the generalized spherical functions (see Equations (20) and (21)) we arrive at the following equations for the coefficients entering into Equations (A1)-(A3) and the corresponding eigenvalues: Solutions to Equations (A4)-(A6) are found numerically with the use of truncating at some l max . In such a way, we find a set of the eigenvalues and the eigenfunctions for each mode, I, W and V.
where [. . .] coincides with the right-hand-side of Equation (24). Taking into account only the first term in the expansion (24) and using the Gaussian ansatz for corresponding eigenfunction, we obtain