Method of spherical phase screens for the

modeling of propagation of diverging beams

in inhomogeneous media

Michael E. Gorbunov

1,2

, Oksana A. Koval

1,

*

, Victor A. Kulikov

1

,

and Alexey E. Mamontov

1

1

A.M. Obukhov Institute of Atmospheric Physics, Russian Academy of Sciences, Pyzhevsky per., 3,

Moscow 119017, Russia

2

Hydrometeorological Research Center of Russian Federation, Bolshoy Predtechensky per., 11-13,

Moscow 123242, Russia

Abstract. 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 screen-

to-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.

1 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–5] and the decimeter

waves propagation during radio occultation sounding of the atmosphere [6–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

*

Corresponding author: kov.oksana20@gmail.com

,

(2019)

https://doi.org/10.1051/itmconf /201930

ITM Web of Conferences

30

CriMiCo'2019

1

150

50

2

27

7

© The Authors, published by EDP Sciences. This is an open access article distributed under the terms of

the

Creative

Commons

Attribution

License 4.0 (http://creativecommons.org/licenses/by/4.0/).

screens are used. In modeling of radio occultation experiments, however, significant

optimization of computational costs is achieved by employing the approximation of one-

dimensional 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.

2 Conclusion of basic relations

2.1 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

,,x y z

and the spherical coordinates

, , ,r

interconnected as follows:

cos cos ,

cos sin ,

sin .

xr

yr

zr

(1)

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

0x

as

0

,u y z

. Hereinafter, we use the Fourier

transform in the following normalization:

, , exp ,

2

, , exp ,

2

k

f f y z ik y ik z dy dz

k

f y z f ik y ik z d d

(2)

where

2/k

is the wave number. Then the wave field in vacuum is written as follows:

0

, , , exp , , ; , ,

2

k

u r u ikS r d d

(3)

,

(2019)

https://doi.org/10.1051/itmconf /201930

ITM Web of Conferences

30

CriMiCo'2019

1

150

50

2

27

7

2

where

, , ; ,Sr

is phase function, having the following form:

22

, , ; , 1 cos cos cos sin sin .S r r

(4)

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:

2 2 2 2

, , ; , 1 .

22

S r r

(5)

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:

, , exp , , ,

, ; , exp

2

, , exp ,

2

u r r ik r r

kr

d d P r r ikr

kr

u r ikr d d

(6)

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:

22

, ; , exp 1 .

2

rr

P r r ik r

r r r r

(7)

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 multiplier

exp ik r

can be omitted because it provides a constant phase addition

in each screen. In this case, the vacuum propagator can be written as follows:

22

, ; , exp .

2

r ikr r

P r r

r r r r

(8)

We rewrite the formula (6) in the operator form:

ˆ

, , , , , ,u r r P r r u r

(9)

where

ˆ

,P r r

is the operator describing the screen-to-screen propagation of the field and

possessing the natural group property:

,

(2019)

https://doi.org/10.1051/itmconf /201930

ITM Web of Conferences

30

CriMiCo'2019

1

150

50

2

27

7

3

2 23 1 12 1 13

ˆ ˆ ˆ

, , , .P r r P r r P r r

(10)

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–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 self-

focused beams. Propagator (7) at describes not a diverging, but a converging beam.

2.2 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

,,N x y z

. Then

, , ,rr

is expressed as follows:

, , , cos cos , cos sin , sin .

rr

r

r r N r r r dr

(11)

To calculate this function, we will use the paraxial approximation. Since

2

r

= 0.675

mm, the values of the second order can be neglected. Thus, we arrive at the following

approximate expression:

, , , , , .

rr

r

r r N r r r dr

(12)

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

NN

κ

, where

,,

x y z

κ

is the three-dimensional vector of spatial

frequencies, and

κ

. For the correlation function, the standard relation holds:

*3

1 2 1 2

exp .

N

N N i d

rr κ r r κ κ

(13)

Expressing explicitly the correlation of the phase path

, , , , , ,r r r r

and inserting

y

r

and

z

r

, we arrive at the expression for the spectral density

of the phase path:

2

2

2

sinc , , .

2

x

N x x

r

r

d

rr

r

μ

(14)

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

,

(2019)

https://doi.org/10.1051/itmconf /201930

ITM Web of Conferences

30

CriMiCo'2019

1

150

50

2

27

7

4

the external scale, the formula can be approximately rewritten like the corresponding

formula for a flat layer [1,2]:

2

2

0, , .

N

r

rr

r

μ

(15)

3 Numerical modeling

3.1 Numerical modeling method

For turbulent fluctuations of the refractive index, we assume the Kolmogorov–Karman

spectrum:

11/6

2 2 2 2 2

0

( ) 0.033 exp / ,

N n m

C

(16)

where

2

n

C

is the structural constant,

00

2/L

,

0

L

is the external scale,

0

5.92 /

m

,

and

0

is the internal turbulence scale.

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]:

2

2

2

22

22

exp ,

44

j l jl

jl

jl

jj ll

ii

NN

N N N N

(17)

where

,

j l j l

and

jl

are discrete Fourier transform of

jl

:

22

exp .

jl j l

jl

jj ll

ii

NN

(18)

Thus,

jl

are uncorrelated random variables with random phases and rms values:

2

2

4

.

jl

jl

NN

(19)

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.

,

(2019)

https://doi.org/10.1051/itmconf /201930

ITM Web of Conferences

30

CriMiCo'2019

1

150

50

2

27

7

5