Method of Spherical Phase Screens for Modeling the Propagation of Diverging Beams in Inhomogeneous Media

The phase-screen (split-step) method is widely used for modeling wave propagation in inhomogeneous media. The method of plane phase screens is best known. However, for modeling a 2D problem of radio occultation sounding of the Earth’s atmosphere, the method of cylindrical phase screen was developed many years ago. In this paper, we propose a further generalization of this method for the 3D problem on the basis of spherical phase screens. In the paraxial approximation, we derive the formula for the vacuum screen-to-screen propagator. We also infer the expression for the phase thickness of a thin layer of aisotropic random media. We describe a numerical implementation of this method and give numerical examples of its application for modeling a diverging laser beam propagating on a 25-km-long atmospheric path.


INTRODUCTION
The phase-screen method has long been widely used for the numerical modeling of the propagation of waves of various natures in inhomogeneous media. As examples, we indicate the problems of modeling the propagation of optical (laser) radiation in a turbulent atmosphere [1][2][3][4][5] and the problems of modeling the propagation of decimeter waves during radio occultation sounding of the atmosphere [6][7][8]. This method in the English literature is called split-step, a method of stepby-step splitting. Accordingly, it is called the method of splitting by physical parameters in the Russian literature. This name reflects the fact that the entire inhomogeneous medium in this method is represented as a sequence of thin layers, and the propagator describing the propagation of a wave through each layer is approximately written as a composition of an infinitely thin layer forming phase distortions of the wave and a vacuum propagator describing diffraction.
The fundamental limitation of the phase-screen method is that it does not take into account backscattering. The phase-screen method can also be considered a finite-dimensional approximation of the path integral describing wave propagation in an inhomogeneous medium [9].
2D phase screens are used in modeling laser radiation in turbulent media, especially to describe the effect of isotropic turbulence. However, significant optimization of computational costs is achieved through the use of 1D phase screens in the modeling of radio occulta-tion experiments. This is due to the fact that, in most cases, atmospheric inhomogeneities with vertical scales from hundreds of meters to kilometers are taken into account, the horizontal scales of which significantly exceed the horizontal size of the Fresnel zone.
In the classical version of the method, plane phase screens are used. This leads to unnecessary computational costs in the case of describing a diverging wave, since at the edges of each screen the angle between the screen and the wave front increases, which results in the need of reduction of the discretization step. In the 2D modeling of radio occultation experiments, it turned out to be quite simple to write down a solution for cylindrical 1D phase screens [7], which allows one to take into account the shape of the phase front of the incident wave.
In this paper, we generalize this approach and develop the method of spherical phase screens. Since it is not possible to derive a simple exact solution in this case, we use a paraxial approximation. In the section Derivation of Basic Relations, we apply the technique of angular spectra [10] and derive the formula for the propagator in spherical coordinates in the paraxial approximation. In the Numerical Modeling section, we show that the numerical implementation of this method is no more complicated than the case of plane phase screens, and we give examples of numerical modeling using the isotropic turbulence model.

DERIVATION OF BASIC RELATIONS
Vacuum Propagator As an example, we will use the beam parameters specified for the DELICAT project (DEmonstration of LIdar based Clear Air Turbulence detection [11,12]: wavelength λ = 354.84 nm and divergence 2α = 0.3 mrad. The method itself can be used for other beam parameters, within the conditions of applicability of the approximation used. We will consider the Cartesian coordinates (x, y, z) and the spherical coordinates , interconnected as follows: (1) Let us consider a wave field in the space with the imposed radiation condition specifying waves propagating in the direction of the x axis. Let the boundary condition be given for the wave field in the plane x = 0. We will use the Fourier transform in the following normalization: (2) where is the wave number. Then the wave field in vacuum is written as follows: where is a phase function with the following form: The maximum absolute values of , , , and ϕ are estimated by the half-beam divergence α = 0.15 × 10 −3 . The propagation distance r is estimated at 30 km. In the future, we will use the paraxial expansion of the phase function; therefore, we estimate the terms of which degrees should be held: (5) Thus, in expansions, one should hold terms up to and including the 3rd order. However, in fact, terms of the third order are absent in the expansions, so it suffices to take into account the 2nd-order terms. This allows us to write the following approximate expression for the phase function: cos cos , cos sin , sin .
x r y r z r cos cos cos sin sin .  (6) Using this approximation, we write the field in space as a Fourier transform: (7) To numerically find the field in an inhomogeneous medium, we will use the phase-screen method. In the framework of this method, we will consider each screen a fragment of a sphere with some (current) radius r. We will calculate the field on the next screen with a radius r + ∆r using the vacuum propagator and multiplying the resulting field by a factor approximately taking into account the medium. It is easy for a vacuum propagator to derive an explicit expression. To do this, using the field on the initial screen, we will calculate the effective spatial spectrum : The effective spectrum is equal to the spatial spectrum of the field in the x = 0 plane, at which the field would be formed on this phase screen in the absence of a medium. Further, using formula (7) with the substitution r + ∆r and and taking into account the influence of the medium, we derive the field in the following screen: (9) where Ψ(r, Δr, θ, ϕ) is the addition of the optical path along direct rays in the layer between r and r + Δr due to the medium.
We rewrite formula (8) in the following form: , , ; ,  Then the field on the next phase screen can be written similarly: Substituting formula (10) here, we get (12) Here we select the integration over the intermediate variables : Performing integration over η and ζ, we infer (14) Calculating this integral for a spherical wave we derive (15) Using the expression ikr r r r ikr r ikr r r r r r ikr d d formula (14) can be rewritten as follows: where is a vacuum propagator in the paraxial approximation: This expression corresponds to the standard splitstep method. The screen-to-screen conversion is done in three steps. At the first step, the spatial spectrum of the field in the initial phase screen is calculated. In the second step, the spectrum is multiplied by the vacuum propagator and, using the inverse Fourier transform, the field on the next phase screen is calculated without taking into account the medium. In the third step, this field is multiplied by a phase factor that takes into account the phase incursion in the medium between the screens.
The factor can be omitted because it gives a constant phase addition in each screen. In this case, the vacuum propagator can be written as follows: We rewrite the formula (17) in the operator form: where is the operator describing the propagation of the field from screen to screen. Consider 3 consecutive phase screens r 1 , r 2 and r 3 . Consider the composition of the operators where and Moreover, we assume that the medium is a vacuum; i.e.,    where is the distance between the first and third screens. Given this identity, from (25) we obtain (27) This proves the natural group property of vacuum operators : It follows from the accuracy of the group property that this propagator is an exact solution of the approximate form of the parabolic equation, in which y/r and z/r are used as transverse coordinates, and the regular decrease in the amplitude of the spherical wave as 1/r is also taken into account. An alternative approach to accounting for wave divergence is possible using lens coordinates [13][14][15]. Nevertheless, in the literature cited, lens coordinates are used to describe beams in a homogeneous nonlinear dispersive medium. The possibility of applying this transformation to an arbitrary inhomogeneous medium requires additional research.
The spherical phase-screen method can also be used to study focused beams, including self-focusing conditions. Propagator (18) for ∆r < 0 describes not a diverging, but a converging beam.

PHASE INCURSION IN THE THIN SPHERICAL LAYER
Each layer of a randomly inhomogeneous medium of thickness Δr between the phase screens r and r + Δr is described by the function which is understood as the realization of a 2D random field as a function of for given r and Δr. Consider the corresponding implementation of the 3D field of the refractive index . Then we can write We will use the paraxial approximation to calculate this function. Since = 0.675 mm, second-order values can be neglected. Thus, the following approximate expression can be written: We will neglect the regular part of the refractive index, considering it a constant, giving in each screen some known insignificant constant phase shift. Field will be considered a 3D statistically homogeneous and isotropic random field with a spectral density where is a 3D vector of spatial frequencies and . For the correlation function, the standard relation is    and , either not in the form of differences, or with the factor ∆r, both in the numerator and in the denominator, can be approximately excluded. In this case, ∆r is considered small in comparison with r. Then expression (34) can be rewritten in the following approximate form: Using the function, this expression can be rewritten in the following form: This allows us to transform the integral (32) as follows: We obtain an expression for the spectral density of the phase path , where is the spatial frequency vector. The following relation is (38) If we put and then we come to the following expression for the spectral density of the phase path: This formula is written in the approximation of a thin layer of a medium; therefore, it does not take into account the fact that, due to beam broadening, not only transverse frequencies µ/r contribute to the angular frequencies µ. For a layer of a medium with a thickness on the order of the external scale of inhomogeneities and more, it can be written by analogy with the corresponding formula for a plane layer [1,2]: ( )

Numerical Modeling Method
For turbulent fluctuations of the refractive index, we will take the Kolmogorov-Karman spectrum as a model: where is the structural constant, L 0 is the external scale, and λ 0 is the internal turbulence scale.
The numerical algorithm is similar to the case of plane phase screens. Since the vacuum propagator (19) is written as a multiplier in the Fourier-transformed space, its numerical implementation reduces to a fast Fourier transform, multiplying by the propagator, and the inverse Fourier transform.
To generate realizations of random uncorrelated phase screens, we use the discrete form of the relation (38) [16]: where and is the discrete Fourier transform from : Thus, are uncorrelated random variables with random phases and RMS data:  The realization of the optical screen thickness is equal to the inverse discrete Fourier transform from [5]: To take into account the fact that the same sampling interval over angular coordinates at different r corresponds to different spatial scales, we use adaptive sampling. In each phase screen, the required sampling step is estimated in angular coordinates; if the current sampling step exceeds this estimate, then the resolution is doubled. For this purpose, the Kotelnikov interpolation is used: in the Fourier space, the frequency grid is supplemented so that, at the same frequency resolution, the maximum frequency becomes twice as high. In the added grid nodes, the Fourier image of the field is set to 0. After the inverse Fourier transform, we get a field interpolated to a spatial grid with half the initial sampling step.
Consider the issue of numerical efficiency of the method. For both plane and spherical screens, the sampling interval should not be below the characteristic scale of the medium (of the order of or greater than the internal scale for turbulence). However, for plane phase screens it is also necessary to take into account that, at the edges, there is a regular phase incursion on the order of which in our example is about 5000 radians. The necessary sampling interval is where is the allowable phase variation between discrete samples of the field.
If we take then mm. With the divergence indicated in the article, the characteristic beam size reaches 5 m, which will require a grid dimension of the order of points, which is a very large value that requires not only a lot of time, but RAM, which is not always available on modern r r r processors. In fact, for the calculation, a step of about 3 mm is sufficient, which reduces the grid dimension by more than 100 times, and the calculation speed of the fast Fourier transform, proportional to increases even more.
The size of the phase screen is usually close to the size of the beam cross section with some margin to account for possible beam broadening. However, this log , N N approach to modeling does not take into account the influence of the low-frequency part of the spectrum of random inhomogeneities, in particular, the external scale, if it significantly exceeds the beam size. Figures 1-6 show examples of a numerical calculation of the amplitude distributions of a Gaussian laser      beam with a wavelength of λ = 354.84 nm with a divergence of 2α = 0.3 × 10 -3 rad at different points of the 24 km long path for = 2.5 × 10 -15 m -2/3 . The figures show the formation of speckles. The step between the phase screens was 100 m.

CONCLUSIONS
In this paper we have developed a new modification of the phase-screen method which reduces com-2 n C putational costs when calculating diverging wave beams. This is achieved through the use of spherical phase screens that repeat the shape of unperturbed wave fronts. In the paraxial approximation, we obtained an expression for the vacuum propagator describing screen-to-screen propagation. We also obtained a formula relating the spectral density of isotropic fluctuations of the refractive index with a phase incursion in a spherical layer of a random medium. This formula is similar to the plane-layer formula. The numerical implementation of the method of spherical phase screens is close to the numerical implementation of the method of plane screens. The paper presents numerical examples of the calculation of fluctuations of a laser beam on a 24 km long path in the regimes of weak and strong fluctuations.