Propagation of low-quality laser beams in turbulent air: new fast simulation methods and experimental results

Novel inherently fast methods basing upon fundamental Gaussian modes are presented to numerically simulate the propagation behavior of a low-quality laser beam, i.e., a beam characterized by a high quality factor M2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M^2$$\end{document}, in turbulent air. Actually, for a given initial diameter, beam wandering and spread at the target position is calculated. Test results are checked against respective simulations using higher-order mode beams having clearly defined quality factors and nearly perfect agreement is demonstrated. Experiments performed on the German Aerospace Center (DLR) laser test range at Lampoldshausen near Stuttgart satisfactorily accord with the predictions of the proposed simulation techniques.


Introduction
High-power lasers for long-range applications are often still seriously affected by inherently poor beam quality characterized by beam quality factors M 2 beyond 5. Typically, even though the overall beam (near field) intensity profiles appear similar to (super) Gaussians, they are commonly joined with very high quality factors. Modeling the propagation of such a beam with given diameter at the laser output plane and quality factor, one is faced with the fact that the angular 1 3 271 Page 2 of 8 scaling of the turbulence strength by a certain power of the beam quality factor is needed, it can only be used for beams described by a single M 2 -value. In the second procedure called method of effective focal length, a defocusing lens is introduced in order to account for the additional diffractive spread of a beam with quality factor greater than unity. As there is no need for any manipulation of the turbulence module in the computer code, this method is most flexible and can be used for any beams represented by two quality factors M 2 x and M 2 y . Finally, as suggested by the numerical results of these two methods, a third, (semi-analytical) "super-fast" procedure for the identification of the final beam diameter is introduced.
The paper is organized as follows: in "Theoretical basis of laser beam propagation in turbulent atmosphere" section, a compact introduction to the theoretical basis of laser beam propagation in a turbulent atmosphere is given, followed by the analytical foundation of the two methods for laser beam quality emulation. "Computations with reference to the propagation of Laguerre-Gaussian beams with distinct quality factors" section presents the results of test computations with the Code TALAP (developed at the Stuttgart DLR Institute of Technical Physics) for both methods with reference to propagation simulations of higher-order Laguerre-Gaussian beams with distinct quality factors. Results of the semi-analytical (third) method conclude this section. Finally, results of propagation experiments done at the DLR laser test range (135 m) at Lampoldshausen near Stuttgart together with conclusions from respective simulations are presented in "Experimental results and numerical simulations" and "Conclusions" sections.

Theoretical basis of laser beam propagation in turbulent atmosphere
The laser beam propagation code TALAP ("Turbulent Atmosphere and Laser Beam Propagation") solves the paraxial wave equation for the (scalar) laser field amplitude u(x, y, z), see e.g. Ref. [8], k = 2 ∕ denoting the laser wave number, by a Fast Fourier Transformation (FFT) method. For the propagation of laser radiation (cw or pulsed) of medium power density in air, thermal blooming plays a minor role. Hence, a stationary, homogeneous and isotropic turbulence flow field may be assumed. Here, the Kolmogorov turbulence model applies stating that energy is fed into the flow at large scales and cascades down at a constant rate per unit mass to small scales where it is dissipated by molecular viscosity, e.g. see Ref. [9]. Hence, between a so-called outer scale length L 0 and inner scale length l 0 inertial forces dominate viscosity and refractive index fluctuations due to the turbulent flow field are characterized by the von Karman power spectrum where C 2 n is the refractive index structure constant with unit [m −2∕3 ]. The spatial wave number is given by 2 = 2 x + 2 y + 2 z , while 0 = 2 ∕L 0 and m = 5.92∕l 0 . Consequently, for a given laser wavelength and propagation distance L, the turbulence problem is governed by the three free parameters C 2 n , L 0 and l 0 . An important quantity is the Fried parameter r 0 (coherence length) [10], scaling as 6∕5 ∶ Physically, the Fried parameter may be regarded as the aperture of a fictive telescope over which the rms phase differences between any two points remain in the order of one radian. That means, a telescope with aperture r 0 in the absence of turbulence (i.e., diffraction-limited) has the same image resolution as any telescope under turbulence conditions having an aperture arbitrarily larger than r 0 . In our code, the numerical turbulence evaluation is based upon the phase-screen approach. Here, the propagation range is divided into a number of segments Δz > L 0 and the cumulative phase shift of any segment is considered at its upper boundary. The cumulative phase shift kΓ(x, y) is given by the perturbation of the optical path length where the local fluctuation of the refractive is denoted by n(x, y, z). As shown in Ref. [11], the optical path length perturbation is derived as with Φ n ( x , y ) = Φ n ( x , y , 0). The complex random function a( x , y ) must have the property a( x , y ) = a * (− x , − y ) because Γ(x, y) is a real quantity. Real and imaginary parts of a( x , y ) are independent random numbers with zero mean value and variance 1/2. For the numerical realization of the phase screen, a Fast Fourier Transform (FFT) method is applied to evaluate Eq. 5. To account for very different orders of magnitude of outer and inner scales (e.g. L 0 in the meter-and l 0 in the millimeter range), 0 = 2 ∕L 0 is not explicitly considered in Eq. 2. As indicated in Fig. 1, two phase screens are computed instead, one with dimensions of the computational area for the beam propagation several centimeters squared and the other with large dimensions L 0 × L 0 [12]. The screens are constructed in a way that the lowest spatial wave number of the small screen matches the highest spatial wave number of the large screen. Phases of the large screen, after interpolation (colored square area in Fig. 1), are then added to the phases of the small screen. This procedure is a kind of technique called "addition of subharmonics" [13]. The method has been successfully tested with reference to the theoretical phase structure function , r = |⃗ r| and r 0 denoting the distance between any two points in the flow field and the Fried parameter, respectively. The brackets indicate an average over an ensemble of phase screens.
In this paper, we concentrate on untwisted ("simple astigmatic") laser beams with principal axes being unchanged along the direction of propagation. Their propagation behaviors are completely defined by the 4 -diameters 2W 0x and 2W 0y based upon the second moments of the intensity distributions I(x, y) at respective waist positions z 0x and z 0y , as well as by the respective beam quality factors M 2 x and M 2 y [7]. The second moment 2 x (variance) of the intensity profile is given by where the first moment of the intensity profile represents the center of gravity position (x-centroid) of the beam; analogous relations are valid for 2 y . The definition of the second moment (variance) based beam diameter aims at the short-term intensity profile at the target position, i.e., the effect of beam wandering is eliminated. In contrast, the long-term diameter follows as the square root of the mean value for an ensemble of propagation runs of the expression . The beam quality factors are defined via the waist parameters as M 2 i = W 0i i ∕ , where i denote the respective far-field divergence half angles and i = x, y.
In the following, two methods for the emulation of lowquality laser beams in a turbulent atmosphere are presented.

Extended method of embedded Gaussian
As is well known, the second moment based half diameters W i (i = x, y) of a beam with quality factors M 2 i are rigorously governed by the propagation equations Defining new half diameters as w i = W i ∕M i , the propagation equations now read This means that a fundamental Gaussian beam (M 2 i = 1) with half diameters w i behaves in the same optical fashion like a beam with quality factors M 2 i having the diameters 2W i = 2w i M i [7]. Hence, to emulate the propagation of a beam with quality factors M 2 i , one has to start at the laser output plane with a Gaussian intensity distribution and the respective lateral dimensions downscaled by the factors M i . After propagation, the lateral dimensions have to be rescaled again at the target plane. Combining the method of embedded Gaussian with the turbulence phase screen approach, one has to consider that the ensemble averaged turbulence induced part of the variance of the focal (far field) intensity profile scales as L 2 ( ∕r 0 ) 2 . This is easily understood by the physical interpretation of the Fried parameter given above (Eq. 3). In order to obtain the correct turbulence spread contribution for the beam diameter when rescaling the lateral dimensions by M at the end of the propagation distance, according to Eq. 3, the refractive index structure constant C 2 n (7) has to be downscaled before propagation by the factor M −5∕3 . However, as the mean square of the beam centroid displacements ⟨x 2 c + y 2 c ⟩ (beam wander) is proportional to C 2 n R −1∕3 0 L 3 (R 0 denoting the initial beam radius) [14,15], the computed rms value of the beam centroid displacement at the target position must be rescaled by M 2∕3 and not by M, thus taking into account the additional effect of R −1∕3 0 being scaled up by M 1∕3 before propagation. A more detailed outline of this scaling relations can be found in the appendix. In order to successfully apply this procedure, we must be able to characterize the beam by only one M 2 , i.e., M 2 x ∼ M 2 y = M 2 .

Method of effective focal length
Within the framework of the complex Gaussian beam parameter formalism [8,16] it is easy to demonstrate that the variance of the far-field divergence of a beam with initial diameter 2w 0 and quality factor M 2 equals the one of a fundamental Gaussian beam (M 2 = 1) with the same initial diameter being modified by a defocusing lens with a certain focal length F < 0. The 4 -diameter of the actual laser beam (being represented by the respective embedded Gaussian, cf. "Extended method of embedded Gaussian" section) at the output plane is denoted as 2w 0 . With zero curvature of the phase front, its beam parameter q 0 at z = z 0 = 0 follows as or For the fundamental Gaussian beam emulating the actual laser beam we have or As in general q(z) = q(0) + z, at the beam position z ≫ z Rayleigh = w 2 0 ∕ , the inverse complex radii of curvature 1∕q and 1 / q, respectively, follow as and (10) For equal far-field divergence angles and ̃, we have to demand, because ≃ w∕z and ̃≃w ∕z, that w(z) ≃ w(z). The two squared beam half diameters are found by multiplying the inverse imaginary parts 1∕ℑ(1∕q) and 1∕ℑ(1∕q) by M 2 ∕ and ∕ , respectively. Equating w 2 ∕z 2 and w 2 ∕z 2 leads to or Now, if a focused (focal length f) propagation of the two beams covering the distance L = f is considered, the spot diameters at the target will be identical, because the focusing lens (mirror) creates the far-field image. Hence, at the target, a fundamental Gaussian beam and a focusing device with optical power D eff = 1∕f eff = 1∕F + 1∕f emulate the 4 -spot-diameter of a beam having the same initial diameter and a quality factor M 2 . The effective focal length follows as In the case of an untwisted beam with principal axes in xand y-direction and respective beam quality factors M 2 x and M 2 y , the procedure can be separately applied to either direction using an astigmatic defocusing lens with two respective radii of curvature. As this method does not involve any scaling of the lateral dimensions, there is no need for any intervention in the turbulence module. Hence, it can be applied to any simple astigmatic laser beam as defined above. Numerical simulations on the basis of Gaussian beam theory or 3D-solutions to the paraxial wave equation confirm this result with very high accuracy, cf. "Computations with reference to the propagation of Laguerre-Gaussian beams with distinct quality factors" section.

Computations with reference to the propagation of Laguerre-Gaussian beams with distinct quality factors
In this section, the results of numerical simulations of the propagation of a focused laser beam (λ = 1.03 μm) over a distance L = 2000 m are presented. In Fig. 2 the short-term beam diameter (rms of the 2 x 2 c + y 2 c in the ensemble) are shown at the target position as functions of the (initial) beam quality factor M 2 for the two proposed methods with reference to pure Laguerre-Gaussian (TEMp0) beams with quality factors M 2 = 2p + 1. T h e t u r b u l e n c e p a r a m e t e r s a r e C 2 n = 5.0 × 10 −13 m −2∕3 , l 0 = 8 mm and L 0 ≈ 3.6 m. All propagation runs start out with initial intensity distributions having the same 4 -diameter of 2R 0 = 8 cm. The computational area of 120 cm × 120 cm is sampled with a resolution of 2 11 × 2 11 mesh points. Beam diameters and respective centroid standard deviations versus refractive index structure constant C 2 n for the same parameters but a fixed value M 2 = 9 are displayed in Fig. 3. It is evident that the results of both novel computational methods coincide with those of the reference Laguerre-Gaussian simulations in a nearly perfect manner.

Super-fast semi-analytical method
The numerical results in Fig. 2, indicating a slightly parabolic behavior of the diameter and a nearly constant centroid deviation, encourage for trying the following semianalytical approach for the evaluation of the beam diameter D (short-term or long-term) after focused propagation over a distance L in turbulent air: Here, D 2 1 denotes the computed square of the second moment-based (short-term or long-term) diameter of a fundamental Gaussian beam at the target position. That means, D 2 1 is proportional to the sum of the Gaussian beam variance due to diffraction and the turbulence-induced variance including or not including beam wander. The term 4( L∕ R 0 ) 2 (M 4 − 1) represents the spread due to diffraction of the higher order beam modes; again, R 0 denotes the initial beam radius. The relation (19) is easily transformed to M 4 total = 1 + M 4 ab + M 4 turb , where M 4 ab and M 4 turb are due to initial higher order modes in the beam (aberration) and turbulence, respectively [6].
With a beam diameter D 1 = 440.0 mm for M 2 = 1 (cf. Fig. 2), the numerical evaluation of Eq. 19 is shown in Fig. 4, again, with reference to pure Laguerre-Gaussian beams. We conclude from the nice agreement that, in this turbulence regime of medium strength (Rytov parameter 2 R = 0.63( L∕r 2 0 ) 5∕6 ≈ 57), inherent diffraction of the beam due to poor initial quality and turbulence-induced spread are statistically fairly independent processes, i.e., the respective variances simply add. In practice, for a given turbulence state, this semi-analytical third method provides "super-fast" information about the relevant beam properties in the focal plane for any value of initial M 2 .

Experimental results and numerical simulations
For the experimental analysis of the dependence of the beam propagation on weather conditions, the DLR operates a free transmission laser test range at Lampoldshausen, Germany. It consists of two stations, the transmitting (TS) and the receiving (RS) one confining a 135-m-long pathway with a beam path 1 m above asphalt ground. Optical turbulence has been measured by a surface layer scintillometer (SLS 20-A, Scintec AG, Germany). A more comprehensive description of the optical test range including installed sensors continuously monitoring the local atmosphere can be found in Ref. [17]. The simultaneous characterization of the laser beam is performed with sensors inside of the TS and RS. These measurements address power, intensity distribution, and jitter of the beam.
The laser system used is a TruDisk 6001 (4C) disk laser (Trumpf, Germany). It operates at a wavelength of 1.03 μm and achieves an adjustable output power (continuous wave) between 180 and 6000 W. The beam quality factor M 2 of the laser system is measured as about 10 at 180 W and increases with power. The collimated beam is guided into a telescope setup consisting of two off-axis parabolic mirrors with respective focal lengths of 227 and 2272 mm and a magnification of 10. The distance between the telescope mirrors is adjusted to obtain the focus point at 135 m distance.
In Figs. 5 ( = 0.532 μm, M 2 = 3) and 6 ( = 1.03 μm, M 2 = 18), comparisons between experiment and simulation are shown with respect to the 4 -beam diameter and the beam centroid standard deviation as functions of the refractive index structure constant C 2 n . Clearly, particularly in Fig. 6, the very poor initial beam quality dominates the turbulence-induced effect on the beam diameter at the target position. In both cases, however, the increase of the 4 -beam diameters with growing C 2 n -values as predicted by the simulation is reflected by the measurement. The experimental and computational absolute values of the diameters also agree very well. Additionally, the measured beam centroid standard deviations exhibit satisfactory consistency with the theoretically expected values.

Conclusions
In this paper, we concentrated on the problem of numerically emulating laser beams propagating in the turbulent atmosphere with initial quality factors M 2 far beyond unity. With the extended method of embedded Gaussian and the method of effective focal length, two numerical techniques were introduced relying on (focused) Gaussian beams which yield the laser intensity distribution at the target position in a fast and flexible manner. As far as the simulation of the diffractive spread due to the initial beam quality is concerned, we presented rigorous analytic foundations for both methods. Hence, with atmospheric conditions admitting a description within the framework of a Kolmogorov turbulence model, both procedures generally apply. A semi-analytical procedure for a first (approximate) determination of the laser focal spot diameter completes our theoretical proposal. These novel numerical tools are ready now to provide support, for example, to investigations of the relative importance of beam quality and air turbulence for laser beam propagation. As a first application, we simulated propagation experiments done at the DLR 135 m laser test range and found good agreement with respect to beam diameter and centroid standard deviation. Experiments with two or three paths along the 135 m laser test range would further increase the turbulence-induced effects and are intended for the future together with applications of the presented tools to more experimental data.