Phase Velocity Analysis of Projected Wave Motion Along Oblique Radar Beams-A Numerical Study of Type-1 Radar Echoes

The nonlinear features of type-1 radar echoes were studied by a two-dimensional (2D) simulation of saturation the Farley-Buneman (FB) wave. The behavior of this FB wave in a plane perpendicular to the Earth’s magnetic field was simulated with a two-fluid code in which electron inertia was discounted while ion inertia was retained. It showed that the appearance of secondary waves propagating vertically and obliquely as the primary horizontal FB wave saturates. The secondary waves originating from nonlinear saturation process will construct the evolution of 2D modes which can be observed by oblique radar beams. We carried out the statistical analysis of projection phase velocities of 2D modes along oblique radar beam at different radar elevation angles. The result revealed that a likely density gradient effect of type-1 radar echoes for the wavelength dependence of phase velocity would appear at a larger radar elevation angle while short wavelength waves would approach isotropic speeds close to ion acoustic speed. This interesting result is primarily attributed to the spectral features of 2D modes.


INTROducTION
Type-1 radar echoes revealed the presence of plasma wave in the earth's equatorial and auroral electrojets.There are three primary sources of type-1 radar echoes.First, in the upper electrojet where there is no appreciable plasma density gradient, a pure FB wave traveling along the electron drift flow will develop when the electron drift speed exceeds the threshold of the Farley-Buneman (FB) instability factor (Kudeki et al. 1987;Pfaff et al. 1987a, b).Second, in the lower electrojet, plasma density gradients can develop kilometer large-scale gradient-drift (GD) waves at an electron drift speed below the threshold of FB instability wherein the large electric field originating from these horizontally propagating large scale primary waves could excite vertically traveling FB waves (Sudan et al. 1973;Kudeki et al. 1987;Pfaff et al. 1987a, b).Third, in the mid-latitude E region, a plasma structure with sharp lateral boundaries will induce an intensely polarized electric field to generate strong electron drift along radar beam direction such that type-1 waves can be excited (Haldoupis et al. 1996).In principle, the linear theory of FB instability can explain many of the observed behaviors of type-1 radar echoes, but there are many nonlinear features taken from other observations which linear theory cannot account for.For example, the analogous phase velocity spectral peak at different elevation angles (Cohen and Bowles 1967;Hanuise and Crochet 1981), the spectral asymmetry (Cohen and Bowles 1967;Farley et al. 1978;Kudeki et al. 1985Kudeki et al. , 1987;;Swartz 1997a, b), and echo power drops 0.3 dB/degree as the elevation angle from 0° (horizontal) (Ierkic et al. 1980;Swartz 1997b) demonstrates such a discrepancy.To explain these features, one must understand the behavior of nonlinear saturated FB waves, the relevant subject has been studied for four decades (Weinstock and Sleeper 1972;Rogister and Jamin 1975;Sudan 1983a, b;Kudeki et al. 1987;Kudeki and Farley 1989;Hamza and St-Maurice 1993a, b;Oppenheim and Otani 1996;Otani andOppenheim 1998, 2006;Fern andKuo 2001, 2009a).In the last decade, numerical simulation has become an important tool in exploring the nonlinear features.
However, owing to the limitations of computer resource, the past numerical simulations mostly focused on small-scale, pure, FB waves.The result of 2D saturation simulations of pure FB waves (Oppenheim and Otani 1996;Fern and Kuo 2001) revealed that, vertically propagating secondary waves emerged as horizontally traveling primary Farley-Buneman wave saturated; The secondary waves evolved from mode-mode coupling; The spectral asymmetry was closely related to the ambient electric field, and the phase velocity of these longer wavelength waves was below the prediction of linear theory, and so on.Oppenheim et al. (2008) firstly used supercomputer and Electrostatic Parallel PIC (EPPIC) code to carry out large-scale simulations of 2D fully kinetic FB turbulence, in a box size reaching about 160 m × 160 m.And their result showed that the evolution of long wavelength waves and the maximum of turbulent energy at saturated stage always shifted toward the longest wavelength waves allowed by the simulation box size.
In addition, if the density gradients exist along the ambient electric field which is slightly weaker than the FB instability threshold electric field, FB waves can also be excited.Farley and Fejer (1975) attributed the type-1 spectrum broadening to the effect of the density gradient term of the FB-GD instability theory.The observed evidence of a density gradient effect on type-1 waves was presented by Hanuise and Crochet (1981) who found the lower type-1 Doppler velocities as their radar observation shifted from higher to lower radar frequencies.Haldoupis et al. (2005) argued against a density gradient effect on short-scale FB waves for that they found no gradient effects on lowering the FB instability threshold and the wavelength dependence of phase velocity by the observation of two radar frequencies at 144 and 50 MHz.They also pointed out that some statistical studies using HF radar in the aurora region whose researchers also reached the same conclusion (e.g., Hanuise et al. 1991;Milan et al. 1997Milan et al. , 2001;;Lacroix and Moorcroft 2001).To clarify density gradient effects on FB waves, Fern and Kuo (2009b) performed 2D simulations of plasma density gradient effect on saturated FB waves and found that density gradient effects were related to the thickness of an unstable layer and ambient vertical electric field (which was parallel to vertical density gradient distribution).
Therefore, a series of 2D simulations of saturation FB waves revealed that the evolution of 2D modes is an important factor in studying the nonlinear features of type 1 radar echoes.Actually, the radar data observed by oblique radar beams can involve the information of 2D or 3D wave motion, because Kuo et al. (2007) performed an analysis of projected wave packet motion along oblique radar beams to study QP echoes.In this paper, we will use the similar analysis of projected wave motion along oblique radar beams to study type-1 radar echoes.Type-1 radar echoes observed by oblique radar beams come mostly from the saturation process of pure FB waves which involve the evolutions of 2D modes, thus the evolution characteristics of 2D modes will be studied in detail and a statistical analysis of projection phase velocities of 2D modes along radar beams will also be performed.It is interesting to note that the statistical results are compatible with the features of type-1 radar echoes.
For example, the statistical results for a larger radar elevation angle can display a likely density gradient effect on the wavelength dependence of phase velocity while shortwavelength waves still approach isotropic phase velocity.The detail results and theoretical analysis will be presented in the following sections.

Numerical Model and data Analysis
Past 2D simulations of pure FB waves such as particle code (Janhunen 1994) or hybrid code (Oppenheim and Otani 1996) focused mostly on smaller scale waves.Recently, Oppenheim et al. (2008) developed larger scale waves of particle code, but the simulations need enormous computer resources including use of a supercomputer.In fact, Pfaff et al. (1987b) compared the growth rate of FB waves calculated by Schmidt and Gary (1973) with the spectrum of electric field fluctuations measured in the CONDOR rocket flight and found that the fluid theory is accurate for FB waves with a wavelength larger than 3 m, and the kinetic theory should be taken into consideration for the waves with a wavelength smaller than 2 m.Based upon the consideration of larger scale wave and computer resources our numerical model uses a fluid code described in a previous paper (Fern et al. 2001).For the convenience of busy readers, we put the numerical model in Appendix A. The numerical computations were performed on a two-dimensional Cartesian mesh using grid points Nx × Nz in the x-z direction.The lengths X (east-west) and Z (vertical) of the simulation box matching grid points were assigned case by case to study the waves of different scales in a fixed resolution.A periodic boundary condition was imposed on both electron density n and electric potential Φ in both x-and z-directions.At t = 0, the electrons were set to move uniformly at a drift velocity of e V (i.e., the E B # drift velocity D V ); the ions assumed a constant velocity , where i X and v in were the ion-gyro frequency and ion-neutral collision frequency respectively 0 E was the background driving electric field and B 0 was the magnetic field strength.Then a density perturbation with amplitude sin n n 0 , d = x ^h was superposed on the background density where , is the wave-number to be assigned.A complete simulation will produce a data set of perturbation density , , n x z t d = ^h for analysis.Via a series of computations, the perturbation density , , n x z t d = ^h at each grid point at every time step was obtained.In addition to plotting the 2D density contour maps at some simulation time step to reveal the 2D evolution of plasma density, the temporal-and spatial-Fourier analysis of plasma density was also performed to obtain the information of 2D modes such as a power spectral distribution in frequency domain.The detailed 2D Fourier analysis was described by Fern et al (2001).

The Phase Velocity Analysis Method for Wave Motion Along Radar Beams
According to a power spectral distribution of every 2D spatial mode, we can evaluate the mean frequency f km which is close to the spectral peak.Depending on the horizontal and vertical wave modes k and m of every 2D mode (k, m), the horizontal phase velocity v px and vertical phase velocity v pz can be estimated roughly as follows: where k k X 2 represent the horizontal and vertical wave number respectively; θ is the polar angle of phase velocity of the corresponding mode.
If the radar beam pointing obliquely to elevation angle {, the way is more convenient to describe the wave propagation in a 2D Cartesian coordinate system perpendicular to magnetic field.Following the previous introduction of 2D simulation geometry in section 2.1, the schematic plot of 2D Cartesian coordinate system associating with the radar beam can be shown as Fig. 1, where the oblique radar beam is always perpendicular to the magnetic field and the aspect angle problem can be discounted.Let the horizontal wave number vector be in the positive x-direction the wave number can be expressed in both the Cartesian coordinate system and the polar coordinate system, where k r is the wave number component along the radar beam, { is the elevation angle of the oblique radar beam.The transformations of the wave number components between the two coordinate systems are given by and inversely, Since the phase velocities are inversely proportional to the wave numbers, i.e., v k , the transformations among the phase velocity components can be obtained from Eqs. ( 7) and (8).
The radial phase velocities v pr along the radar beam are directly observed by the oblique radar beam.If λ is a radar wavelength and the projection wavelength λ ij of the (k i , m j ) mode along the radar beam is given by we can take all of the wave modes with wavelength in the range λδλ < λ ij < λ + δλ and calculate the mean values of projection phase velocity v pr along the radar beam at elevation angle { using the following equation cos cos where θ ij is the propagation angle of the (k i , m j ) mode, P ij is the mean power of the (k i , m j ) mode with a wavelength in the range λδλ < λ < λ + δλ for a interval time and P ij cos(θ ij -{) corresponding to the projection power of the (k i , m j ) mode along the radar beams.Based on the limited 2D modes resulting from simulation space size and space resolution, we consider a small statistical range of 2D modes with a wavelength in the range λδλ < λ < λ + δλ for a given radar wavelength λ which means that a radar system can receive the echoes from the high power waves with a wavelength adjacent to λ.It is shown in section 3.3 that using δλ = 0.1λ or 0.2λ did not have significant different results.

The 2d Mode Analysis for Saturation Simulations of Pure FB Waves
In order to understand the features of type-1 radar echoes such as 2D wave motions from oblique radar beams, a detailed analysis of 2D mode evolution derived from the saturation process of pure FB waves is necessary.First, a simulation box with length X = 27 m and length Z = 108 m combined with limited grid points with Nx = 121 points and Nz = 121 points was designed for simplifying the 2D mode analysis.We then introduced a perturbation of 0.1% [i.e., δn = 0.001N 0 sin(kx)] with single wave mode k = 9 corresponding to a horizontal wavelength of 3 m and considered a downward nighttime electric field Ez = -12.8mV m -1 (corresponding to the electron drift velocity V e = 460 m s -1 ), where this perturbation type is called an initial-single waveperturbation.Other background parameters are listed in Table 1, where the threshold drift velocity V e th for driving FB instability is 437 m s -1 obtained from V e th = Cs(1 + ψ) and Cs = 356 m s -1 .In other words, our simulation conditions can drive pure FB waves.From the time evolution of standard deviation of density as shown in Fig. 2, we can see the perturbation grew and finally reached saturation.Via the analysis of the 2D gray scale map of density variation at different time stages, 2D nonlinear development was displayed in Fig. 3. Figure 3 showed that the primary wave gradually bent and nonlinearly developed secondary waves traveling at vertical direction and finally reached disordered structure.This disordered structure appeared at saturation stage, where the primary wave stopped growing but the secondary waves generated by nonlinear saturation process continued to develop and resulted in a disordered structure.In order to analyze the detailed development of secondary waves, we performed a 2D Fourier analysis of density variation within a time interval of 10 τ 1 ≈ 0.73 s (τ 1 is the period of first wave mode k = 1).It can be seen from Fig. 4, which was the gray scale maps of 2D mode powers at different time stages, that the energy of horizontal primary wave k = 9 would gradually convey to other secondary waves and finally the maximum energy was shifted toward the longest waves allowed by the simulation box size.The result is similar to the results of large-scale simulations conducted by Oppenheim et al. (2008).We also perform the analysis of phase velocity spectra of 2D modes.Considering that nighttime conditions mostly develop rightward waves traveling in the first quadrant (k > 0, m > 0) and the forth quadrant (k > 0, m < 0) as shown in Fig. 4, in Fig. 5 we only show the phase velocity spectra of a few modes in the first quadrant and the forth quadrant obtained from the time series in the time interval  heim et al. 1996;Fern et al. 2001Fern et al. , 2009a)).Considering the corresponding wavelength of 2D modes, + , the phase velocity spectra can be transformed into corresponding frequency spectra.Figure 6 shows the frequency spectra of a several 2D modes.It can be seen that 2D modes with the same horizontal wavelength have the analogous = + , where k x = 2πk /27 and k z = 2πm /108, in unit of the sound speed C s = 356 m s -1 , while the vertical axis shows the power on a linear scale in arbitrary units.The labels across the top give k (in mode number) for that column, while the labels on the left give m (in mode number) for that row.The top number in the upper corner of each plot shows the reduction in that mode's maximum density perturbation with respect to the mode with the largest density perturbation found in the simulation in decibels.The bottom number in the upper corner shows the total power contained in the mode compared to the total power of the mode containing the most power in decibels.
spectral peaks close to their mean frequencies or threshold frequency derived from their horizontal wavelength and a fixed threshold speed Cs.In other words, frequency spectral peaks vary mainly with the horizontal wave number (k) and are almost independent of the vertical wave number (m); therefore, their phase velocity spectral peaks with the same k vary roughly with m.However, the low frequency characteristic of large-scale 2D modes with small k can cause an apparent shift to a lower speed at large m so the lower speed feature of large-scale 2D modes with a larger m will appear in phase velocity spectra.Actually, according to Fig. 3 for density variation at different time stages and Fig. 4 for the time evolutions of 2D mode power, it can be shown that the evolutions of 2D modes originate from horizontal waves (m = 0) and the spectral contents of 2D modes cover the spectrum of horizontal waves.The horizontal waves with higher power from strong FB instability can dominate power spectra of 2D modes, so the frequency spectral peaks of 2D modes are similar to horizontal waves.Meanwhile owing to FB instability strongly driven horizontal waves, Fig. 5 shows different horizontal wavelength waves and have the analogous peaks of phase velocity spectra which satisfy the characteristic of pure FB waves.
In addition, we introduced a 0.1% perturbation with multiple horizontal wave modes k = 1 ~ 10 under the same background condition as shown in Table 1.This perturbation type is called an initial-multiple waves-perturbation.The time evolution of standard deviation of density given in Fig. 7 showed an analogous saturation development of perturbation as in Fig. 2. From the gray scale maps of 2D spatial mode powers at different time stages (see Fig. 8) the maximum energy shifted to shorter wavelength modes, then gradually shifted toward longer wavelength modes at saturated stage.This is understandable; a wave with a shorter wavelength has a higher growth rate; thus it will reach saturation earlier than a wave with longer wavelength waves.
Owing to the development of 2D mode powers at a saturated stage is similar to previous case; phase velocity spectra (not shown) have feature similar to Fig. 5.

Phase Velocity Analysis of Wave Motion Along Radar Beams
Using the analytical method demonstrated in section 2.2, we carried out a statistical analysis of projection phase velocity of 2D modes along oblique radar beams at different elevation angles.Firstly, for the simulation cases of section  3.1, we can obtain the results of phase velocity analysis of different wavelength waves at an elevation angle of 0° ~ 30° for the time interval t = 2.18 ~ 4.37 s (corresponding to 30 τ 1 ) and statistical mode range 0.8λ ≤ λ ≤ 1.2λ, as shown in Fig. 9a, where case A represents the simulation case of initial-single wave-perturbation and case B for the simulation case of initial-multiple waves-perturbation.It can be seen that these two cases have similar variation with elevation angle; the wavelength dependence of phase velocity at a larger radar elevation angle was remarkable and has shorter wavelength modes tended to be more isotropic.Since the simulation box size and space resolution will affect the evolutions of 2D modes, our statistical analysis of a projection phase velocity of 2D modes can only be performed at a limited elevation angle range.However, it is interesting to note that the saturation pure FB waves can also form an analogous gradient effect on phase velocities of different wavelength waves.
In order to increase the elevation angle range of phase velocity analysis, a larger box size and the optimum space resolution need to be considered.Thus, simulated space size and grid points are redesigned as X × Z = 54 m × 108 m and Nx × Nz = 241 Pts.× 241 Pts.corresponding to a space resolution dx × dz = 0.225 m × 0.45 m.For the convenience of comparison, the simulated time of every step (corresponding to time resolution) is invariable, the other background conditions follow Table1, and a 0.1% perturbation with a single wave perturbation corresponding to a horizontal wavelength  of 3 m was introduced.Figure 9b shows the phase velocity analysis of different wavelength waves at an elevation angle 0° ~ 55° for the statistical mode range 0.8λ ≤ λ ≤ 1.2λ, where case C with red symbols for the time interval t = 2.18 ~ 4.37 s (corresponding to 15 1 xl, 1 xl= 2τ 1 derived from horizontal space size which is twice as large as the cases of section 3.1 for the invariant time resolution) is the same as previous time interval, and case D with green symbols for the time interval t = 2.18 ~ 8.74 s is thirty times of the period of the first horizontal wave mode at k = 1 (i.e., 30 1 xl).In other words, for longer wavelength waves as k = 1, case D with the longer time interval is able to analyze more frequency modes that approach the sum of frequency modes in previous case A. We can see, the analysis results for both time intervals were very similar and displayed an analogous statistical feature as in previous cases where the remarkable wavelength dependence of phase velocity can extend to an elevation angle of 55°.Owing to the larger horizontal space size will favor the evolution of longer wavelength modes near the horizontal direction, the analysis of a small elevation angle like as 5° seems to be similar to the characteristic of pure FB waves, namely, different wavelength waves have the analogous phase velocity.The direction of a smaller elevation angle mainly depends on 2D modes near horizontal direction driven by FB instability, so the analysis results can display the characteristic of pure FB waves.
On the other hand, referring to the dispersion equation associated with a linear instability theory (Hanuise and Crochet 1981) given by where X and i X are the electron-neutral collision frequency, ion-neutral collision frequency, electron gyro-frequency and ion gyro-frequency respectively, k is the plasma wave vector, C s is the ion acoustic speed, V e is the electron drift velocity, and L N is the scale length of electron density gradient.
The threshold condition for plasma waves is obtained for 0 c = as shown in Eq. ( 13), The solution of ( 13) at φ = 0 was investigated by Farley and Fejer (1975) as a function of electron density gradient L N .If e k V $ z = ^h corresponds to a radar elevation angle which the electron drift velocity V e is along horizontal ground, the solution of (13) leads to the theoretical phase velocities of different scale waves at different gradient length L N = ^h for some elevation angle φ, where the gradient length L N = ^h is perpendicular to horizontal ground.Here, considering our simulation parameters of Table 1 and fixed elevation angle, then the phase velocity ratio of two different wavelength waves for the different gradient length L N = ^h can be evaluated.Figure 10 shows the theory phase velocity ratio of 5 to 10 m waves for different gradient length L N = ^h where solid line for φ = 0° and dotted line for φ = 55° can roughly cover the possible elevation angle range of our phase velocity analysis.Meanwhile, the phase velocity ratio of 5 to 10 m waves at different elevation angles from simulation case C are also shown in Fig. 10.In principle, the larger phase velocity ratio, which corresponds to the stronger gradient effect, will happen in the shorter gradient length L N = ^h .Thus, it is clear that the phase velocity ratio of larger elevation angles (> 20°) from pure FB simulation can be compared with the theoretical values with a shorter gradient length L N = ^h about 3 -6 km.That is to say, our phase velocity analysis can reach the gradient length of remarkable gradient effect on phase velocity from the HF radar observations (Hanuise and Crochet 1981).

A large-Scale Simulation for Saturation FB Waves
From past studies (Pfaff et al. 1985(Pfaff et al. , 1987a, b;, b;Chu and Wang 2002), the scale of electrojet layers or mid-latitude irregularities for generating type 1 radar echoes can reach a kilometer in order.Oppenheim et al. (2008) also performed a large-scale simulation of FB waves and showed nonlinear saturated process can develop large-scale waves with wavelengths approaching to 160 m (the maximum length allowed by the simulation box size).However, they used a super computer to support massive calculations for a large-scale simulation with a kinetic particle-in-cell code.Here, our simulation model based on fluid code is also available for large-scale simulations.In addition, the demanded computer resource is small.We considered a simulation box size with X × Z = 360 m × 600 m and grid points Nx × Nz = 1201 Pts.× 1201 Pts., and used only a PC computer for running the numerical simulation.The relevant background conditions still follow Table1, and introduced a 0.1% perturbation with a single wave corresponding to a horizontal wavelength 4 m (i.e., k = 90).Figure 11 shows the gray scale maps of 2D spatial mode powers at different time stages of saturation, where the evolutions of 2D spatial mode power are similar to previous cases, i.e., the maximum energy shifts toward the longest waves allowed by the simulation box size and the larger power focus on the vicinity of horizontal modes.Figure 12 shows a phase velocity analysis of different wavelength waves at different elevation angles for the time interval t = 9.70 ~ 19.41 s where red symbols represent the statistical results of 2D modes with a wavelength range of 0.8λ < λ < 1.2λ and green symbols for 2D modes with a wavelength range of 0.9λ < λ < 1.1λ.Evidently, the data for the two statistical ranges of different wavelength modes are very similar.It means that the statistical range of different wavelength modes for a radar wavelength λ can be reduced.However, the similar statistical results are attributed to that the large-scale simulations with more grid points can generate enough 2D spatial modes with a projected wavelength matching the reduced statistical wavelength regime.In other words, if the sum of 2D spatial mode is sufficient, it can be expected that the analogous data can also be obtained from the statistical analysis of 2D modes with the radar wavelength λ.In addition, it is again evident that the remarkable wavelength dependence of phase velocity appears at a larger elevation angle, the feature of pure FB waves with analogous phase velocity displays at a smaller elevation angle and shorter wavelength waves such as 5 m wave approach isotropic speed.Moreover, large-scale simulation can favor the evolution of much longer wavelength wave such as 40 m wave so the phase velocity analysis can roughly compare with the observation results of multiple frequencies HF radar (Hanuise and Crochet 1981) which revealed the phase velocities of 10.8, 20.3, and 39.3 m wavelength waves are respectively about 295 -345, 215 -275, and 80 -185 m s -1 .As seen in Fig. 12, phase velocities at larger elevation angles derived from the analysis of 2D simulations of pure FB waves did roughly approach the observed phase velocities of different wavelength waves.

SuMMARy ANd dIScuSSION
Via a detailed analysis of 2D spatial modes for longterm saturation simulations of pure FB waves, one sees that a maximum energy will shift to the longest wavelength mode allowed by a simulation space size, the larger power comes from the modes near horizontal direction (i.e E × B direction) ascribes to FB instability, and large-scale 2D modes will tend to have a lower phase velocity smaller than a threshold speed.From the analysis of frequency spectra of 2D modes, the peak power frequencies (or approximately mean frequencies) vary mainly with their horizontal wavelengths that can be attributed to a horizontal mode with the larger power.For the 2D mode of any horizontal wavelength, its phase velocity is smaller if its vertical wavelength is shorter.In general, the simulations with a fixed space size and resolution favor the evolutions of large-scale modes at a larger elevation angle corresponding to the development of large-scale modes with a shorter vertical wavelength.Hence, the shorter vertical wavelength combined with the lower frequency characteristic of large-scale modes will cause the lower phase velocity, so that the past 2D simulations of pure FB waves with fixed space size and space resolution (Oppenheim et al. 1996;Fern andKuo 2001, 2009a) could reveal the analogous features of phase velocity spectra in which spectral peaks of large-scale modes obviously shift to lower phase velocity.
In addition, from the phase velocity of different wavelength waves at different radar elevation angles, one sees that the remarkable wavelength dependence of phase velocity appears at a larger elevation angle and the shorter wavelength waves approach isotropic speed.Meanwhile, different wavelength waves can approach the constant threshold speed at a smaller elevation angle that displays the characteristic of pure FB waves.Actually, the smaller elevation angles that they are smaller than threshold angle of 18.2° derived from cos V V  drift velocity, can satisfy the threshold condition of FB instability and 2D modes traveling in smaller elevation angles mainly driven by FB instability, so the statistical characteristic is similar to pure FB waves.As for the larger elevation angles dissatisfying FB instability, the relevant wavelength dependence of phase velocity is primarily attributed to the features of phase velocity spectra of 2D mode: the spectral peaks of large-scale modes with a shorter vertical wavelength shift to lower phase velocity.Because of largescale modes with shorter vertical wavelength corresponding to longer wavelength waves traveling in larger elevation angle, the statistical results of longer wavelength waves at a larger elevation angle will display the lower phase velocities.In other words, the remarkable wavelength dependence of phase velocity at a larger elevation angle mainly depends on the speed variations of longer wavelength waves which can roughly retain the isotropic speed after reaching a lower speed at larger elevation angles.However, it is interesting to note that the statistical analysis for projection phase velocity of 2D modes along radar beams is compatible with radar observations.From the analysis of previous sections, the wavelength dependence of phase velocity at a larger elevation angle are roughly consistent with the observations of multi-frequency HF radar for density gradient effects on FB waves (Hanuise and Crochet 1981), and the phase velocities of shorter wavelength waves approach an isotropic speed near ion acoustic speed that is consistent with saturation features of small-scale FB waves observed by radars (Cohen and Bowles 1967;Kudeki et al. 1987;Swartz 1997b).Actually, the radar elevation angle means the angle of radar beam to horizontal ground.If the electron drift velocity V e is not in the direction of horizontal ground, the radar elevation angle can not be regarded as e k V $ z = ^h corresponding to the simulation elevation angle.Considering the possible error between radar elevation angles and e k V $ z = ^h, the gradient effect on phase velocity observed by an oblique radar beam can also come possibly from the statistical characteristic of a projection phase velocity of 2D modes along radar beams.Furthermore, we see that the analysis of smaller elevation angles satisfying the threshold condition of FB instability tends to display the feature of pure FB waves.If one to consider the larger electron drift velocity (i.e., a stronger electric field), the elevation angle range satisfying the threshold condition of FB instability will be increased.In other words, the elevation angle range for displaying the feature of pure FB waves will be expanded and the remarkable gradient effect on phase velocity possibly moves to much larger elevation angles where the mode powers are very weak, so the observation for a remarkable gradient effect on phase velocity is rare.In fact, recent observations of two different frequency radars (Tiwari et al. 2003;Haldoupis et al. 2005) also tend to show no remarkable gradient effect on phase velocities of Type-1 radar echoes, where Tiwari et al. (2003) indicated Type-1 velocities are quite close to each other at altitude of the peak electrojet (i.e., strong electrojet), and Haldoupis et al. (2005) argued the gradient effect of type-1 radar echoes originating from high-latitude of strong electrojet.However, our simulation discussion on a stronger electric field (i.e., stronger electrojet) can explain the recent observations.It may be that the study of phase velocity analysis of projected wave motion along oblique radar beams will be helpful in clearing up the argument of gradient effect.
In sum, recent studies (Sahr and Farley 1995;Otani andOppenheim 1998, 2006;Dimant 2000;Fern andKuo 2001, 2009a;Oppenheim et al. 2008) showed that the nonlinear development of different wavelength 2D modes can demonstrate the saturation mechanism of pure FB waves.It is plausible that the statistical analysis of projection phase velocities of different 2D modes along oblique radar beams can cover the saturation characteristics of pure FB waves.This paper has completed the relevant statistical analysis for different wavelength waves at different radar elevation angles and showed that the results can explain radar observations.

APPENdIx A: NuMERIcAl MOdEl OF SIMulA-TION
Our numerical model was described in a previous paper (Fern et al. 2001).The x-axis of the rectangular coordinate system points to the east, the y-axis points to the north and the z-axis points upward.The governing equations consist of the continuity equation of plasma, and the equation for quasi-neutrality condition, ^ĥ h is the current density and n is the number density of the plasma; K B is the Boltzmann constant; v in , v en , i X , e X , T i , T e , B , E l , VD < , ve < , vi < are the ion-neutral collision frequency, the electron-neutral collision frequency, ion gyro-frequency, electron gyrofrequency, ion temperature, electron temperature, magnetic field, perturbation electric field, mean flow electron drift velocity, perturbation electron and ion velocity respectively; Vi < and Ve < can be taken as ; Substituting all these physical quantities into (A4), we obtain (A5): t is a constant magnetic field B 0 of 0.28 G pointing to the north.The ion-neutral collision frequency v in and electron-neutral collision frequency v en are assumed to be constant in our simulation range with v in = 2.5 × 10 3 s -1 and v en = 4.0 × 10 4 s -1 .The ion and electron temperatures are assumed to be 230 K.The background plasma can be regarded as uniform with number density of 1.0 × 10 11 m -3 .These background parameters are similar to that used in the past studies of equatorial electrojet (McDonald et al. 1974).
The numerical computations were performed on twodimensional Cartesian mesh using 121 points in the x direction (east-west) and 121 points in the z direction (vertical).A periodic boundary condition is imposed on both electron density n and electric potential Φ in the x and z directions.Flux-corrected transport (FCT) technique (e.g., Boris et al. 1973;Zalesak 1979) has been applied to carry out the time integration of the continuity Eq. (A1).A detailed discussion of the application of FCT technique to study ionospheric irregularities can be referred to (Chou et al. 1996).At t = 0, the electrons are set to move uniformly at drift speed VD < , and the ions assume a constant velocity , which is the steady state velocity of ions under the combined force from electric field E 0 < and the ion-neutral collision (note that v in >> i X ).Then a density perturbation with amplitude δn = n 0 sin(kx) is superposed on the background density N 0 , where k is the wave-number to be assigned.At each time step of computation, the electron velocity at each grid point is calculated from Eq. (A2), and the ion velocity at each grid is obtained by solving the differential Eq. (A3) using 2 nd order Runge-Kutta scheme.These velocities and densities at each grid are substituted into Eq.(A5) to solve for the electric potential Ф(x, z) using the successiveover-relaxation (SOR) technique.Then, the plasma density distribution n(x, z) at time t + δt is calculated by the FCT scheme to complete one cycle of the computation.In order to guarantee the numerical accuracy, we set the absolute error limit in the potential solver as small as 10 -4 .The simulation is called to stop whenever the relative error of any grid fails to converge within this error limit within 10000 steps of SOR iteration.This criterion sets a limit for density gradient, beyond which the off-diagonal terms become so much larger than the diagonal terms in Eq. (A5) that the SOR calculation fails to converge.That means the simulation will be stopped when the density gradient becomes so sharp that the small waves with scale-length comparable to the grid size starts to grow.

Fig. 2 .
Fig. 2. The time evolution of the standard deviation of the spatial variation of density perturbation for introducing single wave, n n 0 2 d .

Fig. 5 .
Fig. 5.The power spectra of phase velocity for 2D spatial modes at saturation stage.The horizontal axis shows the phase velocity, ph v k k x z 2 2

Fig. 6 .
Fig. 6.The power spectra of frequency derived from phase velocity spectra of Fig. 5 where the red lines represent the mean frequency.The horizontal axis shows frequency = ph v k k 2 x z 2 2 # r + _ i, in unit f = C s /27 m, while the vertical axis shows the power on a linear scale in arbitrary units.The labels across the top give horizontal mode k (in mode number) for that column, while the labels on the left give vertical mode m (in mode number) for that row.

Fig. 7 .
Fig. 7.The time evolution of the standard deviation of the spatial variation of density perturbation for introducing multiple waves with k = 1 ~ 10.

Fig. 8 .
Fig. 8. Gray scale maps of the time average of 2D spatial mode power over 10 oscillation periods of the maximum mode (1, 0) at different consecutive time stages for the simulation case of multiple waves, but the first map of time = 0.04 s from the analysis of one oscillation period displaying the initial evolution.The time center of each time stage is denoted on the top of each panel.The indication of gray scale bar is the same as Fig. 4.

Fig. 9 .
Fig. 9. (a) Phase velocities of different wavelength waves at different elevation angles for case A with red symbols and case B with green symbols, where crosshair represents 5 m wave, open circle for 10 m wave and open square for 20 m wave.(b) Phase velocity analysis is the same as (a), where red symbols represent case C and green symbols for case D.

Fig. 11 .
Fig. 11.Gray scale maps of the time average of 2D spatial mode power over 3 oscillation periods of the maximum mode (1, 0) at two different consecutive time stages for large-scale simulation.The time center of each time stage is denoted on the top of each panel.The indication of gray scale bar is the same as Fig. 4.

Fig. 10 .
Fig. 10.Theoretical predictions of the phase velocity ratio at threshold of 5 to 10 m waves for different plasma density gradient length, where solid line for φ = 0°, dotted line for φ = 55°.The phase velocity ratio of 5 to 10 m waves at different elevation angles for simulation case C with open circle.
where V e = 460 m s -1 represents the electron drift velocity along the +x direction and

Fig. 12 .
Fig. 12. Phase velocity analysis is the same as Fig. 9 for two different statistical wavelength ranges of 0.8λ ≤ λ ≤ 1.2λ shown as red symbols and 0.9λ ≤ λ ≤ 1.1λ as green symbols, where crosshair represents 5 m wave, open circle for 10 m wave, open square for 20 m wave and close square for 40 m wave.
x, z) is the potential function for the perturbed electric field E l , i.e.,