Characteristic properties of Nu whistlers as inferred from observations and numerical modelling

The properties of Nu whistlers are discussed in the light of observations by the MAGION 5 satellite, and of numerically simulated spectrograms of lightning-induced VLF emissions. The method of simulation is described in full. With the information from this numerical modelling, we distinguish the characteristics of the spectrograms that depend on the site of the lightning strokes from those that are determined mainly by the position of the satellite. Also, we identify the region in the magnetosphere where Nu whistlers are observed most often, and the geomagnetic conditions favouring their appearance. The relation between magnetospherically reflected (MR) whistlers and Nu whistlers is demonstrated by the gradual transformation of MR whistlers into Nu whistlers as the satellite moves from the high-altitude equatorial region to lower altitudes and higher latitudes. The magnetospheric reflection of nonducted whistler-mode waves, which is of decisive importance in the formation of Nu whistlers, is discussed in detail.


Introduction
Of the several natural sources of VLF waves in the magnetosphere, lightning strokes are the most familiar.According to the commonly accepted notion, a lightning stroke emits electromagnetic waves into the Earth-ionosphere waveguide.While propagating in this waveguide, some of them leak through its upper boundary and penetrate into the magnetosphere, giving rise to whistlers observed in the opposite hemisphere.The investigation of ionospheric and magnetospheric wave phenomena related to lightning strokes began from classical research by Eckersley (1935) and Storey (1953a), among others, and it has continued up to Correspondence to: D. R. Shklyar (david@izmiran.rssi.ru) the present.The first and most profound summary of the research in this field was given in a book by Helliwell (1965), which was a superlative contribution to whistler studies.
The investigation of nonducted whistler-mode waves in the magnetosphere, in particular of MR whistlers and Nu whistlers, which are the subjects of this paper, also has a long history.We will mention only some work that is directly related to − or especially important for − the present study.An unexpected possibility for whistler-wave reflection when the ions are taken into account in the dispersion relation, and the visualisation of this effect by ray tracing, were first demonstrated by Kimura (1966).In a sense, this finding predicted magnetospherically reflected (MR) whistlers, which were found in the spectrograms of wave data from OGO 1 and 3 (Smith and Angerami, 1968).In their study, mainly devoted to MR whistlers, Smith and Angerami also pointed out that the spectrogram of an MR whistler observed far from the equator may have the shape of the Greek letter ν.They called this phenomenon Nu whistler and suggested its basic mechanism.According to these authors, the minimum frequency on a Nu-whistler spectrogram, where the two branches merge, corresponds to the wave that undergoes magnetospheric reflection at the observation point (see also Edgar, 1976).Magnetospheric reflection occurs when the waves reach some point where their frequency is less than the local lower-hybrid-resonance (LHR) frequency, so it is also referred to as LHR reflection.
Although the close relation between MR whistlers and Nu whistlers was established in the initial work by Smith and Angerami (1968), in later work MR whistlers were studied much more than Nu whistlers; see, for instance, Sonwalkar and Inan (1989), Draganov et al. (1992), Thorne and Horne (1994), and Jasna et al. (1990).Some comments on these studies may be found in the paper by Shklyar and Jiȓícek (2000), where an analysis of MR whistlers observed by MAGION 4 and 5 was supplemented by an extensive numerical simulation of MR-whistler spectrograms.Since then, several authors have used numerical simulation of spectrograms in their studies of MR whistlers.Lundin and D. R. Shklyar et al.: Characteristic properties of Nu whistlers Krafft (2001) demonstrated the similarity of MR-whistler spectrograms that appears in a certain range of L-shells and latitudes when the frequency scale of the spectrogram is normalized with respect to the equatorial electron cyclotron frequency on the L-shell of observation.Jiȓícek et al. (2001) investigated the influence of the plasmapause on MR-whistler spectrograms, concluding that the presence of a pronounced plasmapause renders the traces on the spectrograms indistinct, so the "classical" pattern should be observed only under quiet geomagnetic conditions.An essential contribution to the numerical modelling of MR-whistler spectrograms was made by Bortnik et al. (2003), who included wave intensity in the frequency-time plots, thus making them more like real spectrograms.
A further step in the numerical modelling of spectrograms was taken by Chum et al. (2003), who showed that numerical simulations can be used to model spectrograms not only on a short time scale of the order of 10 s, the so-called detailed spectrograms, but also to model overview spectrograms of data taken along a satellite path during tens of minutes, provided that lightning-induced whistlers are the main emission in the region traversed by the satellite.In this case, whistler emission, trapped in the magnetosphere by LHR reflection, evolves into oblique noise bands above the local LHR frequency; these are qualitatively reproduced by numerical simulations of overview spectrograms.We should mention that LHR reflection also plays an important role for several other types of VLF emission in the magnetosphere.Besides MR whistlers, Nu whistlers, and the LHR noise bands, where LHR reflection is the governing factor, it is also important for chorus waves, as was pointed out recently by Parrot et al. (2003) from their analysis of CLUSTER data.
In this paper we concentrate on Nu whistlers, and proceed as follows.Section 2 is devoted to an analytical description of nonducted whistler-wave propagation, with attention focussed on wave reflection at or well below the LHR frequency.The key points on how ray tracing in the framework of geometrical optics can reproduce the spectrograms of the observed electromagnetic field are discussed in Sect.3. Section 4 presents experimental data on Nu whistlers from the MAGION 5 satellite and compares them with computer simulations.Using the information that may be apparent on the modelled spectrograms, but cannot be seen on real ones, the main properties of Nu whistlers are explained.Our findings and conclusions are summarised in Sect. 5.

Some features of nonducted whistler-wave propagation in the magnetosphere
In this section, we discuss some aspects of nonducted whistler-wave propagation in the plasmasphere that are essential for understanding the phenomena discussed in this paper.

Dispersion relation and group velocity
The equations of geometrical optics, for the ray position r and the wave normal vector k of a wave packet with the frequency ω, can be expressed in Hamiltonian form as where the Hamiltonian H (k, r) is given by the local dispersion relation and v g is the group velocity.
Equations (1) are written above in their general form.We now specify the dispersion relation for whistler-mode waves and the resultant expression for the group velocity, which govern the wave propagation in the approximation of geometrical optics, and which we use in our computer simulations.The dispersion relation, which expresses the wave frequency as a function of wave-normal vector and plasma parameters, can be obtained from the general equation for the wave refractive index in a cold, magnetized plasma (see, e.g.Ginzburg and Rukhadze, 1972).In the frequency band ω ci ω ∼ <ω c , which is that of the whistler mode (ω ci is the ion cyclotron frequency and ω c is the magnitude of the electron cyclotron frequency), and in places where the plasma is dense (ω p ω c , where ω p is the electron plasma frequency) which it is in most of the Earth's plasmasphere, the dispersion relation may be written in the approximate form: where the lower hybrid resonance (LHR) frequency ω LH is given by ; (n e , m e are the electron concentration and mass, respectively, while n α , m α are the same quantities for ions of the species where k and k ⊥ are the components of the wave-normal vector parallel and perpendicular to the ambient magnetic field, respectively, θ =cos −1 (k /k), and where c is the speed of light.From Eq. ( 3) one can see that the characteristic value of the wave number in the whistler frequency band is q≡ω p /c, and that for a given wave-normal angle θ , the dependence of the wave frequency on k involves only the ratio k/q.For the so-called quasilongitudinal waves (Ratcliffe, 1959;Helliwell, 1965) k q, whereas the inequality k q corresponds to quasi-resonance waves (see, for example, Walker, 1976;Alekhin and Shklyar, 1980).Some features of quasi-longitudinal and quasiresonance wave propagation, useful for understanding the results of numerical simulations based on the dispersion relation Eq. ( 3), were discussed by Shklyar and Jiȓícek (2000).
The expressions for the longitudinal (parallel to the ambient magnetic field) and transverse (perpendicular to the ambient magnetic field) components of the group velocity that follow from Eq. ( 3) are: From Eq. ( 6) one can see that v g has the same sign as k , hence both quantities change sign at cosθ =0; obviously, from Eq. ( 3), this can happen only when ω<ω LH , in which case v g =0 for ω 2 =ω 2 LH /(1+q 2 /k 2 ) (cf.Eq. ( 3)).As for v g⊥ , it has the same sign as k ⊥ for k<q, while for k>q it is directed opposite to k ⊥ for ω 2 >ω 2 LH /(1−q 4 /k 4 ) and vice versa.The last statement becomes apparent if we rewrite the expression Eq. ( 7) for v g⊥ , eliminating k 2 with the help of the dispersion relation Eq. (3): Also, when k =0 and v g =0, it follows from Eqs. ( 8) and (3) that v g⊥ is always parallel to k ⊥ , since, under these conditions, ω 2 <ω 2 LH .Concerning v g , it is easy to see that the first term in the expression Eq. ( 6) for v g is always much less than the second term.Indeed, for k 2 /q 2 1, the ratio of the first term to the second is less or of the order of ω 2 LH /ω 2 c 1.For k 2 /q 2 <1, the first term is of the order of k ω 2 LH ω q 2 while the second one is of the order of k ω 2 c ω q 2 k 2 q 2 , thus they become comparable only when However, such small values of k 2 /q 2 are outside the range of validity of the approximate dispersion relation Eq. (3), since they correspond to frequencies of the order of the ion cyclotron frequency.

Magnetospheric reflection of whistler-mode waves
The possibility that whistler waves might be reflected within the magnetosphere was suggested and studied by Kimura (1966).In the one-dimensional case, wave reflection corresponds to a change in sign of the group velocity.In the two-dimensional case the situation is more complicated.If, for example, the longitudinal component of the group velocity v g greatly exceeds the transverse one v g⊥ , then the wave reflection corresponds to the point where v g changes its sign, and thus v g =0.However, in the case where v g ∼v g⊥ , a wave may be reflected with respect to one coordinate but continue propagating in the same direction with respect to the other coordinate.For example, in a dipolar magnetic field, the wave reflection with respect to the height (or the radial distance r from the Earth's center) takes place when where λ is the geomagnetic latitude; when condition ( 9) is satisfied, the radial component of the group velocity dr/dt vanishes.
Reflection with respect to both coordinates would require v g =v g⊥ =0 which is impossible for whistlers.Indeed, from Eqs. (3) and ( 6) it follows that v g =0 for while v g⊥ =0 implies (see Eq. ( 8)) Thus, strictly speaking, the reflection of a whistler wave can never happen: magnetospheric reflection is in fact a reversal of the group velocity in a small region of space as the result of refraction (Kimura, 1966).
To find the conditions for magnetospheric reflection, we first describe this phenomenon more rigorously, contrasting it with regular refraction.First of all, we shall speak of reflection as being a property of a ray, regardless of time.If, in a region that is small compared to the characteristic scale of plasma inhomogeneity, a major variation of the direction of the group velocity takes place, whereas before entering and after leaving this region the direction of the group velocity varies relatively slowly, then we shall call this event wave reflection.The foregoing conditions can be expressed as where v 1 and v 2 are the group velocities at the entrance to and at the exit from the reflection region, respectively.The group velocity v g is a function of both k and r, v g =v g (k, r); however, we can neglect the variation of r in a small reflection region.Hence, the amount by which the group velocity varies in passing through the reflection region may be estimated with the help of Eq. (1) as Contours of (dV ⊥ / dk ⊥ )(ω / <V g > V max ) where L is the characteristic scale of the plasma inhomogeneity, l is the length of the part of the ray in the reflection region, <v g > is the average magnitude of the group velocity such that l/<v g > is the duration of the reflection process, and the subscript "max" denotes the maximum value.Using Eq. ( 13) and the notation we rewrite the reflection conditions of Eq. ( 12) in the form: Thus, the wave reflection takes place for those k and k ⊥ , and the corresponding wave frequencies, for which the quantity on the left-hand side in Eq. ( 14) is much larger than unity.This quantity has been calculated numerically using Eq. ( 3), together with the expressions Eqs. ( 6), ( 7) for the group ve-locity.The results are shown in Fig. 1, where we see that the reflection is determined by the parameter which at k →0 greatly exceeds unity, while the other quantities proportional to ∂v g /∂k ⊥ and ∂v g⊥ /∂k ⊥ are ∼ <1 over the whole plane (k ⊥ , k ).We should emphasize that the value of the parameter Eq. ( 15) depends on how we define the size of the reflection region, so it is not determined uniquely.Its only important property is that it is much larger than unity, which ensures that the direction of the group velocity varies rapidly along the ray in the reflection region, compared with its behaviour on other parts of the ray.
As we have seen above, for k =0, the parallel component of the group velocity v g =0, and the wave frequency is determined by Eq. ( 10).From this equation it follows in particular that in the quasi-resonance regime k 2 q 2 , the wave reflection takes place at frequencies close to the LHR frequency ω LH , whereas for k 2 ∼ <q 2 the wave frequency ω may be well below ω LH .
In Fig. 1, the contours of normalized frequency and of the reflection parameters are shown on the (k ⊥ , k )-plane.Although the wave frequency remains constant when the wave propagates in a stationary inhomogeneous medium, this does not mean that it remains on the same contour line of the normalized frequency, since the normalizing LHR frequency may change.Obviously, instead of (k ⊥ , k ), two other quantities may be chosen as the independent variables determining the wave characteristics.In particular, it is convenient to analyze the features of the wave reflection and possible types of ray trajectories in the reflection region with the help of a diagram on the (k/q, ω/ω LH )-plane, as shown in Fig. 2.
In this analysis, we will assume that k ⊥ >0, i.e. that the wave-normal vector is directed towards higher L-shells.The solid line in the figure is determined by Eq. ( 10) and corresponds to k =v g =0.According to Eq. ( 3), the same line defines the minimum possible wave frequency as the function of k/q, so the dispersion relation has no roots below this line.As we have seen above (cf.Fig. 1.), large values of the reflection parameter Eq. ( 15), typical of wave reflection, are attained in the vicinity of k =v g =0.Thus, on the diagram in Fig. 2, the reflection region is represented by the narrow region above the solid line.The dotted line is determined by Eq. ( 11) and corresponds to v g⊥ =0.In the region to the left of this line, v g⊥ has the same sign as k ⊥ (positive under our assumption), while in the region to the right of this line it has the opposite sign.Clearly, when the wave approaches the reflection region, it always moves from higher towards lower values of ω/ω LH on the diagram in Fig. 1.Another important point is that in the reflection region v g⊥ is always positive, while v g changes its sign.Thus, if before and after reflection the wave remains in the shaded region to the left of the dotted line, which is typical of k/q ∼ <1, then v g⊥ does not change sign, and the ray has the shape of an arc.On the contrary, if before and after reflection the wave remains to the right of the dotted line, which is typical of k/q 1, then v g⊥ changes sign, and the ray has the shape of a loop.These features of the wave reflection are illustrated by Figs. 3 and  4. In Fig. 4, L 3 (rather than the more natural quantity L) is chosen as one of the coordinates in order to make the loops in the ray more obvious.Figures 3 and 4 correspond to a 5-KHz wave starting vertically at 15 • geomagnetic latitude, at the height 500 km; the plasmasphere is smooth and the propagation time is set to 3 s.

Spectrogram modelling by means of ray tracing
VLF data from Magion 5 will be presented below in the form of spectrograms.These were made with a sampling frequency of 44 100 Hz and an integration time of 23.22 ms; thus, each instantaneous spectrum was evaluated from 1024 data points.The corresponding resolution in frequency is ∼ <100 Hz.Each spectrogram comprises about 300 instantaneous spectra and covers a time interval of 7 s.It is a representation of spectral intensity in the frequency-time plane, with time along the x-axis, frequency along the y-axis, and the intensity indicated by the degree of darkness on blackand-white spectrograms, or by the use of colour.If the spectral intensity is appreciable only along some curves in the (f, t)-plane, as is the case for MR and Nu whistlers, the problem of spectrogram modelling consists of two parts: firstly, constructing the frequency-time plot, which may have many branches, of course; and secondly, attributing the corresponding intensity to each curve.Here we discuss how it is done by means of ray-tracing calculations based on the equations of geometrical optics.Since very many rays must be calculated in order to reproduce the main features of Nuwhistlers on a model spectrogram, we use relatively simple models for the geomagnetic field and for the distributions of plasma density and LHR frequency, all given by analytical expressions (see Shklyar and Jiȓícek (2000) for details).

Constructing the frequency-time plot
We assume that a thin layer in the upper ionosphere is illuminated by waves from a lightning stroke, and that this process is effectively instantaneous on the time scale of wave propagation to the satellite.We also assume that initially all waves have their normal vectors directed vertically, due to refraction by the ionosphere.Similar assumptions have been used in all of the work on spectrogram modelling cited above.Since the vertical dimension of the illuminated layer (which, in turn, plays the role of an illuminating region for the magnetosphere) is much smaller than its dimension in the horizontal plane, and since the vertical direction of the wavenormal vectors implies that the waves propagate in the meridian plane, computation of the rays is now a two-dimensional problem, with initial conditions given on some line that approximates the thin layer.As such a line, we take a part of the arc at the height of 500 km above the Earth's surface.For numerical modelling of spectrograms, Storey suggested using the notion of a group front.Here we reproduce, with his permission, his definition and physical explanation of this notion.
For any particular frequency, consider all the possible rays that can be traced upwards from the illuminating region, with initial conditions as defined above.Imagine that along each ray, starting at the instant of the lightning stroke, a point moves away from the illuminating region at the local group velocity.Then, at any later instant, the set of all such points defines a surface: this is the group front for the frequency concerned.
A more physical way of visualizing the group front is to imagine that the lightning stroke emits a narrow-band impulse instead of a wide-band one, thus giving rise to a quasi-monochromatic disturbance that propagates through the magnetosphere in the form of a thin sheet.The group front is the surface at the center (i.e.half-way through) of this sheet.Within any given ray tube, the disturbance is a wave packet moving along the tube at the group velocity, and within this packet, the point of maximum amplitude lies on the group front (Storey, 2003, private communication).
In the case under discussion, initially, the group fronts for all frequencies coincide with the illuminating region, i.e. the part of the arc extending over a range of latitudes at 500 km height.With increasing time after the lightning stroke, the group fronts separate due to the different group velocities of waves with different frequencies, while every group front is deformed due to plasma inhomogeneity, and also due to the different initial conditions for the waves of the same frequency starting vertically at different latitudes.
To plot a point in the (f, t)-plane of a spectrogram of data from a satellite, we should find the time at which the group front crosses the satellite position.This procedure can be formalized as follows.The equations of geometrical optics in their general form have been written in the previous section (see Eqs. (1), ( 2)).As is well known (see, for example, Landau and Lifshitz, 1976), when the Hamiltonian H does not depend on time, it is a constant of the motion.Thus, according to Eq. (2), Eq. ( 1) describe a wave packet with constant frequency.We should emphasise that in the 2-D case the wave frequency alone does not determine the wave packet uniquely.
To solve Eq. ( 1), it is most convenient to use canonically conjugate variables.However, once the solution has been found, it can be expressed in terms of any variables that are uniquely related to the canonical ones.The general solution of the Eq. ( 1) has the form while ω(k, r)=ω(k 0 , r 0 ).In the 2-D case considered, both k and r are two-dimensional vectors.Moreover, since we start all rays from a single altitude with the wave normals vertical, there are in fact only two independent initial variables, and we may choose the wave frequency to be one of them.As the second initial variable we choose the initial geomagnetic latitude λ 0 , as is usual in computer simulations of this kind.Then, taking the McIlwain parameter L and the geomagnetic latitude λ as two coordinates, we can rewrite the solution Eq. ( 16) in the form When the solution in the form Eq. ( 17) is known, all of the local wave characteristics, such as the group velocity, the refractive index and the wave-normal angle may be found; the wave frequency is constant along the ray, at the value chosen initially.
The relations Eq. ( 17), being the solution of the equations of motion, are unique functions of their independent variables.The first two relations in Eq. ( 17), which define in a parametric way the time t and the initial latitude λ 0 as functions of ω, λ, and L, can, in principle, be solved for t and λ 0 : These functions, however, may have many branches, that is to say, they may be multi-valued.As we shall see below, the different branches of the solution Eq. ( 18) correspond to different numbers of hops across the equator that the wave packets perform in the magnetosphere.The first function in Eq. ( 18) defines the time when the group front for the frequency ω crosses the satellite position λ, L; thus, it yields the time-frequency curves on the spectrogram.The second function determines the initial latitude for the frequency ω on each branch.This latitude can easily be displayed on a model spectrogram, which, of course, is impossible for real ones.Moreover, since t and λ 0 are now functions of ω, λ, and L, the same is true for k.Thus, all the characteristics of the wave packets that contribute to the spectrogram at a given satellite position become functions of ω and of the hop number, and can be displayed if desired.

Calculation of spectral intensity
A thorough discussion of the rigorous ways of displaying intensity on spectrograms would lead us too far away from the main topic of the present paper: it will be presented elsewhere.Here we discuss only the main aspects of this problem.We regard the wave field as the sum of a set of wave packets propagating in the magnetosphere with their group velocities.The central frequency of each wave packet is conserved, while its wave-normal vector varies along the ray, satisfying a local dispersion relation at each point.In a sense, the wave packet is determined by a bunch of close trajectories in the phase space (k, r) whose projection onto the coordinate space represents the ray tube.The ray itself, and the variation of the wave-normal vector along it, are described by the equations of geometrical optics Eq. (1).
As is well known (see, for instance, Fermi, 1968), the wave packet in geometrical optics is an analog of the mass point in mechanics; as such, it is characterized by the initial coordinates of its amplitude maximum and its wave vector at this point.However, in contrast to a mass point, a wave packet is of finite size; it is also characterized by its width in k-space and the corresponding reciprocal dimension in coordinate space.Thus, in the general 2-D case, a wave packet is characterized by four parameters.However, in the case under consideration, when all rays start vertically at 500 km altitude, only two parameters are needed to characterize a wave packet, and these can be chosen, as above, to be the wave frequency ω and the initial latitude λ 0 .Nevertheless, there are many wave packets with the same frequency, starting at different initial latitudes λ 0 .The receiver on the satellite measures the total field, and even after spectral analysis, it is possible that more than one wave packet may contribute to this field.(We remind the reader that, although different rays never intersect in phase space, their intersection in coordinate space is not forbidden by the equations of geometrical optics.)The key point that simplifies spectrogram modelling in our case is that, as the ray tracing shows, the rays with the same frequency that start at different latitudes never intersect in coordinate space.This means that the satellite, which may be considered as a fixed point, never receives more than one wave packet at any time.Thus, in calculating the spectral intensity at a given frequency, we need to consider only one wave packet, provided that the duration of the time interval over which the spectrum is evaluated is much less than the typical bounce period of the wave packets in the magnetosphere, which we have always found to be the case.(Obviously, different wave packets may come to the satellite on different hops.) Let t − and t + be the times at which the group fronts for the frequencies ω−δω and ω+δω, respectively, cross the satellite position.We can then state that the wave field received by the satellite during the time interval |t + −t − | is that of some wave packet with central frequency ω and bandwidth ω=2δω.If the time of spectral evaluation t is less than |t + −t − |, then the spectral intensity in the frequency band (ω−δω, ω+δω) will be nonzero over the whole interval between t − and t + ; in the opposite case, the spectral intensity will be nonzero throughout some interval of duration t that includes (t + +t − )/2.Here we assume that the frequency resolution is ∼ ω, and thus ω t ∼ > 2π.In calculating the spectral intensity for display on a spectrogram, we need to take into account the variation of the wave-packet amplitude along the ray caused by geometrical factors.To do this, we proceed as follows.Consider, for the frequency ω, the ray that passes through the position of the satellite.Let s be the distance along this ray to any point on it, and let E ω (t, s) be the wave-field component (the waveform) in the frequency band (ω−δω, ω+δω) that the satellite would measure if it were at this point.As has been argued above, for a given s and hop number, this field belongs to one particular wave packet, characterized by ω and λ 0 .The energy-conservation law for this wave packet has the form: This and later equations concern the particular wave packet characterized by the two parameters ω and λ 0 .Henceforth, however, the second parameter will be omitted for shortness.Equation ( 19) can be rewritten as where s is the coordinate along the ray considered, σ is the cross section of a thin ray tube centered on this ray and v g =(v 2 g⊥ +v 2 g ) 1/2 .Obviously, the wave packet passes over the point s during the time interval from t 1 (s) to t 2 (s), where We then integrate Eq. ( 20) over t from t 1 (s) to t 2 (s).Since the integrand tends to zero at both limits of integration, the contribution from the first term vanishes.For the same reason, the integral over t can be shifted into the argument of the derivative with respect to s.As a result we obtain Thus, the quantity U ω (t , s)dt is conserved along the ray.This quantity is the total energy of the given wave packet.
The wave energy density U is related to the electric field E as follows: where E ω (t, s) is the electric-field component of the wave packet measured by the satellite, while the factor w(s, ω, θ ) depends on frequency, wave-normal angle, and the local plasma parameters, and also on the wave mode, of course.At this point we assume that the satellite measures the component of the electric field perpendicular to the Earth's magnetic field B 0 in the (k, B 0 )-plane.Then, in the same range of parameters where the dispersion relation Eq. ( 3) is valid, the expression for w has the form where The quantity w may be regarded as a function of s, ω, and θ since in a cold magnetoplasma the refractive index N is a function of ω and θ .With this notation, the conserved quantity takes the form On the other hand, from the well-known theorem in spectral analysis that relates the field component E ω (t , s) of a received wave packet to its time-dependent spectral amplitude E(ω, s, t), we have If t 2 −t 1 < t, then the integrals in Eqs. ( 25) and ( 26) are equal, while in the opposite case the integrals are proportional to the intervals of integration.Taking these facts into account, we obtain from Eqs. ( 25) and ( 26): We see that, apart from the quantities directly determined by the equations of geometrical optics, an additional quantity that needs to be calculated is the cross section of the ray tube.We will assume that there are no gradients in the azimuthal direction, so the waves propagate in meridional planes.Then the width of the ray tube in the azimuthal direction is where x is the Cartesian coordinate in the meridional plane orthogonal to the dipolar axis of the Earth's magnetic field, L and λ are, as before, the McIlwain parameter and the magnetic latitude, respectively, and is the range of azimuthal angles over which the wave packet extends, which is a constant of its motion.Thus, a non-trivial part of the cross section is its width ξ in the meridional plane, which defines the ray-tube cross section σ according to To find the quantity ξ , let us consider two neighbouring rays.Let (x 1 , z 1 ) and (x 2 , z 2 ) be two neighbouring points on the first and second ray, respectively, and let ψ be the angle between the group velocity and the x-axis.Then the width of the ray tube in the meridional plane, ξ , is Since the vector (−sinψ 2 , cosψ 2 ) is orthogonal to v g and, thus, to the ray, this result does not depend on the particular choice of the point (x 2 , z 2 ), provided that it is close enough to the point (x 1 , z 1 ) where the cross section is calculated.According to Eq. ( 27), the spectral intensity |E(ω, s)| 2 at the observation point s is determined by the factor (σ v g w) s and the conserved value W of the wave packet energy.Thus, to include the spectral intensity in model spectrograms, we need to supplement the ray-tracing calculation with the evaluation, for each ray, of the quantities v g and w determined by Eqs. ( 6), (7), and ( 23), respectively, and with the calculation of a neighbouring ray, which enables us to find ξ Eq. ( 30) and thus the cross section of the ray tube Eq. ( 29).The corresponding data base, which is similar to the one described by Shklyar and Jiȓícek (2000) but supplemented with the relative-intensity parameters, has been computerized.Spectrograms calculated with the help of this data base are presented in the next section.
3.3 Comparison with the approach of Bortnik et al. (2003).
As was mentioned in the Introduction, Bortnik et al. (2003) made an important step in numerical modelling of MR whistlers by including spectral intensity into simulated spectrograms.Like the spectrograms simulated by Bortnik et al. (2003), ours now display spectral intensity.However, the method we use differs from that of Bortnik et al. (2003) in several respects.Firstly, we deal from the outset with wave packets of finite spectral width f , corresponding to the frequency resolution on real spectrograms.In this case, the time interval during which the frequency band f − f/2, f + f/2 is received on the satellite is determined by the group-front crossings of the satellite position, as suggested by Storey ((2003), private communication).This time is determined unambiguously, with no uncertainty; it does not use the notion of detection area, the extent of which is difficult to define consistently due to the continuous merging of different rays at the same frequency.Secondly, as has been shown by Storey (1953b), when a dispersed signal is passed through a bank of narrow band filters, the temporal variation of its instantaneous frequency is measured most accurately when the bandwidth of the filter equals , where f i is the instantaneous frequency.As the wave phenomena that we model are characterized by a rate of frequency variation of the order of a few kHz per second, the filter bandwidth should be of the order of 50 Hz for the sharpest output, so we choose this value as the frequency step in our calculations.Thus, we consider that the interpolation procedure used by Bortnik et al. (2003), which yields a frequency resolution of ∼1 Hz, is superfluous in this respect, all the more so because, finally, they set the width of their frequency bin to 50 Hz.And thirdly, another difference between their approach and ours lies in the way we evaluate the spectral intensity: instead of computing millions of interpolated rays, each weighted with a measure of wave energy, and then calculating the energy carried by those rays that cross the detection area, we calculate the variation of the ray-tube cross-section, then use energy conservation and Parseval's relation to translate the energy in each wave packet, of bandwidth 50 Hz, into spectral intensity displayed on a spectrogram.As for initial distribution of the wave energy among wave packets, we use the following model.We assume that each wave packet is determined by its frequency f and initial latitude λ 0 , and that all wave packets have the same frequency width f and occupy the same spatial width λ 0 at the beginning.Since initially all wave packets have vertical direction of their wave normal vectors, and negligible dimension in radial direction, these parameters determine the wave packet uniquely.The total energy of each wave packet, which, of course, is conserved, is modelled as where ϕ(λ 0 ) is a smooth function, which decreases with the distance from the center of illuminating region; and η(f ), which describes the frequency dependence of the wave energy distribution, is adopted from the paper by Lauben et al. (2001): where f is the wave frequency in kHz.The same frequency dependence has been used by Bortnik et al. (2003).On the other hand, we do not take into account wave growth or damping, so in our case the variation of energy density along the ray is due only to geometrical factors.

Nu whistlers from the MAGION 5 satellite and their modelling
As was noticed in the earliest satellite experiments, VLF data from a satellite exhibit a much richer variety of wave phenomena than data from ground-based observations.The reason is that on satellites, phenomena related to the quasiresonance (or nonducted) type of whistler-wave propagation are observed, as well as those related to quasi-longitudinal (ducted) propagation, while ground-based data are mainly limited to the latter.
In this section, we present examples of MR whistlers and Nu whistlers observed on board the MAGION 5 satellite.The data are available from June 1998 to July 2001.Since they were transmitted in analogue form to the ground station in Panska Ves (50.53 • N, 14.57 • E) in real time, it was the radio visibility of the satellite that limited the parts of the orbits from which data could be obtained.A graphical illustration of the parts on which VLF data were recorded is given in Fig. 5; the smaller parts on which MR and Nu whistlers were observed are marked by asterisks.One can see from this figure that MR whistlers, i.e. the waves that have undergone magnetospheric reflection, are mainly observed on L-shells from L∼1.8 to L∼3.The lines of constant altitude and magnetic latitude are shown along with the magnetic field lines (constant L-shells) for convenience.Usually the range of radio-visible longitudes from 10 • W to 70 • E was covered.We note that MAGION 4, and also MAGION 5 on the descending parts of its orbits, observed MR whistlers in the equatorial region at altitudes from about 1.3 to 2 Earth radii, which is far from the regions where the waves are reflected.Such MR whistlers have been discussed in detail by Shklyar and Jiȓícek (2000); one example is presented in Fig. 6, with its simulated counterpart shown in Fig. 7. On the contrary, over the ascending parts of the MAGION 5 orbits, whistlers could be observed in the regions of their magnetospheric reflection.Before showing examples of spectrograms taken in these regions, which exhibit the ν-shaped patterns characteristic of Nu whistlers, we recall some features of ducted and nonducted whistler wave propagation.As there is no clearcut boundary between these two types of propagation, it is sometimes hard to distinguish between them in satellite data, particularly in the case where the waves propagate from the Earth and are received on a satellite before crossing the equator (fractional-hop whistlers).The degree of dispersion D of the fractional-hop whistlers is very small (∼(10−20) s 1/2 ) due to the short distance of propagation, and it is almost impossible to distinguish between ducted and nonducted propagation in this case; (see, for instance, in the spectrogram of bottom panel in Fig. 8, the single trace at the time ∼3 s).The nature of the propagation in this case can be determined only from the fact that each whistler is followed by Nu whistlers, with almost the same delay in all such events observed during tens of seconds.From time to time, the spectrograms show subsequent traces of reflected ducted whistlers, indicating ducted propagation.When analysing the first magnetospheric reflection on spectrograms, one should also take care to distinguish between Nu whistlers and the traces from double or multiple lightning strokes.For example, the traces shown in Fig. 8, in the second panel from the top, which resemble those of Nu whistlers, are in fact those of normal whistlers originating from multiple lightning strokes in the opposite hemisphere.This can be established from the fact that the traces on the spectrogram all have exactly the same form, while the time delay between successive traces varies randomly in time.
The events shown on the third and fourth panels in Fig. 8 are different.Here we see Nu whistlers in which the first trace corresponds to waves propagating downwards, whereas the second one is formed by waves propagating upwards after MR reflection.Note that the traces are not parallel in this case.These examples show that certain wave phenomena observed on satellites can be identified only by following their evolution and recurrence in the data.
Figure 9 demonstrates how the spectrograms with MR whistler traces change their character along the ascending parts of the MAGION 5 orbits.As the altitude and latitude of the satellite increase (cf.Fig. 5), the time intervals between the traces of successive hops increase also, evidently due to  the lengthening of the ray paths.Note that the initial traces in Fig. 9, the second panel from the top, originate from double lightning strokes.Moreover, the traces of higher-order hops become more nearly horizontal and their range of frequency decreases, while the upper limit of this range approaches the LHR frequency.Simultaneously, the merging frequency at which the traces of the downward and upward propagating waves join one another (i.e. the frequency of the wave reflected at the observation point) increases with the number of hops.This implies that the frequency at the point of the first magnetospheric reflection may be well below the local LHR frequency, which is typical of lower-frequency waves starting at low latitudes.Those waves propagate towards higher L-shells due to the directional properties of their group velocity.This picture is consistent with the reasoning in Sect. 2. We should stress that the key to fitting the simulated spectrogram to the observed one is to choose the illuminating region correctly, since the position of the satellite is precisely known.By this means, computer simulations of the spectrograms may serve as a tool for locating the illuminating region, and thus determining the refractive properties of the ionosphere.
Here we do not discuss the features of MR whistlers, as these are not the subject of the present study (see the literature on MR whistlers cited in the Introduction).We will only mention that pairs of lightning strokes, situated symmetrically with respect to the equator at ground level, will produce MR whistlers with similar spectrograms at the equator in the magnetosphere.The situation is quite different for Nu whistlers since these are observed in the magnetosphere far from the equator.When the lightning stroke is in the same hemisphere as the satellite, the spectrogram starts with the triggering whistler with relatively low dispersion, followed after a significant delay by the Nu whistler, which is identifiable by the characteristic divergence of its two branches from their merging frequency just below the local LHR frequency (see, e.g.Fig. 8).Depending on the position of the satellite with respect to the illuminating region, the triggering whistler may or may not be seen on the spectrogram.
When the source is in the opposite hemisphere with respect to the satellite, the series begins with the Nu whistler, but its two branches diverge less than in the previous case, and their merging frequency is well below the local LHR frequency.Nevertheless, the merging frequency still corresponds to the wave that undergoes magnetospheric reflection at the satellite position; this fact can be proved with the help of simulation, by finding the initial latitude for this frequency and computing the corresponding ray.As was mentioned above, the initial latitudes can be visualized readily on the model spectrogram, and then it is quite easy to find the initial latitude for any point on the Nu-whistler trace.Thus, we find that, in this case, the magnetospheric reflection occurs well below the local LHR frequency, as was discussed in Sect. 2. A series of simulated spectrograms that illustrates the dependence of Nu-whistler shape on the satellite position and on the illuminating region is presented in Fig. 10.
We conclude this section with a discussion of some overview spectrograms from MAGION 5, which represent the VLF spectrum on a time scale of the order of 20 minutes.Since during this time interval the satellite crosses an extended region of the magnetosphere (see the captions below the upper panels of Figs. 8 and 11), the variation of spectral intensity plotted on (f −t)-plane is mainly due to spatial variation of the spectral intensity distribution.If the satellite moves towards higher L-shells and higher latitudes, as it is the case on ascending parts of MAGION 5 orbits (cf.Fig. 5), the pattern of the overview spectrograms may be that of oblique noise bands above the LHR frequency, merging into the LHR noise band when the satellite reaches high altitudes and L-shells (Chum et al., 2003).Examples of the oblique noise bands observed by MAGION 5 are shown on the upper panels of Figs. 8 and 11, the lower panels giving further examples of Nu whistlers.The following considerations may help in understanding the formation of oblique noise bands on overview spectrograms.
We suggest that the oblique noise bands represent lightning-induced VLF emissions that propagate in the magnetosphere in the nonducted mode under quiet magnetospheric conditions.The characteristic features of such propagation are: -transition to the quasi-resonance regime of propagation; -multiple magnetospheric reflections in the regions where the wave frequency is close to the local LHR frequency; and -ray focusing, and merging of the rays that start on different latitudes with the same frequency (see, e.g.Shklyar and Jiȓícek (2000) for details).
We should mention that if the electromagnetic energy radiated by lightning leaks into the magnetosphere at middle latitudes, the resulting whistler mode waves spreading in the nonducted mode attain the quasi-resonance regime of propagation already on the first hop.The L-shell on which waves of a given frequency settle down decreases with an increase of the number of hops, but for a given hop number, the higher the frequency, the lower the corresponding L-shell.Thus, for each hop number there is a rough correspondence between the L-shell and the frequency of the quasi-resonance waves, with higher frequencies corresponding to lower L-shells.Consequently, the satellite, which receives the local spectrum at each instant, observes wave energy in quite narrow frequency ranges, typical of the current L-shell and latitude.These ranges form oblique bands as the satellite moves across L-shells.The center frequency of each band decreases if the satellite moves towards higher Lshells and vice versa, as it is clear from the relationship between frequencies and L-shells mentioned above.We remind the reader that in the quasi-resonance regime, low-frequency waves (ω ω c ) propagate almost along the ambient magnetic field, so the picture described above depends only slowly on latitude.This permits us to use the L-shell as the main parameter governing the spectrum.
In the analysis given above we assumed that on most parts of the trajectories, the waves propagate in the quasiresonance mode, which is true when the lightning activity, and hence the illuminating region, are at middle latitudes.When the waves start at low latitudes, they propagate in the quasi-longitudinal regime during a significant number of hops, especially the lower frequency waves, and the picture described above breaks down.Lightning activity at high latitudes is also unfavorable to the formation of oblique noise bands.In this case the wave energy density decreases so fast due to geometrical factors that these waves contribute very little to the spectrum at the L-shells between 1.8 and 3.7 where the oblique noise bands are observed (for more details, see Chum et al., 2003).In the present paper, we will not go into the details of how we model overview spectrograms, and we limit ourselves to one example of a simulated spectrogram corresponding to MAGION 5 orbit 7102, as shown in Fig. 12.The following lightning statistics have been assumed: lightning strokes appear randomly in time, with a maximum interval of 5 s between strokes; the illuminating  regions from the lightning strokes are in the Northern Hemisphere, at the height of 500 km, and are randomly distributed over latitudes from 15 • to 50 • , with an average width of 15 • .

Summary
Extensive data from MAGION 4 and 5 on VLF phenomena in the plasmasphere have inspired new efforts to understand lightning-related emissions by means of computer simulations.The main focus of the present study is Nu whistlers, which were first reported by Smith and Angerami (1968).We have presented many examples of Nu whistlers from MA-GION 5 measurements, using numerical simulations as a tool for understanding their main features.A short summary of our results now follows.
-MR whistlers, which are observed near the magnetic equator and have spectrograms that are symmetrical with respect to the hemisphere of the source, gradually change into Nu whistlers as a satellite moves towards higher latitudes.The main region where Nu whistlers were observed on MAGION 5 is located between L=2 and L=3, on invariant latitudes from 18 • to 30 • .
-The spectrogram of a Nu whistler is not symmetrical with respect to the source location.If the illuminating region and the satellite are in the same hemisphere, then the spectrogram starts from a single trace of a nonducted sferic, followed, after a certain delay, by a ν-shaped trace with widely diverging branches; the merging frequency in this case is close to the local LHR frequency.
If the illuminating region and the satellite are situated in opposite hemispheres, the spectrogram starts from a ν-shaped trace, usually in a wide frequency band, with narrowly diverging branches; the merging frequency is well below the local LHR frequency.
-The merging frequency of Nu-whistler traces increases with the number of hops and approaches the local LHR frequency, while the traces themselves become flatter and more diffuse.
All these features seen in real spectrograms are readily reproduced and interpreted by simulations.

Fig. 1 .
Fig. 1.Contours of normalized frequency and reflection parameters.Note different colour bars associated with different subplots.Only the parameter shown on the upper right panel greatly exceeds unity and, thus, determines the wave reflection.

Fig. 3 .
Fig. 3. Ray trajectory and wave propagation characteristics for 5-kHz wave.Comparison of the latitude variation and the reflection parameter shown on the upper and middle right panels, respectively, clearly shows that the reflection takes place when the reflection parameter has a peak.

Fig. 4 .
Fig.4.Zoomed-in view of the ray trajectory showing how the arc-type of the trajectory in the reflection region changes to the loop-type as the wave propagation regime changes from a quasilongitudinal to quasi-resonance one (cf.bottom right panel in Fig.3).

Fig. 12 .
Fig. 12. Example of simulated overview spectrogram.Magenta dashed line indicates the LHR frequency along the satellite path.The yellow dashed line corresponds to a quarter of the equatorial electron cyclotron frequency on the current L-shell and serves as an additional reference frequency on the spectrogram.