AnnalesGeophysicae(2004)22:3589–3606 Annales SRef-ID:1432-0576/ag/2004-22-3589 Geophysicae ©EuropeanGeosciencesUnion2004 Characteristic properties of Nu whistlers as inferred from observations and numerical modelling D.R.Shklyar1,J.Chum2,andF.Jir˘´ıc˘ek2 1IZMIRAN,Troitsk,MoscowRegion,142190,Russia 2InstituteofAtmosphericPhysics,Acad. Sci. CzechRepublic,Boc˘n´ıII,14131Prague4,CzechRepublic Received: 16January2004–Revised: 8June2004–Accepted: 11June2004–Published: 3November2004 Abstract.ThepropertiesofNuwhistlersarediscussedinthe thepresent. Thefirstandmostprofoundsummaryofthere- lightofobservationsbytheMAGION5satellite,andofnu- searchinthisfieldwasgiveninabookbyHelliwell(1965), mericallysimulatedspectrogramsoflightning-inducedVLF whichwasasuperlativecontributiontowhistlerstudies. emissions. The method of simulation is described in full. The investigation of nonducted whistler-mode waves in Withtheinformationfromthisnumericalmodelling,wedis- the magnetosphere, in particular of MR whistlers and Nu tinguish the characteristics of the spectrograms that depend whistlers, which are the subjects of this paper, also has a onthesiteofthelightningstrokesfromthosethataredeter- long history. We will mention only some work that is di- minedmainlybythepositionofthesatellite. Also,weiden- rectlyrelatedto−orespeciallyimportantfor−thepresent tifytheregioninthemagnetospherewhereNuwhistlersare study.Anunexpectedpossibilityforwhistler-wavereflection observedmostoften,andthegeomagneticconditionsfavour- when the ions are taken into account in the dispersion rela- ing their appearance. The relation between magnetospher- tion,andthevisualisationofthiseffectbyraytracing,were ically reflected (MR) whistlers and Nu whistlers is demon- first demonstrated by Kimura (1966). In a sense, this find- strated by the gradual transformation of MR whistlers into ing predicted magnetospherically reflected (MR) whistlers, Nu whistlers as the satellite moves from the high-altitude which were found in the spectrograms of wave data from equatorial region to lower altitudes and higher latitudes. OGO 1 and 3 (Smith and Angerami, 1968). In their study, The magnetospheric reflection of nonducted whistler-mode mainly devoted to MR whistlers, Smith and Angerami also waves, which is of decisive importance in the formation of pointedoutthatthespectrogramofanMRwhistlerobserved Nuwhistlers,isdiscussedindetail. farfromtheequatormayhavetheshapeoftheGreekletter ν. TheycalledthisphenomenonNuwhistlerandsuggested Key words. Magnetospheric physics (plasmasphere) – Ra- its basic mechanism. According to these authors, the min- dioscience(radiowavepropagation)–Spaceplasmaphysics imum frequency on a Nu-whistler spectrogram, where the (numericalsimulationstudies) twobranchesmerge,correspondstothewavethatundergoes magnetospheric reflection at the observation point (see also Edgar, 1976). Magnetospheric reflection occurs when the 1 Introduction waves reach some point where their frequency is less than the local lower-hybrid-resonance (LHR) frequency, so it is Of the several natural sources of VLF waves in the magne- alsoreferredtoasLHRreflection. tosphere,lightningstrokesarethemostfamiliar. According AlthoughthecloserelationbetweenMRwhistlersandNu to the commonly accepted notion, a lightning stroke emits whistlers was established in the initial work by Smith and electromagneticwavesintotheEarth-ionospherewaveguide. Angerami (1968), in later work MR whistlers were stud- While propagating in this waveguide, some of them leak ied much more than Nu whistlers; see, for instance, Son- through its upper boundary and penetrate into the magne- walkar and Inan (1989), Draganov et al. (1992), Thorne tosphere, giving rise to whistlers observed in the opposite and Horne (1994), and Jasna et al. (1990). Some com- hemisphere. The investigation of ionospheric and mag- ments on these studies may be found in the paper by Shkl- netospheric wave phenomena related to lightning strokes yar and Ji˘r´ıc˘ek (2000), where an analysis of MR whistlers began from classical research by Eckersley (1935) and observedbyMAGION4and5wassupplementedbyanex- Storey (1953a), among others, and it has continued up to tensive numerical simulation of MR-whistler spectrograms. Correspondenceto: D.R.Shklyar Since then, several authors have used numerical simulation ([email protected]) ofspectrogramsintheirstudiesofMRwhistlers.Lundinand 3590 D.R.Shklyaretal.: CharacteristicpropertiesofNuwhistlers Krafft (2001) demonstrated the similarity of MR-whistler 2.1 Dispersionrelationandgroupvelocity spectrogramsthatappearsinacertainrangeofL-shellsand latitudeswhenthefrequencyscaleofthespectrogramisnor- The equations of geometrical optics, for the ray position r malizedwithrespecttotheequatorialelectroncyclotronfre- andthewavenormalvectork ofawavepacketwiththefre- quencyontheL-shellofobservation. Ji˘r´ıc˘eketal.(2001)in- quencyω,canbeexpressedinHamiltonianformas vestigatedtheinfluenceoftheplasmapauseonMR-whistler dr ∂H(k,r) dk ∂H(k,r) spectrograms,concludingthatthepresenceofapronounced = ≡v ; =− , (1) g plasmapause renders the traces on the spectrograms indis- dt ∂k dt ∂r tinct,sothe“classical”patternshouldbeobservedonlyunder wheretheHamiltonianH(k,r)isgivenbythelocaldisper- quiet geomagnetic conditions. An essential contribution to sionrelation the numerical modelling of MR-whistler spectrograms was madebyBortniketal.(2003), whoincludedwaveintensity H(k,r)=ω(k,r), (2) inthefrequency-timeplots,thusmakingthemmorelikereal spectrograms. andvg isthegroupvelocity. Equations(1)arewrittenaboveintheirgeneralform. We Afurtherstepinthenumericalmodellingofspectrograms nowspecifythedispersionrelationforwhistler-modewaves wastakenbyChumetal.(2003),whoshowedthatnumerical and the resultant expression for the group velocity, which simulations can be used to model spectrograms not only on govern the wave propagation in the approximation of geo- ashorttimescaleoftheorderof10s,theso-calleddetailed metrical optics, and which we use in our computer simu- spectrograms, but also to model overview spectrograms of lations. The dispersion relation, which expresses the wave datatakenalongasatellitepathduringtensofminutes,pro- frequency as a function of wave-normal vector and plasma videdthatlightning-inducedwhistlersarethemainemission parameters, can be obtained from the general equation for intheregiontraversedbythesatellite. Inthiscase,whistler thewaverefractiveindexinacold,magnetizedplasma(see, emission, trapped in the magnetosphere by LHR reflection, e.g. Ginzburg and Rukhadze, 1972). In the frequency band evolves into oblique noise bands above the local LHR fre- quency;thesearequalitativelyreproducedbynumericalsim- ωci(cid:28)ω∼<ωc, which is that of the whistler mode (ωci is the ioncyclotronfrequencyandω isthemagnitudeoftheelec- ulationsofoverviewspectrograms. Weshouldmentionthat c troncyclotronfrequency),andinplaceswheretheplasmais LHRreflectionalsoplaysanimportantroleforseveralother dense(ω (cid:29)ω ,whereω istheelectronplasmafrequency) types of VLF emission in the magnetosphere. Besides MR p c p which it is in most of the Earth’s plasmasphere, the disper- whistlers, Nu whistlers, and the LHR noise bands, where sionrelationmaybewrittenintheapproximateform: LHR reflection is the governing factor, it is also important faol.r(c2h0o0r3u)sfrwoamveths,eiarsawnaalsyspiosionfteCdLoUuStTreEcRendtalytab.y Parrot et ω2 =ω2 k2 + ω2 kk2k2 LH k2+q2 c (k2+q2)2 InthispaperweconcentrateonNuwhistlers,andproceed ω2 ω2cos2θ as follows. Section 2 is devoted to an analytical descrip- ≡ LH + c , (3) tionofnonductedwhistler-wavepropagation,withattention 1+q2/k2 (1+q2/k2)2 focussed on wave reflection at or well below the LHR fre- where the lower hybrid resonance (LHR) frequency ω is LH quency. Thekeypointsonhowraytracingintheframework givenby ofgeometricalopticscanreproducethespectrogramsofthe otibosner4vepdreesleencttsroemxpaegrniemtiecnfitaellddaatraeodniscNuusswedhiisntlSeresctf.r3o.mStehce- ω2 = 1 ωp2ωc2 ; 1 = me X nα , (4) LH M (ω2 +ω2) M n m MAGION5satelliteandcomparesthemwithcomputersim- eff p c eff e ions α ulations. Usingtheinformationthatmaybeapparentonthe (n ,m aretheelectronconcentrationandmass,respectively, modelledspectrograms,butcannotbeseenonrealones,the e e whilen ,m arethesamequantitiesforionsofthespecies mainpropertiesofNuwhistlersareexplained. Ourfindings α α andconclusionsaresummarisedinSect.5. α), k2=kk2+k⊥2, wherekk andk⊥ arethecomponentsofthe wave-normal vector parallel and perpendicular to the ambi- entmagneticfield,respectively,θ=cos−1(kk/k),and q2 =ω2/c2 , (5) p 2 Some features of nonducted whistler-wave propaga- where c is the speed of light. From Eq. (3) one can see tioninthemagnetosphere that the characteristic value of the wave number in the whistler frequency band is q≡ω /c, and that for a given p In this section, we discuss some aspects of nonducted wave-normalangleθ,thedependenceofthewavefrequency whistler-wave propagation in the plasmasphere that are es- on k involves only the ratio k/q. For the so-called quasi- sential for understanding the phenomena discussed in this longitudinal waves (Ratcliffe, 1959; Helliwell, 1965) k.q, paper. whereastheinequalityk(cid:29)q correspondstoquasi-resonance D.R.Shklyaretal.: CharacteristicpropertiesofNuwhistlers 3591 waves (see, for example, Walker, 1976; Alekhin and Shkl- 2.2 Magnetosphericreflectionofwhistler-modewaves yar, 1980). Some features of quasi-longitudinal and quasi- resonancewavepropagation,usefulforunderstandingthere- The possibility that whistler waves might be reflected sults of numerical simulations based on the dispersion rela- within the magnetosphere was suggested and studied by tionEq.(3),werediscussedbyShklyarandJi˘r´ıc˘ek(2000). Kimura (1966). In the one-dimensional case, wave reflec- Theexpressionsforthelongitudinal(paralleltotheambi- tion corresponds to a change in sign of the group velocity. entmagneticfield)andtransverse(perpendiculartotheam- In the two-dimensional case the situation is more compli- bient magnetic field) components of the group velocity that cated. If, for example, the longitudinal component of the followfromEq.(3)are: group velocity vgk greatly exceeds the transverse one vg⊥, then the wave reflection corresponds to the point where vgk vgk ≡ ∂∂kωk = ωkqk2(1+ωkL22H/q2)2 cvhgka∼ngvegs⊥i,tsaswiganv,eanmdatyhubsevrgekfl=ec0t.eHdowwitehverer,spinecthtetocaosneewchoeorre- + kk ωc2 k2+kk2 + k2k⊥2 ! (6) dreisnpaetectbtuot tchoentointhueerpcrooopradginataitneg. iFnotrheexsaammpeled,irienctaiodnipwoiltahr ωq2(1+k2/q2)3 q2 q4 magneticfield,thewavereflectionwithrespecttotheheight (ortheradialdistancer fromtheEarth’scenter)takesplace vg⊥ ≡ ∂∂kω⊥ = ωkq⊥2(1+ωkL22H/q2)2 vwgh⊥en−2vgk tanλ =0, (9) + k⊥ ωc2kk2 (1−k2/q2). (7) where λ is the geomagnetic latitude; when condition (9) is ωq4(1+k2/q2)3 satisfied, the radial component of the group velocity dr/dt vanishes. From Eq. (6) one can see that vgk has the same sign as Reflectionwithrespecttobothcoordinateswouldrequire kk,hencebothquantitieschangesignatcosθ=0;obviously, from Eq. (3), this can happen only when ω<ω , in which vgk=vg⊥=0whichisimpossibleforwhistlers. Indeed,from case vgk=0 for ω2=ωL2H/(1+q2/k2) (cf. Eq.L(H3)). As for Eqs.(3)and(6)itfollowsthatvgk=0for visg⊥di,reitctheadsotphpeossaitmeetosikg⊥nfaosrkω⊥2>foωrL2kH</(q1,−wqh4i/lek4f)orankd>vqiciet ω2=ωL2H k2k+2q2 < ωL2H , (10) versa. Thelaststatementbecomesapparentifwerewritethe expression Eq. (7) for vg⊥, eliminating kk2 with the help of whilevg⊥=0implies(seeEq.(8)) thedispersionrelationEq.(3): k4 ω2 =ω2 > ω2 . (11) vg⊥ =−ωk2(kk2⊥+q2)2[k4(ω2−ωL2H)−q4ω2]. (8) Thus,stLriHctkly4s−peqa4king,thLeHreflectionofawhistlerwavecan neverhappen: magnetosphericreflectionisinfactareversal Also,whenkk=0andvgk=0,itfollowsfromEqs.(8)and(3) ofthegroupvelocityinasmallregionofspaceastheresult that vg⊥ is always parallel to k⊥, since, under these condi- ofrefraction(Kimura,1966). tions,ω2<ω2 . LH To find the conditions for magnetospheric reflection, we Concerningvgk,itiseasytoseethatthefirsttermintheex- first describe this phenomenon more rigorously, contrasting pressionEq.(6)forvgk isalwaysmuchlessthanthesecond it with regular refraction. First of all, we shall speak of re- term. Indeed, for k2/q2&1, the ratio of the first term to the flectionasbeingapropertyofaray,regardlessoftime. If,in second is less or of the order of ω2 /ω2(cid:28)1. For k2/q2<1, LH c aregionthatissmallcomparedtothecharacteristicscaleof thefirsttermisoftheorderof plasmainhomogeneity, amajorvariationofthedirectionof kk ωL2H the group velocity takes place, whereas before entering and ωq2 after leaving this region the direction of the group velocity varies relatively slowly, then we shall call this event wave whilethesecondoneisoftheorderof reflection. Theforegoingconditionscanbeexpressedas kk ωc2 k2 , |v2−v1|∼max(|v1|, |v2|), (12) ωq2 q2 where v and v are the group velocities at the entrance to 1 2 thustheybecomecomparableonlywhen and at the exit from the reflection region, respectively. The groupvelocityv isafunctionofbothkandr,v =v (k, r); kk2 < k2 ∼ ωL2H . however,wecangneglectthevariationofr inagsmalglreflec- q2 q2 ωc2 tionregion. Hence,theamountbywhichthegroupvelocity varies in passing through the reflection region may be esti- However,suchsmallvaluesofk2/q2areoutsidetherangeof matedwiththehelpofEq.(1)as validityoftheapproximatedispersionrelationEq.(3),since they correspond to frequencies of the order of the ion cy- 1v ’−∂vgi ∂ωδt ∼(cid:18)∂vgi(cid:19) ω l , (13) clotronfrequency. gi ∂k ∂x ∂k L<v > j j j max g 3592 D.R.Shklyaretal.: CharacteristicpropertiesofNuwhistlers Contours of normalized frequency f/f Contours of (dV / dk ) ( ω / <V > V ) LHR ll ll g max 0.5 6 0.5 90 80 5 70 4 60 q k/ll 0 3 0 50 40 2 30 1 20 −0.5 −0.5 10 −5 0 5 −5 0 5 Contours of (dV / dk )( ω / <V > V ) Contours of (dV / dk )(ω / <V > V ) ⊥ ll g max ⊥ ⊥ g max 0.5 0.2 0.5 0.2 0.15 0.15 q k/ll 0 0 0.1 0.1 −0.5 0.05 −0.5 0.05 −5 0 5 −5 0 5 k /q k /q ⊥ ⊥ Fig.1. Contoursofnormalizedfrequencyandreflectionparameters. Notedifferentcolourbarsassociatedwithdifferentsubplots. Onlythe parametershownontheupperrightpanelgreatlyexceedsunityand,thus,determinesthewavereflection. where L is the characteristic scale of the plasma inhomo- locity. TheresultsareshowninFig.1,whereweseethatthe geneity,l isthelengthofthepartoftherayinthereflection reflectionisdeterminedbytheparameter region,<v >istheaveragemagnitudeofthegroupvelocity (cid:12) (cid:12) g (cid:12)∂vgk ω (cid:12) such that l/<v > is the duration of the reflection process, (cid:12) (cid:12) , (15) g (cid:12) ∂kk <vg >vgmax(cid:12) andthesubscript“max”denotesthemaximumvalue. Using Eq.(13)andthenotation whichatkk→0greatlyexceedsunity,whiletheotherquanti- tiesproportionalto∂vgk/∂k⊥and∂vg⊥/∂k⊥are∼<1overthe max(|v |, |v |)=v whole plane (k⊥, kk). We should emphasize that the value 1 2 gmax oftheparameterEq.(15)dependsonhowwedefinethesize werewritethereflectionconditionsofEq.(12)intheform: ofthereflectionregion,soitisnotdetermineduniquely. Its only important property is that it is much larger than unity, (cid:12)(cid:12)(cid:12)(cid:18)∂vgi(cid:19) ω (cid:12)(cid:12)(cid:12)∼ L (cid:29)1. (14) whichensuresthatthedirectionofthegroupvelocityvaries (cid:12) ∂k <v >v (cid:12) l rapidlyalongtherayinthereflectionregion,comparedwith j max g gmax itsbehaviouronotherpartsoftheray. Thus,thewavereflectiontakesplaceforthosekkandk⊥,and Aswehaveseenabove,forkk=0,theparallelcomponent the corresponding wave frequencies, for which the quantity ofthegroupvelocityvgk=0,andthewavefrequencyisdeter- on the left-hand side in Eq. (14) is much larger than unity. minedbyEq.(10).Fromthisequationitfollowsinparticular ThisquantityhasbeencalculatednumericallyusingEq.(3), that in the quasi-resonance regime k2(cid:29)q2, the wave reflec- togetherwiththeexpressionsEqs.(6), (7)forthegroupve- tion takes place at frequencies close to the LHR frequency D.R.Shklyaretal.: CharacteristicpropertiesofNuwhistlers 3593 contrary, if before and after reflection the wave remains to Diagram determining the type of wave reflection 5 therightofthedottedline,whichistypicalofk/q(cid:29)1,then 4.5 vg⊥changessign,andtherayhastheshapeofaloop. These featuresofthewavereflectionareillustratedbyFigs.3and 4 4. InFig.4, L3 (ratherthanthemorenaturalquantityL)is 3.5 chosenasoneofthecoordinatesinordertomaketheloops intheraymoreobvious. Figures3and4correspondtoa5- 3 KHzwavestartingverticallyat15◦ geomagneticlatitude,at H ωω/L2.5 theheight500km;theplasmasphereissmoothandtheprop- agationtimeissetto3s. 2 1.5 1 3 Spectrogrammodellingbymeansofraytracing 0.5 VLF data from Magion 5 will be presented below in the 00 1 2 3 4 5 6 7 8 9 10 form of spectrograms. These were made with a sampling k/q frequencyof44100Hzandanintegrationtimeof23.22ms; thus, eachinstantaneousspectrumwasevaluatedfrom1024 Fig.2.Reflectiondiagram. data points. The corresponding resolution in frequency is ∼<100Hz. Each spectrogram comprises about 300 instanta- neous spectra and covers a time interval of 7s. It is a rep- ωLH,whereasfork2∼<q2 thewavefrequencyωmaybewell resentationofspectralintensityinthefrequency-timeplane, belowω . LH with time along the x-axis, frequency along the y-axis, and InFig.1,thecontoursofnormalizedfrequencyandofthe the intensity indicated by the degree of darkness on black- reflection parameters are shown on the (k⊥,kk)-plane. Al- and-white spectrograms, or by the use of colour. If the thoughthewavefrequencyremainsconstantwhenthewave spectral intensity is appreciable only along some curves in propagatesinastationaryinhomogeneousmedium,thisdoes the (f,t)-plane, as is the case for MR and Nu whistlers, not mean that it remains on the same contour line of the theproblemofspectrogrammodellingconsistsoftwoparts: normalizedfrequency,sincethenormalizingLHRfrequency firstly,constructingthefrequency-timeplot,whichmayhave maychange. Obviously,insteadof(k⊥,kk),twootherquan- manybranches,ofcourse; andsecondly,attributingthecor- titiesmaybechosenastheindependentvariablesdetermin- responding intensity to each curve. Here we discuss how it ingthewavecharacteristics. Inparticular,itisconvenientto is done by means of ray-tracing calculations based on the analyzethefeaturesofthewavereflectionandpossibletypes equationsofgeometricaloptics. Sinceverymanyraysmust of raytrajectories in the reflection regionwiththehelp of a becalculatedinordertoreproducethemainfeaturesofNu- diagramonthe(k/q,ω/ωLH)-plane,asshowninFig.2. whistlers on a model spectrogram, we use relatively simple In this analysis, we will assume that k⊥>0, i.e. that the modelsforthegeomagneticfieldandforthedistributionsof wave-normalvectorisdirectedtowardshigherL-shells. The plasma density and LHR frequency, all given by analytical solid line in the figure is determined by Eq. (10) and corre- expressions(seeShklyarandJi˘r´ıc˘ek(2000)fordetails). sponds to kk=vgk=0. According to Eq. (3), the same line defines the minimum possible wave frequency as the func- 3.1 Constructingthefrequency-timeplot tionofk/q,sothedispersionrelationhasnorootsbelowthis line. Aswehaveseenabove(cf. Fig.1.),largevaluesofthe We assume that a thin layer in the upper ionosphere is illu- reflectionparameterEq.(15),typicalofwavereflection,are minatedbywavesfromalightningstroke,andthatthispro- attainedinthevicinityofkk=vgk=0. Thus, onthediagram cess is effectively instantaneous on the time scale of wave in Fig. 2, the reflection region is represented by the narrow propagation to the satellite. We also assume that initially region above the solid line. The dotted line is determined all waves have their normal vectors directed vertically, due byEq.(11)andcorrespondstovg⊥=0. Intheregiontothe to refraction by the ionosphere. Similar assumptions have leftofthisline,vg⊥ hasthesamesignask⊥ (positiveunder beenusedinalloftheworkonspectrogrammodellingcited ourassumption), whileintheregiontotherightofthisline above. Sincetheverticaldimensionoftheilluminatedlayer ithastheoppositesign. Clearly,whenthewaveapproaches (which, in turn, plays the role of an illuminating region for the reflection region, it always moves from higher towards themagnetosphere)ismuchsmallerthanitsdimensioninthe lower values of ω/ω on the diagram in Fig. 1. Another horizontalplane,andsincetheverticaldirectionofthewave- LH importantpointisthatinthereflectionregionvg⊥ isalways normalvectorsimpliesthatthewavespropagateinthemerid- positive,whilevgkchangesitssign. Thus,ifbeforeandafter ianplane,computationoftheraysisnowatwo-dimensional reflection the wave remains in the shaded region to the left problem, with initial conditions given on some line that ap- ofthedottedline,whichistypicalofk/q∼<1,thenvg⊥ does proximates the thin layer. As such a line, we take a part of notchangesign,andtherayhastheshapeofanarc. Onthe thearcattheheightof500kmabovetheEarth’ssurface. 3594 D.R.Shklyaretal.: CharacteristicpropertiesofNuwhistlers Trajectory of the wave 2 1 1 0.5 λ )e n R 0 a 0 Z ( = t 2 Z −1 3 −0.5 L=4 −2 −1 0 1 2 3 4 0 5 10 15 X (R ) e 0.4 200 0.2 100 dl Vg ll 0 V/g ll d 0 −0.2 −0.4 −100 0 5 10 15 0 5 10 15 0.06 8 0.04 6 V⊥g 0.02 22/q4 k 0 2 −0.02 0 0 5 10 15 0 5 10 15 l (length along the ray) l (length along the ray) Fig.3.Raytrajectoryandwavepropagationcharacteristicsfor5-kHzwave.Comparisonofthelatitudevariationandthereflectionparameter shownontheupperandmiddlerightpanels,respectively,clearlyshowsthatthereflectiontakesplacewhenthereflectionparameterhasa peak. For numerical modelling of spectrograms, Storey sug- front is the surface at the center (i.e. half-way through) of gestedusingthenotionofagroupfront. Herewereproduce, this sheet. Within any given ray tube, the disturbance is a with his permission, his definition and physical explanation wavepacketmovingalongthetubeatthegroupvelocity,and ofthisnotion. within this packet, the point of maximum amplitude lies on thegroupfront(Storey,2003,privatecommunication). Foranyparticularfrequency,considerallthepossiblerays thatcanbetracedupwardsfromtheilluminatingregion,with Inthecaseunderdiscussion,initially,thegroupfrontsfor initialconditionsasdefinedabove. Imaginethatalongeach allfrequenciescoincidewiththeilluminatingregion,i.e.the ray, starting at the instant of the lightning stroke, a point partofthearcextendingoverarangeoflatitudesat500km movesawayfromtheilluminatingregionatthelocalgroup height. With increasing time after the lightning stroke, the velocity. Then,atanylaterinstant,thesetofallsuchpoints group fronts separate due to the different group velocities defines a surface: this is the group front for the frequency ofwaveswithdifferentfrequencies,whileeverygroupfront concerned. is deformed due to plasma inhomogeneity, and also due to thedifferentinitialconditionsforthewavesofthesamefre- A more physical way of visualizing the group front is quencystartingverticallyatdifferentlatitudes. to imagine that the lightning stroke emits a narrow-band impulse instead of a wide-band one, thus giving rise to a Toplotapointinthe(f,t)-planeofaspectrogramofdata quasi-monochromatic disturbance that propagates through froma satellite, weshouldfind thetimeat which thegroup the magnetosphere in the form of a thin sheet. The group front crosses the satellite position. This procedure can be D.R.Shklyaretal.: CharacteristicpropertiesofNuwhistlers 3595 formalized as follows. The equations of geometrical optics Wave trajectory for 5 kHz wave on (L3, λ)−plane in their general form have been written in the previous sec- 40 tion(seeEqs.(1),(2)). Asiswellknown(see,forexample, Landau and Lifshitz, 1976), when the Hamiltonian H does 30 not depend on time, it is a constant of the motion. Thus, according to Eq. (2), Eq. (1) describe a wave packet with 20 constant frequency. We should emphasise that in the 2-D λ casethewavefrequencyalonedoesnotdeterminethewave de 10 u packetuniquely. atit To solve Eq. (1), it is most convenient to use canonically ant l 0 conjugate variables. However, once the solution has been vari−10 n found, itcanbeexpressedintermsofanyvariablesthatare I uniquelyrelatedtothecanonicalones. Thegeneralsolution −20 oftheEq.(1)hastheform −30 r =r(r , k , t); k =k(r , k , t) (16) 0 0 0 0 −40 1 2 3 4 5 6 7 8 9 whileω(k, r)=ω(k0, r0). Inthe2-Dcaseconsidered,both L3 k and r are two-dimensional vectors. Moreover, since we start all rays from a single altitude with the wave normals Fig. 4. Zoomed-in view of the ray trajectory showing how the vertical, there are in fact only two independent initial vari- arc-type of the trajectory in the reflection region changes to the ables, and we may choose the wave frequency to be one of loop-type as the wave propagation regime changes from a quasi- them. Asthesecondinitialvariablewechoosetheinitialge- longitudinal to quasi-resonance one (cf. bottom right panel in omagneticlatitudeλ ,asisusualincomputersimulationsof Fig.3). 0 thiskind.Then,takingtheMcIlwainparameterLandthege- omagneticlatitudeλastwocoordinates, wecanrewritethe and L, the same is true for k. Thus, all the characteristics solutionEq.(16)intheform of the wave packets that contribute to the spectrogram at a givensatellitepositionbecomefunctionsofωandofthehop L=L(λ , ω, t);λ=λ(λ , ω, t); number,andcanbedisplayedifdesired. 0 0 k =k(λ , ω, t). (17) 0 3.2 Calculationofspectralintensity When the solution in the form Eq. (17) is known, all of the Athoroughdiscussionoftherigorouswaysofdisplayingin- localwavecharacteristics,suchasthegroupvelocity,there- tensityonspectrogramswouldleadustoofarawayfromthe fractiveindexandthewave-normalanglemaybefound;the main topic of the present paper: it will be presented else- wavefrequencyisconstantalongtheray,atthevaluechosen where. Herewediscussonlythemainaspectsofthisprob- initially. lem. We regard the wave field as the sum of a set of wave TherelationsEq.(17),beingthesolutionoftheequations packets propagating in the magnetosphere with their group of motion, are unique functions of their independent vari- velocities.Thecentralfrequencyofeachwavepacketiscon- ables. The first two relations in Eq. (17), which define in a served, while its wave-normal vector varies along the ray, parametricwaythetimet andtheinitiallatitudeλ asfunc- 0 satisfyingalocaldispersionrelationateachpoint.Inasense, tions of ω, λ, and L, can, in principle, be solved for t and thewavepacketisdeterminedbyabunchofclosetrajectories λ : 0 in the phase space (k,r) whose projection onto the coordi- t =t(ω; λ, L); λ =(ω; λ, L). (18) nate space represents the ray tube. The ray itself, and the 0 variation of the wave-normal vector along it, are described These functions, however, may have many branches, that is bytheequationsofgeometricalopticsEq.(1). to say, they may be multi-valued. As we shall see below, As is well known (see, for instance, Fermi, 1968), the thedifferentbranchesofthesolutionEq.(18)correspondto wave packet in geometrical optics is an analog of the mass different numbers of hops across the equator that the wave pointinmechanics;assuch,itischaracterizedbytheinitial packetsperforminthemagnetosphere. Thefirstfunctionin coordinates of its amplitude maximum and its wave vector Eq. (18) defines the time when the group front for the fre- at this point. However, in contrast to a mass point, a wave quency ω crosses the satellite position λ, L; thus, it yields packetisoffinitesize;itisalsocharacterizedbyitswidthin the time-frequency curves on the spectrogram. The second k-spaceandthecorrespondingreciprocaldimensionincoor- function determines the initial latitude for the frequency ω dinate space. Thus, in the general 2-D case, a wave packet on each branch. This latitude can easily be displayed on a is characterized by four parameters. However, in the case modelspectrogram, which, ofcourse, isimpossibleforreal underconsideration,whenallraysstartverticallyat500km ones. Moreover, since t and λ are now functions of ω, λ, altitude, only two parameters are needed to characterize a 0 3596 D.R.Shklyaretal.: CharacteristicpropertiesofNuwhistlers wave packet, and these can be chosen, as above, to be the v =(v2 +v2 )1/2. Obviously, the wave packet passes over g g⊥ gk wave frequency ω and the initial latitude λ . Nevertheless, thepointsduringthetimeintervalfromt (s)tot (s),where 0 1 2 therearemanywavepacketswiththesamefrequency,start- ingatdifferentinitiallatitudesλ . Thereceiveronthesatel- t1(s)=min[t−(s),t+(s)]; t2(s)=max[t−(s),t+(s)]. (21) 0 litemeasuresthetotalfield,andevenafterspectralanalysis, We then integrate Eq. (20) over t from t (s) to t (s). Since 1 2 itispossiblethatmorethanonewavepacketmaycontribute the integrand tends to zero at both limits of integration, the to this field. (We remind the reader that, although different contribution from the first term vanishes. For the same rea- raysneverintersectinphasespace,theirintersectionincoor- son, the integral over t can be shifted into the argument of dinatespaceisnotforbiddenbytheequationsofgeometrical thederivativewithrespecttos. Asaresultweobtain optics.) Thekeypointthatsimplifiesspectrogrammodelling ! inourcaseisthat,astheraytracingshows,therayswiththe d Z t2(s) σ(s)v (s) U (t0, s)dt0 =0. (22) samefrequencythatstartatdifferentlatitudesneverintersect ds g ω t1(s) incoordinatespace. Thismeansthatthesatellite,whichmay beconsideredasafixedpoint,neverreceivesmorethanone Thus,thequantity wave packet at any time. Thus, in calculating the spectral Z t2(s) intensityatagivenfrequency,weneedtoconsideronlyone W(ω, λ )≡σ(s)v (s) U (t0, s)dt0 0 g ω wave packet, provided that the duration of the time interval t1(s) over which the spectrum is evaluated is much less than the isconservedalongtheray. Thisquantityisthetotalenergy typical bounce period of the wave packets in the magneto- ofthegivenwavepacket. sphere, which we have always found to be the case. (Ob- The wave energy density U is related to the electric field viously,differentwavepacketsmaycometothesatelliteon Easfollows: differenthops.) Lett−andt+bethetimesatwhichthegroupfrontsforthe U=w(s,ω,θ)|Eω(t,s)|2 , frequencies ω−δω and ω+δω, respectively, cross the satel- where E (t,s) is the electric-field component of the wave lite position. We can then state that the wave field received ω packetmeasuredbythesatellite,whilethefactorw(s,ω,θ) by the satellite during the time interval |t+−t−| is that of depends on frequency, wave-normal angle, and the local some wave packet with central frequency ω and bandwidth plasma parameters, and also on the wave mode, of course. 1ω=2δω. If the time of spectral evaluation 1t is less than Atthispointweassumethatthesatellitemeasuresthecom- |t+−t−|, then the spectral intensity in the frequency band ponentoftheelectricfieldperpendiculartotheEarth’smag- (ω−δω,ω+δω)willbenonzerooverthewholeintervalbe- neticfieldB inthe(k,B )-plane. Then,inthesamerange tween t− and t+; in the opposite case, the spectral intensity 0 0 of parameterswhere the dispersion relation Eq. (3) is valid, willbenonzerothroughoutsomeintervalofduration1t that theexpressionforwhastheform includes(t++t−)/2.Hereweassumethatthefrequencyres- olutionis∼1ω,andthus1ω1t ∼>2π. 1 ω2ω2 Incalculatingthespectralintensityfordisplayonaspec- w = p c 8π (ω2−ω2)2 trogram, we need to take into account the variation of the c " # wave-packetamplitudealongtheraycausedbygeometrical × 1+ ωp2 (ωc2+ω2) 1 + ε22 , (23) factors. Todothis,weproceedasfollows. Consider,forthe ω2 (ω2−ω2)(N2−ε ) (N2−ε )2 c 1 1 frequency ω, the ray that passes through the position of the satellite. Lets bethedistancealongthisraytoanypointon where fito,ramn)dinletthEeωfr(et,quse)ncbyebthaendw(aωv−e-δfiωel,dωc+omδωp)ontheanttth(tehesawteallviete- N2 = kω2c22 ; ε1 = ω2ω−p2ω2 ; ε2 =−ω(ωω2p2−ωcω2). (24) would measure if it were at this point. As has been argued c c above, for a given s and hop number, this field belongs to Thequantityw mayberegardedasafunctionofs, ω,and oneparticularwavepacket,characterizedbyω andλ . The θ sinceinacoldmagnetoplasmatherefractiveindexN isa 0 energy-conservationlawforthiswavepackethastheform: functionofωandθ. Withthisnotation,theconservedquan- ∂U (t, s) titytakestheform ω,λ0 +div[v U (t, s)]=0. (19) ∂t g ω,λ0 Z t2(s) W ≡σ(s)v (s)w(s) |E (t0,s)|2dt0 =const. (25) This and later equations concern the particular wave packet g ω characterized by the two parameters ω and λ . Henceforth, t1(s) 0 however,thesecondparameterwillbeomittedforshortness. On the other hand, from the well-known theorem in spec- Equation(19)canberewrittenas tralanalysisthatrelatesthefieldcomponentEω(t0,s)ofare- ceivedwavepackettoitstime-dependentspectralamplitude ∂Uω(t, s) + 1 d (cid:2)σvgUω(t, s)(cid:3)=0, (20) E(ω, s, t),wehave ∂t σ ds where s is the coordinatealong the rayconsidered, σ is the Z t+1t/2 1ω |E (t0,s)|2dt0 =|E(ω, s, t)|2 . (26) crosssectionofathinraytubecenteredonthisrayand t−1t/2 ω 2π D.R.Shklyaretal.: CharacteristicpropertiesofNuwhistlers 3597 If t −t <1t, then the integrals in Eqs. (25) and (26) are 3.3 ComparisonwiththeapproachofBortniketal.(2003). 2 1 equal, while in the opposite case the integrals are propor- AswasmentionedintheIntroduction, Bortniketal.(2003) tionaltotheintervalsofintegration. Takingthesefactsinto made an important step in numerical modelling of MR account,weobtainfromEqs.(25)and(26): whistlersbyincludingspectralintensityintosimulatedspec- | E(ω, s, t) |2 = trograms. Like the spectrograms simulated by Bortnik et  al. (2003), ours now display spectral intensity. However, 2π 1t W ,t −t > 1t 1ω(t2−t1)σ(s)vg(s)w(s) 2 1 themethodweusediffersfromthatofBortniketal.(2003) = (27) in several respects. Firstly, we deal from the outset with  2π W , t −t < 1t . wave packets of finite spectral width 1f, corresponding 1ωσ(s)vg(s)w(s) 2 1 to the frequency resolution on real spectrograms. In this case, the time interval during which the frequency band Weseethat,apartfromthequantitiesdirectlydeterminedby f−1f/2,f+1f/2isreceivedonthesatelliteisdetermined the equations of geometrical optics, an additional quantity bythegroup-frontcrossingsofthesatelliteposition,assug- thatneedstobecalculatedisthecrosssectionoftheraytube. gestedbyStorey((2003),privatecommunication). Thistime Wewillassumethattherearenogradientsintheazimuthal is determined unambiguously, with no uncertainty; it does direction,sothewavespropagateinmeridionalplanes. Then not use the notion of detection area, the extent of which is thewidthoftheraytubeintheazimuthaldirectionis difficulttodefineconsistentlyduetothecontinuousmerging ofdifferentraysatthesamefrequency.Secondly,ashasbeen x18≡R Lcos3λ18, (28) E shownbyStorey(1953b),whenadispersedsignalispassed throughabankofnarrowbandfilters,thetemporalvariation where x is the Cartesian coordinate in the meridional plane of its instantaneous frequency is measured most accurately orthogonal to the dipolar axis of the Earth’s magnetic field, whenthebandwidthofthefilterequals Landλare,asbefore,theMcIlwainparameterandthemag- neticlatitude,respectively,and18istherangeofazimuthal anglesoverwhichthewavepacketextends, whichisacon- (cid:12)(cid:12)dfi(cid:12)(cid:12)1/2 (cid:12) (cid:12) , stantofitsmotion. Thus,anon-trivialpartofthecrosssec- (cid:12) dt (cid:12) tion is its width 1ξ in the meridional plane, which defines where f is the instantaneous frequency. As the wave phe- theray-tubecrosssectionσ accordingto i nomena that we model are characterized by a rate of fre- σ =x181ξ ≡R Lcos3λ181ξ . (29) quencyvariationoftheorderofafewkHzpersecond,thefil- E terbandwidthshouldbeoftheorderof50Hzforthesharpest To find the quantity 1ξ, let us consider two neighbouring output,sowechoosethisvalueasthefrequencystepinour rays. Let(x ,z )and(x ,z )betwoneighbouringpointson 1 1 2 2 calculations. Thus, we consider that the interpolation pro- thefirstandsecondray,respectively,andletψ betheangle cedure used by Bortnik et al. (2003), which yields a fre- betweenthegroupvelocityandthex-axis.Thenthewidthof quency resolution of ∼1Hz, is superfluous in this respect, theraytubeinthemeridionalplane,1ξ,is all the more so because, finally, they set the width of their (z −z )cosψ −(x −x )sinψ frequency bin to 50Hz. And thirdly, another difference be- 1ξ = 2 1 2 2 1 2 . (30) tween their approach and ours lies in the way we evaluate cos(ψ −ψ ) 1 2 the spectral intensity: instead of computing millions of in- Sincethevector(−sinψ , cosψ )isorthogonaltov and, terpolated rays, each weighted with a measure of wave en- 2 2 g thus,totheray,thisresultdoesnotdependontheparticular ergy, and then calculating the energy carried by those rays choiceofthepoint(x ,z ),providedthatitiscloseenough thatcrossthedetectionarea,wecalculatethevariationofthe 2 2 tothepoint(x ,z )wherethecrosssectioniscalculated. ray-tubecross-section,thenuseenergyconservationandPar- 1 1 AccordingtoEq.(27),thespectralintensity|E(ω,s)|2 at seval’s relation to translate the energy in each wave packet, theobservationpoints isdeterminedbythefactor(σv w) of bandwidth 50Hz, into spectral intensity displayed on a g s andtheconservedvalueW ofthewavepacketenergy. Thus, spectrogram. As for initial distribution of the wave energy to include the spectral intensity in model spectrograms, we among wave packets, we use the following model. We as- needtosupplementtheray-tracingcalculationwiththeeval- sume that each wave packet is determined by its frequency uation, for each ray, of the quantities v and w determined f and initial latitude λ , and that all wave packets have the g 0 by Eqs. (6), (7), and (23), respectively, and with the calcu- samefrequencywidth1f andoccupythesamespatialwidth lation of a neighbouring ray, which enables us to find 1ξ 1λ at the beginning. Since initially all wave packets have 0 Eq.(30)andthusthecrosssectionoftheraytubeEq.(29). verticaldirectionoftheirwavenormalvectors, andnegligi- The corresponding data base, which is similar to the one bledimensioninradialdirection,theseparametersdetermine described by Shklyar and Ji˘r´ıc˘ek (2000) but supplemented the wave packet uniquely. The total energy of each wave with the relative-intensity parameters, has been computer- packet,which,ofcourse,isconserved,ismodelledas ized. Spectrogramscalculatedwiththehelpofthisdatabase arepresentedinthenextsection. W ∝ϕ(λ )η(f), 0 3598 D.R.Shklyaretal.: CharacteristicpropertiesofNuwhistlers magnetospheric reflection, are mainly observed on L-shells fromL∼1.8toL∼3. Thelinesofconstantaltitudeandmag- netic 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 de- scending parts of its orbits, observed MR whistlers in the equatorialregionataltitudesfromabout1.3to2Earthradii, whichisfarfromtheregionswherethewavesarereflected. SuchMRwhistlershavebeendiscussedindetailbyShklyar andJi˘r´ıc˘ek(2000); oneexampleispresentedinFig.6,with its simulated counterpart shown in Fig. 7. On the contrary, overtheascendingpartsoftheMAGION5orbits,whistlers couldbeobservedintheregionsoftheirmagnetosphericre- flection. Beforeshowingexamplesofspectrogramstakenin theseregions,whichexhibittheν-shapedpatternscharacter- isticofNuwhistlers, werecallsomefeaturesofductedand Fig.5.RadiovisiblepartsoftheMAGION4andMAGION5satel- nonducted whistler wave propagation. As there is no clear- litetrajectories. cut boundary between these two types of propagation, it is sometimeshardtodistinguishbetweentheminsatellitedata, particularlyinthecasewherethewavespropagatefromthe whereϕ(λ )isasmoothfunction,whichdecreaseswiththe 0 Earthandarereceivedonasatellitebeforecrossingtheequa- distance from the center of illuminating region; and η(f), tor(fractional-hopwhistlers). ThedegreeofdispersionDof which describes the frequency dependence of the wave en- the fractional-hop whistlers is very small (∼(10−20)s1/2) ergy distribution, is adopted from the paper by Lauben et duetotheshortdistanceofpropagation,anditisalmostim- al.(2001): possibletodistinguishbetweenductedandnonductedprop- f2 agationinthiscase;(see,forinstance,inthespectrogramof η(f)= , (f2+0.63)(f2+253) bottompanelinFig.8,thesingletraceatthetime∼3s). The natureofthepropagationinthiscasecanbedeterminedonly wheref isthewavefrequencyinkHz. Thesamefrequency fromthefactthateachwhistlerisfollowedbyNuwhistlers, dependence has been used by Bortnik et al. (2003). On the withalmostthesamedelayinallsucheventsobserveddur- other hand, we do not take into account wave growth or ing tens of seconds. From time to time, the spectrograms damping,soinourcasethevariationofenergydensityalong show subsequent traces of reflected ducted whistlers, indi- therayisdueonlytogeometricalfactors. cating ducted propagation. When analysing the first mag- netosphericreflectiononspectrograms,oneshouldalsotake 4 Nu whistlers from the MAGION 5 satellite and their caretodistinguishbetweenNuwhistlersandthetracesfrom modelling doubleormultiplelightningstrokes. Forexample,thetraces shown in Fig. 8, in the second panel from the top, which As was noticed in the earliest satellite experiments, VLF resemble those of Nu whistlers, are in fact those of normal data from a satellite exhibit a much richer variety of wave whistlers originating from multiple lightning strokes in the phenomenathandatafromground-basedobservations. The opposite hemisphere. This can be established from the fact reason is that on satellites, phenomena related to the quasi- thatthetracesonthespectrogramallhaveexactlythesame resonance(ornonducted)typeofwhistler-wavepropagation form, whilethetimedelaybetweensuccessivetracesvaries are observed, as well as those related to quasi-longitudinal randomlyintime. (ducted) propagation, while ground-based data are mainly TheeventsshownonthethirdandfourthpanelsinFig.8 limitedtothelatter. are different. Here we see Nu whistlers in which the first Inthissection,wepresentexamplesofMRwhistlersand tracecorrespondstowavespropagatingdownwards,whereas Nu whistlers observed on board the MAGION 5 satellite. thesecondoneisformedbywavespropagatingupwardsaf- The data are available from June 1998 to July 2001. Since terMRreflection. Notethatthetracesarenotparallelinthis theyweretransmittedinanalogueformtothegroundstation case.Theseexamplesshowthatcertainwavephenomenaob- inPanskaVes(50.53◦N,14.57◦E)inrealtime,itwasthera- servedonsatellitescanbeidentifiedonlybyfollowingtheir diovisibilityofthesatellitethatlimitedthepartsoftheorbits evolutionandrecurrenceinthedata. fromwhichdatacouldbeobtained. Agraphicalillustration Figure 9 demonstrates how the spectrograms with MR of the parts on which VLF data were recorded is given in whistler traces change their character along the ascending Fig.5;thesmallerpartsonwhichMRandNuwhistlerswere partsoftheMAGION5orbits.Asthealtitudeandlatitudeof observedaremarkedbyasterisks. Onecanseefromthisfig- thesatelliteincrease(cf. Fig.5), thetimeintervalsbetween ure that MR whistlers, i.e. the waves that have undergone thetracesofsuccessivehopsincreasealso, evidentlydueto