Global Variations in the Time Delays Between Polar Ionospheric Heating and the Neutral Density Response

We present results from a study of the time lags between changes in the energy flow into the polar regions and the response of the thermosphere to the heating. Measurements of the neutral density from the Challenging Mini‐satellite Payload (CHAMP) and Gravity Recovery and Climate Experiment (GRACE) missions are used, along with calculations of the total Poynting flux entering the poles. During two major geomagnetic storms in 2003, these data show increased densities are first seen on the dayside edge of the auroral ovals after a surge in the energy input. At lower latitudes, the densities reach their peak values on the dayside earlier than on the night side. A puzzling response seen in the CHAMP measurements during the November 2003 storm was that the density at a fixed location near the “Harang discontinuity” remained at unusually low levels during three sequential orbit passes, while elsewhere the density increased. The entire database of measurements from the CHAMP and GRACE missions were used to derive maps of the density time lags across the globe. The maps show a large gradient between short and long time delays between 60° and 30° geographic latitude. They confirm the findings from the two storm periods, that near the equator, the density on the dayside responds earlier than on the nightside. The time lags are longest near 18–20 hr local time. The time lag maps could be applied to improve the accuracy of empirical thermosphere models, and developers of numerical models may find these results useful for comparisons with their calculations.

, Oliveira and Zesta (2019), and others. Lu et al. (2014) had also reported that traveling atmospheric disturbances (TADs) spread the energy around the globe, and that "the peak of neutral temperature and mass density enhancement lags behind the peak of Joule heating by 3-5 hr." In this paper, we present results from a detailed study of the time delays between auroral heating in the polar ionospheres and the response of the upper atmosphere neutral density. We start with an examination of the variations during two particularly large heating events in October and November 2003, using density measurements taken two satellites. Graphs of density as a function of both time and the "argument of latitude" (https://en.wikipedia. org/wiki/Category:Orbits) in the orbit plane are shown first. These graphs are followed with a comparison with maps of the modeled electric potential patterns in the polar regions, with the measured densities shown on superposed orbit tracks. Finally, maps are derived, from a database spanning decades, that illustrate how the lag times between ionospheric heating and density response varies between different locations around the globe. These results can be used to improve both empirical and numerical models of the thermosphere.

Density Variations Observed During Extreme Heating Events
The time period in October and November 2003 is well known for two exceptionally large geomagnetic storms that had significant heating in the polar aurora, as well as strong magnetic variations (Gopalswamy et al., 2005). We start with the second storm starting on 20 November 2003, with density measurements shown as a function of time in Figure 1. These data are derived from accelerometer measurements of the orbital drag forces taken on the Challenging Mini-satellite Payload (CHAMP) satellite (Bruinsma et al., 2004). The density values had been recalculated by Mehta et al. (2017), and are available with temporal cadence of 10 s for CHAMP. It is evident that significant variations in the density start around 8 UT on the 20th.
For the purpose of studying the response times of the density variations, it is more useful to graph the data in the format shown in Figure 2. For reference, Figure 2a at the top shows the total Poynting flux entering the Northern (red) and Southern (blue) hemispheres as a function of time. These values are calculated from the model described by Weimer (2005b), referred to as W05, which described an improved version of the one that first described how the Poynting flux values can be obtained with an empirical model (Weimer, 2005a). The input values for this model were obtained from the interplanetary magnetic field (IMF) and solar wind velocity measurements obtained with the Advanced Composition Explorer (ACE) spacecraft, that were propagated to the Earth's bow shock using the time delay calculations described by Weimer and King (2008). After propagation, the inputs for the W05 model at a given time were obtained from averages of each IMF vector component and solar wind velocity over the previous 20 min. This averaging introduces a delay between changes in the IMF and the output of the Poynting flux calculation to account for the reconfiguration time of the electrodynamic patterns. Normally the heating calculations would be delayed an additional 20-35 min to account for the propagation of the signal from the bow shock to the polar regions (Weimer et al., 2010). This uncertain propagation delay was not added in this case so that the propagation times will be included in the time lags that we obtain. Figure 2b shows the density values color coded along the orbit track in reference to the color bar shown at the bottom. The densities are graphed as a function of both time and the argument of latitude, which is an orbital parameter that starts at zero when the spacecraft crosses the geographic equator while ascending (increasing in latitude), and smoothly increases until reaching 360 when it returns to the equator on the next orbit. (Geographic latitude is not used since it is not a continuous function for orbits that do not pass exactly over the poles, and the latitude alone does not uniquely identify the position.) The local time of the equatorial crossings are indicated on the right axis. The widths of the color-coded density measurements are greatly expanded for clarity. The thick blue lines show where the satellites pass through 60° magnetic latitude, the approximate, lower boundary of the auroral oval with moderate geomagnetic activity. The actual location of this boundary moves to higher and lower latitudes as conditions vary. As indicated with the labels, the Northern oval passes are at the bottom of the graph, while the Southern oval passes are at the top. Other labels indicate where the CHAMP satellite is on the dayside and nightside. At the left side of the graph the dayside is evident by the higher densities with the lighter color.
Very prominent in Figure 2 are the density increases that take place when the auroral heating increases, followed by a gradual decrease, or cooling, after the heating is reduced. Examining how the densities vary in more detail, it is evident that the effects of the increased heating are first seen on the dayside edge of both auroral ovals, particularly in the South before 1430 UT on the 20th. On the dayside, the densities appear to reach a peak between 19 and 2030 UT, as indicated by the label "A." On the nightside, the peak density is reached later, as indicated with the "B" around 20:30 UT; placing this mark at one orbit later is possible, as in this and subsequent graphs the placement of these labels is a matter of judgment.
Density variation in the auroral oval are more erratic. As mentioned earlier, density increases are first seen on the dayside edges of the auroral oval, but they are not persistent. The largest density peak occurred in Northern oval after 1900 UT, but was not seen on the next orbit. One pass through the Southern oval before 0000 UT on the 21st has a higher density than during the prior orbit. Later on in the Southern oval and on the nightside, density values decreasing sooner than elsewhere are seen before 1200 UT on the 21st.
One of the most curious features of these measurements on CHAMP are the persistent depressions in the density that are seen in Figure 1 between 12 and 16 UT, during a time when the density is increasing elsewhere. Figure 2 shows that these dips are located near the nightside edge of the auroral oval in the Southern hemisphere, and other depressions are seen within the oval. We will return to this topic in Section 3. During this time period density measurements in a different orbit plane were obtained with the Gravity Recovery and Climate Experiment (GRACE) satellites (Tapley et al., 2004). The GRACE orbit is approximately 100 km higher than CHAMP at this time. There are two GRACE satellites that orbit very close together and measured nearly the same density values, so only data from the "A" satellite are shown here. We use the density values calculated by Mehta et al. (2017), that are available with temporal cadence of 5 s for GRACE. Figure 3 shows the density measurements from the GRACE A satellite for the same time period in November 2003, in the same format as Figure 2.
GRACE reaches the ascending node (equator) on the post-noon dayside at 14.50 hr local time, rather than on the nightside at 23.21 hr local time for CHAMP. These orbits result in differences in the overall appearance between Figures 2 and 3, but the general features are the same. Both the dayside and nightside densities increase when the auroral heating is strong, and decrease afterward. The densities on the dayside and nightside reach their peak values on the same orbit, around 1930 UT (label "A") to 2000 UT ("B"). Exact timings are difficult to determine due to the fact that each location is sampled only once per orbit. The auroral oval densities detected with GRACE fluctuate less than in the CHAMP data, and also lack the deep, repeated troughs.
Figures 4 and 5 show the CHAMP and GRACE density measurements in the same format for the earlier storm that occurred on 29-30 October 2003. This earlier storm had a more complex time variation, consisting of a sudden impulse in the heating after 0600 UT that exceeded 2.5 TW in both hemispheres, followed by a lull before an extended period of heating between 1800 UT on the 29th and 0300 UT on the 30th. The CHAMP results in Figure 4 show sudden increases in the density on the dayside boundary of the auroral oval in both hemispheres in the orbit that begins right after the power spike at 0600 UT, and a density spike is also seen at the equator just before 1000 UT, at 1.26 local time. Figure 5 shows that GRACE detected a density spike on the dayside oval boundary in the southern hemisphere around 0700 UT, immediately after the spike in the heating had occurred, and a similar density increase was seen on the next orbit. The first density spike was not seen in the northern hemisphere, as the GRACE orbit had already passed through the northern oval before the heating peaked, but then a density bump was seen in the northern, dayside oval on the following orbit. On the orbit after that, a density spike was seen near the equator at 4.17 local time.
The CHAMP measurements in Figure 4 show that the densities on the dayside peaked around 2200 UT ("A") at 13.26 local time, which the peak on the nightside near 1.26 local time did not reach their maximum until 0100 UT on the next day ("B"). On the other hand, the GRACE measurements in Figure 5 show that the dayside and nightside densities reached their peaks on the same orbit. In this case, the orbit plane passed through the dayside later in the afternoon, at 16.18 local time, and the nightside part of the orbit was in the early morning at 4.17 local time. In both Figures 4 and 5 the density variations within the auroral ovals were as erratic as in Figures 2 and 4, with unusually low or high density values appearing seemingly at random.

Polar Cap Passes
The variability of the density within the polar regions has puzzling features. To see if there is any relationship between these variations and the auroral electric fields and plasma convection patterns, the densities along the orbit paths were plotted over contour maps of the electric potentials from the W05 model. An example of CHAMP measurements is shown in Figure 6 for CHAMP measurements for part of the event that begins on 20 November 2003. The remainder of this event is included in the accompanying Supporting Information S1, which includes similar graphs with both CHAMP and GRACE densities for the two geomagnetic storms featured in Figures 1-5. In these figures, the electric potentials and the orbit tracks are drawn in corrected geomagnetic apex coordinates  (Emmert et al., 2010;Richmond, 1995;VanZandt et al., 1972). The label above each contour plot indicates the date and time when the satellite reaches the highest (absolute value) latitude in apex coordinates. The electric potential patterns were calculated for this same point in time. Level 2 science data from the ACE satellite were used to obtain the IMF and solar wind velocity values required for the W05 model, using time-shifted, 20-min averages. The solid blue line around each contour map marks the location of the lower boundary where the electric potential reaches zero. The configuration of the convection pattern may change while the satellites pass over the pole. Ionospheric plasma in the ionosphere follows the contour lines, moving from noon to midnight at the center, and flowing in a sunward direction at lower latitudes.
It is important to keep in mind that, while the electric potentials and field aligned currents follow the magnetic field lines and map to the same corrected geomagnetic apex coordinates at different altitudes, the neutral atmosphere is not similarly constrained. So it is not expected that upward expansion of the neutral atmosphere will align exactly with the electric potential (plasma convection) patterns. On the other hand, it is known that neutral winds in the polar regions are normally dragged along in the same direction as the ions, which follow the equipotential field lines . It is seen in Figure 6 that the initial response of the neutral density in the polar regions is seen to occur near the dayside "convection throat" (Siscoe & Huang, 1985). Examples are in the plots at times after 08:36 UT and especially at 14:05 UT, but this cusp enhancement is not seen in every case.
Returning to the topic of the three sequential density depressions seen in the CHAMP measurements on 20 November 2003, these depressions appear in Figure 6 in the three Southern hemisphere polar graphs labeled 12:31, 14:04, and 15:37 UT. Since the CHAMP spacecraft was moving from the dayside toward nightside in these graphs, the depressions occur shortly after the times indicated on the graph labels. Curiously, these persistent density minima are found near the nightside "Harang discontinuity" (Gjerloev & Hoffman, 2001) where the 7 of 15 plasma flow changes from Eastward to Westward near local midnight. It seems that colder neutrals could be dragged in from the night side here, but since similar perturbations are not seen near midnight in the Northern hemisphere nor in other graphs in the Supporting Information S1, then this hypothesis does not seem likely.
There are other cases, such as in the map at 10:58 UT, where a lower density is found within the boundaries of a closed convection cell on the dawn side, as well as at other times illustrated in the Supporting Information S1, but there are other cases where this does not happen. Could the apparent density depressions and enhancements be caused by the effects of neutral "tail" and "head" winds, driven by the ion convection? Again, the results cannot confirm this as there is no consistent behavior indicated in these maps to support that possibility. In the example just mentioned near 10:58 UT the low density reported for the CHAMP satellite was located in sunward flow, while in the other cases the density depression was in anti-sunward flow.

Mapping the Time Delays
In Figures 2 through 5, we had seen cases where the neutral density on the dayside responds to increased polar heating faster than on the nightside, and one case where day and night densities reach their peaks at nearly the same time. There appears to be a local time dependence. As the satellites only sample each location once per orbit, exact timings are difficult to measure. So next we turn to statistical means to see how quickly the densities around the globe respond to the heating.
Weimer et al. (2020) described a method by which neutral density measurements from multiple satellites taken over a long period of time are sorted according to location on a geodesic grid. This grid has approximate equal areas in each cell. As measurements taken at different altitudes are combined together, the variation of density  with altitude needs to be accounted for. Using the technique described by Weimer et al. (2016), the density measurements were converted into exospheric temperature values. This conversion was done through use of the thermosphere density model, known as the Naval Research Laboratory Mass Spectrometer and Incoherent Scatter radar Extended (NRLMSISE-00) model (Hedin, 1991;Picone et al., 2002). In the Weimer et al. (2020) paper the accumulated temperature values in each grid cell were used to derive an expression for the exospheric temperature as a function of different quantities, such as solar indices and polar heating. Six different formulas were presented.
In the analysis results presented here the same technique was used except that the newer thermosphere model named NRLMSIS 2.0  was used to convert the satellite densities to exospheric temperatures. This newer model matched the satellite densities calculated by Mehta et al. (2017) better than before, so prior to the exospheric temperature calculation the CHAMP data did not require multiplication by a correction factor (Weimer et al., 2016).

Joule Heating Time Lags
One method used for finding the time lag between the heating and exospheric temperatures was to look at the correlation at different lag times, and select the lag that produced the highest correlation. The results showed promising results, with very short lags at high latitudes on the dayside (around 40 min), as expected, and longer lags at lower latitudes (over 200 min). However, a number of grid cells had extreme lags, over 400 min, and much larger than neighboring cells, particularly in the polar region right next to cells with the lowest lags. In taking a closer look at the data within these particular cells, it was found that graphs of correlation versus lag time had a bimodal response; the correlation peaked at a short time lag, but then the correlation decreased before increasing again to a slightly higher correlation at much longer delays. In these cases, the lower lag should be chosen as the response time.
Rather than looking at just a simple correlation, a better measurement time lags between the Joule heating and exospheric temperatures was accomplished while compensating for the effects of varying solar radiation, the day of year (subsolar latitude and other effects), and universal time. The same technique described by Weimer et al. (2020) is employed, in which the exospheric temperature data within each cell is fit to other parameters using a multiple linear regression fit. The main difference here is that variable time lags are incorporated in the fit. The formula used in the regression is: This formula is basically the same as (6) by Weimer et al. (2020), except for a rearrangement of the terms and incorporation of the variable time lags, as indicated with the (t − t lag ) notation. The symbols S 10 and M 10 represent the solar indices described by Tobiska et al. (2008), that are used in the 2008 Jacchia-Bowman (JB2008) thermosphere model . The S T values in this formula are the total, area integrated Poynting flux values from the W05 model, that were calculated using the solar wind and IMF data from ACE at 4-min steps. The S T values were smoothed using a 60-min boxcar function. The t lag values are referenced to the center of the smoothing window. Use of smaller smoothing windows reduced the correlations.
To generate a map of the delay lag times, the regression fits were repeated within every grid cell using a set of lag times that span a wide range. The least square errors were calculated for each lag. So that lower lag times could be found in the cases that had a double-peaked response, the search in each cell initially went up to only 96 min at 32-min steps. The bisection method (also known as the interval halving or binary search method) was used to resolve the lag time to 2 min: after identification of the lag corresponding to the lowest error, the search process was repeated with a range going from the best lag minus 32 min up to plus 32 min, at 16-min steps, and so on. If the lowest error was found at less than 96 min, then the search ended there, after the finer resolution search with the bisection method. If the best lag was found at the maximum end of the search range, then the next search went from there upward to 296 min, and if the lowest error was not found below 296 min, then the range was extended upward to 504 min. Figure 7a shows the resulting map of the optimal lag times. It is seen that in the polar regions, at latitudes over 60°, the lags are generally less than 60 min, and below 30° they are in the range of 180-300 min. There remain a few places around 30° latitude having optimal lags of 400-458 min, but not near local noon. As mentioned earlier, these lags include the time required for the electrodynamic patterns to respond after the signal in the IMF reaches the magnetopause. The mean correlation over the entire map, between the left and right sides of Equation 1, is 0.96.
There are statistical fluctuations in this map, and vertical striations caused by variations in solar activity while the CHAMP and GRACE orbits moved slowly in local time. For example, there could be flares and geomagnetic storms occurring during multiple orbits while a satellite orbit plane was near one particular local time. The activity could decrease after the satellite moved on to a slightly different local time, then ramp up again (often at 27-day intervals) after the orbit moved a little more. The striations in the maps are reduced using a spatial smoothing. Our smoothing filter averaged the contents of each grid cell with the neighboring grid cells that have a common vertex point. Three passes through the smoothing remove most of the vertical stripes, but with some artifacts still remaining. The result of the smoothing is shown in Figure 7b. As expected, the lowest lags are found in the polar regions, near local noon. There is a sharp gradient near 30° north and south latitude, and the largest lags are found here at more than ±5 hr away from local noon. Near the equator the lags around pre-midnight at 20 hr (approximately 280 min) appear to be larger than the post-midnight lags at 4 hr.

ΔT Time Lags
Another way to calculate the time lags is to use an integrated form of the polar heating referred to as ΔT, that rises in proportion to the heating and then falls at an exponential rate after the heating ceases. The use of this ΔT parameters is based the work by Burke et al. (2009), that was inspired by the findings of Wilson et al. (2006). of ΔT c uses only the geomagnetic Dst index. Weimer et al. (2011Weimer et al. ( , 2015 later showed that it is possible to use the W05 Poynting flux values to calculate the ΔT c values for the JB2008 model, improving the performance using IMF values rather than the geomagnetic Dst index. A similar ΔT is used in the model described by Weimer et al. (2020). The values as a function of time are calculated with the following numerical difference equation: In each time step ΔT increases in proportion (α) to the total Poynting flux in both hemispheres (S T ) and decays at an exponential rate with time constant τ c . The last term in Equation 2 simulates the simulates additional cooling from nitric oxide emissions, represented by P NO , that causes a more rapid decrease in the exospheric temperature following geomagnetic storms (Mlynczak et al., 2003). Referring to equations (10) and (11) by Weimer et al. (2020) for details, the levels of the simulated P NO increases in proportion to the Joule heating decays exponentially (Weimer et al., 2015). Figure 8 shows an example of ΔT that is calculated with Equation 2 for the 20-21 November 2003 time period, along with the total Poynting flux values in the Northern and Southern hemisphere.
The technique used to calculate the time lags between the ΔT parameter and exospheric temperatures is the same as in the prior section, except that ΔT is substituted for S T in Equation (1). In early results, it was found that in some grid cells the lowest error could be encountered at time lags or zero or less. At first it seems nonsensical that negative lag times could be obtained, but Figure 8 explains how this possible. As ΔT reaches a peak in this graph approximately 3 hr behind the peak in the total Poynting flux, the peak response to the Joule heating may occur before the peak in ΔT. Changes in the initial minimum and maximum search times were required before the start of the binary search method to accommodate the possibility of negative time lags.
The final result had a global mean value of 0.97 for the correlation between the exospheric temperatures and the regression fits. Figure 9a shows the resulting map of the optimal lag times for the ΔT parameter that were derived for each grid cell. Figure 9b shows results after running the map through the smoothing filter three times. This map shows that the response in the temperatures to the polar heating is most rapid near the polar cusp at 12 hr local time. The lags increases as the absolute value of the latitude decreases, becoming much longer at latitudes under 30°, with lags times over 140 min found at 20 hr local time near the equator, and a lag near 60 min at 4 hr local time.

Discussion
Figures 2-5 showed examples of density as a function of both time and latitude during two major geomagnetic storms that occurred in October and November 2003, from measurements on the CHAMP and GRACE A satellites, along with the total ionospheric heating on the same timeline. As expected, the density increases when extreme auroral heating occurs, followed by a gradual decrease after the heating ceases. The density increases are first seen on the dayside edge of the auroral ovals. At lower latitudes, the density usually reach their peak values on the dayside earlier than on the night side. Within the auroral oval, the variations are more irregular.
In the 29-30 October 2003 time period the event initially begins with an extreme level of heating that is both sudden and of short duration. Not surprisingly, the density response seen here was also fluctuating and uneven, undoubtedly due to the presence of TADs (Lu et al., 2014), coupled with the orbital characteristics of the satellites that prevents continuous measurements at each location. The initial impulse on the 29th was followed by a period of more continuous heating. At lower latitudes, the response was similar to the other event, with density reaching the peak value on the dayside earlier than on the night side.
One prominent result that is puzzling is the occurrence of a region of lower density that remains near the lower, pre-storm level while increasing densities are seen everywhere else. This behavior is illustrated in different ways in Figures 1, 2, and 6. The behavior of the density during the storms along the orbit path in relation to the electric potential/plasma convection patters are shown in Figure 6. While the results are generally inconsistent, there does seem to be a trend in that the prolonged density minima appears near the pre-midnight Harang discontinuity. As the neutral winds often match the motion of the plasma ions (Killeen & Roble, 1988), the colder neutrals could have been carried into this region from the pole.
There have been previous reports of low-density regions observed at high latitudes near dawn and dusk at around 200 km altitude (Schoendorf & Crowley, 1995;Schoendorf et al., 1996a). Using simulations with the National Center for Atmospheric Research (NCAR) Thermosphere/Ionosphere General Circulation Model (TIGCM), Schoendorf et al. (1996b) have reproduced such density holes, finding that "the high latitude density structure is located completely inside the circulation cells. High densities are inside the anticyclonic flow and low densities are inside the cyclonic flow." Under some conditions "there is a low density region inside the anticyclonic flow. The heating causes the densities outside the flow regions to remain high, while the neutral circulation depletes the density from within the flow" (Schoendorf et al., 1996b). Through use the Coupled Thermosphere-Ionosphere Model (CTIM) (Fuller-Rowell et al., 1987), Fuller-Rowell et al. (1999 found this explanation: "The neutral density holes are dynamically driven, created by the high velocity neutral winds from ion-drag. During strong geomagnetic activity momentum forcing resulting from ion-drag drives neutral winds of several hundreds of m/s. Inertial forces on the gas parcels cause a divergence in the vortex-like circulation. Back pressure induced by the divergent wind field does not completely balance the flow, and air rises within the vortex core, to maintain continuity, and cools adiabatically to produce a low temperature. Due to the cooling that takes place, in spite of the strong Joule heating that must accompany the ion-drag, the cold vortex core appears as a density hole" (Fuller-Rowell et al., 1999). Schlegel et al. (2005) reported on density maxima and minima in the higher altitude CHAMP data, finding elevated densities near the cusp, with minima found above 75° magnetic latitude between 4 and 7 magnetic local time. It is not clear if these minima persisted between multiple orbits. The previous publications reported that the density holes are found near dawn or dusk local times at high latitudes, while the results graphed here from the CHAMP data show a persistent minima near local midnight just outside of the auroral oval. It seems that a different explanation is required. On the other hand, these results do confirm the prior observations and simulations in that lower densities can be found during the 20 October event within the counter-clockwise flow in the dawn-side convection patterns. Examples are found in Figure 6 in the orbit paths indicated for 10:58 UT, 12:31 UT, and 14:04 UT, and at additional times in the Supplemental Information figures.
A global map of the time lags between the polar heating and the thermosphere's response was obtained from the database of CHAMP and GRACE density measurements, that spans several years. The most obvious feature in this map, shown in Figure 7, is that the shortest time lags are found in the polar regions, particularly near noon local time, while the time lags at lower latitudes are much longer. This behavior is not at all surprising. There is a sharp gradient in the lag times between 60 and 30° geographic latitude. While these results generally agree with the previous findings (Fuller-Rowell et al., 1994;Lu et al., 2014;Oliveira & Zesta, 2019;Oliveira et al., 2017), these maps show for the first time measured values of the lag times on a global scale. We did not make any attempt to sort these data according to the magnitude of the heating that preceded each density measurement during some arbitrary time period. Wang et al. (2020) had examined the response of the thermosphere during multiple geomagnetic storms of varying magnitude. They found that larger storms resulted in longer time lags between thermospheric density and Joule heating.
One particularly new and noteworthy result is the lags between the heating and the increasing density are shorter near noon than on the night, with Figure 7 map being in agreement with the preliminary results shown in Figure 2 through 5. The lags are particularly long (up to 350 min) at 30° latitude around 18 hr local time. This extra long delay could perhaps be due to circulation of colder neutrals from the nightside to dayside, following the sub-auroral ion flow that drags the neutrals.
Slightly different results are obtained when the time lags are calculated with respect to the derived ΔT value, as shown in Figure 9. As ΔT is obtained using an integration of the auroral heating over time (followed by exponential cooling when heating ceases), the peak in the ΔT value itself lags behind the peak in the heating by up to 3 hr (Figure 8). As a result, these time lags are smaller by a similar amount, and even negative in places, where the neutral density can reach the maximum prior to the ΔT value. The maps in Figure 9 have less variability than those in Figure 7, but also show a more pronounced minimum near the dayside polar cusps in both hemispheres. The time lags near the equator are considerably longer at 20 hr local time than at 4 hr local time. It is possible that the influx of solar radiation delays the cooling of the thermosphere on the dayside, causing it to lag farther behind the decreasing part of the ΔT curve.
In retrospect, an earlier paper by Licata, Mehta, Weimer, Drob, et al. (2022) suggests similar results. They had constructed three different models using Machine Learning (ML), which were used to study post-storm thermospheric cooling. The CHAMP neutral density measurements were used to build one ML model and also for comparison with all model results. The geomagnetic indices SYM-H (Iyemori, 1990) were used in the CHAMP-ML model to indicate the level of activity and ionospheric heating. Their Figure 2 showed that inputing the SYM-H indices without a time lag had the largest response at the dayside pole, followed by the dayside equator and nightside pole. The nightside equator had the strongest response when driven by a 3-hr lagged SYM-H, hinting that there is a 3-hr delay between the energy input and the response in density at this location. Comparable findings for geomagnetic storm and post-storm modeling were obtained in a related paper , in which a ML model was created to provide uncertainty estimates for the NRLMSIS 2.0 model.

Summary and Conclusions
The results presented here demonstrate how the thermosphere's neutral density responds to auroral energy changes with time lags that vary across the globe. In general, the neutral density changes first in the polar regions, as expected, since this is where the heating occurs. Near the equator, the dayside responds earlier than on the nightside. The time lags are longest near 18-20 hr local time.
The maps that were derived, showing the time delays as a function of location on the globe, can be incorporated into empirical models of the thermosphere in order to improve their accuracy. One example is the EXospheric TEmperatures on a PoLyhedrAl gRid (EXTEMPLAR) model (Weimer et al., 2020), that uses the same grid as the maps presented here. Developers of numerical models of the thermosphere, such as those mentioned by Shim et al. (2012) andKalafatoglu Eyiguler et al. (2019), may find it useful to compare their neutral density calculations with these findings.

Data Availability Statement
A companion data archive is available at https://doi.org/10.5281/zenodo.7667515. This archive contains the numerical values that are mapped in Figures 7 and 9, in the Hierarchical Data Format (HDF5). Included in this archive is a Python program that can read these data files and produce similar illustrations using the Python matplotlib package. The archive also contains the exospheric temperature values used to derive these results, along with the polar heating values and ΔT as a function of time for the years 2000 through 2019. The original CHAMP and GRACE neutral density measurements are available online at https://tinyurl.com/densitysets. IMF measurements from the Advanced Composition Explorer (ACE) satellite are available from the NASA archives at https://cdaweb.gsfc.nasa.gov/pub/data/ace. Refer to Emmert et al. (2020) for how to access the NRL MSIS 2.0 model. The solar indices are available from Space Environment Technologies at http://sol.spacenvironment.net/ JB2008/indices.