On the Electron Temperature in the Topside Ionosphere as Seen by Swarm Satellites, Incoherent Scatter Radars, and the International Reference Ionosphere Model

: The global statistical median behavior of the electron temperature ( T e ) in the topside ionosphere was investigated through in-situ data collected by Langmuir Probes on-board the European Space Agency Swarm satellites constellation from the beginning of 2014 to the end of 2020. This is the ﬁrst time that such an analysis, based on such a large time window, has been carried out globally, encompassing more than half a solar cycle, from the activity peak of 2014 to the minimum of 2020. The results show that Swarm data can help in understanding the main features of T e in the topside ionosphere in a way never achieved before. T e data measured by Swarm satellites were also compared to data modeled by the empirical climatological International Reference Ionosphere (IRI) model and data measured by Jicamarca (12.0 ◦ S, 76.8 ◦ W), Arecibo (18.2 ◦ N, 66.4 ◦ W), and Millstone Hill (42.6 ◦ N, 71.5 ◦ W) Incoherent Scatter Radars (ISRs). Moreover, the correction of Swarm T e data recently proposed by Lomidze was applied and evaluated. These analyses were performed for two main reasons: (1) to understand how the IRI model deviates from the measurements; and (2) to test the reliability of the Swarm dataset as a new possible dataset to be included in the underlying empirical dataset layer of the IRI model. The results show that the application of the Lomidze correction improved the agreement with ISR data above all at mid latitudes and during daytime, and it was effective in reducing the mismatch between Swarm and IRI T e values. This suggests that future developments of the IRI T e model should include the Swarm dataset with the Lomidze correction. However, the existence of a quasi-linear relation between measured and modeled T e values was well veriﬁed only below about 2200 K, while for higher values it was completely lost. This is an important result that IRI T e model developers should properly consider when using the Swarm dataset.


Introduction
The ionosphere is a plasma medium consisting of ions and electrons with also a neutral component that, for most of the cases, are not in thermal equilibrium. In fact, since the early satellites and rockets missions in the 1950s and 1960s [1], a bulk of evidence has been accumulated demonstrating the non-thermal equilibrium between ions, electrons, and neutrals in the ionosphere. Plasma temperature (ions plus electrons) usually exceeds the neutral one by a factor strongly dependent on altitude, time, and location; moreover, large differences among ion and electron temperatures are also observed. As a consequence, when talking about the temperature of the ionosphere, we need to differentiate between ion (T i ), electron (T e ), and neutral (T n ) temperature.
The non-thermal equilibrium of the ionosphere has a number of physical and chemical implications for the ionospheric dynamics and composition. For example, the temperature of charged particles deeply influences the distribution of electron density by driving the recombination rate in the bottomside ionosphere, and affects the diffusive equilibrium state in the topside ionosphere [2,3]. The main source of energy heating the ionosphere is the Sun: either directly, through the EUV and X-ray radiations illuminating the Earth, or indirectly, at high latitudes, through the precipitation of charged particles induced by the interactions between the solar wind and the magnetosphere-ionosphere system and the consequent Joule heating [1,[4][5][6]. In the ionization processes occurring in the ionosphere, a certain amount of energy is released as an excess of kinetic energy of photoelectrons. Such photoelectrons dissipate their energy through collisions with charged and neutral particles in a cascade process that first raises the temperature of the electrons (because the transfer of energy via elastic collisions is more efficient between particles of similar mass). In a second stage, the established temperature difference drives the transfer of energy between electrons and the surrounding medium. Above the F2 region, most of the energy transfer is from electrons to ions because ion-electron Coulomb collisions are more efficient than neutral-electron ones; then, the temperature of the ions also increases above the neutrals one. As a further step, the ions lose their excess energy through collisions with neutrals and other ions, thereby raising the temperature of neutrals. As a result of these complex interactions, the entire photoelectron energy dissipates in the neutrals and it is finally transferred down in the lower layers of the atmosphere by conduction and radiation.
Among the three mentioned temperatures, T e is the most variable because of the very low heat capacity of electrons [4]; moreover, electrons are the most affected by the solar radiation energy input and by the interaction with charged particles entering the Earth's ionosphere following the geomagnetic field lines at auroral latitudes. The close coupling with the geomagnetic field lines and the plasmasphere produces very distinct spatial patterns in the T e distribution [6]. As a consequence, T e exhibits very distinctive diurnal, seasonal, spatial, and solar activity trends. This large variability makes it difficult to obtain a precise and reliable description of T e through in-situ or remote sensing measurements. In fact, in-situ measurements taken by Langmuir Probes (LPs) and Retarding Potential Analyzers (RPAs) [1,7] on-board Low-Earth-Orbit (LEO) satellites provide a global spatial description for different local times, but they are highly dependent on the satellite orbital configuration and altitude and are limited to the satellite mission lifetime. Conversely, ground-based remote sensing techniques such as Incoherent Scatter Radars (ISRs) [7,8] can better describe the altitudinal, diurnal, and seasonal variability but for a fixed geographic location and for the very few existing ISR facilities. Moreover, the solar activity variability can be described when long time series become available. Furthermore, the cost of this measurement technique strongly limits the number and distribution of ISRs worldwide. As a consequence, for modeling purposes, both in-situ satellites and remote sensing ISRs observations are very important and complement each other.
An outstanding example is represented by the International Reference Ionosphere (IRI) model [9], which is the empirical climatological model of the ionosphere used as reference by the ionospheric community, and recently recognized as the official ISO (International Organization for Standardization) standard for the ionosphere [10]. Through the exploitation of different measurement datasets of ionospheric parameters collected by ionosondes, ISRs, LEO satellites, Global Positioning System (GPS) radio occultation, and rockets, IRI describes the hourly monthly median behavior of the electron density, electron and ion temperature, and ion composition in the ionospheric altitude range, on a global basis, for different levels of solar and magnetic activity. Due to its empirical nature, IRI evolves and improves when new datasets become available and are used to derive the empirical coefficients of its analytical formulations. In this framework, the ability to model the changes of the ionospheric ground state as a function of the different climatological conditions of both the interplanetary medium and the Earth's seasons is also fundamental in relation to space weather science.
Despite a long series of measurements acquired by both satellites and ISRs, the major current gap in modeling and studying ionospheric T e consists in the lack of an accurate description of its dependence on the most important geographic and geophysical parameters. There is still insufficient knowledge of the global distribution and relevant description of T e longitudinal variations as well as of its seasonal and hemispheric asymmetries, despite some attempts [11][12][13]. A difference between long-term and short-term T e variations with F10.7 has also been noted [14], which further complicates the modeling of T e . Empirical models still struggle to describe small-scale structures such as the morning overshoot [15][16][17] and its dependence on solar and magnetic activity or the mid-latitude T e enhancement [18,19]. The dependence on the magnetic activity and the choice of an appropriate magnetic activity proxy for T e modeling is still an open problem; for example, Brace and Theis [20] found no clear dependence of T e on the ap index, with exception of perhaps a slight systematic increase of T e with ap at polar latitudes. Apart from modeling, there are still gaps in the fundamental theory. For instance, Pavlov et al. [21] introduced an anomalous conductivity to obtain an agreement between theory and EXOS-D data in the plasmasphere. Even very advanced physics-based models still encounter problems with reproducing observational data (e.g., the post-sunset T e enhancement) and must be driven by measured data [22]. Another tricky point concerns the existing problems associated with the frequent T e differences among satellite data from different missions, and between in-situ data and ground-based ISR measurements [23,24], which require in many cases a cross calibration and/or a correction of these data. Because of this, studying and global modeling of T e has become a very challenging problem.
From this point of view, T e observations collected by LPs on-board the European Space Agency (ESA) Swarm satellites [25] are very promising for a future inclusion in the IRI T e model [26]. In fact, the three Swarm satellites launched at the end of 2013 and still operating provide T e data in the topside ionosphere region on a global basis, for different local times, at a 2-Hz rate. As a consequence, the present dataset of Swarm LP T e measurements is that large to cover different latitudes, local times, seasons, and also solar activity levels. Contrary to the electron density, T e data have not been extensively exploited so far; only a few and rather narrowly focused studies have been published, such as [17,24]. In this context, we aimed to investigate the main climatological trends of T e in the topside ionosphere by taking advantage of seven full years (from the beginning of 2014 to the end of 2020) of Swarm T e observations by comparing the corresponding spatial, diurnal, and seasonal trends with those modeled by IRI. Recently, to improve the accuracy of Swarm T e data, Lomidze et al. [24] proposed a correction to observations, based on the comparison with three ISRs located at different latitudes. In this study, we applied the Lomidze correction to Swarm data to verify its performance when compared against IRI-modeled T e data and T e data observed by Jicamarca (12. The paper is organized as follows. In Section 2, the Swarm and ISRs T e datasets, along with the IRI T e model and the Lomidze correction, are discussed. In Section 3, an overall comparison between Swarm-measured (original and after applying the Lomidze correction) and IRI-modeled T e values is first given based on the entire dataset composed by the three Swarm satellites; in addition, Swarm B and IRI T e data are climatologically compared with the corresponding values recorded by the Jicamarca, Arecibo, and Millstone Hill ISRs. In the same section, we also describe and discuss the geographic, diurnal, and seasonal trends shown by both Swarm B and IRI data. Section 4 is devoted to the final discussion of our results and is supported by the conditioned probability density function analysis of Swarm-measured and IRI-modeled data. The conclusions and the future developments are the subject of Section 5.

ESA Swarm In-Situ Electron Temperature Data and Application of the Lomidze Correction
The ESA Swarm mission is a constellation of three LEO satellites launched at the end of 2013, and still operating, with the aim of studying the geomagnetic field, the electric currents in the magnetosphere and ionosphere, and the impact of the solar wind on the dynamics of the upper atmosphere [25]. The three Swarm satellites are labeled A, B, and C, and share the same design and on-board instrumentation. All three satellites were put in a circular near-polar orbit. Swarm A and C have the same orbit configuration (inclination of 87.35 • , initial altitude of about 460 km, east-west separation of about 1-1.5 • in longitude), while Swarm B has a different one (inclination of 87.75 • , initial altitude of about 510 km). As a consequence, Swarm A and C fly in tandem, while Swarm B moves away from the couple A and C by covering different local times. Each satellite takes about 130-140 days to cover all the local times.
Among the instruments carried by Swarm satellites, we focused on LPs data [27]. LPs provide measurements of in-situ electron density (N e ) and T e along the satellite orbit with an original sampling frequency of 2 Hz in the harmonic mode [24,28]. Moreover, 1-Hz data are also provided by interpolation of the 2-Hz dataset; it is this dataset containing time series of in-situ N e and T e data at 1-Hz rate that is used in this work. Specifically, LP data collected by Swarm A, B, and C from 1 January 2014 to 31 December 2020 (seven full years) are taken into account. Swarm's data are freely downloadable at ftp://swarm-diss.eo.esa.int (accessed on 24 September 2021).
Swarm's LP data are provided with three flags [28]: Flags_LP contains information on the source of LP measurements, and Flags_Ne and Flags_Te characterize N e and T e measurements, respectively. Detailed information about Swarm data quality is available at https:// earth.esa.int/web/guest/swarm/data-access/quality-of-swarm-l1b-l2cat2-products (accessed on 24 September 2021). In this study, only the most reliable data were considered; specifically, those recorded by LPs in High-Gain mode. LP data were then selected with the following choice of the flags: Flags_LP = 1, Flags_Ne ≤ 29, and Flags_Te = 10 or 20.
Lomidze et al. [24] developed a correction of both N e and T e Swarm LP data based on the analysis of Swarm measurements matching coincident values retrieved by ISRs, COSMIC/FORMOSAT-3 GPS radio occultation, and ionosondes. For what concerns the Swarm's T e data correction, Lomidze et al. [24] separately considered the measurements collected when the LPs are in High-Gain mode (HG in the following) from those collected with the LPs in Low-Gain mode. Since only HG data were considered in this study, we applied only the HG Lomidze correction to T e data. As a consequence, the Equations (4)- (6) in Lomidze et al. [24] were applied on the Swarm A, B, and C T e datasets, namely: T e,Lom = 1.2815 · T e,LP − 1167 + 7.293 · N e,LP /10 4 Swarm A, HG T e,Lom = 1.2248 · T e,LP − 1047 + 8.548 · N e,LP /10 4 Swarm B, HG T e,Lom = 1.1334 · T e,LP − 762 + 4.088 · N e,LP /10 4 Swarm C, HG In Equations (1)-(3) T e,Lom is the Lomidze-corrected T e , T e,LP is the original Swarmmeasured value, N e,LP is the corresponding N e value measured by Swarm LP (i.e., unadjusted). In Equations (1)-(3), and throughout the paper, T e values are in K and N e values are in cm −3 . Equation (1) is valid for the dataset of HG Swarm A LP T e dataset, Equations (2) and (3) for the corresponding Swarm B and C datasets, respectively.
As outlined by Equations (1)-(3), the Lomidze correction is based on a regression analysis of Swarm-measured T e,LP with data retrieved by the three ISRs located at Jicamarca, Millstone Hill, and Arecibo from December 2013 to June 2016 (for more information on the dataset selection on which Equations (1)-(3) are based see [24]). The shortness of the dataset did not allow a characterization of possible diurnal, seasonal, geographic, and solar activity dependences of the corrections introduced in Equations (1)-(3). As a consequence, Equations (1)-(3) were applied on the entire Swarm datasets although some residual diurnal, seasonal, geographic, and solar activity dependences were somewhat expected. The application of these ISRs-based corrections allowed greatly improving the accuracy of Swarm T e data when compared to ISRs ones (stepping from median differences between 300 and 400 K to near zero) without affecting the corresponding precision and correlation [24]. As outlined by Lomidze et al. [24], the inclusion of the term N e,LP in Equations (1)-(3) produced a worsening in the precision of the correction (less than 60 K) but with the advantage of removing some certain non-physical patterns from Swarm T e data.
In this work, both Swarm T e,LP original data and Lomidze-corrected T e,Lom data were employed. In Section 3, a statistical comparison between T e data measured by Swarm, calibrated after applying the Lomidze correction, and the corresponding ones modeled by IRI is presented. In the same section, a statistical comparison with T e data recorded by Jicamarca, Arecibo, and Millstone Hill ISRs is also performed.

Electron Temperature Observations by Incoherent Scatter Radars
ISRs exploit the Thomson backscatter from ionospheric electrons to retrieve the value of several ionospheric parameters including T e [8]. Specifically, T e profiles are recorded by ISRs in the ionospheric altitude range from about 100 to 1000 km.
T e observations collected by Jicamarca, Arecibo, and Millstone Hill ISRs were used in this work for validation purposes. Jicamarca ISR is located right above the magnetic equator at a Quasi-Dipole (QD, [29]) latitude of 0.2 • N, Arecibo (QD latitude 27.0 • N) is a low-latitude station, while Millstone Hill (QD latitude 51.8 • N) is a mid-latitude/subauroral station. Jicamarca data are from 1996 to 2020, Arecibo ones from 1974 to 2015, and Millstone Hill ones from 1976 to 2020. ISRs data are not continuous in time but the length of the datasets here used allow for a quite uniform representation of the different diurnal and seasonal trends.
The ISRs T e data (including corresponding T e error) were downloaded from the http://cedar.openmadrigal.org website (accessed on 24 September 2021). The altitudes of Swarm satellites fall into standard remote sensing altitude ranges of all three ISRs, which provide reliable data for decades. A similar dataset referring to T i , for Millstone Hill was recently employed in Pignalberi et al. [30].
To perform a comparison with Swarm B-measured and IRI-modeled values, we selected T e values recorded by ISRs in the range 510 ± 20 km. Only the most reliable ISRs' data were selected, i.e., those with a percentage error lower than 10%. This filtering mostly affected data recorded during the night and early-morning, which are those characterized by the largest dispersion.

Electron Temperature Description by the IRI Model
In the last version of IRI (IRI-2016, [9]), two different models for T e are implemented; they are conventionally named BIL-1995, and TBT-2012 (http://irimodel.org/, accessed on 24 September 2021). The BIL-1995 is the oldest option and is based on the works by Bilitza [31][32][33] who, through the application of the Booker approach [34] to the Brace and Theis [20] and Spenner and Pluge [35] models, describes the T e vertical profile for the ionospheric altitude range between 60 and 2000 km. Since the IRI-2001 version [36], IRI has included an additional option based on Intercosmos data [37], that since the IRI-2012 version [38] was updated with the most recent TBT-2012 model [26] which incorporates all available satellite T e data. The solar activity dependence, also described in Truhlik et al. [26], is mainly based on results from Truhlik et al. [39]. In the current IRI version, this additional option is called TBT-2012+SA and represents the option related to the TBT-2012 model with the solar activity term switched on. Since the TBT-2012+SA is the current default option for the T e modeling in IRI, it has been applied in this work to perform a comparison with Swarm and ISRs measured data and will be briefly described here.
In the TBT-2012 model [26], T e values at five fixed anchor points at 350, 550, 850, 1400, and 2000 km of altitude are modeled through a spherical harmonic expansion in a system of associated Legendre polynomials (up to the 8th order) in terms of the magnetic local time and a latitudinal coordinate based on a combination of invariant and dip latitude (invdip, [40]). Through the use of T e measurements from LPs on-board several satellite missions (for a detailed description of the dataset see Bilitza et al. [23] and Truhlik et al. [26]), the coefficients of the spherical harmonic expansions were calculated for summer and winter solstices, and for the combined equinoxes (independently of the hemisphere); then, interhemispheric differences in magnetic coordinates were not considered in the model. The solar activity dependence was then included by selecting three solar activity ranges and describing the corresponding variability as a function of the PF 10.7 solar index [26].
The complete vertical T e profile was obtained by applying the Booker profile function formalism [33,34]. First, a piece-wise function was built making a linear interpolation of modeled T e values between anchor points, which produced a vertical profile divided in regions with a constant T e gradient. Second, Epstein-step functions were applied to realize the transition from regions with different gradients and obtain the skeleton profile, which produced a continuous analytical representation of the T e vertical profile. Finally, the skeleton profile was integrated to obtain the Booker profile function of T e , which output a continuous and smooth description of the entire modeled T e vertical profile. For a detailed mathematical description of the Booker profile function formalism, refer to the Appendix section of Pignalberi et al. [30].
In IRI, T e is constrained to be equal to T n described by the NRLMSISE-00 model [41] for altitudes lower than 120 km. An example of IRI (TBT-2012+SA default option) modeled vertical T e profile is given in Figure 1 (green curve) where the anchor points that IRI uses to describe the entire profile are highlighted. An anchor point at a non-fixed altitude around 270 km (T e,m in Figure 1) was introduced on the basis of Jicamarca and Arecibo ISRs data to describe a local daytime maximum in T e characterizing mostly the low latitudes [33]. In this work, the IRI model (TBT-2012+SA option) was run on the entire Swarm satellites dataset as if it were collocated and synchronized with LPs on-board Swarm satellites. As a consequence, 1-Hz time series of IRI-modeled T e values were produced for the time and locations probed by Swarm satellites. In this way, a one-to-one comparison between Swarm-measured and IRI-modeled T e values was guaranteed.

Overall Statistical Comparison between Measured and Modeled Electron Temperature Values
First of all, we aimed to investigate the effect of the application of the Lomidze correction on Swarm T e measurements through a comparison with IRI-modeled values for the period 2014-2020 for all the three Swarm satellites.
An overall comparison between Swarm-measured (original and Lomidze-corrected) and IRI-modeled T e values is given in Figures 2 and 3 as joint histograms and distributions of residuals between measured and modeled data. From the residuals between measured and modeled values, the following statistical metrics were calculated: mean and standard deviation of the residuals, Root Mean Square Error (RMSE, Equation (4)), and Relative Root Mean Square Error (RRMSE, Equation (5)).
In Equations (4) and (5), T e,measured is the T e measured by Swarm (i.e., T e,LP ) or corrected according to Lomidze et al. [24] (i.e., T e,Lom ), while T e,modeled is the corresponding value modeled by IRI. Temperatures are in K. N is the number of values on which the statistics are calculated. Figure 2 represents the joint histograms between T e values measured by Swarm and modeled by IRI and gives information on the distribution of measured and modeled data as a function of the corresponding magnitude. Most of the values cluster in the reddish region, i.e., where more than 10 4 values are counted inside each bin (the scale is logarithmic). It is interesting to note that IRI-modeled data never go beyond about 4300 K, while Swarm data (both original and Lomidze-corrected ones) show a long tail well beyond 5000 K (not shown in the figure), even though the application of the Lomidze correction reduces their magnitude. To characterize the accordance between modeled and measured values and the global effect of the application of the Lomidze correction, corresponding residuals were studied. Figure 3 represents the statistical distributions of the residuals between values measured by Swarm and modeled by IRI. By looking at the statistical metrics reported in the upper left corner of each panel, it emerges how the application of the Lomidze correction to Swarm data improved the accuracy when compared to IRI-modeled data. In fact, the mean of the residuals decreased from about 560-650 K to about 150-250 K, with a 400 K average improvement. Conversely, the application of the Lomidze correction did not affect the precision of Swarm data, as highlighted by the standard deviation values that remained in the range 380-470 K. Nevertheless, the application of the Lomidze correction generally improved the agreement with IRI data, as demonstrated by the RMSE and RRMSE values that decreased, respectively, from about 690-780 K to about 440-530 K, and from about 27-29% to about 20-23%.  [24]. From top to bottom, the analysis refers to Swarm A, B, and C. In each panel, the number of total counts is reported in the lower right corner.
To sum up, the results shown in Figure 3 demonstrate how the application of the Lomidze correction to Swarm T e data can improve the agreement with IRI-modeled values, without impacting the corresponding precision. However, we are aware that the comparison with a model, IRI, cannot guarantee the reliability and efficacy of the Lomidze correction; the only comparison that can truly tell whether this correction improves real T e observations is that with other independent measurements. This is why we performed a comparison against T e data observed by ISRs.

Statistical Comparison against ISRs Data
We implemented a statistical comparison among Swarm (original and corrected values), IRI, and ISRs T e data for the Jicamarca, Arecibo, and Millstone Hill locations to verify the ability of Swarm and IRI to describe the diurnal and seasonal trends of T e shown by ISRs data, and their accuracy.
Unlike Lomidze et al. [24], we did not look for correspondences (in both time and space) between ISRs and Swarm data because of the limited number of correspondences achievable (of the order of some hundreds from 2014 to 2020, for each satellite). Instead, we applied a statistical approach by comparing the diurnal and seasonal trends exhibited by both T e values collected by ISRs and by Swarm B for the ISRs locations. In this way, we could investigate how much Swarm, and also the IRI model, can describe the diurnal and seasonal T e trends at equatorial, low, and mid latitudes by making a comparison with Jicamarca, Arecibo, and Millstone Hill data, respectively. This comparison and those of the following section were limited to the Swarm B satellite because it is the nearest to the IRI vertical profile anchor point at 550 km of altitude. Moreover, as follows from Figures 2 and 3, the statistical analysis for all three satellites was very similar, with slightly larger differences for Swarm B.
Specifically, T e values collected at Jicamarca, Arecibo, and Millstone Hill in the range (510 ± 20) km were binned as a function of the magnetic local time (MLT, in QD coordinates) and season, and corresponding boxplots were generated. The binning in 15-minutes wide (in MLT) bins allowed describing the diurnal trend of T e at ISRs location. The seasonal trend was studied by binning data in four periods (expressed in day of the year, doy) around solstices and equinoxes and covering the entire year: • No solar activity data selection was performed because the quantity of both ISRs and Swarm data does not allow further split, and because the investigation of T e dependence on solar activity is outside of the scope of this study. As a consequence, the results in    In summary, every ISR dataset was split in 384 (96 diurnal x 4 seasonal) bins and the following statistical metrics were calculated for the boxplot representation: • 5th percentile, i.e., the lower whisker; • 25th percentile, i.e., the first quartile; • 50th percentile, i.e., the second quartile, representative of the median; • 75th percentile i.e., the third quartile; • 95th percentile, i.e., the upper whisker.
For this kind of comparison, the median is the more appropriate statistical metric because it cuts off the outliers of the distribution. The difference between 75th and 25th percentiles defines the inter-quartile range (IQR), which provides a measure of the data dispersion around the median value, while the whiskers highlight the tails of the distribution. Figures 4-6 show the boxplots for the Jicamarca (Figure 4), Arecibo ( Figure 5), and Millstone Hill (Figure 6) sites. For the comparison with ISRs results, both Swarm Bmeasured (original and corrected with Lomidze) and IRI-modeled T e data were selected and binned as done for ISRs data. Specifically, from the Swarm B dataset, we selected those measurements falling inside a region centered at ISRs' geographic coordinates and extending ±5 • in latitude and ±20 • in longitude. We tested regions of different extension in both latitude and longitude around the ISRs' location. The choice we made is the one balancing the amount of data needed for the statistical comparison and the preservation of the characteristic trends at the ISRs locations. Then, the selected Swarm B (original and corrected with Lomidze) and IRI data were binned as ISRs data and the corresponding median values are superposed to ISRs' boxplots in Figures 4-6. Figure 4 shows how equatorial values at Jicamarca are characterized by the presence of a remarkable morning peak between 5 and 8 MLT, with median values from ISR exceeding 3000 K; this phenomenon is often referred as morning overshoot [15][16][17]42]. Overall, daytime values stand at about 1000-1300 K, while nighttime ones at about 700-1000 K, in a quite steady fashion. Differently from the dawn solar terminator hours, the ones related to the dusk solar terminator do not exhibit a peak but only a small ramp with T e smoothly stepping from the daytime to the nighttime behavior. Of course, the MLTs of both the morning overshoot and the dusk ramp slightly change according to the season. The MLTs characterized by the morning overshoot exhibit the highest dispersion of data as described by the IQRs and whiskers values. The seasonal variability is really weak at Jicamarca due to the illumination conditions that do not change significantly over the year. Both Swarm and IRI can catch the main diurnal features of T e for each season but with remarkable differences in the absolute values. IRI represents reliably nighttime values, while it slightly overestimates daytime values and underestimates values in the morning sector when the overshoot occurs. Differently, Swarm (both original and corrected with Lomidze) performs a general overestimation independently of the hour. Moreover, Swarm values do not exhibit differences between day and night, with the exception of the morning peak. The application of the Lomidze correction reduces the original overestimation of Swarm measurements, but not enough to make them match with ISR data. When comparing with ISRs data, we have to take into account that the IRI TBT-2012+SA option is based on data recorded from satellites [26], that are often limited in both space and time. In fact, the differences between IRI and Jicamarca ISR values might be related to the limited spatial and temporal resolution of the model (20 degrees in latitude and 80 minutes in MLT), resulting in the inability of the model to accurately describe the equatorial daytime dip in T e and thus resulting in an overestimation of daytime T e at the magnetic equator by IRI. At MLTs of the narrow and high morning peak, this limited resolution is most likely at the base of both the underestimation and some mismatches of the position of the morning peak maximum made by the IRI model. Another reason that could explain the observed differences is the fact that in-situ probes can produce higher T e values in a non-equilibrium plasma [43] or if some sort of probe contamination occurs [44]. A possible combination of both effects thus could result in both an overestimation of about 30% during daytime and an underestimation of about 20% of the peak maximum made by the IRI model (as evidenced by Figure 4). During nighttime, when plasma is in thermal equilibrium and the limited spatial and temporal resolution of the model is not so important, the agreement between IRI and ISR is very good. Figure 5 shows the results for Arecibo. Compared to Jicamarca, at Arecibo the seasonal variability is more marked above all for MLT hours from noon to dusk. Like Jicamarca, Arecibo values also show a remarkable morning increase with temperatures abruptly increasing from about 1000 to 2500-3000 K. However, Arecibo daytime values do not decrease as for Jicamarca but remain at around 2000 K (except for the December solstice) with significant differences in the sector 14-19 MLT, depending on the season. In fact, March and September equinoxes show a very similar trend where the first hump corresponding to the morning increase is followed by a small second hump before sunset, while for the June solstice (summer), the second hump does not appear and values steadily decrease from daytime to nighttime. The December solstice (winter) is instead characterized by quite constant values (at about 2500 K like the morning peak one) during daytime and a very distinct decrease at sunset. IRI can reliably represent the nighttime T e at Arecibo, while some departures are evident during daytime, above all in the morning for the June solstice (underestimation) and in the central hours of the day at equinoxes (overestimation). Compared to Jicamarca, at Arecibo, the accordance between Swarm and ISR data is better, especially in the daytime sector, while nighttime values are still affected by a remarkable overestimation (even after applying the Lomidze correction). A very good accordance is found for the winter season (December solstice) during daytime. Swarm and IRI can describe both the diurnal and seasonal trends of T e at Arecibo. Figure 6 summarizes the results for Millstone Hill, where T e manifests a very distinct diurnal behavior characterized by a sharp division between nighttime and daytime values. Nighttime values are of the order of 1100-1500 K, while the daytime values range between 2500 and 3000 K. Differently from Arecibo, in this case, daytime values are always quite constant with a slight tendency to maximize at noon during summer (June solstice). Winter (December solstice) nighttime values are remarkably higher than summer ones due to conjugate heating [23,45]. Except for June solstice, IRI can reliably describe the diurnal behavior of T e at Millstone Hill. It is the same for Swarm after applying the Lomidze correction; in this case the agreement with ISR's data is remarkable for both the representation of the diurnal and seasonal trends and the accuracy of the absolute values; only a slight overestimation is visible during nighttime particularly at December solstice.
It has to be considered that the Lomidze correction, aside from being based on the same three ISRs here studied, is heavily biased toward Millstone Hill ISR data that were more represented in the Lomidze et al. [24] dataset than Arecibo data and above all Jicamarca. As a consequence, since Lomidze et al. [24] cumulated the data of the three ISRs to derive a single calibration curve, it should be expected that the correction is more effective at Millstone Hill than at Arecibo and Jicamarca. This agrees very well with the results shown in Figures 4-6. Moreover, the Lomidze correction does not discriminate between daytime and nighttime or between different seasons. Our results show that the Lomidze correction has the best performance for daytime hours and the worst for nighttime hours, with no distinct seasonal difference. This suggests that a latitudinal and day/night dependence should be included in the Lomidze correction when enough data become available in the future.
As we can see from Figures 4-6, the daily T e profile varies with the ISR geographic location. For example, at Jicamarca, which is at the magnetic equator, T e shows a single sharp maximum around dawn that can be associated with the well-known morning overshoot. This feature is characteristic of low latitudes, and in fact is present also in the daily profiles acquired at Arecibo. Here, and at Millstone Hill (middle and sub-auroral latitudes), the increase of T e in the dayside is due to the increased sunlight and a second hump at dusk marks the decrease of T e in passing from day to night. At Millstone Hill, an increased T e during nighttime is recorded at the December solstice, mainly due to the decrease of plasma density in the F layer and the following reduction of the collisional cooling effect, and also to the conjugate heating. To better highlight and describe all of these different spatial features exhibited by the topside T e , in the next section, the main spatial trends of T e will also be investigated, in both geographic and magnetic coordinates.

Statistical Trends of Swarm-Measured Electron Temperature Values and Comparison with IRI-Modeled Ones
The results of Section 3.1 and 3.2 demonstrate that the application of the Lomidze correction improves the agreement between Swarm data and those modeled by IRI and measured by different ISRs. This is why, in the following analyses, we focus only on Swarm B data corrected according to Lomidze et al. [24].
Swarm B data from 2014 to 2020, and corresponding IRI-modeled values, were binned as a function of the geographic coordinates with bins 2.5 • -wide in latitude and 5 • -wide in longitude. For each bin, the median value was considered, as representative of the main statistical trend, after discarding the outliers through the Median Absolute Deviation statistical procedure (see Section 3.3 by Pignalberi et al. [46] for a detailed description). The binning was first performed on the entire dataset and then by selecting MLT (in QD coordinates) sectors representative of the dawn (5-7 MLT), daytime (12)(13)(14)(15), and nighttime (0-3 MLT) conditions. This was done because at a fixed altitude, the T e diurnal variations are known to be the most important along with the latitudinal ones [20,26,39]. Figure 7 shows the corresponding results considering all the MLTs and the different MLT sectors, for both Swarm-measured and IRI-modeled data, and the corresponding percentage of the normalized residuals between them. The figure highlights how the T e spatial distribution follows the configuration of the Earth's magnetic field lines. In fact, the characteristic shape of the magnetic equator and of the auroral regions stand out, and this is somewhat expected because the plasma in the topside ionosphere is generally constrained by the geomagnetic field [47]. This behavior has also been recently confirmed by Giannattasio et al. [48,49] that, by using Swarm data, identified patterns in the parallel electrical conductivity (which is crucially affected by T e ) related to features of the geomagnetic field such as regions R1 and R2 [50] or the magnetic cusp [51]. They also identified variations in the parallel electrical conductivity due to MLT, QD latitude, season, solar, and geomagnetic activity. The MLTs around dawn highlight the fact that the morning peak is a phenomenon affecting only the low latitudes, which follows the magnetic equator shape and reinforces in the longitudinal sector between 80 • W and 80 • E (including South America and Africa) where the departure between the magnetic equator and the geographic one is the largest. At the morning peak, Swarm B values exceed 3000 K in many longitudinal sectors, well above those recorded at mid latitudes and even higher than high latitude ones. In contrast, during daytime and nighttime, the high latitudes experience the highest values; in both cases, the increase of T e from mid to high latitudes is limited to a quite narrow latitudinal range where the spatial gradient is concentrated: at mid latitudes during daytime and at the sub-auroral and auroral boundaries during nighttime. When compared to IRI-modeled values, the largest differences emerge at magnetic equator latitudes at dawn, and at low and mid latitudes during daytime and nighttime, respectively. Overall, Swarm measurements exceed IRI values for most cases and this is in line with the results of Figures 2 and 3. The excess is particularly clear at the morning peak and at low and mid latitudes at nighttime. If we recall the analysis shown in the previous section, based on the comparison with ISRs data (Figures 4-6), we can say that Swarm overestimates the nighttime values; at the morning peak, Jicamarca ISR results (Figure 4) point out that the difference between Swarm and IRI is mainly due to a joint action of both the Swarm overestimation and the IRI underestimation. As highlighted by Figure 7, the Earth's magnetic field has a strong influence on the diurnal trend of Swarm and IRI values. Therefore, we binned T e data also as a function of QD latitude and MLT with bins 2.5 • -wide and 15 minutes-wide, respectively. The binning was first performed on the entire dataset and then by selecting different days of the year representative of the March equinox, June solstice, September equinox, and December solstice conditions, as done in Section 3.2. for ISRs data. Figure 8 shows that the overall diurnal trend of T e is characterized by different well-defined features at different latitudes. At low latitudes, the morning peak is the main feature and small day-night differences are visible, which confirms the results from Jicamarca ISR (Figure 4). At mid latitudes (30-60 • QD latitude), no distinct morning peak is present and the main feature is the remarkable difference of absolute values between daytime and nighttime, which confirms the results from Millstone Hill ISR (Figure 6). At high latitudes, the diurnal trend is less pronounced with high values throughout the day that intensify during daytime at the auroral boundaries. The agreement between Swarm and IRI is very good for what concerns the overall diurnal pattern. The main differences characterize the morning peak hours and the nighttime hours for low and mid latitudes. It should be noted that IRI cannot describe the very well-located morning peak maximum just over the magnetic equator, and the T e increase at the auroral boundaries during daytime. For both cases the reason has to be found in the limited spatial resolution of the model. As expected, no remarkable differences characterize the equinoxes, while the comparison between the solstices highlights seasonal differences and a different description of these by Swarm and IRI. In fact, while IRI shows a clear seasonal dependence for daytime hours, with higher T e values in the winter hemisphere than in the summer one, Swarm data show a less marked seasonal difference. The diurnal behavior of T e shown in Figures 4-8 is also consistent with previous results in the literature, where the observed trend was recognized to be mainly due to solar UV radiation (in the dayside) and the precipitation of particles (both in the dayside and nightside) caused by magnetosphere-ionosphere coupling [52]. Apart from the obvious contribution due to the excess of energy of photoelectrons in the dayside especially at mid latitudes, electrons, protons and alpha particles of solar origin may be injected into the ionosphere during open magnetosphere conditions, i.e., when reconnection at the magnetopause occurs, at QD latitudes of about 80 • and around 12:00 MLT [53,54]. The precipitation of such particles is expected to enhance T e in the cusp regions [51] and at lower latitudes in the pre-noon sector [55]. T e also increases on the nightside due to particle precipitation. In more detail, plasma acceleration occurs in response to reconnection phenomena in the far-away geomagnetic tail regions, generating a three-dimensional current system modeled by the current wedge paradigm [56] that connects the tail to the nightside ionosphere in correspondence of the R2 region [50]. Interestingly, the peak of such enhancement is observed in correspondence with the boundary between R2 and the main ionospheric trough in the pre-dawn sector due to the joint action of particle precipitation and the decrease of N e [49,57].
The diurnal variation of T e also shows a clear dependence on the local season, which emerges from Figure 8 when considering the six-months shift of local seasons in Northern and Southern hemispheres in correspondence with solstices and equinoxes. Specifically, the increase of T e is remarkable above all in the nightside winter, from about 20 MLT to predawn at QD latitudes of about 60 • and especially in the Southern hemisphere, and in the dayside summer. Such a behavior can be easily interpreted in light of the previous literature. In fact, previous studies concluded that field-aligned currents density changes from winter to summer by about a factor of 2 in the dayside due to the increase of photoelectron flux and the following increase of T e [58][59][60][61][62]. In this case, the increase of T e is driven by both the increase of EUV photoionization and, at very high latitudes (~80 • ), by precipitations in the cusp region around noon. In the nightside, particle precipitation is the dominant process at auroral and sub-auroral latitudes and becomes more and more important when proceeding from summer to winter. In particular, during winter the enhancement of T e reaches the maximum expansion from dusk to about 08:00 MLT, especially in the Southern hemisphere. The explanation of this behavior is the simultaneous collapse of both electron and ion density in the ionospheric F layer, and the following collapse of the contribution of collisional cooling to T e . Therefore, in the absence of mechanisms that reduce the energy of precipitating particles T e is free to increase due, e.g., to energy exchanges with the nightside magnetosphere. Not surprisingly, there is an anticorrelation between N e and T e , and T e is significantly enhanced in regions of depleted N e , such as between region R2 and the main trough in the nightside winter; while, in contrast, it decreases in regions of enhanced N e due to energy loss to the ions [19,48,49,57,63]. Some studies in the past tried to establish a link between N e and T e , since N e is a major factor in electron energy loss due to the ionospheric plasma quasi-neutrality (density of ions is nearly equal to N e ). Generally, an anti-correlation was found. Brace and Theis [20] established an interrelationship between N e and T e using Atmosphere Explorer C data. Kakinami et al. [64] found a U shape in the T e -N e dependence using Hinotori satellite data. Their results suggest that an additional heat source has to be considered, even though the mechanism is unknown. Su et al. [65] dealt with the investigation of T e -N e correlation using data from the DEMETER satellite, by highlighting the variability of this correlation depending on the latitude, local time, and season. Hu et al. [66] used N e measured from COSMIC as a driver for T e in their neural network T e model. Generally, they found a better agreement between their model and ISR data than between ISR data and IRI values. Thus, it is evident that the T e -N e relationship could help in the investigation of T e behavior in the ionosphere but it is a complex problem itself. This is why, in this study, we did not investigate the T e -N e correlation using the Swarm data and we focus only on T e . Compared to the spatial and diurnal trends, the interpretation of the seasonal trends exhibited by Swarm data is more difficult and the differences with IRI are significant. Overall, the main seasonal differences concern the summer and winter seasons for specific MLT sectors. Nighttime T e values at high latitudes have been found to be deeply influenced by the solar illumination conditions [67]; in fact, high values are found in the summer hemisphere lit by the Sun. In contrast, low latitudes are not affected by seasonal changes at nighttime. T e at mid latitudes during local winter can be affected by strong heating due to conjugate photoelectrons if the opposite hemisphere is already sunlit.

Discussion
The results of Section 3.1. highlighted how the application of the Lomidze correction to Swarm-measured data reduces the mismatch with IRI-modeled data but, as pointed out by Section 3.3. results and the comparison with ISRs data (Section 3.2), several differences between measured and modeled data are still present for what concerns both T e magnitude and the corresponding spatial, diurnal, and seasonal variations.
To better highlight the relation between measured and modeled data, the conditioned probability density function between them has been calculated. Specifically, we calculated the probability density of measured T e values, both original and Lomidze-corrected, conditioned by IRI-modeled ones in 50 K-wide range around a given value T 0 (ranging from 25 to 4975 K in steps of 50 K): In Equation (6), T e,measured is the T e measured by Swarm (i.e., T e,LP ) or corrected according to Lomidze et al. [24] (i.e., T e,Lom ), while T e,modeled is the corresponding one modeled by IRI and p is the conditioned probability density. Equation (6) was calculated for each T 0 equal to the T e,modeled values used to define the ranges for the conditioned probability density calculation. Conditioned probability density values are represented in Figure 9 as colored bins, while bins with p < 10 −6 are in white. For each [T 0 − 25; T 0 + 25) K range, the median is calculated and represented as a black circle, and corresponding error bars are the median absolute deviation values [46].
Median values of Figure 9 were used to perform a linear regression analysis and to obtain the corresponding slope, intercept, and Pearson correlation coefficient, namely: where cov is the covariance matrix and σ the variance. Linear regression best fit lines are represented in Figure 9 in magenta, with corresponding slope, intercept, and R values reported in each panel. The analysis was restricted to T e,modeled ≤ 2200 K because the linear relation between measured and modeled data fully breaks above this threshold, as clearly highlighted by the median values of Figure 9.
The results of Figure 9 confirm and complement those shown by Figures 2 and 3. Specifically, the slopes of the linear best fit, which range between 0.80 and 0.89 for the original Swarm data, after applying the Lomidze correction range between 0.97 and 1.00, while the intercept values decrease from 920-990 K to 140-320 K. This confirms the accuracy improvement provided by the Lomidze correction. Differently, the correlation coefficient (R) was slightly affected by the Lomidze correction because this correction is substantially linear, apart from the term dependent on N e,LP (see Equations (1)-(3)); anyway, the correlation is very high in both cases (R > 0.94). However, we have to take into account that these results are valid only in the range T e,modeled ≤ 2200 K, where a linear dependence with a slope near 1 between measured and modeled data is evident. A non-linear dependence is clearly visible for T e,modeled > 2200 K and in principle could be indifferently due to the Swarm LPs instrumentation, possible shortcomings of the Lomidze correction, and the IRI representation of the highest T e values. Recently, Pignalberi et al. [68] studied the occurrence of very high Swarm-measured T e values (T e,measured ≥ 6000 K), and corresponding spatial, diurnal, and seasonal varia-tions. As expected, very high T e values were recorded above all at high latitudes but with well-defined trends. Specifically, the occurrence of very high T e,measured values maximizes in winter and is minimum in summer, with intermediate values at both equinoxes. Moreover, the occurrence is asymmetric between hemispheres, being highest in the Southern hemisphere, where it can reach 40% of the total observations in winter. These extremely high T e values are not represented at all by IRI, which is limited to values below about 4300 K (see Figures 2 and 9). As a consequence, the departure from the linear dependence between T e,measured and T e,modeled for the highest values cannot be related to the application of the Lomidze correction. The fact that the occurrence of T e,measured ≥ 6000 K values is quite remarkable at high latitudes and that it usually happens at small spatial and temporal scales [68] drives us to the conclusion that the median climatological representation given by the IRI model cannot properly represent such very fast and steep T e variations. This is mainly due to the limited spatial (20 degrees in invdip latitude) and temporal (80 minutes in MLT) resolution of its current spherical harmonic expansion approach [26]. Nevertheless, we point out that the Lomidze correction is based on a dataset of T e measurements from Arecibo, Jicamarca, and Millstone Hill ISRs, which are all low-and mid-latitude stations. This means that no high-latitude measurement was considered and, most importantly, only values in the range T e,measured < 3500 K (in HG mode) were used to derive the calibration curves (see Figure 7 in [24] 16.0 • E)) in the Lomidze dataset to also characterize the high T e values usually recorded at auroral latitudes. Finally, the possible overestimation of T e values measured by Swarm LPs due to instrumental/environmental effects cannot be completely ignored. This is why, a more in-depth analysis, based on different datasets and measurement techniques, is needed to clarify the origin and reliability of the very high T e values recorded by Swarm satellites.

Conclusions
In this study, the main statistical median features of the topside ionosphere T e were investigated on a global basis through data collected by LPs on-board Swarm satellites from the beginning of 2014 to the end of 2020. Swarm T e measured data were compared with corresponding values modeled by IRI, and with data measured by Jicamarca, Arecibo, and Millstone Hill ISRs, to evaluate the reliability of both Swarm data and IRI model in the description of the spatial, diurnal, and seasonal trends of T e in the topside ionosphere.
We also investigated the effectiveness of the correction to Swarm T e values proposed by Lomidze et al. [24] through comparison with IRI-modeled data and ISRs-measured ones. Specifically, we found that the application of the Lomidze correction: • improves the agreement between Swarm data and corresponding ones modeled by IRI when the entire dataset is considered, for every Swarm satellite. This is attested by the average 400 K improvement in the mean residuals between Swarm and IRI after the Lomidze correction application; • does not alter either the dispersion of Swarm data around the mean and the correlation between Swarm and IRI, due to the linear character of the correction; • reduces the Swarm data RMSE from about 690-780 K to about 440-530 K, and RRMSE from about 27-29% to about 20-23%, when compared to IRI data; • generally improves the agreement between Swarm and ISR data. The improvement is particularly evident at Millstone Hill while it is lower at Arecibo and even lower at Jicamarca. Moreover, the correction is more effective during daytime than nighttime; • does not alter the linear relation trend between measured and modeled T e values in the range below 2200 K, but it improves the corresponding slope (closer to 1) and intercept (closer to 0) values.
The results obtained with ISRs data suggest that the calibration of Swarm T e data made by Lomidze et al. [24] is heavily biased towards mid latitudes and daytime hours due to the ISRs dataset they used, and to the fact that they cumulated the data of the three ISRs to derive a single calibration curve. As a consequence, the accuracy of Swarm T e data would benefit from the inclusion of the latitudinal and diurnal dependence in the calibration procedure. Swarm T e data are very valuable to study the relative changes of T e , but a more detailed correction is needed to obtain more reliable absolute values of T e especially at low and equatorial latitudes. Moreover, to better describe the high T e values, the inclusion of data recorded by high-latitude ISRs in the Lomidze dataset is advisable.
The orbital configuration of Swarm satellites and the extension of the dataset are very tailored for the study of the global features of T e with a spatial and time resolution hardly accessible in the past. Because of this, we were able to study the median spatial, diurnal, and seasonal trends of T e in the topside ionosphere in a very accurate fashion. Moreover, this allowed implementing an important statistical comparison between measured data and corresponding ones modeled by IRI. This comparison is beneficial both to validate the IRI model and to unveil the most important departures between measured and modeled data for a future inclusion of the Swarm T e dataset in the IRI model. The statistical comparison between Swarm and IRI evidenced that: • the largest differences emerge at magnetic equator latitudes at dawn, and at low and mid latitudes during daytime and nighttime, respectively; • IRI needs to be improved in the description of the morning peak at low latitudes. This can be achieved by increasing the order of spherical harmonics underlying the IRI Te description; • IRI needs to be improved in summer daytime for both hemispheres, and in the description of the high Te values characteristic of daytime auroral oval latitudes; • IRI data never go beyond about 4300 K, while Swarm data show values well beyond 5000 K, and in general Swarm values are higher than IRI ones.
The main statistical features of T e in the topside ionosphere have been pointed out and discussed in the paper and linked to the physical phenomena affecting the T e behavior at these altitudes, and to the existing literature on the topic. One of the most relevant results of this study is the characterization of the high T e values recorded by Swarm LPs, which are customary above all at high latitudes. These high T e values are not represented by the climatological description given by the IRI model. However, IRI has been developed to represent the median behavior of T e at climatological scales, and extreme (off-median) values or very fast and steep T e variations are out of the scope of the model. Nevertheless, a careful characterization of the high T e values recorded by Swarm LPs deserves particular attention for both modeling and scientific purposes.
In this paper, we did not deal with T e variations with the solar and magnetic activity. Since Swarm satellites are still in orbit and the mission will be supported for the next years, the continuous growth of the dataset will allow also studying the effects of the solar and magnetic activity variations on T e in the years to come. Moreover, we focused on the largescale spatial and temporal variations of T e in order to make a fair comparison with the IRI model that describes the climatological behavior of the ionosphere. However, the very-high resolution of Swarm LPs data also allows the study of small-scale T e features [68,69] which complement the results of this study and can be of help in the comprehension of the T e behavior in the topside ionosphere.