- Email: [email protected]
Error sources and travel time residuals in plasmaspheric whistler interpretation
Journal of Atmospheric and Terresrrial Physics. Vol. 51, No. 4. pp. 24’%258, 1989. Printed in Great Britain.
0021-9169/X9 $3.00+ .I0 Pergamon Press plc
Error sources and travel time residuals in plasmaspheric whistler interpretation GY. TARCSAI Geophysical Department, Eotvos University, Budapest, Kun B. ter 2, Hungary, H-1083 H. J. STRANGEWAYS Department of Electrical and Electronic Engineering, The University of Leeds, Leeds LS2 9JT, U.K. and
M. J. RYCROFT British Antarctic Survey, NERC, Madingley Road, Cambridge CB3 OET, U.K. (Received in final form I November 1988)
Abstract-In the interpretation of observed whistlers by curve fitting, systematic travel time residuals appeared which were studied by extensive simulations using ray-tracing, numerical integration and curve fitting. The residuals were found to originate from the commonly used approximations in the refractive index and ray path of whistler mode waves, which result in travel time increments or decrements, not accounted for in whistler interpretation. These approximations and the assumed form of the electron density distribution also lead to systematic errors in the diagnostics of plasmaspheric electron density by whistlers. In addition, the effects of other error sources, including random measurement errors, are also reviewed briefly. It is shown that the fine structure of residual curves is connected to propagation conditions. Thus, their study may yield a new research tool for studying whistler trapping, ducting structures and other features of whistler propagation. The application of residual analysis in conjunction with digital matched filtering of whistlers seems to be especially promising for further whistler studies.
1. INTRODUCTION
by numerical integration, ray-tracing and curvefitting. Emphasis is laid on the effect of these approximations on whistler travel times and on whistler and magnetospheric parameters obtained by curve fitting. In addition, the effect of the assumed form of electron density distribution on curve fitting results is also studied.
It has been previously shown (TARCSAI, 1975) that in the routine interpretation of whistlers by curve fitting, systematic travel time residuals appear. It has been suggested that these residuals can be attributed to various approximations common to routinely used whistler analysis techniques. The residuals were successfully reproduced later by curve fitting to theoretical whistlers computed by ray tracing (STRANGEWAYS, 1978) and by numerical integration (B-R, 1979), and found repeatedly on measured data (TARCSAI, 1981; STRANGEWAYSet al., 1983 ; TARCSAI
et al., 1987). It has been found that the observed systematic residuals appear as a result of three commonly used approximations. These are the omission of a ‘+ 1’ term in, and the contribution of ions to, the plasma refractive index, and also the assumption of purely longitudinal propagation rather than the real ray path, which executes snake-like excursions back and forth across the ducting structure. The present work investigates the importance of these approximations 249
2.
COMPUTATIONOF WHISTLER TRAVELTIMES
2.1. Ray-tracing In order to study the effect of approximations, theoretical frequency-time profiles for whistlers were computed by numerical integration and ray-tracing. The ray-tracing calculations were performed using a computer program developed in its original form by P. D. Alexander (ALEXANDER,1971) ; which uses the full refractive index p of a cold collisionless plasma containing both electrons and ions as given by STIX (1962) (in a number of runs the refractive index was determined for electrons only in order, by comparison, to determine the additional effect of ions on the
250
Gu. TARCSAIel al.
travel times). The travel times were computed for 1.516 kHz rays (in steps of 0.5 or 1.O kHz) assuming field aligned ducts of enhanced electron density centred on nine different L-values in the range L = l&4.5. The highest frequency computed for a given whistler was 16 kHz or, if lower, about half the minimum equatorial electron gyrofrequency along the given field line. The rays were started in the ducts with initial wave normal angles of O”, IO”, +20” and 30” to the geomagnetic field. A starting altitude of 1000 km was employed because the curve fitting procedure of TARCSAI (1975, 1977) was developed for a 1000 km1000 km altitude propagation path. The calculations were performed for six diffusive equilibrium models of plasmaspheric plasma distribution. These are the WN and SD models (ALEXANDER, 1971; DENBY et a/., 1980) and the DE-l-DE-4 models given by PARK (1972). In these models electrons and hydrogen, helium and oxygen ions are distributed according to diffusive equilibrium theory (ANGERAMIand THOMAS, 1964). Although in the majority of computations isothermal field lines were assumed, models including the electron temperature variation along geomagnetic field lines were also used in a modified ray-tracing model (STRANGEWAYS,1984, 1986). The ducts were modelled by Gaussian density enhancements superimposed on the background electron density (STRANGEWAYS and RYCROFT, 1980). The maximum fractional density enhancement 6 (at the centre of the duct) was taken to be 8, 10, 15, 20, 60 and 100%. The effective width of the ducts in the equatorial plane (where the enhancement is reduced to a tenth of its value at the duct centre) was chosen to be 430 km, although other ducts widths were sometimes used. To model the Earth’s magnetic field the same centred dipole approximation (PARK, 1972) was used in ray-tracing as in the numerical integration and curve fitting.
index p reduces to :
mode waves, the phase refractive p2=l+c
ffi 1
f(fHj
+
(3)
c,f) ’
where the summation is for the electrons and all the ionic species. f, are the plasma frequencies, fH,the gyrofrequencies and ci the sign of the charge of the electrons and all the ionic species which constitute the medium. In all of the methods used for whistler analysis the term ‘1’ and the contribution of ions on the right side of (3) are neglected. The plasmaspheric travel times were taken as :
s
1
a-) = 7 2cf
fbf"
(4)
(fH-fyds,
path
where f.and f,, are the plasma frequency and gyrofrequency for electrons, respectively, and the integral is taken for the path above 1000 km altitude. A series of f(t) pairs were computed for each simulated whistler by computing (1) numerically with p’ obtained from the full form of (3) or by neglecting the ‘l’, the ionic effects or both. In order to facilitate comparisons with the ray-tracing results, the plasma densities along the field lines of propagation were adjusted to be the same as those at the centre of the ducts. 2.3. Effect of approximations
on travel times
Comparison of the t(f) values obtained by raytracing and numerical integration for various cases
-
1
effectof
ims
2.2. Numerical integration Assuming strictly field aligned propagation, the whistler travel time, t, at frequency f can be computed by:
t(f) = i
c
s
/id&s,
(1)
path
a For exactly field aligned propagation
-30 g
I
w -40
where c is the speed of light in vacuum, & is the path element along the geomagnetic field line and p’ is the group refractive index given by :
P’ = q w-f).
---. ‘6
-50
j
/
I
,.A
effect of
and for whistler-
my path
“f_2._._._
/
2
I
4
I
6
I
6 Freqwncy
(2)
ducted
I
IO
I
12
I
14
I
16
(kHz)
Fig. 1. The effect of commonly made approximations on whistler travel times at L = 1.6, 2.4 and 3.5 as a function of frequency for the DE-I model.
251
Plasmaspheric whistler interpretation
2
I 2
I 4
I 6
I 6 Frequency
I IO
I 12
I 14
I 16
(kHz)
Fig. 2. The combined effect of the various approximations shown in Fig. 1 on whistler travel times.
permitted the separation of the effects of different approximations. These results are shown in Fig. 1 for the DE-l model and for three whistlers propagating at L = 1.6,2.4 and 3.5, respectively. The effects of the terms omitted from the refractive index were obtained by comparison of travel times (mainly obtained by numerical integration) with and without the respective terms included. Similar calculations by ray-tracing gave very similar results to numerical integration. The effect of the curved and snake-like ducted path was obtained by comparison of travel times computed by numerical integration with the wave normal always parallel to the geomagnetic field direction, and by raytracing, using the full refractive index formula (HINES, 1957) and Gaussian cross-section ducts. It is seen that the effect of ions is to increase travel times, especially at lower frequencies, as was shown by STOREY(1956) and BARRINGTONand NISHIZAKI (1960). The term ’ 1’ in equation (3) also leads to travel time increases of the order of several ms which are proportional to the wave frequency. The total increases in travel times caused by these two factors are 2.7-3.9% at L = 1.6 ; 1.6-2.2% at L = 2.4 and 0.7-1.4% at L = 3.5. On the other hand, the calculations show that along the true snake-like ray paths, the waves travel faster than they do along the axis of ducts in the purely longitudinal mode. This results in travel time decreases as can be seen in the lower part of Fig. 1, which effectively compensate for the effects of the approximations made in the refractive index. Figure 2 shows the combined effect on travel time of the approximations mentioned above, i.e. the difference between the exact travel times computed by ray-tracing using the full refractive index formula and those obtained by numerical integration for field aligned propagation with the approximate refractive index formulae [equation (4)]. It follows from Fig. 2 that the approximate
theoretical travel times used in routine whistler analysis may be lower (for lower L-values) or higher (for L-values greater than about 3) than the true values because of the commonly used approximations. The neglected effects in travel times amount to 2.1-3.6% ; 0.7-1.4% and -(O&1.3)%, at L = 1.6, 2.4 and 3.5, respectively, depending on the wave frequency. These figures are rather similar for different plasma distribution models and ducts and are practically independent of the initial wave normal angle. The consequence of these approximations upon whistler analysis is two-fold. Firstly, part of the neglected travel time increments is absorbed into the procedure of analysis and results in errors in the computed whistler and magnetospheric parameters. Secondly, the non-absorbed contributions appear as systematic travel time residuals in curve fitting which cannot be eliminated by any changes in the fitted whistler parameters. These points are treated in the following sections.
3. CURVE FITTING TO CALCULATEDWHISTLERS 3.1. The FIT procedure The computed whistlers, given by 8-16 f-t pairs, were analysed by a curve fitting method and computer code called FIT (TARCSAI, 1975, 1977). This method proved to be a reliable tool for the analysis of nose or non-nose whistlers. It is based on BERNARD’S (1973) formula for the dispersion o(f) = l.f’/2, which approximates the travel time integral (4) by an analytic expression. The method uses a least squares estimation of the zero frequency dispersion D,, the equatorial electron gyrofrequency fHE in Bernard’s approximation, and a third parameter, T, locating the time of the causative ‘spheric’. The FIT procedure minimises the sum of the squares of the differences between the observed and calculated propagation times. After obtaining these three parameters, the Lvalue of propagation, the electron density at 1000 km (n,) and at the equator (n,,) and the tube electron content (NT) were also calculated using the quasiconstants given by PARK (1972) and extended to lower latitudes by B~~TTNER(1974). The quasi-constants were determined by model calculations using equation (4) for the travel time, and assuming the DE-l model of plasma distribution (for routine whistler analysis this model is used within the plasmasphere). Together with whistler and plasma parameters, the FIT program also yields the estimated standard deviations and the final travel time residuals defined as the travel time of a whistler component at any frequency minus
252
GY. TARC~AI et
al.
the computed travel time at the same frequency calculated with the estimated whistler parameters.
Systematic errors in computed pommeters
:_?J--
3.2. Systematic errors in curve-fitting In order to check the internal accuracy of the FIT procedure, longitudinal whistlers were first analysed which had been computed from equation (4) and with model DE-l (TARCSAI,1981). The calculations demonstrated that the systematic errors in the estimated D0 and fHE (and also in nose frequency f, and the corresponding travel time t,,, which are derived from Do and fHE) are less than 2-3% in the range L = 1.6 5.0 (for L > 3 the bias in fn remains below 0.1%). The time shift between the true time origin and the estimated time of the originating ‘sferic’, T, was less than 2-3% of t, or for L < 2.5 less than 4 ms. The error in L, n,, neq and N, turned out to be less than 1, 3, 5 and l%, respectively, and for L 3 2.5 they remained below 1%. (Note that this is after constant corrections of 4.5%, 2% and 4% in n,, neq and N, respectively, were introduced to account for the approximations in Bernard’s formula and in the quasiconstants.) In conclusion, the internal accuracy of the FIT method is satisfactory and seems to be superior to that of other whistler analysis methods (SMITH et al., 1975 ; CORCUFF, 1977 ; STUART, 1977a,b). It should be noted, however, that it might be worthwhile replacing Bernard’s formula with the more accurate dispersion formula derived by DANIELL (1986) to help to obtain a value of neq that was not biased by an assumed form of electron density distribution along the field lines. The effects of the approximations discussed in Part 2 have been studied by analysing a great number of ray-traced whistlers by the FIT method. The results are summarised in Fig. 3 in the form of error bands, which represent the possible overall systematic errors corresponding to the commonly used approximations in the refractive index and ray path and resulting also from the assumed (fixed) form of electron density distribution and the small internal errors of the FIT procedure. Because the six plasma distributions used in the ray-tracing represent a rather broad range of plasmaspheric conditions representative of the variations under real circumstances, the obtained error bars are fairly realistic. Furthermore, the results obtained are characteristic (although perhaps an underestimate) of errors in other methods of whistler analysis, which are influenced by the same (and also, as a rule, by certain other) approximations as used in FIT. As may be seen in Fig. 3, for L 3 2 the FIT procedure gives the L-value of propagation of whistlers and the tube electron content, N, with errors of
I I.6
I
I
I
2
2.4
2.6
I 32
I
I
I
36
4
4.4
L- value
Fig. 3. Bands of possible systematic errors in the results of whistler interpretation, appearing as a consequence of common approximations in the refractive index, ray path and plasma density distribution (top, L-value; next equatorial electron density ; next total electron content along the geomagnetic field line; next electron density at 1000 km altitude ; bottom, time of causative ‘sferic’). several per cent. The magnitude of Arzeqvalues is less than f20% for 2.5 < L < 3.5. At higher L-values neq tends to be overestimated (eventually by as much as 30%) while below L = 2.2 it is strongly underestimated. The electron density at 1000 km will generally be strongly overestimated at L > 2.3, sometimes by as much as 120%, while in the range L = 1.6-2.3, An, varies in the range + 60%. Lastly, between L = 3 and 4.5, the time shift of the causative ‘sferic’ T may have a systematic error of the order of + 20 ms depending on the actual plasma distribution, with a tendency for T to be underestimated above about L = 3.5. Below L = 3.2, ATwill always be positive indicating that the FIT method places the causative ‘sferic’ earlier than it actually occurs. 4. DISCUSSION The behaviour of AT is the direct consequence of the neglected travel time increments, caused by the
Plasmaspheric whistler interpretation approximations in the ray path and refractive index. These approximations are responsible for the general trends in the error bands shown in Fig. 3, while the width of the error curves at a given L-value is caused by the difference between the fixed electron density distribution (DE-l model) and those distributions which may actually occur [as also indicated by ROTH (1975)]. The relative importance of these error sources depends on the L-value and on which parameter is being considered. The overestimation of the L-value and the strong underestimation of equatorial plane electron densities at low L-values is caused partly by neglect of the sinuous ray path but mainly by approximations in the refractive index. At higher latitudes, however, the unknown shape of the field aligned plasma density distribution is the main source of uncertainties in neq, NT and n,. It is worth noting here that the inversion of the integral equation (4) for electron densities along the field line (TARCSAIand DANIELL,1979) is completely free of any assumptions about the shape of the electron density distribution and is remarkably stable for noisy, measured f-t data. Its successful application is hindered, however, by the approximations in the refractive index and ray path discussed above which are very difficult to avoid. With respect to the errors shown in Fig. 3, it is important to note that at low L-values the errors in L and in neq are correlated, i.e. an overestimation of L is accompanied by an underestimation of neq. This means that the points on the neq vs L curve tend to remain on the profile, i.e. the difference between the true n,,(L’) value and the estimated n&(L’) value both of which are considered at the biased estimated L’ is smaller than the neq bias shown in Fig. 3 at the true (unknown) L or even at L’. This favourable circumstance only holds below about L = 2.3. At Lvalues above 2.3, AL and An,, are rather weakly (and, for various plasma distribution, oppositely) correlated. It is to be noted that the simulated whistlers used in this study covered the full range of frequencies between 1.5 kHz and f,,J2 or 16 kHz (if fHe> 32 kHz). The absence of higher frequency components in the whistler traces at low and mid-latitudes decreases the reliability of nose-extension methods and magnifies the systematic errors of the estimated parameters. For example, at L = 2 the reduction of the highest frequency from 16 kHz to 12.5 or 10 kHz increases AL from 0.08 to 0.12 and 0.16, respectively, with a parallel decrease in neq. This increase in AL counteracts the relation between AL and An,, discussed above. However, at higher L-values, e.g. at L = 3.5, the omission of the whistler components
253
above the nose frequency does not influence the results of curve fitting appreciably. In discussing systematic errors in curve-fitting, there are three additional factors which should be considered. The first is the approximate character of the centred dipole field used in routine whistler analysis, the effects of which have been studied by a number of researchers (LIKHTER and MOLCHANOV,1968; MATHIJRet al., 1972 ; SAGREDOand BULLOUGH,1972). These studies concentrated mainly on ring current effects which are important under disturbed conditions and outside the plasmasphere. For the present work, which refers to the study of whistlers inside the plasmasphere, it is pertinent to refer to the results of SEELY (1977). He constructed a hybrid quiet-time magnetospheric magnetic field model that combined the IGRF spherical harmonic series representation of the main field with the contributions from ring, tail and magnetopause current systems (OLSON and PFITZER, 1974). He showed that using the dipole model in whistler analysis results in L and NT values being overestimated by 24% and less than 15%, respectively, below L = 4.5 ; neqand n, values are correspondingly underestimated by 5-17% and l-15%, respectively. These errors, which are to be added to those discussed above, increase somewhat the systematic error in neq and n, at L-values higher than about 2.2-2.5. The plasma temperature should also be considered as it may influence the results of whistler analysis in a two-fold manner. Firstly, considering the existence of high thermal plasma temperatures in the magnetosphere, it should be determined whether it is appropriate to use the cold plasma approximation for the refractive index in whistler analysis. Although thermal effects may become significant close to the cyclotron resonance (SCARF, 1962; GUTHART, 1965) or beyond the plasmapause (MOREIRA,1982), for the usual whistler frequencies and inside the plasmasphere, the use of warm plasma theory (GUTHART, 1965 ; SAZHIN and SAZHINA,1982) in computing whistler ray paths or travel times produces insignificant changes compared to the cold plasma approximation (GUTHART, 1965 ; HASHIMOTO et aI., 1977 ; STRANGEWAYS, 1986) and, therefore, does not have any noticeable effect on whistler analysis (SEXLY,1977). Secondly, field aligned temperature gradients modify significantly the plasma distribution, a fact neglected in routine whistler analysis where isothermal distributions are used. SEELY(1977) has shown that, in general, a temperature increasing monotonically up the field line, and having its maximum gradient near 1000 km altitude, does not significantly change whistler values of equatorial parameters (L, n,,) with
GY. TARCSAI et al.
254
respect to values obtained using a constant temperature model. However, 1000 km altitude electron densities (n,) may be increased by a factor of 2 or more depending on the equatorial temperature. On the other hand, L and neq values inferred from whistlers are significantly affected by temperature gradients near the equatorial region. STRANGEWAYS (1984, 1986) conducted ray-tracing studies in diffusive equilibrium plasma density models including realistic temperature variations along geomagnetic field lines. He has shown that the equatorial electron density at L - 4.5 can be overestimated by as much as a factor of 2 if isothermal field lines are assumed in curve fitting when, in actual fact, a temperature gradient of the order of 1.6 K km- ’ at 1000 km (and 0.8 K km- ’ at 2000 km) exists. However, contrary to these results, CORCUFF et al. (1972), CARPENTERand CHAPPEL.L (1973), CARPENTERet al. (1981) and TARCSAI(1981) reported good agreement between in situ satellite observations of neq and determinations from curve fitting to whistlers assuming diffusive equilibrium along isothermal field lines. Therefore this question should be studied further, including the rigorous consideration of asymmetric conditions in the conjugate hemispheres (cf. STRANGEWAYS, 1982a). A third factor, leading to systematic errors in whistler analysis and not treated here in detail, is the travel time below 1000 km altitude, which is best subtracted from measured time delays before a curve-fitting analysis. This has two components. The first, the subionospheric propagation time, which may amount to several ms, is absorbed in the computed time-shift parameter T and may result in l-2% errors in the electron densities. The travel time in the ionosphere below 1000 km is generally computed by the formula given by PARK (1972) : At(f) = y, where At is given in s, the critical frequency of the F2layer, averaged for both hemispheres,f,F2, is in MHz and the wave frequency is in Hz. Ar may amount to several tens or hundreds of ms. Our experience and model computations show that an overestimate of f,F2 by 0.5 MHz results in about 0.5, -6 and -3% errors at mid-latitudes in L, neq and NT, respectively. These figures are lower at greater L-values and may be much higher at low L-values where ionospheric and magnetospheric travel times bcomes commensurate. In addition to systematic errors, the results of whistler analysis also contain unpredictable random errors. Figure 4 shows the variation of the median
\
I
errars In computed porarr&ef.s
Fkmdm
I
I
I
I
I
I
I
I
I5
2
2.5
3
35
4
4.5
5
L-value Fig. 4. Median values of the standard deviations of L, neq, n, and Nr as a function of the L-value of propagation. The curves were obtained from curve fitting to about 1500 measured whistlers and are normalised to 20 ,f-f pairs per trace.
value of normalised relative standard deviations of L, neq, n, and NT obtained from whistler curve fitting, as
a function of the L-value of propagation (TARCSAI, 1981). These curves were obtained from data consisting of more than 1200 mid-latitude whistlers (TARCSAI,197.5; TARCSAIet al., 1988) and 300 high latitude whistlers scaled from sonagrams and processed by the FIT method. The normalisation of the standard deviations was performed by reducing them to 20 points (f-t pairs) per whistler trace. Simulations on computed whistlers supported fully the characteristics of error propagation found from the measured data. Figure 4 demonstrates that, for most cases, the standard deviations of the L-value obtained and of the plasma densities remain below 10%. Thus, the systematic errors shown in Fig. 3 and discussed above are comparable to, or sometimes much greater than, the more unpredictable errors associated with whistler curve fitting. 5. SYSTEMATTC TRAVELTIME RESIDUALS 5.1. Measured and computed residuals In curve fitting applications, the residuals represent a very useful by-product which permits the checking of the validity of the theoretical model fitted to actual data. When the residuals are distributed in a nonrandom manner as a function of an independent variable, it indicates the need for an improvement of the model. That is the case with whistler curve fitting (TARCSAI,1975), which led to the simulation studies reported here. The travel time residuals are usually less than the random errors when conventional f-t
Plasmaspheric whistler interpretation
,-
L- 27SD 2 ducts --L-2fSD3hcta
2
4
6
8
10
12
14
255
to a ray-traced whistler (WN model, L = 2.4, 10% enhancement duct with a width parameter of 430 km). Here, no effort has been made to maximise the fit by variation of the ray-tracing plasmaspheric model or duct parameters. Further analysis of simulated whistlers has shown that the greatest part of the residuals is caused by the omission of ionic effects from the refractive index. Whistlers computed with ionic effects gave residual curves with peak-to-peak amplitudes from 4 to 8 times larger than those with the effect of ions neglected for both ray-traced and purely longitudinal simulated whistlers. This supports previous findings by BARRINGTONand NISHIZAKI (1960) and BEGHIN (1966) who, by using special techniques, demonstrated increased dispersion due to ions at the lowest frequencies of several low latitude whistlers. The same study was also attempted by GUPTA and SINGH (1977) with not entirely conclusive results. The term ‘1’ and the curved ray-path lead to residuals of only several tenths of a ms. It is to be noted that the narrower the frequency range of whistlers the lower will be the amplitude of the corresponding residual curves, while an extension of whistler analysis to frequencies as low as 0.7-1.0 kHz results in peakto-peak residual amplitudes of IO-12 ms.
I6
Frequency (kHz)
Fig. 5. Travel time residuals vs frequency for ducted whistlers, computed by ray-tracing for various plasma distributions and at different L-vatues. In panel (b) a measured cnrve constructed from 40 whistlers is afso shown.
scaling techniques are used. In this case averaging at given frequencies for a larger number of whistlers propagating under similar conditions leads to an effective suppression of noise, and systematic residuals appear as shown in Fig. 5b. The overall residual pattern is of the form fl with increasing frequency. Departures from the mean have a characteristic run of signs (+, -, + , -). The peak-to-peak amplitude is generally several ms, sometimes with multiple peaks on the smooth p-shaped curve. The curve fitting to longitudinally propagated whistlers computed by equation (4) has shown that, for such whistlers, the magnitude of residuals remains below about 0.1 ms, i.e. Bernard’s approximation is fairly good. However, the analysis of ray-traced whistlers, computed with the full refractive index, produced residual curves which agree very well with the observed ones. This is seen in Fig. 5b, which shows average residuals obtained from 40 measured whistlers (L = 1.9-2.5) and the residual curve belonging
5.2. Possible applications of residuals The smooth theoretical residual curves, shown in Figs. 5a, 5b and SC for various L-values and plasma distributions, may be considered as being ‘normal’ residuals produced by the known common approximations in whistler curve fitting. These curves vary somewhat with the parameters of the ray-tracing but for normal ducted whistlers always show the same smooth pattern. On the other hand, measured residual curves tend to exhibit a fine structure superposed on the normal pattern, which gives additional information on the propagation conditions. Extra peaks appear rather frequently, especially near the low and high frequency limits of the residual curves, as seen at about 3 and at 14.5 kHz on the measured curve of Fig. 5b (other minor peaks in this curve are due to random measurement errors). According to ray-tracing and curve fitting simulations, the anomalous residuals may perhaps be explained by the following factors. (1) Trapping and untrapping of whistler waves depends on the frequency. It has been shown (STRANGEWAYS, 1978) that the larger the altitude range at which different frequency rays enter and leave the duct, the more complex will be the fine structure in the residual pattern. Thus the altitude of the bases of the ducts and the altitude range of trapping and loss
256
Cu. TARCSAIet al.
might effect the residual pattern of the whistlers trapped in them. Since BERNHARDTand PARK (1977) suggest that the altitude of the end of the ducts will vary between day and night, it would be interesting to investigate whether there is any diurnal variation in the residual pattern of whistlers. (2) It is rather probable that whistlers do not propagate in simple, ‘ideal’ ducts, described by Gaussian or other simple cross-sections, but in field aligned irregularities which have a more complicated structure, forming a conglomeration of ducts, spaced more or less closely to each other. It turns out (STRANGEWAYS, 1978) that the residual patterns for whistlers trapped in ducts containing fine structure will be more complicated than for those trapped in ducts with simple structure. This is illustrated in Fig. 5d where the two residual patterns were obtained for the superposition of two and three Gaussian ducts, respectively. The elemental ducts were characterised by 6 = 15% and a width of 430 km, and were separated from each other by AL = 0.015 (- 100 km at the equator). Both curves are more complicated than the normal, smooth pattern, especially the one for the more complex duct. These residuals are, however, only illustrative of the more complex residual patterns that can occur for more complicated duct structures. Accurate residuals could only be determined by taking the average travel time at each frequency of a large number of rays started from a range of initial latitudes and with a range of initial wave normals. This is because when there is more than one closely spaced duct present, ray paths are strongly influenced by these initial parameters. It would be interesting to investigate whether the whistlers propagating in what is thought to be the same ducting structure have the same residual pattern or not, and whether there is any time-variation of this pattern. Further ray-tracing computations could be performed to investigate whether any predictions can be made as to the differences between the residual patterns of whistlers arising from second degree trapping at high altitude in different fine structure peaks of the same main duct structure as described by STRANGEWAYS(1982b). (3) Extra peaks might occur on the residual curves for whistlers leaving and then being trapped again in ducts or double-ducts (STRANGEWAYS,1982b), for super-whistlers (BERNHARDT, 1979), or for mixedpath modes (SMITH and CARPENTER, 1982) or following more complex paths as in the combined ducted plus sub-protonospheric propagation (STRANGEWAYS et al., 1983). (4) Propagation in the Earth-ionosphere waveguide will also increase the travel times of the rays and will sometimes have an effect on the residual pattern.
Normally, however, this has no real influence on the curve fitting or residuals, as demonstrated by simulations using the dispersion curves of atmospherics measured by TIXIER (1976). This is because the subionospheric travel time depends only weakly on the frequency and is mostly absorbed into the T parameter of the curve fitting. In some conditions, however, Earth-ionosphere waveguide mode propagation effects might appear in whistlers in the vicinity of the first or even the second order mode cut-off frequencies (SHIMAKURAet al., 1987), and may result in large positive residuals at these frequencies. Such residual peaks around 2 and 4 kHz were in fact noted on measured curves, but have not been studied in detail. (5) Various wave-wave or wave-particle interactions might also modify the whistler curves. From this point of view special note is to be taken of VLF emissions appearing in conjunction with whistlers, which indicate the presence of these processes. The foregoing discussion demonstrates that the measured residual curves may reflect the combined effects of various sources. Therefore, their interpretation requires a combination of curve fitting and raytracing studies together with other considerations. The analysis of repeatedly occurring whistler groups, consisting of whistler components which propagated along stable paths, may be particularly useful in studying the fine details of ducted propagation. To avoid averaging in obtaining the residuals, sophisticated digital methods are required in scaling the f-t pairs accurately. The matched-filter approach (HAMAR and TARCSAI, 1982) proved to be very efficient for accurate scaling. It also enables a study to be made of the fine structure of whistler traces, and thus yields additional information on the propagation conditions (TARCSAI et al., 1987).
6. CONCLUSIONS
Accurate curve fitting to a large number of whistlers, computed by numerical integration and raytracing in various magnetospheric plasma distributions, and the analysis of measured whistlers has demonstrated the following four points. (1) Commonly made approximations to the refractive index and ray path lead to whistler travel time errors of the order of several per cent. It is these approximations which hinder the accurate solution of the whistler integral equation for the field aligned electron density distribution. (2) These errors, combined with the use of a fixed plasma distribution in the curve fitting procedure,
Plasmaspheric
whistler interpretation
result in systematic errors of less than +20% in the computed equatorial electron densities in the range 2.5 < L < 3.5. At higher L-values, neq tends to be overestimated, and at lower L-values, underestimated. The corresponding systematid errors in the L-value of propagation and in the tube electron content will be of the order of several per cent for L > 2. The use of the centred dipole approximation in whistler analysis instead of more realistic magnetic field models (SEELY, 1977) increases the overall systematic errors in L and N, somewhat and tends to decrease the systematic error in neq and n, at L-values higher than about 2.22.5. Field aligned temperature gradients, which are not accounted for in routine whistler analysis, may result in a strong overestimation (even by a factor of 2) of neq, although further studies are needed in this respect. The ionospheric and sub-ionospheric contribution to travel times can be taken into account satisfactorily at middle and high latitudes. However, at low latitudes (below about L = 1.5), the effect of the ionosphere on travel times represents a serious problem in whistler diagnostics of the magnetosphere. (3) The standard deviations of the L-value and plasma densities obtained, resulting from random
251
scaling errors of whistler traces, remain below 10% in most cases when at least 20 f-t pairs are used in curve fitting. (4) The commonly used approximations give rise to systematic residuals in the whistler travel times, with a good agreement between measurements and theory. Departures from the smooth, normal residual curves yield a new research opportunity to study the propagation conditions of whistlers and other VLF signals. The further exploitation of travel time residual curves requires the use of digital techniques in whistler recording and scaling. The application of residual analysis in conjunction with digital matched filtering seems to be particularly useful, because the latter, in addition to yielding precise travel times, also gives information on whistler fine structure and thus on the detailed propagation conditions. Ackno&dgemmts-Part of this work was carried out at the Physics Department of Southampton University and the authors are grateful for the facilities extended to them there. Dr A. J. SMITH, Dr I. D. SMITHand Dr K. BULLOUGH are thanked for the provision of scaled whistler data from Halley, Antarctica. GY. TARCSAI acknowledges financial support from the Anglo-Hungarian Cultural Exchange Programme and from the National Fund for Scientific Research (Hungary) through grant No. 1292/86.
REFERENCES ANGERAMI J. J. and THOMAS J. 0. BARRINGTON R. E. and NISHIZAK~ T. BEGHIN C.
1964 1960 1966
BERNARD L. C. BERNHARDT P. A. BERNHARDTP. A. and PARK C. G. CARPENTER D. L. and CHAPPELL C. R. CARPENTER D. L., ANDERSON R. R., BELL T. F. and MILLER T. R. CORCUFF P. CORCUFF P., CORCUFF Y., CARPENTER D. L., CHAPPELLC. R., VIGNERONJ. and KLEIMENOVAN. DANIELLG. J
1973 1979 1977 1973 1981
J. geophys. Res. 69,4537. Can. J. Phys. 38, 1642. Electron Density Projiiles in Ionosphere and Exosphere, p. 587. North-Holland, Amsterdam. J. atmos. terr. Phys. 35, 871. J. geophys. Res. 84,5 13 1. J. geophys. Res. 82, 5222. J. geophys. Res. 18, 3062. Geophys. Res. Lett. 8, 1107.
1977 1972
Ann. Geophys. 33,443. Ann. Geophys. 28,679.
1986 1980
J. atmos. terr. Phys. 48, 27 1. J. atmos. terr. Phys. 42, 51.
1965 1977 1982 1977 1957 1968 1972 1982 1974 1975 1972 1982 1962 1987 1975
Radio Sci. 69D, 14 17. Planet. Space Sci. 25, 1087. Ann. Geophys. 38, 119. Planet. Space Sci. 25,871. J. atmos. terr. Phys. 11, 36. Geomagn. Aeron. 8, 899. Nature 237, 508. J. atmos. terr. Phys. 44,869. J. geophys. Res. 79,3739. Ann. Geophys. 31,69. Planet. Space Sci. 20, 73 1. J. Plasma Phys. 27, 199. Phys. Fluicis 5, 6. J. atmos. terr. Phys. 49, 1081. J. atmos. terr. Phys. 31, 1179.
DENBY M., ALEXANDERP. D., BULLOUGHK. and RYCROFTM. J. GUTHART H GUPTA G. P. and SINGH R. N. HARMARD. and TARCSAIGY. HASHIMOTOK., KIMURA I. and KUMAGI H. HINESC. 0. LIKHTERYA. I. and MOLCHANOV0. A. MATHUR A., RYCROFTM. J. and SAGREDOJ. L. MOREIRAA. OLXIN W. P. and PFITZERK. A. ROTH M SAGREDOJ. L. and BULLQUGHK. SAZHIN S. S. and SAZHINAE. M. SCARFF. L. SHIMAKURAS., TSUBAKIA. and HAYAKAWA M. SMITHA. J., SMITHI. D. and BULLOUGHK.
GY. TARCSAIet al.
258
1982 1962 1956 1982a 1982b 1986 1980 1983
J. atmos. terr. Phys. 44,973. Theory of Plasma Waves. McGraw Hill, New York. Can. J. Phys. 34,1153. J. atmos. terr. Phys. 44, 889. J. atmos. terr. Phys. 44,901. J. atmos. terr. Phys. 48, 67 1. J. atmos. terr. Phys. 42,983. J. atmos. terr. Phys. 45, 387.
1977a 1977b 1975 1985 1979 1988
J. atmos. terr. Phys. 39,415. J. atmos. terr. Phys. 39,427. J. atmos. terr. Phys. 37, 1447. Adv. Space Res. 5, 155. J. atmos. terr. Phys. 41, 967. J. atmos. terr. Phys. 50, 607.
B~~TTNER GY. B~~TTNER GY.
1971 1974 1979
PARKC. G.
1972
SEELYN. T.
1977
STRANGEWAYS H. J. STRANGEWAYS H. J.
1978 1984
TARCSA~ GY.
1977
TARCSAIGY.
1981
TARCSAIGY., HAMARD., SMITHA. J. and YEARBYK. H.
1987
TIXIERM.
1976
Ph.D. thesis, Southampton University. Private communication. Ionoszfera es magnetoszfera fizika VI., pp. 51-58. MTESZ KASZ, Budapest. Stanford University, Technical Report No. 34454-1, Stanford Electronics Labs. Stanford, California. Stanford University, Technical Report No. 3472-1, Stanford Electronics Labs. Stanford California. Ph.D. thesis, Southampton University. Ray-tracing in a warm magnetospheric plasma. Paper presented in session Hl at the 21st URSI Assembly, Florence, Italy. Fortran program for complete whistler analysis (a brief description for users). Circulated to members ofjoint URSI/IAGA WG on passive electromagnetic probing of the magnetosphere. Candidate of Science thesis. Hungarian Academy of Sciences, Budapest. Hyperfine structure of whistlers recorded digitally at Halley, Antarctica. Paper presented in session H2 at the 22nd URSI Assembly, Tel Aviv, Israel. These de Doctorat d’Etat, Universite de Poitiers.
SMITHA. J. and CARPENTER D. L. STIXT. H. STOREYL. R. 0. STRANGEWAYS H. J. STRANGEWAYS H. J. STRANGEWAYS H. J. STRANGEWAYS H. J. and RYCROFTM. J. STRANGEWAYS H. J., MADDENM. A. and RYCROFTM. J. STUARTG. F. STUARTG. F. TARCSAIGY. TARCSAIGY. TARCSAIGY. and DANIELLG. J. TARCSAIGY., SZEMEREDY P. and HEGYMEGI L. Reference is also made to the following unpublished material: ALEXANDER P.
D.
0021-9169/X9 $3.00+ .I0 Pergamon Press plc
Error sources and travel time residuals in plasmaspheric whistler interpretation GY. TARCSAI Geophysical Department, Eotvos University, Budapest, Kun B. ter 2, Hungary, H-1083 H. J. STRANGEWAYS Department of Electrical and Electronic Engineering, The University of Leeds, Leeds LS2 9JT, U.K. and
M. J. RYCROFT British Antarctic Survey, NERC, Madingley Road, Cambridge CB3 OET, U.K. (Received in final form I November 1988)
Abstract-In the interpretation of observed whistlers by curve fitting, systematic travel time residuals appeared which were studied by extensive simulations using ray-tracing, numerical integration and curve fitting. The residuals were found to originate from the commonly used approximations in the refractive index and ray path of whistler mode waves, which result in travel time increments or decrements, not accounted for in whistler interpretation. These approximations and the assumed form of the electron density distribution also lead to systematic errors in the diagnostics of plasmaspheric electron density by whistlers. In addition, the effects of other error sources, including random measurement errors, are also reviewed briefly. It is shown that the fine structure of residual curves is connected to propagation conditions. Thus, their study may yield a new research tool for studying whistler trapping, ducting structures and other features of whistler propagation. The application of residual analysis in conjunction with digital matched filtering of whistlers seems to be especially promising for further whistler studies.
1. INTRODUCTION
by numerical integration, ray-tracing and curvefitting. Emphasis is laid on the effect of these approximations on whistler travel times and on whistler and magnetospheric parameters obtained by curve fitting. In addition, the effect of the assumed form of electron density distribution on curve fitting results is also studied.
It has been previously shown (TARCSAI, 1975) that in the routine interpretation of whistlers by curve fitting, systematic travel time residuals appear. It has been suggested that these residuals can be attributed to various approximations common to routinely used whistler analysis techniques. The residuals were successfully reproduced later by curve fitting to theoretical whistlers computed by ray tracing (STRANGEWAYS, 1978) and by numerical integration (B-R, 1979), and found repeatedly on measured data (TARCSAI, 1981; STRANGEWAYSet al., 1983 ; TARCSAI
et al., 1987). It has been found that the observed systematic residuals appear as a result of three commonly used approximations. These are the omission of a ‘+ 1’ term in, and the contribution of ions to, the plasma refractive index, and also the assumption of purely longitudinal propagation rather than the real ray path, which executes snake-like excursions back and forth across the ducting structure. The present work investigates the importance of these approximations 249
2.
COMPUTATIONOF WHISTLER TRAVELTIMES
2.1. Ray-tracing In order to study the effect of approximations, theoretical frequency-time profiles for whistlers were computed by numerical integration and ray-tracing. The ray-tracing calculations were performed using a computer program developed in its original form by P. D. Alexander (ALEXANDER,1971) ; which uses the full refractive index p of a cold collisionless plasma containing both electrons and ions as given by STIX (1962) (in a number of runs the refractive index was determined for electrons only in order, by comparison, to determine the additional effect of ions on the
250
Gu. TARCSAIel al.
travel times). The travel times were computed for 1.516 kHz rays (in steps of 0.5 or 1.O kHz) assuming field aligned ducts of enhanced electron density centred on nine different L-values in the range L = l&4.5. The highest frequency computed for a given whistler was 16 kHz or, if lower, about half the minimum equatorial electron gyrofrequency along the given field line. The rays were started in the ducts with initial wave normal angles of O”, IO”, +20” and 30” to the geomagnetic field. A starting altitude of 1000 km was employed because the curve fitting procedure of TARCSAI (1975, 1977) was developed for a 1000 km1000 km altitude propagation path. The calculations were performed for six diffusive equilibrium models of plasmaspheric plasma distribution. These are the WN and SD models (ALEXANDER, 1971; DENBY et a/., 1980) and the DE-l-DE-4 models given by PARK (1972). In these models electrons and hydrogen, helium and oxygen ions are distributed according to diffusive equilibrium theory (ANGERAMIand THOMAS, 1964). Although in the majority of computations isothermal field lines were assumed, models including the electron temperature variation along geomagnetic field lines were also used in a modified ray-tracing model (STRANGEWAYS,1984, 1986). The ducts were modelled by Gaussian density enhancements superimposed on the background electron density (STRANGEWAYS and RYCROFT, 1980). The maximum fractional density enhancement 6 (at the centre of the duct) was taken to be 8, 10, 15, 20, 60 and 100%. The effective width of the ducts in the equatorial plane (where the enhancement is reduced to a tenth of its value at the duct centre) was chosen to be 430 km, although other ducts widths were sometimes used. To model the Earth’s magnetic field the same centred dipole approximation (PARK, 1972) was used in ray-tracing as in the numerical integration and curve fitting.
index p reduces to :
mode waves, the phase refractive p2=l+c
ffi 1
f(fHj
+
(3)
c,f) ’
where the summation is for the electrons and all the ionic species. f, are the plasma frequencies, fH,the gyrofrequencies and ci the sign of the charge of the electrons and all the ionic species which constitute the medium. In all of the methods used for whistler analysis the term ‘1’ and the contribution of ions on the right side of (3) are neglected. The plasmaspheric travel times were taken as :
s
1
a-) = 7 2cf
fbf"
(4)
(fH-fyds,
path
where f.and f,, are the plasma frequency and gyrofrequency for electrons, respectively, and the integral is taken for the path above 1000 km altitude. A series of f(t) pairs were computed for each simulated whistler by computing (1) numerically with p’ obtained from the full form of (3) or by neglecting the ‘l’, the ionic effects or both. In order to facilitate comparisons with the ray-tracing results, the plasma densities along the field lines of propagation were adjusted to be the same as those at the centre of the ducts. 2.3. Effect of approximations
on travel times
Comparison of the t(f) values obtained by raytracing and numerical integration for various cases
-
1
effectof
ims
2.2. Numerical integration Assuming strictly field aligned propagation, the whistler travel time, t, at frequency f can be computed by:
t(f) = i
c
s
/id&s,
(1)
path
a For exactly field aligned propagation
-30 g
I
w -40
where c is the speed of light in vacuum, & is the path element along the geomagnetic field line and p’ is the group refractive index given by :
P’ = q w-f).
---. ‘6
-50
j
/
I
,.A
effect of
and for whistler-
my path
“f_2._._._
/
2
I
4
I
6
I
6 Freqwncy
(2)
ducted
I
IO
I
12
I
14
I
16
(kHz)
Fig. 1. The effect of commonly made approximations on whistler travel times at L = 1.6, 2.4 and 3.5 as a function of frequency for the DE-I model.
251
Plasmaspheric whistler interpretation
2
I 2
I 4
I 6
I 6 Frequency
I IO
I 12
I 14
I 16
(kHz)
Fig. 2. The combined effect of the various approximations shown in Fig. 1 on whistler travel times.
permitted the separation of the effects of different approximations. These results are shown in Fig. 1 for the DE-l model and for three whistlers propagating at L = 1.6,2.4 and 3.5, respectively. The effects of the terms omitted from the refractive index were obtained by comparison of travel times (mainly obtained by numerical integration) with and without the respective terms included. Similar calculations by ray-tracing gave very similar results to numerical integration. The effect of the curved and snake-like ducted path was obtained by comparison of travel times computed by numerical integration with the wave normal always parallel to the geomagnetic field direction, and by raytracing, using the full refractive index formula (HINES, 1957) and Gaussian cross-section ducts. It is seen that the effect of ions is to increase travel times, especially at lower frequencies, as was shown by STOREY(1956) and BARRINGTONand NISHIZAKI (1960). The term ’ 1’ in equation (3) also leads to travel time increases of the order of several ms which are proportional to the wave frequency. The total increases in travel times caused by these two factors are 2.7-3.9% at L = 1.6 ; 1.6-2.2% at L = 2.4 and 0.7-1.4% at L = 3.5. On the other hand, the calculations show that along the true snake-like ray paths, the waves travel faster than they do along the axis of ducts in the purely longitudinal mode. This results in travel time decreases as can be seen in the lower part of Fig. 1, which effectively compensate for the effects of the approximations made in the refractive index. Figure 2 shows the combined effect on travel time of the approximations mentioned above, i.e. the difference between the exact travel times computed by ray-tracing using the full refractive index formula and those obtained by numerical integration for field aligned propagation with the approximate refractive index formulae [equation (4)]. It follows from Fig. 2 that the approximate
theoretical travel times used in routine whistler analysis may be lower (for lower L-values) or higher (for L-values greater than about 3) than the true values because of the commonly used approximations. The neglected effects in travel times amount to 2.1-3.6% ; 0.7-1.4% and -(O&1.3)%, at L = 1.6, 2.4 and 3.5, respectively, depending on the wave frequency. These figures are rather similar for different plasma distribution models and ducts and are practically independent of the initial wave normal angle. The consequence of these approximations upon whistler analysis is two-fold. Firstly, part of the neglected travel time increments is absorbed into the procedure of analysis and results in errors in the computed whistler and magnetospheric parameters. Secondly, the non-absorbed contributions appear as systematic travel time residuals in curve fitting which cannot be eliminated by any changes in the fitted whistler parameters. These points are treated in the following sections.
3. CURVE FITTING TO CALCULATEDWHISTLERS 3.1. The FIT procedure The computed whistlers, given by 8-16 f-t pairs, were analysed by a curve fitting method and computer code called FIT (TARCSAI, 1975, 1977). This method proved to be a reliable tool for the analysis of nose or non-nose whistlers. It is based on BERNARD’S (1973) formula for the dispersion o(f) = l.f’/2, which approximates the travel time integral (4) by an analytic expression. The method uses a least squares estimation of the zero frequency dispersion D,, the equatorial electron gyrofrequency fHE in Bernard’s approximation, and a third parameter, T, locating the time of the causative ‘spheric’. The FIT procedure minimises the sum of the squares of the differences between the observed and calculated propagation times. After obtaining these three parameters, the Lvalue of propagation, the electron density at 1000 km (n,) and at the equator (n,,) and the tube electron content (NT) were also calculated using the quasiconstants given by PARK (1972) and extended to lower latitudes by B~~TTNER(1974). The quasi-constants were determined by model calculations using equation (4) for the travel time, and assuming the DE-l model of plasma distribution (for routine whistler analysis this model is used within the plasmasphere). Together with whistler and plasma parameters, the FIT program also yields the estimated standard deviations and the final travel time residuals defined as the travel time of a whistler component at any frequency minus
252
GY. TARC~AI et
al.
the computed travel time at the same frequency calculated with the estimated whistler parameters.
Systematic errors in computed pommeters
:_?J--
3.2. Systematic errors in curve-fitting In order to check the internal accuracy of the FIT procedure, longitudinal whistlers were first analysed which had been computed from equation (4) and with model DE-l (TARCSAI,1981). The calculations demonstrated that the systematic errors in the estimated D0 and fHE (and also in nose frequency f, and the corresponding travel time t,,, which are derived from Do and fHE) are less than 2-3% in the range L = 1.6 5.0 (for L > 3 the bias in fn remains below 0.1%). The time shift between the true time origin and the estimated time of the originating ‘sferic’, T, was less than 2-3% of t, or for L < 2.5 less than 4 ms. The error in L, n,, neq and N, turned out to be less than 1, 3, 5 and l%, respectively, and for L 3 2.5 they remained below 1%. (Note that this is after constant corrections of 4.5%, 2% and 4% in n,, neq and N, respectively, were introduced to account for the approximations in Bernard’s formula and in the quasiconstants.) In conclusion, the internal accuracy of the FIT method is satisfactory and seems to be superior to that of other whistler analysis methods (SMITH et al., 1975 ; CORCUFF, 1977 ; STUART, 1977a,b). It should be noted, however, that it might be worthwhile replacing Bernard’s formula with the more accurate dispersion formula derived by DANIELL (1986) to help to obtain a value of neq that was not biased by an assumed form of electron density distribution along the field lines. The effects of the approximations discussed in Part 2 have been studied by analysing a great number of ray-traced whistlers by the FIT method. The results are summarised in Fig. 3 in the form of error bands, which represent the possible overall systematic errors corresponding to the commonly used approximations in the refractive index and ray path and resulting also from the assumed (fixed) form of electron density distribution and the small internal errors of the FIT procedure. Because the six plasma distributions used in the ray-tracing represent a rather broad range of plasmaspheric conditions representative of the variations under real circumstances, the obtained error bars are fairly realistic. Furthermore, the results obtained are characteristic (although perhaps an underestimate) of errors in other methods of whistler analysis, which are influenced by the same (and also, as a rule, by certain other) approximations as used in FIT. As may be seen in Fig. 3, for L 3 2 the FIT procedure gives the L-value of propagation of whistlers and the tube electron content, N, with errors of
I I.6
I
I
I
2
2.4
2.6
I 32
I
I
I
36
4
4.4
L- value
Fig. 3. Bands of possible systematic errors in the results of whistler interpretation, appearing as a consequence of common approximations in the refractive index, ray path and plasma density distribution (top, L-value; next equatorial electron density ; next total electron content along the geomagnetic field line; next electron density at 1000 km altitude ; bottom, time of causative ‘sferic’). several per cent. The magnitude of Arzeqvalues is less than f20% for 2.5 < L < 3.5. At higher L-values neq tends to be overestimated (eventually by as much as 30%) while below L = 2.2 it is strongly underestimated. The electron density at 1000 km will generally be strongly overestimated at L > 2.3, sometimes by as much as 120%, while in the range L = 1.6-2.3, An, varies in the range + 60%. Lastly, between L = 3 and 4.5, the time shift of the causative ‘sferic’ T may have a systematic error of the order of + 20 ms depending on the actual plasma distribution, with a tendency for T to be underestimated above about L = 3.5. Below L = 3.2, ATwill always be positive indicating that the FIT method places the causative ‘sferic’ earlier than it actually occurs. 4. DISCUSSION The behaviour of AT is the direct consequence of the neglected travel time increments, caused by the
Plasmaspheric whistler interpretation approximations in the ray path and refractive index. These approximations are responsible for the general trends in the error bands shown in Fig. 3, while the width of the error curves at a given L-value is caused by the difference between the fixed electron density distribution (DE-l model) and those distributions which may actually occur [as also indicated by ROTH (1975)]. The relative importance of these error sources depends on the L-value and on which parameter is being considered. The overestimation of the L-value and the strong underestimation of equatorial plane electron densities at low L-values is caused partly by neglect of the sinuous ray path but mainly by approximations in the refractive index. At higher latitudes, however, the unknown shape of the field aligned plasma density distribution is the main source of uncertainties in neq, NT and n,. It is worth noting here that the inversion of the integral equation (4) for electron densities along the field line (TARCSAIand DANIELL,1979) is completely free of any assumptions about the shape of the electron density distribution and is remarkably stable for noisy, measured f-t data. Its successful application is hindered, however, by the approximations in the refractive index and ray path discussed above which are very difficult to avoid. With respect to the errors shown in Fig. 3, it is important to note that at low L-values the errors in L and in neq are correlated, i.e. an overestimation of L is accompanied by an underestimation of neq. This means that the points on the neq vs L curve tend to remain on the profile, i.e. the difference between the true n,,(L’) value and the estimated n&(L’) value both of which are considered at the biased estimated L’ is smaller than the neq bias shown in Fig. 3 at the true (unknown) L or even at L’. This favourable circumstance only holds below about L = 2.3. At Lvalues above 2.3, AL and An,, are rather weakly (and, for various plasma distribution, oppositely) correlated. It is to be noted that the simulated whistlers used in this study covered the full range of frequencies between 1.5 kHz and f,,J2 or 16 kHz (if fHe> 32 kHz). The absence of higher frequency components in the whistler traces at low and mid-latitudes decreases the reliability of nose-extension methods and magnifies the systematic errors of the estimated parameters. For example, at L = 2 the reduction of the highest frequency from 16 kHz to 12.5 or 10 kHz increases AL from 0.08 to 0.12 and 0.16, respectively, with a parallel decrease in neq. This increase in AL counteracts the relation between AL and An,, discussed above. However, at higher L-values, e.g. at L = 3.5, the omission of the whistler components
253
above the nose frequency does not influence the results of curve fitting appreciably. In discussing systematic errors in curve-fitting, there are three additional factors which should be considered. The first is the approximate character of the centred dipole field used in routine whistler analysis, the effects of which have been studied by a number of researchers (LIKHTER and MOLCHANOV,1968; MATHIJRet al., 1972 ; SAGREDOand BULLOUGH,1972). These studies concentrated mainly on ring current effects which are important under disturbed conditions and outside the plasmasphere. For the present work, which refers to the study of whistlers inside the plasmasphere, it is pertinent to refer to the results of SEELY (1977). He constructed a hybrid quiet-time magnetospheric magnetic field model that combined the IGRF spherical harmonic series representation of the main field with the contributions from ring, tail and magnetopause current systems (OLSON and PFITZER, 1974). He showed that using the dipole model in whistler analysis results in L and NT values being overestimated by 24% and less than 15%, respectively, below L = 4.5 ; neqand n, values are correspondingly underestimated by 5-17% and l-15%, respectively. These errors, which are to be added to those discussed above, increase somewhat the systematic error in neq and n, at L-values higher than about 2.2-2.5. The plasma temperature should also be considered as it may influence the results of whistler analysis in a two-fold manner. Firstly, considering the existence of high thermal plasma temperatures in the magnetosphere, it should be determined whether it is appropriate to use the cold plasma approximation for the refractive index in whistler analysis. Although thermal effects may become significant close to the cyclotron resonance (SCARF, 1962; GUTHART, 1965) or beyond the plasmapause (MOREIRA,1982), for the usual whistler frequencies and inside the plasmasphere, the use of warm plasma theory (GUTHART, 1965 ; SAZHIN and SAZHINA,1982) in computing whistler ray paths or travel times produces insignificant changes compared to the cold plasma approximation (GUTHART, 1965 ; HASHIMOTO et aI., 1977 ; STRANGEWAYS, 1986) and, therefore, does not have any noticeable effect on whistler analysis (SEXLY,1977). Secondly, field aligned temperature gradients modify significantly the plasma distribution, a fact neglected in routine whistler analysis where isothermal distributions are used. SEELY(1977) has shown that, in general, a temperature increasing monotonically up the field line, and having its maximum gradient near 1000 km altitude, does not significantly change whistler values of equatorial parameters (L, n,,) with
GY. TARCSAI et al.
254
respect to values obtained using a constant temperature model. However, 1000 km altitude electron densities (n,) may be increased by a factor of 2 or more depending on the equatorial temperature. On the other hand, L and neq values inferred from whistlers are significantly affected by temperature gradients near the equatorial region. STRANGEWAYS (1984, 1986) conducted ray-tracing studies in diffusive equilibrium plasma density models including realistic temperature variations along geomagnetic field lines. He has shown that the equatorial electron density at L - 4.5 can be overestimated by as much as a factor of 2 if isothermal field lines are assumed in curve fitting when, in actual fact, a temperature gradient of the order of 1.6 K km- ’ at 1000 km (and 0.8 K km- ’ at 2000 km) exists. However, contrary to these results, CORCUFF et al. (1972), CARPENTERand CHAPPEL.L (1973), CARPENTERet al. (1981) and TARCSAI(1981) reported good agreement between in situ satellite observations of neq and determinations from curve fitting to whistlers assuming diffusive equilibrium along isothermal field lines. Therefore this question should be studied further, including the rigorous consideration of asymmetric conditions in the conjugate hemispheres (cf. STRANGEWAYS, 1982a). A third factor, leading to systematic errors in whistler analysis and not treated here in detail, is the travel time below 1000 km altitude, which is best subtracted from measured time delays before a curve-fitting analysis. This has two components. The first, the subionospheric propagation time, which may amount to several ms, is absorbed in the computed time-shift parameter T and may result in l-2% errors in the electron densities. The travel time in the ionosphere below 1000 km is generally computed by the formula given by PARK (1972) : At(f) = y, where At is given in s, the critical frequency of the F2layer, averaged for both hemispheres,f,F2, is in MHz and the wave frequency is in Hz. Ar may amount to several tens or hundreds of ms. Our experience and model computations show that an overestimate of f,F2 by 0.5 MHz results in about 0.5, -6 and -3% errors at mid-latitudes in L, neq and NT, respectively. These figures are lower at greater L-values and may be much higher at low L-values where ionospheric and magnetospheric travel times bcomes commensurate. In addition to systematic errors, the results of whistler analysis also contain unpredictable random errors. Figure 4 shows the variation of the median
\
I
errars In computed porarr&ef.s
Fkmdm
I
I
I
I
I
I
I
I
I5
2
2.5
3
35
4
4.5
5
L-value Fig. 4. Median values of the standard deviations of L, neq, n, and Nr as a function of the L-value of propagation. The curves were obtained from curve fitting to about 1500 measured whistlers and are normalised to 20 ,f-f pairs per trace.
value of normalised relative standard deviations of L, neq, n, and NT obtained from whistler curve fitting, as
a function of the L-value of propagation (TARCSAI, 1981). These curves were obtained from data consisting of more than 1200 mid-latitude whistlers (TARCSAI,197.5; TARCSAIet al., 1988) and 300 high latitude whistlers scaled from sonagrams and processed by the FIT method. The normalisation of the standard deviations was performed by reducing them to 20 points (f-t pairs) per whistler trace. Simulations on computed whistlers supported fully the characteristics of error propagation found from the measured data. Figure 4 demonstrates that, for most cases, the standard deviations of the L-value obtained and of the plasma densities remain below 10%. Thus, the systematic errors shown in Fig. 3 and discussed above are comparable to, or sometimes much greater than, the more unpredictable errors associated with whistler curve fitting. 5. SYSTEMATTC TRAVELTIME RESIDUALS 5.1. Measured and computed residuals In curve fitting applications, the residuals represent a very useful by-product which permits the checking of the validity of the theoretical model fitted to actual data. When the residuals are distributed in a nonrandom manner as a function of an independent variable, it indicates the need for an improvement of the model. That is the case with whistler curve fitting (TARCSAI,1975), which led to the simulation studies reported here. The travel time residuals are usually less than the random errors when conventional f-t
Plasmaspheric whistler interpretation
,-
L- 27SD 2 ducts --L-2fSD3hcta
2
4
6
8
10
12
14
255
to a ray-traced whistler (WN model, L = 2.4, 10% enhancement duct with a width parameter of 430 km). Here, no effort has been made to maximise the fit by variation of the ray-tracing plasmaspheric model or duct parameters. Further analysis of simulated whistlers has shown that the greatest part of the residuals is caused by the omission of ionic effects from the refractive index. Whistlers computed with ionic effects gave residual curves with peak-to-peak amplitudes from 4 to 8 times larger than those with the effect of ions neglected for both ray-traced and purely longitudinal simulated whistlers. This supports previous findings by BARRINGTONand NISHIZAKI (1960) and BEGHIN (1966) who, by using special techniques, demonstrated increased dispersion due to ions at the lowest frequencies of several low latitude whistlers. The same study was also attempted by GUPTA and SINGH (1977) with not entirely conclusive results. The term ‘1’ and the curved ray-path lead to residuals of only several tenths of a ms. It is to be noted that the narrower the frequency range of whistlers the lower will be the amplitude of the corresponding residual curves, while an extension of whistler analysis to frequencies as low as 0.7-1.0 kHz results in peakto-peak residual amplitudes of IO-12 ms.
I6
Frequency (kHz)
Fig. 5. Travel time residuals vs frequency for ducted whistlers, computed by ray-tracing for various plasma distributions and at different L-vatues. In panel (b) a measured cnrve constructed from 40 whistlers is afso shown.
scaling techniques are used. In this case averaging at given frequencies for a larger number of whistlers propagating under similar conditions leads to an effective suppression of noise, and systematic residuals appear as shown in Fig. 5b. The overall residual pattern is of the form fl with increasing frequency. Departures from the mean have a characteristic run of signs (+, -, + , -). The peak-to-peak amplitude is generally several ms, sometimes with multiple peaks on the smooth p-shaped curve. The curve fitting to longitudinally propagated whistlers computed by equation (4) has shown that, for such whistlers, the magnitude of residuals remains below about 0.1 ms, i.e. Bernard’s approximation is fairly good. However, the analysis of ray-traced whistlers, computed with the full refractive index, produced residual curves which agree very well with the observed ones. This is seen in Fig. 5b, which shows average residuals obtained from 40 measured whistlers (L = 1.9-2.5) and the residual curve belonging
5.2. Possible applications of residuals The smooth theoretical residual curves, shown in Figs. 5a, 5b and SC for various L-values and plasma distributions, may be considered as being ‘normal’ residuals produced by the known common approximations in whistler curve fitting. These curves vary somewhat with the parameters of the ray-tracing but for normal ducted whistlers always show the same smooth pattern. On the other hand, measured residual curves tend to exhibit a fine structure superposed on the normal pattern, which gives additional information on the propagation conditions. Extra peaks appear rather frequently, especially near the low and high frequency limits of the residual curves, as seen at about 3 and at 14.5 kHz on the measured curve of Fig. 5b (other minor peaks in this curve are due to random measurement errors). According to ray-tracing and curve fitting simulations, the anomalous residuals may perhaps be explained by the following factors. (1) Trapping and untrapping of whistler waves depends on the frequency. It has been shown (STRANGEWAYS, 1978) that the larger the altitude range at which different frequency rays enter and leave the duct, the more complex will be the fine structure in the residual pattern. Thus the altitude of the bases of the ducts and the altitude range of trapping and loss
256
Cu. TARCSAIet al.
might effect the residual pattern of the whistlers trapped in them. Since BERNHARDTand PARK (1977) suggest that the altitude of the end of the ducts will vary between day and night, it would be interesting to investigate whether there is any diurnal variation in the residual pattern of whistlers. (2) It is rather probable that whistlers do not propagate in simple, ‘ideal’ ducts, described by Gaussian or other simple cross-sections, but in field aligned irregularities which have a more complicated structure, forming a conglomeration of ducts, spaced more or less closely to each other. It turns out (STRANGEWAYS, 1978) that the residual patterns for whistlers trapped in ducts containing fine structure will be more complicated than for those trapped in ducts with simple structure. This is illustrated in Fig. 5d where the two residual patterns were obtained for the superposition of two and three Gaussian ducts, respectively. The elemental ducts were characterised by 6 = 15% and a width of 430 km, and were separated from each other by AL = 0.015 (- 100 km at the equator). Both curves are more complicated than the normal, smooth pattern, especially the one for the more complex duct. These residuals are, however, only illustrative of the more complex residual patterns that can occur for more complicated duct structures. Accurate residuals could only be determined by taking the average travel time at each frequency of a large number of rays started from a range of initial latitudes and with a range of initial wave normals. This is because when there is more than one closely spaced duct present, ray paths are strongly influenced by these initial parameters. It would be interesting to investigate whether the whistlers propagating in what is thought to be the same ducting structure have the same residual pattern or not, and whether there is any time-variation of this pattern. Further ray-tracing computations could be performed to investigate whether any predictions can be made as to the differences between the residual patterns of whistlers arising from second degree trapping at high altitude in different fine structure peaks of the same main duct structure as described by STRANGEWAYS(1982b). (3) Extra peaks might occur on the residual curves for whistlers leaving and then being trapped again in ducts or double-ducts (STRANGEWAYS,1982b), for super-whistlers (BERNHARDT, 1979), or for mixedpath modes (SMITH and CARPENTER, 1982) or following more complex paths as in the combined ducted plus sub-protonospheric propagation (STRANGEWAYS et al., 1983). (4) Propagation in the Earth-ionosphere waveguide will also increase the travel times of the rays and will sometimes have an effect on the residual pattern.
Normally, however, this has no real influence on the curve fitting or residuals, as demonstrated by simulations using the dispersion curves of atmospherics measured by TIXIER (1976). This is because the subionospheric travel time depends only weakly on the frequency and is mostly absorbed into the T parameter of the curve fitting. In some conditions, however, Earth-ionosphere waveguide mode propagation effects might appear in whistlers in the vicinity of the first or even the second order mode cut-off frequencies (SHIMAKURAet al., 1987), and may result in large positive residuals at these frequencies. Such residual peaks around 2 and 4 kHz were in fact noted on measured curves, but have not been studied in detail. (5) Various wave-wave or wave-particle interactions might also modify the whistler curves. From this point of view special note is to be taken of VLF emissions appearing in conjunction with whistlers, which indicate the presence of these processes. The foregoing discussion demonstrates that the measured residual curves may reflect the combined effects of various sources. Therefore, their interpretation requires a combination of curve fitting and raytracing studies together with other considerations. The analysis of repeatedly occurring whistler groups, consisting of whistler components which propagated along stable paths, may be particularly useful in studying the fine details of ducted propagation. To avoid averaging in obtaining the residuals, sophisticated digital methods are required in scaling the f-t pairs accurately. The matched-filter approach (HAMAR and TARCSAI, 1982) proved to be very efficient for accurate scaling. It also enables a study to be made of the fine structure of whistler traces, and thus yields additional information on the propagation conditions (TARCSAI et al., 1987).
6. CONCLUSIONS
Accurate curve fitting to a large number of whistlers, computed by numerical integration and raytracing in various magnetospheric plasma distributions, and the analysis of measured whistlers has demonstrated the following four points. (1) Commonly made approximations to the refractive index and ray path lead to whistler travel time errors of the order of several per cent. It is these approximations which hinder the accurate solution of the whistler integral equation for the field aligned electron density distribution. (2) These errors, combined with the use of a fixed plasma distribution in the curve fitting procedure,
Plasmaspheric
whistler interpretation
result in systematic errors of less than +20% in the computed equatorial electron densities in the range 2.5 < L < 3.5. At higher L-values, neq tends to be overestimated, and at lower L-values, underestimated. The corresponding systematid errors in the L-value of propagation and in the tube electron content will be of the order of several per cent for L > 2. The use of the centred dipole approximation in whistler analysis instead of more realistic magnetic field models (SEELY, 1977) increases the overall systematic errors in L and N, somewhat and tends to decrease the systematic error in neq and n, at L-values higher than about 2.22.5. Field aligned temperature gradients, which are not accounted for in routine whistler analysis, may result in a strong overestimation (even by a factor of 2) of neq, although further studies are needed in this respect. The ionospheric and sub-ionospheric contribution to travel times can be taken into account satisfactorily at middle and high latitudes. However, at low latitudes (below about L = 1.5), the effect of the ionosphere on travel times represents a serious problem in whistler diagnostics of the magnetosphere. (3) The standard deviations of the L-value and plasma densities obtained, resulting from random
251
scaling errors of whistler traces, remain below 10% in most cases when at least 20 f-t pairs are used in curve fitting. (4) The commonly used approximations give rise to systematic residuals in the whistler travel times, with a good agreement between measurements and theory. Departures from the smooth, normal residual curves yield a new research opportunity to study the propagation conditions of whistlers and other VLF signals. The further exploitation of travel time residual curves requires the use of digital techniques in whistler recording and scaling. The application of residual analysis in conjunction with digital matched filtering seems to be particularly useful, because the latter, in addition to yielding precise travel times, also gives information on whistler fine structure and thus on the detailed propagation conditions. Ackno&dgemmts-Part of this work was carried out at the Physics Department of Southampton University and the authors are grateful for the facilities extended to them there. Dr A. J. SMITH, Dr I. D. SMITHand Dr K. BULLOUGH are thanked for the provision of scaled whistler data from Halley, Antarctica. GY. TARCSAI acknowledges financial support from the Anglo-Hungarian Cultural Exchange Programme and from the National Fund for Scientific Research (Hungary) through grant No. 1292/86.
REFERENCES ANGERAMI J. J. and THOMAS J. 0. BARRINGTON R. E. and NISHIZAK~ T. BEGHIN C.
1964 1960 1966
BERNARD L. C. BERNHARDT P. A. BERNHARDTP. A. and PARK C. G. CARPENTER D. L. and CHAPPELL C. R. CARPENTER D. L., ANDERSON R. R., BELL T. F. and MILLER T. R. CORCUFF P. CORCUFF P., CORCUFF Y., CARPENTER D. L., CHAPPELLC. R., VIGNERONJ. and KLEIMENOVAN. DANIELLG. J
1973 1979 1977 1973 1981
J. geophys. Res. 69,4537. Can. J. Phys. 38, 1642. Electron Density Projiiles in Ionosphere and Exosphere, p. 587. North-Holland, Amsterdam. J. atmos. terr. Phys. 35, 871. J. geophys. Res. 84,5 13 1. J. geophys. Res. 82, 5222. J. geophys. Res. 18, 3062. Geophys. Res. Lett. 8, 1107.
1977 1972
Ann. Geophys. 33,443. Ann. Geophys. 28,679.
1986 1980
J. atmos. terr. Phys. 48, 27 1. J. atmos. terr. Phys. 42, 51.
1965 1977 1982 1977 1957 1968 1972 1982 1974 1975 1972 1982 1962 1987 1975
Radio Sci. 69D, 14 17. Planet. Space Sci. 25, 1087. Ann. Geophys. 38, 119. Planet. Space Sci. 25,871. J. atmos. terr. Phys. 11, 36. Geomagn. Aeron. 8, 899. Nature 237, 508. J. atmos. terr. Phys. 44,869. J. geophys. Res. 79,3739. Ann. Geophys. 31,69. Planet. Space Sci. 20, 73 1. J. Plasma Phys. 27, 199. Phys. Fluicis 5, 6. J. atmos. terr. Phys. 49, 1081. J. atmos. terr. Phys. 31, 1179.
DENBY M., ALEXANDERP. D., BULLOUGHK. and RYCROFTM. J. GUTHART H GUPTA G. P. and SINGH R. N. HARMARD. and TARCSAIGY. HASHIMOTOK., KIMURA I. and KUMAGI H. HINESC. 0. LIKHTERYA. I. and MOLCHANOV0. A. MATHUR A., RYCROFTM. J. and SAGREDOJ. L. MOREIRAA. OLXIN W. P. and PFITZERK. A. ROTH M SAGREDOJ. L. and BULLQUGHK. SAZHIN S. S. and SAZHINAE. M. SCARFF. L. SHIMAKURAS., TSUBAKIA. and HAYAKAWA M. SMITHA. J., SMITHI. D. and BULLOUGHK.
GY. TARCSAIet al.
258
1982 1962 1956 1982a 1982b 1986 1980 1983
J. atmos. terr. Phys. 44,973. Theory of Plasma Waves. McGraw Hill, New York. Can. J. Phys. 34,1153. J. atmos. terr. Phys. 44, 889. J. atmos. terr. Phys. 44,901. J. atmos. terr. Phys. 48, 67 1. J. atmos. terr. Phys. 42,983. J. atmos. terr. Phys. 45, 387.
1977a 1977b 1975 1985 1979 1988
J. atmos. terr. Phys. 39,415. J. atmos. terr. Phys. 39,427. J. atmos. terr. Phys. 37, 1447. Adv. Space Res. 5, 155. J. atmos. terr. Phys. 41, 967. J. atmos. terr. Phys. 50, 607.
B~~TTNER GY. B~~TTNER GY.
1971 1974 1979
PARKC. G.
1972
SEELYN. T.
1977
STRANGEWAYS H. J. STRANGEWAYS H. J.
1978 1984
TARCSA~ GY.
1977
TARCSAIGY.
1981
TARCSAIGY., HAMARD., SMITHA. J. and YEARBYK. H.
1987
TIXIERM.
1976
Ph.D. thesis, Southampton University. Private communication. Ionoszfera es magnetoszfera fizika VI., pp. 51-58. MTESZ KASZ, Budapest. Stanford University, Technical Report No. 34454-1, Stanford Electronics Labs. Stanford, California. Stanford University, Technical Report No. 3472-1, Stanford Electronics Labs. Stanford California. Ph.D. thesis, Southampton University. Ray-tracing in a warm magnetospheric plasma. Paper presented in session Hl at the 21st URSI Assembly, Florence, Italy. Fortran program for complete whistler analysis (a brief description for users). Circulated to members ofjoint URSI/IAGA WG on passive electromagnetic probing of the magnetosphere. Candidate of Science thesis. Hungarian Academy of Sciences, Budapest. Hyperfine structure of whistlers recorded digitally at Halley, Antarctica. Paper presented in session H2 at the 22nd URSI Assembly, Tel Aviv, Israel. These de Doctorat d’Etat, Universite de Poitiers.
SMITHA. J. and CARPENTER D. L. STIXT. H. STOREYL. R. 0. STRANGEWAYS H. J. STRANGEWAYS H. J. STRANGEWAYS H. J. STRANGEWAYS H. J. and RYCROFTM. J. STRANGEWAYS H. J., MADDENM. A. and RYCROFTM. J. STUARTG. F. STUARTG. F. TARCSAIGY. TARCSAIGY. TARCSAIGY. and DANIELLG. J. TARCSAIGY., SZEMEREDY P. and HEGYMEGI L. Reference is also made to the following unpublished material: ALEXANDER P.
D.