Numerical study of the generation and propagation of ultralow-frequency waves by artificial ionospheric F region modulation at different latitudes

Powerful high-frequency (HF) radio waves can be used to efficiently modify the upper-ionospheric plasmas of the F region. The pressure gradient induced by modulated electron heating at ultralow-frequency (ULF) drives a local oscillating diamagnetic ring current source perpendicular to the ambient magnetic field, which can act as an antenna radiating ULF waves. In this paper, utilizing the HF heating model and the model of ULF wave generation and propagation, we investigate the effects of both the background ionospheric profiles at different latitudes in the daytime and nighttime ionosphere and the modulation frequency on the process of the HF modulated heating and the subsequent generation and propagation of artificial ULF waves. Firstly, based on a relation among the radiation efficiency of the ring current source, the size of the spatial distribution of the modulated electron temperature and the wavelength of ULF waves, we discuss the possibility of the effects of the background ionospheric parameters and the modulation frequency. Then the numerical simulations with both models are performed to demonstrate the prediction. Six different background parameters are used in the simulation, and they are from the International Reference Ionosphere (IRI-2012) model and the neutral atmosphere model (NRLMSISE-00), including the High Frequency Active Auroral Research Program (HAARP; 62.39 N, 145.15W), Wuhan (30.52 N, 114.32 E) and Jicamarca (11.95 S, 76.87W) at 02:00 and 14:00 LT. A modulation frequency sweep is also used in the simulation. Finally, by analyzing the numerical results, we come to the following conclusions: in the nighttime ionosphere, the size of the spatial distribution of the modulated electron temperature and the ground magnitude of the magnetic field of ULF wave are larger, while the propagation loss due to Joule heating is smaller compared to the daytime ionosphere; the amplitude of the electron temperature oscillation decreases with latitude in the daytime ionosphere, while it increases with latitude in the nighttime ionosphere; both the electron temperature oscillation amplitude and the ground ULF wave magnitude decreases as the modulation frequency increases; when the electron temperature oscillation is fixed as input, the radiation efficiency of the ring current source is higher in the nighttime ionosphere than in the daytime ionosphere.


Introduction
The ultralow-frequency (ULF) waves with a frequency in the range of 0.1-10 Hz, which exist extensively in the terrestrial space, are associated with numerous intriguing space physical problems, including magnetosphere-ionosphereatmosphere coupling and radiation belt modeling.As a transitional region from atmosphere to magnetosphere and also an anisotropic medium with background parameters changing rapidly with height, the ionosphere has significant effects on the ionospheric propagation of ULF waves, which therefore is frequently investigated to reveal the mechanism of relevant problems.Tepley and Landshoff (1966) first proposed the waveguide theory for ionospheric propagation of ULF waves, which assumes that ULF waves propagate as shear Alfvén waves along the magnetic field line from low to high latitudes or propagate as compressional waves from high to low latitudes.Greifinger (1972) and Greifinger andGreifinger (1968, 1973) developed the theory of ionospheric waveg-X.Xu et al.: Numerical study of the generation and propagation of ULF waves uide, investigating the coupling, transmission, reflection and cutoff of ULF waves and the effects of propagation direction.The shear Alfvén waves and compressional waves propagate independently in the magnetosphere, while in the ionosphere they are coupled through Hall conductivity, which can affect the reflection and penetration of ULF waves (Yoshikawa andItonaga, 1996, 2000;Yoshikawa et al., 1999).Since the ionosphere is not an absolutely perfect conductor itself, ULF waves can penetrate through the ionosphere into neutral atmosphere then propagate as electromagnetic waves which can be observed on the ground.The 90 • phase rotation of magnetic perturbation from the magnetosphere to the atmosphere is called Hughes rotation and was thoroughly studied by Hughes (1974) and Hughes and Southwood (1976a, b).Following their previous research, Sciffer and Waters (2002) and Sciffer et al. (2004) presented an analytic solution to the problem of the propagation of ULF waves from the magnetosphere to the ground in the oblique background magnetic fields for a thin sheet ionosphere, neutral atmosphere and perfectly conducting ground.In order to investigate the temporal and spatial evolution of ULF waves, Lysak (1997Lysak ( , 1999Lysak ( , 2004) built a numerical model assuming a vertical magnetic field and uniform plasma to study the propagation features of ULF waves in the high-latitude ionosphere with different background parameters and another model including dipole geometry to analyze magnetosphere-ionosphere coupling by Alfvén waves at midlatitudes by performing a three-dimensional simulation.Also, a new model featuring finite ionospheric conductivity and capable of calculating the ground magnetic field of ULF waves is developed by Lysak et al. (2013) to study the ionospheric Alfvén resonator (IAR) in a dipolar magnetosphere.
Apart from research on the propagation of naturally excited ULF waves, the generation of ULF waves by ionospheric modulated heating is considered a physical problem of great importance and interest as well.A strong horizontal electric current driven by an atmospheric dynamo electric field and a magnetospheric electric field flows in the D or E region of the polar ionosphere, which is called the auroral electrojet.Similarly, the equatorial electrojet flows in the lower equatorial ionosphere, which is associated with Cowling conductivity and a tidal electric field.Modulated heating of the lower-ionospheric region where these currents flow with powerful high-frequency (HF) pump waves makes the ionospheric conductivity of the region change periodically, which in turn modulates the preexisting currents in the heated area (Moore, 2007).In the meantime, these oscillating currents form an electric dipole antenna radiating low-frequency waves in the ionosphere.This hypothesis was first suggested by Willis and Davis (1973) and was then proved when artificially excited low-frequency signals were detected in an experiment for the first time (Getmantsev et al., 1974).A series of experiments of generating low-frequency waves following this mechanism were carried out at the European Incoherent Scatter Scientific Association (EISCAT) and the High Frequency Active Auroral Research Program (HAARP; Cohen et al., 2012;Agrawal and Moore, 2012;Moore et al., 2007;Ferraro et al., 1984;Papadopoulos et al., 2003), and this mechanism was named PEJ (polar electrojet) (Stubbe and Kopka, 1977;Barr et al., 1991).The method of artificially generating low-frequency waves by modulating the lower ionosphere is completely dependent on the existence of quasi-stationary ionospheric currents, which to some extent limits the location of the heating facility to the high and equatorial latitudes and also makes the generation of ULF waves more unpredictable (Papadopoulos et al., 2011b).Moreover, some experimental observations such as ULF artificial excited signals at frequencies of 3.0, 5.0 and 6.25 Hz at Arecibo in 1985 still cannot be explained by classic ionospheric current modulation mechanism (Ganguly, 1986), which makes it necessary to develop new theories and experimental methods.
In a series of experiments conducted from 2009 to 2010 at HAARP, artificial ULF and lower ELF signals generated by the modulated heating ionospheric F region in the absence of electrojets were received on the ground far away from the heating facility, and the dependence on the heating conditions differs from the low-frequency waves generated by modulating the ionospheric currents (Papadopoulos et al., 2011a, b;Eliasson and Papadopoulos, 2012).This artificial generation of ULF waves relates to the oscillating diamagnetic drift current in the upper ionosphere due to the modulated heating, which is based on the ICD (ionospheric current drive) theory proposed by Papadopoulos et al. (2007).By modifying the model developed by Lysak (1997), Papadopoulos et al. (2011b) built a new model to study the ICD in the polar ionosphere by F region heating in cylindrical geometry.Eliasson et al. (2012) performed a theoretical and numerical study of the ICD based on a numerical model of the generation and propagation of ULF and ELF waves.The simulation results agree with the HAARP experimental measurements.Utilizing the ICD method, similar experiments conducted at SURA also received artificially generated ULF waves by modulated heating the F region.A comprehensive investigation of ULF wave properties and their dependence on modulation frequencies, polarization, beam inclination, receiving location and the geomagnetic activity was carried out by Kotik et al. (2013Kotik et al. ( , 2015)).
In this paper, we focused on the study of the effects of the background ionospheric parameters and the modulation frequencies on the process of modulated HF heating in the F region and the following generation and propagation of artificial ULF waves by using two mathematical models describing the above physical process.In Sect.2, we first briefly introduce the physical mechanism of the ULF wave generation in the F region by modulated HF heating proposed by Papadopoulos et al. (2011a) and Eliasson et al. (2012).Then we introduce the HF heating model and the following model of ULF wave generation and propagation.In Sect.3, firstly we discuss the possibility of the effects of the background ionospheric parameters and the modulation frequency, based on a relation among the radiation efficiency of the ring current source, the size of the spatial distribution of the modulated electron temperature and the wavelength of ULF waves; secondly, we run the simulation with the HF heating model to investigate the electron temperature response to the modulated heating under different ionospheric conditions with a modulation frequency sweep; thirdly, we run the simulation with the model of ULF wave generation and propagation in the same way to study the radiation efficiency of the ring current source, the propagation loss and the ground ULF wave magnitude.In Sect.4, conclusions based on the simulation results and the corresponding analysis are summarized.
2 Numerical model

Mechanism of artificial generation of ULF waves in ionospheric F region
The radial electron pressure gradient caused by the HF heating in ionospheric F region can drive a local diamagnetic ring current perpendicular to the ambient magnetic field given by where B 0 is the background magnetic field and ∇p e is the electron pressure gradient in the heated F region (Spitze, 1967).When the HF pump wave is amplitude modulated with the oscillation frequency f in the ULF band, the ring current will oscillate with the same frequency.The empirical evidence of the modification of the F region reveals that the change in electron temperature is larger than the electron density and that the response time of the electron density to the HF heating is much longer than that of the electron temperature (Robinson, 1989;Hansen et al., 1992a).Therefore, the electron pressure gradient mainly originates from the electron temperature gradient due to HF heating, and the pressure gradient is given by where n e is the electron density, k B is the Boltzmann constant and T e is the electron temperature.The oscillatory ring current can be expressed as (3) The ring current integrated over the heated volume will create an oscillating magnetic dipole moment parallel to the ambient magnetic field, acting as a virtual antenna which radiates ULF waves at the modulation frequency in the ionosphere.Some of ULF waves will penetrate into the Earthionosphere waveguide and will be received on the ground.

HF heating model
The plasma transport model for the F region ionospheric heating model can be described as follows (Bernhardt and Duncan, 1982;Shoucri et al., 1984;Hansen et al., 1992b): with the subscript α denoting electron (e) and ions (i), n α and m α the number density and mass of electron and ions, respectively, T e and T i the electron and ion temperature, k B the Boltzmann constant, ∂ ∂s the directional derivative along the ambient geomagnetic field, v e = v e • B 0 /B 0 the field-aligned flow velocity, and g = g • B 0 /B 0 the projection of gravity acceleration along the field line.
Equation ( 4) is the steady-momentum equation under the assumption of ambipolar diffusion and quasi-neutrality, and electron inertia can be neglected for the time-and space scales under consideration.The diffusion coefficient D can be expressed as (Hinkel et al., 1992) with the ionic species index i designating the ions involved in our simulation domain, i.e., O + , O + 2 and NO + .υ en and υ in are electron-neutral and ion-neutral collisional rates, respectively, both of which are calculated from Banks and Kocharts (1973) and Schunk and Nagy (1978).
Equation ( 5) is the continuity equation along the direction of the geomagnetic field line for ionospheric electrons, where Q e is the local ionization source term producing the equilibrium density profile without the HF heating and β e is the electron recombination rate.β e is given by (Hinkel et al., 1992;Schunk and Walker, 1973) with k 1 = 4.2 × 10 −7 (300/T e ) 0.85 cm 3 s −1 and k 2 = 1.6 × 10 −7 (300/T e ) 0.55 cm 3 s −1 .The ion densities of O + , O + 2 and NO + are obtained and updated according to the important photochemical equilibrium with the dissociative reaction rates and the source term of ion production rates dependent upon the background densities of neutrals, electron and ions (Park and Banks, 1974;Schunk and Walker, 1973;Schunk and Nagy, 1980;Hinkel et al., 1992).Q e is calculated from Eqs. ( 4) and ( 5) by setting ∂ ∂t = 0 with the initial background ionospheric density profiles (Shoucri et al., 1984), and it is used as a constant during our simulation.Equation ( 6) is the electron energy conservation equation along the geomagnetic field line, which includes the effects of convection and pressure flux and heat conduction.Also note that our simulation is based on the assumption that ion temperatures of O + , O + 2 and NO + and neutral temperatures of N 2 , O 2 , O and H e are the same and the temperature remains invariant during the simulation.
K e is the parallel coefficient of thermal conductivity which takes the form of (Banks and Kocharts, 1973) with N n being the densities of neutrals and Q n the cross section of the mean neutral-electron momentum transfer, which is calculated from Schunk and Nagy (1978).L e is the rate of electron cooling, the physical mechanism of which mainly contains the translational electron-neutral interactions, the rotational and vibrational excitation of N 2 and O 2 , the fine structure excitation of O, and the electronic excitation of O and O 2 .The expression of L e can be found in Schunk and Nagy (1978).
S 0 represents the steady-state source term in the absence of the heating and can be estimated by Eqs. ( 4) and ( 6) by ∂ ∂t = 0 and S HF = 0, which takes the form of The terms in Eq. ( 10) with an extra subscript "0" denote the value at t = 0 in our simulation, which is exactly calculated from the initial background profiles.S 0 is also kept invariant during our simulation like Q e .The profile of S 0 used in the simulation is shown in Fig. 6a.S HF is the energy absorption from heating waves.Since only thermal process is taken into account, S HF provides the external heating source for ionospheric changes.For the F region, the heating source term S HF includes both ohmic loss and anomalous loss due to the nonlinear wave-wave interaction and wave-particle interaction (Meltz et al., 1974;Perkins et al., 1974).The absorption of HF pump power can be expressed as (Shoucri et al., 1984;Gustavsson et al., 2010) with E ± (z) being the HF wave electric field, σ ± (z) the electric conductivity tensor, and the subscript "+" and "−" representing the right-hand circularly polarized O mode and lefthand circularly polarized X mode, respectively.
The wave field at altitude of z is calculated as (Shoucri et al., 1984;Gustavsson et al., 2010) with z 0 being the bottom altitude of the ionosphere, ε ± (z) the ionospheric dielectric tensor, k 0 the vacuum wave number and N ± (z) the ionospheric refractive index expressed as The real part µ ± (z) and the imaginary part χ ± (z) are related to σ ± (z) and ε ± (z) by and the exact expression of σ ± (z) and ε ± (z) can be found in Shoucri et al. (1984).
If ohmic heating is considered as the only mechanism when calculating the absorption of HF pump power by the ionosphere, the energy deposition of the ionosphere when the X mode is used as the pump wave is approximately 4 times larger than that of the O mode (Löfås et al., 2009).This is because the absorption of HF power is dependent on the imaginary part of the refractive index N ± (z) according to Eqs. ( 12) and ( 13), and χ − (z) of the X mode is larger than χ + (z) of the O mode.However, near the O mode reflection point of the F2 peak, huge anomalous absorption can be induced by parametric instability and thermal self-focusing instability when the O mode wave is used as the pump wave, which usually brings about much larger density and temperature enhancement than ohmic heating (Perkins and Valeo, 1974;Gurevich, 1986;Istomin and Leyser, 1997;Kuo, 2015).
In order to achieve a better HF heating effect, the pump wave of the O mode is utilized in our simulation and the frequency of the HF wave is adjusted so that the reflection point is near the F2 peak.At the reflection point in our simulation, the absorption of HF power is calculated as (Meltz et al., 1974;Perkins and Valeo, 1974) with S ohmic being ohmic loss, S anom the anomalous absorption caused by collisional and Landau damping (Perkins et al., 1974), ε 0 the vacuum dielectric constant, ω pe the electron plasma frequency at the reflection point, and ω 0 the frequency of the pump wave.The electron collisional frequency in S ohmic is calculated as with υ ei being the electron-ion collisional frequency at the reflection point; it is derived from Banks and Kocharts (1973) and Schunk and Nagy (1978), while the exact expression of υ e2 in S anom can be found in Perkins et al. (1974).The peak electric field amplitude E 0 is estimated by an empirical formula: with ERP the effective radiated power as a input parameter thus the propagation loss of the HF pump wave is not taken into account.The threshold field value E t for inducing the anomalous absorption is calculated from Kuo (2015).For simplification of the calculation, we adopt a 2-D Gaussian absorption model taking the form of with S HFmax being the HF power absorption at the reflection point calculated by Eq. ( 14), (x r , z r ), the coordinate of the reflection point in our simulation domain; σ x the horizontal half-width; and σ z the vertical half-width.According to the artificial generation of ULF waves in ionospheric F region experiment conducted at HAARP and SURA (Kotik and Ryabov, 2012;Papadopoulos et al., 2011b), the source term S HF is modulated by ULF square waves with 50 % duty cycle in our simulation.

Model of ULF wave generation and propagation
In our simulation, the generation model of ULF waves proposed by Eliasson et al. (2012) is used to calculate the following propagation of ULF waves excited by modulated HF heating.The model is built under the following assumption: firstly, the frequency of the artificially generated ULF wave is far less than the electron cyclotron frequency; secondly, O + is considered as the only ion dominating the propagation of the plasma wave and the condition of quasi-neutrality n 0 = n i = n e for ion and electron number densities is satisfied; thirdly, the electron-neutral and ion-neutral collisions are the main collision mechanism when calculating the wave dynamics.
In the ionosphere, the propagation of artificially generated ULF waves is calculated by the following equations: with E being the electric field, B the magnetic field, A the vector potential introduced by Eq. ( 20), E = ∇ϕ − ∂A ∂t with the gauge ϕ = 0, c the speed of light, ω pe = e 2 n 0 ε 0 m e the electron plasma frequency, v 2 A = B 2 0 µ 0 n 0 m i the Alfvén speed, e the electron charge, n 0 the number density of electron and µ 0 the permeability of free space.
H 1 , H 2 and H 3 are the background parameter matrices taking the form of the following expressions: with ω ce = eB 0 m e and ω ci = eB 0 m i being the ion and electron cyclotron frequencies, θ the dip angle of the ambient geomagnetic field, and υ en and υ in the electron-neutral and ion-neutral collisional rates, respectively, both of which can be calculated according to Banks and Kocharts (1973) and Schunk and Nagy (1978) as in the HF heating model.
In this model, the modulated HF heating effect is estimated by Eq. ( 21), in which δT e is the electron temperature deviation from the background value.When the geomagnetic field is assumed to be vertical, δT e is calculated as with δT emax being the oscillation amplitude of the electron temperature, ω the modulation frequency of HF heating, and D x and D z the half-width of the spatial distribution of the δT e in the horizontal and vertical direction, respectively.Note that Eq. ( 17) and Eq. ( 25) are similar in form yet their physical meanings are completely different.According to Eliasson et al. (2012), the term ∂∇δT e ∂t mainly generates the waves of non-propagating modes which are not studied in this paper.In the meantime, the magnetosonic wave is directly excited by the diamagnetic ring current due to the perpendicular electron temperature pressure gradient, while the Alfvén wave is generated by mode conversion from the magnetosonic wave through Hall conductivity.Thus, in our simulation, only the y component of Eq. ( 21) is considered, and it is expressed as At the boundary between the ionosphere and atmosphere, the boundary condition proposed by Eliasson et al. (2012) is applied.Also, the electromagnetic perturbation in the neutral atmosphere is estimated by analytic expression, and the inversely spatial Fourier transform is applied to obtain the real space wave field (Eliasson and Papadopoulos, 2009;Eliasson et al., 2012).

Parameter setting in numerical simulation
In this paper, simulations with both models introduced previously in Sects.2.2 and 2.3 are conducted in a twodimensional computational domain, with x as the horizontal direction and z as the vertical direction., the temperatures of electron, ions and neutrals, and the dip angle θ of the geomagnetic field.The simulation time is set to 02:00 and 14:00 LT, 10 March 2006, to study the effect of the daytime and nighttime ionosphere.The intensity of the background geomagnetic field is set to be 4.0 × 10 5 T.
The heating parameters for the modulated HF heating are listed in the Table 1, including foF2, the electron plasma frequency at the peak of the F2 layer; f HF , the frequency of pump wave; h r the reflection height of the pump wave; the electron temperature, T e ; and the number density, N e , at h r .Figure 1 shows the electron density and temperature profiles against height from 90 to 1000 km in the ionosphere used in the simulation.Now we introduce the setting of the simulation domain.In the simulation with the HF heating model, the spatial grid size in the x and z direction is 1.0 km, and the time resolution is 1.0 × 10 −3 s.The spatial range is from −100 to 100 km in the x direction and from 180 to 400 km in the z direction.In the simulation with the wave generation and propagation model, the spatial grid size in the x and z direction is 10.0 km, and the time resolution is 1.0 × 10 −4 s.The spatial range is from −1000 to 1000 km in the x direction and from 0 to 1000 km in the z direction, with the ionosphere from 90 to 1000 km and neutral atmosphere from 0 to 90 km. 3 Simulation results and discussion

ULF wave generation and ring current radiation source
According to the physical mechanism of the generation of ULF waves in the ionospheric F region, a diamagnetic ring current source of ULF radiation is driven by the modulated HF heating, which is of crucial importance in the whole wave generation process.So, the relation between the wave generation and the source of ULF radiation due to HF heating should be qualitatively discussed before presenting the simulation results.
Since in the simulation with the model of ULF wave generation and propagation, the electron temperature response to HF heating is estimated with Eq. ( 25), it is natural that both the amplitude of the electron temperature oscillation and the spatial distribution of the δT e are expected to affect the generation of ULF waves.An analytical mathematical deduction by Vartanyan (2015) based on Eqs. ( 18) and ( 19) has revealed the relation between the horizontal half-width of the spatial distribution of the δT e and the radiation efficiency of the ring current source to some extent.This deduction is carried out under the following assumptions: firstly, the geomagnetic field is in the z direction and the ionosphere is uniform, with parameters such as electron density and Alfvén speed constant; secondly, the wave frequency is much lower than the ion cyclotron frequency; thirdly, the ring current source driven by modulated HF heating only radiates the ULF waves in the x direction, which is perpendicular to the ambient geomagnetic field.Based on these assumptions, we find that only magnetosonic waves are considered in this model since the mode conversion between the magnetosonic waves and shear Alfvén waves cannot exist in a uniform ionosphere.Also, the deduction finally turns to looking for a one-dimensional function relating the y component of the vector potential to the shape of the distribution of δT e after simplification, which can be achieved by solving Eq. ( 27): The source term δT e can be estimated by Eq. ( 25) as well but with z = 0. Using the method of the Green function and simply considering the far-field magnitude of the ULF radiation (x D x ), the solution to this equation can be expressed as with k being the wave number of the ULF wave and A y designating the ground far-field magnitude of the ULF radiation.
After substituting k with k = 2π λ , in which λ is the wavelength of the ULF wave, we can obtain the relation among A y , λ and D x , which takes the form of Eq. ( 29) indicates that the ground magnitude A y primarily depends on the ratio of λ to D x .Moreover, A y is less than an order of 10 −5 of the intensity of ULF radiation source when λ ≤ D x , and the magnitude begins to increase with the ratio λ/D x when λ > D x .This magnitude growth is almost linear, in the range of 1 < λ/D x < 5. Similarly, we can conduct the deduction in the vertical direction and find that the vertical radiation efficiency is dependent on the ratio of the wavelength λ to the vertical half-width D z of the δT e spatial distribution.Judging from the above analysis, it is clear that the efficiency of the ULF radiation is higher when the size of the δT e spatial distribution is smaller.Also, the intensity of the ULF ring current radiation source depends on the energy the whole modulated heating electron contains, which can be given by W e = n e k B δT e dxdz.When the amplitude of the electron temperature oscillation is constant, the intensity of the ULF radiation source monotonously increases with the area the modulated heating electron covers.Since both the intensity of the radiation source and the far-field radiation efficiency contribute to the ground wave field magnitude, it is seemingly tricky to directly discuss the dependence of the amplitude of the ground wave field on the size of the spatial distribution of δT e .As a matter of fact, if we think about this problem in terms of engineering, a larger amplitude of the ground wave field can be obtained by adjusting the heating parameters of the HF transmitter with a certain ionospheric profile, which is not much of a difficulty.
However, what is truly uncontrollable is the effects of the background ionospheric profiles on the generation and propagation of artificial ULF waves.Firstly, during the process of modulated HF heating, both the modulation frequency and background ionospheric profile are expected to have an impact on δT emax , the amplitude of electron temperature oscillation, and the shape of the δT e spatial distribution, based on the HF heating model.Secondly, since the wavelength of the ULF wave is decided by λ = 2πv A ω and the Alfvén speed is associated with the ionospheric electron density, the ionospheric profile and the modulation frequency are both expected to affect the radiation efficiency of the diamagnetic ring current source according to Eq. ( 29).Thirdly, the propagation characteristics and propagation loss of the artificial ULF wave, after being generated in the ionosphere, are related to both the modulation frequency and the background ionospheric profile, based on the generation model of the ULF waves.
These factors jointly contribute to the distribution of ground magnetic field intensity, which makes it difficult to find an exact function to describe the relation among them.However, the effects of the background ionospheric parameters and the HF heating modulation frequency can be investigated with the following simulation results.

Electron temperature response to modulated heating
This section is dedicated to studying the effects of the background ionospheric parameters and the modulation frequency on the electron temperature response due to modu-  lated heating.In this simulation, the parameters of the HF power absorption Eq. ( 17) are given by σ x = 20 km and σ z = 40 km, and the ERP is set to be 800 MW.
To begin with, we focus on the spatial distributions of electron temperature change T e with different background ionospheric profiles when the modulation frequency is invariant.The contours of electron temperature changes due to modulated HF heating with a modulation frequency of 0.5 Hz in the nighttime ionosphere are presented in Fig. 2, and the results with the same heating conditions in the daytime ionosphere are presented in Fig. 3.The top, middle and bottom row correspond to the cases where the locations of HF heaters are HAARP, Wuhan and Jicamarca, respectively.
The left, middle and right column represent the T e spatial distribution at the heating time of 0.01, 1.0 and 2.0 s.
Figure 2d, e and f show the change in the spatial distribution of T e in the nighttime ionosphere at Wuhan during a complete modulation period.The HF transmitter is turned on from 0 to 1.0 s and turned off from 1.0 to 2.0 s.At the time of 0.01 s, the energy of T e concentrates completely around the reflection point.Unlike the low-ionosphere HF heating, the transport process becomes quite important during the modulated heating of the ionospheric F2 region plasma, as illustrated in the HF heating model in Sect.2.1.As time proceeds, the spatial distribution of T e manifests a distinct extension along the tilted geomagnetic field due to the transport process, at the time of 2.0 s.However, the half-width of the T e spatial distribution across the field does not show evident change since the conduction and diffusion of heat is much faster and stronger along the field than across the field.Comparing the simulation results in nighttime conditions in Fig. 2 with the daytime condition ones in Fig. 3, we find that the "extension" feature of the T e spatial distribution is relatively less obvious in daytime conditions than in nighttime conditions.Another interesting phenomenon in the HAARP case is that the extension of the T e spatial distribution is faster below the reflection point in the nighttime ionosphere, while it is faster above the reflection in the daytime ionosphere.This is because K e , the parallel coefficient of thermal conductivity per electron, is higher below the reflection point in the nighttime ionospheric HF heating, while the situation is exactly the opposite in the daytime condition, as illustrated in Fig. 6b.Now, we extend the discussion in Sect.3.1 about the relation between the radiation efficiency of the ULF ring current source and the size of the T e spatial distribution to the cases where dip angles of the geomagnetic field are included.The vertical and horizontal half-widths correspond to longitudinal half-width along the field and the transversal half-width across the geomagnetic field, respectively, in this section.Following this conclusion in Sect.3.1, the increase in the longitudinal half-width of the T e spatial distribution will compromise the radiation efficiency along the geomagnetic field.In terms of this explanation, the radiation efficiency of the ULF ring current source is supposed to decline when the latitude of the heating location is increasing since the longitudinal half-width of the T e spatial distribution increases with the latitude according to the results in Figs. 2 and 3. Also, modulated HF heating in the daytime ionosphere can drive a ULF radiation source with larger radiation efficiency than in the nighttime because of the smaller size of the T e spatial distribution.
Since δT emax , the amplitude of the electron temperature oscillation, is a key index measuring the intensity of the ULF radiation source due to the modulated heating, the effects of the ionospheric parameters and the modulation frequency on δT emax should also be studied, which can be achieved eas- ily by investigating the electron temperature at the reflection point.A series of results of the temporal evolution of the electron temperature change at the reflection point are presented in Fig. 4, and the variation of the average electron temperature oscillation amplitude with different modulation frequencies is presented in Fig. 5.The modulation frequency sweep of 0.5, 1.0, 2.0, 4.0 and 8.0 Hz is used in the simulation.
Since in the ionospheric F2 layer the electron cooling time is tens of seconds due to the dramatically decreasing collision rates and the dominating transport process (Gurevich, 1986), the electron temperature deviation from the initial background value cannot return to zero in the cooling phase of a modulation period.Moreover, the oscillation of T e will finally reach a saturation state, as the results from Fig. 4 indicate.Also, the amplitude of the oscillation of T e is not a fixed value due to the drastic increase at the first modulation period.So we make an average of the amplitude of the oscillation of T e with the value of the first period excluded, the results of which are shown in Fig. 5.When the ionospheric profile is fixed, we find that the average amplitude of the oscillation of T e is approximately inversely proportional to the modulation frequency.Especially, when the modulation frequency is high in a ULF band, such as 8.0 Hz in the simulation, the oscillation amplitude becomes so small that 1 2 3 4 5 6 7 8 K e /N e (10 -14 J m -1 s -1 K -1 ) 0 0.5 1 2 3 1.5 2.5 (S HF +S 0 -L e )/N e (10 -20 J s -1 ) 0 0.5 1 2 3 1.5 2.5 L e /N e (10 -21 J s -1 ) HAARP/02:00 LT HAARP/14:00 LT Wuhan/02:00 LT Wuhan/14:00 LT Jicamarca/02:00 LT Jicamarca/14:00 LT the ULF wave generated due to the modulated HF heating is probably undetectable on the ground.This dependence of the average amplitude of the oscillation of T e on the modulation frequency is expected since more pump wave energy will be absorbed by the ionosphere, and the electron temperature can take enough time to cool down when the modulation frequency is lower.According to Sect.2.1, the HF heating model is a selfconsistent simulation model.Also, the background parameters such as plasma number density and temperature affect the electron temperature change jointly.So it is not advisable to directly relate the evolution of T e to a single factor from the background parameters.However, we can analyze the amplitude of the oscillation of T e based on the source terms of the equations in the HF heating model.It is without doubt that the change in the electron temperature at the reflection point in a certain modulation period primarily depends on the source terms of Eq. ( 6).What is more, the change in electron density can be negligible in a modulation period due to its much longer relaxation time than the electron temperature.So the term ∂T e ∂t representing the changing rate of the electron temperature is related to the term 2(S HF +S 0 −L e ) 3n e k B in the heating phase of a modulation period and related to the term 2L e 3n e k B in the cooling phase according to Sect.2.1.The total electron heating rate profile, S HF + S 0 − L e , per electron and the electron cooling rate profile, L e , per electron are shown in Fig. 6c and d, at the modulation frequency of 0.5 Hz.The values of the term (S HF + S 0 − L e )/N e at the re-flection point in the nighttime ionosphere are 1.45 × 10 −20 , 2.22×10 −20 and 2.60×10 −20 J s −1 , which correspond to the HAARP case, the Wuhan case and the Jicamarca case, respectively.In the daytime ionosphere the corresponding values are 1.70 × 10 −20 , 1.43 × 10 −20 and 1.24 × 10 −20 J s −1 .Comparing these values, we find that in the nighttime modulated heating, the total heating rate per electron increases with the latitude at which the HF heating is conducted, while in the daytime it decreases with latitude.Moreover, in the midlatitudes and the equator region, the total heating rate per electron is larger in the nighttime ionosphere, while in the high-latitude region it is larger in the daytime ionosphere.The cooling rate per electron manifests the same characteristics as shown in Fig. 6.Since the amplitude of the oscillation of T e is related to ∂T e ∂t , we can conclude that δT emax is affected by the ionospheric profiles in the same way as the total heating rate per electron.This conclusion is obviously supported by the results in Fig. 5.

Propagation and generation of ULF waves vs. modulation frequency and ionospheric parameters
In this section, we discuss the effects of modulation frequency and ionospheric parameters on the propagation and generation of ULF waves based on the simulation results from the model of ULF wave generation and propagation introduced in Sect.2.2.In our simulation, the parameters of the electron temperature deviation Eq. ( 25) are δT emax = 500 K, D x = 40 km and D z = 40 km, which are fixed in all the cases of the simulation.The magnitude of the total magnetic field vector of the artificially generated ULF waves at t = 2.0 s in the simulation presented in Fig. 7, in which the magnitude of the magnetic field is calculated as |B| = B 2 x + B 2 y + B 2 z .The top, middle and bottom row represent the HAARP case, Wuhan case and Jicamarca case, respectively.The left column corresponds to the generation and propagation of ULF waves in the nighttime ionosphere, while the right column corresponds to the daytime case.The modulation frequency is 4.0 Hz in the simulation of Fig. 7.
Comparing the nighttime case and the daytime case, we find that the amplitude of the magnetic field of the ULF waves in the heating region is larger in the daytime ionosphere than in the nighttime ionosphere, despite the effect of the dip angle.Since all the cases contain the same δT e input according to our parameters setting, we can expect that the efficiency of the energy conversion from the electron temperature oscillation to the ULF waves is higher in the daytime ionosphere in the heating region.However, when it comes to the comparison of the radiation efficiency between cases in the nighttime and daytime ionosphere, the situation turns out to be the other way around.As illustrated in Fig. 7, more than 90 % of the energy of the artificially excited ULF waves in the heating region can spread in the whole simulation domain Figure 7.The magnitude of total magnetic field vector of the artificially generated ULF wave at t = 2.0 s in a simulation.The modulated HF heating is conducted in the ionosphere at HAARP, Wuhan and Jicamarca at 02:00 LT and 14:00 LT, and the modulation frequency is 4.0 Hz.
through wave propagation in the nighttime ionosphere, while in the daytime ionosphere, the wave energy is severely constrained in the heating region, with less than 50 % of the energy getting out of the region.This can be partially explained by the conclusion from Sect.3.1.According to the definition of ULF wave wavelength in Sect.3.1 and the Alfvén speed in Sect.2.2, we can expect that λ ∝ 1 n e when the modulation frequency is invariant.As depicted in Fig. 1a, the electron density is larger in the daytime ionosphere, so the wavelength of ULF wave is relatively shorter.Based on the conclusion on the relation between the radiation efficiency and the wavelength in Sect.3.1, the radiation efficiency is supposed to be lower in the daytime ionosphere than in the nighttime ionosphere.
Now we focus on how the background parameters and the modulation frequency affect the propagation of ULF waves and their penetration to the ground.Although we can draw some conclusion from the simulation with the realistic geomagnetic dip angles and the ionospheric profiles varying together, there is no doubt that the joint effects of these two factors can make the problem more complicated and puzzling.In order to simplify the problem without omitting the parameters we are interested in, we can focus on the HAARP case but with the dip angle assumed to be 90 • .The following discussion is based on the simulation with this assumption.Also, the modulation frequency sweep of 2.0, 4.0 and 8.0 Hz is used in this simulation.
In order to investigate the effects of the ionospheric profiles on the propagation loss of the ULF waves, we introduce the magnetic energy of the ULF waves at a certain height within the range from −100 to 100 km in the horizontal direction, which is calculated as W B = |B| 2 2µ 0 dx.Usually we attribute the energy loss of ULF waves in the high-latitude ionosphere to Joule heating because of the existence of the Pedersen current J P = σ P E 0 (Kelley, 2009), with σ P the Pederson conductivity and E 0 the stable ionospheric electric field at high latitude.So the energy loss can be estimated by J P •δE, with δE the electric field of the ULF waves.Both the daytime and nighttime ionospheric Pederson conductivity profiles calculated from the parameters used in the simulation are perceptible mainly from 200 km in the F region to the bottom of the ionosphere at 90 km, as shown in Fig. 8.So we use the magnetic energy of ULF waves at 90 and 200 km to study the propagation loss.The temporal evolution of the magnetic energy of the artificially generated ULF wave at a height of 90 km (solid line) and 200 km (dashed line) in the daytime (blue) and nighttime (red) ionosphere is presented in Fig. 9, in which the modulation frequency in this case is 2.0 Hz.The left y axis shows the value at 200 km, while the right y axis shows the value at 90 km.We find that the magnetic energy at 200 km increases from zero and reaches a steady oscillation at the modulation frequency, while the magnetic energy at 90 km is not stable compared with that at 200 km.Since 90 km is the boundary of the ionosphere and atmosphere where the mutual conversion of ULF waves and electromagnetic waves in free space happens, this unstable oscillation of the wave energy at 90 km is to be expected.What is more important is the difference between the wave energy in the daytime and nighttime ionosphere.As depicted in Fig. 9, the magnetic energy of ULF waves at 200 km in the daytime ionosphere is much larger than in the nighttime ionosphere, while the wave energy at 90 km in the daytime ionosphere is nearly the same as in the nighttime ionosphere.This feature indicates that a larger propagation loss is expected in the daytime ionosphere.To present the energy loss intuitively, we make a time average of the magnetic energy at 90 and 200 km, denoting them, respectively, as W B90 km and W B200 km , and calculate the ratio of W B90 km /W B200 km .In this way, a smaller ratio indicates a larger magnetic energy loss.The results are shown in Ta-    2 is that the wave energy losses in the daytime and nighttime ionosphere both decrease when the modulation frequency increases.Finally, we discuss the feature of the artificially excited ULF waves penetrating into the ground.The ground magnitudes of the ULF waves at the modulation frequency of 2.0, 4.0 and 8.0 Hz in the daytime and nighttime ionosphere are demonstrated in Fig. 10.We find that both in the daytime and nighttime ionosphere, the amplitude of the magnetic field of the ULF wave decreases when the modulation frequency is enhanced from 2.0 to 8.0 Hz.According to the analysis in Sect.3.1, the wavelength is shorter when the modulation fre- quency is raised, which makes the radiation efficiency of the ring current source lower.Moreover, the higher modulation frequency causes more propagation loss in the ionosphere as indicated in Table 2.Both factors contribute to the smaller ground magnetic field amplitude at the higher modulation frequency.When the modulation is fixed, the ground magnitude of the ULF magnetic field is larger when the wave generation and propagation happen in the nighttime ionosphere than in the daytime ionosphere.Also, this difference of the ground amplitude of the ULF wave due to the nighttime and daytime ionospheric profiles expands when the modulation increases.

Conclusions
Based on the HF heating model and the model of ULF wave generation and propagation, we investigate the effects of the background profiles and the modulation frequencies on the process of modulated HF heating in the ionospheric F region and the subsequent generation and propagation of the artificially generated ULF waves.Some conclusions can be drawn, as follows: 1.The magnitude of the artificially generated ULF wave is related to the intensity of the ring current source driven by F-region-modulated HF heating and its radiation efficiency.The source intensity depends on both the spatial distribution size and the oscillation amplitude of the modulated electron temperature, δT e , due to HF heating.
The radiation efficiency mainly depends on the spatial distribution size of δT e with a fixed ULF wavelength, and a smaller size predicts a higher radiation efficiency.
2. The size of the spatial distribution of δT e is larger in the nighttime ionosphere than in the daytime ionosphere, and it is smaller in the ionosphere at low latitudes than at high latitudes, which indicates a higher radiation efficiency of the ring current source in the daytime ionosphere at low latitudes with a fixed ULF wavelength.
3. The background ionospheric profiles can affect the absorption of the pump wave and the cooling rate of the electron temperature during the modulated HF heating, thus determining the oscillation amplitude of δT e at a fixed modulation frequency.The oscillation amplitude of δT e increases with latitude in the nighttime ionosphere, while it decreases with latitude in the daytime ionosphere.Also, the oscillation amplitude of δT e is approximately inversely proportional to the modulation frequency.
4. The radiation efficiency of the ULF ring current source is larger in the nighttime ionosphere than in the daytime ionosphere regardless of different geomagnetic field dip angles, while the energy conversion efficiency from electron temperature oscillation to ULF waves is lower in the nighttime ionosphere with a fixed δT e input and a fixed modulation frequency.
5. The daytime ionosphere produces more energy dissipation during the propagation of artificially generated ULF waves due to Joule heating, and this propagation loss is larger when the modulation frequency is raised.
6.The ground magnitude of the magnetic field of the artificial ULF wave is larger in the nighttime ionosphere than the daytime ionosphere, and the difference between daytime and nighttime conditions expands as the modulation frequency increases.Also the ground ULF wave amplitude decreases when the modulation frequency increases.

Data availability
Background parameters of our numerical simulation used in this paper are from the International Reference Ionosphere (IRI-2012) model and the neutral atmosphere model (NRLMSISE-00).These data can be accessed at the following websites: http://omniweb.gsfc.nasa.gov/vitmo/iri2012_vitmo.html and http://ccmc.gsfc.nasa.gov/modelweb/models/nrlmsise00.php.

Figure 1 .
Figure 1.The background ionospheric electron density and temperature profile used in our numerical simulation.

Figure 2 .Figure 3 .
Figure 2. Contours of electron temperature changes with modulation frequency of 0.5 Hz at heating time of 0.01, 1.0 and 2.0 s in the nighttime ionosphere.

Figure 4 .
Figure 4. Temporal evolution of electron temperature change at the reflection point.

Figure 5 .
Figure 5.The variation of the average electron temperature oscillation amplitude with different modulation frequencies.

X
. Xu et al.: Numerical study of the generation and propagation of ULF waves

Figure 6 .
Figure 6.(a) The steady-state source term S 0 profile per electron;(b-d) profiles of the parallel coefficient of thermal conductivity K e , the total electron heating rate S HF + S 0 − L e and the electron cooling rate L e per electron at heating time of 1.0 s with modulation frequency of 0.5 Hz.

Figure 8 .
Figure 8.The Pederson conductivity profiles in the daytime and nighttime ionosphere at HAARP used in our simulation.

Figure 9 .
Figure 9. Temporal evolution of the magnetic energy of the artificially generated ULF wave at the modulation frequency of 2.0 Hz at a height of 90 and 200 km in the daytime and nighttime ionosphere.

Figure 10 .
Figure10.The ground magnitude of magnetic field of the artificially generated ULF waves at the modulation of 2.0, 4.0 and 8.0 Hz in the daytime and nighttime ionosphere.
The background profiles used in the simulation are one-dimensional functions with height in z direction and remain constant in x direction.It should be noted that the HF heating model and the ULF wave generation and propagation model share the same background profiles in this paper.In our simulation with both models, HAARP (62.39 • N, 145.15 • W), Wuhan (30.52 • N, 114.32 • E) and Jicamarca (11.95 • S, 76.87 • W) are chosen as the locations of the HF heater.The according ionospheric and atmospheric background profiles are given by the International Reference Ionosphere (IRI-2012) model and the neutral atmosphere model (NRLMSISE-00).The profiles include the number densities of electron, O + , O + 2 , NO + , N 2 , O 2 , O and H e

Table 1 .
Heating parameters used in the model.

Table 2 .
Ratio between the magnetic energy of the artificial ULF wave at 90 and 200 km in the nighttime and daytime ionosphere at the modulation frequency of 2.0, 4.0 and 8.0 Hz.