Inter-hemispheric asymmetries in high-latitude electrodynamic forcing and the thermosphere during the October 8–9, 2012, geomagnetic storm: An integrated data–Model investigation

Inter-hemispheric asymmetry (IHA) in Earth’s ionosphere–thermosphere (IT) system can be associated with high-latitude forcing that intensifies during storm time, e.g., ion convection, auroral electron precipitation, and energy deposition, but a comprehensive understanding of the pathways that generate IHA in the IT is lacking. Numerical simulations can help address this issue, but accurate specification of high-latitude forcing is needed. In this study, we utilize the Active Magnetosphere and Planetary Electrodynamics Response Experiment-revised fieldaligned currents (FACs) to specify the high-latitude electric potential in the Global Ionosphere and Thermosphere Model (GITM) during the October 8–9, 2012, storm. Our result illustrates the advantages of the FAC-driven technique in capturing high-latitude ion drift, ion convection equatorial boundary, and the storm-time neutral density response observed by satellite. First, it is found that the cross-polar-cap potential, hemispheric power, and ion convection distribution can be highly asymmetric between two hemispheres with a clear By dependence in the convection equatorial boundary. Comparison with simulation based on mirror precipitation suggests that the convection distribution is more sensitive to FAC, while its intensity also depends on the ionospheric conductance-related precipitation. Second, the IHA in the neutral density response closely follows the IHA in the total Joule heating dissipation with a time delay. Stronger Joule heating deposited associated with greater high-latitude electric potential in the southern hemisphere during the focus period generates more neutral density as well, which provides some evidences that the high-latitude forcing could become the dominant factor to IHAs in the thermosphere when near the equinox. Our study improves the understanding of storm-time IHA in high-latitude forcing and the IT system.

Inter-hemispheric asymmetry (IHA) in Earth's ionosphere-thermosphere (IT) system can be associated with high-latitude forcing that intensifies during storm time, e.g., ion convection, auroral electron precipitation, and energy deposition, but a comprehensive understanding of the pathways that generate IHA in the IT is lacking. Numerical simulations can help address this issue, but accurate specification of high-latitude forcing is needed. In this study, we utilize the Active Magnetosphere and Planetary Electrodynamics Response Experimentrevised fieldaligned currents (FACs) to specify the high-latitude electric potential in the Global Ionosphere and Thermosphere Model (GITM) during the October 8-9, 2012, storm. Our result illustrates the advantages of the FAC-driven technique in capturing high-latitude ion drift, ion convection equatorial boundary, and the storm-time neutral density response observed by satellite. First, it is found that the cross-polar-cap potential, hemispheric power, and ion convection distribution can be highly asymmetric between two hemispheres with a clear B y dependence in the convection equatorial boundary. Comparison with simulation based on mirror precipitation suggests that the convection distribution is more sensitive to FAC, while its intensity also depends on the ionospheric conductance-related precipitation. Second, the IHA in the neutral density response closely follows the IHA in the total Joule heating dissipation with a time delay. Stronger Joule heating deposited associated with greater high-latitude electric potential in the southern hemisphere during the focus period generates more neutral density as well, which provides some evidences that the high-latitude forcing could become the dominant factor to IHAs in the thermosphere when near the equinox. Our study improves the understanding of storm-time IHA in high-latitude forcing and the IT system. KEYWORDS inter-hemispheric asymmetry, ion convection equatorial boundary, ionosphere-thermosphere system, high-latitude electrodynamic forcing, IMF B y dependence, storm phase response, thermospheric neutral density

Introduction
Earth's ionosphere-thermosphere (IT) system can be highly affected by high-latitude forcing, which plays a significant role in the energy and momentum transfer between solar wind and the magnetosphere-ionosphere (MI) system (Richmond, 2011). This forcing strongly intensifies during geomagnetic storms and is often characterized by large-scale anti-sunward convection flows and equatorward expansions of the auroral oval (Cowley and Lockwood, 1992;Fuller-Rowell et al., 1994). As a fundamental element in storm-time magnetospheric energy transport, the field-aligned currents (FACs) are associated with the ion convection patterns (Shi et al., 2020), along with an increase in Joule heating both locally and globally (Deng et al., 2018). As a global phenomenon, storm-time high-latitude forcing and its effects on the IT system can extend from the polar region to lower latitudes, with characteristic responses at different local times (LTs) and hemispheres (Pi et al., 1997;Jakowski et al., 2005;Basu et al., 2008;Astafyeva et al., 2014). One of the principal expansions is that when the interplanetary magnetic field (IMF) B z suddenly turns negative, the enhanced magnetospheric convection cannot be fully shielded by the ring current in the magnetosphere and the region-2 FAC (Blanc et al., 1983;Goldstein et al., 2005). Different stormrelated IT responses have been reported, such as electron density variations including polar patches and total electron content (TEC), traveling ionospheric disturbances (TIDs), and penetration electric fields (Tsurutani et al., 2008;Tulasi Ram et al., 2009;Jin and Xiong, 2020).
Due to Earth's seasonal dipole tilt, the geomagnetic field configuration, and asymmetric high-latitude forcing (Hong et al., 2021), the storm-time IT responses can be significantly different between the Northern and Southern hemispheres, known as the IT system inter-hemispheric asymmetry (IHA). While all the causes are important, this study focuses on high-latitude forcing since the seasonal effect is considered to be reduced due to the data interval being near the equinox. It should be noted that the influence of the geomagnetic field is already included in high-latitude forcing since the forcing used in this study is the ultimate response of Earth's upper atmosphere. We cannot separate the asymmetry of the original magnetosphere sources from the geomagnetic field effects regardless of which season it is (Förster and Cnossen, 2013). Previous studies have shown that high-latitude forcing strongly manifests asymmetries with complex temporal and spatial changes during geomagnetic storms Reiff and Burch, 1985;Sandholt and Farrugia, 2007;Cousins and Shepherd, 2010), which could drastically affect the IT system in the two hemispheres. It has been emphasized by previous studies that the general FACs deduced from observations also exhibit hemispheric differences (Anderson et al., 2008;Green et al., 2009;Coxon et al., 2016;Workayehu et al., 2020). The IMF B y component in the geocentric solar magnetospheric (GSM) coordinate system, i.e., dawn-dusk direction, is thought to be one main cause of a number of asymmetric features in the magnetosphere and IT system (Walsh et al., 2014). Studies also showed that IMF B y causes obvious IHAs in auroral precipitation (Sandholt and Farrugia, 2007) and high-latitude Joule heating (McHarg et al., 2005). A full understanding of the IT system relies on knowledge of IHA. Despite these extensive studies, the global consequences and causes of the geomagnetic storm-related IHAs remain unknown.
General circulation models (GCMs) are widely used to study the IT system during geomagnetic storms. In order to specify the highlatitude forcing in GCMs, the most common approach is to utilize empirical models, e.g., using Weimer (2005) to specify convection patterns and using Fuller-Rowell and Evans (1987) or Newell et al. (2009) to specify auroral particle precipitation. However, since empirical models mainly describe the average conditions for a given state, they have issues in representing abrupt temporal changes or spatial distributions during a specific storm (Heelis and Maute, 2020). Typically, the two hemispheres are assumed to be mirror images considering a switch in IMF B y and dipole tilt, and the data from the two hemispheres have been combined without considering differences in data coverage. Certainly, more realistic specifications of external forcing are required to accurately simulate the storm-time global IT responses. For example, magnetohydrodynamic (MHD) models can also be used to specify the high-latitude forcing of GCMs, such as the connection between the Lyon-Fedder-Mobarry (LFM) MHD model and the Thermosphere-Ionosphere-Electrodynamics General Circulation Model (TIEGCM), and the coupled magnetosphere-ionosphere-thermosphere (CMIT) model Wiltberger et al., 2004). In addition, the Assimilative Mapping of Ionospheric Electrodynamics (AMIE) technique (Richmond and Kamide, 1988) can provide high-latitude electric potential and electron precipitation patterns used for driving GCMs (Lu et al., 2014;Lu et al., 2020). The Active Magnetosphere and Planetary Electrodynamics Response Experiment (AMPERE) gives another alternative way to calculate the high-latitude electric potential by using the derived FAC and pre-defined electron precipitation or ionospheric conductance patterns (Maute et al., 202l;Robinson et al., 2021;Zhu et al., 2022), which shows improved agreement between simulations and experimental data in driving GCMs for storms (Maute et al., 202l;Maute et al., 2022;Zhu et al., 2022b).
Motivated by the FAC-driven approach and AMIE electron precipitation patterns, the causes and consequences of IHAs of highlatitude forcing and IT responses during the October 8-9, 2012, storm have been investigated systematically with both data and models from over the high, middle, and low latitudes. The remaining part of this paper is organized as follows: Section 2 introduces the methodology of this study. Section 3 overviews the geophysical conditions for the October 8-9, 2012, geomagnetic storm in brief. Section 4 presents the main results of this study through comprehensive data-model and model-model comparisons. Section 5 summarizes the paper.

Methodology
The present study utilizes integrated data from multiinstrument observations and models. The satellite observations are given in Section 2.1. Section 2.2 presents the GITM model. The newly developed FAC-driven procedure is given in Section 2.3. Section 2.4 illustrates how we determine the ion convection equatorial boundary according to DMSP observation.

DMSP ion drift
Data from three Defense Meteorological Satellite Program (DMSP) satellites (i.e., F16, F17, and F18) are used in this study. They flew in sun-synchronous orbits at an altitude of~840 km with an inclination angle of~98.8° (Rich and Hairston, 1994) primarily in the dawn-dusk direction during the October 8-9, 2012, geomagnetic storm. The cross-track ion drift (V y ) is measured by using the onboard Special Sensor for Ions, Electrons, and Scintillation (SSIES), which has a temporal resolution of 1 s. In this study, the data with the quality flag of 1 (i.e., most reliable) are used. A linear baseline correction is applied to the original V y data to remove some co-rotation effects and to ensure V y is zero at 45°|magnetic latitude| (|MLAT|). Afterward, a 13-point sliding window is applied to the corrected data in order to reduce very-highfrequency fluctuations and extract the large-scale ion convection at high-latitude regions.

GOCE neutral density data
The Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) satellite was launched on 17 March 2009. The satellite was a polar-orbiting satellite (inclination angle: 96.5°) and flew in a near-circular orbit with an altitude from 250 km to 280 km . The solar local times (SLTs) of the ascending and descending nodes of the GOCE were about 19 and 07 h, respectively. The neutral mass density (hereafter, neutral density for simplicity) was measured by using an onboard accelerometer with a temporal resolution of 10 s . To reduce the variation caused by the satellite orbit altitude change, the neutral density data are normalized to a constant altitude of 270 km using the NRLMSIS 2.0 model (Emmert et al., 2021). More details of the density normalization technique can be found in Bruinsma et al. (2006).

AMPERE FAC data
The AMPERE high-latitude FAC densities are derived from the horizontal magnetic field perturbations measured by 66 Iridium satellites which are distributed along six longitudinally equally spaced orbital planes (Anderson et al., 2002). Each Iridium satellite flies in a near-polar orbit at an altitude of 780 km with an orbital period of 104 min. The radial current components are calculated by fitting the magnetic perturbations measured within a 10-min time window using spherical cap harmonic basis functions (Waters et al., 2001;Green et al., 2006;Waters et al., 2020). Specifically, the patterns are fitted with a longitude order of 5 and latitude order of 20 (i.e., 3°) between colatitude 0°and 60°in the Altitude-Adjusted Corrected Geomagnetic (AACGM) coordinates (Baker and Wing, 1989). The temporal resolution of the FAC data is up to 2 min, and the spatial resolution of the FAC data is 1°in MLAT and 1 h in the magnetic local time (MLT) . In this study, AMPERE FAC densities with a magnitude below 0.2 μ A/m 2 mainly distributed in the polar cap region and below 50°|MLAT| are removed, defined by being smaller than the noise level. In addition, due to differences in the geomagnetic field magnitude and the curvature of the magnetic field, the FAC pattern is mapped from 780 km to the APEX reference height of 110 km based on a factor of 1.343. Using a bilinear interpolation, the 10-min resolution AMPERE FAC patterns are spatially interpolated to the electrodynamic solver's grids and then temporally interpolated into a 2 s cadence at which the simulation is performed using linear interpolation.

GITM
The Global Ionosphere and Thermosphere Model (GITM) is a three-dimensional first-principle general circulation model for Earth's thermosphere and ionosphere system . The GITM solves continuity, momentum, and energy equations in a spherical coordinate framework to calculate the density, velocity, and temperature of neutrals, ions, and electrons. The GITM has a flexible grid size and can use a stretchable grid in latitude and altitude. Moreover, the GITM relaxes the hydrostatic assumption in the vertical direction, which allows for the evaluation of nonhydrostatic impacts on the IT system (Lin et al., 2017;Deng et al., 2021). The global ionospheric electrodynamic solver in the GITM is the NCAR-3D electrodynamic model (hereafter NCAR-3D, Maute and Richmond, 2017), which is coupled into the GITM by Zhu et al. (2019). More details about the GITM can be found in Ridley et al. (2006).
In this study, two main GITM simulations were run with different high-latitude electrodynamic forcing: driven by empirical models and by more realistic patterns. A summary of the high-latitude forcing settings of these runs is given in Table 1.
For run 1, the high-latitude electric potential is specified by the Weimer 2005 model (hereafter W05, Weimer, 2005), and electron precipitation is specified by the Auroral energy Spectrum and High-Latitude Electric field variabilitY (ASHLEY) model . Both W05 and ASHLEY are driven by IMF and solar wind data. For run 2, the high-latitude electric potential is calculated using AMPERE FAC data and the electron precipitation is specified by the AMIE patterns. The data inputs to AMIE patterns for this storm event include the horizontal magnetic perturbations from 261 ground magnetometers (with 205 in the NH and 56 in the SH) and the electron precipitation measured by the Special Sensor of Ultraviolet Spectrographic Imager (SSUSI) onboard DMSP F16-18 satellites. Details about the AMIE can be found in Lu (2017). More realistic patterns which can resolve additional dynamic and spatial structures than the statistically averaged empirical models are used in run 2 as compared to run 1.
In each simulation, the GITM was run with a spatial resolution of 5°in geographic longitude, 2.5°in geographic latitude, and 1/ 3 scale height in altitude. The time step of both simulations is 2 s. Realistic IMF B y and B z , solar wind, and F 10.7 data from the CDAWeb OMNI data product are used as model inputs. Specifically, in this study, the GITM simulations run 1 and run 2 are carried out with a 3-day pre-run (00UT (universal time) 10/01-00UT 10/03), a 5-day quiet time (00UT 10/03-00UT 10/08), and a 2day event time (00UT 10/08 -00UT 10/10) based on their corresponding high-latitude forcing. The purpose of a pre-run is to ramp up the model from an initial condition to a diurnally reproducible state during the time of interest , and the outputs during that period typically are not used for scientific study.
Frontiers in Astronomy and Space Sciences frontiersin.org

FAC-driven procedure
This subsection gives a brief overview of the FAC-driven procedure used in this study, and details can be found in Maute et al. (2021) and Zhu et al. (2022). The first step is to calculate the high-latitude electric potential (Φ R ) based on the current continuity equation: where J N/S mr is the FAC input in the Northern (N) and Southern (S) hemispheres.
λ m , ϕ m are the magnetic latitude and longitude, respectively, in the modified apex coordinate.
R ( R E + h R ) is the radius to the conducting ionospheric layer with h R = 110 km.
σ R , Φ R are the reference conductivity and high-latitude potential, respectively. p c is the ratio factor changing with magnetic latitude. Σ is the integrated ionospheric conductivities term along a field line.
K D is the integrated neutral dynamo current term integrated along a field line.
The superscripts N and S represent quantities in the Northern and Southern hemispheres, respectively. The specific definitions of the ionospheric conductivities Σ ij (i, j ϕ, λ) and neutral dynamo current K D mj (j ϕ, λ) can be found in Richmond (1995). To specify the boundary of the high-latitude electric potential (Φ R ), a ratio factor p c varying with magnetic latitude λ m is used for both hemispheres: p c is 1 poleward of | λ m | = 45°and is 0 equatorward of | λ m | = 40°and linearly changes between 0 and 1 for the region 40°< | λ m | < 45°. This setting ensures that the high-latitude electric potential (Φ R ) is zero equatorward of | λ m | = 40°.
Once the high-latitude electric potential (Φ R ) in each hemisphere is calculated, the second step is to calculate the global electric potential (Φ) based on Φ R : where Σ T is the total integrated conductivities of both hemispheres. K DT is the total neutral dynamo of both hemispheres. p is the ratio factor changing with magnetic latitude.
Here, the superscript T means the sum of quantities from both hemispheres. p is 1 equatorward of | λ m | = 50°and is 0 poleward of | λ m | = 55°and linearly changes between 1 and 0 for the region 50°< | λ m | < 55°, which means that only the neutral dynamo potential is used equatorward of the lower boundary, e.g., | λ m | = 50°and only the high-latitude potential is used poleward of, e.g., | λ m | = 55°and a combination in between. In this study, similar to the settings in TIEGCM, the NH high-latitude electric potential Φ R is used to determine Φ in midand low latitudes, which is hemispherically symmetric since the field lines are assumed to be equipotential. After that, the high-latitude electric potential in the SH is replaced by (Φ R ) in the SH.

Identification of the ion convection equatorial boundary
In this study, the latitudinal expansion of the high-latitude ion convection to mid and low latitudes was investigated based on the cross-track ion drift V y (nearly zonal direction below (|geographic latitude (GLAT)| = 70). The equatorial boundary is identified as the latitude at which the zonal ion drift suddenly gives a threshold. Heelis and Mohapatra (2009) suggested that the DMSP-observed zonal ion drift flow can be used to obtain the convection pattern equatorward expansion and contraction boundary, which is approximately consistent with the equatorward edge of the region-2 FAC. To quantitatively determine the equatorial boundary, a threshold of 15 m/s per degree in the latitudinal gradient of the zonal ion flow is adopted in this study (Heelis and Mohapatra, 2009). In addition, a criterion of 300 m/s greater than the quiet-time background flow is also applied in order to distinguish the enhancement from the quiettime ion flow background (Hairston et al., 2016). Frontiers in Astronomy and Space Sciences frontiersin.org Figure 1 shows an example of the NH dusk-side ion convection equatorial boundary determined from the DMSP F16-measured zonal ion flow during the October 8, 2012, storm day. Specifically, the radius represents the GLAT, and the polar angle represents the geographic longitude (GLON) at the DMSP F16 fixed dusk-side LT,~18 SLT. The zonal ion flows were binned into 1°in GLAT between 10°and 80°GLAT and 2/3 h intervals (i.e., 10°in GLON). Figure 1A shows the quiet-time zonal ion flow background based on 7-day DMSP dusk-side measurements during October 1-7, 2012; Figure 1B shows the storm period zonal flow (around  Frontiers in Astronomy and Space Sciences frontiersin.org 14 trajectories per day) on October 8. Subtracting the zonal ion flow background 1A from the storm time 1B, the difference is shown in Figure 1C. By combining these two criteria described in the aforementioned paragraph, the ion convection boundary can be determined with black dots consistently for all longitudes, i.e., UT (UT = GLON/15-18 LT). Furthermore, we have also checked that using a larger or smaller threshold such as 400 m/s or 200 m/s only alters the boundary in the GLAT within 2°for this specific storm event.
3 Geophysical conditions Figure 2 shows the IMF components B y and B z , solar wind velocity V x , interplanetary electric field (IEF, E y = -V x x B z ), SYM/H index, and auroral electrojet (AE) index for the moderate geomagnetic storm during October 8-9, 2012. As shown in Figure 2D, the storm sudden commencement (SSC) occurred at 05:16 UT on October 9 2012 and the SYM/H index reached~−100 nT around 10:20 UT on October 8 during the first storm main phase (05:30 UT to 10:30 UT on October 8). The second main phase occurred between 18:00 UT on October 8 and 08:30 UT on October 9 during which the SYM/H index decreased to a minimum of~−116 nT at 02:07 UT on October 9. During the second main phase, the IMF B z component ( Figure 2A) and IEF ( Figure 2B) remained at around −15 nT and 5 mV/m, respectively. Meanwhile, the AE index showed rapid disturbances with basically above 600 nT and can reached up to 1400 nT ( Figure 2C). The IMF B y component stayed at around +8 nT before 01:20 UT and then quickly reversed to~−8 nT. This study mainly focuses on the second main phase which offers us a great opportunity to examine the IMF B y effects on the IHAs of high-latitude electrodynamic forcing and its impacts on the thermosphere.

Results and discussion
In Section 4.1, results from the empirical model-based run 1 (W05 and ASHLEY) and realistic pattern-driven run 2 (FACdriven and AMIE electron precipitation) are compared with DMSP F16-18-observed cross-track ion drift V y and the ion convection equatorial boundary. Section 4.2 illustrates the IHAs in highlatitude forcing, i.e., ion convection and auroral electron precipitation. The asymmetric storm phase response and B y dependence of the ion convection equatorial boundary in both hemispheres are also investigated. The storm-time global evolutions of the thermospheric neutral density simulated by the GITM are compared with the GOCE satellite-measured result. The IHA and physical mechanism are examined in Section 4.3.

Data-model comparisons of the ion convection 4.1.1 Data-model comparison of the ion drift at high latitudes
To determine how well the model output represents the ionospheric ion convection, a data-model comparison of the cross-track ion drift V y has been conducted along DMSP trajectories during the focus period. As shown in Figure 3, six examples from DMSP F16-18 and the corresponding GITM runs are compared with the top and bottom panels representing the Northern and Southern hemispheres, respectively. For each plot, the black dots show the observed cross-track ion drift, and the blue and red curves represent the empirical model-based run 1 (W05) and realistic pattern-driven run 2 (FAC-driven and AMIE electron precipitation), respectively. In general, the cross-track V y from both GITM runs is consistent with the DMSP results. When comparing the two sets of runs, the FAC-driven V y from run 2 can better capture the realistic ion drifts in high latitudes than the W05-based results from run 1, especially inside the auroral zones.
To investigate how GITM simulations capture the highlatitude ion drifts V y from a statistical perspective, Figure 4 compares the simulated and measured cross-track ion drift at the region poleward of 45°|MLAT| from all DMSP trajectories during the focused period for runs 1 and 2. Linear fitting is performed on the data, and the slope (k) and y-intercept (b) are calculated. In addition, the Pearson correlation coefficient (r) and the root mean square error (RMSE) are also calculated. Basically, the cross-track ion drifts V y from both GITM simulations have smaller magnitudes than the observed crosstrack V y since the slope of the best-fit line of the fitted line is much smaller than 1. This underestimation in the simulation results could be due to the underestimation of high-latitude forcing such as the FAC magnitude. This is a known property of AMPERE FAC data in the dusk-side R2 current. However, the slope of the best-fit line is slightly higher in run 2 than in run 1. Meanwhile, compared with run 1, run 2 shows a greater correlation (e.g., for DMSP F16 run 2 = 0.81 vs. run 1 = 0.68) and a smaller RMSE between the simulated and observed ion drifts (336.20 compared to 426.33) than run 1. Overall, Figure 4 indicates that the FAC-driven GITM simulation can better reproduce the high-latitude ion drift measured by DMSP satellites than the W05.

Data-model comparison of the ion convection equatorial boundary
The data-model comparisons of high-latitude ion convection given in Figure 3 and Figure 4 show a general consistency between the measured and simulated ion drift on the dawn and dusk sides, especially for run 2. Moreover, we would like to see how GITM simulations perform in response to the ion convection equatorial boundary. Figure 5 shows a scatterplot of the ion convection equatorial boundary in the GLAT for both hemispheres measured by DMSP F16 versus those from runs 1 and 2 during the focused storm period. In general, the simulations reproduce the DMSP-measured equatorward boundary latitudes well, with a better representation in run 2. Specifically, the correlation coefficients to DMSP F16 are 0.81 in the NH ( Figure 5A) and 0.91 in the SH ( Figure 5C) for run 1 and 0.92 in the NH ( Figure 5B) and 0.95 in the SH ( Figure 5D) for run 2. Moreover, as shown in Figures 5A, C, W05 has a clear underestimation of the convection expansion below 60°for the NH and at all the latitudes for the SH. This underestimation of the ion convection equatorial boundary may have some effects on other electrodynamical processes such as the ion-neutral interactions in that region. Overall, both runs show good agreement with the DMSP-observed convection equatorial Frontiers in Astronomy and Space Sciences frontiersin.org  Frontiers in Astronomy and Space Sciences frontiersin.org 07 boundary. The comparison between runs 1 and 2 is to give the community a rough idea about how well the FAC-driven technique could improve high-latitude electric field specification which is also a justification for using the FAC-driven technique in this study. For later scientific studies in Section 4.2 and Section 4.3, we would like to focus on the GITM run 2 since it has advantages over run 1 as discussed. A summary of the statistical parameters from the data-model comparisons of the two runs can be found in Table 2.
4.2 IHA of high-latitude forcing and the ion convection equatorward expansion 4.2.1 IHA in the high-latitude electrodynamic forcing Figure 6 shows the temporal variations of the cross-polar-cap potential (CPCP, top panel) and hemispheric power (HP, bottom panel) during October 8-9, 2012, from run 1 (left) and run 2 (right). The vertical dashed line marks the time when IMF B y is reversed (00:  20 UT on October 9) during the storm main phase we focused on in this study. All the red and blue lines correspond to the NH and SH, respectively. Following Hong et al. (2021), the asymmetry index (AI) is used to quantify the temporal changes of IHA for a given quantity: where Y NH and Y SH stand for the quantity in the NH and SH, respectively. As shown in Figure 6A, CPCPs from run 1 are almost the same in the NH and SH during the focused period and the CPCP stays around 145 kV. As for run 2, significant asymmetries between the NH and SH can be seen in Figure 6B even for the recovery phase (e.g., after 12 UT on October 9). An evident asymmetry index (AI up to −50%, see Figures 6C, D) occurs in the CPCP during the focused period, especially before IMF B y is reversed. In addition, the CPCP also exhibits more pronounced temporal variations than run 1. Similarly, as shown in the bottom two panels of Figures 6E, F, the hemispheric power in run 2 exhibits more dynamic variations and greater IHA than in run 1. The AIs can reach 75% and −22% in run 2 ( Figure 6H) and run 1 ( Figure 6G), respectively. Different from the CPCP, the HP tends to have greater AI in run 2 when the IMF B y is negative (after 00: 20 UT on October 9). Additionally, the HP in the NH is generally larger than that in the SH in run 2 while the opposite is true in run 1. This is probably due to the fact that the empirical model is mainly dependent on the given IMF, solar wind, and input seasonal conditions. Another reason is that, as described in Section 2.2, the AMIE patterns have better data coverage from the NH than that in the SH.
To our knowledge, the IHAs in high-latitude forcing could be due to the combined action of the season-related dipole tilt and the nonzero B y component, which contribute to the asymmetric interaction in the dayside reconnection (Reistad et al., 2021).    , 2009). This is particularly true when there is an east-west component of the IMF (IMF B y ) . The lobe reconnection due to IMF B y causes the magnetic flux to build up asymmetrically (Lu et al., 1994). For example, the responding speed to the solar wind conditions in the two hemispheres can be different, which can result in different temporal variations of the magnetic flux between the two hemispheres (Milan et al., 2020). Therefore, this phase mismatch may result in some instantaneous asymmetry of the CPCP while it should be reduced significantly when the response in both hemispheres is fully developed. In addition to the reconnection for the IMF B y , the field line potential drops across the hemispheres, and the asymmetries in the current systems produced by differing ionospheric conductivities and neutral dynamos in two polar regions contribute to the possible asymmetry in the CPCP as well. Data-model comparisons are still absolutely crucial to validate the simulation results. While the GITM simulations are able to capture the fundamental features of the ion drifts at high latitudes and the convection boundaries at mid-latitudes, the potential strong asymmetry of the CPCP presented in our simulation has yet to be corroborated with realistic observations. We will make careful comparisons with SuperDARN and other measurements in the follow-on study.

IHA in the ion convection equatorial boundary and its IMF B y dependence
An enhanced geomagnetic storm can lead to large perturbations in the high-latitude ion convection and then extend to lower latitudes. From Section 4.1.2, we have demonstrated that the ion convection equatorial boundary from DMSP observations and GITM simulations shows a fairly good agreement on both the dawn and dusk sides. To have a detailed picture of the UT variations during this storm, Figures 7A, B directly show the temporal variations of DMSP F16-measured ion convection boundary in GLAT on dawn and dusk sides during the focused IMF B y reversed period (18 UT on October 8-10 UT on October 9). Over this period, the NH dawn-side convection boundary (red dashed line) only undergoes a small variation (<5°) around 63°N GLAT, while the dusk-side convection boundary can expand to as far as 45°N (red solid line). As for the SH, both the dawn and dusk sides expand to lower GLATs during this interval of time, with the observed lowest boundary at 45°S and 40°S, respectively. One of the most interesting features is that these two sides reach their lowest | GLAT| under opposite IMF B y conditions: the dawn side first expands under positive B y and then contracts during negative B y (blue dashed line), while the dusk side shows the opposite (blue solid line). This result could be explained by the anti-correlated responses between the dawn and dusk sides to B y (e.g., Cowley, 1981;Reiff and Burch, 1985;Kabin et al., 2003;Tenfjord et al., 2015), the SH dawnside convection cell is much stronger than the dusk-side cell under positive B y , while the dusk side is more favorable under negative B y in the SH.
In terms of the ion convection boundary comparisons between the two hemispheres, the temporal expansion in the NH dusk side (7A) follows a similar tendency to changes in the SH dawn side (7B),  (Cousins and Shepherd, 2010), a stronger convection cell on average expands to lower latitudes. On the other hand, in contrast to the expansion in the SH dusk side (7b), the NH dawn-side convection expansion (7a) seems to be less important, implying a relatively stable convection cell on the dawn side. Moreover, this response is much stronger in the SH than in the NH, suggesting a stronger B y dependence of the two cells in the SH . A possible explanation to the dawn-dusk asymmetry and IHA in the expansion is that the ion convection is subject to asymmetries in FACs and ionospheric conductivities. However, these line plots represent a limited view of ion convection pattern variation at two fixed LTs. The aforementioned hemispheric asymmetries between the NH and SH on both dusk and dawn sides will be further discussed in the following section.  Figures 7C, D, a stronger dusk-cell with minima value −51 kV (d) corresponds to a shrinking convection boundary (a and d, red circle). Furthermore, the dawn-cell increased from 29 kV (c) to 40 kV (d) with subtle changes in the boundary location. Again, the NH displays a clear stretch of the dawn-side cell in the noon-midnight direction instead of an equatorward expansion, which could explain the static phenomena along the DMSP F16 trajectory on the NH dawn side as shown in Figure 7A (red dashed line). This result suggests that it might be insufficient to determine the storm period ion convection expansion based on only the dawn and dusk LTs. Overall, the ion convection pattern simulated by the GITM for these two specific periods is in general agreement with the observed convection boundary expansions. Moreover, a stronger convection cell does not necessarily associate with a lower equatorial boundary.

Causes of the IHAs in the high-latitude convection and its equatorial boundary
In addition, we also conducted data-model comparisons of the equatorial boundary expansion at specific times. As shown in both the linear and polar plots, the red and blue triangles refer to the dawn side for the NH and SH, respectively. Similarly, the red and blue circles refer to the NH and SH but for the dusk side. The snapshot times T1 and T2 represent the corresponding positive B y case at~22 UT on October 8 and negative B y case~04 UT on October 9, respectively, which are marked in the left line plots. The subscripts N and S refer to the NH and SH, respectively. The boundary expansion features in the 2D contour from FAC-driven GITM run 2 are in a good agreement with the line plots from the DMSP F16 observations. For example, the blue triangle T 2S ( Figure 7B, dawn side) almost overlaps the expansion of the dawn-side potential cell ( Figure 7F, blue color) at around 67 o S GLAT.
To identify the possible mechanisms for the observed IHAs in the ion convection patterns, as well as their equatorial boundaries, a separation of the FAC and auroral electron precipitation (associated with the ionospheric conductance) is needed. In reality, the highlatitude FAC and auroral precipitation are connected and respond simultaneously to the changes in geophysical conditions, so it is impossible to completely separate the effects caused by these two variables. However, numerical simulations could help solve this issue by replacing the real forcing with symmetric patterns. We first examine the GITM run 2 (simulation setups in Table 1), which the high-latitude forcing is based on the FAC-driven potential and AMIE electron precipitation patterns. In GITM run 3, the AMIE electron precipitation patterns in the SH are replaced by the mirror precipitation patterns in the magnetic coordinate from the NH accordingly. Therefore, the differences in the SH between run 2 and run 3 illustrate the contributions from the auroral precipitation. Figure 8 shows the Southern hemispheric FAC-driven convection patterns and corresponding auroral energy flux from GITM run 2, run 3, and their differences at 04 UT on October 9 as an example. Specifically, due to the asymmetric electron precipitation, the auroral electron energy flux between the NH and SH is significantly different in both distributions and magnitudes: the SH has a continuous and complete aurora arc structure ( Figure 8D) while the NH is more like a discrete aurora ( Figure 8E). As shown in Figure 8F, the maxima difference in energy flux can vary from −45% to 55%. This result could be explained by the larger hemispheric power captured at that specific time (see Figure 6F); even the hemispheric power in the NH can be overwhelmingly superior to that in the SH during the focused period. However, the SH and NH electron energy flux-based ion convection ( Figures 8A, B) illustrates quite similar distributions, as well as the maxima and minima values for the dawn (33 vs. 34 kV) and dusk (−69 vs. −67 kV) cells. Compared to Figure 8C, the differences caused by the asymmetric auroral precipitation can be ignored at this specific time. Moreover, by comparing the convection contour lines, one can find that the distributions, especially the equatorial boundaries of the convection, are more likely determined by the FAC patterns (R2 FAC edge) while its strength may also depend on the magnitude of the ionospheric conductance-associated precipitation. With this discussion, the effects from the FAC and auroral precipitation can be roughly isolated. It should be noted that this conclusion depends critically on the assumption of the FACdriven procedure. Similar plots for all the UT times during the focus period indicate the same conclusion as the example time shown in Figure 8.

IHAs in the storm-time thermosphere: Neutral density response dρ
To investigate the IHAs in the thermosphere during this storm, the neutral density observed by the GOCE satellite was first divided into ascending (dusk side) and descending (dawn side) trajectories. After that, the daily ascending trajectory data are binned into UT Frontiers in Astronomy and Space Sciences frontiersin.org and GLAT grids with a bin size of 1 h and 1°in GLAT. We then applied a linear temporal interpolation in UT to project the data into uniformly distributed time series for both quiet-time background days and the storm-time days. Similarly, the approach is also used for GITM run 2-simulated neutral density extracted along the GOCE satellite trajectories for both quiet and storm days. Figure 9A shows the observed neutral density response, dρ, on the dusk side as a function of UT and GLAT. The quiet-time neutral density background (given by the 2-day average of October 4 and 5) is subtracted from the storm-time observations. During the period of elevated geomagnetic activity, as indicated by the Ap and SYM/H indices shown in Figure 9D, it is clear that the observation generally shows density enhancements which last for several hours, and the neutral density increases after October 9. During the focus period (i.e., green shaded area), the NH shows dρ enhancements between 60 o −80°N GLAT between 19:30 UT on October 8 and 00 UT on October 9 while the SH polar region represents stronger enhanced dρ at around 21:45 UT on October 8. The roughly 2-h difference in the responses between the NH and SH might be caused by asymmetric energy deposition such as the Joule heating in the two auroral regions. Another possibility is that the GOCE satellite trajectories are in different LTs in the NH and SH above 70°| GLAT|. Shortly after this period, both the neutral density in the NH and SH undergoes significant large-scale enhancements and the enhancements propagate equatorward and finally arrive at 30°S at 8 UT on October 9. Basically, the GOCE satelliteobserved dusk-side dρ has some properties associated with IHA: the density-enhanced time, the density magnitudes, and the latitudinal propagation. Figure 9B is similar to Figure 9A but for the dρ from GITM run 2 (FAC-driven potential and AMIE electron precipitation). Specifically, the run 2-simulated neutral density over the same quite-time during October 4 and 5 is used as the quiet-time background subtracting for GITM as done for GOCE data. In general, the simulated dρ is consistent with the GOCE-measured results though some differences exist. For example, the much greater dρ on October 9 than on October 8 can be captured. In agreement with Figure 9A, Figure 9B shows an equatorward propagation of the neutral density enhancement. It is also worth noting that there are some differences in the data-model comparison: neutral density from GITM run 2 is smeared and has less mesoscale structures than those in the GOCE observations ( Figure 9A), which may be due to three major reasons: 1) the latitudinal resolution in mid-and low latitudes of the observation (~0.6 deg) is much higher than the simulation (2.5 deg); 2) the mesoscale forces may not be included in the simulations Sheng et al., 2021); and 3) large-scale forcing is not perfect.
To quantify the asymmetric density features between the two hemispheres, Figure 9C shows the temporal changes of the AI for the averaged dρ at > 40°|GLAT| on the dusk side for both GOCE (black dot line) and GITM (red dot line). Apparent density asymmetries between the NH and SH can be identified, with the AI varying between −36% and 20% for the observations and between −18% and 20% for GITM results. As expected, the temporal changes of AI show some similarities and differences Frontiers in Astronomy and Space Sciences frontiersin.org between GOCE and GITM results. Both the observational and simulation results show that the AI in dρ reaches its first maxima and minima at around 11 UT and 13 UT on October 8, respectively. During the storm period slightly after B y is reversed, the asymmetry index reached its minimum value at-36% and −18% for the observational and simulation results, respectively. Moreover, a significant positive peak can be found during the recovery phase for both results but the magnitude in the GITM simulation seems to be much greater. As for the shaded focus period, the AI tends to be positive values (i.e., the NH has more density than the SH) under positive B y rather than being significantly negative after B y reversed into negative polarization, especially for the simulated results shown in the red line, implying that dρ might have some IMF B y dependence associated with high-latitude forcing and energy dissipation.
To help understand the hemispheric difference in neutral density, the asymmetry index of the total (hemisphericintegrated) Joule heating >40°|GLAT| between the two hemispheres from GITM run 2 is also plotted in Figure 9C as a reference for the energy dissipation at high latitudes (blue dashed curves). The most remarkable feature is that negative AI (i.e., the SH has more total Joule heating than the NH) occurs shortly after the focus period and lasts almost to the end of the day, October 9, illustrating that the SH dusk side has more energy input during this period. Moreover, as shown in Figure 9C, between 12-20 UT on October 8, more Joule heating was deposited into the NH (positive AI), resulting in increased density in the NH in density dρ during 16-22 UT on October 8. Immediately after this, greater Joule heating was dissipated into the SH with the minima of AI (~-60%) near 05 UT on October 9 in Joule heating. In general, the total Joule heating can lead to changes in neutral density enhancement, with 2-3 h of delay for the thermosphere to respond to the energy input change during storm times (Wang et al., 2020). On the basis of the aforementioned features, it can be summarized that the variations and IHA in the dρ could be significantly due to the contributions from the Joule heating dissipations.
The storm-time IT system responses are largely controlled by two factors: external high-latitude forcing from the solar wind-magnetosphere-ionosphere interaction and the local conditions in the thermosphere and its embedded ionosphere. Therefore, the IHAs in the thermosphere could also come from Frontiers in Astronomy and Space Sciences frontiersin.org several other mechanisms as discussed in the introduction. As described by several previous papers (Laundal and Østgarrd, 2009;Ohtani et al., 2009;Reistad et al., 2015;Laundal et al., 2017), the intrinsic North-South differences in the geomagnetic field could introduce asymmetric interaction of MI coupling in the two hemispheres. It should be pointed out that the IHA due to asymmetric tilt angles and displacements between the two hemispheres is daily UT differences. For example, at any given UT, the two hemispheric auroral zones will be one near nighttime and the other located in the dayside, which implies different solar fluxes would also be concerned during this period. Certainly, the absolute neutral density variation is also subject to the background thermospheric density before the storm, which provides important preconditions (McGranaghan et al., 2014). Previous studies have shown that both the season and geomagnetic field can contribute to the IHAs in the thermospheric neutral density (Emmert, 2015). Further examination shows that no obvious IHAs or enhancements could be found in the neutral density and Joule heating during the quiet-time background, i.e., October 4 and 5 from both GOCE and GITM results; given its proximity to the equinox, the result reveals that the storm-time intensified high-latitude forcing and associated Joule heating could be the dominant factor causing the IHAs to the neutral density dρ in the thermosphere. So far, we have discussed the storm-time responses in highlatitude forcing, Joule heating dissipation, and the thermosphere, all of which have significant IHAs as illustrated by measured and simulated results. To complete the pathways of IHAs from highlatitude forcing to Joule heating, similar to the inter-hemispheric comparisons in Figure 6, we conducted the same examination for the total FAC and total Joule heating (from run 2). As shown in the left part of Figures 10A, C, the total FAC indicates small positive asymmetry during the focused shaded period between the two hemispheres as expected for the equinox condition, with AI no more than 25%. In terms of the total Joule heating shown on the right, an overall negative AI (up to −75%, Figure 10D) can be found during the same period, indicating that more heating is dissipated into the SH. As a conversion between the electromagnetic energy and mechanical form, Joule heating is conventionally taken in the form J · E in the neutral frame of reference where E is the electric field and J is the current density (Vasyliūnas and Song, 2005). Therefore, the negative IHAs in the total Joule heating could be due to the negative IHAs shown in the CPCP ( Figure 6B and Figure 10B). Bringing all the results together allows us to generate a global and systematic perspective of the IHA performance during this storm event.

Summary
The inter-hemispheric asymmetries in the global IT system have been investigated during the October 8-9, 2012, geomagnetic storm for a focused IMF B y reversal period. The IHA in the high-latitude convection and auroral precipitation, the convection equatorial boundary, and the thermospheric neutral density have been investigated by combining data analysis and GITM simulations. The major conclusions are summarized as follows: 1) This study demonstrates the advantage of the FAC-driven technique in the global IT system study. Simulations based on realistic, spatial FAC patterns can better capture the high-latitude ion drift and the convection equatorial boundary than the empirical model-driven results. A summary of statistical parameters for the data-model comparisons of the two simulations is given in Table 2. 2) Significant IHAs have been identified in the FAC-driven simulation in both the high-latitude forcing and the IT system, with AI up to −50% in CPCP, 75% in hemispheric power, and −36% in neutral density. The IHAs show a B y dependence, especially for the convection boundary. 3) By comparing FAC-driven convection patterns based on precipitation from the SH itself or the mirror pattern from

FIGURE 10
Inter-hemispheric comparisons of the total field-aligned current (FAC, (A), total Joule heating (B), and their corresponding asymmetry index (AI, (C, D) during the October 8-9, 2012, geomagnetic storm. All the results are from run 2 (driven by data-based realistic patterns) since run 2 has obvious advantages in reproducing the high-latitude forcing. The vertical purple dashed line marks the IMF B y reversal time during the storm main phase.
Frontiers in Astronomy and Space Sciences frontiersin.org the NH, our results suggest that the distribution of the ion convection, especially the equatorial boundary, is more likely to be determined by the R2 FAC edge while its intensity may also depend on the ionospheric conductance associated with aurora precipitation energy flux. 4) The FAC-driven case well-captures the storm-time thermospheric neutral density response and its IHA on the dusk side. Basically, the density response has some properties associated with IHA: the start time of the enhancement, the relative magnitudes, and the latitudinal propagations between the NH and SH. It is found that the IHA in the neutral density response follows the variations in total Joule heating and has a time delay of~3 h. 5) Previous studies have shown that both the season and geomagnetic field can determine IHAs in the thermosphere.
Our study suggests that Joule heating associated with highlatitude forcing can become the dominant factor even during a moderate storm period when near the equinox.

Data availability statement
We acknowledge the use of NASA/GSFC's Space Physics Data Facility's OMNIWeb (or CDAWeb or ftp) service, and OMNI data. The IMF and solar wind data can be found at (https://omniweb.gsfc. nasa.gov/). The AMPERE FAC data can be found at (https://ampere. jhuapl.edu/). The GOCE neutral density data can be found at (https://earth.esa.int/eogateway/missions/goce/data). The DMSP SSIES ion drift data are available at http://cedar.openmadrigal. org/. The data used in this study along with AMIE outputs and GITM simulation outputs can be found at https://zenodo.org/ record/7475704#.Y6UlZhNKhR4.

Author contributions
YH conducted the data analysis and modeling studies. YD guided each step of this work and contributed to individual discussions of results. QZ coupled the NCAR-3D model into the GITM, added the FAC-driven capability to the NCAR-3Dynamo model, and created the AMIE electron precipitation pattern. AM developed the NCAR-3Dynamo module and the FAC-driven procedure. MH developed the identification of determining the equatorial boundary of ion convection. CW processed the Iridium satellite data and provided advice on differences between the hemisphere data fits. CS provided interpretation for the thermospheric neutral density performances. DW and RL assisted in the analysis of the storm event as experts in the research field of IHA and magnetosphere physics. All the co-authors contributed to the discussions, comments, and improvements to the manuscript.

Funding
This research conducted at the University of Texas at Arlington was supported by AFOSR through award FA9559-16-1-0364 and NASA through awards 80NSSC20K0195, 80NSSC20K1786, 80NSSC20K0606, 80GSFC22CA011, and 80NSSC22K0061. YH was also supported by the National Center for Atmospheric Research (NCAR) Advanced Study Program (ASP) Graduate Visitor Program Fellowship. QZ was supported by the NCAR ASP Postdoctoral Fellowship and NASA GOLD ICON Guest Investigators Program under grant 80NSSC22K0061 through the subaward 2021GC1619. YH and QZ were also supported by the NCAR Early Career Scientist Assembly Visitor Fund. AM was supported through AFOSR award FA9559-17-1-0248 and NASA award 80NSSC20K1784. This material is based on the work supported by the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under Cooperative Agreement No. 1852977.