## 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):

where the optical distance *c* is the velocity of light in a vacuum, *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):

where

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 ^{16} electrons/m^{2}). 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

in units of seconds where *c* is in m/s, TEC is in electrons/m^{2}, 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 |

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]:

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 |

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):

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_{1}, f_{2}). The remaining terms can be written in units of seconds as (e.g., [12])

where Ψ_{1} and Ψ_{2} are the measured pseudoranges and Φ_{1} and Φ_{1} are the measured carrier phases on *f*_{1} and *f*_{2}. The terms Δ*t*_{2} and Δ*t*_{3} 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*_{1} path will be different from that along the *f*_{2} 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 *f*_{1} and *f*_{2} 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^{−5} Tesla for fg = 1.4 MHz and considered constant throughout the propagation. Therefore, assuming the worst case condition with *B* = 5×10^{−5} Tesla and Θ = 0, the maximum estimates of the dual-frequency second-order term (Eq. (15)) for group delay can be written as

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

**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°.

**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 ^{−5} Tesla and Θ = 0, Eq. (16) can be simplified as

where Δt_{3} is measured in seconds, c is in meters, TEC is in electrons/ m^{2}, f in Hz, and the maximum ionization Nm is measured in electrons/m^{3}. The maximum ionization Nm can be estimated for a Chapman profile using the expression given by [12]

where VTEC is the TEC in vertical direction and H is the atmospheric scale height. The VTEC will be measured in electrons/m^{2} when H is in meters and Nm is in electrons/m^{3}. 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*_{3} at pair of frequencies (2, 15), (8, 30), and (8, 60) GHz for VTEC = 300 and 5 TECU is plotted in **Figure 3**.

**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°.

**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]

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**.

**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°.

**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.

where d_{1} = 1.4563 and d_{2} = 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

where b_{1} = 2.495×10^{8}, b_{2} = 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** 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.

**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].

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*_{2} 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*_{2} + Δ*t*_{3} + Δ*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**.

**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°.

Similarly the elevation angle dependence of the total dual-frequency residuals in the phase advance (Δ*t*_{TEC} + Δ*t*_{2}/2 + Δ*t*_{3}/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**.

**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°.

**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):

Dual frequency (deduced from Eq. 15):

In which

where *r*_{1}, *r*_{2}, and *y*_{1} 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.

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):

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

The third-order term Δ*t*_{3} will be measured in meters when *f* is measured in Hertz and the maximum ionization *N*_{m} and TEC in electrons/m^{3} and electrons/m^{2}, 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.

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.

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^{−16} 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^{−18} 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.