Feynman Path Integrals - Application to Wave Propagation

In this section we introduce an alternative representation for a Green function based on Feynman path integrals. This approach is most suitable for the analysis of the wave propagation and scattering in a line-of-sight region. This section provides some basics for study in an unbounded medium, while later in Wave Field Fluctuations in Random Media over a Boundary Interface“Understanding Wave Field Fluctuations in Random Media” we introduce a path integral representation in the presence of a boundary.

The basic Feynman postulate states that the probability amplitude for the transition of a particle-wave from point

{\vec{r}}_{0}

to point

\vec{r}

can be represented by contributions made by individual trajectories over which the particle-wave can propagate between points

{\vec{r}}_{0}

and

\vec{r}

. The contribution of each trajectory is proportional to exp {jS/ħ}, where S is a classical action and ħ is Planck’s constant. As is well known, the quantum-mechanical probability amplitude of transition is equivalent to the Green function

G (\vec{r}, {\vec{r}}_{0})

of the Helmholtz equation written in a parabolic approximation, after the time has been replaced by the x-coordinate along the selected direction of propagation and m/ħ has been replaced by the wavenumber k = 2π/λ, where m is a particle’s mass, k is the wave number of the admitted wave, and λ is the wavelength. Turning to the trajectory continuum, we can write:

G (\vec{r}, {\vec{r}}_{0}) = \int D \vec{ρ} (x) e x p \{j k S [\vec{ρ} (x)]\}, E q . 1

where:

$\vec{r} = \{x, \vec{ρ}\}, {\vec{r}}_{0} = \{0, {\vec{ρ}}_{0}\}$
are the coordinates of the observation point and source, respectively;
$\vec{ρ} = \{γ, z\}, \vec{ρ} = \{γ_{0}, z_{0}\}$
are the vectors in the plane orthogonal to the direction of propagation x, and;
$D \vec{ρ} (x)$
is the differential in the space of continuous trajectories.

The action S is given by:

S [\vec{ρ} (x)] = \int_{0}^{x} d x^{'} [L_{0} (x^{'}, \vec{ρ} (x^{'})) + \frac{1}{2} δ ε (x^{'}, \vec{ρ} (x^{'}))], E q . 2

where:

$L_{0} (x, \vec{ρ} (x)) = \frac{1}{2} {(\frac{d ρ}{d x})}^{2} + U_{0} (x, \vec{ρ} (x))$
is an unperturbed (when δε = 0) Lagrangian in the small-angle approximation;
$U_{0} (x, \vec{p} (x))$
is an unperturbed potential;
$U_{0} (x, \vec{p} (x)) = 1 / 2 <ε (x, \vec{ρ} (x)) - 1>$
, the angle brackets < … > here and later denote averaging over the ensemble of dielectric permittivity;
$ε (x, \vec{ρ}), δ ε (x, \vec{ρ})$
is a random component of
$ε (x, \vec{ρ}), <δ ε> = 0 .$

Consider calculation of the first two moments of the Green function (1) in the case of an unbounded stationary medium filled with random inhomogeneities of dielectric permeability

ε (x, \vec{ρ})

. In continuous notation the average Green function

<G (\vec{r}, {\vec{r}}_{0})>

and coherence function:

Γ (x, {\vec{ρ}}_{1}, {\vec{ρ}}_{2}, ρ_{0}^{'}, {\vec{ρ}}_{0}^{''}) = <G ({\vec{r}}_{1}, {\vec{r}}_{0}^{'}) G^{*} ({\vec{r}}_{2}, {\vec{r}}_{0}^{''})>

<G (\vec{r}, {\vec{r}}_{0})> = \int D \vec{ρ} (x) e x p \{j k S_{0} [\vec{ρ} (x)]\} <e x p (j \frac{k}{2} \int_{0}^{x} d x^{'} δ ε (x^{'}, \vec{ρ} (x^{'})))>, E q . 3

Γ (x, {\vec{ρ}}_{1}, {\vec{ρ}}_{2}, {\vec{ρ}}_{0}^{'}, {\vec{ρ}}_{0}^{''}) = <G ({\vec{r}}_{1}, {\vec{r}}_{0}^{'}) G^{*} ({\vec{r}}_{2}, {\vec{r}}_{0}^{''})> = \int D \vec{ρ} (x) \int D \vec{R} (x) e x p \{j k \int_{0}^{x} d x^{'} \frac{d \vec{R}}{d x^{'}} \frac{d \vec{ρ}}{d x^{'}}\} \cdot <e x p \{j \frac{k}{2} \int_{0}^{x} d x^{'} [δ ε (x^{'}, \vec{R} (x^{'}) + \frac{\vec{ρ} (x^{'})}{2}) - δ ε (x^{'}, \vec{R} (x^{'}) - \frac{\vec{ρ} (x^{'})}{2})]\}> . E q . 4

In Eq. (4) we have introduced the sum

\vec{R} (x) = 1 / 2 ({\vec{ρ}}_{1} + {\vec{ρ}}_{2} (x))

and the difference

\vec{ρ} (x) = {\vec{ρ}}_{1} (x) - {\vec{ρ}}_{2} (x)

of the trajectories with the following boundary conditions:

\vec{R} (x^{'} = 0) = \frac{1}{2} ({\vec{ρ}}_{0}^{'} + {\vec{ρ}}_{0}^{''}) \vec{ρ} (x^{'} = 0) = {\vec{ρ}}_{0}^{'} - {\vec{ρ}}_{0}^{''},

\vec{R} (x^{'} = x) = \frac{1}{2} ({\vec{ρ}}_{1} + {\vec{ρ}}_{2}) \vec{ρ} (x^{'} = x) = {\vec{ρ}}_{1} - {\vec{ρ}}_{2} .

We can assume that the fluctuations in δε are statistically uniform:

<δ ε (x_{1}, {\vec{ρ}}_{1}) δ ε (x_{2}, {\vec{ρ}}_{2})> = B_{ε} (x_{1} - x_{2}, {\vec{ρ}}_{1} - {\vec{ρ}}_{2})

and have different correlation scales:

L_x in direction x and;
L_⟂ in the plane (γ, z).

If x ≫ L the integral

\int_{0}^{x} d x^{'} δ ε (x^{'}, \vec{ρ} (x^{'}))

represents a Gaussian random value in accordance with the central-limit theorem. Hence:

<e x p (j \frac{k}{2} \int_{0}^{x} d x^{'} δ ε (x^{'}, \vec{ρ} (x^{'})))> = e x p (- γ_{s}), E q . 5

<e x p \{j \frac{k}{2} \int_{0}^{x} d x^{'} [δ ε (x^{'}, \vec{R} (x^{'}) + \frac{\vec{ρ} (x^{'})}{2}) - δ ε (x^{'}, \vec{R} (x^{'}) - \frac{\vec{ρ} (x^{'})}{2})]\}> = e x p (- D_{s}), E q . 6

where:

γ_{S} = \frac{k^{2}}{8} \int_{0}^{x} d x^{'} \int_{0}^{x} d x^{''} <δ ε (x^{'}, \vec{ρ} (x^{'})) δ ε (x^{''}, \vec{ρ} (x^{''}))> E q . 7

is the variance of the phase fluctuations along the trajectory

\vec{ρ} (x)

and:

D_{S} = \frac{k^{2}}{8} \int_{0}^{x} d x^{'} \int_{0}^{x} d x^{''} <[δ ε (x^{'}, \vec{R} (x^{'}) + \frac{\vec{ρ} (x^{'})}{2}) - δ ε (x^{'}, \vec{R} (x^{'}) - \frac{\vec{ρ} (x^{'})}{2})] \cdot δ ε (x^{''}, \vec{R} (x^{''}) + \frac{\vec{ρ} (x^{''})}{2}) - δ ε (x^{''}, \vec{R} (x^{''}) - \frac{\vec{ρ} (x^{''})}{2})> E q . 8

is a structure function of the phase difference along the trajectories

\vec{R} (x) + \frac{\vec{ρ} (x)}{2}

and

\vec{R} (x) - \frac{\vec{ρ} (x)}{2} .

Introduce a two-dimensional spatial spectrum of the fluctuations de in a plane (γ, z):

F_{ε} (x^{'} - x^{''}, {\vec{κ}}_{⊥}) = \frac{1}{4 π^{2}} \int d^{2} \vec{ρ} B_{ε} (x^{'} - x^{''}, \vec{ρ}) e^{- j \vec{κ} \vec{ρ}} . E q . 9

Then for γ_S and D_s we obtain:

γ_{S} = \frac{k^{2}}{4} \int_{0}^{x / 2} d η \int_{- 2 η}^{2 η} d ξ \int d^{2} {\vec{κ}}_{⊥} F_{ε} (ξ, {\vec{κ}}_{⊥}) e x p (j \vec{κ} ξ \frac{d \vec{ρ}}{d η}), E q . 10

D_{S} [\vec{R}, \vec{ρ}] = \frac{k^{2}}{4} \int_{0}^{x / 2} d η \int_{- 2 η}^{2 η} d ξ \int d ξ \int d^{2} {\vec{κ}}_{⊥} F_{ε} (ξ, {\vec{κ}}_{⊥}) e x p (j \vec{κ} ξ \frac{d \vec{R}}{d η}) \cdot [\cos (\vec{κ} ξ \frac{d \vec{ρ}}{d η}) - \cos (\vec{κ}, \vec{ρ})] . E q . 11

To analyse Eqs. (10) and (11) let us introduce the parameter

θ_{R} = m a x [\frac{d \vec{R}}{d x}],

the characteristic angle of the ray’s trajectory along the x-axis. The value of θ_R is contributed to by the three characteristic parameters:

θ_d, the gradient of the ray trajectory due to either the position of the correspondent or regular refraction (if any), in case of normal refraction θ_d is given by Understanding Parabolic Approximation in Wave Propagation: Analytical Methods and Applicationsthis equation);
$θ_{F} = 1 / \sqrt{k x},$
the angular size of the Fresnel zone;
θ_s, the characteristic angle of scattering on the fluctuations δε, θ_s ~ 1/kL_⟂.

Therefore θ_R = max{θ_d, θ_F, θ_s}.

Let us also introduce α parameter of anisotropy α in the fluctuations of δε: α = L_⟂ / L_x. The effective width of function

F_{ε} (ξ, {\vec{κ}}_{⊥})

over ξ does not exceed a correlation scale L_x. When inequality:

θ_{R} ≪ α E q . 12

holds, we can assume that

{\vec{κ}}_{⊥} ξ d \vec{ρ} / d η, {\vec{κ}}_{⊥} ξ d \vec{R} / d η ≪ d η 1

in a significant region over ξ in Eqs. (10) and (11). As a result, we obtain the following approximations:

γ_{S} (x) = \frac{π k^{2} x}{4} \int d^{2} {\vec{κ}}_{⊥} Φ_{ε} (0, κ_{⊥}); E q . 13

and:

D_{S} [\vec{R}, \vec{ρ}] \equiv D_{S} |\vec{ρ}| = \frac{π k^{2}}{2} \int_{0}^{x} d x^{'} \int d^{2} κ_{⊥} Φ_{ε} (0, κ_{⊥}) [1 - \cos (κ_{⊥} \overset{⇀}{ρ} (x^{'}))] E q . 14

where:

2 π Φ_{ε} (0, {\vec{κ}}_{⊥}) = \int_{- \infty}^{\infty} d ξ F_{ε} (ξ, {\vec{κ}}_{⊥}) E q . 15

and

Φ_{ε} (\vec{κ})

is a three-dimensional spatial spectrum of fluctuations in δε. We can reasonably assume that θ_R ~ θ_s, and the inequality (12) takes the form:

\frac{k L_{⊥}^{2}}{L_{x}} ≫ 1 E q . 16

which means that the longitudinal correlation scale has to be small compared with the distance

k L_{⊥}^{2}

where the diffraction on the irregularities of the scale L_⟂ become significant. The expressions (5) and (14) provide the fundamental solution to the equation of the coherence function Γ in the Markov approximation:

\frac{\partial Γ}{\partial x} - \frac{j}{k} \frac{\partial^{2} Γ}{\partial \vec{R} \partial \vec{ρ}} - \frac{π k^{2}}{4} H (\vec{ρ}) Γ = 0, E q . 17

H (\vec{ρ}) = 2 \int d^{2} κ_{⊥} Φ_{ε} (0, κ_{⊥}) [1 - \cos (κ_{⊥} \overset{⇀}{ρ} (x^{'}))]

and:

D_{S} [\vec{ρ}] = π k^{2} / 4 \int_{0}^{x} d x^{'} H (\vec{ρ} (x^{'})) .

“Local” conditions of the Markov approximation are bounded by inequalities (3) and (4), for an average field we also need smallness of the attenuation of the average field over the distance of L_x:γ_S(L_x) ≪ 1. “Non-local“ conditions lead to the inequality:

R_{c} ≫ λ, E q . 18

where R_c is a coherence radius of the scattered field in a plane (γ, z), which depends on the distance x, R_c ≡ R_c(x). The parameter R_c can be found from:

D_{S} (R_{c}) . E q . 19

To illustrate the path integral approach, we obtain a known solution to Eq. (17) for the unbounded medium filled with statistically uniform fluctuations of

δ ε (x, \vec{ρ}) .

Using Eqs. (4), (16) and (14) we obtain:

Γ (x, {\vec{ρ}}_{1}, {\vec{ρ}}_{2}, {\vec{ρ}}_{0}^{'}, {\vec{ρ}}_{0}^{''}) = \int D \vec{ρ} (x) \int D \vec{R} (x) \cdot e x p \{j k \int_{0}^{x} d x^{'} \frac{d \vec{R}}{d x^{'}} \frac{d \vec{ρ}}{d x^{'}} - D_{S} [\vec{ρ} (x^{'})]\} . E q . 20

The integral over trajectories R(x^′) represents a continuous Fourier transform of the delta-function

δ (d^{2} \vec{ρ} / d x^{2}) .

A Lagrangian in the exponent can be expanded into a functional series up to the second order terms in the vicinity of the trajectory

{\vec{ρ}}_{0} (x^{'}) = \vec{ρ} x^{'} / x + {\vec{ρ}}_{0} (1 - x^{'} / x) .

Integrating then over

\vec{ρ} (x)

, we obtain the well-known solution for an unbounded medium:

Γ (x, {\vec{ρ}}_{1}, {\vec{ρ}}_{1}, {\vec{ρ}}_{0}^{'}, {\vec{ρ}}_{0}^{''}) = \frac{k^{2}}{4 π^{2} x^{2}} e x p \{j k \frac{(\vec{R} - {\vec{R}}_{0}) (\vec{ρ} - {\vec{ρ}}_{0})}{x} - \frac{π k^{2}}{4} \int_{0}^{x} d x^{'} H (\vec{ρ} \frac{x^{'}}{x} + {\vec{ρ}}_{0} (1 - \frac{x^{'}}{x}))\} . E q . 21

The relationship between the average of the path integrals (3), (4) and the average over the Fermat paths was analysed in Refs. The asymptotic evaluation of the path integral representation for the Green function (3) may be performed using the method of stationary phase. The asymptotic method consists of finding the extreme path

{\vec{ρ}}^{*} (x)

which renders the minimum value of the phase

S [\vec{ρ} (x)]

given by Eq. (2). It can be noted that in obtaining the solution for the coherence function in an unbounded medium (21) we actually did not face any difficulties with integration over R(x) due to an isotropy of the medium,

<ε (\vec{r})> = c o n s t .

However, this is not the case for more complex media. The variational problem for determining

{\vec{ρ}}^{*} (x)

gives rise to the Euler equation for the Fermat paths:

\frac{d^{2} \vec{ρ}}{d x^{2}} + \nabla δ ε (x, \vec{ρ} (x)) = 0, E q . 22

where

\nabla = \{\frac{\partial}{\partial γ}, \frac{\partial}{\partial z}\} .

The rigorous mathematical solution to the problem encountered difficulties in a general case of the multi-scale inhomogeneities of δε and, therefore, the basic physics has to be involved for a qualitative analysis.

As mentioned above, there are two characteristic transverse scales of the phase fluctuations in

S [\vec{ρ}] .

The first is the Fresnel zone size resulting from the first term in Eq. (3). The second term is due to a random component of the dielectric permittivity de and is of the order of the coherence scale R_c introduced by Eq. (19).

We assume a small value of the mean-square variation of the fluctuations in a phase difference of the fields approaching any points separated in space by a distance of the order of the Fresnel zone size. Hence, this assumption is rendered by the inequality:

δ S_{1} = \frac{k^{2} π}{4} \int_{0}^{x} d x^{'} H (\sqrt{\frac{x^{'}}{κ}}) ≪ 1 . E q . 23

As known, the meaning of this inequality is that the Fresnel zone size is the characteristic region in the plane x^′ = constant, from which the rays arrive at the receiving point in phase, even in the presence of the δε fluctuations. It also means that the integration in Eq. (3) can be fulfilled along a single canal-ray tube bounded in lateral cross-section by the Fresnel volume. We assume the extreme trajectory

{\vec{ρ}}^{*} (x)

to deviate slightly from the unperturbed path

{\vec{ρ}}_{0} (x)

, defined by the unperturbed Lagrangian

L_{0} [\vec{ρ} (x)],

we should require the fluctuations in the arrival angle of the wave defined by

{\vec{ρ}}^{*} (x)

to be small compared with angular size of the Fresnel zone ~(kx)^-1/2. In a functional representation this requirement is equivalent to:

δ S_{2} = \frac{k^{2} π}{8} \int_{0}^{x} d x^{'} δ {\vec{ρ}}^{2} (x^{'}) \nabla^{2} H (\vec{ρ} (x^{'})) \leq \frac{k^{2} π}{8} x^{2} \int d^{2} {\vec{κ}}_{⊥} ≪ 1 . E q . 24

In Eq. (7) we evaluated the deviation

δ \vec{ρ} = {\vec{ρ}}^{*} - {\vec{ρ}}_{0}

to be of the order of the Fresnel zone size.

Having satisfied inequalities (22) and (23), the integration in Eq. (3) can be performed along the unperturbed ray trajectories

{\vec{ρ}}_{0} (x)

which are the Fermat paths for a free wave particle propagation (in the particular case of an unbounded medium):

{\vec{ρ}}_{0} (x^{'}) = (\vec{ρ} - {\vec{ρ}}_{0}) \frac{x^{'}}{x} . E q . 25

Such a straight-line approximation of trajectories in a path integral is similar to the extended Huygens–Fresnel principle. The applicability of extended Huygens–Fresnel principle shown that this approximation besides yielding the exact solution for the first two moments, provides a qualitatively correct solution for the high-order moments of the wave field.

Numerical Methods of Parabolic Equations

Among many numerical methods used in the problems of applied electromagnetics and wave propagation one may distinguish two methods most widely used in a VHF/UHF propagation in the atmospheric boundary layer, namely: split-step-Fourier and split-step Pade methods. Both methods are based on the parabolic approximation to a wave equation and differ in the method of obtaining the approximation of the exponential operator. Initially, the split-step parabolic approximation was introduced in a problem of applied acoustics where the wave propagation can be described by a scalar wave equation of elliptical type:

\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial z^{2}} + k^{2} n^{2} u = 0 E q . 26

where u is the spectral amplitude of the wave component with frequency ω, k = ω/c is the wavenumber, c is the group velocity of the propagation and n is the refractive index of the medium and may vary with the coordinates, we assume, for now, a constant value of n, n = 1. Equation (26) is written for a two-dimensional case of the propagation since this case is the one most widely considered in applications of the split-step approximation.

Equation (26) can be factorised as follows:

[\frac{\partial}{\partial x} + j \sqrt{k^{2} + \frac{\partial^{2}}{\partial z^{2}}}] \cdot [\frac{\partial}{\partial x} - j \sqrt{k^{2} + \frac{\partial^{2}}{\partial z^{2}}}] u = 0, E q . 27

representing the forward and backward propagating waves respectively. We may also remove fast phase variations with distance by introducing the slow-varying amplitude ũ, u = ũ exp (jkx). The truncated and factorised equation for amplitude ũ takes the form:

[\frac{\partial}{\partial x} + j k (1 + \sqrt{1 + \frac{1}{k^{2}} \frac{\partial^{2}}{\partial z^{2}}})] \cdot [\frac{\partial}{\partial x} + j k (1 - \sqrt{1 + \frac{1}{k^{2}} \frac{\partial^{2}}{\partial z^{2}}})] ũ = 0 . E q . 28

By retaining the only forward propagating wave we obtain the forward parabolic equation:

[\frac{\partial}{\partial x} + j k (1 - \sqrt{1 + \frac{1}{k^{2}} \frac{\partial^{2}}{\partial z^{2}}})] ũ = 0 . E q . 29

which has a formal solution:

ũ (x + ∆ x) = e x p \{j k ∆ x [\sqrt{1 + \frac{1}{k^{2}} \frac{\partial^{2}}{\partial z^{2}} - 1}]\} ũ (x) E q . 30

where Δx is a range increment. As apparent from Eq. (30) the forward propagating wave at range x + Δx can be obtained from values of the wave amplitude at the previous distance x by applying the exponential operator in Eq. (30). In order to actually calculate the wave field, the exponential operator in Eq. (30) should be approximated in a form suitable for computations. Let us introduce a notation Z = 1/k²∂²/∂z² and, depending on the type of approximation, we obtain three known types of the split-step approximation to parabolic equations:

a The standard or narrow angle approximation, obtained by expanding the square root in the exponent into a Taylor series and retaining the linear term:

\sqrt{1 + Z} = 1 + \frac{Z}{2} . E q . 31

This approximation is used in a split-step Fourier method to be discussed later, and is normally valid for narrow angles of propagation relative to axis x, not exceeding 15°–20°.

b The Claerbout approximation in the form:

\sqrt{1 + Z} \approx \frac{1 + \frac{3}{4} Z}{1 + \frac{Z}{4}} E q . 32

which is reported to be valid for wider angles up to 30–40°. And finally:

c The split-step Pade approximation, when exponential is expanded into series:

e x p [j k ∆ \sqrt{1 + Z} - 1] \approx 1 + \sum_{m = 1}^{M} \frac{α_{m} Z}{1 + b_{m} Z} E q . 33

where the coefficients α_m, b_m are determined numerically in the complex plane. The split-step Pade approximation is reported to be valid for angles of propagation relative to the x-axis of up to 90°.

All the above methods have been realized in a very powerful computational tech- nique well suited for parallel computing and widely used in numerous applications. With regard to the problem of radio wave propagation in the earth’s troposphere the most established computational approach is based on the split-step Fourier method equivalent to a standard narrow angle approximation of the square root in the expo- nential operator. Whilst the wide angle approximations, such as Claerbout and Pade approximations have been used in the problems of scattering on objects submerged in a free-space, the author does not know of any systematic derivation of the PadØ approximation in the case of radio wave propagation in a stratified troposphere over terrain or the sea boundary surface. We may notice that while the fundamental elliptic equation of type (22) can be obtained for potentials (Hertz functions) under free-space propagation conditions, the equations for the Debye’s potentials take a different form in the case of a stratified troposphere. In the presence of a random component of the refractive index the main equation for the slow varying envelope is given by Understanding Parabolic Approximation in Wave Propagation: Analytical Methods and Applicationsthis equation.

Consider Understanding Parabolic Approximation in Wave Propagation: Analytical Methods and Applicationsthis equation and retain the term

\frac{1}{r^{2}} \frac{\partial^{2} W_{1}}{\partial θ^{2}}

on the right-hand side of the equation. While this is not very consistent with the arguments used in Understanding Parabolic Approximation in Wave Propagation: Analytical Methods and Applications“Parabolic Approximation to the Wave Propagation”, we retain this term in order to get in line with the approach of other authors to obtaining the parabolic equation. We may notice that to strictly follow the procedure we have to leave the other terms in the right-hand side of Understanding Parabolic Approximation in Wave Propagation: Analytical Methods and Applicationsthis equation in order to be consistent in the accuracy of the approximation. Introducing the same coordinates x = αθ, z = r – α and γ = ϕ · α sin ϑ we obtain:

\frac{\partial^{2} W_{1}}{\partial x^{2}} + j 2 k \frac{\partial W_{1}}{\partial x} + ∆_{⊥} W_{1} + k^{2} [(ε_{m} (z) - 1) + δ ε (\vec{r})] W_{1} = 0 . E q . 34

Then we may introduce the operator:

L = \frac{1}{k^{2}} [∆_{⊥} + k^{2} [(ε_{m} (z) - 1) + δ ε (\vec{r})]] E q . 35

and factorize Eq. (34) in a form similar to Eq. (25):

[\frac{\partial}{\partial x} + j k (1 + \sqrt{1 + L})] \cdot [\frac{\partial}{\partial x} + j k (1 - \sqrt{1 + L})] W_{1} = 0 .

Then, following the above described procedure, we may retain the forward propagating wave obeying the truncated equation:

\frac{\partial W_{1}}{\partial x} = j k (\sqrt{1 + L} - 1) W_{1} E q . 36

with formal solution at the marching step in the form:

W_{1} (x + ∆ x) = W_{1} (x) e x p [j k ∆ x (\sqrt{1 + L} - 1)] . E q . 37

The exponential can then be expanded into a series similar to Eq. (33) and the next range step solution is given by:

W_{1} (x + ∆ x, γ, z) = W_{1} (x) + \sum_{m = 1}^{M} \frac{α_{m} L}{1 + b_{m} L} W_{1} (x, γ, z) . E q . 38

We may introduce the auxiliary function:

f_{m} (x + ∆ x, γ, z) = \frac{α_{m} L}{1 + b_{m} L} W_{1} (x, γ, z) . E q . 39

which can be found by multiplying both parts of Eq. (40) by 1+b_mL and solving a second-order differential equation for f_m. Substituting f_m into Eq. (39) we have:

W_{1} (x + ∆ x, γ, z) = W_{1} (x) + \sum_{m = 1}^{M} f_{m} (x + ∆ x, γ, z) . E q . 40

It is important to notice that all the functions f_m can be solved independently and in parallel at each step along the distance x. The marching equation (40), as well as the forward propagating wave equation (30), shall be appended by both boundary conditions (at z = 0 and

|\vec{r}| \to \infty

) and initial conditions at x = 0.

A somewhat modified approach to the implementation of the Pade method is to substitute a Pade series with a Pade product:

\sum_{m = 1}^{M} \frac{α_{m} L}{1 + b_{m} L} = \prod_{m = 1}^{M} \frac{1 + λ_{m} L}{1 + μ_{m} L} . E q . 41

Now, we concentrate on the split-step Fourier method. Consider Understanding Parabolic Approximation in Wave Propagation: Analytical Methods and Applicationsthis equation and assume the presence of deterministic and stratified refractivity in the medium, the random component δε is set to 0. Assume first that the propagation takes place in an unbounded medium, i. e., we have only requirements on the appropriate decay of the field at infinity. The next assumption normally employed in the split-step method applied to radio wave propagation in the troposphere is to assume the medium to be stratified over the z-coordinate, i. e. a stratified troposphere, thus removing the dependence on the y-coordinate and truncating Understanding Parabolic Approximation in Wave Propagation: Analytical Methods and Applicationsthis equation to a two-dimensional (x, z) parabolic equation.

Understanding Parabolic Approximation in Wave Propagation: Analytical Methods and ApplicationsEquation in this case can be presented in the form:

\frac{\partial W}{\partial x} = j [\frac{k}{2} (ε_{m} (z) - 1) + \frac{1}{2 k} \frac{\partial^{2}}{\partial z}] W . E q . 42

In a vertically stratified medium the formal solution to Eq. (42) is then given by:

W (x + ∆ x) = e x p \{j ∆ x [\frac{k}{2} (ε_{m} (z) - 1) + \frac{1}{2 k} \frac{\partial^{2}}{\partial z^{2}}]\} W (x) E q . 43

where Δx is a range increment. Equation (43) forms the basis for a split-step Fourier implementation. Let us define a Fourier transform of the envelope W(x, z) as follows:

\tilde{W} (x, p) = J [W (x, z)] = \int_{- \infty}^{\infty} d z W (x, z) e x p (- j p z) . E q . 44

Apparently p is a vertical component of the wave vector of the incident field W(x, z). The fundamental assumption in a split-step solution to Eq. (43) is that a Fourier transform of the envelope W(x, z) is performed while treating the term

\frac{k}{2} (ε_{m} (z) - 1)

as a constant, i. e., no dependence on the z-coordinate. With that assumption, we obtain a marching solution:

\tilde{W} (x + ∆ x, p) = \tilde{W} (x, p) e x p [- (k^{2} (e_{m} (z) - 1) - p^{2}) \frac{∆ x}{2 j k}] . E q . 45

The inverse Fourier transform is given by:

W (x + ∆ x, z) = e x p [j \frac{k}{2} (ε_{m} (z) - 1) ∆ x] \cdot J^{- 1} \{\tilde{W} (x, p) e x p (- j p^{2} \frac{∆ x}{2 k})\} . E q . 46

In the above equation the term

\frac{k}{2} (ε_{m} (z) - 1)

is no longer a constant, its variations actually results in variation of the angular spectrum p of the marching field W(x+Δx, z) with the distance.

The boundary condition in the case of ideal reflection (a perfectly conducting boundary surface) can be realized by adding the mirror image of the field W(x, z) in the upper half-space (z > 0) into the lower half space (z < 0) with an appropriate phase to reproduce the odd or even image of the field W relative to z = 0. The odd composition satisfies a condition W(x, z = 0) = 0 for a vertically polarized field, while the even composition results in the boundary condition

\partial W / \partial z |_{z = 0}

for horizontal polarization. In these cases of perfectly conducting boundary the Fourier transform reduces to a one-sided sine or cosine transform, respectively:

{\tilde{W}}_{o} (x, p) = - 2 j J [W_{o} (x, z)] = - 2 j \int_{0}^{\infty} d z W_{o} (x, z) \sin (p z),

{\tilde{W}}_{e} (x, p) = 2 J [W_{e} (x, z)] = 2 \int_{0}^{\infty} d z W_{e} (x, z) \cos (p z) E q . 47

where the subscripts e and o indicate that the field is even or odd, which, in turn, corresponds to horizontal or vertical polarization, respectively. The impedance boundary condition (Understanding Parabolic Approximation in Wave Propagation: Analytical Methods and Applicationsthis equation) was also treated in Ref. where it was shown that such a boundary condition can be realized using a mixed Fourier transform.

{\tilde{W}}_{q} (x, p) = \int_{0}^{\infty} d z W (x, z) [q \sin (p z) - p \cos (p z)]

where

q = \frac{j m}{\sqrt{η + 1}}

(see Understanding Parabolic Approximation in Wave Propagation: Analytical Methods and Applications“Parabolic Approximation to a Wave Equation in a Stratified Troposphere Filled with Turbulent Fluctuations of the Refractive Index”). The marching solution for W(x, z) is then given by:

W (x + ∆ x, z) = e x p [j \frac{k}{2} (ε_{m} (z) - 1) ∆ x] \cdot \{e x p (j q^{2} \frac{∆ x}{2 k} - q z) K (x) + \frac{2}{π} \cdot \int_{0}^{\infty} d p \frac{q \sin (p z) - p \cos (p z)}{q^{2} + p^{2}} e x p (- j \frac{p^{2}}{2 k} ∆ x) \cdot \int_{0}^{\infty} d z^{'} W (x, z^{'}) (q \sin (p z^{'}) - p \cos (p z^{'}))\}

where the term K(x) is defined as:

K (x) = 2 q \int_{0}^{\infty} d z W (x, z) e x p (- q, z); R e (q) > 0, E q . 48

K (x) = 0; R e (q) > 0 .

It is apparent that the implementation of the split-step Fourier algorithm employing the formulas (47) and (48) will be more complicated than the rather simple implementation for a perfectly conducting interface using formula (46). As shown in many other publications (we may refer to a practical implementation of the numerical methods in radio coverage prediction systems), the value of the implementation of the impedance boundary conditions becomes pronounced at vertical polarization at smaller distances and at lower VHF frequencies for long-range propagation, where the impact of the surface wave produced by the distributed images of the source provides a substantial contribution to the received field strength.

The important development in an implementation of the split step Fourier transform method is the capability to treat an irregular terrain profile, as reported in numerous publications. The approach was initially developed in underwater acoustics and then implemented in the studies of radio propagation in the troposphere.

Let us consider the field of horizontal polarization over an irregular terrain with the terrain profile given by the function of the distance T(x). (We may note that with |q| ≫ 1 the boundary condition is W(x, z = 0) = 0 for both polarizations.) The range dependent boundary condition is then given by W(x, z = T(x)) = 0. The approach is then to map the irregular terrain profile to a smooth or more flat one with consequent modifications of the wave equation. The mapping is made by introducing a change of variables. Let the new height and range variables be represented by:

χ = x, ς = z - T (x) . E q . 49

In the new coordinate system the slow varying envelope W(x, z) can be written in the form:

W (x, z) = w (χ, ς) e x p (j θ (χ, ς) E q . 50

where θ(χ, ς) is a phase correction term due to the irregular terrain. Substituting the new variables into Understanding Parabolic Approximation in Wave Propagation: Analytical Methods and Applicationsthis equation and following the procedure we end up with the following equation for a new envelope:

\frac{\partial^{2} w}{\partial ς^{2}} + j 2 k \frac{\partial w}{\partial χ} + k^{2} [(ε_{m} (ς + T (χ)) - 1) + 2 ς T^{''} (χ)] w = 0 E q . 51

and:

θ (χ, ς) = - k ς T^{'} (χ) - k^{3 / 2} \int_{0}^{χ} {[T^{'} (α)]}^{2} d α . E q . 52

As observed, the implementation of the split-step Fourier transform can be done in a way similar to the case of a smooth terrain, the difference from Understanding Parabolic Approximation in Wave Propagation: Analytical Methods and Applicationsthis equation is the presence of the second derivative of the boundary interface, 2ςT^′′(χ), in the exponential operator. Nonetheless, given that the split-step Fourier transform algorithm treats this term as a constant, the addition of the second derivative does not make a difference in the implementation though ensuring the important correction of the wave front due to scattering on the irregular terrain.

We may notice apparent similarities between the split-step Fourier method and the discrete implementation of the path integral method. In the case of propagation over the boundary interface, considered in Wave Field Fluctuations in Random Media over a Boundary Interface“Understanding Wave Field Fluctuations in Random Media”, implementation of the boundary conditions into the path integral approach will result in the appearance of the mirrored source, similar to the above formulas (47) and (48). In fact, numerical calculation of the discrete version of the path integral is practically realised by applying the fast-fourier transform algorithm to the exponential propagator at given step along the range x.

Finally, we can state that the above numerical methods constitute a powerful framework for quantitative analysis of radiowave propagation through the tropospheric boundary layer in the very general case of the refractivity condition. In particular, the methods described above can incorporate either non-uniformity of the refractivity profiles in the horizontal plane as well as in an irregular terrain. The limitation of the above methods, as in fact of all numerical methods, is that while they are quite helpful in a quantitative estimation of the propagation phenomena when the physical mechanism is clearly understood they lack the capability to provide a framework for qualitative analysis, especially in the case of combined effects of refraction and scattering on random irregularities of refractive index or a randomly rough sea surface.

Author

Olga Nesvetailova

Freelancer

Literature

Fock, V.A. Electromagnetic Diffraction and Propagation Problems, Pergamon Press, Oxford, 1965.
Fock, V.A. Distribution of currents excited by a plane wave on the surface of conductor, J. Exp. Theor. Phys., (Nauka), 1945, 15 (12), 693–702.
Leontovitch, M.A. Method of solution to the problems of em wave propagation over the boundary surface, Izv. Acad. Nauk, Ser. Fiz., 1944, 8 (1), 16–22.
Leontovitch, M.A., Fock, V.A. Solution to a problem of em wave propagation over the Earth’s surface using the method of parabolic equation, J. Exp.Theor. Phys., (Nauka), 1946, 16 (7), 557–573.
Fock, V.A. Theory of wave propagation in a non-uniform atmosphere for elevated source, Izv. Acad. Nauk, Ser. Fiz., 1960, 14 (1), 70–94.
Krasnushkin, A.E. The expansion over normal waves in a spherically-stratified medium, Dokl. Acad. Nauk SSSR, 1969, 185 (6), 1262–1265.
Kravtsov, Yu.A., Orlov, Yu.I., Geometric Optics of Non-uniform Media, Nauka, Moscow, 1980, 340 pp.
Rotheram, S. Radiowave propagation in the evaporation duct, The Marconi Rev., 1974, 67 (12), 18–40.
Booker, H.G., Walkinshaw W. The mode theory of tropospheric refraction and its relation to waveguides and diffraction, in Meteorological Factors in Radiowave Propagation, The Royal Society, London, 1946, pp. 80–127.
Andrianov, V.A. Diffraction of UHF in a bilinear model of the troposphere over the Earth’s surface, Radiophys. Quantum Electron., 1977, 22 (2), 212–222.

Footnotes