Probability of High Intensities of the Light Wave Propagating in a Turbulent Atmosphere

We examine statistics of fluctuations of the laser beam intensity at its propagating in turbulent atmosphere. We are interested in relatively large propagating distances and the remote tail of the probability density function. The tail is determined by the stretched exponent, we find its index.

The subject of our work is the theoretical study of the physical properties of a laser beam propagating in a turbulent atmosphere or, more generally, an electromagnetic wave in a turbulent medium.The main effect that is being investigated in the framework of this problem is the diffraction of a laser beam by fluctuations of the refractive index, which are induced by pressure fluctuations of turbulent pulsations.The fluctuations of the refractive index are a random field whose properties are described statistically.Therefore, theoretical predictions of the behavior of the laser beam concern the average values (or the corresponding probability density functions), which are obtained by averaging over the statistics of fluctuations.
The problem of laser beam propagation in the atmosphere has a long history.The basic theoretical results concerning the propagation of a laser beam in a turbulent atmosphere were obtained back in the sixties and seventies of the last century, they are summarized in a number of monographs [1][2][3][4][5].Recently, there has been a revival of interest in this topic [6][7][8][9][10], connected mainly with numerical modeling of the propagation of a laser beam (electromagnetic wave) in a turbulent medium, which allows obtaining detailed information about the process.
The known theoretical results describe mainly the typical behavior of a laser beam in the atmosphere.At the same time, in the equation for the electromagnetic field, the random refractive index is a multiplicative noise, which leads to very non-trivial statistical properties of the electromagnetic field.It can be expected that the probability of rare events when one or another value (say, the intensity of the electromagnetic field) has an abnormally large value will be significantly higher than naive estimates based on the analysis of typical processes.As an example of such behavior, we can cite the statistics of quantum particles in a random potential (see, for example, [11]), which is the multiplicative noise in the Schrödinger equation for a quantum particle.
We theoretically study the propagation of a monochromatic electromagnetic wave in an unlimited turbulent medium.All the characteristic scales of the problem (the size of the wave packet, the length of its propagation) are assumed to significantly exceed the wavelength, so that the description of the wave in terms of the complex envelope Ψ is applicable.Due to the large value of the speed of light, we can assume that the state of the medium does not change during the propagation of the wave and use a stationary approximation to describe its envelope, that is, consider Ψ as a function of the coordinate z in the direction of wave propagation, and a two-dimensional radius vector r in the plane transverse to the direction of propagation.Of course, the state of the turbulent medium varies with time.In the accepted approximation, the envelope Ψ(r, z) adiabatically adjusts to the current state of the medium.
In this paper, we consider the intensity of the electromagnetic wave to be quite small, so that nonlinear effects (such as self-focusing) are insignificant.Then the equation for the envelope Ψ(r, z) in suitable units of measurement has the form of a two-dimensional Schrödinger equation where ∇ is the two-dimensional gradient in the plane transverse to the direction of propagation, and ξ(r, z) is the fluctuating component of the refractive index.By virtue of the equation (1), the envelope Ψ is a field that depends on the implementation of ξ.The latter varies over time, leading to variations in the envelope.Therefore, Ψ can be considered as a random field whose statistical properties can be extracted by averaging over long times.
The distance traveled by the wave is assumed to be greater than the integral turbulence scale, whereas the transverse size of the wave packet is assumed to be smaller than the integral scale.In this case, the refractive index ξ fluctuates rapidly along the direction of wave propagation.We are interested in the integral characteristics associated with ξ, therefore, by virtue of the central limit theorem, the random field ξ can be assumed to have Gaussian probability distribution.It is determined by the pair correlation function of fluctuations in the refractive index: where angle brackets mean averaging over implementa-tions of the state of the medium, r 12 = |r 1 − r 2 | and c is a number.The expression (2) is valid for distances r 12 lying in the inertial interval of turbulence.The constant in the ratio (2) is determined by turbulent fluctuations on the integral scale, and the power correction is determined by fluctuations on the scales ∼ r 12 .Fluctuations of the refractive index ξ in a turbulent medium are determined mainly by pressure fluctuations.For the Kolmogorov spectrum, the exponent of c in the expression ( 2) is c = 5/3.Further, the exponent c is treated as an arbitrary number lying in the interval 1 < c < 2. We focus mainly on the statistical properties of the intensity I of the laser beam (of the electromagnetic wave) at some observation point, which we choose as the origin: The distance z traveled by the beam is assumed to be large, z ≫ 1 in our units.This means that the effects associated with random diffraction are strong.As the initial state, we choose a plane wave, Ψ(r, 0) = 1.A generalization of our calculation scheme to other cases, say, to the initial Gaussian beam, will be published elsewhere.We will be interested in the moments of I n for large values of n.These values determine the probability of events with an abnormally high intensity value.
The average value I n can be expressed as I n = F 2n (0, . . ., 0, z), where F 2n is the 2n-point correlation function of the envelop: This correlation function is represented as a convolution with the Green function G 2n : Here Ψ in (r) = Ψ(r, z = 0) is the initial value of the envelop, and the Green function depends on z and 4n radius-vectors.
The Green function can be represented as a path integral, which is derived from the Schrödinger equation (1).After averaging over the fluctuations of the refractive index in accordance with the expression (2), we arrive at an integral over the variables y j , which are functions of the coordinate ζ, 0 < ζ < z: One takes arguments x j and r j of the Green function as boundary conditions for the trajectories y j (ζ) at ζ = 0 and ζ = z, see Eq. ( 5).The "potential" W (8) is the function of the variables y ij = |y i − y j |.The constant appearing in the expression (2) falls out of consideration, as it should be.It is possible to find an explicit expression for the pair Green function G(r 1 , r 2 , x 1 , x 2 , z): which determines the behavior of the pair correlation function.Here the designations R = (r For the initial state in the form of a plane wave, when Ψ in = 1, the pair correlation function has the simple form [12,13]: Thus, the characteristic distance between the points 1 and 2 at a distance z is r ph = z −1/c ≪ 1.This value makes sense of the phase corruption length due to fluctuations of the refractive index.The characteristic value of the variable X is estimated as z 1/c+1 . The analysis of the 2n-point Green function (5) shows that for z ≫ 1 it has sharp maxima on configurations where there are n pairs of close points x j located at distances of the order of r ph separated by much larger distances [13,14].In this case, the first point of the pair is taken from the set x 1 , . . ., x n , and the second is taken from the set x n+1 , . . ., x 2n .The trajectories y j (ζ) in the path integral (6) starting from close points x j also remain close.In this case, the "potential" W (8) can be approximated by the sum of y c ij for the pairs, and the remaining summands in the sum (8) mutually reduce each other.Then the Green function is factorized.For example, for close couples (x 1 , x n+1 ), . . .(x n , x 2n ) Obviously, there is n! of similar contributions, by the number of splits of x j into pairs.Thus, the space of integration by the initial coordinates in the expression for the moments of intensity (12) for the case of an initial plane wave is divided into n! areas corresponding to different point splits into pairs.The contributions of all such regions are the same and in the main approximation for intensity moments we get I n = n!, which means the exponential probability density function P (I) = exp(−I).Note that P (I) = exp(−I) corresponds to the Gaussian statistics of the Ψ field, natural for a complex field with a random phase.
It is demonstrated in the works [13][14][15] that corrections to the approximation (11) related to the discarded terms in the expression (8) for W are proportional to the parameter α: it is small at z ≫ 1.For the Kolmogorov spectrum of turbulence, the exponent a is equal to a = 11/15.Corrections to the value I n = n! grow as n increases and become essential when αn ∼ 1.In this paper, we determine the value of I n at αn ≫ 1 and demonstrate that they significantly exceed n!.In other words, at I ≫ α −1 , the tail of the probability density function P (I) arises, significantly exceeding exp(−I).
If αn ≫ 1 then the integration space in the integral ( 12) is still divided into n! regions, corresponding to well separated pairs of the trajectories y j (ζ), the regions give identical contributions to I n .Let us take for definiteness the following subsetting into the pairs: (y 1 , y n+1 ), . . .(y n , y 2n ).Introducing the variables Y j = (y j + y j+n )/2 and ρ j = y j − y j+n , where j runs through values from 1 to n, we get a functional representation for the Green function with arguments in the selected area in the form: In the path integral ( 14), some non-zero initial values are implied Y j (0) = X j and zero final values (at ζ = z) are implied for the trajectories Y j and ρ j , the values are dictated by examining I n .As it follows from the representation (6), the action S in Eq. ( 14) is equal to In this expression, the indices j, k run through the values from 1 to n.The quantity U at ρ j , small in comparison with the characteristic values of Y j − Y k , is determined by the formula following from the expression (8): Here the Greek indices α, β designate the components of the vectors in a plane transverse to the direction of wave propagation.
As we will see below, the values of Y j in the integral ( 14) are parametrically large by n.This allows us to calculate the integral over Y j in the expression (14) in the saddle-point approximation.The integral over ρ j cannot be calculated in the approximation.
We first investigate the saddle-point value of Y j , which is determined by the extremum condition δG/δY j = 0. Substituting here the expression (14), we find the equation i 2 Here ρj is the average value of ρ j , it is non-zero because of non-zero initial values X j of the trajectories Y j .When deriving the equation ( 18), we replaced the average of the product of ρ jα ρ kβ with the product of averages.The reason for this is that the fluctuations of ρ j are determined mainly by the second term in the action (15), which is diagonal by j.Therefore, the contribution of fluctuations to the mean ρ jα ρ kβ , where j = k, turns out to be negligible.
In field theory, the average value of the fluctuating field is found from the condition of the extremum of the socalled quantum effective action (see, for example, [16]).
In the problem under study, this condition is equivalent to the relation δG/δ ρj = 0, which determines the average ρj .Utilizing the expression ( 14), we find the equation i 2 At deriving the equation we have neglected the contribution of the second term in the action (15).The reason is that the fluctuations of ρ j (which are determined mainly by this term) turn out to be much larger than ρj , which is why the dependence of the second term on ρj is weak.Further, we will justify this neglect.Note that the extremum condition for the initial value of ρj (0) leads to the condition dY j /dζ(0) = 0, due to the structure of the action (15).The equations (18,19) show that the saddle-point value of the expression ( 14) is written as The last equality in (20) is obtained after integration in parts, taking into account the equation ( 19) and the boundary conditions.The value (20) is negative because Y j are real quantities, whereas the averages ρj are purely imaginary quantities, and the function positive for the values of the exponent c considered here: The system of equations (18,19) enables one to evaluate the saddle-point quantity (20).Assuming that all values of Y j together with their boundary values have the same order of Y j ∼ X, as well as the average ρj ∼ i ρ, we arrive at the estimates Note the absence of dependence on n in the estimation for ln G sp .Now we turn to accounting for the fluctuation contribution to G 2n .Fluctuations of ρ j can be estimated as z −1/c , the value is much larger than the average ρ (21).Therefore after shifting ρ j by its mean value and shifting Y j by its saddle-point value, the effects associated with the "interaction potential" ( 16) can be neglected (estimation of the accuracy see below).As a result, we arrive at a factorized approximation of the type (11).Thus, the expression for the Green function is factorized G 2n = G sp G f l , where G sp is determined by the expression (20) and the fluctuation factor is given by the product of the pair Green functions The expression standing here for the pair Green functions is derived from the general formula (9) after substituting r = R = 0 and taking into account that the initial values of the shifted variables Y j − X j are zero.
Returning to the expression (12) for the moments of intensity, we find The factor G sp depends solely on the coordinates of the centers of the pairs X j = (x j + x j+n )/2, therefore, after integration over differences x j − x j+n we find For a large n, the value of the integral (24) is estimated as: where C is a constant of the order of one.Here X is the characteristic value of variables X j introduced above.
Optimizing the expression (25) by X, we find: The substitution of Eq. (26) into Eq.( 25) gives the desired asymptotics of the high intensity moments: The corresponding asymptotics of the probability density function P (I) has the form For the Kolmogorov spectrum β = 7/13.Since β < 1, we conclude that there is a higher probability of large values of I in comparison with the exponential probability density exp(−I).Formula (28) is in accordance with the results of recent numerical simulations [17].Note that the expression (28) means a significant non-Gaussianity of the statistics of the wave field Ψ, formed by the chaotic scattering on the turbulent fluctuations.The asymptotic law (27) turns to I n = n! at αn ∼ 1, or at z ∼ n c/(4−c 2 ) , as it should be.A similar statement is true for the probability density function, which is matched with the exponent exp(−I) at I ∼ 1/α.Eq. ( 26) enables one to find the characteristic size of the region at the front of the initial wave defining I n at n ≫ 1/α: At αn ≫ 1 the quantity (29) is much larger than the estimate z 1/c+1 for X obtained in the Gaussian regime.
Physically, this means that to create high intensity values, it is required to collect energy from a large area of the original wave.The expression (29) implies an estimate for the average values ρ (21 These averages are small in comparison with the characteristic amplitude of fluctuations: ρ ≪ r ph , which justifies the approximation (22).Note, however, that the saddle-point action (20) is determined precisely by these averages, and in action (20) this smallness is compensated by a large number of terms.The found asymptotic law (27) is bounded from above by n.This limitation is determined by the mean square of the fluctuation component of the interaction of pairs between themselves that we have discarded.The assessment is fair for him .
The smallness of this value in comparison with the main value of the fluctuation action, estimated as ∼ n, leads to the condition of applicability of our approach n ≪ α −4/c .The inverse limit case requires a special analysis.
We have theoretically established that when a laser beam propagates in a turbulent medium at distances where diffraction on random fluctuations of the refractive index plays a significant role, the probability of abnormally large intensity fluctuations is significantly higher than estimates made on the basis of Gaussian statistics of the electromagnetic wave envelope Ψ.This is due to the multiplicativity of the refractive index in the equation for the envelop and it is a universal property of such stochastic systems.We have found the form of the probability distribution function of abnormally large intensity values.Our conclusions are made for the simplest case of the initial plane wave, although the analysis scheme itself is applicable to other initial forms of the laser beam.The results of their research will be published elsewhere.The case of extremely high intensity values requires a sepa-rate analysis.A separate analysis also requires to take into account the nonlinearity (for example, the effect of self-focusing), which may be significant for the large intensities of the electromagnetic wave we are considering.
The work is supported by the scientific program of the National Center for Physics and Mathematics, the project "Physics of high energy densities", stage 2023-2025.The authors thank P.M. Lushnikov for useful discussions.