Discover the intricacies of wave propagation through stratified and random media with our in-depth analysis of parabolic approximation methods. This article covers the application of the wave equation in turbulent tropospheres, highlighting the impact of refractive index fluctuations.
Learn about the Green function for parabolic equations in stratified environments and enhance your knowledge of wave behavior in complex settings.
Analytical Methods in the Problems of Wave Propagation in a Stratified and Random Medium
The analogy between the non-stationary problem of quantum mechanics and the parabolic approximation in a problem wave propagation was likely first explored by Fock. It is well known that the methods of classic wave theory were utilised in quantum mechanics in the early stages of its development. At present, the situation is rather the reverse, and the mathematical methods of quantum mechanics have become widely adopted in the latest developments of wave theory.
The asymptotic solution to the problem of the diffraction of a plane monochromatic wave over a sphere of large radius, compared with the wavelength of radiation, was obtained by Fock in 1945. In obtaining this solution, he developed an asymptotic theory of diffraction on the basis of Leontovich’s boundary conditions and the concept of the “local” field. The approach utilised several large parameters involved in problem formulation, such as |η| ≫ 1, m = (kα/2)1/3 ≫ 1, η = εg+j4πσω, εg, σ, the dielectric permittivity and conductivity of the earth, respectively, ω is the cyclic frequency of the radiated wave, k = 2p/λ is a wavenumber, λ is a wavelength, α is the radius of the earth’s sphere. The same results were obtained by Fock later by means of the parabolic approximation to the wave equation which is based on significantly different parameters, with m ≫ 1, and the characteristic scales of the wave oscillations in two different directions: along the earth’s surface and normal to it. With regards to radio wave propagation at a frequency above 1 GHz the important development was a solution of the parabolic equation for the stratified troposphere, where the refractive index depends solely on the height above the ground surface.
The analysis of the wave field at low altitudes z, which are small compared with the radius of the earth’s curvature α, is convenient to perform by means of the introduction of the modified refractivity:
The solution obtained is suitable for the very general case of the dependence n(z) and impedance boundary conditions at the earth’s surface. The solution is represented by a contour integral which in the shadow region is calculated by a sum over residues in the poles of the integrand, i. e. a series of normal waves. In the line-of-sight (LOS) region the contour integral is commonly calculated using stationary phase methods similar to a geometric optic approximation. Further study of the wave propagation in a stratified tropo- sphere was dedicated to the development of the effective methods of the field representation for various height profiles of the refractivity and combinations of the relative location of the transmitting and receiving antennae. An alternative analytical method to the solution of wave propagation in stratified medium, the method of “invariant submergence”, was developed by Klatskin. Using that method, a boundary problem, which is, in fact, a problem of diffraction, can be transformed into a problem of evolution. A significant advantage of this method is its inherent capability to study the wave propagation through random media, and the method is especially effective when applied to a randomly stratified medium.
Together with a ducting mechanism, the wave scattering on random inhomogeneities of the refractive index plays a significant role in the study of radio wave propa- gation in the atmospheric boundary layer at frequencies above 1 GHz. The straight-forward application of either the contour integral or normal mode series to a problem of wave propagation in a random troposphere will face significant difficulties. The major problem will be related to the divergence of the matrix elements of the normal wave transformation due to a scattering on random fluctuations in the refractive index. The deployment of the method of “invariant submergence” also does not lead to an analytical solution to the problem though numerical calculation is possible.
In this category we use the representation of the Green function in the form of an expansion over the eigen functions of a continuous spectrum which allows one to avoid the above mentioned difficulties of divergence of matrix elements and then to separate the matrix coefficients with a “nearly discrete” and continuum spectrum with normalised height gain functions. This approach led to some advances in solving several problems, in particular, the problem of wave propagation in a tropospheric duct filled with random irregularities of the refractive index in addition to a regular supper-refractive gradient of the refractivity. The disadvantage of the method is that the separation of the discrete spectrum of eigen functions is performed in the “unperturbed”, i. e. regular, part of the Schrödinger operator and, therefore, the spectrum of eigen functions is left unchanged by the random component of the potential term. The study of the spectrum of Schrödinger’s equation with random potential is a complex problem by itself and most advances so far have been achieved in the one-dimensional problem.
The theory of wave scattering in a random medium has been developed during the last decades. The most significant progress has been achieved in the solution to the problems of wave propagation in either a random but unbound medium or a uniform medium with a random boundary interface. In the solution of these problems a statistical approach is commonly used. This approach is targeted at obtaining the statistical characteristics of the scattered field:
- probability function or,
- statistical moments of the field.
Among them the first two moments: average field (coherent component) and either coherence function or field intensity, are of significant interest in many practical applications. There are three known methods of obtaining the closed equations for the field moments: Feynman diagrams, local perturbation and the Markov process approximation.
The study of a coherent component can be done using the modified perturbation theory, a similar result is obtained from the diagram methods. The advances in the solution of the higher order moments of the field have been obtained after application of the parabolic approximation to the wave equation. Neglecting back scattering in the parabolic approximation allows one to employ the concept of the Markov process and all the following analogies from the statistical theory of Markov processes.
The advantage of the functional approach, is that the formal solution to the problem is written via quadrature. The moment of the field of any order can, in principle, be obtained by averaging the respective dynamic solution which is described by a multiple functional integral in the space of the virtual trajectories. The solutions to the first two moments of the field obtained from both the method of closed equations for the moments and the method of the path integral provide exactly the same result, under conditions of the Markov process approximation. The study of the field moments of order higher than two is easier with the use of the path integral, since this approach does not formally require a solution of the equations for the field moments. In fact, this advantage became apparent in the calculation of the moments of the intensity (signal strength), which are not described by a closed equation.
The current limitations of the path integral method could be rather associated with the still limited mathematical methods of analytical calculation of the trajectory integrals, despite the serious studies in that area.
Finally, we comment on studies concerning scattering on a rough sea surface. The basic conclusion is that, until now, there is no adequate theory that takes into account the combined effect of refraction, diffraction and scattering on random inhomegeneities of the refractive index in the volume of the tropospheric layer as well as scattering on a rough sea surface. The impact from a rough surface cannot be regarded as negligible in the general case, at least from a theoretical point of view. Known attempts are limited to the deterministic model of the sea surface or a semi-empirical approach in treating the scattering on the sea surface in the Kirchhoff approximation that is applicable to the coherent component of the scattered field. While the theory of wave scattering on a random surface is an established science by itself, the studies known to the author were concerned with scattering theory on a random surface while treating the propagation media in a very simplistic way, basically, as deterministic and uniform.
In applying the scattering theory to a ducting mechanism, the approach should eventually be modified in order to take into account the multiple effects of refraction, diffraction and scattering. As discussed in Atmospheric Boundary Layer and Basics of the Propagation Mechanisms“Atmospheric Boundary Layer and Key Propagation Mechanisms”, the sea waves could be responsible for the modulation of the turbulence spectrum in the near-surface layer, where the major action of radio propagation takes place. On the other hand, all previous advances in theory have been driven by some unresolved problems in experiments on the radio wave propagation phenomenon. So far, the currently available radio coverage prediction system is reported to be in reasonably good agreement with the results of existing models, at least for measured signal strength.
Parabolic Approximation to a Wave Equation in a Stratified Troposphere Filled with Turbulent Fluctuations of the Refractive Index
Let us consider the vertical electric dipole in a spherical coordinate system (r, ϑ, ϕ) with the spherical axis coming through the dipole. Assume that the dielectric permeability ε is a function of the radius only, e. g. ε = ε(r). The classical solution to the electromagnetic field can be derived by introduction of Debye’s potentials U and V.
The field components are given by the following equations:
We assume the field to be a harmonic function of time – exp (-jωt), ω is a cyclic frequency of the electromagnetic radiation, k = ω/c is a wave number, and c is the speed of light. The notation Δ* represents a Laplace operator on a sphere:
The Maxwell equations:
will be satisfied if potentials U and V obey the following equations:
The approximate Leontovitch’s boundary conditions for U and V have the form:
with r = α, where a is the earth’s radius.
Let us examine the vertically polarised field for which V = 0, U ≠ 0. For high frequencies when kα >> 1, we can select a “dedicated” direction of wave propagation x = αϑ aligned with the direction along the arc of the earth’s radius. With kα >> 1, the radial component of the electric field Er is related to the potential U via: Er ≈ -kα2 U.
Introduce a new function:
which is governed by the equation:
and obeys the boundary condition:
Let us isolate the slow varying complex amplitude W1 in the wave field by:
and introduce coordinate γ = (r sin ϑ)φ. The amplitude W1 is then given by equation:
at the right-hand side of which we have a correction term. Using qualitative arguments, it can be shown that in the case of a stratified medium ε = ε(r) the terms on the left have an order of magnitude of k2W1/m2 while on the right the terms containing derivatives are of the order of k2W1/m4. The right-hand side terms containing sin ϑ in the denominator will also be small if the following inequality holds:
respectively. Having compared the terms in the left- and right-hand sides of Eq. (14) we obtain the inequalities:
which, being satisfied, allow one to treat the terms in the right-hand side of Eq. (14) as negligible. Therefore, with fulfilment of inequalities (17) the right-hand side of Eq. (14) can be replaced with zero.
where cz is a structure constant, κm = 5,92/l0, l0 is an internal scale of turbulence. Substituting Eq. (20) into Eq. (17) we obtain:
Assuming reasonably short distances and high frequencies of radiation, e. g. m ≫ 1 and x ≤ 103km, we can see that the major impact comes from variations of the z-coordinate, and the inequality (22) thus reduces to:
and boundary condition:
The potential U is given by:
The next step is to figure out the correct representation of the field amplitude W1 at small distances from the source in order to obtain an initial condition. In particular, at the distance where the curvature of the earth can be neglected as well as the ray’s refraction we have to obtain a reflection formula in the parabolic approximation. Introducing the coordinates of the source:
- x = 0;
- γ = γ0;
- z = z0;
- and assuming that η ≫ 1 the reflection formula will be satisfied when:
we make a transition in Eq. (27) with x → 0 and obtain an expression for singularity at the source:
It should be noted that in the case of very small heights of the source, another approach can be used and the solution to Eq. (24) should be tailored with the solution of Veil-Van-der-Paul at small distances x.
In conclusion, one should note that there is some ambiguity with the definition of the complex amplitude via Debye’s potential U. Indeed, equations similar to Eq. (24) can be obtained for the other functions W2 and W3 defined by the following expressions:
The area of applicability of the parabolic approximation for W3 is defined by the same set of inequalities (16), (22) while for W3 the inequality (16) can be replaced by a less stringent one:
Let us summarize the major steps undertaken in this section:
- The components of the electromagnetic fieldandcan be represented via Debye’s potentials U and V. Even in the presence of a random turbulent component of the refractive index the de-polarization effects can be neglected and, therefore, the field components are governed by independent equations for U and V.
- The wave equations for both potentials U and V can be approximated by parabolic equations for a slow varying complex amplitude. In this approximation all radiated waves are travelling only in one direction, away from the source. The Leontovich boundary conditions are used along with the parabolic equation. The conditions in the source are also modified in a parabolic approximation.
- During the transformation from the Helmoltz equation to a parabolic equation we actually replaced the earth sphere with a cylinder of the same radius and then flattened the cylinder by means of compensating the earth’s curvature by “ray’s curvature” in the opposite direction. This is a well-known “flat earth” approximation with modified refractivity and locally Cartesian coordinates.;
Green Function for a Parabolic Equation in a Stratified Medium
where:
obeys the boundary condition:
and taking into account the definition of the δ-function, we obtain:
We will seek a solution to Eq. (34) as a composition over functions ΨE(z) which are governed by the equation:
and satisfy the boundary conditions:
Equation (35) is similar to the one-dimensional Schrödinger equation for a particle in a field with potential energy which is proportional to the term (εm(z)-1). As known, in the case of unlimited potential (in Eq. (35) (εm(z)-1)→2z/α with z → ∞)), the motion of the particle in a stationary state is infinite to z → +∞ and the spectrum of the eigen values E is purely continuous.
The eigen functions of a continuous spectrum satisfy the orthogonality conditions:
and “completeness”:
The conditions (38) and (39) are known from quantum mechanics where they are proven for potentials limited at infinite z. Using methods similar to those, one can prove Eqs. (38) and (39) for a potential with a linear increment at infinity. The sign * means a complex conjugate, here and below.
Substituting Eq. (37) into Eq. (34) and taking into account Eqs. (38) and (39) we obtain:
and then performing the integration in Eq. (37) we obtain:
In all practical problems related to radio wave propagation and scattering in the troposphere, the non-uniformities of the refractive index are localised in the boundary layer of the atmosphere and εm(z)-1 can be approximated by a linear function 2z/α with z → ∞.
where contour Γ comes from infinity along the ray arg ξ = -2π/3 to zero and then from zero to infinity ξ → +∞ along the real axis ξ. The function w2(t) is defined by a complex conjugate of Eq. (42).
The S-matrix S(E) has poles En in the upper half-space of E and residues in En completely define the field in the shadow region. The value of ΨE(z) in the pole E = En provides a normalised height function of normal wave with number n:
The normal wave χn(z) with complex value of the propagation constant En(Im En > 0) is unlimited at infinite z. This unlimited growth is a consequence of the exponential decay of the field over the x-coordinate, since the field at z → ∞ is created by radiation propagating from points distant to -∞ over the x-coordinate.
Let us calculate integral (41) for normal refraction when εm(z) = 1+2z/α at z > 0 and the eigenfunctions ΨE(z) the continuous spectrum (43) are defined by:
One can expect that, in a line-of-sight region, expression (41) shall represent a “reflection formula” in accordance with geometric optics, i. e. represent the field as a sum of the direct wave and the wave reflected from the earth’s surface.
The wave propagating in the direction up from the source can be expressed as follows:
Using an asymptotic expression for the Airy function we obtain the trajectory equation by defining the extremum of the phase in Eq. (46):
A similar equation can be obtained for the ray in the direction down from the source:
Differentiating Eqs. (47) and (48) we can find the equation for angle ϑ and its relation to the stationary value of E:
As observed from Eq. (49) the values of E satisfying the inequality E ≤ μ3 z0 determine the real hit-angles ϑ, those, in turn, represent uniform plane waves with angle ϑ/m between the normal to their wavefront and tangential to the radius r = α+z0 at the point of source location. Trajectory (48) has a minimum in:
The ray trajectories described by Eqs. (47) and (48) are shown in Figure.

Let us substitute Eqs. (45)–(47) into Eq. (51) and assume z > z0 and:
The integration in Eq. (51) can be performed using the method of stationary phase. The result is given by:
Substituting Eq. (55) into Eq. (45) we obtain S(E) = 1. Let us introduce the Airy function v(t) by the relationship:
and the variable τ = E/μ2-μz0. The integral (52) can then be transformed into:
Without significant error the low limit in the integral (58) can be extended to -∞. Substituting into Eq. (58) the integral representation for ν(t):
and performing the integration, we obtain:
Equation (60) describes the direct wave in the range of distances:
As a result, we obtain the Green function representation in terms of the normal wave series:
In conclusion, the Green function has been built as an expansion over the set of complete and orthogonal eigenfunctions of the continuous spectrum. In the case of normal refraction, such a representation produces final formulas for the field of the point source that are similar to those introduced by Fock on the basis of the classical approach of the contour integral.
It is worthwhile to emphasise the purpose of the above presentation for a Green function. In a problem of multiple scattering the whole path is involved in the scattering process and one needs to have a representation for the field (or Green function) which is equally applicable in the line-of-sight and shadow regions. Also, in the process of obtaining the closed set of equations for statistical moments of the field, the Green function has to possess some properties of orthogonality and completeness which can be achieved by using the above expansion over the set of eigenfunctions of the continuous spectrum.