© D. Vasylyev et al., Published by EDP Sciences 2022

1 Introduction

The scintillation phenomenon is familiar to everyone who has looked at the sky during cloudless nights. By looking for some period of time, one can observe the remarkable twinkling of stars. These random changes in the apparent star brightness or its position are due to the random distortion of light coming from a star caused by refractive index irregularities in a turbulent troposphere. The fluctuation in brightness or, in other words, in intensity is one example of the scintillation phenomenon. The pulsation of star images (star dancing) is another manifestation of this phenomenon known as angular scintillation.

If we were able to sense the electromagnetic radiation at radio wavelength, we would similarly perceive scintillation of distant radio sources of celestial origin and from artificial satellites. In this case, scintillation is caused by random scattering on irregularities of electron density of the electrons and protons in the Earth ionosphere, solar wind, or in the interstellar medium. As the electron density varies as a random function of location, the associated refractive index changes similarly. Hence, as a plane wavefront of the incoming signal passes through the region where the refractive index is low, the velocity of radio waves will be high, and the wavefront will advance farther. With a higher refractive index, the wavefront is delayed. In combination, these effects will result in a corrugated wavefront and phase fluctuations. Some portions of the wavefront will be convex in the direction of propagation, and some of them will be concave. At some distance from the scattering region, the partial waves with such concave or convex wavefronts will interfere constructively or destructively, producing a complex intensity pattern.

This article will focus on ionospheric scintillation phenomena, i.e., amplitude or intensity and phase fluctuations. Amplitude scintillation is more intense in equatorial regions, while phase scintillation is more frequent and stronger at high latitudes (Aarons, 1982; Liu et al., 2012; Tsai et al., 2017). Amplitude scintillation normally occurs when the Fresnel dimension of the propagating radio signal is of the order of irregularity scales in the ionosphere. This condition is met for the post-sunset equatorial irregularities formed in equatorial F-spread and the equatorial plasma bubbles (Moorthy et al., 1979; Patel et al., 2009; Yokoyama & Stolle, 2017). The phase scintillation is primarily caused by steep ionospheric density gradients and irregularities formed in the cusp and auroral precipitation regions and polar cap patches (Spogli et al., 2009; Jiao & Morton, 2015). Fluctuations in phase in these regions are additionally enhanced due to geometric effects, i.e., the relative position of field-aligned irregularities and the signal path (Forte & Radicella, 2004). The occurrence of ionospheric scintillation is also strongly time-dependent. It varies daily, seasonal, and yearly and is more frequent during post-sunset hours at low latitudes (Basu et al., 1996; Huang et al., 2014a; Jin et al., 2018). During vernal and autumn equinoxes, there is a higher probability of scintillation (Muella et al., 2013; Jiao & Morton, 2015). Finally, the strength and occurrence of phase and intensity scintillation correlate with solar and geomagnetic activity at high latitudes (Abdu et al., 1998; Souto Fortes et al., 2015; Guo et al., 2017; Béniguel, 2019).

Intense phase and intensity fluctuations may have a strong negative impact on modern technical infrastructure, the robustness of some life-critical services, and the quality of scientific data acquired in remote-sensing measurement campaigns. Intensity scintillation causes the signals from a satellite to fade below the average level. If the level of fading exceeds some prescribed margin of a receiver, the signal becomes masked by the noise and cannot be reliably retrieved. Phase scintillation induces a frequency shift that is responsible for cycle slips in GNSS signals (Knight & Finn, 1998; Horvath & Crozier, 2007). When this shift exceeds the bandwidth of the phase lock loop, the satellite tracking might fail, and in some cases, the loss of lock with the satellite might follow (Conker et al., 2003; Aquino et al., 2005; Sreeja et al., 2011). For positioning, navigation, and timing (PNT) applications, such as precise point positioning (PPP), this effect can impair the functionality of this service, putting vital systems at risk.

In astronomical and remote sensing applications, scintillation may degrade the contrast and resolution of images as well as incorporate artifacts and increase the noise level in retrieved data. In radio astronomy, ionospheric and interplanetary scintillations lead to image wander, degradation of image contrast, pulse broadening of pulsar signals, angular broadening of point sources, and thus degradation of the resolving power of radio telescopes (Strom et al., 2001). In applications such as GNSS reflectometry, intensity scintillation has the greatest impact on the degradation of retrieval performance in altimetric and scatterometric applications (Camps et al., 2016). As the algorithms for the sounding of Earth and planetary atmospheres are usually based on ill-posed inversion problems (Leitinger, 2001), the introduction of additional sources of noise due to scintillation may considerably influence the performance of retrieval methods. In occultation measurements, strong refractive scintillation impacts the measured values of bending angle and contributes to the residual ionospheric error (Hinson, 1986; Carrano et al., 2011; Ludwig-Barbosa et al., 2020). The remarkable effect of intensity and phase scintillation on synthetic aperture radar (SAR) missions is the “azimuth streaking” notable as a spatial modulation in the position of optimal azimuth registration (Gray et al., 2000; Carrano et al., 2012; Sato et al., 2021). Other effects of scintillation that influence the performance of SAR are defocusing (Belcher, 2008; van de Kamp et al., 2009), signal decorrelation (Ishimaru et al., 1999), degradation of coherence properties Béniguel & Hamel (2011), and random phase perturbations (Rogers et al., 2014).

However, it should be noted that scintillation has not only negative effects. The knowledge of the scintillation level makes it possible to localize irregularities of the electron density in the ionosphere (Sokolovskiy, 2000), deduce the climatological trends for their occurrence (Siefring et al., 2011), retrieve the properties of ionospheric, interplanetary, and interstellar medium (Rickett, 1977; Bhattacharyya et al., 1992; Shishov et al., 2010), or even infer the structure of the local magnetic fields in planetary ionospheres (Hinson & Tyler, 1982; Hinson, 1984). In combination with other instrumentation, scintillation monitoring becomes a powerful tool for investigating the formation and evolution of ionospheric structures such as spread-F and equatorial plasma bubbles (Kelley et al., 2002).

Due to the aforementioned impact of scintillation on various important services and on remote sensing data products, the theoretical modeling and computer simulation of these phenomena were attracting the attention of researchers for more than seventy years (for review, see e.g. Ratcliffe, 1956; Gurvich & Tatarskii, 1975; Fante, 1975a; Crane, 1977; Yeh & Liu, 1982; Kravtsov, 1992; Priyadarshi, 2015). The first theoretical ideas on how to model scintillation were closely related to the emergent statistical theory of turbulence (Kolmogorov, 1941) and the development of the theory of random functions and corresponding methods of spectral analysis (see e.g. Rice, 1944, 1945). The former theory gave an understanding of how one can describe complex random systems such as turbulence or randomly inhomogeneous media, while the theory of random functions provided the rigorous mathematical apparatus for dealing with signals modulated randomly in phase and amplitude. The major theoretical milestones in such developments were connected with the need for a description of optical light propagation in tropospheric turbulence or sound waves in a turbulent ocean (Tatarski, 2016). To such landmarks belong the introduction of the concept of locally homogeneous random media (Obukhov, 1949), the development of Rytov’s theory of smooth perturbations (Obukhov, 1953), the derivation of field moment equations (Shishov, 1968), and asymptotic methods for solving such equations (Gochelashvily & Shishov, 1975; Yakushkin, 1975; Fante, 1975b). These concepts allow one to establish a rigorous mathematical apparatus for dealing with wave propagation in a random medium and obtain formulas of scintillation indices for multiple cases important in practice.

In connection with the ionospheric propagation, the first theoretical attempts to explain scintillation were associated with the diffraction theory (Ratcliffe, 1956). The propagation of electromagnetic waves in a random ionosphere has been proposed to be modeled via a thin sheet of material that modulates the wave amplitude and phase. The transmitted wave is then represented as a set of uniform infinite plane waves traveling in different directions. The values of corresponding directional cosines are considered to be governed by the specific angular spectrum (Booker & Clemmow, 1950). At some distance from the modulating sheet of material, the interference of plane waves forms a specific diffraction pattern in intensity. If the modulation strength of the thin sheet is a random function of position, the phase shifts of plane waves emanating from it, as well as the formed diffraction pattern are random functions of the position (Booker et al., 1950). In these initial approaches, the randomness of scattering directions, i.e., of the directional cosines values, is the driver of random modulation for both phase and amplitude directly behind the screen. However, it was soon realized that for most cases of radio wave propagation in the ionosphere, the amplitude modulation factor could be neglected, and the spectrum of random phase modulation is sufficient to model scintillation (Hewish, 1951). This concept of the phase-modulated sheet, called simply the phase screen, has become an extremely handy tool for studying and modeling ionospheric phase and amplitude scintillation.

Among further advances and milestones in phase screen modeling of scintillation, one can distinguish the incorporation of spatial anisotropy of field-aligned ionospheric irregularities in phase screen generation (Briggs & Parkin, 1963; Singleton, 1970a), the introduction of multiple-phase screens for a description of extended random media (Fejer, 1953), and the replacement of Gaussian spectra for phase screen generation with more realistic ones with power-law dependence (Rufenach, 1972). Further, Uscinski (1985) has shown that there exists a close resemblance between the methods of moment equations for electromagnetic wave amplitudes and the theory of multiple-phase screens. The popularity of phase screen simulations of scintillation was also growing due to advances in numerical mathematics. After developing the fast Fourier transform (FFT) algorithm (Cooley & Tukey, 1965) and the split-step algorithm for solving wave equations with a variable refractive index (Hardin et al., 1973), the FFT-based phase screen generation became spread out as the numerical method for both intensity and phase scintillation simulations.

There are many models of ionospheric scintillation based on the aforementioned fundamental concepts, either on heuristic ideas or on phenomenology. Many of them are collected and discussed in some detail in the review of Priyadarshi (2015), which briefly discusses scintillation mechanisms and the morphology of high- and low-latitude scintillation. Among analytical models, one should mention those developed by Bramley (1967), Fremouw & Rino (1973), Rino (1979a,b), Aarons (1985), Franke & Liu (1985), Hinson (1986), Retterer (2010). The prominent numerical climatological models are WBMOD¹ (Fremouw & Secan, 1984; Secan et al., 1995, 1997), GISM² (Béniguel, 2002; Béniguel & Hamel, 2011; Béniguel, 2019), SIGMA³ (Deshpande et al., 2014, 2016; Deshpande & Zettergren, 2019), and some related models for GNSS applications such as the Cornel model (Humphreys et al., 2009), the UPC/OE/RDA model (Camps et al., 2017), and the multi-frequency GNSS scintillation model (Rino et al., 2018). The models incorporating in situ data are the equatorial and high-latitude models of Basu et al. (1976, 1981), the WBMGRID⁴ model (Groves et al., 1997), and the Wernik-Alfonsi-Materassi (WAM) model (Wernik et al., 2007), to name just a few. Table 1 summarizes the major development stages of scintillation modeling from a historical perspective. Many mentioned models are used for scientific purposes as well as in services for scintillation prediction, determination of time periods and geographic areas with strong scintillation, warning systems, etc.

The present article reviews the basic methods and theories used in the modeling of ionospheric scintillation. Instead of reviewing the ideas and results for particular scintillation models, we have the intention to provide the basic theoretical framework for describing and modeling electromagnetic field propagation in a random ionosphere that is applicable to the majority of the models listed in Table 1. The principal focus of this article is on the scintillation modeling with phase screens, namely the single- and multiple-phase screen approaches. Section 2 presents the theoretical methods for describing the propagation of electromagnetic waves in a disordered medium. Here we derive the wave equation for an electromagnetic field propagating in the random ionospheric media. The solution to this equation can be obtained by effectively superposing the wave propagation in a vacuum together with the random modulation of its phase. Particularly, the statistical properties of this phase modulation are deduced from the power spectral densities of electron density fluctuations. Several popular models of such spectra are summarized in Section 3. Section 4 outlines the influence of geometric parameters of the communication link and scintillation-causing irregularities on the strength of the scintillation. We discuss in some detail the assumptions of flat- and spherical geometries in the problem of scattering on random ionospheric inhomogeneities. Geometric parameters such as zenith and azimuth angles of the link, shape, and direction of anisotropic ionospheric irregularities and propagation factors are discussed. The analytic results for scintillation indices that include these effects are given in Section 5. In particular, the indices obtained in the flat-geometry approximation serve as a basis of the WBMOD scintillation model. The provided results in the spherical-geometry case can be useful for the consistent extension of the WBMOD applicability to the arbitrary slant communication links. The extension of these results to the case of multiple scattering is given in Section 6, where we discuss the FFT-based multiple phase screen (MPS) simulation algorithm. The Global Ionospheric Scintillation Model (GISM) is used to exemplify the scintillation simulation via the MPS. Finally, some conclusions have been drawn in Section 7.

2 Theoretical models of scintillation

When propagating in a random medium, an electromagnetic wave undergoes chaotic distortions of its phase and amplitude (or, alternatively, its power or intensity). In order to analyze, predict and effectively mitigate the degradation of the transmitted signal encoded in such a wave, the probability distribution functions (PDFs) for random phase and amplitude are needed. In general, the construction of such PDFs requires knowledge of all their moments or cumulants. In practice, it appears that moments up to the second-order are of importance and could alone serve as proxy parameters characterizing affects of random media on the propagating signal. This restriction of the interest to low-order moments is driven by two factors. Firstly, for a strongly irregular medium, one would expect that the PDF of interest approaches the Gaussian distribution as a consequence of the limiting theorems of probability theory (Papoulis & Pillai, 2002). In this case, knowledge of the first and second-order moments is sufficient to determine the total PDF. Secondly, the most important information is generally contained in low-order moments. For example, the mean-field provides information on the way that the incident wave is attenuated on passing through the medium. The second field moments are related to the visibility of fringes of two interfering signals, yielding information on the angular spread and power of the propagating signal. The second intensity moment (fourth field moment) is useful for the characterization of intensity fluctuations of the transmitted wave.

For this reason, the following statistical moments and their combinations are of primary interest in studying the effects of a randomly inhomogeneous ionosphere on the propagation of radio and microwave signals:$$S_4^2=\frac{\langle I^2\rangle -\langle I{\rangle }^2}{\langle I{\rangle }^2},\hspace{1em} I=|E{|}^2,$$(1)$${\sigma }_S^2=\langle S^2\rangle -\langle S{\rangle }^2,$$(2)where E is the field amplitude, I is the field intensity, S is the phase departure of the wave detected by the observer. The statistical averaging (ensemble- or time-averaging) here and in the following is denoted as 〈...〉. Conventionally, the one-minute time averages are used for the scintillation metrics defined by equations (1) and (2). The proxy quantities S₄ and σ_S are called the amplitude (or intensity) and phase scintillation indices, correspondingly. As can be seen from the definition, both indices provide information on how strong fluctuating intensity and phase differ from their regular or mean values.

Restricting our attention to the analysis of the second-order statistical moments (1), (2) is one of the simplifications made while treating scintillation phenomena from the theoretical perspective. In the following sections, we will make a number of further simplifications. We will neglect the bending of signal rays in the ionosphere and consider the depolarization effects caused by the propagation in the ionosphere to be negligibly small. For simplicity, the spherical shape of the Earth is assumed, but the flat Earth approximation is discussed in the appropriate context. The ionospheric plasma would be considered to be collisionless, its random variations would be treated by statistical means that again rely on second-order statistical quantities such as autocorrelation functions. The simplified model for autocorrelation properties of electron density fluctuations would be utilized for incorporating the anisotropy effects in scintillation phenomena due to scattering on field-aligned irregularities. The bulk of hte random ionosphere would be treated effectively as a set of thin random phase screens, whose temporal variation would be incorporated into the spatial variation with the help of the so-called Taylor hypothesis. In view of all these assumptions and simplifications, the resulting theoretical methods are quite capable of gaining insights into the nature of scintillation phenomena and to be used in computer simulations of scintillation fading levels.

In this section, we outline the theoretical foundations of signal wave propagation in a random ionospheric medium and discuss some approximations and assumptions used. In particular, we derive the wave equation for a signal wave propagating in a randomly inhomogeneous medium, which is the central equation in the forthcoming analysis. The solution of this equation is essential for the calculation of the scintillation indices, as will be shown in the following sections.

2.1 Wave propagation in random medium

Let us consider a plane wave propagating in z direction. This wave impinges a random medium that extends over the range 0 ≤ z ≤ Δz, cf. Figure 1. In this section, we restrict our consideration to the normal incidence case that will be generalized to the oblique incidence case in Section 4. The electromagnetic field of this wave we write as$$\mathbf{E}(\mathbf{r},t)=\{\begin{array}{ll}{\mathbf{E}}_0{e}^{i(\mathrm{kz}-\omega t)},& z<0,\\ \mathbf{f}\left(\mathbf{r}\right)e^{-\mathrm{i\omega }t},& 0\le z\le \Delta z,\end{array}$$(3)where r = (x y z)^T is the three-dimensional spatial coordinate vector, ω is the angular frequency, k = ω/c is the free-space wave number, c is the speed of light. The quantities E₀ and f(r) are the time-independent field amplitudes for waves propagating in a vacuum and in a random medium, respectively.

Figure 1 Geometry of propagation through a thin (a) and extended (b) random medium. These regimes correspond to a situation when the medium thickness is comparable to or larger than the typical scale $\mathcal{l}$ of medium inhomogeneities.

If ε is the dielectric permittivity of the random medium, its general expression can be written as$$\epsilon (\mathbf{r},t)=\langle \epsilon (\mathbf{r},t)\rangle \left[1+\delta \epsilon (\mathbf{r},t)\right],$$(4)where 〈ε〉 is the mean dielectric permittivity and δε is the dimensionless random modulation factor possessing the property 〈δε〉 = 0. With the use of Taylor’s hypothesis (Taylor, 1938) the time dependence of the dielectric permittivity can be incorporated into the spatial coordinate. Taylor argued that if the velocity fluctuations of a random medium are small compared to the mean flow velocity, the temporal variation at any fixed point in space can be considered as a passage of a “frozen” spatial structure past the observer at the mean flow velocity v_drift. In this case, the time dependence can be incorporated in the spatial coordinate, i.e., ε(r, t) = ε(r − v_drift t,0). For convenience, in the following, we will simply write the spatial dependence as ε(r). In many practical situations, the Taylor hypothesis can be applied with relatively mild restrictions (Tennekes, 1975) and is suitable for analyzing ionospheric plasma turbulence (Mikkelsen & Pécseli, 1980; Dyrud et al., 2008; Rino, 2011).

For electromagnetic waves with frequencies much larger than the plasma frequency (i.e. larger than 1–1.5 MHz) propagating in the collisionless plasma, the background dielectric permittivity can be written as (Rawer, 1993)$$\langle \epsilon \left(\mathbf{r}\right)\rangle ={\epsilon }_0\left[1-\frac{e^2\langle N_e\left(\mathbf{r}\right)\rangle }{{\epsilon }_{0}m {\omega }^2}\right]={\epsilon }_0\left(1-\frac{{\omega }_p^2}{{\omega }^2}\right),$$(5)where ε₀ is the vacuum permitivity, ω_p is the plasma frequency, e is the elementary charge, N_e is the electron density and m is the electron mass. From this equation follows that the presence of free electrons modifies the free-space dielectric permittivity and makes it dependent on the frequency of the propagating electromagnetic wave. In the following we use the refractive index $n\left(\mathbf{r}\right)=\sqrt{\epsilon \left(\mathbf{r}\right)/{\epsilon }_0}$ instead of the dielectric permittivity.

For a weakly inhomogeneous random ionosphere (Rino, 2011) the depolarization effects can be neglected such that each polarization component f_p (r), p = 1, 2 of the field f(r) can be treated separately. For simplicity in the following, we omit the index p associated with certain polarization components. If we assume that the signal wave has a beam-like form with a definite propagation direction, which in our case is along the z axis, we can separate the fast and slowly varying field components as f(r) = u(r)e^ikz. The fast varying component e^ikz will not be of much interest in the following. Indeed, this term is eliminated in the definition of the intensity scintillation index, and in phase measurements it is usually treated as the phase reference; thus, it is also not relevant for a calculation of the phase scintillation index. The wave equation for the slowly varying amplitude u(r) can then be written in the so-called paraxial approximation as$$2\mathrm{ik}\frac{\mathrm{\partial }u({\mathbf{r}}_{\perp },z)}{\mathrm{\partial }z}+{\Delta }_{\perp }u({\mathbf{r}}_{\perp },z)+2k^2\delta n({\mathbf{r}}_{\perp },z) u({\mathbf{r}}_{\perp },z)=0,$$(6)where Δ_⊥ = ∂²/∂x² + ∂²/∂y² is the component of the Laplace operator transversal to the propagation direction and r_⊥ = (x y)^T is the corresponding transversal coordinate vector. In the paraxial approximation, we have neglected the term with second-order derivatives with respect to z since for the directed beam, it is much smaller than the term k∂u/∂z. The fluctuating part of the refractive index in equation (6) is$$\delta n=-\frac{r_e{\lambda }^2}{2\pi }\delta N_e,$$(7)where r_e ≈ 2.818 × 10⁻¹⁵ m is the classical electron radius and λ = 2π/k is the wavelength of the electromagnetic wave.

The wave equation (6) is the starting point for a derivation of statistical quantities defined in equations (1) and (2). In this connection, various approximative schemes have been developed for the solving of (6), such as the geometric optics approximation (Tatarskii, 1971; Wheelon, 2004), the method of smooth perturbations, or Rytov approximation (Clifford, 1978; Wheelon, 2003), the Huygen–Kirchhoff method (Ratcliffe, 1956; Banakh & Mironov, 1979), the method of field moments (Tatarskii, 1971; Liu et al., 1974; Prokhorov et al., 1975), and the method of functional integrals (Tatarskii et al., 1992; Zavorotny et al., 1992), to name just a few.

2.2 Split-step approach and the concept of a random phase screen

The popular solution method of equation (6) is the so-called split-step algorithm (Hardin et al., 1973; Tappert & Hardin, 1975). It is based on the observation that the wave equation is a combination of the terms responsible for wave propagation in a vacuum [the first and second terms in the L.H.S of Eq. (6)] and in a random medium with the fluctuating refractive index δn [the first and the third terms in the L.H.S of Eq. (6)]. Now, if the spatial scale along the z axis, over which the refractive index fluctuation changes considerably, is denoted by $\mathcal{l}$, two situations might occur. If the thickness of the random medium $\Delta z\approx \mathcal{l}$, cf. Figure 1a, the medium slab can be effectively replaced with a thin randomly phase changing screen. In the case of an extended random medium with $\Delta z>l$, cf. Figure 1b can be effectively represented as a set of thin phase screens placed in a vacuum. We refer to the first case as the weak scattering regime and the second case as the multiple scattering regime. The combination of the split-step method with the model of single or multiple phase screens is a popular tool to solve equation (7), which we discuss in some detail in Section 6.1.

In order to see why the concept of the random phase screen appears while solving equation (6), let us consider the propagation along the infinitesimally small distance $\delta z\ll \mathcal{l}$ in the medium. The distance is chosen in such a way that δn(r_⊥, z + δz) ≡ δn(r_⊥, z), i.e., the fluctuating part of the refractive index, does not change sufficiently. While propagating this distance, the field gains a phase increment k δn δz and can be written as$$u({\mathbf{r}}_{\perp },z+\delta z)=u({\mathbf{r}}_{\perp },z)e^{\mathrm{ik} \delta n \delta z}\approx u({\mathbf{r}}_{\perp },z)\left(1+\mathrm{ik} \delta n \delta z-\frac{1}{2}(k \delta n \delta z{)}^2+...\right).$$(8)

By forming the combination$$2\mathrm{ik}\frac{\mathrm{\partial }u}{\mathrm{\partial }z}=2\mathrm{ik}\underset{\delta z\to 0}{\lim }\frac{u({\mathbf{r}}_{\perp },z+\delta z)-u({\mathbf{r}}_{\perp },z)}{\delta z}=-2k^2 \delta n u,$$(9)we see that it is exactly the third term in the L.H.S. of equation (6) up to the minus sign. In this way, the medium-related part of the derivative ∂u/∂z cancels the third term in (6), and only the vacuum propagation wave equation is left over. We can thus interpret the paraxial equation as the propagation in a vacuum superimposed with the random phase modulating effect due to the spatially variable refractive index fluctuation δn.

The cumulative phase gain while propagating in a random medium of thickness Δz is$$\delta \phi \left({\mathbf{r}}_{\perp }\right)=k\underset{z}{\overset{z+\Delta z}{\int }}\mathrm{}\delta n({\mathbf{r}}_{\perp },z^{\mathrm{\prime }}) \mathrm{d}{z}^{\mathrm{\prime }}.$$(10)

The second-order statistical properties of random phase fluctuations are inferred from the autocorrelation function, which can be written for a spatially homogeneous medium as$$\begin{array}{c}{\rho }_{\delta \phi }\left({\mathbf{r}}_{1,\perp }-{\mathbf{r}}_{2,\perp }\right)=k^2\underset{z_1}{\overset{z_1+\Delta z}{\int }}\mathrm{}\underset{z_2}{\overset{z_2+\Delta z}{\int }}\mathrm{}\langle \delta n\left({\mathbf{r}}_{1,\perp },z^{\mathrm{\prime }}\right)\delta n\left({\mathbf{r}}_{2,\perp },z\mathrm{\text{'}\text{'}}\right)\rangle \mathrm{d}z\mathrm{\prime } \mathrm{d}{z}^{\mathrm{\text{'}\text{'}}},\\ =k^2\Delta z\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\mathrm{}{\rho }_{\delta n}({\mathbf{r}}_{1,\perp }-{\mathbf{r}}_{2,\perp },z)\mathrm{d}z,\end{array}$$(11)where for any (real) quantity X the autocorrelation is usually defined as$${\rho }_X({\mathbf{r}}_1,{\mathbf{r}}_2)=\langle X\left({\mathbf{r}}_1\right)X\left({\mathbf{r}}_2\right)\rangle -\langle X\left({\mathbf{r}}_1\right)\rangle \langle X\left({\mathbf{r}}_2\right)\rangle .$$(12)

The second line in equation (11) is obtained by changing the integration of coordinates according to z = z′ − z″, 2Z = z′ + z″ and the subsequent integration over Z. The limits of the integration were extended to infinite values since the correlation function ρ_δn decays quickly, as the argument exceeds the correlation length $\mathcal{l}$, which in turn does not exceed the medium thickness Δz.

The relation of autocorrelation functions for phase and refractive index fluctuations (or, equivalently, of the electron density fluctuations) provided by equation (11) suggests that the knowledge of statistical properties of the random medium is useful for generating physically sound phase screens. It appears that for phase screen generation, it is more convenient to operate with the Fourier images of the autocorrelation functions. The latter quantities, known as the power spectral densities (PSDs), and underlying analytical models are discussed in the next section.

3 Statistical description of ionospheric random medium

The autocorrelation function for an arbitrary real quantity X, as defined in equation (12), has the Fourier image$${\mathrm{\Phi }}_X\left(\mathbf{\kappa }\right)=\int \mathrm{}{\rho }_X\left(\mathbf{r}\right)e^{-i\mathbf{\kappa }\cdot (\mathbf{r}-{\mathbf{r}}^{\mathrm{\prime }})} \frac{{\mathrm{d}}^3\mathbf{r}}{(2\pi {)}^3}=\int \mathrm{}{\rho }_X\left(\mathbf{r}\right)\cos \left(\mathbf{\kappa }\cdot \mathbf{r}\right) \frac{{\mathrm{d}}^3\mathbf{r}}{(2\pi {)}^3},$$(13)which is the three-dimensional PSD of process X. The dimensionality of the PSD refers to the three-dimensional vector of spatial frequencies κ = (κ_x κ_y κ_z)^T. Examples of quantity X that are of interest in the context of wave propagation in the ionosphere are the electron density fluctuations, δN_e, or, equivalently, the reflective index fluctuations, δn, the random phase increment due to the transmission through a medium slab, δφ, the phase departure of the wave, S, observed by a ground observer, and the fluctuations of the amplitude δE/〈E〉 or of the log-amplitude χ = log(1 + δE/〈E〉). In the following, we will substitute the corresponding quantity instead of X in the specified contexts.

In practical situations, not all of the information on all components of the vector κ is usually available. For example, if two or more ground stations record the diffraction pattern of a signal from a satellite, the correlation analysis can be performed to infere the properties of ionospheric irregularities (Briggs, 1984; Costa et al., 1988; Kriegel et al., 2017). The statistical quantity inferred is then defined in a plane rather than in the three-dimensional space. Similarly, the satellite or rocket measurements usually yield one-dimensional information along their flight track (Inhester et al., 1990). In all these situations, the full three-dimensional spectrum (13) is not accessible, but rather its two- or one-dimensional analogies can be inferred. One therefore defines the two-dimensional PSD as (Tatarskii, 1971)$$F_X({\mathbf{\kappa }}_{\perp },z)=\int \mathrm{}{\mathrm{\Phi }}_X({\mathbf{\kappa }}_{\perp },{\kappa }_z)\cos \left({\kappa }_{z}z\right) \mathrm{d}{\kappa }_z, {\mathrm{\Phi }}_X\left(\mathbf{\kappa }\right)=\frac{1}{2\pi }\int \mathrm{}{F}_X({\mathbf{\kappa }}_{\perp },z)\cos \left({\kappa }_{z}z\right) \mathrm{d}z$$(14)and the one-dimensional one as (Kovasznay et al., 1949)$$W_X\left(\kappa \right)=\frac{1}{\pi }\underset{0}{\overset{\mathrm{\infty }}{\int }}\mathrm{}{\rho }_X\left(r\right)\cos \left(\kappa r\right) \mathrm{d}r, {\mathrm{\Phi }}_X\left(\kappa \right)=-\frac{1}{2\pi \kappa }\frac{\mathrm{d}{W}_X\left(\kappa \right)}{\mathrm{d}\kappa }.$$(15)

The last relation in equation (15) is valid if the quantity X is spatially isotropic, i.e. its variation does not depend on some preferable direction in space.

The isotropy condition tends to be very restrictive and cannot be fully met in the ionosphere, where the electron density irregularities tend to align along the geomagnetic lines of force. The spatial homogeneity condition used in equation (11) is less restrictive and means that the correlation functions depend on the relative distance between two spatial points, where the quantities of interest are recorded, rather than on the exact point coordinates. In general cases, this holds if we are interested in spatial separations smaller than the characteristic sizes of ionospheric irregularities. Since the medium thickness, Δz, is larger than such sizes, it is useful to introduce the concept of local homogeneity. This concept follows from the observation that although the quantity X(r) might not be spatially homogeneous, the difference X(r) − X(r′) as the detrended process might possess this property. So instead of the autocorrelation functions, the locally spatially homogeneous process X is characterized by the structure function (Wheelon, 2004; Rino, 2011; Tatarski, 2016)$$D_X(\mathbf{r}-{\mathbf{r}}^{\mathrm{\prime }})=\langle {\left[X\left(\mathbf{r}\right)-X\left({\mathbf{r}}^{\mathrm{\prime }}\right)\right]}^2\rangle =2\left[{\rho }_X\left(0\right)-{\rho }_X(\mathbf{r}-{\mathbf{r}}^{\mathrm{\prime }})\right].$$(16)

In particular, by using equations (11) and (14) the following relationship can be established between the structure functions of phase and refractive index fluctuations:$$\begin{array}{c}{D}_{\delta \phi }\left({\mathbf{r}}_{\perp }-{\mathbf{r}}_{\perp }^{\mathrm{\prime }}\right)=4\pi k^2\Delta z\int \mathrm{}{\mathrm{\Phi }}_{\delta n}\left({\mathbf{\kappa }}_{\perp },{\kappa }_z=0\right)\left\{1-\cos \left[{\mathbf{\kappa }}_{\perp }\cdot \left({\mathbf{r}}_{\perp }-{\mathbf{r}}_{\perp }^{\mathrm{\prime }}\right)\right]\right\} {\mathrm{d}}^2{\mathbf{\kappa }}_{\perp },\\ =2k^2\Delta z\int \mathrm{}\int \mathrm{}{F}_{\delta n}({\mathbf{\kappa }}_{\perp },z)\left\{1-\cos \left[{\mathbf{\kappa }}_{\perp }\cdot ({\mathbf{r}}_{\perp }-{\mathbf{r}}_{\perp }^{\mathrm{\prime }})\right]\right\} {\mathrm{d}}^2{\mathbf{\kappa }}_{\perp } \mathrm{d}z.\end{array}$$(17)

Thus, the local homogeneity of refractive index fluctuations implies the local homogeneity of random phase increments for the propagating wave.

In the following, we will primarily use the structure functions instead of the autocorrelation functions. In this way, we emphasize that the random medium under consideration is not restricted to being spatially homogeneous, and the obtained formulas are valid in less restricted cases of local homogeneity. Additionally, as will be discussed in Section 6.1, the phase structure function (17) is a convenient metric to check if the relevant parameters in a computer simulation of scintillation with phase screens have been properly chosen.

We now list the most widely used spectral models for modeling ionospheric irregularities. For now, we will concentrate merely on the functional dependence of these models and discuss the inclusion of anisotropy effects in the next section.

3.1 Gaussian model

The simplest idea is to assume that the refractive index (electron density) fluctuations are correlated according to Gaussian law with some characteristic correlation radius r₀. The corresponding spectra are$$\begin{array}{c}{\mathrm{\Phi }}_{\delta n}\left(\mathbf{\kappa }\right)=\frac{r_0^3\langle \delta n^2\rangle }{8\pi \sqrt{\pi }}\exp \left(-r_0^2|\mathbf{\kappa }{|}^2/4\right),\\ F_{\delta n}({\mathbf{\kappa }}_{\perp },z)=\frac{r_0^2\langle \delta n^2\rangle }{4\pi }\exp \left(-\frac{z^2}{r_0^2}-r_0^2{\kappa }_{\perp }^2/4\right),\\ W_{\delta n}\left(\kappa \right)=\frac{r_0\langle \delta n^2\rangle }{2\sqrt{\pi }}\exp \left(-r_0^2{\kappa }^2/4\right).\end{array}$$(18)

The spectra are properly normalized such that$$\int \mathrm{}{\mathrm{\Phi }}_{\delta n}\left(\mathbf{\kappa }\right) {\mathrm{d}}^3\mathbf{\kappa }=\int \mathrm{}{F}_{\delta n}({\mathbf{\kappa }}_{\perp },0) {\mathrm{d}}^2{\mathbf{\kappa }}_{\perp }=\int \mathrm{}{W}_{\delta n}\left(\kappa \right) \mathrm{d}\kappa =\langle \delta n^2\rangle$$(19)holds true. The Gaussian correlation function and corresponding spectral models were popular in the early stages of theoretical descriptions of scintillation in the ionosphere (Bramley, 1954; Bowhill, 1961; Briggs & Parkin, 1963; Liu et al., 1974) and had their validity for characterizing the random ionosphere under highly disturbed conditions (Fejer, 1953; Bramley, 1954; Mercier, 1962; Uscinski, 1978; Ishimaru, 1978; Priyadarshi, 2015).

3.2 Power-law models

After the demonstration by Rufenach (1972) that the ionospheric F region irregularities exhibit power-law spectra, multiple measurements have confirmed that this dependence is quite universal for various ionospheric conditions and geographic latitudes (Singleton, 1974; Crane, 1976; Umeki et al., 1977; Franke & Liu, 1983; Bhattacharyya & Rastogi, 1985; Basu et al., 1988a). The three-dimensional power-law PSD in the general form can be written as$${\mathrm{\Phi }}_{\delta n}\left(\mathbf{\kappa }\right)=\frac{f\left(\kappa \right)}{\left(Q\right[\mathbf{\kappa }]+c{)}^{p^{\left(3\right)}/2}},$$(20)where Q[κ] is a quadratic form of spatial frequency, p⁽³⁾ is the three-dimensional spectral index, c is an optional scalar parameter, and f(κ) is some function to assure the convergence of the integrated spectrum (Voitsekhovich, 1995). The quadratic form incorporates the spatial anisotropy of ionospheric irregularities placed in the geomagnetic field, the subject to be addressed in Section 4.

The power law of type (20) signifies that there exists some nonlinear process that allows the inflow of the energy, e.g. due to solar heating or wind shears, to be redistributed between the inhomogeneities. This process is of a cascade type where the energy is redistributed among the different scales while larger inhomogeneities are split into smaller ones and are called the Richardson energy cascade.⁵ An illustration of such a process is the evolution of turbulent flows. In the developed turbulence, the energy-containing inhomogeneities of the largest size L₀, referred to as the outer scale, consist of stretching and bending under the action of inertial forces that results in their splitting into less energetic and smaller inhomogeneities. This process continues until the inner scale ${\mathcal{l}}_0$, i.e., the smallest inhomogeneity scale, is reached, the viscous forces start to dominate the inertial ones, and the energy is drained by dissipation.

Since the spectra (14) and (15) are of reduced dimensionality compared to (20), the corresponding two- and one-dimensional spectral indices would differ from the tree-dimensional one according to$$p^{\left(3\right)}=p^{\left(2\right)}+1=p^{\left(1\right)}+2.$$(21)

It follows that depending on the type of experiment used for obtaining the spectral characteristics of a random medium, the specific spectral index can be deduced, and the other two indices can be calculated by virtue of equation (21). In the following, we adopt the convention of using the one-dimensional spectral index and denote it as p = p⁽¹⁾. This choice is dictated by the convenience in comparison with the vast majority of experimental data and results of a numerical simulation. However, one should keep in mind that only the three-dimensional spectra have a direct physical interpretation as they exhibit a true scaling behavior of the underlying nonlinear cascade process. Table 2 summarizes the typical values of spectral indices for particular ionospheric instability phenomena.

3.2.1 von Karman model

The popular spectral models for phase screen simulation are the von Karman PSDs$$\begin{array}{c}{\mathrm{\Phi }}_{\delta n}\left(\kappa \right)=\frac{\langle \delta n^2\rangle }{{\pi }^{3/2}}\frac{\mathrm{\Gamma }\left(\frac{p+2}{2}\right)}{\mathrm{\Gamma }\left(\frac{p-1}{2}\right)}\frac{{\kappa }_0^{p-1}}{({\kappa }^2+{\kappa }_0^2{)}^{\frac{p+2}{2}}},\\ F_{\delta n}({\kappa }_{\perp },z)=2\frac{\langle \delta n^2\rangle }{\pi \mathrm{\Gamma }\left(\frac{p-1}{2}\right)}{\kappa }_0^{-2}{\left(\frac{{\kappa }_0^{2}z}{2\sqrt{{\kappa }_{\perp }^2+{\kappa }_0^2}}\right)}^{\frac{p+1}{2}}{K}_{\frac{p+1}{2}}\left(z\sqrt{{\kappa }_{\perp }^2+{\kappa }_0^2}\right),\\ W_{\delta n}\left(\kappa \right)=\frac{\langle \delta n^2\rangle }{\sqrt{\pi }}\frac{\mathrm{\Gamma }\left(\frac{p}{2}\right)}{\mathrm{\Gamma }\left(\frac{p-1}{2}\right)}\frac{{\kappa }_0^{p-1}}{({\kappa }^2+{\kappa }_0^2{)}^{\frac{p}{2}}},\end{array}$$(22)where κ₀ = 2π/L₀ is the spatial frequency associated with the outer scale L₀, Γ(x) is the gamma function, and $K_n\left(x\right)$ is the modified Bessel function of the second kind (Macdonald function). The normalization constants in (22) are obtained from equation (19). The corresponding structure functions for the refractive index and phase fluctuations are$$D_{\delta n}\left(r\right)=2\langle \delta n^2\rangle \left[1-\frac{1}{\mathrm{\Gamma }\left(\frac{p-1}{2}\right)}\frac{({\kappa }_{0}r{)}^{\frac{p-1}{2}}}{2^{\frac{p-3}{2}}}{K}_{\frac{p-1}{2}}\left({\kappa }_{0}r\right)\right], \hspace{1em}r\ge 0,$$(23)$$D_{\delta \phi }\left(r_{\perp }\right)=4\langle \delta n^2\rangle \frac{k^2\Delta z}{{\kappa }_0}\frac{\sqrt{\pi }}{\mathrm{\Gamma }\left(\frac{p-1}{2}\right)}\left[\mathrm{\Gamma }\left(\frac{p}{2}\right)-2{\left(\frac{{\kappa }_0{r}_{\perp }}{2}\right)}^{p/2}{K}_{\frac{p}{2}}\left({\kappa }_0{r}_{\perp }\right)\right].$$(24)

Expression (24) is obtained by substituting (22) in equation (17) and integrating.

3.2.2 Shkarofsky model

The von Karman spectra do not depend on the spatial frequency ${\kappa }_m\propto 1/{\mathcal{l}}_0$ associated with the inner scale and, thus, do not account for the change of spectral properties in the high-frequency dissipation region, κ > κ_m. The spectra proposed by Shkarofsky (1968) are useful if this behavior is of importance, see e.g. (Liu & Yeh, 1978; Knepp, 1983; Rino, 2011). The normalized Shkarofsky spectra are$$\begin{array}{c}{\mathrm{\Phi }}_{\delta n}\left(\kappa \right)=\frac{\langle \delta n^2\rangle }{(2\pi {)}^{3/2}}{\left(\frac{{\kappa }_0}{{\kappa }_m}\right)}^{\frac{p-1}{2}}{\kappa }_m^{-3}{t}^{-\frac{p+2}{2}}\frac{K_{\frac{p+2}{2}}\left(t\right)}{K_{\frac{p-1}{2}}\left(\frac{{\kappa }_0}{{\kappa }_m}\right)},\\ F_{\delta n}({\kappa }_{\perp },z)=\frac{\langle \delta n^2\rangle }{2\pi }{\left(\frac{{\kappa }_0}{{\kappa }_m}\right)}^{\frac{p-1}{2}}{\kappa }_m^{-2} {\tilde{t}}^{-\frac{p+1}{2}}\frac{K_{\frac{p+1}{2}}\left(\tilde{t}\right)}{K_{\frac{p-1}{2}}\left(\frac{{\kappa }_0}{{\kappa }_m}\right)},\\ W_{\delta n}\left(\kappa \right)=\frac{\langle \delta n^2\rangle }{(2\pi {)}^{1/2}}{\left(\frac{{\kappa }_0}{{\kappa }_m}\right)}^{\frac{p-1}{2}}{\kappa }_m^{-1}{t}^{-\frac{p}{2}}\frac{K_{\frac{p}{2}}\left(t\right)}{K_{\frac{p-1}{2}}\left(\frac{{\kappa }_0}{{\kappa }_m}\right)},\end{array}$$(25)where$$t=\frac{1}{{\kappa }_m}\sqrt{{\kappa }^2+{\kappa }_0^2}, \tilde{t}=\frac{1}{{\kappa }_m}\sqrt{{\kappa }_{\perp }^2+{\kappa }_0^2}\sqrt{1+{\kappa }_m^2{z}^2}.$$(26)

The corresponding structure functions of refractive index and phase fluctuations read as$$D_{\delta n}\left(r\right)=2\langle \delta n^2\rangle \left[1-(1+{\kappa }_m^2{r}^2{)}^{\frac{p-1}{4}}\frac{K_{\frac{p-1}{2}}\left(\frac{{\kappa }_0}{{\kappa }_m}\sqrt{1+{\kappa }_m^2{r}^2}\right)}{K_{\frac{p-1}{2}}\left(\frac{{\kappa }_0}{{\kappa }_m}\right)}\right],$$(27)$$D_{\delta \phi }\left(r_{\perp }\right)=2\sqrt{2\pi }\langle \delta n^2\rangle \frac{k^2\Delta z}{{\kappa }_m}\sqrt{\frac{{\kappa }_m}{{\kappa }_0}}\frac{K_{\frac{p}{2}}\left(\frac{{\kappa }_0}{{\kappa }_m}\right)}{K_{\frac{p-1}{2}}\left(\frac{{\kappa }_0}{{\kappa }_m}\right)}\left[1-(1+{\kappa }_m^2{r}_{\perp }^2{)}^{\frac{p}{4}}\frac{K_{\frac{p}{2}}\left(\frac{{\kappa }_0}{{\kappa }_m}\sqrt{1+{\kappa }_m^2{r}_{\perp }^2}\right)}{K_{\frac{p}{2}}\left(\frac{{\kappa }_0}{{\kappa }_m}\right)}\right].$$(28)

In deriving equation (28) the relation (17) has been used.

3.2.3 Two-component spectral models

The last spectral model we would like to mention is the two-component spectrum typical for complex systems with two cascade processes involved in forming the refractive index structures. Examples of such behavior are the direct and inverse cascades in the two-dimensional turbulence (Davidson, 2004). The two-component behavior of the spectrum might also be related to the formation of specific geometric forms of irregularities, such as the shocklike structures in equatorial bottom side spread-F (Hysell et al., 1994b). In any case, there exists an additional break at frequency κ_b such that κ₀ < κ_b < κ_m. The two-component spectra have been studied by many authors, see e.g. Refs. Shishov (1974, 1976); Carrano & Rino (2015), and here we derive the normalized spectra following the ideas of Shkarofsky (1968) as$${\mathrm{\Phi }}_{\delta n}\left(\kappa \right)=\frac{1}{{\pi }^{\frac{3}{2}}}\langle \delta n^2\rangle \frac{\mathrm{\Gamma }\left(\frac{p_2+2}{2}\right)}{\mathrm{\Gamma }\left(\frac{p_2-1}{2}\right)}\frac{{\kappa }_0^{p_1-1}{\kappa }_b^{p_2-p_1}{\left({\kappa }^2+{\kappa }_0^2\right)}^{-\frac{p_1+2}{2}}{\left({\kappa }^2+{\kappa }_b^2\right)}^{-\frac{p_2-p_1}{2}}}{{}_2{\mathrm{F}}_1\left(\frac{3}{2},\frac{p_2-p_1}{2},\frac{p_2+2}{2},1-\frac{{\kappa }_0^2}{{\kappa }_b^2}\right)},$$(29)$$F_{\delta n}\left({\kappa }_{\perp },0\right)=\frac{1}{\pi }\langle \delta n^2\rangle \frac{\mathrm{\Gamma }\left(\frac{p_2+1}{2}\right)}{\mathrm{\Gamma }\left(\frac{p_2-1}{2}\right)}\frac{{\kappa }_0^{p_1-1}{\kappa }_b^{p_2-p_1}{\left({\kappa }_{\perp }^2+{\kappa }_0^2\right)}^{-\frac{p_1+1}{2}{\left({\kappa }_{\perp }^2+{\kappa }_b^2\right)}^{-\frac{p_2-p_1}{2}}}}{{}_2{F}_1\left(1,\frac{p_2-p_1}{2};\frac{p_2+1}{2};1-\frac{{\kappa }_0^2}{{\kappa }_b^2}\right)},$$(30)$$W_{\delta n}\left(\kappa \right)=\frac{1}{\sqrt{\pi }}\langle \delta n^2\rangle \frac{\mathrm{\Gamma }\left(\frac{p_2}{2}\right)}{\mathrm{\Gamma }\left(\frac{p_2-1}{2}\right)}\frac{{\kappa }_0^{p_1-1}{\kappa }_b^{p_2-p_1}({\kappa }^2+{\kappa }_0^2{)}^{-\frac{p_1}{2}{\left({\kappa }^2+{\kappa }_b^2\right)}^{-\frac{p_2-p_1}{2}}}}{{}_2{\mathrm{F}}_1\left(1,\frac{p_2-p_1}{2};\frac{p_2}{2};1-\frac{{\kappa }_0^2}{{\kappa }_b^2}\right)}.$$(31)

Here ₂F₁(a, b; c; x) is the hypergeometric function and p₁ and p₂ are two spectral indices that correspond to spatial frequency regions κ₀ < κ ≤ κ_b and κ_b < κ < κ_m, respectively. Table 2 provides some typical spectral indices for ionospheric processes characterized by two-component spectra.

It is worth noting that, in contrast to spectra in equation (25), the spectra of different dimensionality (29), (30), and (31) do not relate to each other by means of relations (14) and (15). Each spectrum is derived by assuming a certain dependence on the spatial frequency with exponents following from the dimensionality reasoning. The one-dimensional spectrum (31) was proposed by Shkarofsky (1968), and the two-dimensional spectrum (30) was used Franke et al. (1984), Franke & Liu (1985) for modeling scintillation in equatorial regions. The expressions for the corresponding structure functions are barely useful for practical calculations, and their expressions will be omitted.

To summarize, each spectral model considered above has its specific advantage from a practical perspective. The Gaussian spectrum is attractive since it allows one to reduce the computation complexity when dealing with the computation of scintillation indices, cf. Ref. Uscinski (1982). Moreover, as has been noted above, this spectrum is applicable to a highly disturbed random ionosphere. The Shkarofsky spectrum has universal applicability in the whole range of spatial frequencies of interest. The two-component spectrum can be considered as the generalization of the Shkarofsky model to the case when the electron density exhibits dependence on two spectral indices. In situations when the two-component spectra relevant for scintillation modeling might occur in structures associated with plasma bubbles (Basu et al., 1983; Muralikrishna et al., 2007; Yokoyama, 2017b), they are also connected with the nonlinear convection (Hysell et al., 1994a), reversed flow events Jin et al. (2019), and polar cap patches Spicher et al. (2015), to name just a few. Finally, the von Karman model is a handy alternative to the Shkarofsky spectrum in the phase screen computer simulations. The corresponding algorithms are based on the discretization of both the spatial domain, where the phase screen is defined, and the spatial frequency domain of the corresponding PSD. The dependence of the spectrum on the inner scale ${\mathcal{l}}_0$, i.e., on the wave frequency κ_m, is not important if the grid spacing of the corresponding phase screen is chosen to be greater than ${\mathcal{l}}_0$. The latter condition is usually met in the practice where the speed of computations is essential.

The discussed Gaussian and power-law spectra and the related structure functions are given in Figure 2. From Figure 2a, it is evident that the von Karman (22), Shkarofsky (25), and two-component spectra (29)–(31) behave very similarly at small spatial frequencies and start to deviate at the spectral break, κ_b, and at the inner scale, κ_m, frequencies. For a two-component spectrum, the change of the spectral slope p₁ in the range κ₀ < κ < κ_b to the slope p₂ in the range κ_b < κ < κ_m is evident. Figure 2b shows the structure functions for phase and refractive index fluctuations for several spectral models. The general behavior of the structure functions is to saturate at spatial separations exceeding the largest irregularity scale given by the outer scale L₀ = 2π/κ₀. The interesting observation is that although the structure functions D_δn behave differently for the models considered, the phase structure functions D_δφ for the Gaussian, von Karman, and Shkarofsky spectra are almost identical. This result is quite surprising considering that the spectral characteristics of the Gaussian spectrum are quite distinct from those of the power-law spectra, cf. Figure 2a. Thus, although the specific spectral components of the Gaussian PSD contribute to fluctuations δn in a distinct way in comparison to the corresponding spectral components of the power-law PSDs, the statistical properties of the resulting phase fluctuations are almost identical. As a consequence, the simulation of the scintillation indices via the random phase screen technique, discussed in Section 6, yields similar results for power-law and Gaussian phase screens, provided that the correlation radius of the Gaussian PSD is set to be equal to the outer scale, i.e. when r₀ = L₀ = 2π/κ₀.

Figure 2 Power spectral densities (a) and the structure functions (b) for several models discussed in the text. The spectral indices are chosen to be: p = p1 = 1.3 and $p_2=3.8$. The one-, two-, and three-dimensional spectra are normalized such that $\langle \delta n^2\rangle =1$. The different power-law models show similar behavior at the outer scale frequency, ${\kappa }_0$, but deviate at the spectral break, ${\kappa }_b$, and the inner scale, ${\kappa }_m$, spatial frequencies. The correlation radius in the Gaussian spectrum is $r_0=2\pi /{\kappa }_0$.

4 Geometric parameters in the modeling of scintillation

The scintillation strength depends on the relative positions of the observer and sender as well as on the geometric form of scintillation-producing irregularities crossing their field of view (Parkin, 1967). Expressed in terms of local coordinates, the first aspect corresponds to the dependence of scintillation indices on the zenith and azimuth angles. Additionally, after scattering within the ionospheric layer, the signal propagates some distance in free space until being detected by the ground observer. This free space propagation should also be accounted for when calculating the scintillation. The second aspect is considered in the modeling of scintillation by incorporating the anisotropy of ionospheric irregularities into the power spectral density of the refractive index fluctuations.

4.1 Dependence on zenith and azimuth angles

To characterize the position of ionospheric irregularities in space, we use the zenith and the azimuth angles in the following considerations. We should differentiate these angles for some spatial point along the communication path where the incoming wave is scattered by irregularity and the corresponding variables of the ground observer. We denote the zenith and azimuth angles for the former case as θ and ϕ, correspondingly, while for the latter case, they are denoted as θ_o and ϕ_o. The line element ds along the slant path at the point with height z above the ground is related to the corresponding vertical length element dz as$$\mathrm{d}s=\sec \theta \mathrm{d}z,$$(32)where θ is the zenith angle with respect to the local vertical direction of the given point. If one neglects the finite curvature of the Earth or the finite curvature of the ionospheric shell, the zenith angle θ would coincide with the observer’s zenith angle θ_o for all points along the straight propagation path.

If the finiteness of the curvature is taken into account, the following relation can be established under the assumption of spherical geometry:$$\mathrm{d}s=\sec \theta \mathrm{d}z=M\left({\theta }_o\right)\mathrm{d}z,\hspace{1em} M\left({\theta }_o\right)=\frac{1}{\sqrt{1-{\left(\frac{R_{\oplus }}{R_{\oplus }+z}\right)}^2{\sin }^2{\theta }_o}}.$$(33)

Here R_⊕ is the Earth radius, z is the altitude of the given point. The distinction between the zenith angles θ and θ_o is important as some authors prefer to use the coordinates associated with the point along the propagation path, referring to equation (32) (Butler, 1951; Ellison & Seddon, 1952; Booker, 1958; Rino, 1979a), while some other authors prefer to perform recalculations according to equation (33) in order to get the scintillation index in terms of the observer’s coordinates (Briggs & Parkin, 1963; Wand & Evans, 1975). The flat geometry approximation θ = θ_o can be used under the condition of small values of zenith angles.

In the approximation of flat geometry, the zenith angle does not depend on the position along the propagation path. Thus, the integration of equation (32) for the slant path yields ${\int }_0^z\mathrm{}\mathrm{d}s=z\sec \theta =z\sec {\theta }_o$. For the spherical geometry, the corresponding path length is obtained from equation (33)$$\underset{0}{\overset{z}{\int }}\mathrm{}\mathrm{d}s=z \cdot S\left({\theta }_o\right)=z \cdot S^{\mathrm{\prime }}\left(\theta \right),$$(34)where the slant path correction factors expressed in terms of the dimensionless variable $\xi =R_{\oplus }/z$ are$$\begin{array}{c}S\left({\theta }_o\right)=\sqrt{1+2\xi +{\xi }^2{\cos }^2{\theta }_o}-\xi \cos {\theta }_o, \hspace{1em}{\theta }_o\in [0,\pi /2]\\ S^{\mathrm{\prime }}\left(\theta \right)=\left(1+\xi \right)\cos \theta -\sqrt{{\left(1+\xi \right)}^2{\cos }^2\theta -1-2\xi },\hspace{1em} \theta \in [0,{\theta }_m],\\ {\theta }_m=\arccos \left[\frac{\sqrt{1+2\xi }}{1+\xi }\right].\end{array}$$(35)

Depending on which coordinates are adopted, either the factor S(θ_o) or the factor S′(θ) should be used in the calculation of the slant path length.

In similar manner we should distinguish the azimuth angle of the observer, ϕ_o, and the corresponding angle at some point along the slant path. The relationship between these two quantities involves the geographic coordinates, namely the latitude Φ and the longitude Λ, of the corresponding points and reads as$$\tan \varphi =\tan {\varphi }_o \frac{\cos {\mathrm{\Phi }}_o\cos (\mathrm{\Lambda }-{\mathrm{\Lambda }}_o)}{\mathrm{cos\Phi }-\tan {\varphi }_0\cos {\mathrm{\Phi }}_o\mathrm{sin\Phi }\sin (\mathrm{\Lambda }-{\mathrm{\Lambda }}_o)}.$$(36)

This equation suggests that in general the azimuth angles ϕ, ϕ_o differ from each other and only in the case of Φ ≈ Φ_o and Λ ≈ Λ_o one can set ϕ ≈ ϕ_o.

The usage of the coordinate sets {θ, ϕ} or {θ_o, ϕ_o} is dictated merely by practical convenience. For example, in the case of the GNSS applications when the signal is recorded by the ground receivers, the appropriate coordinates are θ_o and ϕ_o. For remote sensing applications with SAR satellites, the choice of θ, ϕ is more appropriate for studying the signal scattering in the ionospheric layer. In this case it is also convenient to express these variables in terms of the local coordinates of the transmitting satellite using relations similar to equations (33) and (36).

4.2 Dependence on geometric parameters of irregularities

The dependence of scintillation indices on zenith and azimuth angles is closely related to the geometrical parameters of the anisotropic scintillation-produced irregularities (Singleton, 1970b, 1973). Moreover, as the irregularities are usually field-aligned, the additional rotation on the geomagnetic dip and declination angles should be performed for proper accounting of the relative position of a scintillation-causing irregularity and propagation path.

Some information on the geometrical structure of ionospheric irregularities can give, for example, the classical correlation analysis of diffraction patterns observed by spatially separated ground observers (Briggs, 1968). This analysis reveals that the diffraction patterns could be analyzed if one assumes that the surfaces of constant correlation of refractive index fluctuations are the ellipsoids with the major semi-axes a, b, c, cf. Figure 3a. Following the Singleton convention (Singleton, 1970a), we use the coordinates r, s, t such that the minor ellipsoid semi-axis c is along the t-axis and the semi-axes a and b are along the axes s and r respectively, cf. Figure 3b. Finally, it is convenient to associate the minor semi-axis with the reference correlation radius r₀. Then other semi-axes are expressed through r₀ via the dimensionless scaling factors α, β ≥ 1. Specifically, the equation (Briggs & Parkin, 1963; Singleton, 1970a; Moorcroft & Arima, 1972)$$\frac{t^2}{r_0^2}+\frac{s^2}{(\alpha r_0{)}^2}+\frac{r^2}{(\beta r_0{)}^2}=\mathrm{const}$$(37)defines the surface of constant correlation, i.e. the surface on which ρ_δn and D_δn attain the constant values. Alternatively, the surfaces (37) with α → 1/α, β → 1/β, r₀ → 1/r₀ represent the surfaces in the reciprocal domain where Φ_δn attain constant values.

Figure 3 The ellipsoidal surface with semi-axes a, b, c, is the surface of the constant autocorrelation function of refractive index fluctuations. For expressing the equation of the surface in local geographic coordinates, the transformation from North–East-Down coordinate system x,y,z to the Singleton coordinates $r,s,t$ is used. The transformation consists of: (a) the rotation on the magnetic declination angle ϕ followed by the rotation on the magnetic inclination or dip angle ψ; (b) the rotation on angle $\gamma$ that characterizes the tilt of the ionospheric layer relative to the observer plane. In order to visualize the orientation of the plane within which the geomagnetic field has constant magnitude the multiple field lines are shown in the sketch (b).

In the spatial domain, it is convenient to relate the Singleton coordinates r, s, t to those in the conventional coordinate system x, y, z. If we choose x, y, z to be a NED (North, East, Down) coordinate system, the Singleton and NED coordinates can be related by a series of proper rotations, as shown in Figure 3. The rotation angles are: the magnetic declination δ and inclination (dip) ψ angles as well as the tilt angle γ that characterizes the tilt of a bunch of magnetic field lines relative to the observation plane. These geometric considerations result in the following definition of the quadratic form of spatial frequencies that enters the general definition of three-dimensional PSD in equation (20) (Rino & Fremouw, 1977; Rino, 1979a, 1982a, 2011):$$Q\left[\mathbf{\kappa }\right]={\mathbf{\kappa }}^T\widehat{\bar{C}}\mathbf{\kappa }.$$(38)

Here $\widehat{\bar{C}}$ is a symmetrical matrix with the following elements:$$\begin{array}{c}{\widehat{\bar{C}}}_{11}=r_0^2[{\alpha }^2{\cos }^2\psi {\cos }^2\delta +{\sin }^2\psi {\cos }^2\delta \left({\cos }^2\gamma +{\beta }^2{\sin }^2\gamma \right)\\ +{\sin }^2\delta \left({\sin }^2\gamma +{\beta }^2{\cos }^2\gamma \right)+\frac{1}{2}\sin 2\gamma \sin \psi \sin 2\delta \left({\beta }^2-1\right)],\end{array}$$(39)$$\begin{array}{c}{\widehat{\bar{C}}}_{22}=r_0^2[{\alpha }^2{\cos }^2\psi {\sin }^2\delta +{\sin }^2\psi {\sin }^2\delta \left({\cos }^2\gamma +{\beta }^2{\sin }^2\gamma \right)\\ +{\cos }^2\delta \left({\sin }^2\gamma +{\beta }^2{\cos }^2\gamma \right)-\frac{1}{2}\sin 2\gamma \sin \psi \sin 2\delta \left({\beta }^2-1\right)],\end{array}$$(40)$${\widehat{\bar{C}}}_{33}=r_0^2\left[{\alpha }^2{\sin }^2\psi +{\cos }^2\psi \left({\cos }^2\gamma +{\beta }^2{\sin }^2\gamma \right)\right],$$(41)$$\begin{array}{c}{\widehat{\bar{C}}}_{12}={\widehat{\bar{C}}}_{21}=\frac{r_0^2}{2}\sin 2\delta \left[{\sin }^2\gamma +{\beta }^2{\cos }^2\gamma -{\alpha }^2{\cos }^2\psi -{\sin }^2\psi \left({\cos }^2\gamma +{\beta }^2{\sin }^2\gamma \right)\right]\\ +\frac{r_0^2}{2}\cos 2\delta \sin 2\gamma \sin \psi \left({\beta }^2-1\right),\end{array}$$(42)$${\widehat{\bar{C}}}_{13}={\widehat{\bar{C}}}_{31}=\frac{r_0^2}{2}\sin 2\psi \cos \delta \left({\alpha }^2-{\cos }^2\gamma -{\beta }^2{\sin }^2\gamma \right)-\frac{r_0^2}{2}\sin 2\gamma \cos \psi \sin \delta \left({\beta }^2-1\right),$$(43)$${\widehat{\bar{C}}}_{23}={\widehat{\bar{C}}}_{32}=-\frac{r_0^2}{2}\sin 2\psi \sin \delta \left({\alpha }^2-{\cos }^2\gamma -{\beta }^2{\sin }^2\gamma \right)-\frac{r_0^2}{2}\sin 2\gamma \cos \psi \cos \delta \left({\beta }^2-1\right).$$(44)

In the limiting case δ = 0 these formulas coincide with the results of Rino & Fremouw (1977).

If the statistical properties of phase fluctuations are of interest, for example, in connection with the phase screen modeling, then, in the view of definition (17), the quadratic form to be inserted in the spectrum ${\mathrm{\Phi }}_{\delta n}$ is$$Q[{\mathbf{\kappa }}_{\perp },{\kappa }_z=0]={\widehat{\bar{C}}}_{11}{\kappa }_x^2+2{\widehat{\bar{C}}}_{12}{\kappa }_x{\kappa }_y+{\widehat{\bar{C}}}_{22}{\kappa }_y^2.$$(45)

This result is valid if the propagation direction coincides with the z axis of the NED coordinate sysem, i.e., for zero zenith angle. For small zenith angles, the flat geometry approximation can be used, i.e., one can assume that the ionospheric layer is plan-parallel and the Earth’s surface is approximately flat. In this case, the corrections to equation (45) due to the finiteness of the zenith angle can be incorporated by substituting ${\kappa }_z=-\tan \theta {\mathbf{a}}_{k,\perp }\cdot {\mathbf{\kappa }}_{\perp }$ as has been shown in Rino & Fremouw (1977). The quadratic form can be written as$$Q[{\mathbf{\kappa }}_{\perp },{\kappa }_z=-\tan \theta {\mathbf{a}}_{k,\perp }\cdot {\mathbf{\kappa }}_{\perp }]=\mathcal{A}{\kappa }_x^2+2\mathcal{B}{\kappa }_x{\kappa }_y+\mathcal{C}{\kappa }_y^2,$$(46)where$$\begin{array}{c}\mathcal{A}={\widehat{\bar{C}}}_{11}+{\widehat{\bar{C}}}_{33}{\tan }^2\theta {\cos }^2\varphi -2{\widehat{\bar{C}}}_{13}\tan \theta \cos \varphi ,\\ \mathcal{B}={\widehat{\bar{C}}}_{12}+{\widehat{\bar{C}}}_{33}{\tan }^2\theta \sin \varphi \cos \varphi -\tan \theta \left({\widehat{\bar{C}}}_{13}\sin \varphi +{\widehat{\bar{C}}}_{23}\cos \varphi \right),\\ \mathcal{C}={\widehat{\bar{C}}}_{22}+{\widehat{\bar{C}}}_{33}{\tan }^2\theta {\sin }^2\varphi -2{\widehat{\bar{C}}}_{23}\tan \theta \sin \varphi \end{array}$$(47)are the coefficients of the quadratic form expressed in terms of the observer azimuth ϕ and zenith θ as long as the approximation (32) is applicable. The vector ${\mathbf{a}}_{k,\perp }=(\cos \varphi \sin \varphi {)}^T$ is the unit vector along the transverse direction to the wave vector k such that ${\mathbf{a}}_{k,\perp }\cdot {\mathbf{\kappa }}_{\perp }={\kappa }_x\cos \varphi +{\kappa }_y\sin \varphi$. The relevant geometrical aspects of the model (46) are shown in Figure 4a. We also note that the assumption of the flat geometry is equivalent to the assumption that θ = θ_o as discussed in Section 4.1. Some improvement of this approximation could serve the hybrid model, in which one uses the quadratic form (46) but defines the dependence between the zenith angles of the observer and the point along the slant path via the relation (36). Such modification allows one to avoid the divergencies associated with secθ at large zenith angles.

Figure 4 Cross-sections (hatched areas) of ellipsoids of constant autocorrelation of electron density fluctuations with the plane associated either with the slab of irregularities (a) or with the wave propagation direction (b). Here the local coordinate system x,y,z (North, East, Down) originates at height z above the observation plane XY. The coordinate system x′, y′, z′ is chosen in such a way that the z′ axis is pointed along the wave vector k while the x′ and z′ axes lie in the plane formed by k and the down axis z. In turn, the direction of the wave vector is defined in terms of the zenith angle θ and the azimuth angle ϕ. The models (a) and (b) also assume the flat and spherical geometries, respectively. This results in different definitions of the slant ranges and transformation rules for θ, ϕ into the observer’s zenith, θo, and azimuth ϕo, angles (see the text for more details).

Finally, the following quadratic form can be derived for an arbitrary value of the zenith angle:$$Q[{\mathbf{\kappa }}_{\perp }^{\mathrm{\prime }},{\kappa }_z^{\mathrm{\prime }}=0]={\mathcal{A}}^{\mathrm{\prime }}{{\kappa }_x^{\mathrm{\prime }}}^2+2{\mathcal{B}}^{\mathrm{\prime }}{\kappa }_x^{\mathrm{\prime }}{\kappa }_y^{\mathrm{\prime }}+{\mathcal{C}}^{\mathrm{\prime }}{{\kappa }_y^{\mathrm{\prime }}}^2,$$(48)where$$\begin{array}{c}{\mathcal{A}}^{\mathrm{\prime }}={\widehat{\bar{C}}}_{11}{\cos }^2\theta {\cos }^2\varphi +{\widehat{\bar{C}}}_{22}{\cos }^2\theta {\sin }^2\varphi +{\widehat{\bar{C}}}_{33}{\sin }^2\theta +{\widehat{\bar{C}}}_{12}{\cos }^2\theta \sin 2\varphi -{\widehat{\bar{C}}}_{13}\sin 2\theta \cos \varphi -{\widehat{\bar{C}}}_{23}\sin 2\theta \sin \varphi ,\\ {\mathcal{B}}^{\mathrm{\prime }}=\frac{1}{2}\left({\widehat{\bar{C}}}_{22}-{\widehat{\bar{C}}}_{11}\right)\cos \theta \sin 2\varphi +{\widehat{\bar{C}}}_{12}\cos \theta \cos 2\varphi +{\widehat{\bar{C}}}_{13}\sin \theta \sin \varphi -{\widehat{\bar{C}}}_{23}\sin \theta \cos \varphi ,\\ {\mathcal{C}}^{\mathrm{\prime }}={\widehat{\bar{C}}}_{11}{\sin }^2\varphi +{\widehat{\bar{C}}}_{22}{\cos }^2\varphi -{\widehat{\bar{C}}}_{12}\sin 2\varphi .\end{array}$$(49)

Here the spatial frequencies ${\kappa }_x^{\mathrm{\prime }}$, ${\kappa }_y^{\mathrm{\prime }}$ are reciprocal coordinates of the coordinates $x^{\mathrm{\prime }}$, $y^{\mathrm{\prime }}$ defined in the plane transversal to the propagation direction, i.e., to the wave vector k. The model (49) assumes the spherical symmetry of the Earth and ionospheric shell. This assumption, to which we refer as the spherical geometry approximation, is not biased by the projection transformation of the kind used in (46), and, thus, is consistent for arbitrary values of zenith angles.

4.3 Propagation factors

So far, we have considered the propagation through a random ionospheric medium of thickness Δz, cf. equation (3). In practical situations, there is a finite distance from the layer of such medium and the observation plane placed at z > Δz in the coordinate system chosen, cf. Figures 1a and 1b. The propagation over such a distance can be considered to happen in a vacuum. The propagation factors take this effect into account (Clifford, 1978; Tatarski, 2016).

Here we consider two models corresponding to equations (46) and (48) with the note that the obtained results for the zero zenith angle are the same as obtained within the model (45). We also consider a weak scattering case that corresponds to placing a single phase screen at z = 0. For the spectrum based on the model (46) the correlation functions for the phase departure, S, and for its log-amplitude, χ = log(1 + δE/〈E〉), can be obtained as$$\begin{array}{c}\left\{\begin{array}{l}{\rho }_S\left({\mathbf{R}}_{\perp }\right)\\ {\rho }_{\chi }\left({\mathbf{R}}_{\perp }\right)\end{array}\right\}=k^{2}z{\sec }^2\theta \underset{0}{\overset{z}{\int }}\mathrm{} \int \mathrm{}{\mathrm{\Phi }}_{\delta n}\left(\mathbf{\kappa }\right)\left\{\begin{array}{l}{\mathcal{P}}_S\left[\frac{{\kappa }^{2}z\sec \theta }{k}\right]\\ {\mathcal{P}}_{\chi }\left[\frac{{\kappa }^{2}z\sec \theta }{k}\right]\end{array}\right\}\cos ({\mathbf{\kappa }}_{\perp }\cdot {\mathbf{R}}_{\perp })\cos \left[\right({\kappa }_z+\tan \theta {\mathbf{a}}_{k,\perp }\cdot {\mathbf{\kappa }}_{\perp }\left)z^{\mathrm{\prime }}\right]\mathrm{d}{z}^{\mathrm{\prime }} {\mathrm{d}}^3\mathbf{\kappa }\\ =2\pi k^{2}z{\sec }^2\theta \int \mathrm{}{\mathrm{\Phi }}_{\delta n}({\mathbf{\kappa }}_{\perp },{\kappa }_z=-\tan \theta {\mathbf{a}}_{k,\perp }\cdot {\mathbf{\kappa }}_{\perp })\left\{\begin{array}{l}{\mathcal{P}}_S\left[2h\right({\mathbf{\kappa }}_{\perp }\left)Z\right]\\ {\mathcal{P}}_{\chi }\left[2h\right({\mathbf{\kappa }}_{\perp }\left)Z\right]\end{array}\right\}\cos ({\mathbf{\kappa }}_{\perp }\cdot {\mathbf{R}}_{\perp }){\mathrm{d}}^2{\mathbf{\kappa }}_{\perp },\end{array}$$(50)where R_⊥ = (X Y)^T is the vector defined in the observer plane such that ${\mathbf{R}}_{\perp }={\mathbf{r}}_{\perp }-z \tan \theta {\mathbf{a}}_{k,\perp }$, cf. Figure 4 and$$h\left({\mathbf{\kappa }}_{\perp }\right)={\kappa }_{\perp }^2+{\tan }^2\theta ({\mathbf{a}}_{k,\perp }\cdot {\mathbf{\kappa }}_{\perp }{)}^2, Z=\frac{z}{2k}\sec \theta .$$(51)

In a similar way for the spectral model based on equation (48) we obtain$$\begin{array}{c}\left\{\begin{array}{l}{\rho }_S\left({\mathbf{R}}_{\perp }\right)\\ {\rho }_{\chi }\left({\mathbf{R}}_{\perp }\right)\end{array}\right\}=2\pi k^2 z \left[S^{\mathrm{\prime }}\right(\theta ){]}^2 \int \mathrm{}{\mathrm{\Phi }}_{\delta n}({\mathbf{\kappa }}_{\perp }^{\mathrm{\prime }},0\left)\left\{\begin{array}{l}{\mathcal{P}}_S\left[2\right({\mathbf{\kappa }}_{\perp }^{\mathrm{\prime }}{)}^2{Z}^{\mathrm{\prime }}]\\ {\mathcal{P}}_{\chi }\left[2\right({\mathbf{\kappa }}_{\perp }^{\mathrm{\prime }}{)}^2{Z}^{\mathrm{\prime }}]\end{array}\right\}\cos \right({\mathbf{\kappa }}_{\perp }^{\mathrm{\prime }}\cdot {\mathbf{R}}_{\perp }^{\mathrm{\prime }}){\mathrm{d}}^2{\mathbf{\kappa }}_{\perp }^{\mathrm{\prime }},\\ Z^{\mathrm{\prime }}=\frac{z}{2k}{S}^{\mathrm{\prime }}\left(\theta \right) {\mathbf{R}}_{\perp }^{\mathrm{\prime }}=\widehat{\mathbf{T}}{\mathbf{R}}_{\perp },\end{array}$$(52)where $S^{\mathrm{\prime }}\left(\theta \right)$ is given by equation (35) and the variable ${\mathbf{\kappa }}_{\perp }^{\mathrm{\prime }}$ is defined in the plane transversal to the propagation direction while ${\mathbf{R}}_{\perp }$ is transversal to the observer’s vertical direction, see Figure 4. Here $\widehat{\mathbf{T}}$ is the transformation matrix from ${\mathbf{R}}_{\perp }$ to the coordinates that are transversal to the propagation direction and have the origin at the ground. The explicit form of $\widehat{\mathbf{T}}$ is not relevant for further discussions.