Temporal analysis of laser beam propagation in the atmosphere using computer-generated long phase screens

Temporal analysis of the irradiance at the detector plane is intended as the first step in the study of the mean fade time in a free optical communication system. In the present work this analysis has been performed for a Gaussian laser beam propagating in the atmospheric turbulence by means of computer simulation. To this end, we have adapted a previously known numerical method to the generation of long phase screens. The screens are displaced in a transverse direction as the wave is propagated, in order to simulate the wind effect. The amplitude of the temporal covariance and its power spectrum have been obtained at the optical axis, at the beam centroid and at a certain distance from these two points. Results have been worked out for weak, moderate and strong turbulence regimes and when possible they have been compared with theoretical models. These results show a significant contribution of beam wander to the temporal behaviour of the irradiance, even in the case of weak turbulence. We have also found that the spectral bandwidth of the covariance is hardly dependent on the Rytov variance. ©2008 Optical Society of America OCIS codes: (010.1300) Atmospheric and oceanic optics: Atmospheric propagation; (010.1330) Atmospheric and oceanic optics: Atmospheric turbulence


Introduction
Temporal turbulence-induced phenomena in optical propagation through the atmosphere were first studied by Tatarskii for plane waves [1].An excellent summary of his work can be found in a review paper by Lawrence and Strohbehn [2].Later, Clifford extended Tatarskii's theory to spherical waves [3], and Ishimaru applied it to obtain the frequency spectrum of the temporal amplitude fluctuations on the beam axis of a Gaussian wave [4].Shelton gave an accurate model for the calculation of the power spectrum in realistic conditions, obtaining a general good agreement with the experiments made on a ground-to-satellite optical link [5].Finally, Andrews and Phillips obtained general and more simplified expressions for the temporal covariance of the irradiance fluctuations and its corresponding power spectrum, both on the optical axis and in points located at a certain distance from it [6].
The analysis of the temporal behaviour of the irradiance is the primary step to provide significant merit figures in an optical link.Actually, without the temporal information, the overall performance of a laser link remains essentially incomplete.In a communication system that suffers any kind of fading due to multipath interference, rain, fog or turbulence, the mean signal-to-noise ratio or the mean bit error rate are not the relevant parameters, so we have to find the instantaneous values of these figures.Otherwise, the burst error rate and the overall availability of the system cannot be known [7].Moreover, the temporal characterization of the system is necessary to develop optimal schemes for detection and coding [8][9][10].
However, all these models suffer from the fact that they rely on the conventional Rytov analysis, which is limited to weak turbulence conditions.Alternatively to theoretical models, numerical experiments have proved to be an important tool in free-space optical propagation research.The use of simulation techniques has been useful to put in evidence the importance of the beam wander effect into the performance of ground-to-satellite uplinks [11][12][13][14].Nevertheless, the practical difficulties encountered in the generation of large phase-screens have prevented the extension of these techniques to the study of temporal phenomena.
Numerical experiments on wave propagation through the atmosphere customarily rely on the parabolic or Fresnel approximation of the wave equation.The turbulent medium is represented by series of independent phase screens, and the split-step technique is used to propagate the beam.The screens are generated mainly by using two different techniques.The first of them is the FFT-based or spectral method, first proposed by Macaskill and Ewert [15] and later improved by Lane et al. [16] and Frehlich [17] with the addition of subharmonics.The second one is the covariance matrix method, developed by Harding et al. [18].This second one is more accurate.Nevertheless, the huge amount of computational effort required to evaluate large covariance matrices limits the practical size of the generated screens, thus making necessary the use of fractal interpolation techniques [16,18].Though some attempt has been reported in the field of adaptive optics to extend these methods to the generation of the long phase screens required in the simulation of temporal phenomena [19], no results have been provided concerning Gaussian beam propagation in optical links.
In the present work we extend the applicability of the numerical split-step technique to the temporal analysis of Gaussian beams propagating through the atmosphere.With this goal in mind, we have implemented an algorithm to generate and handle long rectangular phase screens.The screens generated with this method are used to simulate the temporal behaviour of the irradiance in terrestrial horizontal links, in which we consider the refractive index structure function C n 2 as a constant.The paper is organised as follows: in section II we provide the theoretical background for the temporal analysis of the optical irradiance of a wave propagating into the atmospheric turbulence, paying special attention to some approximations commonly accepted in the theory and their limitations.In section III the computational method developed for temporal analysis by using large moving phase screens is explained.Finally, the results obtained are discussed in section IV.

Classical procedure
Analytical methods in atmospheric optical propagation routinely apply the frozen atmosphere hypothesis, or Taylor's hypothesis, to extend the spatial statistical properties of a travelling wave to the study of the temporal behaviour of the irradiance at the detector plane.The conversion from the spatial covariance to the temporal one is performed by means of a variable change.The spatial coordinate ρ representing the position in a plane perpendicular to the wave propagation direction is transformed into the time variable τ by expressing it as a function of the transverse component of the wind velocity V ⊥ .The specific conversion formula depends on the type of wave that we are considering.Following Andrews and Phillips [6], we have τ ρ ⊥ = V (1) for a plane wave, and for a spherical wave, where ξ = 1 -z/L, with z being the wave propagation direction from the emitter to receiver and L being the path length.In the case of a Gaussian beam the conversion formula is being a complex parameter of the Gaussian beam, with k being the wave number in the free-space, and W 0 and R 0 being, respectively, the beam width and radius of curvature at z = 0.
The temporal power spectrum S I is straightforwardly derived from the temporal covariance B I by means of the Fourier transform The most relevant parameter in order to assess the availability of the link is the mean fade time, which is defined as the average continuous period of time in which the received power remains below a specified level.This time period is expressed mathematically as follows where F is the fading level below the average power and F 0 is some specific fading value.〈n(F 0 )〉 is the expected number of fades per second below the F 0 level and P{F < F 0 } is the probability of fading, i.e. the probability of the received power being under the F 0 level.The probability of fade can be obtained directly from the probability density function (PDF) of the received irradiance, without the need to resort to the temporal analysis of the wave propagation.Currently, different PDF models proposed either for the irradiance or for the power collected by a finite aperture may be found in the literature [6].The calculation of the expected number of fades per second 〈n(F 0 )〉, however, requires the knowledge of the temporal statistics of the signal, although even with this information this quantity is not easy to obtain analytically.When the analysis is performed by using the irradiance as the working variable, it is necessary to assume that the PDF of its temporal derivative is a Gaussian function to obtain a closed analytical expression for 〈n(F 0 )〉.If the log-amplitude of the wave is used instead, the log-amplitude PDF itself has to be assumed to be Gaussian.Otherwise only approximated formulae can be obtained [6,[20][21].The above described procedure followed to obtain the mean fade time in a free-space laser link works well under the weak turbulence regime.However, there are several steps in which significant divergences between theory and experiment may be introduced in horizontal links, when the turbulence level is moderate or high.In the following subsections we will address to two of these issues

Loss of the Gaussian shape along propagation
The first key assumption made in the theory is given by Eq. (3).At first glance, the conversion from spatial to temporal phenomena obtained from this equation is adequate only in those conditions in which the perturbation of the wave due to the turbulence is small, i.e. the kind of situations that are well described by the Rytov approximation.In the case of laser beams, this means that the beam is expected to remain approximately Gaussian during the whole path, so that the effect of the turbulence can be considered essentially as a perturbation affecting the main wave.Unfortunately, this is not the case in many situations of interest, and particularly in low altitude horizontal links, where typical values of the structure constant of the air refractive index C n 2 varies, roughly speaking, between 5×10 -13 and 5×10 -14 at sea level, and the laser beam may lose its Gaussian shape after propagating a few hundred meters under these conditions, depending on the wavelength and the initial values of the beam parameters W 0 and R 0 .
This phenomenon occurs as the accumulated phase perturbation due to the turbulence distorts and eventually breaks the Gaussian wave, giving rise to several bright centers in the instantaneous field profile.The breakdown of the beam can be observed both experimentally and in a numerical simulation, by representing the field of the wave propagated over a long enough distance.Nevertheless, the important point to make here is that even the long-term average irradiance may lose its expected Gaussian shape.
In order to measure the degree of departure from the Gaussian shape, we have performed a series of propagation simulations over a fixed distance and taking different C n 2 values.The beam parameters used in the simulations are λ = 1.064 μm for the wavelength and W 0 = 2 cm for the beam waist radius.In order to apply the split-step method, the ray path has been divided into 30 propagation steps, each of them covering a 200 m propagation distance.At each step, the turbulence-induced perturbation in the wave phase is simulated by a phase screen containing 512x512 points on a square grid of 1.2 mm spacing, which has been generated assuming a pure Kolmogorov turbulence.For a pure Gaussian wave, the beam width W is usually defined as the distance from the beam center at which the field amplitude decays by 1/e from the peak value.In a practical case, however, it is more reliable to obtain the beam width by measuring the power encompassed by a circle centered at the peak of irradiance.The amount of power that is located within the circle area is directly related to the circle radius.On a pure Gaussian wave, the 1/e radius circle contains the 86.5% of the total power carried by the beam.However, other values may be used instead, as shown in Table 1.Provided that the beam is perfectly Gaussian, all these criteria lead to the same beam width W. On the contrary, a non-Gaussian beam can be easily detected by the lack of convergence in the results obtained.98.17 % In our case, we have measured the W value of the long-term spot size at the receiver by using the three following expressions ), the beam remains nearly Gaussian during the whole path, as shown by the fact that all three criteria given by Eqs. ( 6)-( 8) lead to the same result.In the other two cases (C n 2 = 5•10 -15 and C n 2 = 10 - 14 m -2/3 ), however, each criterion provides a different value for the estimated beam width.It becomes apparent from these curves that the averaged field deviates more and more from the Gaussian profile as the turbulence accumulated along the wave path increases.This departure is shown in Fig. 2, where we have plotted the average field profiles of the beam, measured at a distance L = 5 km in the same three cases.
As a conclusion, the transformation given in Eq. ( 3) should be considered only as an approximation.This is true not only for the analytical models, but also for the conventional numerical simulations, in which the calculations are made without paying attention to the temporal correlation between the phase screens in successive repetitions.Hence the problem of how to transform the spatial analysis in a temporal one still holds.).The deviation of the field profile from the Gaussian shape increases with the accumulated turbulence along the beam path.Dashed black lines have been obtained by using Eq. ( 6), dot-dashed blue lines have been obtained by using Eq. ( 7), and red solid lines have been obtained by using Eq. ( 8).The grid parameters correspond to an effective inner scale value l 0 = 2.4 mm.Fig. 2. Simulated average field profile (blue lines) obtained at a distance of 5000 m from the transmitter, for the same values of the turbulence strength as in the previous figure.Red and magenta lines represent the theoretical Gaussian shape for each case.The numerical results again match well with the theoretical Gaussian profile for a structure constant C n 2 =10 -16 m -2/3 , whereas for higher values the difference becomes more and more apparent.

Statistics of non-homogeneous fields
Another critical point in the analysis of the temporal behaviour of the irradiance at the detector plane is the fact that the stochastic process under consideration is not homogeneous in space, as usually assumed.On the contrary, it has been shown that the scintillation index for a Gaussian beam is a function of the distance from the beam axis [22].This phenomenon can be easily understood by considering that the diffracted light will have greater influence over the points at which the field takes low values, as it happens in the boundary of the beam, rather than in the proximity of the beam center, where the field maintains relatively high values.To address this problem, it is a common practice to express the spatial covariance of the field in the analytical approach as follows The circular symmetry of the beam allows to consider the particular case , which is intended to be representative enough of the whole beam statistics.This mathematical treatment is commonly accepted, although it contains the implicit assumption of an isotropic turbulence.Against this approach it could be argued that the turbulence may present a certain degree of anisotropy, firstly due to the specific direction of the wind (which is not radial at all), and secondly, introduced by the proximity to the ground in horizontal links [23].
All the points mentioned in sections 2.2 and 2.3 justify the interest of looking for an alternative tool, capable of giving an insight about the dynamics of the laser beam propagation in turbulent media.In our case, we have resorted to a specifically designed numerical technique to simulate the wind effect in the atmospheric turbulence.This technique is described in the following section.

Practical generation of large phase screens
Under the Taylor hypothesis the movement of the atmosphere due to the wind can be simulated in a numerical experiment merely by displacing the phase screens in a transverse direction.The technique would then be as simple as repeating the complete propagation algorithm, from the transmitter to the receiver, and shifting the phase screens one or several grid points each time a propagation sequence is performed.However, due to the large computational requirements in using large phase screens, it is not clear in practice how to implement this movement while keeping realistic conditions.This is so because a minimum screen length is required to obtain statistically relevant results.To give an idea about the practical difficulty of this direct approach, we can imagine a not very ambitious numerical experiment in which we were using a 512×512 grid of 0.8 m width, which is supposed to be enough for describing both the laser beam profile and the turbulence along the complete path.To simulate the wind motion by moving an equivalent rectangular screen in the transverse direction, 512×32000 grid points would be needed to generate a 10 second measurement event under a wind of 5 m/s.If we suppose that 50 of such phase screens are required to simulate the whole propagation length, 6.5Gbytes of memory would then be necessary only to allocate the screens.A further restriction is that the minimum screen length should at least double the outer scale value L 0 to properly comprise the largest eddies of the turbulence.Typical measured values of L 0 in experimental campaigns extend between 10 and 20 meters [24].Hence, a reasonable screen length should be between 20 and 40 meters in the shift direction.
To overcome this problem, we have partitioned the long rectangular screens into several (square) ones as follows.First, a low-resolution rectangular screen is generated, by means of the covariance method [18].This initial screen must cover the whole time interval being simulated and will be stored in the computer memory during the whole simulation process.Then, the initial screen is divided into several blocks, and a high resolution screen is constructed from each block by means of an interpolation technique [16].To achieve the desired resolution, the interpolation is done in several consecutive steps.The inner grid points are obtained by using a nearest four point interpolation, whereas the points at the edges are interpolated from the nearest two edge points.To minimize the error introduced at each interpolation step, we have followed the procedure proposed by Recolons et al. [25].Finally, the screen is moved in a direction transverse to the beam propagation, thus simulating the wind.When the beam spot reaches the edge of the screen, a new high-resolution screen is generated by interpolating the next low-resolution block.The method used here differs essentially from that proposed by Assémat et al. [9], in the sense that we do not perform extrapolations of successive columns from an initially exact screen, but rather make a series of interpolations over a basic low-resolution frame.Although in our case it is not possible to create an arbitrarily long screen, due to the huge memory requirements inherent to the covariance method, we have preferred this approach, in order to ensure that a proper correlation is maintained along the whole phase screen in both dimensions.Fig. 3. From the exact low-resolution phase-screen (above in the figure) a certain number of columns is taken to create a smaller high-resolution screen (below).L y is the length of the initial screen in the moving direction.
The whole procedure is sketched in Figs. 3 and 4. Figure 3 shows the initial low-resolution screen.The stars represent the exactly calculated points, by using the covariance method.From this screen we select the first columns on the left, as shown in the figure.After several interpolation steps, we obtain the high-resolution screen of Fig. 4. The final high-resolution block may be considered as a series of square screens (three in the example of the figure) tied together.The main difficulty of the method resides in the fact that the points laying in the vicinity of the right end of the block (we are assuming here that the screen moves to the left) cannot be properly interpolated, for at this moment we are ignoring the information from the first columns of the next block, so they must be calculated later.When the beam spot reaches the third square screen on the right side, the first square screen on the left is chopped off, and a second section from the low-resolution screen is added to the right.The new section is interpolated in the same way as before, together with the points that had been left uncalculated in the previous stage.This procedure can be repeated as many times as necessary, until the end of the low-resolution screen is reached.By using this procedure, only a small part of the entire high-resolution screen is kept in the computer memory at a time.The saving in memory storage allows dealing with phase screens that are large in size, with the only limitation of the maximum number of low-resolution points that can be calculated by the covariance method.In our case, we departed from a 5×200 low-resolution grid and performed between seven and nine consecutive interpolations to achieve the desired resolution, depending on the case.
Figure 5 shows the effectiveness of the proposed technique in generating possible continuations for a particular screen.The two different phase screens in the figure can be considered as a combination of three square screens, tied together, which have been generated departing from the same low-resolution grid.The leftmost section is identical in both screens, whereas the central and right sections have been interpolated independently.We can see that in addition to the correlation conditions, continuity between sections is fully satisfied.In Fig. 6 we plot the structure function obtained in the horizontal direction from 512×1024point screens generated with this technique.The results are compared with the theoretical Kolmogorov curve.We can see that the match is reasonably good, being the differences attributable to the error introduced by the interpolation process.We must note that even a small departure from the expected behaviour may affect the statistical behaviour of the simulated turbulence.Therefore minimizing the error committed is of crucial importance.In our case, we used the optimization technique of Ref. 25. Increasing the density of the lowresolution grid, thus reducing the number of interpolation steps needed, also results in a better accuracy, with the only limitation being the computing resources required by the covariance method.

Results
For the sake of simplicity and in order to avoid the not yet well known influence of beam wander in horizontal links, we have limited ourselves to the analysis of the temporal evolution of the irradiance at single points in the reception plane.A further advantage of assuming a point detector is that there are theoretical models available to compare with our simulations.
To investigate the temporal behaviour of the laser beam irradiance, by means of the temporal covariance analysis, we have selected three different situations that correspond to low, moderate and high turbulence strength, respectively.The numerical experiments have been worked out for a Gaussian laser beam with W 0 = 2 cm and λ = 1.064 microns.The turbulence structure constant considered is C n 2 = 10 -13 m -2/3 in all three cases, whereas the difference between them comes from the different path lengths chosen.The results obtained are compared with the previously published theoretical Gaussian beam model for the case of weak turbulence [6].To this date, no other theoretical model has been proposed for Gaussian beam propagation in the cases of moderate and strong turbulence.
Several covariance functions and power spectra for the normalized irradiance have been calculated in each case.First of all, we have obtained the covariance of the irradiance at the optical axis (r = 0), i.e., the center of the unperturbed beam, and at a distance from the optical axis equal to half of the long-term beam width (r = W LT /2).We have also considered a second reference point in the instantaneous field, the centroid or 'center of mass' of the irradiance.Both reference points are eventually the same in the long-term beam, but due to wander effects, they may differ in the short-term beam.Consequently, we have also examined the temporal evolution of the irradiance at the centroid and at a W LT /2 distance from it.
The effective or long-term beam width is defined as the radius at which the average intensity decays by a factor 1/e 2 with respect to the maximum, according to the expression which holds at least for weak turbulence.In the simulations, the long-term beam radius has been set as W LT = r 86 , i.e., the radius of the beam containing the 86% of the total power.

Gaussian-beam case: Weak turbulence
The path distance chosen in the first simulation is L = 500 m, which yields a value σ R 2 = 0.867 for the Rytov variance.(We take here the usual criterion σ R 2 = 1 as the limit between weak and moderate turbulence).To simulate the wind motion, a set of low-resolution rectangular phase screens have been generated, covering each of them a distance of 30 m in the transverse direction.By using the partition technique described in section 3, these screens have been decomposed into 256×256 square phase screens, each of them corresponding to a working window with dimensions L x = L y = 0.174 m.The grid interval is 0.68 mm, thus giving an effective inner scale l 0 = 1.36 mm.The wind motion was simulated by moving the screens horizontally 4 grid points at each step from the previous position.The whole simulation has been repeated ten times under the same conditions.The temporal covariances obtained for the irradiance at the optical axis and at the beam centroid are shown in Fig. 7(a).The numerical results are compared with the theoretical results given by the approximated analytical expressions in Ref. 6, chap.8, which can be written in a simplified way for collimated beams, as follows where In the figure, the time variable τ has been normalized by multiplying it by the angular frequency ω t = V•(k 0 /L) 1/2 , which corresponds, in first approximation, to the bandwidth of a plane wave affected by weak turbulence [3], thus giving a frequency f t = ω t /2π = 17.3 Hz for a wind velocity V = 1 m/s in the transverse direction.As seen in the figure, the covariance measured at the beam centroid and at the optical axis lead to somewhat different results in the weak turbulence regime.The discrepancies observed may be attributed to the beam wander effect.There is also a good agreement with the theoretical prediction given by Eq. (11).
The corresponding spectra are shown in Fig. 7(b).The numerical results have been obtained directly by Fourier transforming the covariances of Fig. 7(a), whereas the theoretical spectrum has been obtained by using the approximate expression given in ref. 6.The plots in Fig. 7(b) show again an acceptably good agreement between theory and simulations at low frequencies.However, the spectral decay at higher frequencies in the theoretical spectrum is faster than in those obtained in the simulations.This discrepancy may be attributable to the approximations made in Ref. 6. Nevertheless, it is worth noting that the cut-off frequency of the irradiance spectrum still agrees with the theoretical prediction.In the case considered in Fig. 7, this frequency moves between 1.5ω t and 2ω t .
The covariance and spectra for the normalized irradiance obtained at a distance r = W LT /2, both from the optical axis and from the beam centroid are plotted in Fig. 8.The numerical values for the irradiance covariance have obtained by averaging the covariance measured at four different points, all of them placed at a distance 0.5 W LT from the reference point considered, and located at the extremes of an imaginary cross.If there were some nonhomogeneity present in the wind-generated turbulence, this would lead to different covariance values, depending on whether the point of measurement is located in a direction parallel or perpendicular to that of the wind.Although in some cases there seemed to exist some difference in the results obtained, it was not significant enough to give definitive conclusions on this subject.For this reason, we omitted these discrepancies, and restricted ourselves to the study of the ensemble average covariance.The obtained results are plotted in Fig. 8(a).Here, the effect of beam wander is more noticeable than in Fig. 7, which is reflected by the lower covariance values obtained when taking the centroid as the reference point.Nevertheless, the covariance values measured at the optical axis still agree with those predicted by theory.
The corresponding spectra are plotted in Fig. 8 (b).We must note that the theoretical spectrum is plotted only for angular frequency values ω > ω t .The reason is found in the fact that the analytical expression of Eq. ( 11) for the irradiance covariance includes the Kummer hypergeometric function 1 F 1 , which is difficult to evaluate for large values of the argument.
Thus, it is difficult to obtain good theoretical values for τ • ω t > 4. As a consequence, the Fourier transform obtained numerically from the theoretical covariance gives no values at very low frequencies.Unfortunately, there is no approximate analytical expression for the theoretical spectrum that could overcome this problem for the case r ≠ 0. Nevertheless, the agreement between theory and simulation is still good within the cut-off region, as seen in the figure.

Gaussian-beam case: Moderate turbulence
The path length chosen to simulate the moderate turbulence case was L = 1500 m, which corresponds to a Rytov variance σ R 2 = 6.49.The effective grid size was 512×24576 points for the rectangular phase screens, with a total screen length of 30 m in the shift direction.The grid interval now is 1 mm, thus resulting in an equivalent inner scale l 0 = 2.0 mm.The temporal covariances and the power spectral densities of the normalized irradiance, obtained both at the optical axis and at the beam centroid, are shown in Fig. 9.The normalized irradiance spectra obtained at a distance r = 0.5 W LT from the two reference points is also depicted in Fig. 10.The value of the normalizing frequency at the horizontal axis is now f t = 10 Hz for a wind velocity V = 1 m/s.As one might expect, the contribution of beam wander effects to the on-axis covariance irradiance is much greater than in the weak turbulence case, which can be seen by comparing the plots in Fig. 9(a) with those in Fig. 7(a).For the sake of comparison, the theoretical covariance and power spectrum, obtained by assuming that the weak turbulence model was still valid, has also been plotted.The difference in behaviour between the theoretical prediction and the numerical results is now more apparent than in the weak turbulence regime (see Fig 7), and is affecting the shape of the curves.Here, the departure of the numerical results from the theoretical prediction given by the Rytov theory is not only quantitative but also qualitative.Furthermore there is no theoretical expression for the temporal covariance at the centroid.Oppositely to the weak turbulence case, the theoretical prediction for the cut-off frequency gives poor results, as seen in Figs.9(b) and 10.Fig. 10.Irradiance spectral power density curves obtained for the case of moderate turbulence, at points located at a W LT /2 distance, both from the optical axis and from the beam centroid.

Gaussian-beam case: Strong turbulence
The last case analysed for laser beam propagation is that of a strong turbulence.The results plotted in Fig. 11 have been obtained by setting σ R 2 = 23.1.This value of the Rytov parameter corresponds to a propagation distance L = 3000 m, keeping the same value C n 2 = 10 -13 m -2/3 for the structure constant as in the previous cases.The overall size of the rectangular phase screens was 1024×29696 grid points, corresponding to a 30 m length in the wind direction.Figure 11 shows the power spectra of the irradiance obtained at the two reference points considered, i.e., the optical axis and the beam centroid, and at a distance W LT /2 from them.The normalizing frequency is now f t = 7 Hz.Unfortunately, no theoretical spectra can be provided for comparison in this case.As a reference, we have re-plotted instead the curves of Fig. 8 for moderate turbulence.It should be noted that the difference between the spectrum at the optical axis and at the beam centroid in Fig. 11(a) tends to be smaller in the case of strong turbulence.We attribute this effect to the breaking of the beam that occurs at such turbulence conditions.Another interesting point is that the normalized bandwidth of the spectra does not change significantly from the weak turbulence case, giving a value between 3ω t and 4ω t for the spectrum of the irradiance at the optical axis and somewhat lower at the centroid.
Moreover, since the normalizing frequency ω t decreases, the real bandwidth remains nearly constant.In the three cases considered, the cut-off frequency may be estimated to be around 25-35 Hz assuming a transverse wind velocity of 1 m/s.

Plane wave: Strong turbulence
Finally, we have studied the case of a wide super-Gaussian beam, propagating in conditions of high turbulence, in order to simulate the behaviour of a plane wave.The covariance and spectra obtained are plotted in Fig. 12.The super-Gaussian beam launched is of the form with a 20 cm initial width.The propagated distance was L = 5000 m and the Rytov variance σ R 2 = 59.The grid size is 1024×15360 points for the rectangular phase screens, with the total length being 40 m.We performed twelve independent experiments with the same conditions.
The plots in Fig. 12(a) show that the behaviour of the irradiance at the optical axis (r = 0) and at the beam centroid is nearly identical, as expected on a plane wave (for there is no wander at all), with a covariance approaching that of an impulse function.The saturation effect in the scintillation is made evident by the low covariance value at τ = 0.For the sake of comparison we have also plotted in the figure the theoretical values derived from the expressions in Ref. 6, Chap.9. Noticeably, the shapes of the curves differ appreciably from those obtained in the simulations.

Conclusions
The loss of the Gaussian shape in the long-term irradiance under moderate and strong turbulence poses an additional difficulty in the temporal analysis of a laser beam propagating into the atmosphere.Hence the need of a numerical simulation of the temporal phenomena.To this end, we have adapted the covariance method introduced by Harding et al. to generate long phase screens.In spite of the error introduced by the interpolation process, the phase screens obtained exhibit a reasonably degree of accuracy.With these phase screens we have simulated the temporal behaviour of a laser beam propagating along a horizontal path into the atmosphere under weak, moderate and strong turbulence conditions.The error committed in the temporal covariance amplitude has been estimated to be 8% in the worst case.The results obtained in the simulations show a good agreement with the theory in the weak turbulence case.The effect of beam wander in the covariance irradiance of a Gaussian beam has also been put on evidence.This effect is greater in the moderate turbulence case, but reduces again in the strong turbulence case.The contribution of beam wander to the temporal covariance is also noticeable in the weak turbulence regime if we place the detector at a distance W LT /2 from the beam peak.In general, the effect of beam wander in the irradiance behaviour is always greater at a distance W LT /2 from the beam peak than at the peak itself, no matter the value of the Rytov variance.Our results also show that for a Gaussian beam the real bandwidth of the power spectra at the optical axis remains nearly constant, independently of the Rytov variance.In the three cases considered the cut-off frequency may be estimated to be around 25-35 Hz for a transverse wind velocity of 1 m/s.We have also seen that the bandwidth of the irradiance fluctuations increases only slightly at points located at a certain distance outside the center, with respect to the value obtained in the peak.
We have also simulated the propagation a plane wave under strong turbulence conditions.The discrepancies found between theory and simulations are remarkable.
Additionally, in some cases we have observed slight discrepancies in the covariance values measured at points placed in a direction perpendicular to the wind, with respect to the values obtained at points placed along the wind direction.These small discrepancies may be attributed to the inhomogeneities caused by the wind.Unfortunately, we were not able to provide significant results at this moment, with the sample sizes employed.Further research is needed to give definitive conclusions on this issue.Also, future work is planned to include the use of a finite aperture in the reception plane.
results are shown in Fig. 1.In the case of weak turbulence (C n 2 = 10 -16 m -2/3

Fig. 1 . 2 = 10 -
Fig. 1.Progressive loss of the Gaussian shape in the long term analysis of a laser beam propagating in the atmosphere.This effect is negligible under weak turbulence conditions (C n 2 = 10 -16 m -2/3).The deviation of the field profile from the Gaussian shape increases with the accumulated turbulence along the beam path.Dashed black lines have been obtained by using Eq.(6), dot-dashed blue lines have been obtained by using Eq.(7), and red solid lines have been obtained by using Eq.(8).The grid parameters correspond to an effective inner scale value l 0 = 2.4 mm.

Fig. 4 .
Fig. 4. Generation of the high-resolution screen by interpolating the low-resolution positions (black stars).The red points on the right cannot be properly calculated.The blue ones are correct and they will have to be taken into account when interpolating the next block to ensure continuity (below).

Fig. 5 .
Fig. 5. Example of two different phase screens interpolated from the same low-resolution grid.The leftmost section is identical in both screens.The figure shows two possible continuations in the middle and right sections, which have been generated independently.

Fig. 6 .
Fig. 6.Structure function D(y) in the horizontal direction Y for screens calculated with the method presented in this work (blue line).The black line represents the theoretical Kolmogorov structure function D = 6.884 (y/r 0 ) 5/3 , where y is the coordinate in the wind direction, and r 0 is the Fried parameter

Fig. 7 (
Fig. 7 (a) Temporal irradiance covariance and (b) spectral power density of the irradiance at the reception plane obtained at a propagation distance L = 500 m from the transmitter, in the weak turbulence regime.The curves in black represent the theoretical covariance and spectrum at the optical axis.In both figures, the horizontal axis has been normalized with respect to the theoretical cut-off frequency ω t for plane waves.

Fig. 8 .
Fig. 8. (a) Temporal covariance and (b) spectral power density of the irradiance in the same case as Fig. 7 but for points situated at a distance W LT /2 from the two centers considered in the beam.

Fig. 9 .
Fig. 9. (a) Temporal covariance and (b) power spectrum of the irradiance at the reception plane, after a propagation distance L = 1500 m, obtained both at the optical axis and at the beam centroid, in the moderate turbulence regime.The normalizing frequency is f t = 10 Hz for a 1 m/s wind velocity.

Fig. 11 .
Fig. 11.(a) Irradiance spectral power density obtained at the theoretical center of the beam and at the centroid in the case of strong turbulence.Gray and green curves correspond to the moderate turbulence case.(b) Id. at a distance W LT /2 from the two previous points.

Fig. 12 .
Fig. 12.(a) Temporal covariance and (b) power spectrum of the irradiance for a plane wave, simulated by a super-Gaussian beam, in regime of strong turbulence.The difference between the behaviour at the center (r = 0) and at the centroid becomes negligible.

Table 1 :
Several radius definitions for a Gaussian beam as a function of the power fraction encompassed by circle centered at the peak of irradiance