3-D global hybrid simulations of magnetospheric response to foreshock processes

It has been suggested that ion foreshock waves originating in the solar wind upstream of the quasi-parallel (Q-||) shock can impact the planetary magnetosphere leading to standing shear Alfvén waves, i.e., the field line resonances (FLRs). In this paper, we carry out simulations of interaction between the solar wind and terrestrial magnetosphere under radial interplanetary magnetic field conditions by using a 3-D global hybrid model, and show the properties of self-consistently generated field line resonances through direct mode conversion in magnetospheric response to the foreshock disturbances for the first time. The simulation results show that the foreshock disturbances from the Q-|| shock can excite magnetospheric ultralow-frequency waves, among which the toroidal Alfvén waves are examined. It is found that the foreshock wave spectrum covers a wide frequency range and matches the band of FLR harmonics after excluding the Doppler shift effects. The fundamental harmonic of field line resonances dominates and has the strongest wave power, and the higher the harmonic order, the weaker the corresponding wave power. The nodes and anti-nodes of the odd and even harmonics in the equatorial plane are also presented. In addition, as the local Alfvén speed increases earthward, the corresponding frequency of each harmonic increases. The field-aligned current in the cusp region indicative of the possibly observable aurora is found to be a result of magnetopause perturbation which is caused by the foreshock disturbances, and a global view substantiating this scenario is given. Finally, it is found that when the solar wind Mach number decreases, the strength of both field line resonance and field-aligned current decreases accordingly.


Introduction
Planetary magnetic fields provide effective obstacles to the solar wind which exists ubiquitously in astrophysical environments (Chapman and Ferraro 1930). Since solar wind flows are supersonic, collisionless planetary bow shocks emerge to occupy an exclusive position in front of each magnetized planet (Russell 1985). The study of interactions between planetary bow shocks and magnetic fields, therefore, lends important insight not only into the associated phenomena of collisionless shocks (Lu et al. 2009), but also into the nature of the magnetized planetary obstacles in response to the solar wind impact (Blanco et al. 2006). Among those, the Earth's bow shock as well as its geomagnetic field response is of particular interest, as so much of the planetary astrophysics was learnt from it.
It is well known that the Earth's bow shock not only primarily mediates the supermagnetosonic solar wind flow forming the magnetosheath, but also reflects a significant portion of energetic particles forming suprathermal backstreaming particles in the quasi-parallel (Q-||) shock region where the angle between the shock normal and interplanetary magnetic field (IMF) θ B n 45 • (Asbridge et al. 1968;Tsurutani and Rodriguez 1981;Eastwood et al. 2005). Those reflected suprathermal ions in turn interact with the incoming supersonic solar wind plasma, resulting in kinetic instabilities, generating ultralow-frequency (ULF) waves, and forming the Earth's foreshock (Omidi and Winske 1990). The foreshock disturbance phenomena, such as hot flow anomalies (Schwartz et al. 1985), foreshock cavities (Thomas and Brecht 1988) and foreshock bubbles (Omidi et al. 2010), can create significant magnetospheric impacts through the magnetosheathmagnetopause system and generate magnetospheric ULF waves (Archer et al. 2015), substantiated by recent observations from spacecrafts and ground stations (Eastwood et al. 2011;Hartinger et al. 2013;Zhao et al. 2017). Note that it is termed as foreshock disturbances in this paper for general foreshock processes and waves.
The ubiquitously pre-existing magnetospheric ULF waves are important to the planetary magnetosphere (Chian and Oliveira 1994;Zong et al. 2017), and can help us better understand the global processes between the solar wind and the magnetosphere (Mcpherron 2005;Claudepierre et al. 2008). These magnetospheric ULF waves have been suggested to be driven by a variety of mechanisms through either internal or external processes (Claudepierre et al. 2010). For the former, it has been shown that magnetospheric plasma with unstable configurations can lead to a local generation of ULF waves, such as the drift-mirror instability due to the pressure anisotropy in the Earth's inner magnetosphere (Hasegawa 1969). Another example regarding the internal processes is the generation of Alfvén waves by an inverted plasma energy distribution in the magnetosphere (Hughes et al. 1978). On the other hand, examples of externally driven magnetospheric ULF waves are related to Kelvin-Helmholtz instability caused by shear flows at the magnetopause (Southwood 1968;Eriksson et al. 2006), fluctuations in the solar wind dynamic pressure (Kepko et al. 2002), as well as foreshock processes (Hartinger et al. 2013;Russell et al. 1983;Takahashi et al. 1984Takahashi et al. , 2016Menk 2011;Bier et al. 2014;). In this paper, we focus on one of the external mechanisms-ion foreshock processes.
One of the most frequently observed planetary magnetospheric ULF waves is the standing Alfvén waves, also known as toroidal Alfvén waves (with azimuthal magnetic and velocity perturbations) or field line resonances (FLRs) (Takahashi and Ukhorskiy 2007;James et al. 2019;Rusaitis et al. 2021). These toroidal Alfvén waves play important roles in the magnetospheric physics. For instance, it has been suggested that electrons can be adiabatically accelerated by a drift-resonance through the interaction with toroidal mode ULF waves (Elkington et al. 1999;Hudson et al. 2000), leading to radial transport of energetic electrons in the outer radiation belt (Ukhorskiy et al. 2005). The FLRs are resonant Alfvénic oscillations of closed geomagnetic field lines whose foot points reside in the ionosphere. Moreover, the FLRs are of signatures of discrete frequencies, and such discrete frequencies arise from the restriction to an integral number of half wavelengths along the closed field lines (Dungey 1955;Cummings et al. 1969). This may be schematically illustrated by the classical figure of harmonic field line oscillations [e.g., Fig. 2 in Southwood and Hughes (1983)].
Theoretical endeavors have been made to describe the generation of the FLRs, which can be so far ascribed to four proposed mechanisms: resonant cavity, waveguide, surface waves and mode conversion. In the resonant cavity model (also called "box" model), the ionosphere and the magnetopause can serve as reflecting boundaries to reflect signals propagating in the magnetosphere either across (e.g., compressional waves) or along (e.g., Alfvén waves) the magnetic field lines, so that the compressional wave frequencies are quantized the same as the shear Alfvén wave frequencies (Kivelson and Southwood 1986;. It has also been proposed in the "box" model and observed by spacecrafts that poloidal Alfvén waves can evolve into toroidal Alfvén waves through phase-mixing process Sarris et al. 2009a). The waveguide model is a revised model for cavity modes, with anti-sunward propagation and reflection of discrete frequency waves at the magnetopause and at the turning points on dipolar field lines (Samson et al. 1992;Wright 1994). The Kelvin-Helmholtz instability model assumes the shear flows on the magnetopause producing surface waves whose amplitude decays from the surface, while the compressional wave modes couple the boundary to the resonant magnetic field lines (Southwood 1974). By contrast, in the mode conversion model, the incoming compressional waves are mode-converted to shear Alfvén waves at the locations where Alfvén resonance condition is satisfied (Chen and Hasegawa 1974a), but only the standing shear Alfvén waves along the closed field lines can be finally excited inside the magnetosphere due to the reflection boundary conditions. Therefore, the mode conversion is regarded as a direct driving mechanism, which will be focused in this paper.
Recent satellite and ground-based observations have shown that the transient and localized disturbances in the foreshock, a previous neglected source providing the compressional energy, are also able to generate standing Alfvén waves in the magnetosphere Wang et al. 2018Wang et al. , 2019Wang et al. , 2020. In addition to the global effects on the magnetopause, foreshock processes may also play an important role both in the longitudinal localization of magnetospheric ULF waves and in the plasma precipitation into cusp region. Obtaining a global view of FLRs as magnetospheric response to foreshock processes from spacecraft observations is challenging, since it requires coordinated multi-point observations in the foreshock and the magnetosphere, which have been only achieved hitherto in quite a few case studies Wang et al. 2020). On the other hand, although 2-D kinetic models can be used to study the entry of compressional perturbations and their propagation into the outer magnetosphere (Takahashi et al. 2020),fully 3-D global kinetic models are needed to study the magnetospheric response as a result of ion foreshock processes, as the bow shock and the magnetosheath continuously expand sunward in the 2-D approximation (see comments made in Lin and Wang (2005)) and the mode conversion at the magnetopause is of 3-D physics (Shi et al. 2013). So far, no self-consistent global simulations have been done to demonstrate such defined processes. Therefore, to fill this gap, we carry out 3-D dayside global hybrid simulations to study the FLRs in response to foreshock disturbances.
In this paper, we use the radial IMF and solar wind conditions similar to THEMIS observations by Wang et al. (2018). In addition, since the high solar wind Mach number conditions ( M A > 3 ) are prevalent at the Earth's bow shock (Wilkinson 2003), we present two supercritical foreshock cases associated with M A = 8 and M A = 5 to qualitatively describe the solar wind Mach number dependence. The remainder of this paper is organized as follows. Section "Model" briefly describes the hybrid model that has been used in this work. Section "Results and discussion" demonstrates the simulation results of the toroidal mode of the FLRs generated by the foreshock disturbances, and the field-aligned current (FAC) to the ionosphere produced by the FLRs. In the end, Sect. "Summary" concludes this paper.

Model
The dayside global-scale hybrid model employed in this paper was developed by Swift (1996) and then implemented in 3-D by Lin and Wang (2005). More detailed numerical discussion of this hybrid model can also be found in Shi et al. (2017). Note that in this model, ions are treated as fully kinetic particles, while electrons are regarded as a massless fluid.
The motion of ions is described by their equation of motion (in the simulation units): where v i is the ion particle velocity, E is the electric filed, B is the magnetic field, ν is the collision frequency, and V i and V e are the bulk flow velocities of ions and electrons, respectively.
The electric field is obtained from the momentum equation of the electron fluid: where P e is the thermal pressure of the electrons, and N is the electron number density (equal to ion number density based on the assumption of quasi-charge neutrality), respectively.
Ampère's law is used to evaluate the electron bulk flow velocity: where α = 4πe 2 /m i c 2 , e is the electron charge, and m i is the ion mass. Note that the constant α is related to the ion inertial length d i = 1/ √ αN 0. Finally, the magnetic field is updated with Faraday's law: For the results presented in this paper, the physical quantities are normalized as follows. The magnetic field B is scaled by the IMF B 0 , ion number density N by the solar wind density N 0 , time t by the inverse of the solar wind ion gyrofrequency −1 i0 , plasma flow velocity by the solar wind Alfvén speed V A0 , temperature by m i V 2 A0 , and length by the Earth radius R E .
Spherical coordinates are utilized in our 3-D simulation, the same as in Lin and Wang (2005) and Shi et al. (2013). The simulation domain contains the global system of the bow shock, magnetosheath, and magnetosphere in the dayside region with GSM x > 0 and a geocentric distance 3 R E ≤ r ≤ 33 R E . The Earth is located at the origin x, y, z = (0, 0, 0) as shown in Fig. 1a. The inner boundary at r = 3 R E is assumed to be perfectly conducting, and the outer boundary at r = 33 R E is applied with the solar wind inflow boundary conditions, while all other boundaries are utilized with outflow boundary conditions with the constant IMF. In addition to the fully kinetic ions which dominate in the outer magnetosphere, a cold and incompressible ion fluid is assumed to be dominant and coexist with particle ions in the inner magnetosphere, mainly occupying r < 6.5 R E (Swift 1996;Lin and Wang 2005), and its number density is given by N f = 1000/r 3 [1 − tanh(r − 6.5)] normalized to N 0 (Lin et al. 2014). The total ion bulk flow speed is determined by the subscripts p and f stand for the discrete particles and the fluid, respectively. The fluid ion velocity is updated by dV f /dt = E + V f × B (Swift 1996). The constant N f condition corresponds to the approximation �N /N ≪ 1 , where N is the density perturbation in the inner magnetosphere. Note that in this paper the collision frequency ν (equivalent to the resistivity) is set to be zero (Shi et al. 2013). Initially the uniform solar wind and IMF are assumed to occupy the region of r > 15 R E , and a centered 3-D dipole field combined with an image dipole field occupies the geomagnetic region within r < 15 R E . A transition layer (with a half-width of 0.5 R E ) exists between the two regions centered at r = 15 R E . A uniform solar wind plasma flows into the system along the −x direction with an isotropic drifting Maxwellian distribution and with ion beta value β i = 0.5 , and a total number of ∼ 5 × 10 8 particles are utilized in the simulations. The time step interval to advance the particles is 0.05 −1 i0 . A total grid of N r × N θ × N φ = 250 × 140 × 140 is used. Nonuniform grids are used in the radial direction, with a higher resolution of r ∼ 0.05R E through the magnetopause, magnetosheath and bow shock, while r ∼ 0.1 − 0.25R E in the upstream solar wind regions. Under radial IMF conditions, the ion foreshock disturbances are generated and propagate mostly in the radial direction, and the foreshock is located in the subsolar region. Note that such a resolution is enough to resolve the ion foreshock disturbances (spatial size up to several R E , ref. Fig. 1) in our hybrid simulations (Lin and Wang 2005).
We have performed runs with doubling spatial resolution and macro-particle number density. The results are found to be similar to those shown in this paper, indicating a good convergence of simulations. The hybrid model is valid for low-frequency physics with ω ∼ � i as well as where ω is the wave frequency, k the wave number, i the ion gyrofrequency, and ρ i the ion Larmor radius. For this range of wave frequency and wavelength, the ion kinetic physics in the near-Earth instabilities are resolved with grid sizes ∼ ρ i or ion inertial length d i . The finite ion gyroradius effects are resolved with particle time steps t much smaller than the gyroperiod. It is worth noting that to save the computational resources while retaining ion kinetic effects on a global scale, i.e., while being able to resolve the particle scales of ρ i and i , the ion inertial length of the solar wind in our simulations is chosen to be d i0 = 0.1 R E , a few times larger than that in reality. This choice can be justified due to (but not limited to) a number of facts as follows. (I) The wavelengths of the ion foreshock disturbances in our simulations can be longer than those in reality since the physics is scaled with ρ i , but they still satisfy the relation d i ≪ R E . (II) The Q-|| shock width is determined by such wavelengths. (III) The ion foreshock extends to large scale regions in length ( ≫ d i ). (IV) Our hybrid model and simulations have been proved to be able to self-consistently capture the ion kinetic physics in the global system of solar wind-magnetosphere interactions (Lin and Wang 2005). For instance, Wang et al. (2009) have shown that the characteristics of wave power intensity at the Q-|| shock are consistent with spacecraft observations, and Lee et al. (2020) have found a good agreement by comparing the magnetosheath wave powers under the same solar wind conditions between the magnetospheric multiscale (MMS) observations and our hybrid simulations.

Results and discussion
Global view and identification of foreshock structures for Case I Figure 1a presents the global view of ion density N at t = 100 −1 i0 in the noon-meridian and equatorial planes showing the 3-D structures of self-consistently generated global system. a Presents the global view of N in the noon-meridian and equatorial planes, and the Earth is also shown at the origin in our GSM coordinate system. The white dashed lines roughly mark the magnetopause and bow shock positions, respectively. b, c Give the zoom-in equatorial views of Q-|| foreshock. Typical magnetic field lines with arrows showing the field line direction are displayed in a and b. Note that the 3-D field lines in a are partly hidden by the contour planes which are completely opaque. The wiggled field lines in the foreshock region are caused by the ion foreshock disturbance. Typical foreshock cavities and compressional pressure pulses are illustrated in c obtained from the simulation together with typical magnetic field lines (arrows showing the field line direction), which shows the 3-D self-consistently generated dayside global system with ion foreshock, bow shock, magnetosheath, magnetopause and magnetosphere. Figure 1b, c zooms in the equatorial B and N contours of the ion foreshock, which can be identified in the upstream of the Q-|| shock by the perturbed/wiggled magnetic field lines as shown in b. The ion foreshock disturbances, including compressional pressure pulses (Shi et al. 2013) and diamagnetic cavities (Lin and Wang 2005), as respectively, marked in c, are generated mainly by the interaction of backstreaming ions and incoming solar wind (Lin and Wang 2005), and other wave-particle interactions and nonlinear evolution of waves which lead to the generation of structures observed in the foreshock (Eastwood et al. 2005;Wilson 2016). For instance, Wang and Lin (2006) found that in the nonlinear stage of the beamplasma interaction, the nonresonant modes decay and the resonant modes grow through a nonlinear wave coupling, and the interaction among the resonant modes leads to the formation of filamentary structures, which are the k⊥B structures of magnetic field, density, and temperature. The pressure pulses of the foreshock disturbance are carried by the super-Alfvénic solar wind into the downstream of the Q-|| (Lin and Wang 2005), through the magnetosheath, and onto the magnetopause (Shi et al. 2013). The bow shock can be identified from the deflection of the IMF, which is originally along the x -direction, and from the enhanced magnetic field and ion density on the earthward side. The magnetopause can be distinguished through the sharply increased (decreased) magnetic field (ion density) toward the Earth. Specifically, the self-consistently generated bow shock is around a standoff radial distance of ∼ 13 R E , and the magnetopause is roughly at ∼ 11 R E in the subsolar region, while the magnetosheath is located in the region between the bow shock and the magnetopause. The average locations of the bow shock and the magnetosheath are marked by the white dashed lines in Fig. 1a, while their actual shapes oscillate in response to the foreshock perturbations. In the following, we provide a detailed tracking of these perturbations and analyze their direct impacts to the magnetosphere. Shi et al. (2013) have shown that the compressional ion foreshock disturbance pulses can generate surface waves, such as kinetic Alfvén waves, at the magnetopause. In the interaction between the foreshock disturbance pulses and the magnetopause, the magnetopause is largely perturbed by the thermal and dynamic pressures due to its pressure balance nature. To see how the pressure pulses generated from the ion foreshock region penetrate through the magnetosheath and perturb the magnetopause, Fig. 2 shows the evolution of the pressure components (i.e., magnetic pressure P B , thermal pressure P th , total pressure P tot , dynamic pressure P dyn and dynamic pressure components) as well as the flow velocity components along the Sun-Earth line at four different moments of time. Note that for simplicity we only consider the total pressure as P tot = P B + P th + P dyn , and P dyn = P dyn,x + P dyn,y + P dyn,z . At each moment of time, the black arrows trace the same series of dynamic pressure pulses originating from the ion foreshock, propagating through the magnetosheath, and finally interacting with the magnetopause. At t = 102 −1 i0 , a series of dynamic pressure pulses in the foreshock region emerge and marked by the black arrow in Fig. 2a, f. At this moment of time, the P tot of the traced pressure pulses is almost dominated by the dynamic pressure, which is mostly determined by the x-component of the dynamic pressure ( P dyn,x ). In other words, at this moment of time, the P tot of traced pulses is predominated by the earthward plasma flows, since P dyn,x is largely determined by the earthward flow speed ( V ix < 0 ) as shown in Fig. 2k. It can be seen that in the foreshock and its immediate downstream, the V ix dominates in the flow velocity components so that the fluctuation in V ix (e.g., around the traced arrow, with a tailward enhancement) results in the dynamic pressure pulses. At t = 105 −1 i0 , the same dynamic pressure pulses move to the downstream of the shock as shown by the black arrows in Fig. 2b, g. The dynamic pressure of the traced pulses increases to its maximum due to the density compression though the shock, shown in Fig. 2l. At t = 112 −1 i0 , the traced dynamic pressure pulses move into the middle of magnetosheath as shown by the black arrows in Fig. 2c, h. It can be seen that, for the traced pulses, although P dyn,x still dominates P dyn , the magnitude of P dyn decreases as the pulses slow down towards the Earth along the Sun-Earth line, and P th becomes comparable to the P dyn,x . At t = 116 −1 i0 , the traced pulses are approaching the magnetopause as shown in Fig. 2d, i. Although P dyn further decreases to be comparable to P B , it can be seen that the pressure pulses are still dominated by P dyn . In the meantime, the traced pulses are greatly expanded in the radial direction, and the expanded leading part at x = 11 R E begins to interact with the magnetopause, resulting in deflection of the magnetosheath flow (i.e., the increase of v z and P dyn,z ). At t = 119 −1 i0 , the traced pulses interact with the magnetopause, as shown by the shaded area in Fig. 2e, j. The P dyn becomes almost trivial on the magnetosheath side, while it increases on the magnetopause side as marked by the black arrow at this moment of time. The dramatic decrease in P dyn,x and increase in P dyn,y shown in (j) indicate that the pulsed magnetosheath flow is deflected by the magnetopause. At the same time, the locally increased P dyn,y contribute to the P tot as shown in (e). Note that the foreshock disturbance pulses are continuously generated and impinge onto the magnetopause. When stronger dynamic pressure pulses approach and interact with the magnetopause, the magnetopause is compressed, as shown at t = 116 −1 i0 and at t = 119 −1 i0 . On the other hand, when weaker dynamic pressure pulses interact with the magnetopause, the a-e Evolution of magnetic pressure ( P B ), thermal pressure ( P th ), dynamic pressure ( P dyn ) and total pressure ( P tot = P B + P th + P dyn ). f-j Evolution of dynamic pressure components ( P dyn,x , P dyn,y , P dyn,z ), as well as dynamic pressure ( P dyn = P dyn,x + P dyn,y + P dyn,z ) shown in the corresponding left column. k-o Evolution of flow velocity components ( V x , V y , V z ). The black arrows trace a series of dynamic pressure pulses originating from the ion foreshock, propagating through the magnetosheath, and finally interacting with the magnetopause. The shaded areas in e, j and o mark the location inside and near the magnetopause, where the deflection of magnetosheath flow occurs. The vertical black dashed lines reference the average positions of the magnetopause and the bow shock, respectively magnetopause is relaxed, i.e., the magnetopause moves toward the Sun, as shown at t = 112 −1 i0 . It can also be seen from Fig. 2a-e that the whole process of the traced pulses perturb the magnetopause by ∼ 0.1 − 0.2 R E . Note that, by definition, such traced high speed dynamic pressure pulses can also be referred to as the magnetosheath jets (Plaschke et al. 2018).
To give a closer look of the propagation of foreshock compressional waves through the magnetosheath, Fig. 3 shows the temporal evolution of spatial profiles of (a) B and (b) N along the Sun-Earth line during the time interval t = 80 − 140 −1 i0 . Again, the magnetopause boundary layers can be seen from the sharp increase in B and the abrupt decrease in N , while the bow shock position is roughly around 13 R E . Figure 3 therefore clearly shows the oscillation of the magnetopause due to the pressure pulses discussed in Fig. 2. The compressional waves can be perceived by the in-phase relation of B and N (Shi et al. 2013), and some of their propagation is marked in Fig. 3. It can be seen that the compressional waves propagate along the Sun-Earth line from the foreshock into the magnetopause boundary layer. For example, at t ∼ 100 −1 i0 , a new packet of compressional waves emerge from the foreshock, and propagate through the magnetosheath. At t ∼ 120 −1 i0 , the compressional pulses arrive at the magnetopause boundary layer, interact with the magnetopause, and transmit into the magnetosphere.
Previous simulations (Lin and Wang 2005;Shi et al. 2013) have shown that when the compressional wave packets interact with the magnetopause, kinetic Alfvén waves with k ⊥ ≫ k || are excited near and inside the magnetopause boundary layer where the Alfvén resonance condition is satisfied (Chaston et al. 2005), consistent with the mode conversion process (Chen and Hasegawa 1974b;Johnson and Cheng 1997;Lin et al. 2012). According to the mode conversion model, as the broadband remainder of compressional waves penetrates into the magnetosphere, standing toroidal Alfvén waves can be excited along the closed geomagnetic field.
It has been shown that the relatively large amplitude FLRs are produced at the source frequency that is matched to a harmonic of the Alfvén shear modes (Lee and Lysak 1991). To check this property, we plot the power spectra of the compressional ( B ) and toroidal ( B φ ) wave components of the magnetic field along the Sun-Earth line from the foreshock into the magnetosphere in Fig. 4. It clearly shows that the self-consistently generated foreshock and toroidal Alfvén wave spectra cover a wide range including ∼ � i0 /2 . The frequency ω normalized by the solar wind i0 at various locations is obtained from the simulation data t = 80 − 200 −1 i0 . During this time interval, the bow shock is oscillating around X ∼ 13 R E , and the overall magnetopause position is around X ∼ 10.5 R E . Around the oscillating magnetopause boundary layer from X = 9.5 − 10.5 R E , locally increasing compressional wave powers appear with a spectral width of ω ∼ 0.2−0.4 because most of the magnetic pressure pulses impinge onto the magnetopause and this enhances local magnetopause perturbations. Meanwhile, some of the compressional wave power is transmitted into the magnetosphere ( X < 9 R E ).
Note that due to the Doppler effects, the frequency of the same wave mode appears to be different in the solar wind and various locations throughout the magnetosheath with the variable background flow velocities (Shi et al. 2013). The Doppler shift of the frequency is �ω = k x V sw ∼ 0.75 � −1 i0 , where k x and V sw are wave number and solar wind speed in the foreshock, respectively. The wave frequency in the solar wind plasma frames is estimated to be where ω is the frequency in the simulation (Earth) frame. The negative sign of the wave frequency ω r corresponds On the other hand, although toroidal wave power is also present in the magnetosheath, the toroidal waves do not propagate into the magnetosphere. This can be seen by the fact that there is a disruption at ∼ X = 9.5 R E , i.e., the toroidal waves are excited locally in the magnetosphere. Also note that in our simulations we do not observe cavity mode waves. Figure 5 shows the time sequence of toroidal mode components B φ and E r along a closed dipole field line through the equatorial radius of L = 8 R E in the noonmeridian plane from north to south; i.e., the starting point is located at the ionospheric boundary in the north, and the ending point is in the south boundary. To emphasize the field line oscillations near the equator, B φ has been multiplied by an envelope factor f (S) = 0.5[tanh( S − 0.2S 0 ) − tanh(S − 0.8S 0 )] , where S and S 0 are the distance and field line length along the dipole field line from the north to south with the starting and ending points, respectively. The vertical axes of each plot represent the time series of the corresponding toroidal modes. Note that the selected field line is approximately the same line at different moments of time, since the magnetospheric field undergoes a systematic low-frequency oscillation in response to the foreshock oscillations. It can be seen that the oscillation of the field lines in B φ and E r bounces between  Case I: toroidal mode components E r and B φ at L = 8 R E for Case I in the noon-meridian plane as a function of S , where S is the distance along the field line from north to south. Note that the magnitudes of E r and B φ are marked on the bottom right north and south. The perturbation corresponds to the combined effects of multiple modes in k || between the boundaries. It can be seen that before t ∼ 20 −1 i0 , the fluctuation along the field line is negligible, but becomes more and more stronger later on. After t ∼ 40 −1 i0 , a nearly standing wave pattern can be perceived, and the wave mode number approximately corresponds to a half wavelength; that is, the fundamental mode is dominant. In addition, the resonance frequency satisfying ω = k || V A around the equator at L = 8 R E (ref. Fig. 9 below) is found to be ∼ 0.2 −1 i0 , consistent with the resonance period of ∼ 35 −1 i0 (ref. Fig. 7 at ∼ 85 − 120 −1 i0 below). Also note that the magnitude of the toroidal wave components E r is ∼ 1 mV/m, which is on the same order of magnitude compared with those in recent observations at similar L shells in the dayside magnetosphere Wang et al. 2019;Sarris et al. 2009b). For instance, in Wang et al. (2019) the amplitude of the toroidal mode component E r at L ∼ 9 R E is ∼ 2 mV/m (ref . Fig 3e), and in Sarris et al. (2009b) the E r observed by THEMIS-D at L ∼ 8 R E is about 2-3 mV/m (ref. Figs. 1, 2 THEMIS-D around 7:00UT). A close comparison, however, requires similar solar wind and IMF conditions between the simulations and observations. In the present study, the solar wind conditions are completely steady and the IMF orientations are radial, while in the observations the solar wind conditions are changing with time and the IMFs are more oblique. Since the amplitude of the FLRs arising from external drivers depends on the strength of the initial disturbance, the dependence of the FLRs on the foreshock disturbances needs to be further investigated under more general solar wind conditions. Figure 6a shows the dominant fundamental harmonic of the eigenstructures of the toroidal component E r between t = 0 − 200 −1 i0 in time sequence. The oscillations near the equator become more and more evident at t > 40 � −1 i0 , and the dominant wave period of the fundamental harmonic appears to be ∼ 30 −1 i0 . Figure 6b shows the root integrated powers (RIPs) of the fundamental, second-and third-order harmonic eigenmodes of E r along the dipole-like field line at L = 8 R E . In this paper, the RIP for a given time series is defined the same way as in Claudepierre et al. (2010), and thus the RIP is a measure of the fluctuation amplitude at the frequency of interest. The fundamental harmonic eigenstructure reaches its maximum near the equator, but not exactly at the equator due to that a much longer simulation time may be needed for a more precise solution at the low frequency. It also can be seen that the second harmonic eigenstructure has two peaks off the equator, and the third harmonic eigenstructure has three peaks along the same magnetic field line. Such results are qualitatively consistent with MHD simulations under prescribed external driving perturbations (Claudepierre et al. 2010), as well as the FLR theories (Cummings et al. 1969). In addition, the fundamental harmonic has the largest wave amplitude along the magnetic field line, consistent with the direct observation in Fig. 5 that the fundamental mode of the FLRs is dominant. For higher order of harmonics, the power is much weaker, and the RIP calculation that includes the effects of nonuniform magnetic field along field lines in the eigenfunction analysis does not produce the results as predicted by the theories. This may be because there are other transient mode perturbations from the foreshock and the system is never in a pure steady state. In the following, to focus on the toroidal Alfvén mode for both the fundamental and the extended higher harmonics, we use a simplified wave power analysis. Specifically, we decompose the toroidal wave modes at the harmonic frequencies by using Fourier analysis into simple plane waves which satisfy the boundary conditions of the standing waves. Figure 7 shows the first 7 harmonics of toroidal mode oscillations (top) and spectra (bottom) of E r along the dipole field line in the noon-meridian plane at L = 8 R E . The oscillations in time sequence have been rescaled to show the oscillations along field lines for each harmonic mode, while their strength can be perceived from wave power spectra in the corresponding bottom panel. The rescaling means the amplitudes of higher harmonics are magnified in order to show the harmonic oscillations, in time sequence. b The RIPs of the fundamental, second and third harmonic mode structures of E r along the dipole-like field line in the noon-meridian plane at L = 8 R E . The vertical dashed line marks the equator position. The relative magnitude of oscillations for the fundamental harmonic is also marked at the right bottom in a as higher harmonics have much smaller amplitudes. It shows in the top panel that the sinusoidal bouncing oscillations become more and more evident at t > 60 � −1 i0 , and the higher the harmonic order, the smaller the period. Specifically, the period of the fundamental mode can be seen as ∼ 30 −1 i0 (e.g., t ∼ 90 − 120 −1 i0 ), while the second harmonic has a period of t ∼ 15 −1 i0 (e.g., t ∼ 80 − 95 −1 i0 ), and so on. From the bottom panel, it is seen that the fundamental mode has the largest wave power but the lowest frequency, while as the order of harmonics goes higher, the wave power decreases but the frequency increases. Besides, the frequency of the fundamental harmonic mode is ∼ 10 mHz, which is consistent with the results ( ∼ 7 − 8 mHz) from spacecraft observations around 8 R E (Sarris et al. 2009b). It also clearly shows the symmetry of E r for the harmonic modes of standing Alfvén waves about the equator. Moreover, since the perturbations in E r are symmetric (antisymmetric) for odd (even) harmonics with respect to the equatorial plane, only odd number toroidal modes for E r of the FLRs can be excited at the equator, while the counterparts for B φ are opposite at the equator (Lee and Lysak 1989;Sarris et al. 2009b). In other words, for the harmonics of toroidal Alfvén mode component E r , the odd (even) harmonic numbers would encounter anti-nodes (nodes) at the magnetic equator; while for B φ the even (odd) harmonic numbers would encounter anti-nodes (nodes). These alternate pattern of wave power is also obvious from the wave power spectra of E r at the equator shown in Fig. 8. It can be seen that after t ∼ 10 −1 i0 , the wave powers of the odd harmonic modes begin to grow, and the growth rates of these modes appear to be more or less the same. After t ∼ 40 −1 i0 , the wave powers of Note that the oscillations in time sequence have been rescaled to show the bouncing of field lines for each harmonic mode, but the strength can be perceived from wave power spectra Fig. 8 Case I: wave power spectra of E r over time for the first 7 harmonics at the subsolar equator and L = 8 R E . The wave power appears not to grow or decay. It also clearly shows the wave power trend at the equator for E r , i.e., the odd harmonics are much stronger than the negligible even harmonics, and the higher the harmonic order, the weaker the corresponding wave power these growing modes almost reach a steady state. On the other hand, however, the wave powers of the even harmonic modes grow marginally, and remain at a negligible level. It can also be seen that the fundamental mode is almost dominant in the wave power of E r . In addition, only the odd harmonics are excited, and the wave power of these modes decreases as harmonic order increases. It is also found that the total wave power of all the harmonic modes reaches almost a steady state after the wave amplitudes remain almost constant.

Magnetospheric response as FLRs for Case I
This feature is clearly seen in Fig. 9, which depicts the radial profiles of power spectral density (PSD) as a function of L for E r in (a) and B φ in (b) at the equator in the noon-meridian. The white dashed lines mark the first 7 orders of n-harmonic modes of the FLRs, and are obtained from the dipole field line approximation of local Alfvén speed in (c) by averaging B and N during t = 50 − 250 −1 i0 in the simulation. The enhanced oddharmonic numbers in E r and even-harmonic numbers in B φ are the piece of evidence for standing toroidal Alfvén waves. As the local Alfvén speed increases, the frequency of excited toroidal mode of the FLRs increases accordingly. At the equator, the wave power of fundamental harmonic mode in E r is the strongest, and the wave power decreases as the excited odd harmonic number increases. The same trend is seen in B φ but with even harmonic numbers. In short, the lower the excited harmonic mode number, the stronger the wave power.

Inner magnetospheric signatures due to FLRs for Case I
Recent aurora observations (Wang et al. , 2019 using the ground-based all-sky imager suggested that the foreshock disturbances can locally perturb the magnetopause, leading to the FLRs in the magnetosphere and the auroral FACs. Figure 10 shows the global view of (a) J || and (b) B φ contours in the equatorial plane and at r = 4 R E with typical closed magnetic field lines near the magnetopause. Note that the azimuthal sizes of the magnetopause perturbation are about 1 − 2 R E . At the high latitude ( ∼ 70 • ), the azimuthal mode with mode number m ∼ 5 of the FACs that flow along closed geomagnetic field lines is dominant. In view of their direction to the polar cap in the northern hemisphere (downward near the dawn and upward near the dusk) and their connection to the region near the magnetopause (indicating that it is mainly due to the magnetopause perturbation), these dominant FACs with m ∼ 5 azimuthal mode are consistent with the Region-1 current. It should be pointed out that fine azimuthal structures, indicating that high azimuthal wavenumber ( m ∼ 140 ) mode waves coupled with the FLRs (Wang 2019), are not observed in our simulations. This may be because we simply use the conducting inner magnetospheric boundary conditions in our model. Since only the FACs of Region-1 pattern are observed, it is unclear how the azimuthal mode number depends on the model assumptions. Our future work will focus on such FAC distribution including the nightside by employing an improved inner boundary conditions coupling to ionosphere. Shi et al. (2013) have shown that the Alfvénic structures that are mode-converted from compressional waves at the magnetopause can propagate along the magnetic field lines into the cusp region (ref. Fig. 12in that paper). In the same scenario, the low m-number J || shown in Fig. 10 can be caused by the perturbations originating from the compressional waves in the magnetosphere; i.e., the mode-converted Alfvén waves propagate along the closed magnetic field lines, resulting in the FACs into the cusp region. As the field line inside the magnetosphere is closed, the toroidal mode of standing shear Alfvén waves (i.e., the FLR) appears. Note that the large azimuthal modes (e.g., m 60 ) at even higher latitudes observed by Wang et al. (2018) are the FACs caused by the perturbations from the night side of the open magnetic field lines, so they are not considered in this paper.
Aurora phenomena, including diffuse and discrete aurora patterns, have been suggested to be an intrinsic feature of ionospheric response to the magnetospheric FLRs (Kozlovsky and Kangas 2002;Tanaka et al. 2012), and the local brightening of both aurora soon after foreshock disturbance events has been observed by allsky imager cameras . It is also suggested that the discrete aurora is generated by electrons accelerated along magnetic field lines, and is associated with the FACs which correspond to energetic particle precipitations . Therefore, to explore these possibly observable auroral features produced by the FLRs, we traced the FAC signatures near the inner simulation boundary. Figure 11 shows a polar contour view of evolution of the FAC power log J 2 || in northern hemisphere at r = 4 R E from t = 160 −1 i0 to t = 190 −1 i0 . The Region-1 type current shown in Fig. 10 are marked by black circles. The evolution of the FAC signatures along the closed field lines corresponding to L ∼ 6 − 8 R E is tracked with purple circles with Fig. 11 Case I: polar contour view of evolution of the FAC power log J 2 || in northern hemisphere at r = 4 R E from t = 160 −1 i0 to t = 190 −1 i0 . The azimuthally tailward moving FAC signatures along the closed field lines corresponding to L ∼ 6 − 8 R E is tracked with purple circles with their peaks marked by cross-dotted targets. Note that the Region-1 type current shown in Fig. 10 are marked by black circles Fig. 12 Case II global view of N at t = 100 −1 i0 in the noon-meridian and equatorial planes, showing the 3-D structures of the global system. The white dashed lines roughly mark the magnetopause and bow shock positions, respectively, and typical magnetic field lines are also displayed. The foreshock disturbances and the magnetosheath density are much weaker compared with Case I their peaks marked by cross-dotted targets. It can be seen that the purple circled FACs move duskward (i.e., azimuthally tailward), consistent with the convection of the driver waves near the magnetopause. The moving speed is estimated to be ∼ 400 m/s when mapped to the ionosphere (at ∼ 1.07 R E ), consistent with poleward moving speed of the discrete aurora in observations (Kozlovsky and Kangas 2002). Figure 12 shows the global view of ion density N at t = 100 −1 i0 in the noon-meridian and equatorial planes with typical magnetic field lines. The magnetopause and the bow shock positions are marked by the white dashed lines. The standoff radial distances of the bow shock and the magnetopause move farther away from the Earth, at 15 R E and 12.5 R E , respectively. Since the lower Mach number results in a smaller fraction of solar wind particles being reflected at the Q-|| shock, the foreshock disturbances are weaker compared with Case I. This in turn results in a lower wave power in the foreshock and less compressive fluctuations (Gary 1993;Turc et al. 2015Turc et al. , 2018. Figure 13 gives the toroidal Alfvén component E r of the fundamental harmonic mode along the closed field line in the noon-meridian plane at L = 8 R E for Case II. The oscillations of the fundamental eigenstructures are shown in the left panel. The period of the fundamental mode in Case II is roughly 50 −1 i0 (e.g., t ∼ 100 − 150 −1 i0 ), smaller than Case I ( ∼ 30 −1 i0 ). From the right panel, it can be seen that the maximum wave power of the fundamental harmonic by using plane wave analysis is ∼ 8 × 10 −3 (V/m) 2 /Hz (which is similar to the equatorial power based on the RIP method, and much smaller than the corresponding wave power ∼ 3 × 10 −2 (V/m) 2 /Hz at the same location in Case I).

Magnetospheric response for Case II
In short, it clearly shows that as the solar wind Mach number becomes lower, the period of the FLR becomes longer and the strength becomes weaker. Figure 14 shows the radial profiles of toroidal Alfvén wave power spectra at the equator for Case II. It can be seen that the wave power decreases dramatically at the same location, while other basic features of the FLRs remain quite similar to those for Case I. Specifically, in addition to the weaker fundamental harmonic wave power, higher harmonic wave powers also become weaker when the solar wind Mach number decreases. For instance, the odd-number harmonic wave powers for E r are almost one order of magnitude smaller than those in Case I. In addition, although similar to Case I, each of the harmonic wave power decreases as L decreases (towards the Earth side), the decay rate in Case II appears to be faster than the corresponding part in Case I. Figure 15 shows the global view of (a) J || and (b) B φ contours in the equatorial plane and at r = 4 R E with typical closed magnetic field lines near the magnetopause. It can be seen that the FAC signatures to the cusp region are similar to Case I. At the high latitudes, the azimuthal mode with mode number m ∼ 5 of the Region-1 type FACs (with m ∼ 5 azimuthal mode) that flow along closed geomagnetic field lines are still present and dominant, while the high m -number azimuthal mode is not observed. Although the correlated J || and B φ can be seen, the magnitude in both J || and B φ become smaller when the Mach number decr eases.
Note that the presented simulations only describe steady radial IMF cases, when the foreshock lies upstream of the subsolar bow shock. Under otherwise the same solar wind conditions, non-radial IMF configurations may result in less significant dayside foreshock disturbances, and thus weaker FLRs. Our future work will also focus on the FLRs and the FACs under more general IMF and solar wind conditions.

Summary
In this paper, self-consistent 3-D global hybrid simulations are carried out to investigate the FLRs which are generated by the foreshock disturbance processes under radial IMF conditions, and the simulation results support the scenario from the observations Wang et al. 2019) that the compressional waves arising from the Q-|| foreshock region excite toroidal Alfvén waves in the magnetosphere. These results provide insights for space environment of magnetized planets, and are summarized as follows: 1. This paper presents the first global simulations of the FLRs driven by foreshock processes. 2. The ion foreshock disturbances are generated self-consistently at the Q-|| shock by the interaction between the backstreaming superthermal and incoming supersonic ions, and it is shown that these foreshock disturbances propagate into the magnetosphere and externally drive the magnetospheric toroidal Alfvén waves, i.e., the FLRs, through the direct mode conversion process. 3. The spectrum of the self-consistently generated foreshock disturbances covers a wide frequency range including � i0 /2 . After considering the Doppler shift effects, it is found that the frequency band of the FLR harmonics match the source frequency range at the foreshock. 4. Although various ULF waves can be generated in the magnetosphere by the foreshock disturbances, the fundamental harmonic mode of the FLRs can still be distinguished. The analysis of harmonics modes of the FLRs shows that the fundamental mode is the strongest, and as the order of the harmonic increases, the corresponding wave power decreases. In addition, the wave power of the toroidal Alfvén waves seems to be not growing or decaying over simulation period of time. Since the solar wind inputs in the simulations are steady, it is suggested that the FLRs are part of the steady-state magnetosphere under quasi-radial IMF conditions. 5. At the equator, the alternating feature of toroidal Alfvén wave components E r and B φ is presented. As the Alfvén speed increases inwards the magnetosphere, the frequencies of toroidal Alfvén harmonic modes increase accordingly.