Method of spherical phase screens for the modeling of propagation of diverging beams in inhomogeneous media

The phase-screen (split-step) method is widely used for the modeling of wave propagation in inhomogeneous media. Most known is the method of flat phase screens. An optimized approach based on cylindrical phase screen was introduced for the 2-D modeling of radio occultation sounding of the Earth’s atmosphere. In this paper, we propose a further generalization of this method for the 3-D problem of propagation of diverging beams. Our generalization is based on spherical phase screens. In the paraxial approximation, we derive the formula for the vacuum screento-screen propagator. We also derive the expression for the phase thickness of a thin layer of an isotropic random media. We describe a numerical implementation of this method and give numerical examples of its application for the modeling of a diverging laser beam propagating on a 25 km long atmospheric path.


Introduction
The method of phase screens has been widely used for the numerical simulation of the wave propagation of various nature in inhomogeneous media, including the modeling of the optical (laser) radiation propagation in a turbulent atmosphere [1][2][3][4][5] and the decimeter waves propagation during radio occultation sounding of the atmosphere [6][7][8]. This method is referred to as split-step. 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 the composition of an infinitely thin layer that forms phase distortions of the wave and a vacuum propagator describing diffraction.
The phase screen method has a fundamental limitation: it does not account for backscattering. The method of phase screens can also be considered as a finite-dimensional approximation of the path integral [9]. In problems of modeling of laser radiation in turbulent media, especially in order to describe the effect of isotropic turbulence, 2-dimensional phase screens are used. In modeling of radio occultation experiments, however, significant optimization of computational costs is achieved by employing the approximation of onedimensional phase screens, because in most cases atmospheric inhomogeneities with vertical scales from hundreds of meters to kilometers are taken into account, while their horizontal scales significantly exceed the horizontal size of the Fresnel zone.
In the classical version of the method, flat phase screens are used. This leads to excessive computational costs when describing a diverging wave: the increasing angle between the screen and the wavefront at the edges of each screen results in oversampling. In the two-dimensional (2D) modeling of radio occultation experiments, it turned out to be quite simple to write down a solution for cylindrical 1-dimensional phase screens [7], which takes 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, using the paraxial approximation. In the Section 2, we derive the basic relations using the technique of angular spectra [10]. The main result here is the formula for the propagator in spherical coordinates in the small angle approximation. In Section 3, we show that the numerical implementation of this method is no more complicated than the case of flat phase screens, and we give examples of numerical modeling using the model of isotropic turbulence. In Section 4, we offer our conclusions.

Vacuum propagator
As an example, we will use the beam parameters specified for the DELICAT project (DEmonstration of LIdar based Clear Air Turbulence detection, Demonstration of clear sky turbulence detection using lidar) [11,12]: wavelength 354.84 nm, divergence 20.3 mrad. The described method can also be used for other beam parameters, within the framework of the applicability conditions of the approximation used.
We will consider the Cartesian coordinates   Let us consider a wave field in space   ,, u x y z with an imposed radiation condition that implies waves propagating in the direction of the axis x . We denote the boundary condition for the field in the plane 0 x  as   0 , u y z . Hereinafter, we use the Fourier transform in the following normalization: where 2/ k    is the wave number. Then the wave field in vacuum is written as follows: The maximum absolute values  ,  ,  , and  are estimated by the value of the half beam divergence of 0.15×10 −3 . The propagation distance r is estimated at 30 km.
We can write the following approximate expression for the phase function: Our numerical solution for the field in an inhomogeneous medium will be based on the split-step method. In the framework of this method, the medium is split in spherical layers, represented by phase screens centered at the wave source. We will calculate the field on the next screen using the vacuum propagator and multiplying the resulting field by the factor that describes the influence of the random medium: where   , , , rr     is the phase thickness of the phase screen, and   , ; , P r r    is the vacuum propagator. In the paraxial approximation, the following expression can be derived: This formula follows the standard split-step method of phase screens. Propagation of the field from screen to screen is performed in the following steps: 1) The spatial spectrum of the field in the initial phase screen is calculated. 2) The spectrum is multiplied by the vacuum propagator. 3) The inverse Fourier transform is taken producing the field on the next phase screen, without taking into account the medium. 4) The field is multiplied by the phase factor that takes into account the phase perturbation in the medium between the screens. The We rewrite the formula (6) in the operator form: where   , P r r  is the operator describing the screen-to-screen propagation of the field and possessing the natural group property: 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 / yr and / zr are used as transverse coordinates, and the regular decrease in the amplitude of the spherical wave as 1/ r is taken into account. An alternative approach to accounting for wave divergence is possible using lens coordinates [13][14][15]. Nevertheless, in the referenced papers, the lens coordinates are used to describe beams in a homogeneous nonlinear dispersive medium. The possibility of applying this transformation to an arbitrary heterogeneous medium requires additional research.
The method of spherical phase screens can also be used to study focused and selffocused beams. Propagator (7) at describes not a diverging, but a converging beam.

Phase raid in a thin spherical layer
Each layer of a random medium thick r  between phase screens r and rr  is described by the function   , , , rr     , which is the realization of a 2D random field as a function of ,  for given r and r  . Consider a realization of the 3-D field of the refractive index We will neglect the regular part of the refractive index, considering it a constant, producing in each screen a constant phase shift. The field   ,, N x y z will be considered a 3-D statistically homogeneous and isotropic random field with spectral density is the three-dimensional vector of spatial frequencies, and  κ . For the correlation function, the standard relation holds: This formula is written in the approximation of a thin layer of a medium with rr  , neglecting the beam broadening. For a layer of a medium with a thickness of the order of  (15) 3 Numerical modeling

Numerical modeling method
For turbulent fluctuations of the refractive index, we assume the Kolmogorov-Karman spectrum: where The numerical algorithm is similar to the case of flat phase screens. Since the vacuum propagator (8) is written as a multiplier in the Fourier-transformed space, it is numerically implemented using the fast Fourier transform.
To generate realizations of random uncorrelated phase screens, we use the discrete form of relation the correlation of the phase path [16]: To take into account the fact that the same discretization step according to angular coordinates for different radiuses corresponds to different spatial scales, we use adaptive discretization. On each phase screen, the required sampling step is estimated in angular coordinates, and if the current sampling step exceeds this estimate, the resolution is doubled. To this end, the Kotelnikov interpolation is used: in the space of Fourier transforms. The frequency grid is enlarged, keeping the frequency resolution, and increasing the maximum frequency twice. 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 the spatial grid with a half-step discretization.

Conclusions
In this work, we developed a new modification of the phase-screen method, which provides a reduction in computational costs when modeling the propagation of diverging wave beams. The method utilizes spherical phase screens that better conform with the shape of unperturbed wave fronts. In the paraxial approximation, we obtained the expression for the vacuum propagator describing the propagation from screen to screen. We also obtained a formula relating the spectral density of isotropic fluctuations of the phase thickness of a spherical phase screen in an isotropic random medium. The formula is similar to the flat layer formula. The numerical implementation of the method of spherical phase screens is very close to the numerical implementation of the method of flat screens. Numerical examples of the calculation of fluctuations of a laser beam propagating along a path 25 km long in the weak and strong fluctuation modes are presented.
This work was carried out with financial support of the Russian Foundation for Basic Research (grant No. 18-35-00368).