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

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 altim-etry. 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 naviga- tion, Earth observation, solar system science, and telecommunications. However, all transionospheric microwave signals are subject to ionospheric refraction and subse- quent 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 computa- tion. 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. 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.


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 transionospheric 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).
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][11][12] derived different approximation formulas to correct the secondand 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 secondorder 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.

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): 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, Þds is the ionospheric time delay or advance, and t len I 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): where p 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 ð n e 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 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.
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.
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.

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 t len 1 and t len 2 are the excess path length in addition to the LoS path length for 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.

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

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 2 e can be simplified assuming that the ionosphere is composed of Chapman layer function [23]. In this case, the so-called shape parameter Nm · TEC where Nm is the maximum ionization can be assumed to be 0.66 [2,24], and thus the integral Ð n 2 e ds in Eq. (16) can be written as 0.66NmTEC. Therefore, assuming the worst case condition with B = 5 + 10 −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] VTEC 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°.  Table 3. Maximum estimates of the second-order residual term in the dual-frequency group delay computation.
Wave Propagation Concepts for Near-Future Telecommunication Systems 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.   Table 4. Estimates of the third-order residual term in the dual-frequency group delay in units of seconds.

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] ΔTEC bend ¼ 1:108 · 10 À3 exp À2:1844β 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)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(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.

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 len I in 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 Eq. (17). To estimate the excess path length, we follow the estimates based on numerical raytracing computations given by [4] t len 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.

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.   Table 7. Estimates of the total residual error in the dual-frequency group delay in units of seconds.  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°.

Higher-order term correction
In the following sections, we have discussed the possibility of higher-order propagation effects correction at very high frequencies.
Single frequency (deduced from Eqs. 2 and 5):  Dual frequency (deduced from Eq. 15): In which 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°. where BcosΘ is the average value of the longitudinal component of the Earth's magnetic field along the ray path. The parameters 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. 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.

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.

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.

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