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

1 Introduction

This paper is the second part of a study on ionospheric scintillation modeling and analysis. It addresses the power spectral analysis (PSD) of ionospheric scintillation in a plane wave (PW) and a spherical wave (SW) approach. The first part focuses on the analysis of scintillation indices and is the topic of a companion paper.

Random variability of the received electromagnetic field due to fluctuations of electronic density in the ionosphere can cause degradation in positioning and communication system performances. Thus, accurate knowledge of the amplitude and phase signal received can help data processing research to mitigate these effects. They are related to the diffusion of electromagnetic waves in the ionosphere, which can present strong turbulent irregularities. Characterizing these ionospheric irregularities and their impact on the propagated field remain a challenge. However, the theoretical interpretation of the variability of the received signal can provide information about the medium through which the waves have propagated.

Nowadays, many two-frequency radionavigation systems allow for Total Electron Content (TEC), amplitude and phase scintillation measurements of the signal. The set of these measurements potentially contains a wide range of information describing the ionospheric medium and its variability. Mastering an analytic formalism to describe signal variations in response to medium fluctuations is essential for tackling an inversion problem effectively. In addition to understanding the current ionospheric conditions, mastering the long-term climatology of amplitude and phase scintillation phenomena allows a better design of electromagnetic systems by associating an operational availability. Statistics of scintillation events are therefore relevant. These are linked to the long-term solar activity and its cycle, the Earth’s magnetic activity, the season, the local time, the position on the Earth and the geometry of the radio-electric link with respect to the Earth’s magnetic field.

In the past, combined observations of incoherent scatter and scintillation observations have led to a better understanding of the morphology of ionospheric irregularities (Aarons, 1982; Livingston et al., 1982; Rino et al., 1983). In this context, the theoretical work of modelling the scintillation indices in amplitude, S₄ and phase, σ_φ, were particularly useful (Rino, 1979). More recently, in Carrano et al. (2016) the authors propose a method to estimate the zonal irregularity drift from measurements of amplitude and phase scintillation indices by a single GNSS station. The Fresnel frequency f_F has also been studied by several authors to characterize ionospheric medium irregularities (McCaffrey & Jayachandran, 2019; Song et al., 2022; Sun et al., 2023). It is given by $f_{\mathrm{F}}=V/\left(\sqrt{2}{r}_{\mathrm{F}}\right)$ with V the relative layer drift (composed by ionospheric piercing point of the LOS and the irregularity drift), $r_{\mathrm{F}}=\sqrt{\lambda L_v}$ the Fresnel distance, where λ is the wavelength and L_v the ionospheric layer-receiver distance. In the literature, the Fresnel frequency is sometimes associated with the frequency of the maximum of the log-amplitude PSD (that is also called rollover frequency) (Forte & Radicella, 2002). In this paper, the limits of this assumption are specifically studied. Going forward, to describe morphological and spectral features of irregularities, inversion of GPS observations has been used. For instance, in Conroy et al. (2021), Vaggu et al. (2023) an inverse method is applied to derive the parameters that describe the ionospheric layer irregularities. The approach estimates the best fit of numerical model outputs to GPS observations. The numerical model is a 3D-PWE/2D-MCPS (Parabolic Wave Equation method with Multiple Curved Phase Screens) resolution of the Helmholtz propagation equation, and the inversion approach is based on spectral analysis. This approach requires the use of a numerical regression algorithm. To optimize computation times and facilitate the inversion of results, the development of an analytic formalism is fundamental.

The modelling of the amplitude and phase of the received signal has been the subject of intensive theoretical research and several approaches can be found in the open literature (Yeh & Liu, 1982; Vasylyev et al., 2022). Notably, one can cite the work of Tatarskii (1971), Rytov et al. (1989) and Wheelon (2003), which led to an asymptotic approach under low disturbance and thin layer hypothesis. One strong feature of these formulations is that they would allow for attempts to access certain parameters, describing the statistics of variability of the ionospheric environment.

The originality of this paper, by comparison to Fabbro & Feral (2012) and Galiègue (2015), is to develop analytic formulations of log-amplitude and phase PSD, which could be used to make inversion of ionospheric parameters describing irregularities. To facilitate the inversion and interpretation of spectra, asymptotic formulations at low and high frequencies are also proposed, as well as an approximate formulation of the maximum of the log-amplitude spectrum. The goal is to provide to the community key formulations or parameters, giving more direct access to the description of the ionospheric turbulence. The analysis is carried out in terms of log-amplitude and phase Power Spectral Densities (PSD). To validate the presented results, both 3D-PWE/2D-MCPS numerical scheme and Rytov’s theory of smooth perturbation results for spherical waves are compared, and differences between PSD in plane and spherical wave formalisms are also studied. In addition, analytic derivations of the frequency of the maximum of the log-amplitude PSD obtained under the weak scattering and thin-layer hypotheses are also proposed.

Knowing the spectral characteristics of the ionospheric medium makes it possible to better describe the dynamics of the scintillation of the received signal, to describe its anisotropy, to understand its cascade of energy in the wavenumer domain by studying its inertial slope, and to know the size of irregularities via its outer scale of turbulence. As the measured scintillation, through the amplitude and phase scintillation indices (S₄ and σ_φ) studied in the first paper, have different dynamics of variation and different statistical occurrence depending on the geophysical conditions, we can then imagine different characteristics of the PSD of the irregularities depending on these conditions. Knowledge of this local climatology can be particularly relevant for system design.

The paper is organized as follows: Section 2 is focused on analytic formulations of the ionospheric scintillation impact in terms of log-amplitude and phase PSD for both incident Plane Wave (PW) and Spherical Wave (SW). In addition to these formulations, a heuristic approach, called Corrected Plane Wave (CPW) is proposed as a modification of the plane wave formalism. The frequency of the maximum of the log-amplitude spectrum is then derived, and its accuracy is studied. In Section 3, the different approaches are compared (i.e., analytic incident PW/SW/CPW formulations), additionally, the sensitivity of the spectra to ionospheric parameters is analyzed. Finally, the different results obtained are summarized and discussed in Section 4 which concludes the paper.

2 PSD analytic formulations under weak scattering assumption

In Galiègue et al. (2016), the spectrum of the electronic density fluctuations $S_{\Delta N_e}$ has been derived in the line of sight coordinate system (u, v, s). It is given by:

$$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},$$(1)

where A, B, and C are the anisotropic coefficients expressed in the line of sight basis (cf. expressions in Galiègue et al., 2016). They depend on the link geometry: γ, α_z and ψ angles (see Fig. 1 of Morel et al., 2025) and anisotropic ratios A_y, A_z. a_x is a dimensionless scaling factor, $K_0=\frac{2\pi }{L0}$ [rad m⁻¹] is the wave number associated with the outer scale of the turbulence and p_m the slope of the Power Spectral Density of the irregularities, finally C_S $\left[{\mathrm{m}}^{-3-p_m}\right]$s the strength of the ionospheric turbulence (Rino, 2011; Galiègue et al., 2016).

From this description, the Rytov’s theory of smooth perturbations (Rytov et al., 1989; Wheelon, 2003), allows to derive analytic formulations of log-amplitude and phase spectra. The electric field E₁(R) can then be decomposed in terms of amplitude $\mathcal{A}\left(\mathbf{R}\right)$ and phase φ(R) with $E_1\left(\mathbf{R}\right)=\mathcal{A}\left(\mathbf{R}\right){\mathrm{e}}^{j({\phi }_0+\phi (\mathbf{R}\left)\right)}$. For more details, the reader should refer to Morel et al. (2025).

The signal first traverses a slant distance L_t from the transmitter to the upper part of the ionospheric layer, then covers a distance R_iono within the ionospheric layer, and finally travels a distance L_v from the lower part of the ionospheric layer to the ground station. The overall distance traveled by the signal is $R=L_t+R_{\mathrm{iono}}+L_v$. In this context, H corresponds to the altitude of the ionospheric layer, and ΔH indicates its thickness as described in more detail in Morel et al. (2025).

2.1 Log-amplitude and phase PSD

The log-amplitude W_χ and phase W_φ spectra are respectively the temporal Fourier transform of the log-amplitude and phase covariances $\langle \chi (\mathbf{R},t)\chi (\mathbf{R},t+\tau )\rangle$, and $\langle \phi (\mathbf{R},t)\phi (\mathbf{R},t+\tau )\rangle$, assuming that both processes are stationary (Fabbro & Feral, 2012):

$$W_{\chi }\left(\omega \right)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }\mathrm{}\langle \chi (\mathbf{R},t)\chi (\mathbf{R},t+\tau )\rangle {\mathrm{e}}^{-\mathrm{j\omega }\tau }\mathrm{d}\tau ,$$(2)

$$W_{\phi }\left(\omega \right)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }\mathrm{}\langle \phi (\mathbf{R},t)\phi (\mathbf{R},t+\tau )\rangle {\mathrm{e}}^{-\mathrm{j\omega }\tau }\mathrm{d}\tau ,$$(3)

where ω [rad s⁻¹] is the pulsation of the propagated wave, t and τ [s] re the time, (τ is the dual variable of ω in the Fourier space). The transformation between space and time is introduced by assuming Taylor’s hypothesis (Taylor, 1938): the medium is considered frozen and advected with a uniform speed V (Yeh & Liu, 1982). This velocity vector V is defined by its coordinates in the Line Of Sight (LOS) basis $\mathbf{V}=(V_u,V_v,V_s)$, which allows defining the drift speed transverse to the line of sight $V_r=\sqrt{V_u^2+V_v^2}$. In such conditions, it can be shown that (Ishimaru, 1972; Wheelon, 2003):

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

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

where F_χ and F_φ are the filtering functions in amplitude and phase respectively (Morel et al., 2025). By combining equations (2) and (4), and equations (3) and (5), it comes (Ishimaru, 1972; Wheelon, 2003; Fabbro & Feral, 2012):

$$W_{\chi }\left(\omega \right)=\frac{\pi r_e^2{\lambda }^2{R}_{\mathrm{iono}}}{V_u}{\int }_{-\infty }^{+\infty }\mathrm{}{S}_{\Delta N_e}\left(\frac{\omega -k_v{V}_v}{V_u},k_v,0\right){\mathcal{F}}_{\chi }\left(\frac{\omega -k_v{V}_v}{V_u},k_v\right)\mathrm{d}{k}_v,$$(6)

$$W_{\phi }\left(\omega \right)=\frac{\pi r_e^2{\lambda }^2{R}_{\mathrm{iono}}}{V_u}{\int }_{-\infty }^{+\infty }\mathrm{}{S}_{\Delta N_e}\left(\frac{\omega -k_v{V}_v}{V_u},k_v,0\right){\mathcal{F}}_{\phi }\left(\frac{\omega -k_v{V}_v}{V_u},k_v\right)\mathrm{d}{k}_v.$$(7)

In these formulations, the incident wave illuminating the ionosphere irregularites can be a spherical wave (SW) or a plane plane (PW). SW or PW specificity can be taken into account by changing the spectral filtering function, named ${\mathcal{F}}_{\chi ,\phi }^{\mathrm{PW}}$ and ${\mathcal{F}}_{\chi ,\phi }^{\mathrm{SW}}$, as realized in (Morel et al., 2025) to compute the amplitude and phase variances. The plane wave formulation of the filtering functions are:

$${\mathcal{F}}_{\begin{array}{l}\chi \\ \phi \end{array}}^{\mathrm{PW}}(k_u,k_v)=1\mp \mathrm{sinc}\left(\frac{(k_u^2+k_v^2){r_F^{\mathrm{PW}}}^2{R}_{\mathrm{iono}}}{4\pi L_v}\right)\cos \left(\frac{(k_u^2+k_v^2){r_F^{\mathrm{PW}}}^2}{2\pi }\left(1+\frac{R_{\mathrm{iono}}}{2L_v}\right)\right),$$(8)

with $r_F^{\mathrm{PW}}=\sqrt{\lambda L_v}$ the Fresnel radius for the incident plane wave. For the spherical wave formulation, (8) becomes:

$$\begin{array}{c}{\mathcal{F}}_{\begin{array}{l}\chi \\ \phi \end{array}}^{\mathrm{SW}}(k_u,k_v)=1\mp \frac{2}{R_{\mathrm{iono}}}\sqrt{\frac{k_0(L_t+R_{\mathrm{iono}}+L_v)}{4(k_u^2+k_v^2)}}\times \\ [\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)]\end{array},$$(9)

where t₁ and t₂ are given in the first paper ans C() and S() are the Fresnel integrals. This way, similarly to variance formulation correction in Morel et al. (2025), it is here proposed to correct the plane wave filtering functions by substitution of the Fresnel distance to get ${\mathcal{F}}_{\chi ,\phi }^{C\mathrm{PW}}$. The filtering functions are then given by:

$${\mathcal{F}}_{\begin{array}{l}\chi \\ \phi \end{array}}^{C\mathrm{PW}}(k_u,k_v)=1\mp \mathrm{sinc}\left(\frac{(k_u^2+k_v^2){r_F^{\mathrm{SW}}}^2{R}_{\mathrm{iono}}}{4\pi L_v}\right)\cos \left(\frac{(k_u^2+k_v^2){r_F^{\mathrm{SW}}}^2}{2\pi }\left(1+\frac{R_{\mathrm{iono}}}{2L_v}\right)\right),$$(10)

with $r_F^{\mathrm{SW}}=\sqrt{\frac{\lambda L_v{L}_t}{L_v+L_t}}$, the Fresnel radius for an incident spherical wave. The impact of the incident plane wave is analyzed in more details in Section 3.

2.2 Log-amplitude and phase PSD asymptotes

Whatever the filtering function (i.e., under SW, PW, or CPW approximation), the analytic expressions of the integrals in equations (6) and (7), are complicated to derive. However, it is possible to obtain asymptotic expressions at low and high frequencies that can be useful for inversion purposes. Assuming a thin turbulent layer compared to the layer-receiver distance (R_iono ≪ L_v) and considering simplified weighting functions, asymptotic formulations can be derived for the log-amplitude and phase spectra. The details of the mathematical derivations are given in Appendix A. They lead to the following high-frequency behavior of log-amplitude and phase PSD (similar for plane wave and spherical wave formalisms):

$$\begin{array}{c}{W}_{\begin{array}{l}\chi \\ \phi \end{array}}^{\mathrm{PW},\mathrm{SW}}(\omega \to +\mathrm{\infty })=\frac{\pi \sqrt{\pi }{r}_e^2{\lambda }^2{R}_{\mathrm{iono}}{a}_x^{3-p_m}{A}_y{A}_z{C}_s\mathrm{\Gamma }\left(\frac{p_m-1}{2}\right)}{\mathrm{\Gamma }(p_m/2)(A{V}_v^2+B{V}_u^2-2C{V}_u{V}_v{)}^{1-p_m/2}}\\ \times {\left({\omega }^2\left(\mathrm{AB}-C^2\right)+\frac{K_0^2}{a_x^2}\left(A{V}_v^2+B{V}_u^2-2C{V}_v{V}_u\right)\right)}^{\frac{1-p_m}{2}}.\end{array}$$(11)

We note that this expression can be simplified by assuming a large outer scale, and so a negligible K₀:

$$W_{\begin{array}{l}\chi \\ \phi \end{array}}^{\mathrm{PW},\mathrm{SW}}(\omega \to +\mathrm{\infty })=\frac{\pi \sqrt{\pi }{r}_e^2{\lambda }^2{R}_{\mathrm{iono}}{a}_x^{3-p_m}{A}_y{A}_z{C}_s\mathrm{\Gamma }\left(\frac{p_m-1}{2}\right)}{\mathrm{\Gamma }(p_m/2)(\mathrm{AB}-C^2{)}^{(p_m-1)/2}(A{V}_v^2+B{V}_u^2-2C{V}_u{V}_v{)}^{1-p_m/2}}{\omega }^{1-p_m}.$$(12)

At low-frequency, it is necessary to assume a large outer scale to derive the asymptotes, and so the asymptotic formulations become:

$$W_{\chi }^{\mathrm{PW},\mathrm{SW}}(\omega \to 0)=\frac{{\pi }^2{r}_e^2{\lambda }^2{R}_{\mathrm{iono}}{a}_x^{3-p_m}{A}_y{A}_z{C}_s{V}_r^{p_m-1}}{2(A{V}_v^2+B{V}_u^2-2C{V}_v{V}_u{)}^{p_m/2}\mathrm{\Gamma }\left(\frac{p_m+1}{2}\right)\cos \left(\pi -\frac{p_m+1}{4}\pi \right)}{\left(\frac{{r_F^{\mathrm{PW},\mathrm{SW}}}^2}{2\pi }\right)}^{(p_m-1)/2},$$(13)

$$W_{\phi }^{\mathrm{PW},\mathrm{SW}}(\omega \to 0)=\frac{2\sqrt{\pi }\pi r_e^2{\lambda }^2{R}_{\mathrm{iono}}{a}_x^2{A}_y{A}_z{C}_s{K}_0^{1-p_m}\mathrm{\Gamma }\left(\frac{p_m-1}{2}\right)}{\mathrm{\Gamma }(p_m/2)\sqrt{A{V}_v^2+B{V}_u^2-2C{V}_u{V}_v}}-W_{\chi }^{\mathrm{PW},\mathrm{SW}}(\omega \to 0),$$(14)

and these low-frequency asymptotes are valid for 1 < p_m < 5.

These relations will be validated by comparison to 3D-PWE/2D-MCPS numerical simulations, in the next subsection. Note that these asymptotic formulations (12), (13), and (14) are useful to derive turbulent parameters from experimental PSD.

2.3 Validation of log-amplitude and phase PSD modeling

In Sections 2.1 and 2.2, log-amplitude and phase PSD for an incident spherical wave illuminating the ionosphere are introduced. These expressions are based on Rytov theory. To validate these formulations, the 3D-PWE/2D-MCPS method is used, considering spherical phase screen as described in Morel et al. (2025), and a Monte-Carlo process is applied to derive statistical moments and average PSD. The polar configuration is used, and the associated parameters are given in Table 1. Note that 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 (ΔH ≪ H). The integrated strength of the turbulence C_kL is set to 10³⁴, which corresponds to moderate scintillation.

A validation example is proposed by comparison of log-amplitude (Fig. 1a) and phase PSD (Fig. 1b) for SW.

Figure 1 Comparison of theoretical and numerical PSD for SW formalism. (a) Log-amplitude PSD for SW formalism. (b) Phase PSD for SW formalism.

The polar configuration is considered (cf. Table 1), and log-amplitude and phase spectra are plotted in Figures 1a and 1b. These graphs show the comparison between theoretical Rytov PSD for incident spherical wave and incident corrected plane wave formalism (in red and red crosses, Eqs. (6), (7), (9), and (10)), asymptotes of incident spherical PSD (in black and magenta dashed lines respectively, Eqs. (12), (13), and (14)) and 3D-PWE/2D-MCPS (in cyan dashed line). The temporal 3D-PWE/2D-MCPS numerical spectra are obtained by applying Taylor’s principle along u-axis on spatial modeling of the propagated field. It must be pointed out that for better comparison and plot clarity, the satellite velocity has not been taken into account in this example, as this would significantly shifts the spectrum toward high frequencies.

Excellent matches are observed between theoretical and numerical curves. The simultaneous use of the two methods allows an intervalidation of the analytic and numerical spectra. For log-amplitude and phase PSD, there is a good match between high-frequency asymptotes and theoretical curves beyond 5 Hz. Analytic curves join low-frequency asymptotes for f < 0.1 Hz, with a small shift on the log-amplitude spectrum because K₀ is neglected to derive this formulation (cf Appendix A.2). These results validate the expression of the asymptotes given by (12), (13), and (14). Finally, the CPW approach gives a really good approximation of the spherical one, thus the CPW can be used to derive another key parameter that can be relevant for inversion: the frequency corresponding to the maximum of the log-amplitude spectrum, also called the rollover frequency. The latter is specifically studied in the following sub-section.

2.4 Maximum of the log-amplitude PSD

Different authors consider that the peak value of the log-amplitude spectrum, corresponding to its maximum value as a function of frequency, occurs for the Fresnel frequency f_F (Yeh & Liu, 1982; Forte & Radicella, 2002; Wernik et al., 2003; McCaffrey & Jayachandran, 2019). However, this is a hypothesis that deserves to be quantified. To do this, an exercise of estimating the difference between the frequency of the PSD maximum and the Fresnel frequency was performed. Note that the Fresnel frequency is approximately the first maximum of the filtering function for log-amplitude. Indeed, for the analysis, it is interesting to compute the value of $k_r=\sqrt{k_u^2+k_v^2}$ corresponding to the first null of the cardinal sine and cosine:

$$k_{\mathrm{sinc}}=\frac{2\pi }{\sqrt{\lambda R_{\mathrm{iono}}}},$$(15)

$$k_{\cos }=\frac{2\pi }{\sqrt{4\lambda L_v(1+\frac{R_{\mathrm{iono}}}{2L_v})}}.$$(16)

With the thin layer assumption R_iono ≪ L_v, it follows that k_cos ≪ k_sinc, so the first maximum value of the filtering functions corresponds approximately to the first −1 value of the cosine term, i.e., for:

$$k_m\simeq \frac{2\pi }{\sqrt{2}{r}_F^{\mathrm{PW},\mathrm{SW}}}.$$(17)

Finally, with the relationship between frequency and wave number, the Fresnel frequency can be expressed as:

$$f_F^{\mathrm{PW},\mathrm{SW}}=\frac{V_r}{\sqrt{2}{r}_F^{\mathrm{PW},\mathrm{SW}}}.$$(18)

For the simulations, two different configurations have been defined, corresponding to a LEO satellite moving over the polar or equatorial region, where ionospheric scintillation phenomena are frequent. The parameters describing these configurations are reported in Table 1. The main difference between polar and equatorial conditions is the orientation of the Earth’s magnetic field parametrized by angles (γ, ψ) and the anisotropy of the ionospheric irregularities, characterized by the anisotropic ratios (A_y, A_z). The drift velocity of the medium is set to 1000 m s⁻¹ for the polar configuration and 100 m s⁻¹ for the equatorial configuration. However, the velocity of the ionospheric piercing point must also be taken into account, as explained in (Morel et al., 2025). Naturally, this velocity depends on the satellite’s altitude. The other parameters have been kept equal between the two configurations, to facilitate the comparison. Note that the satellite altitude is only used for the SW formalism.

Figure 2 shows the relative error between the frequency of the maximum of the temporal log-amplitude spectrum (eq. (6)) estimated numerically and the Fresnel frequency (eq. (18)) for the polar (Fig. 2a) and equatorial (Fig. 2b) configurations. The study is proposed according to the main anisotropy A_z and the slope variations of PSD p_m, corresponding to the abscissa and ordinate axes, respectively, whereas the relative error is reported in a color scale. Because relative error is plotted, results are the same regardless of the approach used, PW or SW.

Figure 2 Relative error in percentage between the frequency of the maximum of the log-amplitude spectrum estimated numerically from (6) and the Fresnel frequency (18), $\frac{|f_{\max }-f_F|}{f_{\max }}$, in color scale for both polar (a) and equatorial (b) configurations according to the anisotropic ratio Az and the slope of the irregularities PSD pm.

Only for values of slope lower than 4 and high values of anisotropy, the Fresnel frequency seems to roughly approximate the frequency of the maximum f_max, though the relative error is still greater than 10% even in the most favorable cases.

So the statement that the Fresnel frequency corresponds to the frequency of the maximum of the log-amplitude spectrum is an approximation and does not really correspond to the frequency of the maximum of the log-amplitude spectrum in all configurations, as shown in Figure 2. A better expression for the frequency of the maximum has been researched here. Initially, we aim to find ω_max such that ${\left[\frac{\mathrm{d}{W}_{\chi }\left(\omega \right)}{\mathrm{d}\omega }\right]}_{\omega ={\omega }_{\max }}=0$, but since expression (6) makes an integral appear, obtaining an analytic expression for ω_max is extremely challenging. The key idea is therefore to first determine k_max such that ∇W_χ(k) = (0,0), where W_χ(k) = W_χ(k_u, k_v) is the 2D spatial spectrum given by:

$$W_{\chi }\left(k_u,k_v\right)=\pi r_e^2{\lambda }^2{R}_{\mathrm{iono}}{S}_{\Delta N_e}\left(k_u,k_v,0\right){\mathcal{F}}_{\chi }\left(k_u,k_v\right),$$(19)

where the demonstration for this derivation can be found in Wheelon (2003), Morel et al. (2024) for an incident plane wave and an incident spherical wave, respectively. The expression of the 2D spatial spectrum does not make an integral appear as the temporal spectrum defined by equation (6). It is therefore easier to calculate the spectrum derivative and its roots to estimate the maximum. The mathematical details are reported in Appendix B.1. Two methods are used, providing two different solutions:

$$k_{\max }=\{\begin{array}{cc}\sqrt{\frac{4K_0^2}{p_m-4}}& \mathrm{if} p_m>4+2.3\Delta \\ \sqrt{\frac{4{\pi }^2}{2r_{\mathrm{F}}^2}-\frac{4\pi p_m}{2r_{\mathrm{F}}^2\left(\pi +\Delta \right)}}& \mathrm{otherwise},\end{array}$$(20)

with $\Delta =\frac{K_0^2{r}_F^2}{2\pi }$, a_x and G do not appear here because we are in the particular case where A_y = A_x = 1 and so a_x = G = 1. The transformation from the spatial domain to the temporal domain is achieved using Taylor’s hypothesis, allowing wave numbers to be converted into frequencies. It is possible to find a relation between the temporal spectrum and the spatial one by using the following relation W_χ(k,ω) = δ(k · V_r − ω)W_χ(k), which can be derived from (Tatarskii, 1971).

$$W_{\chi }\left(\omega \right)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}W(\mathbf{k},\omega ){\mathrm{d}}^2\mathbf{k},$$(21)

$$={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}\delta (\mathbf{k}\cdot {\mathbf{V}}_{\mathbf{r}}-\omega )W_{\chi }\left(\mathbf{k}\right){\mathrm{d}}^2\mathbf{k},$$(22)

$$={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}\delta (k_u{V}_u+k_v{V}_v-\omega )W_{\chi }(k_u,k_v){\mathrm{d}}^2\mathbf{k},$$(23)

$$=\frac{1}{V_u}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}\delta \left(k_u-\frac{\omega -k_v{V}_v}{V_u}\right)W_{\chi }(k_u,k_v){\mathrm{d}}^2\mathbf{k},$$(24)

$$=\frac{1}{V_u}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}{W}_{\chi }\left(\frac{\omega -k_v{V}_v}{V_u},k_v\right)\mathrm{d}{k}_v.$$(25)

This last relationship can be verified by independently computing temporal (Eq. (6)) and spatial (Eq. (19)) spectra leading to equality (25). Let θ be the angle between the u-axis and the principal direction of the 2D spectrum, i.e., where the energy is mainly concentrated, k_max can be decomposed as: k_max = k_maxcos(θ)u+k_max sin(θ)v. Using now the dispersion relationship ω = 2πf = k_uV_u + k_vV_v, the frequency f_max corresponding to the wave number k_max can be derived as:

$$f_{\max }=\left|\frac{1}{2\pi }\left(k_{\max }\cos \right(\theta )V_u+k_{\max }\sin (\theta \left)V_v\right)\right|.$$(26)

However, it is important to keep in mind that this transformation is an approximation, as the frequency corresponding to the maximum of the temporal log-amplitude spectrum cannot be determined analytically due to the integral in equation (6). In particular cases where the Earth’s magnetic field H₀ (z axis) is aligned with u or v, or if one of the two velocity components is zero, the frequency of the maximum is:

$$f_{\max }=\{\begin{array}{cc}\frac{V_r}{2\pi }\sqrt{\frac{4K_0^2}{p_m-4}}& \mathrm{if} p_m>4+2.3\Delta \\ \frac{V_r}{2\pi }\sqrt{\frac{4{\pi }^2}{2r_{\mathrm{F}}^2}-\frac{4\pi p_m}{2r_{\mathrm{F}}^2\left(\pi +\Delta \right)}}& \mathrm{otherwise},\end{array}$$(27)

with V_r = V_u or V_v. An attempt to improve the accuracy was made by applying one iteration of the Newton–Raphson algorithm. The gain is not significant, and the obtained analytic expressions are more complicated, and thus would be difficult to use for inversion parameters. The results of this analysis are, however, reported in Appendix B.1, for information.

2.5 Assessments of log-amplitude PSD maximum modeling

To quantify the relevance of the proposed formulations for the log-amplitude PSD maximum, a specific study has been led.

Figure 3 shows the relative error between the frequency of the maximum of the temporal log-amplitude spectrum (computed numerically from Eq. (6)) and the frequency estimated with equations (27a) and (27b), for the polar (Fig. 3a) and equatorial (Fig. 3b) configurations. The relative error is represented versus the main anisotropy ratio A_z and the PSD slope p_m, reported in the abscissa and ordinate axes respectively. Let’s remark that, in theory, we deal with integral quantities, such as the asymptotic log-amplitude variance, which is valid for 2 < p_m < 6 (cf. Eq. (24) in Galiègue et al., 2016), or the derivation of the low-frequency asymptote of the log-amplitude PSD, valid for 1 < p_m < 5 (eq. (13)). These integrals are defined and finite for a restricted range of the spectral index p_m. In those cases, the real range of p_m is more restricted than that found in Figure 3. For the polar configuration (left graph), a region of high error is visible for low A_z values (below 5) and p_m values between 3 and 4. In this area, the error reaches nearly 100%, indicating a significant divergence between the frequency of the maximum estimated numerically from (6) and the estimated one with (27a) and (27b). From this area, when A_z or p_m increases, the error gradually decreases. For higher values of p_m and A_z, the error becomes very low, indicating a better match between the frequencies. For the equatorial configuration (right graph), the relative errors remain generally lower compared to the polar configuration. The error does not exceed 60%, even in the areas of higher error. The error remains relatively stable as A_z increases, while it gradually decreases with increasing p_m, following a similar trend to the polar configuration but with fewer extreme variations. The region where p_m is close to 5 (4.5 ≤ p_m ≤ 5.5) gives a high error of 60%, when A_z ≥ 11 for polar configuration and A_z ≥ 3 for equatorial configuration. In both configurations, the error is lower for higher values of A_z and p_m, suggesting that the frequency estimation is more accurate for these values. On the other hand, the polar configuration appears more sensitive to variations in A_z and p_m, showing much higher errors than the equatorial one for low values of A_z and p_m. Even if for some combinations of the parameters describing the ionospheric irregularities, the results lead to a non-negligible error, these graphs suggest that the equation applied to estimate the frequency maximum is more reliable than the Fresnel frequency commonly used and could serve in a PSD inversion objective.

Figure 3 Relative error in percentage between the frequency of the maximum of the log-amplitude spectrum estimated numerically from (6) and the estimated one from equations (27a) and (27b) in color scale for both polar (a) and equatorial (b) configurations according to the anisotropic ratio Az and the slope of the irregularities PSD pm.

A last error study is proposed by considering the relative error on the PSD maximum, versus the outer scale of turbulence L₀ and slope p_m, reported in abscissa and ordinate axes, respectively. As shown in Figure 4, the estimated frequency using relations (27a) and (27b) provides very good results for low PSD slope (p_m ≤ 3.5) and equatorial configuration, or for high PSD slope (p_m ≥ 5.5) and polar configuration. However, for slope values p_m close to 4 (for polar configuration) and close to 4.7 (for equatorial configuration), the relative error reaches up to 90% for large outer scale values.

Figure 4 Relative error in percentage between the frequency of the maximum of the log-amplitude spectrum estimated numerically from (6) and the estimated one from equations (27a) and (27b) in color scale for both polar (a) and equatorial (b) configurations according to the outer scale of turbulence L0 and the slope of the irregularities PSD pm.

Recalling that these results were obtained assuming A_y = 1 the tests were conducted for cases where 1 < A_y ≪ A_z showed similar performance.

3 Comparison of the different approaches

In Sections 2.1 and 2.2, different approaches to model the log-amplitude and phase PSD for a wave illuminating the ionosphere are introduced: the incident plane wave (PW), the incident spherical wave (SW) and a heuristic approach consisting of correcting the Fresnel length in the filtering functions: Corrected Plane Wave (CPW). These approaches are compared hereafter, and the sensitivity of the PSD to ionospheric parameters is also studied.

The comparison exercise allows studying the received signal PSD for an incident wave illuminating irregularities that are a plane wave or a spherical wave. Figures 5a and 5b show the log-amplitude and phase spectra computed from equations (6) and (7) for the three analytic approaches, plane wave (PW) (in blue), corrected plane wave (CPW) (in red and yellow crosses), and spherical wave (SW) (in red and yellow), in polar configuration defined by Table 1. The red and the yellow curves corresponds to a LEO satellite at an altitude of 600 km and a MEO satellite at an altitude of 20000 km respectively. The Corrected Plane Wave (CPW) approach gives a good approximation of the spherical one (SW) for both LEO and MEO satellites. The Fresnel frequencies $f_F^{\mathrm{PW}}$, $f_F^{\mathrm{MEO}}$, and $f_F^{\mathrm{LEO}}$ (Eq. (18)) are indicated with vertical black dashed lines on Figure 5a and do not match with the maximum of the log-amplitude spectrum. The frequencies of the maximum of the log-amplitude spectrum $f_{\max }^{\mathrm{PW}}$, $f_{\max }^{\mathrm{MEO}}$ and $f_{\max }^{\mathrm{LEO}}$ obtained with equations (27a) and (27b) are plotted in colored dashdotted vertical lines and give a more accurate approximation of the frequencies of the log-amplitude maximum PSD than the Fresnel frequencies. However, it can be observed that these values do not correspond perfectly to the maximum.

Figure 5 Temporal power spectral densities for different approaches.

Figure 5 allows to study the PW approximation by comparison to SW approach. The spectra for incident plane (PW) and incident spherical wave (SW) approaches present differences for a satellite in Low-Earth Orbit. The shape of log-amplitude and phase PSD is similar. However, a frequency and level shift can clearly be observed on both spectra. The difference is clearly observed in the example plotted in Figure 5a. The low frequency asymptote, the maximum value, and the frequency of the maximum are different between PW and LEO SW. Note that the spectra for the incident plane (PW) and the incident spherical wave (SW) approach converge as the satellite’s altitude increases, i.e. when the SW and PW Fresnel distances converge. From this finding, it is preferable to use the more general formalism of the incident spherical wave, in particular for satellites in Low-Earth Orbit. This result was expected and has also been demonstrated in the first part of this study (Morel et al., 2025).

To go further in the analysis, we tested other configurations whose results are plotted in Figure 6. Figures 6a and 6b show the results obtained for a polar and equatorial configurations respectively for different values of anisotropy A_z. Polar and equatorial simulation parameters are given in Table 1. In the graphs, the maximum of each log-amplitude PSD is indicated by a point in color. Figures 6c and 6d show the results considering the same configurations, but for different values of slope p_m. For the four graphs, the satellite velocity and so the ionospheric piercing point velocity is not taken into account. The impact of the satellite velocity would be to shift the PSD toward higher frequencies, without changing the proposed analysis results. First, we notice that for both configurations (Figs. 6a and 6b), expression (27b), which is $f_{\max }=\frac{V_r}{2\pi }\sqrt{\frac{4{\pi }^2}{2r_F^2}-\frac{4\pi p_m}{2r_F^2(\pi +\Delta )}}$, corresponds better to the frequency of the maximum than the Fresnel frequency, especially for high anisotropy values. Then, in Figures 6c and 6d, we notice that the highest values of the p_m slope (p_m = 5 or 7), the more the frequency of the maximum shifts towards the low frequencies. In this case, the frequency of the maximum could be approximated by expression (27a) which is $f_{\max }=\frac{V_r}{2\pi }\sqrt{\frac{4K_0^2}{p_m-4}}$. The opposite reasoning is also valid, when p_m decreases, the frequency of the maximum moves towards high frequencies, and the Fresnel frequency will then correspond better to the frequency of the maximum.

Figure 6 Log-amplitude PSD of SW approach, for two configurations, polar and equatorial, and different values of anisotropy Az and slope pm.

The high and low frequency asymptotes are also plotted in dashed lines in Figure 6. We notice an excellent match between the theoretical spectra and the high-frequency asymptotes for frequencies beyond 5 Hz. The analytic curves match the low-frequency asymptotes for frequencies lower than 0.1 Hz. These asymptotic elements are particularly relevant for extracting certain parameters characterizing the turbulent ionospheric layer. Indeed, we can observe the dependence of the low frequency asymptote to the anisotropy. Likewise, the peak frequency of the log-amplitude spectrum provides a simple dependence on the speed and the high frequency behavior gives direct access to the inertial slope of the spectrum of inhomogeneities. All these key parameters or formulas could be used in a numerical inversion process of the log-amplitude PSD to retrieve the parameters describing irregularities PSD.

4 Conclusion

This paper presents an analysis of the log-amplitude and phase PSD and their potential for estimating parameters that characterize the crossed turbulent ionospheric layer. The study is carried out using an analytic formalism under Rytov’s hypothesis, assuming a small turbulent layer thickness compared to the link range.

First, the developed integral formulas enable the investigation of errors potentially induced by using the incident plane wave model (PW) instead of the incident spherical wave model (SW) to predict 3D scintillation effects. By defining realistic propagation scenarios, comparisons of PW and SW approaches to predict scintillation effects are proposed. Initially, it is demonstrated that both approaches (PW and SW) yield similar PSD patterns, albeit with some differences in terms of level, low-frequency asymptote, and the maximum log-amplitude PSD position, particularly when the satellite altitude is low. Furthermore, the approximation consists of replacing the Fresnel distance in the PW formalism of PSD with the SW Fresnel distance is studied. Comparisons with the more rigorous spherical incident wave formalism, indicates that this approximation is highly accurate. Ultimately, this study illustrates the variation trends of the PSD as the satellite altitude changes, with typical including LEO, MEO, and GNSS’s typical satellite altitudes.

Additionally, the Rytov formalism allows deriving some asymptotic formulations of the log-amplitude and phase PSD at low and high frequencies. Rytov formalism also allows to study the maximum of the log-amplitude PSD, and several analytic formulations have been proposed. The latter yield satisfactory results in specific configurations with high anisotropy and low value of the slope of the irregularities PSD. The relative error compared to the numerically estimated frequency of the maximum of log-amplitude PSD is calculated to be less than 5% in favorable cases.

This expression of the frequency of the log-amplitude PSD maximum, along with the high- and low-frequency PSD asymptotes, constitute a set of analytic equations useful for addressing the inverse problem. This involves evaluating key parameters describing irregularities PSD, such as the anisotropy ratio along the Earth’s magnetic field, the slope of the irregularities PSD and the drift speed of the ionosphere.

These advancements lay the groundwork for an inversion method adapted to recover key parameters of ionospheric medium irregularities. They can be applied to GNSS signals affected by ionospheric scintillation, at high or low latitudes. Another potential application is Synthetic Aperture Radar measurements performed in the L-band (or at lower frequency), which provide images distorted by ionosphere irregularities.

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: Asmptotic formulations for log-amplitude and phase spectra

A.1 Log-amplitude and phase spectra high frequency asymptote

For ω tends to infinity, ${\lim }_{\omega \to \mathrm{\infty }}{\mathcal{F}}_{\chi }, {\mathcal{F}}_{\phi }=1$, equations (6) and (7) becomes:

$$\begin{array}{c}{W}_{\chi ,\phi }\left(\omega \right)=\frac{\pi r_e^2{\lambda }^2{R}_{\mathrm{iono}}{a}_x^{3-p}{A}_y{A}_z{C}_s}{V_u}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}{\left(A{\left(\frac{\omega -k_v{V}_v}{V_u}\right)}^2+B{k}_v^2+2C{k}_v\left(\frac{\omega -k_v{V}_v}{V_u}\right)+\frac{K_0^2}{a_x^2}\right)}^{-p_m/2}\mathrm{d}{k}_v,\\ =\frac{\pi r_e^2{\lambda }^2{R}_{\mathrm{iono}}{a}_x^{3-p}{A}_y{A}_z{C}_s}{V_u}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}(P{k}_v^2+Q{k}_v+R{)}^{-p_m/2}\mathrm{d}{k}_v,\end{array}$$

with:

$$P=A\frac{V_v^2}{V_u^2}+B-2C\frac{V_v}{V_u},$$

$$Q=2\mathrm{A\omega }\frac{V_v}{V_u^2}-2C\frac{\omega }{V_u},$$

$$R=A\frac{{\omega }^2}{V_u^2}+\frac{K_0^2}{a_x^2}.$$

The integral can be calculated if $P{k}_v^2+Q{k}_v+R>0$ which implies 4PR − Q² > 0, and if p_m > 1. By completing the square, by setting $T^2=\frac{R}{P}-\frac{Q^2}{4P^2}$, and using the change of variable $x=k_v+\frac{Q}{2P}$, the integral becomes:

$${\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}(P{k}_v^2+Q{k}_v+R{)}^{-p_m/2}\mathrm{d}{k}_v=2P^{-p_m/2}{\int }_0^{+\mathrm{\infty }}\mathrm{}(x^2+T^2{)}^{-p_m/2}\mathrm{d}x,$$(A.1)

which is equal to $P^{-p_m/2}{T}^{1-p_m}\sqrt{\pi }\frac{\mathrm{\Gamma }\left(\frac{p_m-1}{2}\right)}{\mathrm{\Gamma }(p_m/2)}$, after some calculations, the expression of equation (11) appears.

A.2 Log-amplitude spectrum low frequency asymptote

For ω tends to zero, K₀ is neglected, and simplified weighting functions are used, i.e., only the cosine term is considered (we consider the thin layer assumption i.e., $\frac{R_{\mathrm{iono}}}{2L_v}\ll 1$). For an incident plane wave, the simplified functions can be expressed as: ${\mathcal{F}}_{\chi }\left(\frac{\omega -k_v{V}_v}{V_u},k_v\right)=1-\cos \left(\frac{L_v}{k_0}\left(\frac{(\omega -k_v{V}_v{)}^2}{V_u^2}+k_v^2\right)\right)$ which is equal to $1-\cos \left(\frac{V_r^2{L}_v}{V_u^2{k}_0}{k}_v^2\right)$ when ω tends to zero.

$$W_{\chi }(\omega \to 0)=\frac{\pi r_e^2{\lambda }^2{R}_{\mathrm{iono}}{a}_x^{3-p}{A}_y{A}_z{C}_s}{V_u}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}{\left(A\frac{V_v^2}{V_u^2}{k}_v^2+B{k}_v^2-2C{k}_v^2\frac{V_v}{V_u}\right)}^{-\frac{p_m}{2}}(1-\cos (a{k}_v^2\left)\right)\mathrm{d}{k}_v.$$

Factoring kv and substituing $u=k_v^2$, the previous equation becomes:

$$W_{\chi }(\omega \to 0)=\frac{\pi r_e^2{\lambda }^2{R}_{\mathrm{iono}}{a}_x^{3-p_m}{A}_y{A}_z{C}_s}{V_u{\left(A\frac{V_v^2}{V_u^2}+B-2C\frac{V_v}{V_u}\right)}^{p_m/2}}{\int }_0^{+\mathrm{\infty }}\mathrm{}{u}^{-(p_m+1)/2}(1-\cos (\mathrm{au}\left)\right)\mathrm{d}u,$$

and ${\int }_0^{+\mathrm{\infty }}\mathrm{}{x}^{-\mu }(1-\cos (\mathrm{bx}\left)\right)\mathrm{d}x=\frac{b^{\mu -1}\pi }{2\mathrm{\Gamma }\left(\mu \right)\cos (\pi -\frac{\pi \mu }{2})}$, for 1 < μ < 3 and b > 0 which is equivalent to 1 < p_m < 5 and $\frac{V_r^2{L}_v}{V_u^2{k}_0}>0$. And so, after some calculations, the expression of equation (13) appears. To get the expression for the incident spherical wave model, just substitute the Fresnel radius for an incident plane wave into the Fresnel radius for incident spherical wave.

A.3 Phase spectrum low frequency asymptote

The low frequency asymptote for the phase spectrum can now be easily expressed with the log-amplitude asymptote.

$$\begin{array}{c}{W}_{\phi }(\omega \to 0)=\frac{\pi r_e^2{\lambda }^2{R}_{\mathrm{iono}}{a}_x^{3-p_m}{A}_y{A}_z{C}_s}{V_u}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}{\left(A\frac{V_v^2}{V_u^2}{k}_V^2+B{k}_v^2-2C{k}_v^2\frac{V_v}{V_u}+\frac{K_0^2}{a_x^2}\right)}^{-p_m/2}(2-F_{\chi }(...\left)\right)\mathrm{d}{k}_v\\ =\frac{2\pi r_e^2{\lambda }^2{R}_{\mathrm{iono}}{a}_x^{3-p_m}{A}_y{A}_z{C}_s}{V_u}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{}{\left(A\frac{V_v^2}{V_u^2}{k}_V^2+B{k}_v^2-2C{k}_v^2\frac{V_v}{V_u}+\frac{K_0^2}{a_x^2}\right)}^{-p_m/2}\mathrm{d}{k}_v-W_{\chi }\left(0\right)\\ =\frac{4\pi r_e^2{\lambda }^2{R}_{\mathrm{iono}}{a}_x^{3-p_m}{A}_y{A}_z{C}_s{K}_0^{-p_m}}{V_u{a}_x^{-p_m}}{\int }_0^{+\mathrm{\infty }}\mathrm{}(\xi k_v^2+1{)}^{-p_m/2}\mathrm{d}{k}_v-W_{\chi }(0)\end{array}$$

with $\xi =\left(A\frac{V_v^2}{V_u^2}+B-2\frac{V_v}{V_u}\right)\frac{a^2}{K_0^2}$, the integral is equal to $\frac{\mathrm{B}(1/2,p_m/2-1/2)}{2\sqrt{\xi }}$, with B being the beta function, And so, after some calculations, the expressions of equation (14) appears.

Appendix B: Maximum of log-amplitude PSD

B.1 Basic solutions

In the ionosphere, in particular at low latitudes, the irregularities are strongly anisotropic. The idea is to reduce the two-dimensional spectrum given by (19) to a one-dimension spectrum by determining the direction of the projection of the anisotropy on the (u, v) plane. This direction, where the energy is concentrated is given by $\theta =\frac{1}{2}\mathrm{atan}\left(\frac{2C}{A-B}\right)+n\frac{\pi }{2}$, with n an integer. Assuming that the anisotropy is predominant on the z axis, A_x, A_y ≪ A_z, in this case, it can be shown that θ ≈ −α_z. Then, we can set k_u = k_r cos(θ) and k_v = k_r sin(θ). For the filtering function, the thin layer approximation allows simplifying the expression (valid for PW and SW formalism):

$${\mathcal{F}}_{\chi }\left(k_r\right)\approx 1-\cos \left(\frac{k_r^2{r}_F^2}{2\pi }\right).$$(B.1)

The log-amplitude spectrum becomes:

$$W_{\chi }\left(k_r\right)=\pi r_e^2{\lambda }^2{R}_{\mathrm{iono}}{a}_x^{3-p_m}{A}_y{A}_z{C}_s{\left(G{k}_r^2+\frac{K_0^2}{a_x^2}\right)}^{-p_m/2}\left(1-\cos \left(\frac{k_r^2{r}_F^2}{2\pi }\right)\right),$$(B.2)

with $G=A{\cos }^2\left(\theta \right)+B{\sin }^2\left(\theta \right)+C\sin \left(2\theta \right)$. Using the previous assumptions on anisotropies i.e θ ≈ −α_z, the expression of G reduces to:

$$G=A_y^2{\cos }^2\left(\psi \right){\sin }^2({\alpha }_y-{\alpha }_z)+\left(\sin \right(\gamma \left)\sin \right(\psi )+\cos (\gamma \left)\cos \right(\psi \left)\cos \right({\alpha }_y-{\alpha }_z){)}^2.$$(B.3)

Then, using the relationship ${\alpha }_y={\alpha }_z+\arccos \left(\frac{\tan \left(\psi \right)}{\tan \left(\gamma \right)}\right)$, which implies 0 < |ψ| ≤ |γ| (cf Appendix of Galiègue et al., 2016), it can be shown that G reduces to:

$$G=A_y^2{\cos }^2\left(\psi \right)-\frac{{\sin }^2\left(\psi \right)}{{\sin }^2\left(\gamma \right)}\left(A_y^2{\cos }^2\right(\gamma )-1)$$(B.4)

Some particular cases give interesting values:

$$G=\{\begin{array}{lll}1,& \mathrm{if} \psi =\gamma \mathrm{or} A_y=1,& \\ A_y^2,& \mathrm{if} \psi =0.& \end{array}$$(B.5)

It follows from the differentiation of (B.2) that it is proportional to:

$$\frac{\mathrm{d}{W}_{\chi }}{\mathrm{d}{k}_r}\left(k_r\right)\propto k_r{\left(G{k}_r^2+\frac{K_0^2}{a_x^2}\right)}^{-p_m/2-1}\left[\frac{r_F^2}{\pi }\sin \left(\frac{k_r^2{r}_F^2}{2\pi }\right)\left(G{k}_r^2+\frac{K_0^2}{a_x^2}\right)+G{p}_m\left(\cos \left(\frac{k_r^2{r}_F^2}{2\pi }\right)-1\right)\right],$$(B.6)

which leads to the following equation system:

$$\{\begin{array}{ll}{k}_r& =0,\\ {\left(G{k}_r^2+\frac{K_0^2}{a_x^2}\right)}^{-p_m/2-1}& =0,\\ \frac{r_F^2}{\pi }\sin \left(\frac{k_r^2{r}_F^2}{2\pi }\right)\left(G{k}_r^2+\frac{K_0^2}{a_x^2}\right)+G{p}_m\left(\cos \left(\frac{k_r^2{r}_F^2}{2\pi }\right)-1\right)& =0.\end{array}$$

We are therefore interested in the third equation only, since the first is trivial and the second has no solution because G > 0. To solve the third equation, we use the trigonometric formula ${\sin }^2\left(x\right)=\frac{1-\cos \left(2x\right)}{2}$, set $X=\frac{k_r^2{r}_F^2}{2\pi }$ and $\Delta =\frac{K_0^2{r}_F^2}{2\pi a_x^{2}G}$. The equation becomes:

$$\sin \left(X\right)(X+\Delta )-p{\sin }^2\left(\frac{X}{2}\right)=0.$$(B.7)

By successively dividing the equation by $\sin \left(\frac{X}{2}\right)$ and $\cos \left(\frac{X}{2}\right)$, it comes:

$$\tan \left(Y\right)=\frac{4}{p}Y+\frac{2\Delta }{p},$$(B.8)

with X = 2Y. We cannot find an exact solution for this type of equation, but it is possible to either do an asymptotic expansion of the n-th root, the first root corresponding to the global maximum of the spectrum, or do a series expansion of the tangent. The expressions obtained with the first and the second method respectively gives:

$$Y^{\left(1\right)}\approx \frac{\pi }{2}-\frac{p_m}{2(\pi +\Delta )},$$(B.9)

$$Y^{\left(2\right)}\approx \frac{2\Delta }{p_m-4},$$(B.10)

and the corresponding wave numbers are:

$$k_{\max }^{\left(1\right)}=\sqrt{\frac{4{\pi }^2}{2r_F^2}-\frac{4\pi p_m}{2r_F^2(\pi +\Delta )}},$$(B.11)

$$k_{\max }^{\left(2\right)}=\sqrt{\frac{4K_0^2}{a_x^{2}G(p_m-4)}},$$(B.12)