© G. Morel et al., Published by EDP Sciences 2025

1 Introduction

This paper is the first part of a study on ionospheric scintillation modeling and analysis. It addresses the analysis of scintillation indices. The second part focuses on power spectral density analysis and is the topic of a companion paper.

The ionospheric scintillation phenomenon is a rapid fluctuation of the amplitude and phase of an electromagnetic wave propagated through the ionosphere. It is related to small-scale turbulent fluctuations of the ionospheric refractive index that are likely to develop (Tatarskii, 1971; Ishimaru, 1997). Despite advances in technology and signal processing, the impact of ionospheric scintillation is still a problem for many radiofrequency systems operating below 10 GHz. Variations in amplitude, linked to constructive or destructive interferences of the propagated wave, are frequently observed at the equator during equinoxes and after sunset. Phase variations, linked to electron density disturbances associated with significant velocities of the ionospheric medium, are observed in polar areas during periods of geomagnetic storms and intense solar activity. These amplitude and phase scintillations create GNSS receiver loss of lock or positioning errors that affect geopositioning system performances at L band (Paziewski et al., 2022). Scintillation effects can also affect low-frequency, L-band and P-band, Synthetic Aperture Radar (SAR) systems (Belcher, 2008; Rogers et al., 2014) or telecommunications.

Since the ionosphere can be a highly fluctuating medium, its effects on propagation are generally statistically quantified, mainly in terms of signal log-amplitude and phase variability (Wheelon, 2001). To model electromagnetic wave propagation, two approaches are commonly used: The Rytov’s theory of smooth perturbations, allowing derivations of analytic and asymptotic formulations of the electromagnetic field moments (Tatarskii, 1971; Rytov et al., 1989; Wheelon, 2003), and a numerical approach, covering all scattering regimes and based on the Parabolic Wave Equation (PWE) method (Tappert, 1977) that can be solved in a stochastic medium (Rino, 2011). The PWE method is then combined with Multiple Phase Screens (MPS) to take into account the turbulent characteristics of the random propagation medium. PWE-MPS method is a numerical approach that has been widely used in the literature to assess the impact of atmospheric turbulence on wave propagation (Knepp, 1983; Fabbro & Feral, 2012). The PWE-MPS method can be applied in three dimensions (i.e., 3D PWE associated with 2D Multiple Phase Screens: 3D-PWE/2D-MPS) or reduced to two dimensions in certain geometrical conditions (i.e., 2D PWE associated with 1D Multiple Phase Screens: 2D-PWE/1D-MPS), taking advantage of the anisotropy property of the ionosphere irregularities (Galiègue et al., 2016). Those methods are often used by considering an electromagnetic plane wave illuminating the ionosphere layer, which is well suited for GNSS applications, i.e., when the satellite is far from the ionosphere, but becomes questionable for lower orbit transmitters. Therefore, in order to extend applications of PWE-MPS methods, the use of an incident spherical wave has been proposed by several authors such as Grimault (1998) or Knepp (2016).

The originality of this paper, by comparison to “Fabbro & Feral (2012), Galiègue (2015), Galiègue et al. (2016)” is to assess the impact of the incident plane wave approximation in ionospheric scintillation modeling by comparing it with the more general electromagnetic incident spherical wave approach. To do so, Rytov’s theory and 3D-PWE/2D-MPS numerical approach are adapted for an incident spherical wave illuminating the ionosphere. To validate the presented results, both methods are compared. The analysis is carried out in terms of amplitude and phase scintillation indices.

The paper is organized as follows: Some reminders on the ionospheric inhomogeneities modeling are presented in Section 2. Section 3 focuses on analytic formulations of the ionospheric scintillation impact in terms of: Log-amplitude and phase variances for both incident Plane Wave (PW) and Spherical Wave (SW). In addition to these formulations, a heuristic approach using a modification of the plane wave formalism, called Corrected Plane Wave (CPW), is proposed. In Section 4, the PWE-MCPS numerical method considering an incident spherical wave is described to be compared in Section 5 with the other approaches presented in this paper, i.e., analytic incident PW/SW/CPW formulations. In Section 6, real system configurations are considered, and the impact of the sphericity of the wave is studied for LEO and MEO satellites. Finally, the different results obtained are summarized and discussed in Section 7, which concludes the paper.

2 Statistical description of ionospheric irregularities

In the ionosphere and for a signal propagated at a frequency, much larger than the plasma frequency, which is usually in the range of a few MHz (Hapgood, 2018), i.e., for microwave signals, the refractive index n is a function of the position r [m] and time t [s]. It is related to the electron density N_e(r, t) by:

$$n^2(\mathbf{r},t)\mathrm{ }=1-\frac{e^2}{{ϵ}_{0}m{\omega }^2}{N}_e(\mathbf{r},t),$$(1)

where e [A s], m [kg] are the charge and the mass of the electron, ϵ₀ [F m⁻¹] the permittivity in free space and ω [rad s⁻¹] the signal pulsation. The refractive index is usually decomposed into a stable component 〈n₀〉 that describes the large-scale refraction effects and a random turbulent component Δn(r, t):

$$n(\mathbf{r},t)=\langle n_0\rangle (1+\Delta n(\mathbf{r},t\left)\right).$$(2)

Attention being focused on the random part of the refractive index, its stable component can be written, without loss of generality, as 〈n₀〉 = 1. Comparing the first order series expansion of n² (r, t) with its previous expression (1), it comes:

$$\Delta n(\mathbf{r},t)=-\frac{e^2}{2{ϵ}_{0}m{\omega }^2}\Delta N_e(\mathbf{r},t)=-\frac{r_e\lambda }{k_0}\Delta N_e(\mathbf{r},t),$$(3)

where r_e = e²/(4πϵ₀mc²) [m] is the electronic radius, λ [m] the signal wavelength, and k₀ [m⁻¹] the free space wave number and ΔN_e (r, t) the electronic density fluctuations.

It is assumed that the turbulent random medium is statistically homogeneous (Wheelon, 2001), with mean 〈ΔN_e (r, t)〉 = 0 and variance ${\sigma }_{∆N_e}^2$ [m⁻⁶]. In such conditions, it follows that the 3-D spatial covariance of electronic density fluctuations $B_{\Delta N_e}(\mathbf{r},\mathbf{r}\mathit{\text{'}})=\langle \Delta N_e\left(\mathbf{r}\right)\Delta N_e\left(\mathbf{r}\mathit{\text{'}}\right)\rangle$ reduces to $B_{\Delta N_e}(\mathbf{r}-\mathbf{r}\mathit{\text{'}})$. Then by applying the Wiener–Khintchine theorem, the following relationship between the 3-D spatial covariance and the 3-D spatial spectrum $S_{\Delta N_e}\left(\mathbf{k}\right)$ follows (Fabbro & Feral, 2012):

$$B_{\Delta N_e}(\mathbf{r}-\mathbf{r}\mathit{\text{'}})={\mathrm{F}}^{-1}\left\{S_{\Delta N_e}\right(\mathbf{k}\left)\right\}={\int }_{-\infty }^{+\infty }\mathrm{}{S}_{\Delta N_e}\left(\mathbf{k}\right)e^{j\mathbf{k}.(\mathbf{r}-\mathbf{r}\mathit{\text{'}})}{\mathrm{d}}^3\mathbf{k},$$(4)

$$S_{\Delta N_e}\left(\mathbf{k}\right)=\mathrm{F}\left\{B_{\Delta N_e}\right(\mathbf{r}\left)\right\}=\frac{1}{(2\pi {)}^3}{\int }_{-\infty }^{+\infty }\mathrm{}{B}_{\Delta N_e}\left(\mathbf{r}\right)e^{-j\mathbf{k}.\mathbf{r}}{\mathrm{d}}^3\mathbf{r},$$(5)

where F{} and F⁻¹ {} are the direct and inverse Fourier transforms. The spectrum of electronic density fluctuations $S_{\Delta N_e}\left(\mathbf{k}\right)$ has now to be defined. Observations show that plasma irregularities are elongated and aligned along the terrestrial magnetic field (Aarons, 1982). Thus, the irregularities are commonly modeled with an ellipsoid whose dimensionless scaling factors are a_x ≤ a_y ≤ a_z, with z the axis chosen along the Earth’s magnetic field vector, (x, y, z) forming a direct Cartesian coordinate system (cf. Fig. 1). By defining anisotropic ratios:

$$A_x=\frac{a_x}{a_x}=1,$$(6)

$$A_y=\frac{a_y}{a_x}\ge 1,$$(7)

$$A_z=\frac{a_z}{a_x}\ge 1,$$(8)

the spectrum can be expressed in the dual magnetic basis (k_x, k_y, k_z) according to (Galiègue et al., 2016):

$$S_{\Delta N_e}\left(\mathbf{k}\right)=a_x^{3-p_m}{A}_y{A}_z{C}_s{\left(k_x^2+A_y^2{k}_y^2+A_z^2{k}_z^2+\frac{K_0^2}{a_x^2}\right)}^{-p_m/2},$$(9)

where $k_0=\frac{2\pi }{L_0}$ [m⁻¹] is the wave number associated with the outer scale of the turbulence, p_m the slope of the power spectral density of the irregularities, and C_s [m^−3−pm] the strength of the ionospheric turbulent structure (Rino, 1979; Galiègue et al., 2016). Then, it is possible to express the spectrum of electronic density fluctuations in the Line Of Sight (LOS) coordinate system (u, v, s) as illustrated in Figure 1, where γ is the angle between the line of sight direction s and the z axis that coincides with the local terrestrial magnetic field H₀. The orthogonal projection of H₀ in the plane (u, v) orthogonal to the line of sight s makes an angle α_z with the direction (Ov). ψ is the angle between (Oy) and its orthogonal projection in the plane (u, v). In this new coordinate system, the spectrum of irregularities becomes (Galiègue et al., 2016):

$$S_{\Delta N_e}(k_u,k_v,k_s)=a_x^{3-p_m}{A}_y{A}_z{C}_s{\left(A{k}_u^2+B{k}_v^2+2C{k}_u{k}_v+f\left(k_s\right)+\frac{K_0^2}{a_x^2}\right)}^{-p_m/2},$$(10)

where A, B and C are the anisotropic coefficients expressed in the line of sight basis, depending on the geometry (i.e., γ, α_z, and ψ angles, see Fig. 1) and anisotropic ratios (A_y, A_z), their expressions are given in Appendix A. The development of f(k_s), which depends on the geometry, is not given here because in the next sections k_s = 0 and f(k_s) = 0.

Figure 1 Geometry of the ionospheric turbulent irregularities in the LOS coordinate system (u, v, s).

Substituting expressions (9) or (10) in the covariance expression (4) and by setting r = 0, one can derive the variance of electron density fluctuations:

$${\sigma }_{\Delta N_e}^2=\frac{{\pi }^{3/2}\mathrm{\Gamma }\left(\frac{p_m-3}{2}\right)}{K_0^{p_m-3}\mathrm{\Gamma }\left(\frac{p_m}{2}\right)}{C}_s.$$(11)

3 3D analytic formulation under weak scattering assumption

In the previous section, ionospheric irregularities have been described by their power spectral density. From this description, based on the Rytov assumption, analytic formulations of variances of log-amplitude and phase are presented. The Rytov method is an approach that allows analytic expressions of log-amplitude and phase variances, scintillation indices, and spectra to be derived under several assumptions, which are: Neglecting the polarization effects, neglecting the induced current density at the transmitter, and considering weak perturbations (Rytov et al., 1989; Wheelon, 2003). Under those assumptions, the Helmholtz propagation equation reduces to:

$${\nabla }^{2}E\left(\mathbf{r}\right)+k_0^2(1+2\Delta n(\mathbf{r}\left)\right)E\left(\mathbf{r}\right)=0,$$(12)

with k₀ the free space wave number of the transmitted wave. Equation (12) can be solved by applying the Rytov approximation, which is valid under weak scattering regime, i.e., for a variance 〈χ²〉 of the log-amplitude χ lower than one (Wheelon, 2003). The idea is to express the solution as the product of two terms, the first one is the electric field E₀(r) which would be measured without ionospheric irregularities, and the second one is the substitution function e^ψ(r): E(r) = E₀(r)e^ψ(r). The function ψ can be expressed as a series of powers of the refractive index variations: $\psi \left(\mathbf{r}\right)={\sum }_{m=1}^{+\infty }\mathrm{}{\psi }_m\left(\mathbf{r}\right)$, with ψ_m ∝ (Δn)^m. By considering just the first order of this series development, and calculating for r = R, which corresponds to the position vector of the receiver, it comes:

$$E_1\left(\mathbf{R}\right)=E_0\left(\mathbf{R}\right){\mathrm{e}}^{{\psi }_1\left(\mathbf{R}\right)},$$(13)

$${\psi }_1\left(\mathbf{R}\right)=-2k_0^2{\int }_{-\infty }^{+\infty }\mathrm{}G(\mathbf{R},\mathbf{r})\Delta n(\mathbf{r},t)\frac{E_0\left(\mathbf{r}\right)}{E_0\left(\mathbf{R}\right)}{\mathrm{d}}^3\mathbf{r}.$$(14)

So, the received electric field is expressed in terms of amplitude $\mathcal{A}\left(\mathbf{R}\right)$ and phase φ(R) such as E₁ $\left(\mathbf{R}\right)=\mathcal{A}\left(\mathbf{R}\right){\mathrm{e}}^{j({\phi }_0+\phi (\mathbf{R}\left)\right)}$, and by definition, the log-amplitude is $\chi =\log \frac{\mathcal{A}}{E_0}$, and finally (Rytov et al., 1989; Wheelon, 2003):

$$\phi \left(\mathbf{R}\right)=-2k_0^2{\int }_{-\infty }^{+\infty }\mathrm{}\Delta n\left(\mathbf{r}\right)\mathfrak{I}\left(G(\mathbf{R},\mathbf{r})\frac{E_0\left(\mathbf{r}\right)}{E_0\left(\mathbf{R}\right)}\right){\mathrm{d}}^3\mathbf{r},$$(15)

$$\chi \left(\mathbf{R}\right)=-2k_0^2{\int }_{-\infty }^{+\infty }\mathrm{}\Delta n\left(\mathbf{r}\right)\mathfrak{R}\left(G(\mathbf{R},\mathbf{r})\frac{E_0\left(\mathbf{r}\right)}{E_0\left(\mathbf{R}\right)}\right){\mathrm{d}}^3\mathbf{r},$$(16)

with $\mathfrak{R}\left(\right)$ and $\mathfrak{J}\left(\right)$, the real and imaginary parts and G(R, r) the Green function. The chosen geometry is described in Figure 2 where ϑ is the angle between the local vertical and the line of sight. The (u, v) plane is the plane transverse to the propagation direction. L_t is the distance traveled by the signal between the transmitter and the top of the ionospheric layer, R_iono is the slant distance traveled by the signal inside the ionospheric layer, and L_v is the distance traveled by the signal between the bottom of the ionospheric layer and the receiver. The total distance is R = L_t + R_iono + L_v. H is the ionospheric layer altitude and ΔH is its thickness. With some trigonometric calculations, it can be shown that: L_v only depends on the Earth radius R_T, ϑ, and H; whereas R_iono depends on R_T, ϑ, H, and ΔH. Indeed, using the law of cosines, expressions for L_v, R_iono, and L_t are the following:

$$L_v=\sqrt{R_T^2{\cos }^2\left(\vartheta \right)+H^2+2R_{T}H}-R_T\cos \left(\vartheta \right),$$(17)

$$R_{\mathrm{iono}}=\sqrt{(R_T+H+\Delta H{)}^2-2R_T{\sin }^2(\vartheta )}-\sqrt{R_T^2{\cos }^2\left(\vartheta \right)+H^2+2R_{T}H},$$(18)

$$L_t=\sqrt{R_T^2{\cos }^2\left(\vartheta \right)+H_{\mathrm{sat}}^2+2R_T{H}_{\mathrm{sat}}}-\sqrt{(R_T+H+\Delta H{)}^2-2R_T{\sin }^2(\vartheta )},$$(19)

we note that only L_t depends on the satellite altitude H_sat.

Figure 2 Description of the satellite-ground link in the Line of Sight (LOS) coordinatesystem (u, v, s). R = Lv + Riono + Lv.

3.1 3D log-amplitude and phase variances

3.1.1 General expression

Starting from equations (15) and (16), the expressions for log-amplitude and phase variances can be derived as (Fabbro & Feral, 2012; Galiègue et al., 2016):

$$\langle {\chi }^2\rangle =\pi r_e^2{\lambda }^2{R}_{\mathrm{iono}}{\int }_{-\infty }^{+\infty }\mathrm{}{\int }_{-\infty }^{+\infty }\mathrm{}{S}_{\Delta N_e}(0,k_u,k_v){\mathcal{F}}_{\chi }(k_u,k_v)\mathrm{d}{k}_u\mathrm{d}{k}_v,$$(20)

$$\langle {\phi }^2\rangle =\pi r_e^2{\lambda }^2{R}_{\mathrm{iono}}{\int }_{-\infty }^{+\infty }\mathrm{}{\int }_{-\infty }^{+\infty }\mathrm{}{S}_{\Delta N_e}(0,k_u,k_v){\mathcal{F}}_{\phi }(k_u,k_v)\mathrm{d}{k}_u\mathrm{d}{k}_v,$$(21)

where λ [m] is the wavelength of the transmitted wave, R_iono [m] is the inhomogeneities thickness, F_χ and F_φ are the spectral weighting functions for amplitude and phase (also called filtering functions). Depending on the nature of the incident wave, plane or spherical, the formulation of these spectral weighting functions is different. This point is detailed in the following sections.

3.1.2 Incident plane wave approach

For an incident plane wave illuminating the ionosphere, the spectral weighting functions for log-amplitude and phase, ${\mathcal{F}}_{\chi }^{\mathrm{PW}}(k_u,k_v)$ and ${\mathcal{F}}_{\phi }^{\mathrm{PW}}(k_u,k_v)$ respectively, can be expressed as the following (Wheelon, 2003):

$${\mathcal{F}}_{\chi }^{\mathrm{PW}}(k_u,k_v)=\frac{2}{R_{\mathrm{iono}}}{\int }_{L_v}^{L_v+R_{\mathrm{iono}}}\mathrm{}{\sin }^2\left(\frac{s{k}^2}{2k_0}\right)\mathrm{d}s,$$(22)

$${\mathcal{F}}_{\phi }^{\mathrm{PW}}(k_u,k_v)=\frac{2}{R_{\mathrm{iono}}}{\int }_{L_v}^{L_v+R_{\mathrm{iono}}}\mathrm{}{\cos }^2\left(\frac{s{k}^2}{2k_0}\right)\mathrm{d}s,$$(23)

where $k^2=k_u^2+k_v^2$. We note that, with the change of variable s = R − s′, these expressions are equivalent to the following:

$${\mathcal{F}}_{\chi }^{\mathrm{PW}}(k_u,k_v)=\frac{2}{R_{\mathrm{iono}}}{\int }_{L_t}^{L_t+R_{\mathrm{iono}}}\mathrm{}{\sin }^2\left(\frac{(R-s\text{'})k^2}{2k_0}\right)\mathrm{d}s\text{'},$$(24)

$${\mathcal{F}}_{\phi }^{\mathrm{PW}}(k_u,k_v)=\frac{2}{R_{\mathrm{iono}}}{\int }_{L_t}^{L_t+R_{\mathrm{iono}}}\mathrm{}{\cos }^2\left(\frac{(R-s\text{'})k^2}{2k_0}\right)\mathrm{d}s\text{'}.$$(25)

Then, after some mathematical developments and by setting:

$$X=\frac{(k_u^2+k_v^2)L_v}{k_0},$$(26)

$$\varsigma =\frac{R_{\mathrm{iono}}}{2L_v},$$(27)

it comes (Galiègue et al., 2016):

$${\mathcal{F}}_{\chi }^{\mathrm{PW}}(k_u,k_v)=1-\mathrm{sinc}\left(\mathrm{X\varsigma }\right)\cos \left(X\right(1+\varsigma \left)\right),$$(28)

$${\mathcal{F}}_{\phi }^{\mathrm{PW}}(k_u,k_v)=1+\mathrm{sinc}\left(\mathrm{X\varsigma }\right)\cos \left(X\right(1+\varsigma \left)\right).$$(29)

Asymptotic formulations can be derived for the log-amplitude and phase variances under several assumptions: A thin turbulent layer compared to the propagation distance (R_iono ≪ L_v), and a high value of the outer scale of turbulence, indicating that K₀ can be neglected (Galiègue et al., 2016):

$$\langle {\chi }^2{\rangle }_{\mathrm{asym}}^{\mathrm{PW}}=\frac{2^{1-p_m/2}{\pi }^{5/2}{\lambda }^2{r}_e^2{R}_{\mathrm{iono}}{a}_x^{3-p_m}{A}_y{A}_z{C}_s\mathrm{\Gamma }(3/2-p_m/4)}{(p_m/2-1)\mathrm{\Gamma }(p_m/4){B\text{'}}^{p_m/2}}{\left(\frac{L_v}{k_0}\right)}^{p_m/2-1}{}_2{\mathrm{F}}_1\left(\frac{p_m}{2},\frac{1}{2};1;1-\frac{A\text{'}}{B\text{'}}\right),$$(30)

$$\langle {\phi }^2{\rangle }_{\mathrm{asym}}^{\mathrm{PW}}=4{\pi }^2{r}_e^2{\lambda }^2{R}_{\mathrm{iono}}{a}_x{A}_y{A}_z{C}_s\frac{K_0^{2-p_m}}{\left(p_m-2\right)\sqrt{\mathrm{AB}-C^2}}-\langle {\chi }^2{\rangle }_{\mathrm{asym}}^{\mathrm{PW}},$$(31)

where ₂F₁ is the hypergeometric function, A′ and B′ are anisotropic coefficients, see (Galiègue et al., 2016). Note that in (30), the Fresnel distance for an incident plane $r_{\mathrm{F}}^{\mathrm{PW}}=\sqrt{\lambda L_v}$ implicitly appears in the term $\frac{L_v}{k_0}$:

$$\frac{L_v}{k_0}=\frac{(r_{\mathrm{F}}^{\mathrm{PW}}{)}^2}{2\pi }.$$(32)

This Fresnel distance has a fundamental impact on amplitude scintillation (Rino, 1979; Carrano et al., 2019; Morel et al., 2024).

3.1.3 Incident spherical wave approach

For an incident spherical wave, the spectral weighting function for log-amplitude and phase is now given by (Wheelon, 2003):

$${\mathcal{F}}_{\chi }^{\mathrm{SW}}(k_u,k_v)=\frac{2}{R_{\mathrm{iono}}}{\int }_{L_t}^{L_t+R_{\mathrm{iono}}}\mathrm{}{\sin }^2\left(\frac{s(R-s)}{2k_{0}R}\left(k_u^2+k_v^2\right)\right)\mathrm{d}s,$$(33)

$${\mathcal{F}}_{\phi }^{\mathrm{SW}}(k_u,k_v)=\frac{2}{R_{\mathrm{iono}}}{\int }_{L_t}^{L_t+R_{\mathrm{iono}}}\mathrm{}{\cos }^2\left(\frac{s(R-s)}{2k_{0}R}(k_u^2+k_v^2)\right)\mathrm{d}s,$$(34)

with the transmitter-receiver distance R = L_v + R_iono + L_v. Then, after some mathematical developments, it comes:

$$\begin{array}{c}\begin{array}{cc}{\mathcal{F}}_{\chi }^{\mathrm{SW}}(k_u,k_v)& =1-\frac{2}{R_{\mathrm{iono}}}\sqrt{\frac{k_0(L_t+R_{\mathrm{iono}}+L_v)}{4(k_u^2+k_v^2)}}\times \end{array}\\ \begin{array}{cc}& \left[\cos \left(\frac{(L_t+R_{\mathrm{iono}}+L_v)(k_u^2+k_v^2)}{4k_0}\right)\left(\mathrm{C}\right(t_2)-\mathrm{C}(t_1\left)\right)+\sin \left(\frac{(L_t+R_{\mathrm{iono}}+L_v)(k_u^2+k_v^2)}{4k_0}\right)\left(\mathrm{S}\right(t_2)-\mathrm{S}(t_1\left)\right)\right],\end{array}\end{array}$$(35)

$$\begin{array}{cc}{\mathcal{F}}_{\phi }^{\mathrm{SW}}(k_u,k_v)& =1+\frac{2}{R_{\mathrm{iono}}}\sqrt{\frac{k_0(L_t+R_{\mathrm{iono}}+L_v)}{4(k_u^2+k_v^2)}}\times \\ & \left[\cos \left(\frac{(L_t+R_{\mathrm{iono}}+L_v)(k_u^2+k_v^2)}{4k_0}\right)\left(\mathrm{C}\right(t_2)-C(t_1\left)\right)+\sin \left(\frac{(L_t+R_{\mathrm{iono}}+L_v)(k_u^2+k_v^2)}{4k_0}\right)\left(\mathrm{S}\right(t_2)-\mathrm{S}(t_1\left)\right)\right],\end{array}$$(36)

with $t_1=\frac{\sqrt{k_u^2+k_v^2}(L_t-(R_{\mathrm{iono}}+L_v\left)\right)}{2\sqrt{k_0(L_t+R_{\mathrm{iono}}+L_v)}}$, $t_2=\frac{\sqrt{k_u^2+k_v^2}(L_t+R_{\mathrm{iono}}-L_v)}{2\sqrt{k_0(L_t+R_{\mathrm{iono}}+L_v)}}$; and S() and C() the Fresnel integrals defined by $S\left(x\right)={\int }_0^x\mathrm{}\sin \left(t^2\right)\mathrm{d}t$ and $C\left(x\right)={\int }_0^x\mathrm{}\cos \left(t^2\right)\mathrm{d}t$. As expected, it can be shown that the spectral weighting functions for the incident spherical approach tend to the spectral weighting functions for the incident plane wave approach as L_t tends to infinity. The details of this demonstration are provided in Appendix B.

3.1.4 Correction of the incident plane wave approach

In the asymptotic formulation of the log-amplitude variance (30), it has been observed that the Fresnel distance r_F appears. So, in order to correct the approximation of an incident plane wave in the variance computations, one empirical solution is to replace the plane wave Fresnel distance by the spherical Fresnel distance $r_{\mathrm{F}}^{\mathrm{SW}}=\sqrt{\frac{\lambda L_v L_t}{L_v+ L_t}}$. This substitution leads to multiplying $r_{\mathrm{F}}^{\mathrm{PW}}$ by $\sqrt{\frac{L_t}{L_v+ L_t}}$. This correction is proposed in (Lawrence et al., 1964). Following the same heuristic approach, another correction can be thought for log-amplitude and phase indices. Indeed, by looking at the spectral weighting functions for plane wave, it is possible to express X with the Fresnel radius: $X=\frac{k^2(r_{\mathrm{F}}^{\mathrm{PW}}{)}^2}{2\pi }$. Thereby, to correct the incident plane wave formalism into an incident spherical wave formalism, the Fresnel distance can be substituted in the argument: $X_{\mathrm{cor}}=\frac{k^2{\left(r_{\mathrm{F}}^{\mathrm{SW}}\right)}^2}{2\pi }=\frac{k^2{L}_v{L}_t}{k_0(L_v+L_t)}$. By correcting the weighting functions, it is possible to use this correction for phase and log-amplitude variances as well as power spectral densities (Morel et al., 2025). In order to validate the analytic results obtained under the incident spherical wave formalism, the Parabolic Wave Equation with Multiple Curved Phase Screens (PWE-MCPS) can be used.

4 Parabolic Wave Equation with multiple phase screens

4.1 Propagation scheme

In this section, the theory allowing the propagation of an incident spherical wave using the (PWE-MCPS) is described (Grimault, 1998; Knepp, 2016). The equation of the numerical scheme, described in Figure 3, is given by:

$$u(r_{i+1},\theta ,\alpha )={\mathrm{e}}^{j{\varphi }_i(\theta ,\alpha ,t)}{\mathrm{F}}^{-1}\left\{{\mathrm{e}}^{-j\frac{1}{2k_0}(\frac{1}{r_i}-\frac{1}{r_{i+1}})(k_{\theta }^2+k_{\alpha }^2)}\mathrm{F}\right\{u(r_i,\theta ,\alpha )\left\}\right\},$$(37)

Figure 3 Scheme of the 3-D PWE model with 2-D multiple curved phase screens (3D PWE/2D MCPS) in spherical coordinates (r, θ, α), with $\theta =\frac{\pi }{2}-{\theta }_c$.

The first term ${\mathrm{e}}^{j{\varphi }_i(\theta , \alpha , t)}$ corresponds to the phase screen at range r_i. This equation characterizes a Split-Step Fourier (SSF) scheme to solve the propagation between r_i and r_i+1 = r_i + δr, associated with phase screen ϕ_i representing the medium variability. The generation of such a random phase screen is shortly described within the next section.

4.2 Phase screen spectrum

Recalling that 〈Δn(r, t)〉 = 0 and ${\varphi }_i(\theta ,\alpha ,t)=k_0{\int }_{r_i}^{r_i+\delta r}\mathrm{}\Delta n(r,\theta ,\alpha ,t)\mathrm{d}r$, it follows that ϕ is a centred 2D random variable which 2D spatial covariance function B_ϕ (θ, α) = 〈ϕ(θ′, α′) ϕ(θ′ + θ, α′ + α)〉 is given by (Fabbro & Feral, 2012):

$$B_{\varphi }\left(\theta ,\alpha \right)=k_0^2{\int }_{r_i}^{r_i+\delta r}\mathrm{}{\int }_{r_i}^{r_i+\delta r}\mathrm{}\langle \Delta n\left(r_1,{\theta }^{\mathrm{\prime }},{\alpha }^{\mathrm{\prime }},t\right)\Delta n\left(r_2,{\theta }^{\mathrm{\prime }}+\theta ,{\alpha }^{\mathrm{\prime }}+\alpha ,t\right)\rangle ,$$(38)

$$B_{\varphi }(\theta ,\alpha )=k_0^2{\int }_{r_i}^{r_i+\delta r}\mathrm{}{\int }_{r_i}^{r_i+\delta r}\mathrm{}{B}_{\Delta n}(r_2-r_1,\theta ,\alpha )\mathrm{d}{r}_1\mathrm{d}{r}_2.$$(39)

Changing the integration variable to sum and difference, equation (39) becomes:

$$B_{\varphi }(\theta ,\alpha )=k_0^2\delta r{\int }_{-\delta r}^{\delta r}\mathrm{}{B}_{\Delta n}(w,\theta ,\varphi )\left(1-\frac{\left|w\right|}{\delta r}\right)\mathrm{d}w.$$(40)

If δr is large with respect to the correlation distance of the turbulent refractive index, the integration can be extended to infinity so that:

$$B_{\varphi }(\theta ,\alpha )=k_0^2\delta r{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}{B}_{\Delta n}(w,\theta ,\alpha )\mathrm{d}w,$$(41)

and using relation (3), it comes:

$$B_{\varphi }(\theta ,\alpha )=r_e^2{\lambda }^2\delta r{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}{B}_{\Delta N_e}(w,\theta ,\alpha )\mathrm{d}w.$$(42)

So, by taking the Fourier transform, the spectrum of the phase screen can be obtained:

$$\begin{array}{cc}{S}_{\varphi }(k_{\theta },k_{\alpha })& =\frac{1}{(2\pi {)}^2}{\int }_0^{2\pi }\mathrm{}{\int }_0^{\pi }\mathrm{}{B}_{\varphi }(\theta ,\alpha ){\mathrm{e}}^{j(k_{\theta }\theta +k_{\alpha }\alpha )}\mathrm{d}\theta \mathrm{d}\alpha ,\\ & =\frac{r_e^2{\lambda }^2\delta r}{(2\pi {)}^2}{\int }_0^{2\pi }\mathrm{}{\int }_0^{\pi }\mathrm{}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}{B}_{\Delta N_e}(w,\theta ,\alpha ){\mathrm{e}}^{j(k_{\theta }\theta +k_{\alpha }\alpha )}\mathrm{d}w\mathrm{d}\theta \mathrm{d}\alpha ,\\ \begin{array}{c}\\ \end{array}& \begin{array}{c}=\frac{r_e^2{\lambda }^2\delta r}{(2\pi {)}^2{r}^2}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}{B}_{\Delta N_e}(w,\theta ,\alpha ){\mathrm{e}}^{j(k_{\theta }\frac{\widehat{u}}{r}+k_{\alpha }\frac{\widehat{v}}{r})}\mathrm{d}w\mathrm{d}\widehat{u}\mathrm{d}\widehat{v},\\ =2\pi \frac{r_e^2{\lambda }^2\delta r}{r^2}{S}_{\Delta N_e}\left(0,\frac{k_{\theta }}{r},\frac{k_{\alpha }}{r}\right).\end{array}\end{array}$$(43)

It can be noted that this expression is similar to that obtained for the incident Plane Wave formalism (Rino, 2011; Galiègue et al., 2016). Equation (43) defines the transverse spectrum S_ϕ (k_θ, k_α) that is used to generate the turbulent random phase screens ϕ(θ, α). The numerical implementation of the phase screen consists in filtering a white Gaussian noise by $\sqrt{S_{\varphi }(k_{\theta },k_{\alpha })}$ as explained in Macaskill & Ewart (1984).

5 Comparison of the different approaches

In the previous sections, we have presented different approaches to model scintillation for a wave illuminating the ionosphere irregularities: The incident plane wave approach (PW), the incident spherical wave (SW), and a heuristic approach consisting of correcting the Fresnel distance in the filtering functions of the plane wave formalism (CPW). These three approaches are based on Rytov theory and provide theoretical formulations of log-amplitude and phase variances in a turbulent ionosphere. The last approach, PWE-MCPS, is numerical, and a Monte Carlo process is used to derive statistical moments. These four approaches are compared in this section.

The parameters used for all simulations are given in Table 1. This configuration corresponds to a LEO satellite moving over the polar region where ionospheric scintillation phenomena are frequent. The thickness and altitude of the ionosphere are set to ΔH = 20 km and H = 350 km, respectively, which allows respecting the thin layer approximation. The integrated strength of the turbulence $C_{k}L=\left(2\pi {)}^3\right(\frac{1000}{2\pi }{)}^{p_m}\mathrm{\Delta }\mathrm{H}{C}_s$ is set to 10³⁴, which corresponds to moderate scintillation.

Figure 4 shows the scintillation indices S₄ (left figures, which are related to the log-amplitude variance: $S_4\approx 2\sqrt{\langle {\chi }^2\rangle }$) and σ_φ (right figures, which are related to the phase variance: ${\sigma }_{\phi }\approx \sqrt{\langle {\phi }^2\rangle }$) obtained by the four methods, whereas the transmitter-inhomogeneous layer distance L_t is varying. The top graphs of Figure 4 represent the indices versus L_t, whereas the bottom graphs represent the relative errors in percent referenced to the results of the incident spherical wave Rytov formalism, chosen to be the reference under a weak scattering, regime i.e.:

$$\frac{|S_4-S_4^{\mathrm{SW}}|}{S_4^{\mathrm{SW}}}=\frac{\Delta S_4}{S_4^{\mathrm{SW}}},$$(44)

$$\frac{|{\sigma }_{\phi }-{\sigma }_{\phi }^{\mathrm{SW}}|}{{\sigma }_{\phi }^{\mathrm{SW}}}=\frac{\Delta {\sigma }_{\phi }}{{\sigma }_{\phi }^{\mathrm{SW}}}.$$(45)

Figure 4 Top: Scintillation indices for an incident plane wave (blue curve), an incident spherical wave (red curve), an incident corrected plane wave (green curve) and PWE-MCPS result (cyan curve). Bottom: Relative error with respect to the incident spherical wave model. The satellite altitude varies between Hsat = H + ΔH = 370 km (on top of the ionospheric layer) and Hsat = 5,000 km.

Scintillation indices are computed from equations (20) and (21) with filtering functions for incident plane wave approach (Eqs. (28) and (29)) and incident spherical wave approach (Eqs. (35) and (36)). For both indices S₄ and σ_φ, the values obtained with the incident spherical wave approach tend, as expected, to the values obtained with the plane wave approach for large values of L_t, this behavior is demonstrated theoretically in Appendix B. As illustrated by the relative errors (bottom graphs of Fig. 4), the impact of the incident wave sphericity is more important on the amplitude than on the phase, with a dynamic of the evolution of the S₄ index stronger than the evolution of the σ_φ index. On the one hand, the relative error on the S₄ index is important, typically for LEO satellites, H_sat ≤ 2,000 km i.e., L_t ≤ 1,670 km, the relative error for the incident plane wave model is greater than 10% and even exceeds 100%, while the incident corrected plane wave model gives much better results with a relative error of less than 1% as soon as the satellite altitude is higher than 540 km. On the other hand, the relative error for the phase scintillation index is less than 3% regardless of the propagation model; nevertheless, the incident corrected plane wave approach gives better results with a relative error below 0.1% even at low altitude. We also notice a good match between the analytic and numerical models for the amplitude and phase scintillation indices. This point gives confidence in the simulation results. Finally, one can remark that the evolution of the scintillation indices with L_t is known analytically. Indeed, looking at equations (30) and (31), the plane wave Fresnel distance can be replaced by the spherical Fresnel distance, and it comes:

$${\langle {\chi }^2\rangle }_{\mathrm{asym}}^{\mathrm{SW}}\propto {\left(r_{\mathrm{F}}^{\mathrm{SW}}\right)}^{p_m-1/2}={\left(\lambda \frac{L_v{L}_t}{L_v+L_t}\right)}^{p_m/2-1},$$(46)

$${\langle {\phi }^2\rangle }_{\mathrm{asym}}^{\mathrm{SW}}\propto \mathrm{Cst}-{\left(r_{\mathrm{F}}^{\mathrm{SW}}\right)}^{p_m-1/2}=\mathrm{Cst}-{\left(\lambda \frac{L_v{L}_t}{L_v+L_t}\right)}^{p_m/2-1},$$(47)

where Cst is a constant depending on ionospheric and geometric parameters and can be derived from (30) and (31). These analytic trends are confirmed by the results in Figure 4.

The propagation channel in point-to-point geometry has the property of being reciprocal. To verify that this property is well respected by the scintillation models, another test case is proposed here.

Figure 5 illustrates the scintillation indices S₄ and σ_φ (top and bottom graphs, respectively) for the PW and SW models, plotted as functions of L_t at a fixed satellite altitude H_sat = 600 km. All other configuration parameters are detailed in Table 1. In Figure 5, the red curves, which correspond to the SW model, display a clear symmetry around L_t = R/2. This symmetry reflects the reciprocity of the propagation channel when using the incident spherical wave approach: Whether the distance from the satellite to the ionospheric layer is L_t or R − L_t, the scintillation indices remain unchanged. This reciprocity results from the symmetry in the filtering functions, as shown in equations (33) and (34), the integrals can be calculated over either the range L_t to L_t + R_iono or L_v to L_v + R_iono, yielding identical results. This is not the case, however, for the incident plane wave model, as seen in equations (22) and (23). In contrast, the blue curves in Figure 5, representing the PW model, do not show the same values for L_t and R − L_t, indicating a non-reciprocal propagation channel. For the incident plane wave approach, the log-amplitude and phase variances vary in opposite ways.

Figure 5 Scintillation indices for incident plane wave (blue curve) and spherical wave (red curve) approaches. The satellite altitude is set at Hsat = 600 km, the distance Lt between satellite and irregularities varies between 0 and R − Riono = 598 km. At Lt = 0 km, the ionospheric layer is very close to the satellite; at Lt = 598 km, the ionospheric layer is very close to the receiver.

6 Impact of real system configuration on ionospheric scintillation

In the previous sections, the simulations performed were not fully representative of a real system. First, the ionospheric medium velocity relative to the LOS axis was identical, regardless of the altitude of the satellite. However, this velocity is the sum of the ionosphere drift velocity and satellite projection velocity on the ionospheric layer. Indeed, the satellite’s speed depends on its altitude in the hypothesis of a circular orbit $V_{\mathrm{sat}}=\sqrt{\frac{\mu }{R_T+H_{\mathrm{sat}}}}$, with μ [m³ s⁻²] the standard gravitational parameter and R_T [m] the Earth radius. Thus, we add to the drift velocity of the medium the scan velocity of the ionosphere. From geometric consideration in Figure 6, the scan velocity of the ionosphere can be expressed as:

$$V_{\mathrm{scan}}=\frac{\beta }{\alpha }\frac{R_T+\mathrm{H}+\mathrm{\Delta }\mathrm{H}}{R_T+H_{\mathrm{sat}}}{V}_{\mathrm{sat}},$$(48)

where α is the angle between the local vertical at the receiver point and the local vertical passing through the satellite, and β the angle between the local vertical at the receiver point and the local vertical passing through the top of the ionospheric layer. In compliance with Figure 6, α and β can be expressed as:

$$\alpha =\arccos \left(\frac{R_T^2+(R_T+H_{\mathrm{sat}}{)}^2-R^2}{2R_T(R_T+H_{\mathrm{sat}})}\right),$$(49)

$$\beta =\arccos \left(\frac{R_T^2+(R_T+H+\Delta H{)}^2-(L_v+R_{\mathrm{iono}}{)}^2}{2R_T(R_T+H+\Delta H)}\right).$$(50)

Figure 6 Description of the satellite-ground link with angles α and β.

In our simulations, the altitude of the satellite H_sat varies from 370 to 20,000 km, so its speed V_sat varies from about 7.5 to 4 km s⁻¹, and the scan speed of the ionosphere varies from about 7.5 to 0.1 km s⁻¹. It is important to specify that the cutoff frequency of the scintillation measurements (with a typical value of 0.1 [Hz] for GNSS) is necessary to take into account.

As a reminder, the scintillation indices can be computed with the integral of the log-amplitude and phase spectra whose expressions are given in (Morel et al., 2024, 2025):

$$S_4^2=4{\int }_{f_{\mathrm{cut}}}^{f_{\max }}\mathrm{}{W}_{\chi }\left(f\right)\mathrm{d}f,$$(51)

$${\sigma }_{\phi }^2={\int }_{f_{\mathrm{cut}}}^{f_{\max }}\mathrm{}{W}_{\phi }\left(f\right)\mathrm{d}f,$$(52)

so the scintillation indices could be impacted by the satellite velocity and the cutoff frequency.

Figure 7 shows the scintillation indices relative to satellite altitude and velocity (top graphs) and the percentage deviation from the value at 375 skm, i.e., (bottom graphs):

$$\frac{|S_4-S_4(375\left)\right|}{S_4\left(375\right)}=\frac{\Delta S_4}{S_4\left(375\right)},$$(53)

$$\frac{|{\sigma }_{\phi }-{\sigma }_{\phi }(375\left)\right|}{{\sigma }_{\phi }\left(375\right)}=\frac{\Delta {\sigma }_{\phi }}{{\sigma }_{\phi }\left(375\right)}.$$(54)

Figure 7 Top: Scintillation indices for an incident spherical wave versus the satellite altitude and speed. Bottom: Percentage deviation from the scintillation indices at Hsat = 375 km.

Firstly, we observe that the evolution of the amplitude scintillation index S₄ is similar to the one shown in Figure 4, meaning that the speed and the cutoff frequency do not influence the amplitude scintillation. For this particular configuration, a LEO satellite amplitude scintillation index S₄ is about 0.085, whereas a MEO satellite amplitude scintillation index S₄ is about 0.14, which means that the farther away the satellite is, the higher the amplitude scintillation. Secondly, the dynamic of evolution for the phase scintillation index is a little bit stronger than the one presented in Figure 4. Indeed, for this configuration, a LEO satellite has a phase scintillation index σ_φ of about 0.292 which must be compared to a value of about 0.263 for a MEO satellite, representing a deviation of about 10%. These results are consistent with (Morton et al., 2022), where it is demonstrated that “the level of scintillation severity is determined by how close the ionospheric irregularities are located to the LEO satellite orbit: The closer they are, the stronger the scintillation”.

7 Conclusion

Errors potentially introduced by applying a hypothesis of a plane wave illuminating the ionosphere instead of a spherical wave to model 3D scintillation effects were evaluated quantitatively for a typical Earth-satellite system configuration at L-band. To carry out this study, different models were considered. On the one hand, starting from the 3D scalar propagation equation, the log-amplitude and phase variances have been derived under Rytov formalism for a Plane Wave and a Spherical Wave illuminating the ionospheric irregularities. On the other hand, as part of a numerical approach, the scintillation effects have been evaluated from a 3D Parabolic Wave Equation (PWE) approach associated with Multiple Curved Phase Screen (MCPS).

Based on realistic scenarios inducing ionospheric scintillation effects, comparisons of different models considering a spherical or a plane incident wave have been conducted.

First, scintillation effects derived from 3D-PWE/2D-MCPS numerical schemes match very well the corresponding 3D analytic formulations. The proposed method to correct the incident plane wave formalism into the incident spherical wave one shows very good results on scintillation indices, even at low altitude. Secondly, it has been shown that at low altitude (L_t < 500 km), when the satellite is close to the ionospheric layer, the relative error of the incident plane wave compared to spherical wave formalism is more than 50% for the amplitude scintillation index, whereas the impact on the phase scintillation is negligible. Thirdly, the propagation channel reciprocity for incident spherical wave has been shown, which is not the case for incident plane wave. Thus, two-way path propagation scintillation effects must be studied with the incident spherical wave model or by applying a correction of the Fresnel distance into a model considering an incident plane wave. Finally, a study with more realistic conditions, taking the satellite velocity into account, showed that the phase scintillation index is about 10% higher for LEO satellites than for GNSS satellites, and the amplitude scintillation index is about 40% lower.

Acknowledgments

The authors would like to thank the Direction Générale de l’Armement (DGA), the Agence de l’Innovation de Défense (AID), for their financial support. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

Funding

Part of this work was funded by the Agence de l’Innovation de Défense (AID).

Conflicts of interest

The authors declare no Conflict of Interest.

Data availability statement

No observational or model data were used for or generated by this work.

References

Appendix A

Anisotropy coefficient

This Appendix has been added for the sake of completeness. The expression anisotropy coefficients in the Line Of Sight coordinates system are (Galiègue, 2015; Galiègue et al., 2016)

$$\begin{array}{c}A=\left(\sin \right(\gamma \left)\cos \right({\alpha }_z\left)\sin \right(\psi )+\cos (\gamma \left)\cos \right({\alpha }_y\left)\cos \right(\psi ){)}^2\\ +A_y^2{\cos }^2\left(\psi \right){\sin }^2\left({\alpha }_y\right)+A_z^2{\sin }^2\left(\gamma \right){\sin }^2\left({\alpha }_z\right),\end{array}$$(A.1)

$$\begin{array}{c}B=\left(\cos \right(\gamma \left)\sin \right({\alpha }_y\left)\cos \right(\psi )+\sin (\gamma \left)\sin \right({\alpha }_z\left)\sin \right(\psi ){)}^2\\ +A_y^2{\cos }^2\left(\psi \right){\cos }^2\left({\alpha }_y\right)+A_z^2{\sin }^2\left(\gamma \right){\cos }^2\left({\alpha }_z\right),\end{array}$$(A.2)

$$\begin{array}{c}C=-\left(\sin \right(\gamma \left)\cos \right({\alpha }_z\left)\sin \right(\psi )+\cos (\gamma \left)\cos \right({\alpha }_y\left)\cos \right(\psi \left)\right)\\ \times \left(\cos \left(\gamma \right)\sin \left({\alpha }_y\right)\cos \left(\psi \right)+\sin \left(\gamma \right)\sin \left({\alpha }_z\right)\sin \left(\psi \right)\right)\\ +A_y^2{\cos }^2\left(\psi \right)\cos \left({\alpha }_y\right)\sin \left({\alpha }_y\right)+A{z}^2{\sin }^2\left(\gamma \right)\cos \left({\alpha }_z\right)\sin \left({\alpha }_z\right).\end{array}$$(A.3)

The coefficients A′ and B′, which appear in the asymptotic formulation of log-amplitude and phase variances, are (Galiègue, 2015; Galiègue et al., 2016):

$$A\mathrm{\prime }=A{\cos }^2\left(ϵ\right)+B{\sin }^2\left(ϵ\right)+2C\cos \left(ϵ\right)\sin \left(ϵ\right),$$(A.4)

$$B\mathrm{\prime }=A{\sin }^2\left(ϵ\right)+B{\cos }^2\left(ϵ\right)-2C\cos \left(ϵ\right)\sin \left(ϵ\right),$$(A.5)

with ϵ the angle between the u axis and u′ axis, which is the major axis of the ellipse.

Appendix B

Asymptotic weighting functions

Starting with asymptotic formulas for the Fresnel integrals for t tends to infinity:

$$\begin{array}{c}S\left(t\right)={\int }_0^t\mathrm{}\sin \left(x^2\right)\mathrm{d}x\sim \sqrt{\frac{\pi }{8}}-\frac{\cos \left(t^2\right)}{2t},\\ C\left(t\right)={\int }_0^t\mathrm{}\cos \left(x^2\right)\mathrm{d}x\sim \sqrt{\frac{\pi }{8}}+\frac{\sin \left(t^2\right)}{2t},\end{array}$$

and setting $k^2=k_u^2+k_v^2$, the cosine and sine terms in brackets in equations (35) and (36) can be expressed as:

$$\begin{array}{c}\cos \left(\frac{(L_t+R_{\mathrm{iono}}+L_v)k^2}{4k_0}\right)\left(C\right(t_2)-C(t_1\left)\right)\approx \cos \left(\frac{(L_t+R_{\mathrm{iono}}+L_v)k^2}{4k_0}\right)\frac{\sin \left(t_2^2\right)}{2t_2}\\ -\cos \left(\frac{(L_t+R_{\mathrm{iono}}+L_v)k^2}{4k_0}\right)\frac{\sin \left(t_1^2\right)}{2t_1},\end{array}$$(B.1)

$$\begin{array}{c}\sin \left(\frac{(L_t+R_{\mathrm{iono}}+L_v)k^2}{4k_0}\right)\left(S\right(t_2)-S(t_1\left)\right)\approx \sin \left(\frac{(L_t+R_{\mathrm{iono}}+L_v)k^2}{4k_0}\right)\frac{\cos \left(t_2^2\right)}{2t_2}\\ -\sin \left(\frac{(L_t+R_{\mathrm{iono}}+L_v)k^2}{4k_0}\right)\frac{\cos \left(t_1^2\right)}{2t_1}.\end{array}$$(B.2)

The sum of (B.1) and (B.2) gives:

$$[...]\approx \frac{1}{2t_2}\sin \left(t_2^2-\frac{(L_t+R_{\mathrm{iono}}+L_v)k^2}{4k_0}\right)-\frac{1}{2t_1}\sin \left(t_1^2-\frac{(L_t+R_{\mathrm{iono}}+L_v)k^2}{4k_0}\right).$$(B.3)