Higher harmonic tweek sferics observed at low latitude: estimation of VLF reflection heights and tweek propagation distance

Lightning generated signals recorded at a lowlatitude station, Suva (18.2 ◦ S, 178.3 E) Fiji, in the South Pacific region, during September 2003–July 2004, are used to study the propagation features and the reflection heights of tweek atmospherics in the waveguide formed by the Earth’s surface and the lower ionosphere. Tweeks are observed only during the local night and the maximum harmonic ( n) recorded is six. The occurrence of tweeks with higher n p ogressively decreases as n increases. The dispersed part of tweeks decreases as n increases. The attenuation factor has been calculated for tweeks with n=1–3. The ionospheric reflection heights obtained assuming the transverse magnetic mode of propagation for tweek signals vary from 83–92 km. A higher harmonic of the same tweek is reflected from about 2.0 km higher than the lower harmonic. For 90% of tweeks, propagation distances are estimated to be between 1000– 5000 km. Tweeks with lower n propagate longer distances than the tweeks with higher n.


Introduction
The return strokes of lightning discharges are natural powerful transmitters of electromagnetic energy over a wide electromagnetic spectrum extending from few Hz to beyond the visible region (Prasad and Singh, 1982). Most of the lightning energy is radiated in the Extremely Low Frequency (ELF, 3-3000 Hz) and the Very Low Frequency (VLF, 3-30 kHz), with the peak spectral density of these atmospherics (or sferics) being in the frequency below 15 kHz (Volland et al., 1987;Ramachandran et al., 2007). ELF and VLF energy Correspondence to: S. Kumar (kumar su@usp.ac.fj) generated by the lightning discharge propagates by multiple reflections in the waveguide bounded by the Earth (ground or sea) and the lower region of the ionosphere (Barr et al., 2000). This guided propagation occurs with a low attenuation rate (a few dB/1000 km), allowing ELF radio sferics to be observed around the world from a single lightning discharge. The conductivity of the boundaries of the waveguide and the total path traveled in it result in appreciable dispersion of the sferics at the lower frequency end, with a lower cut-off at around 1.8 kHz. Such long-delayed sferics are referred to as "tweeks" in the literature (Yamashita, 1978), due to the distinct chirping sound they make when heard on a loudspeaker or earphone.
The properties of tweeks have been investigated and explained by a number of authors (e.g. Ohtsu, 1960;Yano et al., 1989Yano et al., , 1991Hayakawa et al., 1994Hayakawa et al., , 1995Sukhorukov, 1996;Ferencz, 2004;Ferencz et al., 2007). Tweeks have been reported to be left-hand polarized (L) waves (Reznikov et al., 1993). The strong gyrotropy of the lower ionosphere provides total internal reflection for L-waves and transmission of right-hand polarized (R) waves into the outer ionosphere as whistler mode waves which can be recorded by satellites or at ground-based stations in the opposite hemisphere. Tweeks generated by strong lightning can leak out of the Earth-ionosphere waveguide (EIWG) and excite the "Spiky Whistlers" (SpW) even after several thousand of km of subionospheric propagation. SpW whistlers have been recently observed by the DEMETER satellite (Ferencz et al., 2007).
Tweeks have been utilized by researchers to estimate the ionospheric reflection height (h), equivalent electron densities (n e ) at reflection heights, propagation distance (d), and the geographic locations of the source discharge (e.g. Reeve and Rycroft, 1972;Kumar et al., 1994;Shvets and Hayakawa, 1998;Ohya et al., 2003Ohya et al., , 2006. From the analysis of tweeks observed during February-March 1991 at Bhopal, a low-latitude ground station in India, Kumar et    al. (1994) found h and d to vary in the range of 83-89 km and 1500-2500 km, respectively. By analyzing the tweeks up to the 8th harmonics observed during January-April 1991, Shvets and Hayakawa (1998) estimated an increase in the n e values from 28-224 el cm −3 at h=81-83 km. From the firstorder mode cut-off frequency of tweeks, at low-mid latitudes in Japan, Ohya et al. (2003), estimated h and n e to vary in the range of 80-85 km and 20-28 el cm −3 , respectively.
In this work, initial findings on tweeks observed during September 2003-July 2004, at a low-latitude station Suva (geog. lat. 18.2 • S, geog. long. 178.3 • E), Fiji, in the South Pacific region, are presented. The geographic location of the recording station, on an island in the vast ocean, may provide an excellent opportunity to record tweeks with higher harmonics with good dispersion. The sferics data are recorded using the World Wide Lightning Location Network (WWLLN) VLF setup at the University of the South Pacific. There are 28 universities and national institutes all over the world participating in WWLLN (R. L. Dowden, personal communication).

Experimental data and analysis
The transverse magnetic (TM) mode has a vertical electric field component (E Z ) (Wait, 1957;Budden, 1962). The experimental set-up consists of a short (1.5 m) whip antenna to receive E Z of sferics. The pre-amplifier fixed at the bottom of the whip antenna is coupled to the ELF-VLF receiving unit (SU). A Global Position System (GPS) antenna is connected to GPS receiver built into the SU unit. The GPS receiver generates a pulse per second with an accuracy of 100 ns, and is used to measure the arrival time of the sferics. The output of the SU unit is fed to the sound card in a PC which records the data in files of 11 MB per minute using the lightning software. The details of the system have been described by Dowden et al. (2002).
During September 2003-January 2004, sferics data were recorded for a total of 30 min at different hours in the night and the day as well. From February 2004-July 2004, sferics data were recorded for 2 min at every hour only in the night. The ELF-VLF data files are analyzed using a MATLAB code which produces dynamic spectrograms of one-second duration. Tweeks sferics are visually identified from the spectrograms (see Fig. 1) and then analyzed in two steps. In the first step, for the first harmonic tweek, the cut-off frequency (f c1 ) and two other close frequencies (f 1 and f 2 ) near the cut-off frequency and their times of appearance (t 1 and t 2 ) in the spectrograms are estimated. The spectrograms produced by MATLAB codes are changed into bmp format files and are then uploaded into Microsoft Paint. The Paint software displays spectrograms on an identical coordinate system with arbitrary values from which coordinates corresponding to f and t are noted. In the second step, these coordinates are calibrated to produce the frequency and time using a Microsoft Excel spreadsheet. A program has been developed in Excel to convert these coordinates to the actual frequency and time. The accurate cut-off frequency of higher harmonics (f cn ) is obtained using the above steps and extrapolating up to the time of arrival (t a ) of the cut-off frequency of the first harmonic (see Fig. 1b). From this graphical analysis the resolutions in frequency and time components are 35 Hz and 1 ms, respectively.

Theoretical considerations
The ELF-VLF radio waves propagate by multiple reflections in the EIWG. The electromagnetic field in the waveguide can be decomposed into a sequence of independent field structures (modes) which propagate with different group velocities. Each of the modes is defined by its cut-off frequency (f cn ). For a waveguide having perfectly conducting boundaries, f cn of the nth mode is given by (Budden, 1961;Yamashita, 1978): where c is the velocity of light in free space and h is the height of the waveguide. If f >f cn and close to f cn , the mode propagates with a group velocity v gn given by (Ohtsu, 1960): The v gn approaches zero as f approaches f cn . No propagation occurs if f <f cn , as the wave is strongly attenuated. Using the Eq. (2), total distance d, propagated by the tweek of the nth mode with a perfect conducting, the EIWG can be obtained as (Prasad, 1981): where δt=t 2 −t 1 is the difference in arrival times of the two frequencies, f 2 and f 1 , close to f cn of the tweeks of any mode; v gf 1 and v gf 2 are the corresponding group velocities of the waves centered at frequencies f 1 and f 2 .
The main propagation characteristics of the modes emitted by a vertical dipole in the waveguide with spherical walls of finite conductivities are described by a factor e i k S n d (Wait, 1957), where k is the wave number and S n is the "sin" of the complex angle between the wave vector and the normal at the points of incidences in the EIWG (Wait, 1970, pp. 147-153) which is given by: The value of is governed by the electrical conductivities of the Earth's surface and the lower ionosphere. for the ground ( gi )/sea ( gi )-ionosphere paths is given by: where N is the refractive index of the medium inside the waveguide and N i , N g and N s are the complex refractive indices of the lower layer of the ionosphere, ground, and sea given by: where ω is the angular frequency of tweeks. ε i , ε g , ε s , and ε o are the permittivity of the lower ionosphere, ground, sea, and the free space, respectively. σ i , σ g , and σ s are conductivities of the lower ionosphere, ground, and sea, respectively.
The conductivity and permittivity of the ground and sea are well known. The major problem lies in the assignment of the parameters describing the ionospheric conductivity. The distribution of plasma in the ionosphere depends in a complicated way on latitude, solar zenith angle, season, and solar activity, etc. One of the simplest ionospheric profiles is an exponential variation of conductivity with height. The height dependent conductivity parameter ω r (h), for the Dregion ionosphere, is given by Wait and Spies (1964): (7) where ω o (h) is the electron plasma frequency and ν (h) is the effective electron-neutral collision frequency, both being functions of height h. The ω r (h) varies exponentially with h, at a rate determined by H , the reference height, and β the sharpness factor measured in km −1 . The frequency dependent conductivity of the lower ionosphere is given by If the wave frequency ω is less than ν (h) (∼10 6 Hz at 83 km), which is the present case, the conductivity can be written as: The numerical value of the attenuation rates (α n ) defined in dB per 1000 km path length, is given by (Wait, 1957;Hayakawa et al., 1995): The imaginary part of S n is −( /kh)/Re(S n ), so that where Re(S n ) is real part of S n .

Morphological characteristics of tweeks
The tweeks were recorded during September 2003 to January 2004 both in the local day and the night for durations of 30 min at the beginning of different hours. Analysis of the data showed that tweeks occur only in the nighttime (18:00-06:00 LT). Therefore, from February-July 2004, data were recorded in the nighttime only. The occurrence pattern of tweeks has been determined from the measurements made during September to December 2003 in the nighttime. A total of 2428 tweeks were observed from the recordings in this period in the night. The occurrence of tweeks has been classified into two categories-pre-midnight (18:00-00:00 LT) and post-midnight (00:00-06:00 LT). Table 1 shows the occurrence pattern of tweeks. Tweek sferics occur throughout the night, but there are more in the post-midnight period. Further higher harmonic tweeks (n≥3) occur more often in the postmidnight period than in the pre-midnight period. Figure 1a-e presents typical spectrograms showing the higher harmonic tweeks. Tweeks with second and third harmonics, as shown in panels (a) and (b), have stronger dispersion and longer duration than the tweeks with higher harmonics (panels ce). Most (>90%) of the tweeks are associated with a quasitransverse electromagnetic (QTEM) mode, as shown in the panels (b) and (c). The duration of the tweeks presented here is estimated to be in the range of 15-60 ms. Reznikov et al. (1993) observed the duration of typical tweeks to be in the range 40-50 ms, but may reach to about 100 ms. The distance of propagation, in general, is also larger for lower harmonic tweeks (not shown). It can also be noted that the duration of the dispersed part of any tweek decreases with an increase in harmonic number. Wave propagation in the VLF range is analyzed in terms of waveguide modes, characterized as quasi-transverse electric (QTE) and quasi-transverse magnetic (QTM) modes (Budden, 1962), with a cut-off frequency (f c ) of about 1.8 kHz. The single mode which has no cut-off frequency and propagates in the EIWG below 1.8 kHz is called the QTEM mode. The attenuation of the QTEM mode increases exponentially with frequency, and most of the QTEM frequency components are strongly attenuated above ∼1 kHz and hence do not overlap the QTE and QTM mode frequency components (Sukhorukov and Stubbe, 1997). The QTM mode is similar to the TM mode except that it has a small longitudinal (in the direction of wave propagation) magnetic field component that cannot exist significantly over long distances. For frequencies less than 15 kHz the lower order QTM modes approximate pure TM modes more than the higher order modes (Wood, 2004). We therefore treat the modes in the EIWG as pure TM modes above f c of the first order mode (∼1.8 kHz) and below that a single QTEM mode which has no f c . Thus the propagation of energy below the cut-off frequency of the first order TM mode, as shown in the spectrograms (Fig. 1), must have occurred as QTEM mode waves.
The tweek sferics, as shown in Fig. 1, having ELF (<1.5 kHz) frequency components, suggest that these might be produced by strong lightning discharges, possibly positive cloud-ground (+CG) lightning discharges which make up roughly 10% of all lightning strokes (Uman, 2001). About 50% of ELF sferics with slow tail are associated with sprites (Rodger, 1999). Strong +CG discharges with large return stroke peak currents can be associated with ELF radio atmospherics observed at large distances (∼1500 km) from the discharge (Sukhorukov and Stubbe, 1997;Cummer and Inan, 1997;Pasko et al., 1998;Ohkubo et al., 2005).

Attenuation of sferics propagating in a partially conducting waveguide
To explain the occurrence of tweeks at this station, the attenuation rate (α n ) is calculated for modes n=1-3 using Eq. (11). In this equation, Re(S n ) is a function of the mode number, frequency, and the refractive index of the reflecting boundaries of the waveguide. Equations (5) and (6) clearly show that for a perfectly conducting EIWG, will be zero and α n for each mode would then vanish and higher harmonic tweeks should be observed at all times. In practice, tweeks are not observed always and tweeks with higher harmonics occur less often (Kumar et al., 1994;Singh and Singh, 1996). The ground, sea, and ionospheric boundaries are good conductors at VLF but not perfect conductors. The propagation path of the tweeks to this station is mostly over sea and its larger conductivity than ground offers less attenuation, allowing higher harmonic tweeks to occur more often. In the analysis of attenuation, the values of the parameters are taken from (Prasad, 1981;Westerlund, 1974): ionospheric dielectric constant ε i ε o =0.5, ground dielectric constant ε g ε o =20, sea dielectric constant ε s ε o =81, ground conductivity σ g =10 −4 S/m, and sea conductivity σ s =5 S/m. Daytime VLF reflection heights are reported to be around 70-75 km (Glukhov et al., 1992) and our observations show the nighttime tweek reflection heights mostly in the range of 85-90 km. Hence, for calculation purposes, the ionospheric reflection heights h, for daytime, were chosen as 71 km and 74 km and for nighttime as 86 km and 90 km.   The lower ionosphere can be characterized as the "Wait ionosphere", defined by a reference height, H in km, and the exponential sharpness factor, β in km −1 . The ionospheric conductivity, σ i , has been calculated using the height dependent conductivity parameter ω r (h) given in Eq. (7). For the day and the nighttime conditions H and β are taken to be 70 km, 0.3 km −1 and 80 km, 0.5 km −1 , respectively. The magnitude of ω r (h) calculated using the Eq. (7) is 3.4×10 5 s −1 for the daytime and 3.7×10 7 s −1 for the nighttime. Figure 2a shows a plot of α n versus f calculated using Eq. (11) for the EIWG with ground as the lower surface. It can be seen from Fig. 2a that the attenuation increases sharply as the frequency approaches the cut-off frequency. The attenuation for all the modes also increases when the reflecting layer height is lowered. Figure 2b shows a plot of α n versus f for the waveguide formed by the sea surface and the lower ionosphere. Figure 2 (a and b) shows that for the waveguide considered here the attenuation is larger in the daytime (h=71 km and 74 km) than that in the nighttime (h=86 km and 90 km) and increases sharply as the frequency approaches the cut-off frequency. The geographical location of our station on a small island, Fiji, in the South Pacific region, is such that the tweek sferics propagate a significant portion in the waveguide formed by sea and the lower ionosphere. It can be confirmed from panels (a) and (b) of Fig. 2 that ground-ionosphere path offers larger attenuation by about 3-5 dB/1000 km than the sea-ionosphere path for all modes. This could be the main reason for the more probable occurrence of tweeks with higher harmonics at this station compared to those reported where ground dominates      as the lower wall of the waveguide in the total propagation path. The attenuation is larger for higher harmonics tweeks, therefore, a higher harmonic of tweeks are not observed frequently.

Tweek reflection heights
Using Eq. (1), the ionospheric reflection heights h are calculated for the tweeks shown in Fig. l (a-e) and are presented in Table 2. The value of h varies in the range 83-92 km. Based on the clarity of dispersion seen in the spectrogram, the analysis for h was made using 503 selected tweeks observed during September 2003-July 2004. This selection included a total of 320 tweeks during September-December 2003 and 183 tweeks from January-July 2004. The mean cut-off frequency f cm for the different modes are first calculated using the cut-off frequencies estimated from the spectrograms and the number of sferics considered for each mode. The variation of f cm with mode number n is shown by the solid line (left hand scale) in Fig. 3. The frequency obtained by dividing f cm by n, is defined as the "mean fundamental frequency". The variation of f cm /n with mode number n is also plotted in Fig. 3  This indicates that the cut-off frequencies of higher harmonics (n>1) are not exact multiples of the cut-off frequency for the first harmonic (f c1 ) and n. They are slightly less than f c1 ×n, as shown by the trend line. The negative gradient of f cm /n indicates an increase in the electron density with height in the lower ionosphere. It also indicates that the higher order tweeks penetrate slightly deeper into the lower ionosphere. This result supports that the ionosphere is not    sharply bounded and homogeneous. The cut-off frequency depends on the reflection height, which, in turn, depends on the electron density. The mean cut-off frequency of the different modes can be used to determine the reflection heights of different harmonics. The calculated mean reflection height h m of modes n=1-6 is shown in Fig. 4; h m increases with increase in n and varies in the range 83.4-85.6 km for n=1-3. The number of tweeks with n=4-6 are not significant to consider the variation of h m with n in this range. Shvets and Hayakawa (1998) estimated the steepness of the electron density profile from the fundamental cut-off frequencies of multimode tweek sferics. Their results showed a negative gradient in f cm /n, indicating an increase in the electron density from 28-224 cm −3 in the altitude range of 2 km at the reflection height of about 81 km. Thus, f cm is also a useful parameter for investigating electron density variations in the nighttime D-region ionosphere (e.g. Ohya et al., 2003Ohya et al., , 2006.

Ionospheric reflection height variation with time, location, and solar cycle
The analysis shows that the mean cut-off frequencies of the different harmonics vary with the time of observation.
The temporal variations of f cm and h m at Suva are shown in Fig. 5. The plots have been made only for n=1-3 harmonic tweek sferics due to the small number of tweeks observed with higher harmonics. For the first harmonic, h m obtained from f cm increases up to 23:00 LT and then decreases with a minimum at around 01:00 LT and then gradually increases up to 05:00 LT. The h m obtained from second and third harmonics shows similar variation with maximum at around 04:00 LT. Using a small number of tweeks, Kumar et al. (1994) estimated the h to vary in the range of 83-89 km from the tweek observations taken at a low-latitude station (23.1 • N) in India during the period January-March 1992. The average sun spot number (R z ) for the period  January-March 1992 was about 110. From the tweeks observations carried out in the Southern Hemisphere in the Indian and Atlantic Ocean onboard the research vessel "Academician Vernadsky" during the period January-April 1991 Shvets and Hayakawa (1998) estimated the h m to vary from 81.5-83.3 km for n=1-8. The R z for the period of January-April 1991 was about 140. The ionospheric reflection heights in this work during the period September 2003-July 2004, which fall in the lower (moderate) solar activity period (R z ∼45), are higher by 1-1.5 km than those estimated by Shvets and Hayakawa (1998), indicating the spatial and solar activity dependence of reflection heights. McRae and Thomson (2004), using long wave propagation capability (LWPC) codes to interpret the VLF amplitude and phase of navigational transmitter signals, showed that the ionospheric VLF reflection height is somewhat lower (by about ≤1 km) during solar maximum than during solar minimum. Higher hm values estimated in the present study as compared to those reported by Kumar et al. (1994) and Shvets and Hayakawa (1998) also indicate a decrease in the electron density in the nighttime D-region ionosphere during low solar activity. The dependence of the D-region electron concentration on solar activity is not very clear, in spite of the fact that during recent decades this problem has been attempted by many researchers (see Danilov, 1998, and references therein). Some authors reported an increase in the D-region electron concentration while others found a decrease in the concentration with an increase in solar activity. It is also widely known that the D-region electron concentration measurement meets serious technical difficulties; however, the ELF/VLF method seems to be more reliable at least for nighttime D-region studies. Bremer and Singer (1977) studied the variation of D-  region electron density with solar activity using radio wave propagation in a broad range (10-2500 kHz) of frequencies.
They found that the solar activity effect on D-region electron density to be small but positive at all the altitudes, which seems to support our results. Detailed discussion on solar activity effects in the ionospheric D-region is beyond the scope of this paper; the interested reader is advised to see the paper by Danilov (1998).

Tweek propagation distances
Using Eq. (3), the propagation distances d are calculated for the tweeks shown in Fig. 1a-e and are given in Table 2, which vary from 1000-4500 km. Figure 6 shows the range of d traveled by the different tweeks from the overall analysis of 503 tweeks. The maximum propagation distance calculated from a tweek was about 12 000 km. About 90% of tweeks propagated distances in the range of 1000-5000 km, with maximum in the range of 2000-3000 km. Using a similar method, Kumar et al. (1994) found d to vary from 1500 km to 2500 km. This may be because of the different path conditions in the EIWG where the lower boundary of the waveguide is ground for the Indian observations but largely sea for the observations made here. Christian et al. (2003) found about 78% of lightning on the Earth to occur between the ±30 • around the geographic equator, with a highest rate at the west coast of Australia. It can be said here that most of the tweeks observed at Suva (18.2 • S) probably came from over the sea at the low latitudes (±30 • around the geographic equator) and a small number from other locations, most likely from west coast of Australia. However, a direction finding method can determine the locations of the lightnings associated with tweeks, which is an area of further experimental and theoretical research interest to us.

Summary and conclusions
The lightning generated ELF-VLF tweek sferics observed with the WWLLN VLF system during the period from September 2003 to July 2004, have been analyses and the results are presented. The VLF recording and analysis system employed in this work is relatively new and advanced, which provides precise timing measurements using GPS and highresolution spectrograms using MATLAB code. The ELF-VLF radio wave technique has demonstrated the potential to monitor changes in the bottom part of the lower ionosphere.
Tweek sferics have been utilized to determine: 1) the nighttime variation of ionospheric reflection heights, and 2) their propagation distance in the EIWG. The main findings of this study are concluded as follows: -Tweek sferics occur during the local night between 18:00-06:00 LT. The higher harmonics occur more often in the post-midnight period than in the pre-midnight period. Tweeks up to the 6th harmonic were observed. The occurrence of tweeks with higher harmonic decreases as the harmonic number increases. The dispersed portion of the tweeks reduces as harmonic number increases. Higher harmonic of the same tweek is reflected from about 2.0 km higher than the lower harmonic.
-The waveguide formed by the ground and the lower ionosphere offers more attenuation than that formed by the sea and the lower ionosphere. The attenuation is less at night, which may be the main reason why tweek sferics occur only at night.
-Most of the tweeks are found to have ELF frequency components. The propagation of ELF waves (below cutoff frequency of first order) must occur by the QTEM mode.
-The ionospheric reflection heights calculated from tweeks observed during lower solar activity period of present study are found to vary in the range 83-92 km, which are higher by about 1-1.5 km than those estimated during the high solar activity period.
-About 90% of the tweeks propagated distances between 1000 and 5000 km in the waveguide to the receiving station.