Transionospheric Microwave Propagation: Higher-Order Effects up to 100 GHz

Abstract

Ionospheric refraction is considered as one of the major accuracy limiting factors in microwave space-based geodetic techniques such as the Global Positioning System (GPS), Satellite Laser Ranging (SLR), very-long-baseline interferometry (VLBI), Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS), and satellite altimetry. Similarly, a high-performance ground-to-space and space-to-ground microwave link is considered to be very important for synchronizing clocks in global networks. Moreover, precise time and frequency transfer may lead to new applications in navigation, Earth observation, solar system science, and telecommunications. However, all transionospheric microwave signals are subject to ionospheric refraction and subsequent delays in the travel time. Since the ionosphere is a dispersive medium for radio signals, the first-order propagation effect can be removed by combining signals at two or more frequencies. Anyway, higher-order ionospheric effects remain uncorrected in such combinations. The residuals can significantly affect the accuracy of precise positioning, navigation, as well as the performance of time and frequency transfer. Here, we studied ionospheric propagation effects including higher-order terms for microwave signals up to 100 GHz frequencies. The possible combination between the L, S, C, X, Ku, and Ka band frequencies is studied for the first-order ionosphere-free solutions. We estimated the higher-order propagation effects such as the second- and third-order terms and ray-path bending effects in the dual-frequency group delay and phase advance computation. Moreover, the correction formulas originally developed for global navigation satellite systems (GNSS) L-band frequencies are tested for mitigating residual errors at higher frequencies up to 100 GHz.

1. Introduction

The propagation of a radio wave through the ionospheric plasma can be described by the refractive index of the ionosphere given by the Appleton-Hartree formula [1]. At very high frequencies (>100 MHz), the refractive index depends mainly on the electron density and on the strength and the direction of the geomagnetic field in relation to the ray path. It becomes evident that the spatial distribution of the electron density along the ray path and the corresponding geomagnetic field relationships determine the ionospheric impact on the electromagnetic wave.

Since the ionosphere is a dispersive medium, radio wave propagation is frequency dependent. Therefore, by combining two signals, more than 99% of the ionospheric propagation delay can be corrected. However, higher-order propagation effects such as the second- and third-order terms in the refractive index and ray-path bending errors remain uncorrected in dual-frequency ionosphere-free combination. The range (or travel time) computation using trans-ionospheric signal is affected up to several centimeters (or 1.e−9 s) due to higher-order ionospheric terms. Therefore, higher-order effects should not be neglected in precise time and positioning applications, especially during times of enhanced total electron content (TEC).

Early work was done by Brunner and Gu [2] in computing higher-order ionospheric effects and developing correction formulas for them. Since then higher-order ionospheric effects have been studied by different authors during the last two decades, e.g., Bassiri and Hajj [3], Jakowski et al. [4], Strangeways and Ioannides [5], Kedar et al. [6], Fritsche et al. [7], Hawarey et al. [8], Hoque and Jakowski [10–16], Hernandez-Pajares et al. [17], Datta-Barua et al. [18], Morton et al. [19], etc.

Brunner and Gu [2] considered the second-order term and the curvature correction term for the dual-frequency ionospheric correction of the Global Positioning System (GPS) observations. Similarly, Bassiri and Hajj [3] studied the second-order ionospheric term assuming an Earth-centered tilted dipole approximation for the geomagnetic field. Jakowski et al. [4] studied ionosphere induced ray-path bending effects in precise satellite positioning systems such as GPS (20,000 km altitude) and Navy Navigation Satellite System (NNSS) signals (1000 km altitude). Strangeways and Ioannides [5] considered the ratio of the curved path lengths and the geometric path length and the ratio of the TECs along GPS L1 and L2 ray paths for determining the geometric distance of the Earth-GPS path.

Hawarey et al. [8] investigated the second-order ionospheric term for very-long-baseline interferometry (VLBI) and found that further improvements can be achieved using a more realistic magnetic field model such as the International Geomagnetic Reference Field (IGRF) model. Fritsche et al. [7] found that Global Navigation Satellite Systems (GNSS) satellite positions can also be improved to the centimeter level by applying higher-order ionospheric corrections.

Datta-Barua et al. [18] estimated the higher-order ionospheric errors using data from the federal aviation administration’s Wide Area Augmentation System (WAAS). They found that during ionospheric storms when slant range delays at GPS L1 can be as high as 100 m, the higher-order group errors in the GPS L1–L2 or L1–L5 dual-frequency combination can be tens of centimeters. Morton et al. [19] studied the second-order error analysis based on an extensive collection of electron density profiles measured by the Arecibo incoherent scatter radar and geomagnetic field vectors generated using the IGFR model.

Hoque and Jakowski [10–12] derived different approximation formulas to correct the second- and third-order terms, errors due to excess path length in addition to the free space path length, and TEC difference at two GNSS frequencies for L-band signals. Petrie et al. [21] investigated the potential effects of the bending terms on global GPS network.

Recently, Hernandez-Pajares et al. [17, 20] made a comprehensive summary of the second-order effect in receiver position and clock, tropospheric delay, geocenter offset, and GNSS satellite position and clock. They considered all the relevant higher-order contributions such as the second- and third-order terms, geometric bending, and slant total electron content bending (i.e., the difference between the slant total electron content (STEC) for straight and bent paths).

In the present paper, we have investigated the ionospheric propagation effects including higher-order terms for microwave signals up to 100 GHz frequencies. The possible combination between the L, S, C, X, Ku, and Ka (1–2, 2–4, 4–8, 8–12, 12–18, 27–40 GHz, respectively) band frequencies is studied for the first-order ionosphere-free solutions. We estimated the remaining higher-order propagation effects such as the second- and third-order terms and ray-path bending effects in the dual-frequency group delay and phase advance computation.

In our previous work, we developed correction formulas for mitigating the second-order ionospheric term [10, 11], the third-order term, and errors due to ray-path bending [4, 12, 15, 16] in the dual-frequency ionosphere-free linear combination. The correction formulas are developed and validated mainly for precise (e.g., centimeter and millimeter level) GNSS positioning at L-band frequencies. In the present paper, we investigated the performance of these correction formulas for mitigating remaining errors at higher frequencies up to 100 GHz.

2. First-order ionospheric term

The travel time τ of the signal for traveling the geometric distance between the transmitting satellite S and the receiver R can be written in units of seconds as (by dividing the true range expression given in [15] by the speed of light):

$\tau =\frac{L}{c}+\frac{1}{c}{\displaystyle \underset{S}{\overset{R}{\int }}\left(1-n\right)ds}-t_I^{len}$E1

where the optical distance $L={\displaystyle {\int }_S^{R}nds}$ is the line integral of the refractive index between the transmitting satellite and the receiver along the ray path, c is the velocity of light in a vacuum, $\frac{1}{c}{\displaystyle {\int }_S^R\left(1-n\right)ds}$ is the ionospheric time delay or advance, and $t_I^{len}$ is the excess time delay due to ionospheric ray-path bending in which ds is the ray-path element and n is either group or phase refractive index of the ionosphere. The Appleton-Hartree formula of the refractive index and its expansion for phase advance and group delay can be found in [3, 9, 15] and references therein. The ionospheric group delay t_Igr and phase advance t_I can be written in units of seconds as (by dividing ionospheric delay expressions given in [15] by the speed of light):

$t_{Igr}=\frac{1}{c}{\int }_S^R\left(n_{gr}-1\right)ds=\frac{p}{c{f}^2}\mp \frac{q}{c{f}^3}+\frac{u}{c{f}^4}$E2

$t_I=\frac{1}{c}{\int }_S^R\left(1-n\right)ds=\frac{p}{c{f}^2}\mp \frac{q}{2c{f}^3}+\frac{u}{3c{f}^4}$E3

where

$p=40.3\int n_{e}ds=40.3·TEC$E4

$q=2.2566\times {10}^{12}\int n_e\mathit{Bcos}\mathrm{\Theta }ds$E5

$u=2437\int n_e^{2}ds+4.74\times {10}^{22}\int n_e{B}^2\left(1+{\mathit{cos}}^2\mathrm{\Theta }\right)ds$E6

The upper (−) and lower (+) signs in the expressions (2) and (3) are related to the ordinary and extraordinary waves, respectively. The group refractive index n_gr is greater than unity causing travel time greater than the speed of light, whereas the phase refractive index n is less than unity causing travel time less than the speed of light.

The quantity n_e is the electron density, B is magnetic induction, and Θ is the angle between the Earth’s magnetic field vector and the propagation vector. The presence of the Earth’s magnetic field makes the ionosphere anisotropic that means the refractive index as well as the signal propagation depends on the propagation direction. The integral ${\displaystyle \int {\text{n}}_{\text{e}}\mathrm{ds}}$ along a signal path (i.e., curved path) is defined as the total electron content (TEC) and often measured in TEC units (1 TECU = 10¹⁶ electrons/m²). The parameters p, q, and u given by Eqs. (4)–(6) together with signal frequency f determine the first-, second-, and third-order ionospheric refraction effects, respectively.

The first-order ionospheric effect (group delay or phase advance) contributes above 99% of the total ionospheric effect and can be written as

$t_{Igr}^{\left(1\right)}=-t_I^{\left(1\right)}=\frac{40.3}{c{f}^2}\int n_{e}ds=\frac{40.3}{c{f}^2}TEC$E7

in units of seconds where c is in m/s, TEC is in electrons/m², and f is in Hz. Eq. (7) indicates that the travel time caused by the first-order term is equal in magnitude but opposite in sign for group delay and phase advance.

The first-order term is directly proportional to the total number of electrons TEC, encountered by the signal during its travel through the ionosphere and inversely proportional to the square of the signal frequency. So, if frequency and link-related TEC are known, the first-order propagation effect can be derived by Eq. (7). The first-order ionospheric delay on different frequencies is determined for 1 TECU and given in Table 1.

Frequency [GHz]	Ionospheric delay [s]	Ionospheric delay [m]
2	3.3607e−10	0.1007
8	2.1004e−11	0.0063
15	5.9745e−12	0.0018
30	1.4936e−12	4.4778e−04
60	3.7341e−13	1.1194e−04
100	1.3443e−13	4.0300e−05

Table 1.

Estimates of the first-order term at different frequencies for TEC = 1 TECU.

The vertical TEC may vary between 1 TECU and 300 TECU depending on a number of factors such as local time, geographic/geomagnetic location, season, solar activity level, etc. The frequency dependence of the first-order ionospheric term up to 100 GHz frequency at different levels of ionospheric ionization characterized by the vertical total electron content (VTEC) = 300, 150, 50, and 5 TECU is plotted for elevations 10°, 30°, and 60° in Figure 1. The vertical total electron content of 300 and 150 TECU corresponds to VTEC during extreme space weather conditions, 50 TECU corresponds to midlatitude day time, and 5 TECU corresponds to midlatitude night time VTEC. The following obliquity factor or mapping function is used to convert the vertical TEC to slant TEC [15]:

Figure 1.

Frequency dependence of the first-order term at different levels of ionospheric ionization and elevation angles. Left and right scales correspond to the units of seconds and meters, respectively.

$MF=\frac{\mathit{STEC}}{\mathit{VTEC}}=\frac{h_m+R_E}{\sqrt{{\left(h_m+R_E\right)}^2-{\left(R_h+R_E\right)}^2{\mathit{cos}}^2\beta }}$E8

where h_m is the maximum ionization height (e.g. 350 km), R_E is the Earth’s mean radius (~6371 km), R_h (~0) is the receiver height from the Earth’s surface, and β is the elevation angle.

Table 2 gives estimates of the first-order term at 2, 8, 15, 30, 60, and 100 GHz frequencies for 10° elevation angle and VTEC = 300 and 5 TEC units.

Frequency [GHz]	VTEC [TECU]
	300	5
	First term
	[s]	[s]
2	2.81e−07	4.69e−09
8	1.76e−08	2.93e−10
15	6.00e−09	8.33e−11
30	1.25e−09	2.08e−11
60	3.13e−10	5.21e−12
100	1.13e−10	1.88e−12

Table 2.

Estimates of the first-order term at different frequencies for VTEC = 300 and 5 TEC units and elevation = 10°.

As Figure 1 and Table 2 demonstrate, the first-order ionospheric delay can reach up to the 1.e−09 and 1.e−10 s levels at 30 and 60 GHz frequencies, respectively.

3. Higher-order terms

Although higher-order effects are less than 1% of the total ionospheric effects, they cannot be ignored in precise time and position applications. In precise applications dual-frequency measurements are commonly used to eliminate the major part (first-order effect) of the ionospheric propagation effects. The remaining higher-order terms can be up to 30 cm at L-band frequencies (15, 18). However, they will be much less at C, X, Ku, and K band frequencies. When analyzing ionospheric effects on signal propagation, it is a common practice to consider the code pseudo-range and carrier-phase observation equations for group delay and phase advance computation. Considering extraordinary wave of propagation, the code pseudo-range and carrier-phase expressions can be written in terms of ionospheric effects as (by dividing pseudo-range and carrier-phase expressions given in [12] by the speed of light):

$\mathrm{\Psi }=\mathrm{\tau }+t_{Igr}+t_I^{len}=\mathrm{\tau }+\frac{p}{c{f}^2}+\frac{q}{c{f}^3}+\frac{u}{c{f}^4}+t_I^{len}$E9

$\mathrm{\Phi }=\mathrm{\tau }-t_I+t_I^{len}=\mathrm{\tau }-\frac{p}{c{f}^2}-\frac{q}{2c{f}^3}-\frac{u}{c{f}^4}+t_I^{len}$E10

where τ is the travel time in vacuum. Although the measured travel times Ψ and Φ are biased by satellite and receiver clock errors, instrumental biases and tropospheric effect, and multipath effects, we have ignored them for simplicity in Eqs. (9) and (10). Due to the dispersive nature of the ionosphere, the propagation effect is frequency dependent, and we can eliminate the first-order term by combining signals at two different frequencies (f₁, f₂). The remaining terms can be written in units of seconds as (e.g., [12])

$\frac{f_1^2}{f_1^2-f_2^2}{\mathrm{\Psi }}_1-\frac{f_2^2}{f_1^2-f_2^2}{\mathrm{\Psi }}_2=\tau -\mathrm{\Delta }{t}_{TEC}-\mathrm{\Delta }{t}_2-\mathrm{\Delta }{t}_3-\mathrm{\Delta }{t}_{len}$E11

$\frac{f_1^2}{f_1^2-f_2^2}{\mathrm{\Phi }}_1-\frac{f_2^2}{f_1^2-f_2^2}{\mathrm{\Phi }}_2=\tau +\mathrm{\Delta }{t}_{TEC}+\frac{\mathrm{\Delta }{t}_2}{2}+\frac{\mathrm{\Delta }{t}_3}{3}-\mathrm{\Delta }{t}_{len}$E12

$\mathrm{\Delta }{t}_{TEC}=\frac{40.3\times \left({TEC}_2-{TEC}_1\right)}{c\left(f_1^2-f_2^2\right)}=\frac{40.3\times \left(\Delta {TEC}_{\mathit{bend}2}-\Delta {TEC}_{\mathit{bend}1}\right)}{c\left(f_1^2-f_2^2\right)}$E13

${TEC}_{1,2}=\int n_{e}ds=\left({TEC}_{LoS}+\Delta {TEC}_{\mathit{band}1,2}\right)$E14

$\Delta t_2=\frac{2.2566\times {10}^{12}\int n_e\mathit{Bcos}\mathrm{\Theta }ds}{c{f}_1{f}_2\left(f_1+f_2\right)}$E15

$\mathrm{\Delta }{t}_3=\frac{2437}{c{f}_1^2{f}_2^2}\int n_e^{2}ds+\frac{4.74\times {10}^{22}}{c{f}_1^2{f}_2^2}\int n_e{B}^2\left(1+{\mathit{cos}}^2\mathrm{\Theta }\right)ds$E16

$\mathrm{\Delta }{t}_{len}=\frac{t_2^{len}{f}_2^2-t_1^{len}{f}_1^2}{c\left(f_1^2-f_2^2\right)}$E17

where Ψ₁ and Ψ₂ are the measured pseudoranges and Φ₁ and Φ₁ are the measured carrier phases on f₁ and f₂. The terms Δt₂ and Δt₃ are the dual-frequency second- and third-order ionospheric terms. Due to the dispersive nature of the ionosphere, the bending effects as well as the total electron content (TEC) along the f₁ path will be different from that along the f₂ path. Considering this, the additional terms Δt_len and Δt_TEC are introduced in Eqs. (11) and (12) referring the dual-frequency residual errors due to the excess path and TEC difference, respectively. The quantities ΔTEC_bend1 and ΔTEC_bend2 are the excess TEC due to bending in addition to the straight line of sight (LoS) TEC, and $t_1^{len}$ and $t_2^{len}$ are the excess path length in addition to the LoS path length for f₁ and f₂ signals, respectively.

In the following sub sections, we estimated various dual-frequency residual terms for signal combination at different frequency bands considering different levels of ionospheric ionization.

3.1. Dual-frequency second-order term

The imposed Earth’s magnetic field causes an electron in the plasma to oscillate around the magnetic field line with the gyrofrequency fg = eB/(2πm) which is usually less than 1.4 MHz [22]. The value of B can be derived as ~5×10⁻⁵ Tesla for fg = 1.4 MHz and considered constant throughout the propagation. Therefore, assuming the worst case condition with B = 5×10⁻⁵ Tesla and Θ = 0, the maximum estimates of the dual-frequency second-order term (Eq. (15)) for group delay can be written as

$\Delta t_2=\frac{11.28\times {10}^7}{c{f}_1{f}_2\left(f_1+f_2\right)}\int n_{e}ds=\frac{11.28\times {10}^7}{c{f}_1{f}_2\left(f_1+f_2\right)}TEC$E18

The elevation angle dependence of the Δt₂ at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU is plotted in Figure 2.

Figure 2.

Elevation angle dependence of the Δt₂ at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU. Left and right scales correspond to the units of seconds and millimeters, respectively.

Table 3 gives the maximum estimates of the second-order residual term at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU and elevation = 10°, 30°, and 60°.

	300 TECU			5 TECU
	Elevation [deg]			Elevation [deg]
Freq-comb. (code) [GHz]	10	30	60	10	30	60
2–15	6.17e−12	3.88e−12	2.51e−12	1.03e−13	6.46e−14	4.19e−14
8–30	3.45e−13	2.17e−13	1.41e−13	5.75e−15	3.61e−15	2.34e−15
8–60	9.65e−14	6.06e−14	3.93e−14	1.61e−15	1.01e−15	6.55e−16

Table 3.

Maximum estimates of the second-order residual term in the dual-frequency group delay computation.

Figure 2 and Table 3 show that the dual-frequency second-order residual term is bigger than the 1.e−14 and 1.e−15 s levels for frequency combinations (2–15) and (8–30) GHz, respectively, even at VTEC = 5 TECU. Therefore, the second-order residual term cannot be ignored if 1.e−16/1.e−17 s level accuracy is required in the time transfer.

3.2. Dual-frequency third-order term

Eq. (16) indicates that the third-order residual term is proportional to the square of the electron density as well as geomagnetic induction B and angle between the geomagnetic field vector and propagation direction. The dependency on $n_e^2$ can be simplified assuming that the ionosphere is composed of Chapman layer function [23]. In this case, the so-called shape parameter $\eta =\frac{\int n_e^{2}ds}{N_m\times TEC}$ where Nm is the maximum ionization can be assumed to be 0.66 [2, 24], and thus the integral $\int n_e^{2}ds$ in Eq. (16) can be written as 0.66NmTEC. Therefore, assuming the worst case condition with B = 5×10⁻⁵ Tesla and Θ = 0, Eq. (16) can be simplified as

$\mathrm{\Delta }{t}_3=\frac{1608.4}{c{f}_1^2{f}_2^2}{N}_{m}TEC+\frac{23.7\times {10}^{13}}{c{f}_1^2{f}_2^2}TEC=\left(1608.4N_m+23.7\times {10}^{13}\right)\frac{TEC}{c{f}_1^2{f}_2^2}$E19

where Δt₃ is measured in seconds, c is in meters, TEC is in electrons/ m², f in Hz, and the maximum ionization Nm is measured in electrons/m³. The maximum ionization Nm can be estimated for a Chapman profile using the expression given by [12]

$\mathit{VTEC}=4.13H{N}_m$E20

where VTEC is the TEC in vertical direction and H is the atmospheric scale height. The VTEC will be measured in electrons/m² when H is in meters and Nm is in electrons/m³. The scale height H can be assumed as 70 km for a rough estimation of the third-order ionospheric term. The elevation angle dependence of the Δt₃ at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU is plotted in Figure 3.

Figure 3.

Elevation angle dependence of the Δt₃ at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU. Left and right scales correspond to the units of seconds and millimeters, respectively.

Table 4 gives the maximum estimates of the third-order residual term at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU and elevation = 10°, 30°, and 60°.

	300 TECU			5 TECU
	Elevation [deg]			Elevation [deg]
Freq-comb. (code) [GHz]	10	30	60	10	30	60
2–15	5.25e−13	3.30e−13	2.14e−13	2.66e−16	1.67e−16	1.08e−16
8–30	8.20e−15	5.15e−15	3.34e−15	4.16e−18	2.61e−18	1.69e−18
8–60	2.05e−15	1.29e−15	8.35e−16	1.04e−18	6.53e−19	4.24e−19

Table 4.

Estimates of the third-order residual term in the dual-frequency group delay in units of seconds.

Figure 3 and Table 4 show that the dual-frequency third-order residual term can be bigger than the 1.e−13 and 1.e−15 s levels for frequency combinations (2–15) and (8–30) GHz, respectively, during times of high TEC such as VTEC = 300 TECU.

3.3. Dual-frequency excess TEC term

Due to ray-path bending or curvature effect, the estimate of electron density integrated over a curved path is slightly greater that that along the straight line of sight propagation (TEC_LoS). The excess TEC along a curved path in addition to the TEC_LoS can be written as [12]

$\Delta {TEC}_{\mathit{bend}}=\frac{1.108\times {10}^{-3}\mathit{exp}\left(-2.1844\beta \right){TEC}^2}{f^{2}H{h}_m^{0.3}}$E21

where ΔTEC_bend is measured in TECU, atmospheric scale height H is in km, the maximum ionization height h_m is in km, the signal frequency f is in GHz, TEC is in TECU, and elevation angle β is in radians. Typical values for the parameters H = 70 km and h_m = 350 km are used in the present studies. Thus, knowing ΔTEC_bend at two frequencies, we can calculate Δt_TEC by Eq. (21). The elevation angle dependence of the Δt_TEC at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300, and 5 TECU is plotted in Figure 4.

Figure 4.

Elevation angle dependence of the Δt_TEC at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU. Left and right scales correspond to the units of seconds and millimeters, respectively.

Table 5 gives the maximum estimates of the residual excess TEC term at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU and elevation = 10°, 30°, and 60°.

	300 TECU			5 TECU
	Elevation [deg]			Elevation [deg]
Freq-comb.; (code) [GHz]	10	30	60	10	30	60
2–15	1.95e−12	3.59e−13	4.81e−14	5.42e−16	9.96e−17	1.34e−17
8–30	3.05e−14	5.60e−15	7.51e−16	8.47e−18	1.56e−18	2.09e−19
8–60	7.62e−15	1.40e−15	1.88e−16	2.12e−18	3.89e−19	5.22e−20

Table 5.

Estimates of the residual excess TEC term in the dual-frequency code combinations in units of seconds.

Figure 4 and Table 5 show that the dual-frequency residual excess TEC term can be bigger than the 1.e−12 and 1.e−14 s levels for frequency combinations (2–15) and (8–30) GHz, respectively, at 10° elevation angle during times of high TEC such as VTEC = 300 TECU.

Eq. (21) requires the knowledge of the ionospheric parameters H and hmF2. If actual parameters are not known, the formula may not be used in practical cases. Considering this, Hoque and Jakowski [16] derived the following correction approach considering excess TEC dependency only on TEC and elevation angle.

$\Delta {TEC}_{\mathit{bend}}=\frac{d_1}{f^2}\left(\frac{1}{\sqrt{1-d_2{\mathit{cos}}^2\beta }}-1\right){TEC}^2$E22

where d₁ = 1.4563 and d₂ = 0.8260. The coefficients are derived based on a nonlinear fit with ray-tracing results in least square senses. Although the correction given by the original approach Eq. (21) gives the best performance in comparison with the ray-tracing results, the correction results given by the new approach Eq. (22) is comparable.

3.4. Dual-frequency excess path term

Again the estimate of the total path length along a curved path is slightly larger than the LoS path length, and their difference is defined as the excess path length and denoted by $t_I^{len}$ in Eq. (17). To estimate the excess path length, we follow the estimates based on numerical ray-tracing computations given by [4]

$t_I^{len}=\frac{b_1}{c{f}^4}\left(\frac{1}{\sqrt{\left(1-b_2{\mathit{cos}}^2\beta \right)}}-1\right){TEC}^2$E23

where b₁ = 2.495×10⁸, b₂ = 0.859, and β is the elevation angle. The excess path length will be estimated in milliseconds when β is measured in radians, f is in MHz, TEC is in TECU, and c is in m/s. Now the dual-frequency residual excess path term Δt_len can be calculated by Eq. (17) in conjunction with Eq. (23). The elevation angle dependence of the Δt_len at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU is plotted in Figure 5.

Figure 5.

Elevation angle dependence of the Δt_len at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU. Left and right scales correspond to the units of seconds and millimeters, respectively.

Figure 5 and Table 6 show that the dual-frequency residual excess path term can be bigger than the 1.e−13 and 1.e−14 s levels for frequency combinations (2–15) and (8–30) GHz, respectively, at 10° elevation angle during times of high TEC such as VTEC = 300 TECU.

	300 TECU			5 TECU
	Elevation [deg]			Elevation [deg]
Freq-comb. (code) [GHz]	10	30	60	10	30	60
2–15	9.38e−13	1.73e−13	1.34e−14	2.61e−16	4.80e−17	3.83e−18
8–30	1.47e−14	2.70e−15	2.16e−16	4.07e−18	7.50e−19	5.99e−20
8–60	2.67e−15	4.92e−16	3.93e−17	7.42e−19	1.37e−19	1.09e−20

Table 6.

Estimates of the residual excess path term in the dual-frequency group delay computation in units of seconds.

Table 6 gives the estimates of the residual excess path term at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU and elevation = 10°, 30°, and 60°.

The estimation of excess path length can be improved if additional ionospheric parameters such as scale height H and peak density heights hmF2 are known. Such a formula is derived by [12].

$d_I^{len}=\frac{7.5\times {10}^{-5}\mathit{exp}\left(-2.13\beta \right){TEC}^2}{f^{4}H{\left({hmF}_2\right)}^{1/8}}$E24

where the excess path is measured in meters, TEC in TEC units, frequency f in GHz, atmospheric scale height H and height of maximum ionization hmF₂ in kilometers, and elevation angle β in radians.

3.5. Total dual-frequency residuals in group delay and phase advance

In the previous sections, we estimated different dual-frequency-residual terms separately. It is worthy to estimate their combined effects on group delay and phase advance computation. The elevation angle dependence in group delay (Δt_TEC + Δt₂ + Δt₃ + Δt_len, see Eq. 11) at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU is plotted in Figure 6.

Figure 6.

Elevation angle dependence of the total residuals (absolute values) in group delay computation at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU. Left and right scales correspond to the units of seconds and millimeters, respectively.

Table 7 gives the maximum estimates of the total dual-frequency residuals in the group delay at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU and elevation = 10°, 30°, and 60°.

	300 TECU			5 TECU
	Elevation [deg]			Elevation [deg]
Freq-comb.; (code) [GHz]	10	30	60	10	30	60
2–15	−9.59e−12	−4.74e−12	−2.79e−12	−1.04e−13	−6.49e−14	−4.20e−14
8–30	−3.99e−13	−2.30e−13	−1.45e−13	−5.77e−15	−3.62e−15	−2.35e−15
8–60	−1.09e−13	−6.37e−14	−4.03e−14	−1.61e−15	−1.01e−15	−6.55e−16

Table 7.

Estimates of the total residual error in the dual-frequency group delay in units of seconds.

Similarly the elevation angle dependence of the total dual-frequency residuals in the phase advance (Δt_TEC + Δt₂/2 + Δt₃/3 − Δt_len, see Eq. (12)) at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU is plotted in Figure 7.

Figure 7.

Elevation angle dependence of the total residuals in phase advance measurements at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU. Left and right scales correspond to the units of seconds and millimeters, respectively.

Table 8 gives the maximum estimates of the total dual-frequency residuals in the phase advance at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU and elevation = 10°, 30°, and 60°.

	300 TECU			5 TECU
	Elevation [deg]			Elevation [deg]
Freq-comb. (code) [GHz]	10	30	60	10	30	60
2–15	4.27e−12	2.23e−12	1.36e−12	5.18e−14	3.24e−14	2.10e−14
8–30	1.91e−13	1.13e−13	7.19e−14	2.88e−15	1.81e−15	1.17e−15
8–60	5.39e−14	3.16e−14	2.01e−14	8.06e−16	5.05e−16	3.28e−16

Table 8.

Estimates of the total residual term in the dual-frequency phase advance in units of seconds.

Figures 6 and 7 and Tables 7 and 8 show that the residual error in the dual-frequency group delay and phase advance cannot be ignored if 1.e−16/1.e−17 s level accuracy is required in the time transfer.

4. Higher-order term correction

In the following sections, we have discussed the possibility of higher-order propagation effects correction at very high frequencies.

4.1. Possibilities of second-order term correction

For the second-order term correction, the following correction formula is proposed by [11] for the global navigation satellite systems (GNSS) users in Europe with geographic latitude 30–65° N and longitude 15° W–45° E. Although the formula gives the best performance for the L-band signals, the second-order term correction at S (2–4 GHz), C (4–8 GHz), X (8–12 GHz), Ku (12–18 GHz), K (18–27 GHz), and Ka (27–40 GHz) band signals is also possible by this formula.

Single frequency (deduced from Eqs. 2 and 5):

${\mathrm{t}}_{\mathrm{I}\mathrm{g}\mathrm{r}}^{\left(2\right)}\mathrm{q}=\frac{2.2566\times {10}^{12}}{\mathrm{c}{\mathrm{f}}^3}\int {\mathrm{n}}_{\mathrm{e}}\mathrm{Bcos\Theta ds}=\frac{2.2566\times {10}^{12}}{\mathrm{c}{\mathrm{f}}^3}\overline{\mathrm{Bcos\Theta }}\times \mathrm{T}\mathrm{E}\mathrm{C}$E25

Dual frequency (deduced from Eq. 15):

$\Delta {\mathrm{t}}_2=\frac{2.2566\times {10}^{12}\times \overline{\mathrm{Bcos\Theta }}\times \mathrm{T}\mathrm{E}\mathrm{C}}{\mathrm{c}{\mathrm{f}}_1{\mathrm{f}}_2\left({\mathrm{f}}_1+{\mathrm{f}}_2\right)}$E26

In which

$\overline{\mathit{Bcos}\mathrm{\Theta }=-y_1\mathit{cos\alpha }+\left|\sqrt{r_1^2-y_1^2{\mathit{sin}}^2\alpha }\right|-2r_2\mathit{cos}\stackrel{}{\mathit{ά}}}$E27

where $\overline{\mathit{Bcos\Theta }}$ is the average value of the longitudinal component of the Earth’s magnetic field along the ray path. The parameters r₁, r₂, and y₁ are the functions of the receiver-to-satellite elevation angle β, geographic latitude ϕ, and longitude λ at the receiver position. The quantity α is the receiver-to-satellite azimuth angle and α’ is the modified azimuth angle. For details of formulation, we refer to [11].

To check the performance of the correction formula at very high frequency, e.g., 30 GHz, we estimated the second-order term by a 2D ray-tracing tool [12] and also by the correction formula Eqs. (25) and (27). The effect of the Earth’s magnetic field on the radio wave propagation is taken into account by considering the international geomagnetic field (IGRF) model in the ray-tracing tool. Figure 8 gives their comparisons.

Figure 8.

Estimates of second-order term by ray-tracing simulation and correction formula at 30 GHz for VTEC = 300 and 5 TECU (top and bottom plots, respectively) at elevation = 10°, 30°, and 60°.

We found that during times of high TEC such as VTEC = 300 TECU, the differences between the ray-tracing and correction results are in the level of 1.e−15 s.

4.2. Possibilities of third-order term correction

To correct the third-order ionospheric term, we proposed a correction formula based on analytical integration of the Chapman layer [12].

Single frequency (deduced from Eqs. 2 and 6):

$t_{Igr}^{\left(3\right)}=\frac{2487}{c{f}^4}\int n_e^{2}ds=\frac{1608}{c{f}^4}{N}_{m}TEC$E28

Dual frequency (deduced from Eqs. 15 and 36):

${\mathrm{\Delta }t}_3=\frac{1608}{c{f}_1^2{f}_2^2}{N}_{m}TEC$E29

The third-order term Δt₃ will be measured in meters when f is measured in Hertz and the maximum ionization N_m and TEC in electrons/m³ and electrons/m², respectively.

The above correction formula was validated for the GNSS L-band signals in former studies. Therefore, to check the performance of the correction formula at very high frequency, e.g., 30 GHz, we have estimated the third-order term by a 2D ray-tracing tool [12] and also by the correction formula Eq. (28). Figure 9 gives their comparisons.

Figure 9.

Estimates of the third-order term by ray-tracing simulation and correction formula at 30 GHz for VTEC = 300 and 5 TECU as a function of the elevation.

We see that during times of high TEC such as VTEC = 300 TECU, the differences between the ray-tracing and correction results are in the level of 1.e−18 s.

4.3. Possibilities of excess path term correction

To correct the excess path length, Jakowski et al. [4] proposed a correction formula based on simulation studies of the Chapman layers. The formula is already given by Eq. (23) for single frequency applications.

The correction formula was mainly derived for the GNSS L-band signals. Therefore, to check the performance of the correction formula at very high frequency, e.g., 30 GHz, we estimated the excess path length by the 2D ray-tracing tool [12] and also by the correction formula Eq. (23). Figure 10 gives their comparisons.

Figure 10.

Estimates of the excess path term by ray-tracing simulation and correction formula at 30 GHz for VTEC = 300 TECU as a function of the elevation.

We found that during times of high TEC such as VTEC = 300 TECU, the differences between the ray-tracing and correction results are in the level of 1.e−17 s.

5. Conclusions

From our simulation studies, we have found that the time delay due to the first-order ionospheric term is in the range of 1.e−9–5.e−12 s depending on the ionospheric ionization and elevation angles at 30–60 GHz frequencies. The corresponding second- and third-order terms are in the ranges of 1.e−13–2.e−16 s and, respectively. The corresponding excess path length error due to bending is in the range of 1.e−15–2.e−20 s.

We found that for the second-order ionospheric term, the differences between the ray-tracing and correction results are in the order of 1.e−15 s for a microwave link at 30 GHz when assuming high vertical total electron content (VTEC) such as VTEC = 300 TECU. The corresponding differences are at the level of 1.e−18 and 1.e−17 s for the third-order and ray-path bending terms, respectively.

European Space Agency (ESA) Atomic Clock Ensemble in Space (ACES) mission on board International Space Station (2014–2016) is aimed to transfer clock signals from low earth orbiting (LEO) satellites to the Earth using a time and frequency transfer link in the microwave domain supporting the clock performance at the level of 10⁻¹⁶ s [25]. High-performance frequency comparison of optical clocks between space and ground is essential for supporting future missions in the field of fundamental physics. Such a metrology link must provide frequency and time comparison and dissemination with an uncertainty level of 10⁻¹⁸ and beyond. Our investigation shows that for achieving such an accuracy level in the time and frequency transfer using trans-ionospheric microwave links, the higher-order ionospheric propagation effects must be corrected for.