- Email: [email protected]
Mode theory of whistler ducts: integrated group delay times
Journal of Atmospheric and Terrestrial Physics, Vol. 54, No. 11/12. pp. 1599-1607, 1992. Printed in Great Britain.
0021-9169/92 $5.00+ .00 ~t'2 1992 Pergamon Press Ltd
Mode theory of whistler ducts: integrated group delay times M. J. LAIRD Department of Mathematics, King's College, Strand, London, WC2R 2LS, U.K. (Received in final form 30 March 1992; accepted 27 April 1992)
Abstract--Field-aligned enhancements of plasma density act as waveguides for whistler-mode waves. Calculations are presented of mode characteristics for a planar structure, where the enhancement has a sech-squared profile. Particular attention is paid to the group velocity for different waveguide modes. Application is then made to the magnetosphere by calculating the group delay per unit length for a given mode at different latitudes along a duct (taken to be centred on a dipole field line) and then integrating to obtain the total group delay over the magnetospheric part of the whistler path. Compared with strictly longitudinal propagation, group travel times are less at low frequencies,but greater at frequenciesnear to one-half the equatorial electron gyrofrequency, the magnitude of the residuals at the extremes of the frequency range being greater for higher-order modes. Differences between modes can be of the order of milliseconds.It therefore seems possible that recently reported fine structure seen in whistler spectrograms could be due in part to multi-mode propagation.
I. INTRODUCTION Whistler-mode waves in the magnetosphere may be channelled along magnetic field lines by means of ducts consisting of field-aligned enhancements of plasma density. The refractive index varies with the density and in consequence the ducts act in a similar way to dielectric waveguides. For routine whistler analysis it is usually assumed that propagation in the duct is purely longitudinal. The consequences of this, and of other commonly used approximations, have been studied in detail by TAI~CSAIet aL (1989), who found decreases in travel time when the actual snakelike ray path was taken into consideration. A matched filtering technique applied to whistlers recorded at Halley, Antarctica (HAMARet al., !990) and to whistlers observed on the Active (Intercosmos 24) spacecraft (LICHTENBERGERet al., 1991) has revealed fine structure. Whistlers which appear as a single trace on a normal spectrogram often contain a number of components, which can be individually analysed. Travel time residuals can differ by several milliseconds, and one possible reason for this is that different components may be related to different waveguide modes in the ducted magnetospheric part of the path (LICHTENBERGERet al., 1991). Open waveguides, such as ducts, differ in a number of ways from metal waveguides. For a simple planar structure, there is no cut-off for the lowest order mode; if the enhancement goes to zero, the mode fields become those of a plane wave. For higher order modes, at cut-off the fields extend to infinity, and the group velocity is then that for the ambient medium
(though this is not necessarily true for a cylindrical duct). Thus for a multi-moded duct one would expect significant differences in group velocity for different modes. That this is the case was found by SCARABUCCI and SMlTH (1971) in a study of the mode theory of bell-shaped (Gaussian profile) ducts. For a magnetospheric duct, both the physical and the geometric characteristics vary considerably along its length. For this reason, most studies of whistler propagation for realistic models have employed ray tracing. However WAL~ER (1966a) did propose that an alternative approach would be to evaluate the group refractive index for a given mode at a number of points along the duct, and then calculate the integrated group delay time by numerical quadrature. In this paper we shall carry out this procedure for a simple duct model, using a diffusive equilibrium model given in DENBY et al. (1980) for the ambient plasmaspheric electron density, and a centred dipole model of the geomagnetic field.
2. DUCT MODEL, MODE CONDITION AND GROUP DELAY We consider a plane-stratified duct for which the (constant) magnetic field B is parallel to the z-axis, and the electron number density N is a function of x alone. In their study, SCARABUCCIand SMITH (1971) employed both full wave methods and phase integral techniques. They found that the latter gave good results when checked against the former, and this is the method that we shall use. We shall only consider
1599
M. J. LAIRD
1600
propagation directions parallel to the xz-plane, so that for a wave vector k making an angle 0 with 0z (that is with B), ky is zero and kz is constant. For simplicity we shall consider a symmetrical enhancement centred on x = 0, and then, for a mode, one has the mode condition (WALKER,1966a)
io
k x d x = 2k=
f:o
Y the ratio of the electron gyrofrequency fB ( = c°a/2~) to f As kz is constant, and equal to k/sec 0, this may be written in the form
YsecO-sec2 0 = (~ou/ck~) 2 = F(x).
F(x) as defined here is proportional to the electron n u m b e r density, and so, for a fractional enhancement
tan0dx = (p-l)rc,
(1)
x0
F(0) = F(oo)(1 +6).
where p is a positive integer, and Xo, - x o are the reflection levels, where 0 is zero. To calculate the group delay, for a ray as depicted in Fig. 1, the group velocity is &o/c~k, where k and the angular frequency co ( = 2nf) are related by the dispersion equation. For a mode, the group delay per unit length zg--the reciprocal of the mode group velocity v~--is just the group time taken to go from A, where x = 0, to B, where x = x0, divided by the distance travelled along the duct, AC, and is therefore given by
¢~o zg
=
Jo
/ Oco
dx - - + J o
/ t3kx
1""0 Oo~tOkz dx--
&o/t3kx "
(2)
It can readily be shown that this is equivalent to the formulation used by WALKER (1966a) and SCARABUCCIand SMITH (1971).
3. D I S P E R S I O N
RELATION
(4)
AND WAVE
TRAPPING
We shall adopt the usual expression for the whistlermode refractive index in a cold plasma, so that the dispersion relation is (refractive index) 2 = (ck/o~) 2 = X / ( Y c o s O - 1),
(5)
Provided that F(x) < Y2/4, (4) gives two real values for sec 0, s(x) and S(x), with the properties
s ( x ) + S ( x ) = r,
s(x)S(x) = F(x),
S(x) > ~ Y > s(x).
(6)
The situation for wave trapping in a density enhancement, which requires Y > 2, is then illustrated in Fig. 2, which is a plot of sec 0 against F(x). Note that at the reflection levels _ x 0 , one value of sec 0 must be one, so that F(+_Xo) = Y - 1. For wave trapping we have the relationships ½Y > s(0) >~secO>~s(+_Xo) = 1 > s(oo).
(7)
The wave-normal angle 0 is then constrained to lie between + 0o, where sec 0o = s(0), giving the kind of ray path illustrated in Fig. 1. (Values of sec 0 less than one correspond to 0 imaginary and hence, for waveguide modes, to exponentially decaying fields outside the reflection levels.) The above picture is incomplete in that S(x) corresponds to a wave that can propagate at all values of x and so can leak energy from the duct (ADAcHI, 1965). However, from the work of SCARABUCCIand SMITH (1971) and LAIRD and NUNN (1975) this leakage is expected to be small for mag-
(3) X being the square of the ratio of the electron plasma frequency fN ( = ~ou/2n) to the wave frequency f, and
Y S(ool Y-I SIO)
B Y/2
sCO)
-x o
I F(oo~
I Y-I
I I F(O~' Y=/4
Ftx)
Fig. 1. A ducted ray path showing reflection at x = _+x0. Group delay per unit length is the group time from A to B divided by the distance A C.
Fig. 2. Sec 0 against F(x), where 0 is the wave-normal angle and F(x) is proportional to the electron number density. For a trapped wave, sec 0 is bounded by 1 and s(0).
Mode theory of whistler ducts netospheric ducts. An analytic treatment has been given by KARPMAN and KAUFMAN (1981a,b, 1982a,b, 1984). At lower cut-off, s ( ~ ) = 1 and so F ( ~ ) = Y - 1 . However, near Y = 2 another possibility may arise, namely that F(0) > Y2/4. If this happens, in terms of a ray picture, a ray starting from the reflection level x0 can no longer reach the centre of the duct. Instead, it is turned back where S = s = ~ Y, when the ray direction is parallel to B, and can then escape to infinity as sec0 follows the branch S(x). [Although the wave normal continues to turn towards higher refractive index, the ray direction turns in the opposite sense (SM1vn et al., 1960).] As noted by SCARABUCCl and SMITH (1971), this implies that there can be an upper cut-off. There is therefore a critical value Y,, when the conditions for lower cut-off and upper cut-off are satisfied simultaneously, which satisfies, using (5),
F(O)
=
Y~/4
=
F(~)(l+6)
= ( Y - l ) ( 1 +6), (8)
which on solving for Y gives
Y,. = 2(1 + 6 + (5(1 + 6)) ,/2).
(9)
To find the group velocity, we return to the dispersion relation (3) and write it in the form
k 2-
2
Ykk:+coN/c
~
= 0.
(10)
Taking #/Ok, and remembering that Y is inversely proportional to u~. we get
Ykk_ ~o 2 k - Y(f:k:+kS)+ ~o Ok
-
0.
f
,ru =
~0 j ,
(s 2 -
1)'"2F'(x)/J,
(s2Si?'2g'b)
'
(16) where F'(x) is expressed in terms of s via (4). For a linear or an exponential profile the integrals can be evaluated in terms of elementary functions. However, we may want a differentiable profile with a maximum, and that is the topic of the next section.
4. THE SECH-SQUARED ENHANCEMENT
A profile with the desired properties is obtained if we take the density enhancement to be proportional to sech2(x/a). This profile has been considered by KARPMAN and KAUVMAN(1982b), and is well known in other contexts (see for example BUDDEN, 1985). The parameter a is a measure of the width of the duct ; the half-width (the distance between the two points at which the enhancement is one-half of its maximum value) is 1.76~. We recall that F(x) is proportional to the electron number density N, and so (~oN/ck:) 2 = F(x) = F ( ~ ) ( I + 6 sechZ(x/~r))
= F ( O ) - S F ( ~ ) tanh2(x/a),
(17)
where 1 + 6 = F(O)/F(oo) = N(O)/N(oo). Then, from (4) and (6),
(12)
(Yk:/o~) &~/~k: = Y + cos 0( Ycos 0 - 2) (13)
where s = sec 0. Equation (4) can then be solved for sec 0 in terms of x : sec0 = s(x) = ~(Y--(y2--4F(x))L'2).
( Y - 2 s ) ( s 2 - 1)~,.2ds/F'(x) = ( p - ~)~,
(is)
(Yk:/e~) &o/?k~ = sin 0( Ycos 0 - 2)
= ( y s 2 - 2 s + Y)/s 2,
~(0)
- 2 k : jL
(11)
Dividing by k and taking components we obtain
= ( Y - 2 s ) ( s 2 - 1 ) t : 2 / s 2,
1601
6F(oo) sech2(x/cr) = ( S ~ , - s ) ( s - s ~ ),
(18)
6 F ( ~ ) tanh2(x/a) = ( S o - s ) ( s o - s ) ,
(19)
where S~, s~, So, So denote S(oo), s(vo), S(0), s(0), respectively. Differentiation of F(x) and substitution in (15) and (16) gives, for the mode condition and group delay per unit length zu :
(aO~(~Z)6'/2/C)lm = (p-- ~)~,
(14)
z.
Substitution in (1) and (2) gives the mode condition and group delay per unit length in integral form. These integrals can be evaluated numerically for an arbitrary density profile, but for an analytical treatment it turns out to be more convenient to take s as the variable of integration. When we do so, the mode condition and group delay per unit length are given by
= ( Yk:/~o)l,/l:,
(20) (21)
where
I,, =
f] s° ( Y - 2 s ) ( s z I, =
1) ds/(S,~ - s ) ( s - s . ~ ) Q ~'2, (22)
s2 d s / ( S ~ - s ) ( s - s ~ ) Q
1:2,
(23)
1602
Iz=
M. J. LAIRD
I "o( Y s 2 - 2 s +
,/i
Y ) d s / ( S o - s ) ( s - s ~ ) Q I/2. (24)
The three integrals above all include the square root of the quartic Q defined by
Q = (So -s)(so - s ) ( s - 1)(s÷ 1),
(25)
and so can be expressed in terms of elliptic integrals. These can easily be calculated with great accuracy and compared with tabulated values. Full details of the reduction of the integrals to standard form are given in the Appendix. The mode condition (20) becomes 4aogN(~)61/2
c(So - 1) ~12(So+ 1),/2 × {21-I(n., m) - rI(n~, m) - rI(n3, m)} = ( p - ~)z, (26) where l-I(n, m) denotes the complete elliptic integral of the third kind, and the parameters are given by
n, -= (so-1)/(so+l),
n2 = n , ( S ~ + l ) l ( S ~ - l ) ,
n3 = - n , ( 1 +s~)/(1 - s ~ ) ,
m = n~(So+ 1)/(S0-1).
(27) For given values of 6 and Y, (5) and (6) show that so, So and So are determined by So, and so the mode condition (26) can be regarded for integer p as an equation for So, and so one obtains a discrete set of values for 0o, the wave-normal angle at the duct centre. The advantages of an analytic solution are that relevant quantities are quickly calculated and that limiting cases are readily investigated. A general analysis shows that other profiles leading to elliptic integrals include the square-law (parabolic) profile, and an asymmetric profile where, in (17), an additional term proportional to tanh (x/a) may be added. Further, from examination of the nature of the singularities of the integrands in (15) and (16), any other profile with a single maximum may be expected to lead to broadly similar results.
In the limit 6 ~ 0, n~, n2 and m all vanish and so the other two elliptic integrals in the mode condition (26) take the value 7t/2. So the total number of modes, P, provided that Y exceeds the critical value Yc given by (9), may be expressed in the form (cf. the results of KARPMAN and KAUFraAN (1984) from the approximate treatment of a cylindrical duct) P =
[(GOON(O0)/C)(2(~) '/z(Y--2)-,12A(6, Y) +½], (29)
where the square brackets denote the integral part, and A(0, Y) --- 1 (except at Y = 2, where the value is (4/•)1n(1÷21/2), approximately 1.12). Plots of A(6, Y) against I/Y for several values of 6 are shown in Fig. 3. As one might expect, P is directly proportional to the duct width, as characterized by or, and to the ambient plasma frequency (i.e. to N:/2). The dependence on the size of the enhancement and the normalised frequency A ( = I / Y =fife) is more complex, but for small enhancements P is roughly proportional to 6 u2 and at low frequencies to A ~/2. At upper cut-off, So = So = Y/2. Reference to (22) and (25) shows that I,, can be evaluated in terms of elementary functions, as detailed in the Appendix. The maximum number of modes, P . . . . occurs when Y = Y~ which for small enhancements is approximately 2(I+6~/2). Inserting this in (29), as a first approximation emax =
[1.12tro9~(~)61/4/c+½],
(30)
the error being of order 3 3/4 . We now turn to the group velocity. Denote by V the group velocity for propagation parallel to the magnetic field. Then from (3) and (13) V - 2 ( Y - 1)~o _ 2 c ( Y - 1 ) 3/2
kY
Xi/2y
2c(fn-f)312f 1/2 fNfs ' (31)
1.5
~o 30 1.2 do
I.I
5. M O D E
PROPERTIES:
CUT-OFFS
AND GROUP
VELOCITY
At lower cut-off, s~ = 1, so that the parameter n 3 goes to minus infinity and the associated integral, II(n3, m) vanishes. Then So~ = Y - 1 , so that So and So are the roots of
s z - Ys+ ( r - 1)(1 +6) = 0.
(28)
1.0 0.0
0.I
, 0.5
0.2
0.4
I 0.5
y-I
Fig. 3. The function A(6, Y) against Y- 1for different values of the percentage enhancement ft. The broken curve shows A(~, Y:) where Y~is the upper cut-offlimit for Y. (N.B. for Y> 2, A(0, Y) = 1.)
1603
Mode theory of whistler ducts which is the usual starting point for calculating plasmaspheric travel times. For waveguide modes, the group velocity vg is related to the propagation of energy along the guide, and so for low-order modes, where the fields are confined to the neighbourhood of the duct centre, one would expect it to be close to V0 [corresponding to fN(0)], whereas at lower cut-off, where the fields extend to infinity, it would be V~. Investigation of the limiting values of I,/1. for so = 1 and for s~ = 1 confirms this. At upper cut-off, however, the situation is rather different. Here one finds that the group velocity is equal to the phase velocity for the mode, og/k;. As at this cut-off F(O) = soSo = Y2/4, one finds that the group velocity is just ~'/B/2J~v(O), which is independent of frequency (a well-known result for propagation when Ycos 0 = 2). Hence, at upper cut-off, e~/l/~ = Y2/4( Y - 1) 3'2 = sec 2 00/(2 sec 0 0 - l) 3'2, (32) where 0o is the wave-normal angle at x = O.
Figure 4 is a plot of vJV0 [where Vq is the reciprocal of z~ given by (21)] against 00 for a 10% enhancement and different values of the normalised frequency A ( = I / Y ) . At low frequencies it is seen that vq increases with 00 to V~ ( = V0(1 +6)~/2) at cut-off. As the critical value A,. (0.349...) is approached, vg first drops significantly before rising to the cut-off value. Well above A,., t'u decreases from V0 to the upper cutoff value. Thus for a 10% enhancement the group velocity can be from about 5% greater than Vo to 20% less. Of course, the mode condition only allows for a discrete set of values of 00 for given model parameters. SCARABUCCIand SMITH (1971) noted that
AIlil'ude (Mm) 5I I . O ~ 15 I0I
0.5--
~ 1
fN(MHz) ~ 0.0
I0 20 Distance (Mrn)
Fig. 5. Duct width parameter a and normalised frequency A, as fractions of their equatorial values, against distance in megametres from the equator along the field line L = 4. The scale at the top is altitude. Also shown is the electron plasma frequency f~ (MHz) for the winter night model of RYCROVr and ALEXANDER(1969), described in DENBYet al. (1980).
at higher frequencies the group velocity decreased with increasing mode number p, but also stated that in the low frequency band (A < 0.25) the group velocities for all modes were the same. This is not necessarily the case for our duct profile, though around A = 0.1 the plot is quite flat for 0o < 16". It should be added that group velocities do depend on the profile used ; for example, for a simple slab model the group velocity decreases with increasing 0o until very close to lower cut-off. 6. VARIATION ALONG A DUCT
0.5 0"9C
FC
f 20°
\l
30°
\
i 40*
I
00 Fig. 4. Group velocity vg (normalized to the group velocity for parallel propagation at the duct centre, V0) against wavenormal angle 00 at the duct centre, for an enhancement of 10%, and for different values of the normalised frequency A ( = 1/Y). The broken curves give the group velocity at lower cut-off (A < A,.) and upper cut-off (A > Ac), where A, = 0.349 ....
In our model, the mode properties described in the previous section depend on four quantities: the enhancement factor 6, the duct width as characterised by ~r, the ambient electron plasma frequency, and A, the frequency normalised to the electron gyrofrequency. In applying our results to the magnetosphere we shall suppose that 6 is constant along a duct centred on a dipole field line with L-value L0, though as the integration along a duct is numerical there is no particular problem with ~i varying along the duct, or with different magnetic field models. Then along the duct the other three quantities vary considerably, as shown in Fig. 5, in which ~/(Teq and A/Aeq ( = fs,,q/fB) are plotted as functions of distance from the equator
1604
M. J. LAIRD
measured along the field line with L equal to four. Also shown is the plasma frequency for the winter night model of RYCROFT and ALEXANDER (1969) and described in DENBY et al. (1980). The plots terminate at an altitude of 1000 km, below which the electron density rises rapidly. It can be seen that tr changes by a factor greater than ten. The changes in fN and f e compensate each other to a certain extent : starting at 1000 km altitude, the refractive index for propagation parallel to the magnetic field typically decreases somewhat to a minimum between 1400 and 2200 km altitude (STRANGEWAYS and RYCROFT, 1980), and then increases to the equator to about four times its minimum value. For example, for a frequency of 3 kHz the Appleton-Lassen formula gives values for the refractive index of 8.9 at 1000km altitude, 8.3 at 1600 km, and 31.2 at the equator for the field line with L equal to four. (A correction for ions needs to be made at low altitudes, but the trend is clear.) A major effect of these variations, both geometrical and physical, is on the number of modes that can be supported. For a 10% enhancement with •eq = 100 km, only a few modes exist at low altitudes (P = 4 at 4000 km for a 2 kHz wave) whereas many may theoretically exist at the equator (77 at 2kHz), again for L equal to four. Another important length scale when one seeks to integrate over the length of a duct is the distance along the duct covered by a complete oscillation of the ray path, that is, four times the distance A C in Fig. 1. A C is just the denominator in expression (16) for za, and hence for our model 4 A C = 2a(6F(oo))'/2I~.
(33)
By taking the limit of lz as So ~ 1 (which corresponds to an oscillation with infinitesimal displacement from the duct centre) one can obtain a minimum value z0 given by
(34)
(As Zo depends only on the value of N(x) near x = 0, the second expression applies to any profile with a maximum there, as can be checked by expanding the dispersion relation to find 09 to first order in x 2 and k~.) Along a dipole field line, at magnetic latitude 2, a = aeqCOS32/(l+3sin22)I/2,
(35)
and the element of length, dl, is given by dl= Lrecos2(l+3sin22)l/2d2,
f
dl/zo = (Lre/2naeq)(6/2(1 + 6))1/2(4 tan 2 - 32).
(37) For L = 4 , 2 = 5 0 ° , aeq= 100 km and 6 = 0 . 1 , this gives a value of 18.58. This will decrease with higher frequency, but not by much--for most of its path a whistler is in the low frequency regime. At the equatorial cut-off Y~q = 2, the number of oscillations for the above values of the parameters is 15.82. The results from applying formula (34) for z0 seem to agree reasonably well with published ray-tracing studies, for example STRANGEWAYSand RYCROFT (1980) (though it should be noted that they use an enhancement proportional to exp (-x2/2tr~), so that in making comparisons one should put a = 21/2ffd). For the waveguide modes, one can calculate, from (33), a similar quantity z p ( f ) ( = 4 A C ) for mode number p at frequency f, and integrate numerically along a field line. Generally, z p ( f ) increases (and so the number of oscillations decreases) with p; the behaviour with respect to f is more complex. For the same values as above, for p = 1 the number of oscillations is 17.2 at 2kHz decreasing to 15.7 at 6.5 kHz; for p = 3, the number is 13.4 at 2 kHz rising to 14.3 at 4.5 kHz and then falling to 13.7 at 6.5 kHz.
7. INTEGRATED MODE GROUP DELAY TIMES
The group delay time Tg(20, f ) along a duct centred on L = L0 between latitude 20 and the magnetic equator is taken to be given by Ta(20, f ) =
fl °za dl = Lor e dol
"ca
× cos 2(1 + 3 sin22) ~/2 d2,
z0 = 2rca(2(l + 6 ) ( Y - 1 ) / 6 ( Y - 2 ) ) '/2 = 4rc((Y- I)/(Y-2))~/Z(-N(0)/N"(0))'/2.
such a ray path between the equator and latitude 2 is
(36)
where re is the Earth's radius. Thus, in the low frequency limit (Y>> 1), the number of oscillations of
(38)
where za(2, f ) is the group delay per unit length at latitude 2 for frequency f For selected values of Lo, 6 and 6eq, and given f, one can, at any latitude, calculate P, the total number of modes that can propagate. The mode condition may then be solved for the modes required and the group delays per unit length calculated. After performing these operations at a suitable number of different latitudes along the duct, one can integrate (38) numerically to give the integrated group delay time for each desired mode at the given frequency. Repetition of this procedure for a number of frequencies permits plots such as those shown in Figs 6 and 7 to be drawn. For the purpose of illustration, a 10% enhance-
Mode theory of whistler ducts P 6
5 -r ~4
2
[ -i
I i
o
I 2
~rg I ms)
Fig. 6. Residual group delay time G (ms) against frequency f (kHz) for the lowest three modes for a 10% enhancement centred on L = 4. a,,q was 100km and the integration was carried out between the equator and magnetic latitude 50'. ment centred on L = 4 was chosen. According to STRANGEWAYS and RYCROFT (1980), waves may become trapped in ducts over a range of altitudes between 1000 and 7000km. We have taken 20 to be 50, corresponding to an altitude of about 4000 km. In both figures, the group delay for propagation strictly parallel to the magnetic field has been subtracted, allowing the residual group delay time t o ( f ) to be exhibited. Only modes that exist over the whole latitude and frequency range (2-6.5 kHz) are considered. Figure 6 shows results for the three lowest modes for a,,u = 100 kin, whereas Fig. 7 shows results for the six P
1605
lowest modes for a 200km duct. Common features are that, at low frequencies the residual delay is negative and at higher frequencies positive. These can be understood by reference to Fig. 4 and the accompanying discussion in Section 5. Around 3.5 kHz the mode dispersion is low. At the extremes of the frequency range shown, the absolute value of the residual increases with mode number p. Extension of the integration range to a lower altitude increases the magnitude of the residuals at low frequencies and decreases them at the top end of the scale, and raises the frequency for minimum mode dispersion. (TARCSAI et al. (1989), with a starting altitude of 1000 kin, obtained negative residuals at all frequencies, the greatest magnitude being at the lowest frequency.) Compared with the a,~q = 100km duct, modes 1, 2 and 3 for the 200km duct have smaller residuals: this corresponds to the fact that more modes can propagate in the wider duct, so that the lower order mode fields are confined more closely to the duct centre (as a proportion of the duct width). Calculations for narrower ducts lead to greater residuals. Increasing the enhancement increases the number of modes, which, as explained above, tends to reduce the size of the residuals for the lower order modes, but at the same time it increases the range of possible mode group velocities, which has the opposite effect: computations have to be carried out in individual cases. Finally, it should be noted that the delays shown in the figures relate only to one half of the supposed magnetospheric path ; symmetry between ).0 and - 2 , would double the delays.
8. D I S C U S S I O N
5
4
O
I
I
I
2
tg ( m s )
Fig. 7. Residual group delay time ta (ms) against frequency f (kHz) for the lowest six modes for a 10% enhancement centred on L = 4. ~q was 200 km and the integration was carried out between the equator and magnetic latitude 50'.
The main aim of this paper has been to see if mode theory can give fine structure on the scale of milliseconds; the results presented indicate that it can. Figures 6 and 7 represent what one might expect to see from an impulsive source if the modes shown were sufficiently excited. However, there are a number of factors that have not been taken into consideration, and a number of questions unanswered. Among the factors are the use of the full Appleton-Lassen formula for the refractive index rather than the highdensity limit (3), and the contribution of ions, both discussed by TARCSM et al, (1989) with regard to parallel propagation. There are also warm plasma effects (SAZHIN et al., 1990), and cross-duct gradients in the magnetic field strength and plasma density. All of these could he included, if necessary, by numerical evaluation of the integrals occurring in the mode condition and the expression for group delay. Curvature
1606
M.J. LAIRD
is another important factor--we have essentially been assuming that the duct is straight. The appropriate parameter is a/pB, where ps is the radius of curvature of a field line. Within the plasmasphere, PB does not change much along a field line, and so the most marked effects are near the magnetic equator, where the combined effects of the curvature and cross-duct gradients are to increase the value of 6 needed for g u i d a n c e (LAIRD, 1981; STRANGEWAYS, 1991), or alternatively, for given 3, to reduce the effective duct enhancement and width, and to shift the effective duct centre to higher L-value. These effects are frequencydependent. Although some work has been done (BOOKER, 1962; WALKER, 1966b, 1971), the theory including curvature, especially for wide ducts, needs to be extended. Unanswered questions include the excitation of modes (and also their amplification or damping via interactions with energetic particles). Usually one would expect the low-order modes to predominate; however, scattering by irregularities could excite higher-order modes, and the variation of duct characteristics along the duct, which is appreciable over a
ray oscillation, could be expected to lead to coupling between the modes. Another possible mechanism to explain the observations is propagation in complex duct structures (HAMAR et aL, 1990; LICHTENBERGERet al., 1991). One interesting case, under active consideration, is of a narrow duct or ducts within a broader duct. There are then two kinds of mode, one largely confined to the narrow duct and the other for which the narrow duct is a perturbation. A ray-tracing study by STRANGEWAYS(1982) showed that rays first trapped in the main duct at low altitude may be further trapped within fine structure enhancements at higher altitude. However, to answer these questions one probably needs some form of full-wave theory. What is clear is that whistler analysis may be expected to continue to yield information about plasma distribution within the plasmasphere. Acknowledgements--I am indebted to Dr A. J. Smith for
drawing my attention to this problem, and to him, Dr Gy. Tarcsai and Dr J. Lichtenberger for discussingtheir whistler observations and analysis.
REFERENCES
Radio Sci. J. Res. NBS 69D, 493. J. geophys. Res. 67, 4135.
BUDDENK. G.
1965 1962 1985
BYRDP. F. and FRIEDMANM. D.
1971
DENBYM., BULLOUGHK., ALEXANDERP. D.
1980
Handbook o f Elliptic Integrals for Engineers and Scientists, 2nd edn. Springer, Berlin. J. atmos, terr. Phys. 42, 51.
1990
J. atmos, terr. Phys. 52, 801.
1981a 1981b 1982a 1982b 1984 1981 1975 1991
JETP Lett. 33, 252. Soy. Phys. JETP 53, 956. J. plasma Phys. 27, 225. Soy. J. Plasma Phys. 8, 180 (translation of Fiz. Plazmy 8, 319). Planet. Space Sci. 32, 1505. J. atmos, terr. Phys. 43, 81. Planet. Space Sci. 23, 1649. J. geophys. Res. 96, 21,149.
1990 1971 1960 1982 1991 1980 1989 1966a 1966b 1971
Ann. Geophys. 8, 273. Radio Sci. 6, 65, J. geophys, Res. 65, 815. J. atmos, terr. Phys. 44, 901. J. atmos, terr. Phys. 53, 151. J. atmos, terr. Phys. 42, 983. J. atmos, terr. Phys. 51,249. J. atmos, terr. Phys. 28, 747. J. atmos, terr. Phys. 28, 807. Proc. R. Soc. Lond. A. 321, 69.
ADACHI S. BOOKER H. G.
The Propagation of Radio Waves. Cambridge Uni-
versity Press, Cambridge.
and RYCROFT M. J. HAMAR D., TARCSAIGY., LICHTENBERGERJ., SMITH A. J. and YEARBY K. H. KARPMAN V. I. and KAUFMANR. N. KARPMAN V. I. and KAUFMANR. N. KARPMAN V. I. and KAUFMANR. N. KARPMAN V. I. and KAUFMANR. N.
KARPMANV. I. and KAUFMANR. N. LAIRD M. J.
LAIRDM. J. and NUNND. LICHTENBERGER J., TARCSAIGY., P.~,SZTORS., FERENCZ C., HAMAR n . , MOLCHANOVO. A. and GOLYAVIN A. SAZHIN S., SMITH A. J. and SAZmNA E. M. SCARABUCCIR. R. and SMITH R. L. SMITHR. L., HELLIWELLR. A. and YABROFFI. W. STRANGEWAYSH. J. STRANGEWAYSn . J. STRANGEWAYSH. J. and RYCROFT M. J. TARCSAI GY., STRANGEWAYSH. J. and RYCROFT M. J.
WALKERA. D. M. WALKERA. D. M. WALKERA. D. M.
Reference is also made to the following unpublished material:
RYCROFTM. J. and ALEXANDERP.D.
1969
Paper A.3.11. presented at the 12th Plenary Meeting of COSPAR, Prague, Czechoslovakia.
Mode theory of whistler ducts APPENDIX
To evaluate the integrals I,., I, and/=, defined by (22), (23) and (24), put the rational part of the integrand into partial fractions. Hence
1~ = Y l o + 2 1 , - ( S ~ - l ) I 2 - ( l - s ~ ) I 3 , 1, = (S212+s[13)/(S~--s~)-lo, I= = ti1, + 13 - 12,
( S ~ - 1)12 = 9
~
,2 (1 - n t sin2 4)) d4)/ (1 - n 2 sin 2 4))(1 - m sin 2 4))1,2, (A8)
'n,2
(AI) (A2)
1607
(1 -s~)13 = 9
(A3)
t
(1 --n, sin24)) d4)/
0
(l-n3sin24))(1-msin24)) ''2.
lo= f[°ds/Q ''2,
,%
1,= f[sds/Q ~'2, II(n,m)=
I2= f f ds/(S~-s)Q 1/2, I3= f[°ds/(s-so~)Q':2.(A4) When Q is factored as in (25), the appropriate transformation is (BYRD and FRIEDMAN, 197l, p. 112 et seq.) sin 2 4) = (s-- 1)/n,(s + 1),
(A5)
d4)/(l-nsin24))(1-msin24)) ~,~2. (AI0)
For the complete ray path, the corresponding incomplete elliptic integrals are required. At upper cut-off, so = So = ~ Y, and Im is given by
~Im = fl Y/2 ( s 2 - l ) l ' e d s / ( S ~ - s ) ( s - s ~ ) ,
s~ < 1,
where n~ and other relevant parameters are defined by (27). One then finds that
n2
1o = 9
I do
d4)/(l - m sin 24)) 1,2 = 9K(m),
(A6)
where 9 = 2/(So-1)~Z2(So+ 1)~z2, and K(m) is the complete elliptic integral of the first kind. The other integrals are given by I~ = 9
(A9)
I~, 12 and 13 may then be expressed in terms of K(m) and complete elliptic integrals of the third kind
where
(l +n~ sin24)) d4)/ ( 1 - n l sin24))(1-rosin24)) I/2, (A7)
(All) which, as stated, is readily integrated in terms of elementary functions. At the critical value Y,, given by (9), s~ = 1, S~ = Y,.- 1, and the integral becomes ½1,.= fl Y''2 (s+ 1)"2(s- 1) L'2ds/(Y,.--1--s)
= y,!/2(y,.--2)-~/2 In { y,~/Z(y,,+2),.'2 + Yc+ l } - c o s h '~Y,., (= 26 ~/4In (1 + 2 ~/2)+0(6'/4)).
(A12)
0021-9169/92 $5.00+ .00 ~t'2 1992 Pergamon Press Ltd
Mode theory of whistler ducts: integrated group delay times M. J. LAIRD Department of Mathematics, King's College, Strand, London, WC2R 2LS, U.K. (Received in final form 30 March 1992; accepted 27 April 1992)
Abstract--Field-aligned enhancements of plasma density act as waveguides for whistler-mode waves. Calculations are presented of mode characteristics for a planar structure, where the enhancement has a sech-squared profile. Particular attention is paid to the group velocity for different waveguide modes. Application is then made to the magnetosphere by calculating the group delay per unit length for a given mode at different latitudes along a duct (taken to be centred on a dipole field line) and then integrating to obtain the total group delay over the magnetospheric part of the whistler path. Compared with strictly longitudinal propagation, group travel times are less at low frequencies,but greater at frequenciesnear to one-half the equatorial electron gyrofrequency, the magnitude of the residuals at the extremes of the frequency range being greater for higher-order modes. Differences between modes can be of the order of milliseconds.It therefore seems possible that recently reported fine structure seen in whistler spectrograms could be due in part to multi-mode propagation.
I. INTRODUCTION Whistler-mode waves in the magnetosphere may be channelled along magnetic field lines by means of ducts consisting of field-aligned enhancements of plasma density. The refractive index varies with the density and in consequence the ducts act in a similar way to dielectric waveguides. For routine whistler analysis it is usually assumed that propagation in the duct is purely longitudinal. The consequences of this, and of other commonly used approximations, have been studied in detail by TAI~CSAIet aL (1989), who found decreases in travel time when the actual snakelike ray path was taken into consideration. A matched filtering technique applied to whistlers recorded at Halley, Antarctica (HAMARet al., !990) and to whistlers observed on the Active (Intercosmos 24) spacecraft (LICHTENBERGERet al., 1991) has revealed fine structure. Whistlers which appear as a single trace on a normal spectrogram often contain a number of components, which can be individually analysed. Travel time residuals can differ by several milliseconds, and one possible reason for this is that different components may be related to different waveguide modes in the ducted magnetospheric part of the path (LICHTENBERGERet al., 1991). Open waveguides, such as ducts, differ in a number of ways from metal waveguides. For a simple planar structure, there is no cut-off for the lowest order mode; if the enhancement goes to zero, the mode fields become those of a plane wave. For higher order modes, at cut-off the fields extend to infinity, and the group velocity is then that for the ambient medium
(though this is not necessarily true for a cylindrical duct). Thus for a multi-moded duct one would expect significant differences in group velocity for different modes. That this is the case was found by SCARABUCCI and SMlTH (1971) in a study of the mode theory of bell-shaped (Gaussian profile) ducts. For a magnetospheric duct, both the physical and the geometric characteristics vary considerably along its length. For this reason, most studies of whistler propagation for realistic models have employed ray tracing. However WAL~ER (1966a) did propose that an alternative approach would be to evaluate the group refractive index for a given mode at a number of points along the duct, and then calculate the integrated group delay time by numerical quadrature. In this paper we shall carry out this procedure for a simple duct model, using a diffusive equilibrium model given in DENBY et al. (1980) for the ambient plasmaspheric electron density, and a centred dipole model of the geomagnetic field.
2. DUCT MODEL, MODE CONDITION AND GROUP DELAY We consider a plane-stratified duct for which the (constant) magnetic field B is parallel to the z-axis, and the electron number density N is a function of x alone. In their study, SCARABUCCIand SMITH (1971) employed both full wave methods and phase integral techniques. They found that the latter gave good results when checked against the former, and this is the method that we shall use. We shall only consider
1599
M. J. LAIRD
1600
propagation directions parallel to the xz-plane, so that for a wave vector k making an angle 0 with 0z (that is with B), ky is zero and kz is constant. For simplicity we shall consider a symmetrical enhancement centred on x = 0, and then, for a mode, one has the mode condition (WALKER,1966a)
io
k x d x = 2k=
f:o
Y the ratio of the electron gyrofrequency fB ( = c°a/2~) to f As kz is constant, and equal to k/sec 0, this may be written in the form
YsecO-sec2 0 = (~ou/ck~) 2 = F(x).
F(x) as defined here is proportional to the electron n u m b e r density, and so, for a fractional enhancement
tan0dx = (p-l)rc,
(1)
x0
F(0) = F(oo)(1 +6).
where p is a positive integer, and Xo, - x o are the reflection levels, where 0 is zero. To calculate the group delay, for a ray as depicted in Fig. 1, the group velocity is &o/c~k, where k and the angular frequency co ( = 2nf) are related by the dispersion equation. For a mode, the group delay per unit length zg--the reciprocal of the mode group velocity v~--is just the group time taken to go from A, where x = 0, to B, where x = x0, divided by the distance travelled along the duct, AC, and is therefore given by
¢~o zg
=
Jo
/ Oco
dx - - + J o
/ t3kx
1""0 Oo~tOkz dx--
&o/t3kx "
(2)
It can readily be shown that this is equivalent to the formulation used by WALKER (1966a) and SCARABUCCIand SMITH (1971).
3. D I S P E R S I O N
RELATION
(4)
AND WAVE
TRAPPING
We shall adopt the usual expression for the whistlermode refractive index in a cold plasma, so that the dispersion relation is (refractive index) 2 = (ck/o~) 2 = X / ( Y c o s O - 1),
(5)
Provided that F(x) < Y2/4, (4) gives two real values for sec 0, s(x) and S(x), with the properties
s ( x ) + S ( x ) = r,
s(x)S(x) = F(x),
S(x) > ~ Y > s(x).
(6)
The situation for wave trapping in a density enhancement, which requires Y > 2, is then illustrated in Fig. 2, which is a plot of sec 0 against F(x). Note that at the reflection levels _ x 0 , one value of sec 0 must be one, so that F(+_Xo) = Y - 1. For wave trapping we have the relationships ½Y > s(0) >~secO>~s(+_Xo) = 1 > s(oo).
(7)
The wave-normal angle 0 is then constrained to lie between + 0o, where sec 0o = s(0), giving the kind of ray path illustrated in Fig. 1. (Values of sec 0 less than one correspond to 0 imaginary and hence, for waveguide modes, to exponentially decaying fields outside the reflection levels.) The above picture is incomplete in that S(x) corresponds to a wave that can propagate at all values of x and so can leak energy from the duct (ADAcHI, 1965). However, from the work of SCARABUCCIand SMITH (1971) and LAIRD and NUNN (1975) this leakage is expected to be small for mag-
(3) X being the square of the ratio of the electron plasma frequency fN ( = ~ou/2n) to the wave frequency f, and
Y S(ool Y-I SIO)
B Y/2
sCO)
-x o
I F(oo~
I Y-I
I I F(O~' Y=/4
Ftx)
Fig. 1. A ducted ray path showing reflection at x = _+x0. Group delay per unit length is the group time from A to B divided by the distance A C.
Fig. 2. Sec 0 against F(x), where 0 is the wave-normal angle and F(x) is proportional to the electron number density. For a trapped wave, sec 0 is bounded by 1 and s(0).
Mode theory of whistler ducts netospheric ducts. An analytic treatment has been given by KARPMAN and KAUFMAN (1981a,b, 1982a,b, 1984). At lower cut-off, s ( ~ ) = 1 and so F ( ~ ) = Y - 1 . However, near Y = 2 another possibility may arise, namely that F(0) > Y2/4. If this happens, in terms of a ray picture, a ray starting from the reflection level x0 can no longer reach the centre of the duct. Instead, it is turned back where S = s = ~ Y, when the ray direction is parallel to B, and can then escape to infinity as sec0 follows the branch S(x). [Although the wave normal continues to turn towards higher refractive index, the ray direction turns in the opposite sense (SM1vn et al., 1960).] As noted by SCARABUCCl and SMITH (1971), this implies that there can be an upper cut-off. There is therefore a critical value Y,, when the conditions for lower cut-off and upper cut-off are satisfied simultaneously, which satisfies, using (5),
F(O)
=
Y~/4
=
F(~)(l+6)
= ( Y - l ) ( 1 +6), (8)
which on solving for Y gives
Y,. = 2(1 + 6 + (5(1 + 6)) ,/2).
(9)
To find the group velocity, we return to the dispersion relation (3) and write it in the form
k 2-
2
Ykk:+coN/c
~
= 0.
(10)
Taking #/Ok, and remembering that Y is inversely proportional to u~. we get
Ykk_ ~o 2 k - Y(f:k:+kS)+ ~o Ok
-
0.
f
,ru =
~0 j ,
(s 2 -
1)'"2F'(x)/J,
(s2Si?'2g'b)
'
(16) where F'(x) is expressed in terms of s via (4). For a linear or an exponential profile the integrals can be evaluated in terms of elementary functions. However, we may want a differentiable profile with a maximum, and that is the topic of the next section.
4. THE SECH-SQUARED ENHANCEMENT
A profile with the desired properties is obtained if we take the density enhancement to be proportional to sech2(x/a). This profile has been considered by KARPMAN and KAUVMAN(1982b), and is well known in other contexts (see for example BUDDEN, 1985). The parameter a is a measure of the width of the duct ; the half-width (the distance between the two points at which the enhancement is one-half of its maximum value) is 1.76~. We recall that F(x) is proportional to the electron number density N, and so (~oN/ck:) 2 = F(x) = F ( ~ ) ( I + 6 sechZ(x/~r))
= F ( O ) - S F ( ~ ) tanh2(x/a),
(17)
where 1 + 6 = F(O)/F(oo) = N(O)/N(oo). Then, from (4) and (6),
(12)
(Yk:/o~) &~/~k: = Y + cos 0( Ycos 0 - 2) (13)
where s = sec 0. Equation (4) can then be solved for sec 0 in terms of x : sec0 = s(x) = ~(Y--(y2--4F(x))L'2).
( Y - 2 s ) ( s 2 - 1)~,.2ds/F'(x) = ( p - ~)~,
(is)
(Yk:/e~) &o/?k~ = sin 0( Ycos 0 - 2)
= ( y s 2 - 2 s + Y)/s 2,
~(0)
- 2 k : jL
(11)
Dividing by k and taking components we obtain
= ( Y - 2 s ) ( s 2 - 1 ) t : 2 / s 2,
1601
6F(oo) sech2(x/cr) = ( S ~ , - s ) ( s - s ~ ),
(18)
6 F ( ~ ) tanh2(x/a) = ( S o - s ) ( s o - s ) ,
(19)
where S~, s~, So, So denote S(oo), s(vo), S(0), s(0), respectively. Differentiation of F(x) and substitution in (15) and (16) gives, for the mode condition and group delay per unit length zu :
(aO~(~Z)6'/2/C)lm = (p-- ~)~,
(14)
z.
Substitution in (1) and (2) gives the mode condition and group delay per unit length in integral form. These integrals can be evaluated numerically for an arbitrary density profile, but for an analytical treatment it turns out to be more convenient to take s as the variable of integration. When we do so, the mode condition and group delay per unit length are given by
= ( Yk:/~o)l,/l:,
(20) (21)
where
I,, =
f] s° ( Y - 2 s ) ( s z I, =
1) ds/(S,~ - s ) ( s - s . ~ ) Q ~'2, (22)
s2 d s / ( S ~ - s ) ( s - s ~ ) Q
1:2,
(23)
1602
Iz=
M. J. LAIRD
I "o( Y s 2 - 2 s +
,/i
Y ) d s / ( S o - s ) ( s - s ~ ) Q I/2. (24)
The three integrals above all include the square root of the quartic Q defined by
Q = (So -s)(so - s ) ( s - 1)(s÷ 1),
(25)
and so can be expressed in terms of elliptic integrals. These can easily be calculated with great accuracy and compared with tabulated values. Full details of the reduction of the integrals to standard form are given in the Appendix. The mode condition (20) becomes 4aogN(~)61/2
c(So - 1) ~12(So+ 1),/2 × {21-I(n., m) - rI(n~, m) - rI(n3, m)} = ( p - ~)z, (26) where l-I(n, m) denotes the complete elliptic integral of the third kind, and the parameters are given by
n, -= (so-1)/(so+l),
n2 = n , ( S ~ + l ) l ( S ~ - l ) ,
n3 = - n , ( 1 +s~)/(1 - s ~ ) ,
m = n~(So+ 1)/(S0-1).
(27) For given values of 6 and Y, (5) and (6) show that so, So and So are determined by So, and so the mode condition (26) can be regarded for integer p as an equation for So, and so one obtains a discrete set of values for 0o, the wave-normal angle at the duct centre. The advantages of an analytic solution are that relevant quantities are quickly calculated and that limiting cases are readily investigated. A general analysis shows that other profiles leading to elliptic integrals include the square-law (parabolic) profile, and an asymmetric profile where, in (17), an additional term proportional to tanh (x/a) may be added. Further, from examination of the nature of the singularities of the integrands in (15) and (16), any other profile with a single maximum may be expected to lead to broadly similar results.
In the limit 6 ~ 0, n~, n2 and m all vanish and so the other two elliptic integrals in the mode condition (26) take the value 7t/2. So the total number of modes, P, provided that Y exceeds the critical value Yc given by (9), may be expressed in the form (cf. the results of KARPMAN and KAUFraAN (1984) from the approximate treatment of a cylindrical duct) P =
[(GOON(O0)/C)(2(~) '/z(Y--2)-,12A(6, Y) +½], (29)
where the square brackets denote the integral part, and A(0, Y) --- 1 (except at Y = 2, where the value is (4/•)1n(1÷21/2), approximately 1.12). Plots of A(6, Y) against I/Y for several values of 6 are shown in Fig. 3. As one might expect, P is directly proportional to the duct width, as characterized by or, and to the ambient plasma frequency (i.e. to N:/2). The dependence on the size of the enhancement and the normalised frequency A ( = I / Y =fife) is more complex, but for small enhancements P is roughly proportional to 6 u2 and at low frequencies to A ~/2. At upper cut-off, So = So = Y/2. Reference to (22) and (25) shows that I,, can be evaluated in terms of elementary functions, as detailed in the Appendix. The maximum number of modes, P . . . . occurs when Y = Y~ which for small enhancements is approximately 2(I+6~/2). Inserting this in (29), as a first approximation emax =
[1.12tro9~(~)61/4/c+½],
(30)
the error being of order 3 3/4 . We now turn to the group velocity. Denote by V the group velocity for propagation parallel to the magnetic field. Then from (3) and (13) V - 2 ( Y - 1)~o _ 2 c ( Y - 1 ) 3/2
kY
Xi/2y
2c(fn-f)312f 1/2 fNfs ' (31)
1.5
~o 30 1.2 do
I.I
5. M O D E
PROPERTIES:
CUT-OFFS
AND GROUP
VELOCITY
At lower cut-off, s~ = 1, so that the parameter n 3 goes to minus infinity and the associated integral, II(n3, m) vanishes. Then So~ = Y - 1 , so that So and So are the roots of
s z - Ys+ ( r - 1)(1 +6) = 0.
(28)
1.0 0.0
0.I
, 0.5
0.2
0.4
I 0.5
y-I
Fig. 3. The function A(6, Y) against Y- 1for different values of the percentage enhancement ft. The broken curve shows A(~, Y:) where Y~is the upper cut-offlimit for Y. (N.B. for Y> 2, A(0, Y) = 1.)
1603
Mode theory of whistler ducts which is the usual starting point for calculating plasmaspheric travel times. For waveguide modes, the group velocity vg is related to the propagation of energy along the guide, and so for low-order modes, where the fields are confined to the neighbourhood of the duct centre, one would expect it to be close to V0 [corresponding to fN(0)], whereas at lower cut-off, where the fields extend to infinity, it would be V~. Investigation of the limiting values of I,/1. for so = 1 and for s~ = 1 confirms this. At upper cut-off, however, the situation is rather different. Here one finds that the group velocity is equal to the phase velocity for the mode, og/k;. As at this cut-off F(O) = soSo = Y2/4, one finds that the group velocity is just ~'/B/2J~v(O), which is independent of frequency (a well-known result for propagation when Ycos 0 = 2). Hence, at upper cut-off, e~/l/~ = Y2/4( Y - 1) 3'2 = sec 2 00/(2 sec 0 0 - l) 3'2, (32) where 0o is the wave-normal angle at x = O.
Figure 4 is a plot of vJV0 [where Vq is the reciprocal of z~ given by (21)] against 00 for a 10% enhancement and different values of the normalised frequency A ( = I / Y ) . At low frequencies it is seen that vq increases with 00 to V~ ( = V0(1 +6)~/2) at cut-off. As the critical value A,. (0.349...) is approached, vg first drops significantly before rising to the cut-off value. Well above A,., t'u decreases from V0 to the upper cutoff value. Thus for a 10% enhancement the group velocity can be from about 5% greater than Vo to 20% less. Of course, the mode condition only allows for a discrete set of values of 00 for given model parameters. SCARABUCCIand SMITH (1971) noted that
AIlil'ude (Mm) 5I I . O ~ 15 I0I
0.5--
~ 1
fN(MHz) ~ 0.0
I0 20 Distance (Mrn)
Fig. 5. Duct width parameter a and normalised frequency A, as fractions of their equatorial values, against distance in megametres from the equator along the field line L = 4. The scale at the top is altitude. Also shown is the electron plasma frequency f~ (MHz) for the winter night model of RYCROVr and ALEXANDER(1969), described in DENBYet al. (1980).
at higher frequencies the group velocity decreased with increasing mode number p, but also stated that in the low frequency band (A < 0.25) the group velocities for all modes were the same. This is not necessarily the case for our duct profile, though around A = 0.1 the plot is quite flat for 0o < 16". It should be added that group velocities do depend on the profile used ; for example, for a simple slab model the group velocity decreases with increasing 0o until very close to lower cut-off. 6. VARIATION ALONG A DUCT
0.5 0"9C
FC
f 20°
\l
30°
\
i 40*
I
00 Fig. 4. Group velocity vg (normalized to the group velocity for parallel propagation at the duct centre, V0) against wavenormal angle 00 at the duct centre, for an enhancement of 10%, and for different values of the normalised frequency A ( = 1/Y). The broken curves give the group velocity at lower cut-off (A < A,.) and upper cut-off (A > Ac), where A, = 0.349 ....
In our model, the mode properties described in the previous section depend on four quantities: the enhancement factor 6, the duct width as characterised by ~r, the ambient electron plasma frequency, and A, the frequency normalised to the electron gyrofrequency. In applying our results to the magnetosphere we shall suppose that 6 is constant along a duct centred on a dipole field line with L-value L0, though as the integration along a duct is numerical there is no particular problem with ~i varying along the duct, or with different magnetic field models. Then along the duct the other three quantities vary considerably, as shown in Fig. 5, in which ~/(Teq and A/Aeq ( = fs,,q/fB) are plotted as functions of distance from the equator
1604
M. J. LAIRD
measured along the field line with L equal to four. Also shown is the plasma frequency for the winter night model of RYCROFT and ALEXANDER (1969) and described in DENBY et al. (1980). The plots terminate at an altitude of 1000 km, below which the electron density rises rapidly. It can be seen that tr changes by a factor greater than ten. The changes in fN and f e compensate each other to a certain extent : starting at 1000 km altitude, the refractive index for propagation parallel to the magnetic field typically decreases somewhat to a minimum between 1400 and 2200 km altitude (STRANGEWAYS and RYCROFT, 1980), and then increases to the equator to about four times its minimum value. For example, for a frequency of 3 kHz the Appleton-Lassen formula gives values for the refractive index of 8.9 at 1000km altitude, 8.3 at 1600 km, and 31.2 at the equator for the field line with L equal to four. (A correction for ions needs to be made at low altitudes, but the trend is clear.) A major effect of these variations, both geometrical and physical, is on the number of modes that can be supported. For a 10% enhancement with •eq = 100 km, only a few modes exist at low altitudes (P = 4 at 4000 km for a 2 kHz wave) whereas many may theoretically exist at the equator (77 at 2kHz), again for L equal to four. Another important length scale when one seeks to integrate over the length of a duct is the distance along the duct covered by a complete oscillation of the ray path, that is, four times the distance A C in Fig. 1. A C is just the denominator in expression (16) for za, and hence for our model 4 A C = 2a(6F(oo))'/2I~.
(33)
By taking the limit of lz as So ~ 1 (which corresponds to an oscillation with infinitesimal displacement from the duct centre) one can obtain a minimum value z0 given by
(34)
(As Zo depends only on the value of N(x) near x = 0, the second expression applies to any profile with a maximum there, as can be checked by expanding the dispersion relation to find 09 to first order in x 2 and k~.) Along a dipole field line, at magnetic latitude 2, a = aeqCOS32/(l+3sin22)I/2,
(35)
and the element of length, dl, is given by dl= Lrecos2(l+3sin22)l/2d2,
f
dl/zo = (Lre/2naeq)(6/2(1 + 6))1/2(4 tan 2 - 32).
(37) For L = 4 , 2 = 5 0 ° , aeq= 100 km and 6 = 0 . 1 , this gives a value of 18.58. This will decrease with higher frequency, but not by much--for most of its path a whistler is in the low frequency regime. At the equatorial cut-off Y~q = 2, the number of oscillations for the above values of the parameters is 15.82. The results from applying formula (34) for z0 seem to agree reasonably well with published ray-tracing studies, for example STRANGEWAYSand RYCROFT (1980) (though it should be noted that they use an enhancement proportional to exp (-x2/2tr~), so that in making comparisons one should put a = 21/2ffd). For the waveguide modes, one can calculate, from (33), a similar quantity z p ( f ) ( = 4 A C ) for mode number p at frequency f, and integrate numerically along a field line. Generally, z p ( f ) increases (and so the number of oscillations decreases) with p; the behaviour with respect to f is more complex. For the same values as above, for p = 1 the number of oscillations is 17.2 at 2kHz decreasing to 15.7 at 6.5 kHz; for p = 3, the number is 13.4 at 2 kHz rising to 14.3 at 4.5 kHz and then falling to 13.7 at 6.5 kHz.
7. INTEGRATED MODE GROUP DELAY TIMES
The group delay time Tg(20, f ) along a duct centred on L = L0 between latitude 20 and the magnetic equator is taken to be given by Ta(20, f ) =
fl °za dl = Lor e dol
"ca
× cos 2(1 + 3 sin22) ~/2 d2,
z0 = 2rca(2(l + 6 ) ( Y - 1 ) / 6 ( Y - 2 ) ) '/2 = 4rc((Y- I)/(Y-2))~/Z(-N(0)/N"(0))'/2.
such a ray path between the equator and latitude 2 is
(36)
where re is the Earth's radius. Thus, in the low frequency limit (Y>> 1), the number of oscillations of
(38)
where za(2, f ) is the group delay per unit length at latitude 2 for frequency f For selected values of Lo, 6 and 6eq, and given f, one can, at any latitude, calculate P, the total number of modes that can propagate. The mode condition may then be solved for the modes required and the group delays per unit length calculated. After performing these operations at a suitable number of different latitudes along the duct, one can integrate (38) numerically to give the integrated group delay time for each desired mode at the given frequency. Repetition of this procedure for a number of frequencies permits plots such as those shown in Figs 6 and 7 to be drawn. For the purpose of illustration, a 10% enhance-
Mode theory of whistler ducts P 6
5 -r ~4
2
[ -i
I i
o
I 2
~rg I ms)
Fig. 6. Residual group delay time G (ms) against frequency f (kHz) for the lowest three modes for a 10% enhancement centred on L = 4. a,,q was 100km and the integration was carried out between the equator and magnetic latitude 50'. ment centred on L = 4 was chosen. According to STRANGEWAYS and RYCROFT (1980), waves may become trapped in ducts over a range of altitudes between 1000 and 7000km. We have taken 20 to be 50, corresponding to an altitude of about 4000 km. In both figures, the group delay for propagation strictly parallel to the magnetic field has been subtracted, allowing the residual group delay time t o ( f ) to be exhibited. Only modes that exist over the whole latitude and frequency range (2-6.5 kHz) are considered. Figure 6 shows results for the three lowest modes for a,,u = 100 kin, whereas Fig. 7 shows results for the six P
1605
lowest modes for a 200km duct. Common features are that, at low frequencies the residual delay is negative and at higher frequencies positive. These can be understood by reference to Fig. 4 and the accompanying discussion in Section 5. Around 3.5 kHz the mode dispersion is low. At the extremes of the frequency range shown, the absolute value of the residual increases with mode number p. Extension of the integration range to a lower altitude increases the magnitude of the residuals at low frequencies and decreases them at the top end of the scale, and raises the frequency for minimum mode dispersion. (TARCSAI et al. (1989), with a starting altitude of 1000 kin, obtained negative residuals at all frequencies, the greatest magnitude being at the lowest frequency.) Compared with the a,~q = 100km duct, modes 1, 2 and 3 for the 200km duct have smaller residuals: this corresponds to the fact that more modes can propagate in the wider duct, so that the lower order mode fields are confined more closely to the duct centre (as a proportion of the duct width). Calculations for narrower ducts lead to greater residuals. Increasing the enhancement increases the number of modes, which, as explained above, tends to reduce the size of the residuals for the lower order modes, but at the same time it increases the range of possible mode group velocities, which has the opposite effect: computations have to be carried out in individual cases. Finally, it should be noted that the delays shown in the figures relate only to one half of the supposed magnetospheric path ; symmetry between ).0 and - 2 , would double the delays.
8. D I S C U S S I O N
5
4
O
I
I
I
2
tg ( m s )
Fig. 7. Residual group delay time ta (ms) against frequency f (kHz) for the lowest six modes for a 10% enhancement centred on L = 4. ~q was 200 km and the integration was carried out between the equator and magnetic latitude 50'.
The main aim of this paper has been to see if mode theory can give fine structure on the scale of milliseconds; the results presented indicate that it can. Figures 6 and 7 represent what one might expect to see from an impulsive source if the modes shown were sufficiently excited. However, there are a number of factors that have not been taken into consideration, and a number of questions unanswered. Among the factors are the use of the full Appleton-Lassen formula for the refractive index rather than the highdensity limit (3), and the contribution of ions, both discussed by TARCSM et al, (1989) with regard to parallel propagation. There are also warm plasma effects (SAZHIN et al., 1990), and cross-duct gradients in the magnetic field strength and plasma density. All of these could he included, if necessary, by numerical evaluation of the integrals occurring in the mode condition and the expression for group delay. Curvature
1606
M.J. LAIRD
is another important factor--we have essentially been assuming that the duct is straight. The appropriate parameter is a/pB, where ps is the radius of curvature of a field line. Within the plasmasphere, PB does not change much along a field line, and so the most marked effects are near the magnetic equator, where the combined effects of the curvature and cross-duct gradients are to increase the value of 6 needed for g u i d a n c e (LAIRD, 1981; STRANGEWAYS, 1991), or alternatively, for given 3, to reduce the effective duct enhancement and width, and to shift the effective duct centre to higher L-value. These effects are frequencydependent. Although some work has been done (BOOKER, 1962; WALKER, 1966b, 1971), the theory including curvature, especially for wide ducts, needs to be extended. Unanswered questions include the excitation of modes (and also their amplification or damping via interactions with energetic particles). Usually one would expect the low-order modes to predominate; however, scattering by irregularities could excite higher-order modes, and the variation of duct characteristics along the duct, which is appreciable over a
ray oscillation, could be expected to lead to coupling between the modes. Another possible mechanism to explain the observations is propagation in complex duct structures (HAMAR et aL, 1990; LICHTENBERGERet al., 1991). One interesting case, under active consideration, is of a narrow duct or ducts within a broader duct. There are then two kinds of mode, one largely confined to the narrow duct and the other for which the narrow duct is a perturbation. A ray-tracing study by STRANGEWAYS(1982) showed that rays first trapped in the main duct at low altitude may be further trapped within fine structure enhancements at higher altitude. However, to answer these questions one probably needs some form of full-wave theory. What is clear is that whistler analysis may be expected to continue to yield information about plasma distribution within the plasmasphere. Acknowledgements--I am indebted to Dr A. J. Smith for
drawing my attention to this problem, and to him, Dr Gy. Tarcsai and Dr J. Lichtenberger for discussingtheir whistler observations and analysis.
REFERENCES
Radio Sci. J. Res. NBS 69D, 493. J. geophys. Res. 67, 4135.
BUDDENK. G.
1965 1962 1985
BYRDP. F. and FRIEDMANM. D.
1971
DENBYM., BULLOUGHK., ALEXANDERP. D.
1980
Handbook o f Elliptic Integrals for Engineers and Scientists, 2nd edn. Springer, Berlin. J. atmos, terr. Phys. 42, 51.
1990
J. atmos, terr. Phys. 52, 801.
1981a 1981b 1982a 1982b 1984 1981 1975 1991
JETP Lett. 33, 252. Soy. Phys. JETP 53, 956. J. plasma Phys. 27, 225. Soy. J. Plasma Phys. 8, 180 (translation of Fiz. Plazmy 8, 319). Planet. Space Sci. 32, 1505. J. atmos, terr. Phys. 43, 81. Planet. Space Sci. 23, 1649. J. geophys. Res. 96, 21,149.
1990 1971 1960 1982 1991 1980 1989 1966a 1966b 1971
Ann. Geophys. 8, 273. Radio Sci. 6, 65, J. geophys, Res. 65, 815. J. atmos, terr. Phys. 44, 901. J. atmos, terr. Phys. 53, 151. J. atmos, terr. Phys. 42, 983. J. atmos, terr. Phys. 51,249. J. atmos, terr. Phys. 28, 747. J. atmos, terr. Phys. 28, 807. Proc. R. Soc. Lond. A. 321, 69.
ADACHI S. BOOKER H. G.
The Propagation of Radio Waves. Cambridge Uni-
versity Press, Cambridge.
and RYCROFT M. J. HAMAR D., TARCSAIGY., LICHTENBERGERJ., SMITH A. J. and YEARBY K. H. KARPMAN V. I. and KAUFMANR. N. KARPMAN V. I. and KAUFMANR. N. KARPMAN V. I. and KAUFMANR. N. KARPMAN V. I. and KAUFMANR. N.
KARPMANV. I. and KAUFMANR. N. LAIRD M. J.
LAIRDM. J. and NUNND. LICHTENBERGER J., TARCSAIGY., P.~,SZTORS., FERENCZ C., HAMAR n . , MOLCHANOVO. A. and GOLYAVIN A. SAZHIN S., SMITH A. J. and SAZmNA E. M. SCARABUCCIR. R. and SMITH R. L. SMITHR. L., HELLIWELLR. A. and YABROFFI. W. STRANGEWAYSH. J. STRANGEWAYSn . J. STRANGEWAYSH. J. and RYCROFT M. J. TARCSAI GY., STRANGEWAYSH. J. and RYCROFT M. J.
WALKERA. D. M. WALKERA. D. M. WALKERA. D. M.
Reference is also made to the following unpublished material:
RYCROFTM. J. and ALEXANDERP.D.
1969
Paper A.3.11. presented at the 12th Plenary Meeting of COSPAR, Prague, Czechoslovakia.
Mode theory of whistler ducts APPENDIX
To evaluate the integrals I,., I, and/=, defined by (22), (23) and (24), put the rational part of the integrand into partial fractions. Hence
1~ = Y l o + 2 1 , - ( S ~ - l ) I 2 - ( l - s ~ ) I 3 , 1, = (S212+s[13)/(S~--s~)-lo, I= = ti1, + 13 - 12,
( S ~ - 1)12 = 9
~
,2 (1 - n t sin2 4)) d4)/ (1 - n 2 sin 2 4))(1 - m sin 2 4))1,2, (A8)
'n,2
(AI) (A2)
1607
(1 -s~)13 = 9
(A3)
t
(1 --n, sin24)) d4)/
0
(l-n3sin24))(1-msin24)) ''2.
lo= f[°ds/Q ''2,
,%
1,= f[sds/Q ~'2, II(n,m)=
I2= f f ds/(S~-s)Q 1/2, I3= f[°ds/(s-so~)Q':2.(A4) When Q is factored as in (25), the appropriate transformation is (BYRD and FRIEDMAN, 197l, p. 112 et seq.) sin 2 4) = (s-- 1)/n,(s + 1),
(A5)
d4)/(l-nsin24))(1-msin24)) ~,~2. (AI0)
For the complete ray path, the corresponding incomplete elliptic integrals are required. At upper cut-off, so = So = ~ Y, and Im is given by
~Im = fl Y/2 ( s 2 - l ) l ' e d s / ( S ~ - s ) ( s - s ~ ) ,
s~ < 1,
where n~ and other relevant parameters are defined by (27). One then finds that
n2
1o = 9
I do
d4)/(l - m sin 24)) 1,2 = 9K(m),
(A6)
where 9 = 2/(So-1)~Z2(So+ 1)~z2, and K(m) is the complete elliptic integral of the first kind. The other integrals are given by I~ = 9
(A9)
I~, 12 and 13 may then be expressed in terms of K(m) and complete elliptic integrals of the third kind
where
(l +n~ sin24)) d4)/ ( 1 - n l sin24))(1-rosin24)) I/2, (A7)
(All) which, as stated, is readily integrated in terms of elementary functions. At the critical value Y,, given by (9), s~ = 1, S~ = Y,.- 1, and the integral becomes ½1,.= fl Y''2 (s+ 1)"2(s- 1) L'2ds/(Y,.--1--s)
= y,!/2(y,.--2)-~/2 In { y,~/Z(y,,+2),.'2 + Yc+ l } - c o s h '~Y,., (= 26 ~/4In (1 + 2 ~/2)+0(6'/4)).
(A12)