Time-dependent responses of the neutral mass density to fixed magnetospheric energy inputs into the cusp region in the thermosphere during a period of large IMF BY: a high-resolution two-dimensional local modeling

Thermospheric mass density values around the 400-km altitude in the cusp can be significantly enhanced as compared to regions around the cusp. To gain insights into the extent to which the magnitude of the cusp mass density enhancements can be explained by the static distributions of moderate electric field and electron precipitation typical for a period of large IMF BY, we employed a high-resolution two-dimensional local model that can represent the plasma features that are characteristic of the cusp: azimuthal ion flow and low-energy electron precipitation. We also calculated the thermospheric dynamics with and without neutral–ion drag. We found that in the calculation with this drag the obtained mass density enhancement is 10% at most, indicating that the thermospheric dynamics imposing the moderate static electric field and electron precipitation can only explain about one-third of the typical magnitude of cusp thermospheric mass density, i.e., 33%. We also found that in the calculation without neutral–ion drag the magnitude of the mass density enhancement is slightly larger than the one with the neutral–ion drag. To explain the average magnitude of the cusp mass density enhancements completely, other energy inputs such as Alfvén waves, in addition to the static distributions of electric field and electron precipitation, are needed.


Introduction
The cusp is typically located around 75° magnetic latitude between ~ 1000 and ~ 1400 magnetic local time (MLT) in the altitudes of the ionosphere-thermosphere, where low-energy electrons almost directly come from the dayside magnetosheath. The ionosphere-thermosphere dynamics in the cusp region are strongly dependent on the condition of interplanetary magnetic fields (IMFs) and the condition of solar extreme ultraviolet (EUV) radiation as well as the inputs from the dayside magnetosheath (e.g., solar wind particles). Recent CHAllenging Minisatellite Payload (CHAMP) satellite observations have shown that the neutral mass density around the 400-km altitude in the cusp is remarkably larger than that of ambient regions. After Lühr et al. 's (2004) discovery, the anomalous mass density structure and other related phenomena have been extensively investigated in observational and modeling studies. Kervalishvili and Lühr (2013) have shown that the mass density enhancement is, on average, 33%.
The mass density anomaly is considered to be generated by thermospheric heating processes, such as Joule heating and particle heating, which drive neutral upwelling. Joule heating is caused by Pedersen currents depending on Pedersen conductivities, perpendicular electric fields, and neutral winds. Electron precipitation enhances Pedersen conductivity and, thus, Joule heating by ionization. In addition, electron precipitation directly heats the neutral atmosphere. The variability of the perpendicular electric field is also important for the enhancements of the Joule heating (e.g., Matsuo et al. 2003;Matsuo and Richmond 2008;Deng et al. 2009;Zhu et al. 2018).
Many modeling studies have been conducted to reproduce the mass density anomaly. Under geomagnetically disturbed conditions, previous studies have partially been successful in reproducing the mass density anomaly, while still facing difficulties in reproducing sufficient mass density enhancements under quiet conditions. Crowley et al. (2010) calculated mass density enhancements during strong IMF B y (+ 20 nT ) conditions and managed to generate mass density enhancements of over 200%. Wilder et al. (2012) also employed the Thermosphere-Ionosphere-Mesosphere-Electrodynamics General Circulation Model (TIME-GCM) (Roble et al. 1988;Roble and Ridley 1994) and created a density change of over 100% for a geomagnetic storm. Drawing on Ridley et al. 's (2006) Global Ionosphere-Thermosphere Model (GITM), Deng et al. (2013) imposed an intense Poynting flux of 75 mW/m 2 and soft particle precipitation, and consequently obtained mass density enhancements of more than 50%. Crowley et al. (2010), Wilder et al. (2012), and Deng et al. (2013) all assumed highly strong disturbances. Brinkman et al. (2016), who used the Aerospace Dynamic Model (ADM) (Walterscheid and Schubert 1990), obtained the result consistent with observation as far as the model of moderate energy inputs is concerned. However, for the strongest energy inputs arising from geomagnetic storms, their model's density enhancements were well below the observations. Overall, the complete explanation of the mass density anomaly for various conditions remains to be established.
Many previous studies have used global models such as TIME-GCM and GITM, most of which have a horizontal resolution of around 100 km and a vertical resolution of one-half scale height in TIME-GCM (or one-third scale height in GITM). This horizontal resolution is comparable to the typical latitudinal size of the cusp, and further studies using a higher spatial resolution model need be done when we consider the latitudinal profiles of the ion flow and electron precipitation which show prominent features in the cusp. The recent localized version of GITM can also reach the high resolution of 8 km (Lin et al. 2017). However, few studies still have used such high-resolution models to investigate the neutral mass density anomaly (Lin et al. 2017 focused on acoustic-gravity waves). Furthermore, a recent study has suggested that the time for the neutrals to respond to changes in the ionospheric plasma convection can be significantly reduced, and that the reduced time dampens Joule heating (Billet et al. 2020).
The purpose of this study is to gain insights into the extent to which the magnitude of the cusp mass density enhancements can be explained by the moderate static electric field and electron precipitation by employing a high-resolution (kilometer-scale) two-dimensional local model that considers time-dependent responses of neutrals and ions. We acknowledge the importance of the electric field variability on the energy inputs, but in this study we aim to understand to what extent the static electric field can explain the mass density anomaly.

Neutral dynamics
We consider the neutrals to be composed of N 2 , O 2 , NO , N , O , and He . The continuity equation of neutrals is where n j is the number density of neutral species j , and u n is the neutral flow velocity. The subscript n denotes the neutral. All neutral species are assumed to move at the same velocity.
Including Coriolis force, pressure, collision to ions, gravity, and viscosity, the momentum equation of neutrals can be written as: where Ω , ρ n , p n , ν ni , u i , G , and η are the angular velocity of the Earth's rotation, the neutral mass density, the neutral total pressure, the neutral-ion collision frequency, the ion flow velocity, the gravitational acceleration, and the dynamic viscosity, respectively. The subscript i denotes the ion. The collision frequency is derived by the formulas in Schunk and Nagy (2009). The dynamic viscosity is in Banks and Kockarts (1973).
Considering adiabatic expansion, heat conduction, and external heating, the energy equation of neutrals is given by where T n , R , c v , κ , Q are the neutral temperature, the specific gas constant, the specific heat capacity at constant volume, the heat conductivity, and the heating rate per unit volume, respectively. The heat conductivity is given by Banks and Kockarts (1973). In the auroral region, the heat sources are mainly Joule heating ( Q J ), particle heating ( Q P ), and viscous heating ( Q V ). Thus, Q = Q J + Q P + Q V . Defining E as the electric field and B as the magnetic field, the first Q J is given by where σ P is the Pedersen conductivity. The last Q V is given by where u , v , and w are the components of the neutral velocity of x, y, and z, respectively. The second Q P will be described in "Electron precipitation" section.
In general, there are many other heating and cooling processes, such as neutral-electron frictional heating, inelastic collisions, or infrared cooling (Schunk and Nagy 2009). However, they are very small compared with the three major processes above for typical conditions in the cusp. For example, when the electron temperature is 2500 K as a typical value in the cusp (e.g., Prölss 2006;Taguchi et al. 2017), the magnitudes of Joule heating, particle heating, neutral-electron frictional heating, and N 2 rotation cooling at 300 km altitude are the order of 10 −9 , 10 −10 , 10 −13 , and 10 −14 W/m 3 , respectively [these values are derived by formulas in Schunk and Nagy (2009)]. It is assumed that the initial state of the neutral atmosphere set by NRLMSISE-00 model (Picone et al. 2002) is in equilibrium (Shinagawa and Oyama 2006). That is, in the energy equation solar heating in the thermosphere is considered to be balanced with cooling by heat conduction initially. In the model, thermospheric motion is driven by additional processes such as Joule heating, particle heating and ion-neutral drag force.

Ion dynamics
We consider the ions to be composed of O + 2 , N + 2 , NO + , and O + . The continuity equation of ions is where n k is the number density of ion species k . S k is the net source term by ionization, recombination, and other chemical reactions. This term will be described again in "Electron precipitation" section.
Assuming time derivative, advection, Coriolis force, and viscosity to be zero, the momentum equation of ions is where ρ i , p i , e , m i , and ν in are the ion mass density, the ion total pressure, the elementary charge, the mean molecular mass of ions ( m i = ρ i /n i , where n i is the ion number density), and the ion-neutral collision frequency, respectively. (5) In the perpendicular direction of (7), Lorenz and collisional forces are dominant. Thus, we can derive the component perpendicular to the magnetic field of u i (i.e., u i⊥ ) as follows: The motion of ions to the component parallel to the magnetic field is determined by ambipolar diffusion. Assuming all ion species have the same temperature, u i is given by where k B , T p , D a , and H p are the Boltzmann constant, the plasma temperature, the ambipolar diffusion coefficient, and the plasma scale height, respectively. T e is the electron temperature, T i is the ion temperature, and I is the geomagnetic inclination. In this study, I is set to be 90 • since we consider idealized and simplified situations.
The ion temperature is approximated well by assuming a balance between the frictional heating and the heat (8) Fig. 1 Schematic illustration that shows the geometry of the simulation. The x-axis, y-axis, and z-axis are directed eastward, northward, and upward. The peak of electron energy flux ( Q 0 ) is located at the center ( y = 0 km). The peak of southward electric field ( E y ) is shifted from the center to the north to maximize the upward field-aligned current at the center exchange with neutrals as follows (St.-Maurice and Hanson 1982): where m n is the mean molecular mass of neutrals ( m n = ρ n /n n , where n n is the neutral number density).

Electron precipitation
To describe the effects of electron precipitation, we employed Fang et al. 's (2010) empirical model, which derives the altitude profile of the total ionization rate due to electron precipitation. We assumed the differential number flux of precipitating electrons φ(U ) as a kappa distribution as follows: where Q 0 is the total energy flux, and U 0 is the characteristic energy. The peak altitude of the ionization increases for lower characteristic energies. In the cusp region, the electron precipitation is characterized by "soft" (~ 100 eV) electrons coming almost directly from the magnetosheath. Another effect of electron precipitation is particle heating. Precipitating electrons collide with neutral molecules and transfer energy. Some energy is lost by dissociation and radiation, and the rest eventually heats the molecules. Using the total ionization rate P /m 3 s , the particle heating rate Q P eV/m 3 s is given by where �ε = 35 eV is the mean ionization energy, and C eff is the heating efficiency as an empirical function of height (Rees et al. 1983;Richards 2013).

Numerical implementation
As can be seen in many convection models such as the Heppner and Maynard model (1987, for example their Figures 1, 2, and 3), the east-west dimension of the relatively strong electric field near the cusp is about 3 MLT. This is roughly 1400 km at the altitude of 300 km (assuming at 75 magnetic latitude). Since the typical latitudinal (north-south) size of the region where such a northward/ southward electric field occurs is roughly 200 km, the gradient in the east-west direction is sufficiently small when compared with that in the north-south direction. These indicate that the basic feature of the cusp can be represented by a two-dimensional model (meridian vs. (11) Q P = �εPC eff , altitude). We note that southward or northward electric field can more often occur in the cusp than dawn-to-dusk electric field because IMF tends to form a spiral "gardenhose" configuration.
We have developed a new two-dimensional local model, which is based on Shinagawa and Oyama's (2006) thermospheric neutral model. Figure 1 is a schematic illustration that shows the geometry of the simulation. We set the x-axis, y-axis, and z-axis to be directed eastward, northward, and upward. All physical quantities are assumed to be uniform in the x-direction. The numerical domain ranges from 0 to 700 km in altitude and from − 4000 to 4000 km in meridional distance. We separate the domain into cells with a vertical size of Δz = 5 km and a horizontal size of Δy = 10 km. The time step Δt is set to be 1 ms. We employ the cubic-interpolated pseudoparticle (CIP) method (Takewaki et al. 1985; to obtain the time evolution.
We set boundary conditions as follows: At y = − 4000, 4000 km, ∂f /∂y = 0 for any physical quantity f . At z = 0 km, u n = 0, ∂T n /∂z = 0, and ∂n n /∂z = ∂n i /∂z = 0 . At z = 700 km, ∂u n /∂z = 0, ∂T n /∂z = 0, ∂n n /∂z = −n n /H n , and ∂n i /∂z = −n i /H p . The last two conditions mean diffusion equilibrium for neutrals and ions. The precipitating electron energy flux Q 0 and the southward electric field E y are set to be Gaussian functions as follows: where W is the scale width and set to be 200 km as a typical meridional width of the cusp. The peak of electron energy flux is located at the center ( y = 0 km ). The electric field peak is shifted from the center to the north to maximize the meridional gradient of E y at the center, which reproduces the situation that cusp electron precipitation tends to occur at the upward field-aligned current for IMF B y < 0 (i.e., lower-latitude part of the two-sheet cusp field-aligned currents), as has been shown by Taguchi et al. (1993).
In this study, we set the peak electric field E to be 60 mV/m . The electric field of 60 mV/m is based on the satellite observation of the cusp (Maynard et al. 1982), which may be regarded as a moderate condition. The electric field in the cusp can be variable during highly disturbed interval (i.e., strongly southward IMF) or during the interval when the direction of the IMF often changes. However, when the direction of the IMF is stably eastward or westward, which usually occurs because IMF tends to form a spiral "garden-hose" configuration, Fig. 2 The resulting north-south profiles of the fractional density change and neutral velocity (top), neutral temperature change (middle top), eastward neutral velocity (middle bottom), and specific heating rate (bottom) in each case. In the top, the contour maps show the fractional density change, and vectors show the neutral flow velocity. The right side (positive values) of the horizontal axis is the north. The peak of electron precipitation is located at the horizontal center there should be southward or northward electric field in the cusp depending on IMF B y < 0 or B y > 0 , respectively (IMF B y < 0 is assumed in this study as mentioned above). The simulation conducted in our study is aimed at revealing the effect caused by this sort of the stable electric field being imposed. This is the reason why we adopt the static southward electric field in the simulation. Electron precipitation is imposed with total energy flux Q of 1.6 mW/m 2 and characteristic energy U 0 of 100 eV , which indicates "soft" electrons into the cusp (Newell and Meng 1988).
The initial condition of neutrals is set by the NRLM-SISE-00 model (Picone et al. 2002)  Since the original profile given by IRI-2016 is not in equilibrium, we first run the model for 6 h without any external forcing. The resulting ion profile is used as the initial condition in the following calculations. The electron temperature is fixed to the value of the IRI-2016.

Modeling runs
We performed two modeling runs to investigate the contributions of neutral-ion drags to the neutral atmosphere. We calculated with neutral-ion drags in Case 1 and without them in Case 2 for comparison. In Case 2 the collisional term of the momentum Eq. (2) −ν ni (u n − u i ) was dropped (i.e., ν ni = 0 ). Note that Joule heating, as well as drag forces, is generated by collisions between neutrals and ions, and that the approach of turning off the neutral-ion drag term in the momentum equation but keeping Joule heating term in the energy equation may cause inconsistency.
The outline of the procedure of the calculation is as follows: 1. Two-dimensional profiles of neutrals and ions are set up. 2. Electric field and electron precipitation start at t = 0. 3. The time evolution of neutrals and ions is calculated over 7200 s (2 h).

Comparing the contributions of various ionospheric processes
We define fractional density change as �ρ/ρ 0 , where ρ 0 is the initial neutral mass density, and �ρ is the difference of neutral mass density from ρ 0 . Therefore, �ρ/ρ 0 indicates the relative enhancement of mass density. For instance, �ρ/ρ 0 = 0.1 means a 10% increase from the initial condition. We refer to heating rate per unit volume as volumetric heating rate, and heating rate per unit mass as specific heating rate. What matters for the neutral response is the latter. Figure 2 shows the resulting north-south profiles around the center at t = 7200 s. The right side (positive values) of the horizontal axis is the north. The contour maps in Fig. 2a show the fractional density change, and vectors show the neutral flow velocity. Figure 1b-d shows the neutral temperature change, eastward neutral velocity, and specific heating rate, respectively. In altitudes of 200 km to 400 km, the specific heating rate is maximized, and then neutral air heats, which causes neutral upwelling and mass density enhancements. Figure 2 shows that the mass density enhancement, upward neutral velocity, and specific heating rate of Case 1 are all smaller than those in Case 2. The peak values of mass density changes, neutral temperature changes, and vertical neutral velocity are summarized in Table 1. The peak of mass density is located south from the center in Case 1. When neutral-ion drags are present, the neutral air is pulled into the direction of the eastward E × B drift (Fig. 2c). After that, the Coriolis force moves the neutrals southward, causing the large mass density in the southern region. The differences in peak locations between the electric field and electron precipitation also cause weak asymmetry. Figure 3 shows the neutral atmosphere profiles at the 400-km altitude. Figure 3a-e shows the mass density enhancements, temperature changes, eastward velocity, northward velocity, and upward velocity. Similar to Fig. 2, the fractional density change, neutral temperature change, and vertical neutral velocity in Case 1 are all smaller than those in Case 2. The mass density on the southern side in Case 1 is larger than that on the northern side. Figure 3c shows that the large eastward neutral  . 3 The resulting profiles of the fractional density change (top), neutral temperature change (middle top), eastward neutral velocity (middle), northward neutral velocity (middle bottom), and vertical neutral velocity (bottom) at the 400-km altitude in each case. Positive distance is directed to the northward. The peak of electron precipitation is located at the center velocity is caused by the neutral-ion drag, as expected. Figure 3d shows that the large eastward flow suppresses the northward flow via the Coriolis force in the north region ( y ≥ 250 km) in Case 1. The peak values of the mass density changes at the 400-km altitude are 9.7% in Case 1 and 12.3% in Case 2. The peak values of mean mass density are about one-third of the mean value, i.e., 33% (Kervalishvili and Lühr 2013). Figure 4 shows the time evolution of mass density changes at the 400-km altitude and the specific heating rate at the 300-km altitude, where Joule heating drives neutral upwelling most effectively. Figure 4a shows that the neutral mass density oscillates until about 50 min due to atmospheric gravity waves caused by sudden commencement of heating at 0 min (in our calculations, electric fields and electron precipitation rise as step functions at the beginning). After 50 min, the mass density increases very slowly in both cases. Figure 4b shows the time evolution of the specific heating rate. In both cases, Joule heating initially increases by ionization and then decreases for several tens of minute at least. The specific heating rate at the 300-km altitude is initially 99 W/kg and eventually reaches 140 W/kg in Case 1, but when we do not consider the neutralion drag (Case 2), a much greater specific heating rate ( 245 W/kg ) is produced. We also note that in both cases the Joule heating rate decreases with time, which is why mass density enhancements do not keep on increasing even for longer times of energy inputs. Figure 5 shows altitude profiles of the ion density, specific heating rate, ion temperature, and vertical ion velocity at the center at 40-and 120-min intervals in each case. Assuming u i = E × B/B 2 , the volumetric Joule heating rate (4) can be written as follows:

Effects of each process on the time evolution of Joule heating
Thus, larger velocity differences between neutrals and ions generate larger Joule heating rates. In Case 1, the horizontal neutral velocity is pulled into the E × B drift direction (x-direction) and finally reaches a value where the neutral-ion drag and viscous forces are balanced (In this direction, the Coriolis force is tiny, and the pressure gradient is zero because of the two-dimensional assumption.), showing that the neutral-ion drag can reduce the Joule heating. Figure 4b, c shows that both ion temperature and specific heating rate of Case 2 are larger than those of Case 1 since large velocity differences between neutrals and ions are maintained in Case 2.
In the F layer, the major ion species is O + , and the dominant chemical reactions are: Neutral upwelling brings molecule-rich air to higher altitudes (Fuller-Rowell et al. 1996;Lu et al. 2016). Thus, the three chemical reactions above all act to decrease the O + density in the F layer (Fig. 5a). Figure 5d shows ion down-flow at 120 min, corresponding to the reduction of ions.
The results indicate that both neutral-ion drags and chemical reactions in this model reduce Joule heating (16) The resulting profiles of the ion density (left), specific heating rate (middle left), ion temperature (middle right), and vertical ion velocity (right) at 40 min (top) and 120 min (bottom) intervals in each case rates, and that the enhancement of the neutral mass density is overestimated when the fixed ionosphere condition is adopted.

Comparison with previous studies
In our calculation, we used the input of the electric field of 60 mV/m representing a moderate condition. The direction of the electric field was assumed to be southward, typical for large negative IMF B Y , and the peak of cusp precipitating electron energy flux was set to be located at the lower-latitude part of that electric field which corresponds to the upward field-aligned current in the cusp. One of the unique features of our modeling is to incorporate these distributions as inputs. While these distributions are consistent with the observation result on a few 100-km scale, our numerical simulation has reproduced only about one-third of the average observation.
We also performed several modeling runs, and found that the electric field of approximately 150 mV/m is needed to produce the typical mass density enhancement, i.e., ~ 33%. In contrast, Brinkman et al. (2016) obtained the neutral mass density enhancement of more than 30% on the assumption similar to the one used in our Cases 1 and 2, i.e., the total electric field of 60 mV/m , and the total energy flux of 1.6 mW/m 2 with a characteristic energy of 100 eV . This difference is reasonable because the enhancement of the neutral mass density is overestimated when the fixed ionosphere condition is adopted as mentioned in "Effects of each process on the time evolution of Joule heating" section. Clemmons et al. (2008) have reported the mass density depletions at 123-325 km altitudes by the Streak mission. Figure 6 shows the fractional density change at 250 km altitude in Case 1. The difference between the maximum and minimum is 3.4%. The mass density depression reported by Clemmons et al. (2008) is a few percent, which is consistent with our results.
As has been stated in "Numerical implementation" section, our numerical model has been run at a horizontal resolution of 10 km. This is higher than previous models. ADM model is run at a resolution of 20 km. TIME-GCM is run at a resolution of 5° in longitude and 5° in latitude. GITM model in Deng et al. (2013) is run at a resolution of 2.5° in longitude by 0.375° in latitude. Although Lin et al. (2017) have used the localized version of GITM with a high-resolution of 0.08° in longitude and latitude, their target is the acoustic-gravity wave (not the mass density anomaly). Some previous studies have reported that electric field variability can play an important role in Joule heating (Deng et al. 2009;Zhu et al. 2018). This variability is driven by electromagnetic waves (typically Alfvén waves). Since alternating electric fields keep velocity differences between neutrals and ions large, the mass density's time evolution may differ from what we presented above. Lotko and Zhang (2018) have shown that volumetric Joule heating rates generated by Alfvén waves are maximized at F layer altitudes, with the altitude profile depending on wavelength and frequency. In contrast, quasi-static electric fields maximize volumetric Joule heating rates at E layer altitudes. This means that Alfvén waves can generate very large specific heating rates at F layer altitudes, which would play an important role in the mass density anomaly.
In this study, we have focused on what extent the static electric field can explain the mass density anomaly, and clarified that other energy inputs different from those related to the static distributions of electric field and electron precipitation are absolutely needed. When we consider the above-mentioned recent results on the importance of the volumetric Joule heating rates generated by Alfvén waves, a strong candidate for additional energy sources will be the Alfvén wave heating. For an in-depth discussion, a precise calculation of the height profile of Alfvénic heating needs to be done.

Summary
We used a high-resolution numerical model to investigate to what extent the moderate static southward electric field, which is typical for a period of negative IMF B Y in the cusp, and the soft electron precipitation distribution can explain the mass density anomaly. Contributions of neutral-ion drags were compared using two cases with and without neutral-ion drags in the momentum equation. We found that neutral-ion drag forces decrease velocity differences between neutrals and ions. Chemical reactions reduce ions at F layer altitudes in response to neutral upwelling. The mass density enhancement in the calculation containing the neutral-ion drag process is 10% at most, indicating that the thermospheric dynamics imposing the moderate static electric field and electron precipitation can only explain about one-third of the typical magnitude of cusp thermospheric mass density. To explain the average magnitude of the cusp mass density enhancements completely, other energy inputs such as Alfvén waves, in addition to the static distributions of electric field and electron precipitation, are needed.