Field line distribution of density at L=4.8 inferred from observations by CLUSTER

. For two events observed by the CLUSTER spacecraft, the ﬁeld line distribution of mass density ρ was inferred from Alfv´en wave harmonic frequencies and compared to the electron density n e from plasma wave data and the oxygen density n O + from the ion composition experiment. In one case, the average ion mass M ≡ ρ/n e was about 5 amu (28 October 2002), while in the other it was about 3 amu (10 September 2002). Both events occurred when the CLUSTER 1 (C1) spacecraft was in the plasmatrough. Nevertheless, the electron density n e was signiﬁcantly lower for the ﬁrst the distribution of electron density can also be inferred from measurements by multiple spacecraft.

ferred mass density decreases to a minimum 23% lower than the equatorial value at |MLAT|=15.5 • , and then steeply increases as one moves along the field line toward the ionosphere. For this event we were also able to examine the spatial dependence of the electron density using measurements of n e from all four CLUSTER spacecraft. Our analysis indicates that the density varies with L at L∼5 roughly like L −4 , and that n e is also locally peaked at the magnetic equator, but with a smaller peak. The value of n e reaches a density minimum about 6% lower than the equatorial value at |MLAT|=12.5 • , and then increases steeply at larger values of |MLAT|. This is to our knowledge the first evidence for a local peak in bulk electron density at the magnetic equator. Our results show that magnetoseismology can be a useful growth and propagation of electromagnetic ion cyclotron (EMIC) waves. The field line dependence of electron density affects the growth and propagation of magnetospheric whistler waves, hiss and chorus.
Direct measurement of mass density in the Earth's magnetosphere is difficult because it requires particle detectors with mass discrimination capable of measuring particles with low (∼eV) energies. Because of this, and also because of the potential for remote sensing, toroidal (azimuthally oscillating) Alfvén wave frequencies have increasingly been used to determine magnetospheric mass density both from the ground (Waters et al., 2006) and from space . We sometimes refer to this as magnetoseismology.
Assuming a power law form for the mass density, a number of authors have used toroidal Alfvén wave harmonic frequencies to infer the field line dependence of mass density (see references in Waters et al., 2006). Price et al. (1999) cast the wave equation into finite difference form in order to solve for the mass density at several locations along a field line. Denton et al. (2001Denton et al. ( , 2004 solved for the field line dependence of mass density using a polynomial expansion in a coordinate related to distance along the field line. These last two studies revealed an equatorial peak in mass density. Statistical studies Takahashi and Denton, 2007) have revealed that the tendency for the mass density to peak at the magnetic equator is greatest for large L 6 in the afternoon local time sector, especially during geomagnetically active times. We define L to be the maximum radius to any point on the field line (based on the TS05 model (Tsyganenko and Sitnov, 2005)) divided by R E . In dipole coordinates, this is the ordinary L shell.
For this paper we return to an event study. Although the statistical studies of toroidal Alfvén frequencies have the greatest potential for revealing the typical field line dependence and give some indication of the range of variability, it is not clear exactly how much of the variability in possible solutions is due to actual variation in the field line dependence at different times and how much is related to the uncertainty of the frequencies measured for individual events (Takahashi and Denton, 2007). (Takahashi and Denton ar-gued that both of these contribute to the range of possible solutions.) For instance, we can conclude from Takahashi and Denton's results that for L>6 in the afternoon local time sector, there is on average a peak in mass density at the magnetic equator, but we are not sure how much variability there is in the statistical distribution, that is, to what extent the distribution can be more or less peaked. For this reason, it is still of interest to look at individual events, particularly if the frequencies can be measured with great accuracy. In this study, we use toroidal Alfvén frequencies observed by the CLUS-TER spacecraft on 28 October 2002 in order to find what is probably the most precise field line distribution of mass density determined to date, one for which the errors in the inferred mass density are probably dominated by factors other than the uncertainty of the Alfvén frequencies (e.g., magnetic field model and theoretical wave equation).
The best technique for determining the field line dependence of electron density is probably the active sounding technique of Reinisch et al. (2004), but this technique has so far been useful for detecting the near-equatorial density distribution only in the plasmasphere, where electron density is relatively flat . Here we develop a least squares fitting technique for probing the distribution of electron density using plasma wave observations by the four CLUSTER spacecraft. For an event on 10 September 2002, we find (to our knowledge for the first time) an equatorial peak in bulk electron density at the same time that we find a (larger) equatorial peak in mass density inferred from Alfvén wave frequencies. In Sect. 2, we describe our method for inferring mass density; in Sect. 3, we describe our least squares fitting method for determining the spatial distribution of electron density; in Sect. 4, we describe our results for the field line distribution of the mass density on 28 October 2002; in Sect. 5, we describe our results for the field line distribution of mass and electron density for 10 September 2002; in Sect. 6, we compare the densities to other density measurements and compare the two events to each other; and in Sect. 7, we end with discussion.

Solving for the mass density
Our method for solving for the mass density ρ is described by Denton et al. (2004). In brief, we solve the Alfvén wave equation (Singer et al., 1981) with perfectly conducting ionospheric boundaries using the TS05 magnetic field model (Tsyganenko and Sitnov, 2005) and the following form for the base 10 logarithm of the mass density ρ log 10 ρ = c 0 + c 2 τ 2 + c 4 τ 4 + c 6 τ 6 + . . . , with up to seven terms (up to c 12 τ 12 ). The Alfvén crossing time coordinate τ is where ds is the differential length along the field line and V A is the Alfvén speed. The integral in Eq.
(2) is calculated from the magnetic equator (position of minimum B 0 ) to any position along the field line. The coordinate τ is the most natural choice for evaluating the mass density, since in the WKB approximation, the nodes of an Alfvén wave are evenly spaced with respect to this coordinate. Note that we only use even powers, assuming that the distribution of ρ symmetric about the magnetic equator. This assumption has been functionally required because the solutions are not well constrained if we allow the distribution to be asymmetric about the magnetic equator (see discussion by Denton et al., 2004); however, in this paper, we present evidence that the field line distribution of n e is symmetric about the magnetic equator, at least for the 10 September 2002 event. Because the coordinate τ is itself a function of the field line distribution for ρ that we are seeking, we need to assume a field line dependence for ρ for the purpose of defining τ . Here, we use a power law field line distribution, where R is the geocentric radius, and LR E is the largest value of R at any point on the field line (normally the magnetic equator); at the location corresponding to this largest value of R, the mass density is ρ 0 . The power law coefficient is found from the slope of the observed frequencies as described by Denton and Gallagher (2000). Please note carefully that the power law form here is used only to define the coordinate τ used by our method. The resulting field line distribution using Eq. (1) is not restricted to the form of Eq. (3). To find the coefficients c i yielding theoretical frequencies f th−n best matching the observed frequencies f obs−n , we minimize where n is the harmonic number (number of antinodes between the ionospheric boundaries for the electric field or velocity perturbation) and N freq is the number of harmonics observed. We first do a grid search over a range of values for the coefficients c i , and then use a root finding or minimization routine to zero in on the solution (Press et al., 1997). Alternately we may solve for the density using log 10 ρ = C 0 + C 2 z 2 + C 4 z 4 + C 6 z 6 + . . . , where z= √ 1−R/(LR E ), where R is the geocentric radius, and LR E is defined as the maximum distance to any point along the field line. In dipole coordinates, L is the usual L shell, and z is the sine of the dipole coordinate latitude MLAT. (It was not possible to use the coordinate τ used by Takahashi and Denton (2007) because the usefulness of that coordinate requires that the mass density not have a large variation along the field line, and we found that there was a large variation, especially for the event observed on 28 October 2002.) For the solution of the wave equation, we assume a perfectly conducting boundary at a height of 300 km (120 km) for the 28 October 2002 (10 September 2002) event. (There is no good reason for using different values, but the difference in results is so small that it's not worth redoing the calculations to make these values the same. For the 28 October 2002 event, using the 300 km height leads to values of ρ that are 0.9% higher at the magnetic equator, and 14% higher at a height of 300 km. See the discussion by CLUSTER spacecraft (where C1, C2, C3, and C4 are black, red, green, and blue, respectively) in a meridian with X=R cos(MLAT) and Z=R sin(MLAT). The dotted tracks show the spacecraft orbits from 11:40 UT to 12:07 UT on 10 September 2002, and the position of the spacecraft at 12:07 UT is indicated by the squares. The displacement in the out of plane (azimuthal) direction at 12:07 UT is indicated by the numerals next to the squares.

Method for detecting the distribution of electron density using multiple spacecraft
For the 10 September 2002 event, we have electron density measurements from all four spacecraft, so it should be possible in principle to infer the distribution of density. We started out using calculations of the density gradient (Darrouzet et al., 2006). However, tests revealed that the gradients were not accurate. In particular, if we created simulated density measurements having a dependence (LR E /R) α (5/L) β , where α is a power law coefficient for the parallel dependence (Eq. 3) and β is a power law coefficient for the L shell dependence, then we inferred the coefficients from the gradient based on simulated data, the inferred coefficients varied greatly from the coefficients of the simulated data. For instance, with α=0 and β=4 for the simulated data, we inferred from the density gradient α=2-4 (depending on exactly how we define the coefficient) and β=8 at L=5 with the effective β varying greatly with respect to L. We conclude based on this test that the gradient calculated from the four spacecraft does not accurately represent the actual gradient at the center of mass of the four spacecraft. The gradient method would work if the gradient were constant over the region of space occupied by the four spacecraft, but that was not the case for our data, and the spatial arrangement of the spacecraft was not at all tetrahedral, accentuating the problem. See Fig. 1 We then developed the following least squares fitting method to examine the distribution of electron density for the 10 September 2002 event. If several CLUSTER spacecraft pass through a region and we are able to determine the local electron density n e along the trajectories (this method could be used with any number of spacecraft, but several are needed to constrain the possible solutions), we are able to determine the spatial and temporal distribution of n e under the key assumptions of separability and smoothness. We assume that we can separate the dependence of n e into a product of terms n e,mod = n e,L n e, n e,MLT n e,UT , where n e,L is a function only of L, n e,MLT is a function only of MLT, n e, is a function only of a parallel coordinate (ideally orthogonal to L and MLT), and n e,UT is a function only of the time in hours. In addition, we assume that these separated dependencies are relatively smooth (first few derivatives continuous) so that we can describe them in terms of a superposition of a few simple functions (like polynomials). The problem becomes easier to solve numerically if we take the (natural) logarithm of Eq. (6), ln n e,mod = ln n e,L + ln n e, + ln n e,MLT + ln n e,UT .
Now suppose that we can write the individual logarithmic dependencies in terms of a sum of terms depending on the coordinate. For instance, for the L dependence, we may use ln n e,L = a L,0 + a L,1 ln(5/L d ), where L d is the L value based on a dipole field, that is, where R is the geocentric radius, MLAT (defined in Eq. 11) has values very close to those of the magnetic latitude MLAT, and the number "5" in Eq. (8) is an arbitrary choice. With (8), n e,L is exp ln n e,L , or showing that this is a power law dependence in terms of L d . The adjusted magnetic latitude MLAT is Using the TS05 magnetic field model to map the position of the CLUSTER spacecraft to the position of maximum geocentric radius, we find that this position is offset from MLAT=0 by MLAT 0 (1.63 • for the 10 September 2002 event), and we judge MLAT=MLAT 0 to be a better description of the location of the magnetic equator. Note that the measured values of n e observed by the CLUSTER spacecraft for the 10 September 2002 event peaked at this slightly offset value. Alternately, we may use a more general polynomial expansion, ln n e,L = a L,0 + a L,1 L d + a L,2 L 2 d + a L,3 L 3 d , where in most cases we drop the last (cubic) term. Keep in mind that the a terms are constants that multiply field values (e.g., L d ). With Eq. (12), n e,L is n e,L = exp a L,0 + a L,1 L d + a L,2 L 2 d + a L,3 L 3 d .
For the parallel dependence (along the field lines), we use either ln n e, = a ,1 ln (L d /R) , or ln n e, = a ,1 z + a ,2 z 2 + a ,3 z 3 + a ,4 z 4 +a ,5 z 5 + a ,6 z 6 , where only a limited number of these terms is typically kept. Equation (14) yields n e, = (L d /R) a L,1 , which is the power law form (3) for n e (rather than ρ) with α=a ,1 . For a dipole field, (L d /R) a L,1 =1/ cos 2 (MLAT), a function of MLAT; MLAT is not an orthogonal coordinate to L d , and neither is MLAT . (This is not a terrible problem, but it should be kept in mind that L d and MLAT are not entirely independent.) Equation (15)   of the event. In a number of cases, we use a term proportional to This combination takes into account corotation (for constant MLT , MLT and UT increase together) and leads to a better match between the observations and the model than does MLT alone if explicit UT dependence is not modeled. A full list of the possible terms used in the method appears in the middle range (between the second and third horizontal lines) of the first column of Table 1. (The specific fits will be discussed in Sect. 3.) In schematic form, we have a string of observations of n e , n e (j ), where j is a single subscript representing the variation in time and spacecraft, and we want to approximate the values of ln (n e (j )) by ln n e,mod (j ) = k a(k) t (j, k), where the index k represents one of the terms t (j, k) used in the model (one of the terms listed in Table 1). The coefficients a(k) are determined to minimize the averaged squared difference between the observed values, ln (n e (j )), and the model values, ln n e,mod (j ) . Since the model values come from a linear superposition, this (linear algebra) problem has a unique answer, and the standard difference can be determined, where N j is the number of j values (number of observations). Since the standard difference χ is that of the natural logarithm of the density, in order to get a measure of the agreement of the observed and model densities, we take exp (χ ). The resulting value gives a multiplicative error; the model values are likely to be in the range of the observation multiplied or divided by exp (χ ). If exp (χ ) is  close to unity, then exp (χ) −1 is the fractional error. Both χ and exp (χ) are listed at the bottom of Table 4.

Alfvén wave frequencies for 28 October 2002
Survey plots of power spectra of 4 s resolution magnetic field data from the CLUSTER Flux Gate Magnetometer (FGM) (Balogh et al., 2001) were scanned to look for toroidal Alfvén wave events. Clear events with dominant toroidal (azimuthal) polarization and with harmonics were not common, but there were a few. A Fourier spectrogram for one of these events observed by CLUSTER spacecraft 1 (C1) on 28 October 2001, is shown in Fig. 2. In order to detect toroidal Alfvén frequencies with sufficient accuracy, it is necessary to take a Fourier transform over a sufficiently long time at nearly constant L shell. For a spacecraft like CRRES, the best opportunity for determinations of mass density based on Alfvén frequencies is at spacecraft apogee . For the CLUSTER spacecraft, the orbit apogee is at geocentric radius R∼20 R E (Escoubet et al., 2001). At this distance, magnetic field models are not reliable. Probably the best location for magnetoseismology using CLUSTER data is at spacecraft perigee with R∼4-5 R E . On 28 October 2001, perigee occurred at about 02:50 UT, at which time the C1 spacecraft was 3.7 • north of the SM magnetic equator at magnetic local time MLT=8.7 (dawn). The middle panel of Fig. 2 shows the power spectrum of the toroidal (azimuthal) component of the magnetic field and bands of wave power are clearly evident at around 02:40-02:45 UT for the second harmonic (n=2 at 27 mHz marked by the magenta label "2") and 4th harmonic (n=4 at 55 mHz). The spacecraft crossed the SM magnetic equator at 02:43 UT, which is about the time with the clearest Alfvén wave bands. At earlier and later times (02:00 and 03:05 UT), wave power can be seen at the fundamental frequency (n=1 at 11 mHz). Weaker but still evident at 02:43 UT are other wave bands such as the 3rd harmonic (41 mHz), 5th harmonic (67 mHz), 6th harmonic (80 mHz), and even the 8th harmonic (103 mHz). The harmonic index n is the number of anti-nodes in the electric field or velocity perturbation between the ionospheric boundaries. Modes with odd harmonic number n, the fundamental, 3rd harmonic, 5th harmonic, etc. have a node (antinode) in the wave magnetic field (electric field) at the magnetic equator, while modes with even n, the 2nd harmonic, 4th harmonic, etc. have an antinode (node) in the wave magnetic field (electric field) . Because of this, the 2nd harmonic and 4th harmonic are the strongest modes seen near the magnetic equator.
In order to measure the frequencies of the harmonics, we use a semiautomated procedure. Figure 3 shows a grayscale plot of the power spectrum of the azimuthal component of the magnetic field near the magnetic equator (middle panel of Fig. 2). Here the darker shades correspond to higher wave power. The power spectrum was calculated from 4 s resolution data using 20 min windows with a Welch window (Press et al., 1997) and a moderate amount of frequency whitening (wave power at higher frequencies enhanced by multiplying by the frequency). Using an interactive program, we select wave bands by using mouse clicks (red diamonds in Fig. 3  Then the program finds the best linear fit to the power spectrum in the vicinity of the mouse clicks (red lines in Fig. 3). This is done by computing the integrated wave power (integrated with respect to time and frequency) around the fit weighting the power at each frequency in the power spectrum by a Gaussian in the difference between the wave frequency and the fit frequency with a Gaussian width of 5 mHz. Figure 4 is like Fig. 3, except for the power spectrum of the radial component of the electric field from the Electric-Field and Wave experiment (EFW) on Cluster spacecraft 1 (Gustafsson et al., 1997). The resolution of the data is the same as for the FGM instrument, and the power spectrum is calculated exactly the same way. In the case of the electric field, the strongest waves are observed for the fundamental and 3rd harmonic (because these harmonics have an anti-node in the electric field at the magnetic equator). The fit lines in Figs. 3 and 4 have a similar slope, and we normalized the frequencies to the second harmonic frequency to reduce the effect of changing frequency. These normalized frequency ratios (f n /f 2 ) are listed in Table 2 as a function of the harmonic number n (=1 for the fundamental) for the magnetic field data ((f n /f 2 ) B ) and electric field data ((f n /f 2 ) E ), where the uncertainties are found from the variation of the frequency ratio (f n /f 2 ) from the fit curves over the common frequency range that was fit for both f n and f 2 .
In the past, we have found the uncertainties of observed frequencies by finding the width of the frequency peaks (Denton et al., 2001Takahashi et al., 2004). However, this method probably overestimates the error because the error in the mean of a distribution is less than the width of the distribution (Takahashi and Denton, 2007). Unfortunately, our method of estimating the error listed in Table 2 (error from range of variation in the values of f n /f 2 from the linear fits) yields uncertainties that are too low. We can see this by comparing the values of (f n /f 2 ) B and (f n /f 2 ) E . For modes n=3, 4, and 5, the difference in these values is 0.015, 0.023, and 0.038. These values are larger than the uncertainties listed in Table 2. However, they are smaller than the uncertainties that would result from the frequency resolution of the 20 min time window used for the Fourier analysis (0.8 mHz). If the uncertainty of each harmonic frequency is assumed to be the same value f , the uncertainty of (f n /f 2 ), (f n /f 2 ) would be .031 (using f 2 =26.63 mHz), Eq. (21) yields (f n /f 2 )=0.057, 0.071, and 0.084 for n=3-5, which is larger than the difference between (f n /f 2 ) B and (f n /f 2 ) E (0.015, 0.023, and 0.038).
If we take f/f 2 =0.014, we get (f n /f 2 )=0.026, 0.032, and 0.038 for n=3-5, all of which are at least as great as (f n /f 2 ) B −(f n /f 2 ) E (0.015, 0.023, and 0.038). The uncertainties in Table 2 for (f n /f 2 ) are found using f/f 2 =0.014, using Eq. (21), but the values of (f n /f 2 ) for n=3-5 have been reduced by a factor of 1/ √ 2 to account for the fact that these values come from two measurements ((f n /f 2 ) B and (f n /f 2 ) E ). The unnormalized frequencies are labelled f n−input (mHz) in the next to last column of Table 2. These are found by multiplying (f n /f 2 ) by f 2 =26.63 mHz.

Inferred mass density for 28 October 2002
Given the frequencies of the toroidal Alfvén harmonics in Table 2, we can solve for the mass density distribution along the field line. Besides solving for the field line distribution of mass density consistent with the peak frequencies, we also do a Monte Carlo simulation for the mass density solution using Gaussian distributions of harmonic frequencies with mean and standard deviation equal to the mean values and www.ann-geophys.net/27/705/2009/ Ann. Geophys., 27, 705-724, 2009   uncertainties listed for f n−input in Table 2. In other words, we pick a random set of combinations of harmonic frequencies such that the distribution of the frequencies of each harmonic is consistent with the mean and uncertainty of f n−input for the harmonic number n. In order to determine the range of mass density distributions consistent with the frequency uncertainties, we solve for the mass density using 51 combinations of random frequencies. However, when using all eight frequencies as input to our mass density inversion code, and using the field line distribution (Eq. 1) with 5 polynomial terms (c 0 to c 8 ), our code did not converge on a solution for all combinations of the frequencies we tried. We only consider solutions valid if the theoretical frequencies found by the code are within the uncertainties of the original measurements from the input frequencies (which may include the random variation discussed immediately above). In fact, only 51 out of 209 frequency combinations attempted resulted in a valid solution. (We keep trying random combinations until we get 51 solutions.) The actual mean and standard deviation of the frequencies of the harmonics used in the Monte Carlo simulation (those frequency combinations yielding a valid solution) are listed as f n−solution in Table 2. As can be seen from Table 2, these values are within the range of frequencies listed for f n−input , but have a somewhat smaller standard deviation. The resulting solutions for mass density are displayed in the top panels of Fig. 5. The bold curve is the mass density based on the peak frequencies. There are three thin solid curves plotted in Fig. 5a and b. At every value of MLAT, the middle thin curve shows the mean value from the Monte Carlo simulation of random combinations of harmonic frequencies, while the upper and lower curves show the mean plus or minus one standard deviation. However, these four curves are so close together that it is difficult to see the variation. (Fig. 5c-f are discussed in Sect. 7.3.) In order to better show the variation in the solutions, we show in Fig. 6a the values of mass density from the Monte Carlo simulation divided by the solution for mass density using the peak frequencies (thick curves in Fig. 5a and b). As in our earlier solutions for the field line density of ρ (Denton et al., 2001Takahashi et al., 2004;Takahashi and Denton, 2007), the greatest variation in the solutions occurs at high latitude. However, the solutions here are staggeringly more precise than those from our earlier studies. We have a variation in the solutions of about a factor of 2 ( Fig. 6) with a variation of ρ with respect to MLAT of more than four orders of magnitude ( Fig. 5a and b). There are two reasons for the improvement. First of all, the relative uncertainty of the frequencies is about a factor of six smaller than in the study of Denton et al. (2004). The relative uncertainties are lower partly because of our better method for finding the frequencies (Figs. 3 and 4), but the main reason is that we were Ann. Geophys., 27, [705][706][707][708][709][710][711][712][713][714][715][716][717][718][719][720][721][722][723][724]2009 www.ann-geophys.net/27/705/2009/ able to use a comparable time window for the Fourier analysis (20 min versus 30 min used by Denton et al., 2004) while the frequencies were about a factor of 10 higher (because of the lower L=4.8 shell sampled at the perigee of CLUSTER compared to L∼7 sampled by CRRES). While τ may be in some sense the ideal coordinate (see Sect. 2), it is difficult to interpret the meaning of particular polynomial coefficients. For this reason, we have also solved for the field line distribution using a 5 polynomial fit for ρ with respect to the z coordinate. Figure 6b shows that this solution also lies within the range of solutions from the Monte Carlo simulation using the τ coordinate. This polynomial function is log 10 ρ = 1.57 + 0.35z 2 + 0.83z 4 + 3.95z 6 + 4.12z 8 , (22) for ρ expressed in amu/cm 3 .

Alfvén wave frequencies and inferred mass density for 10 September 2002
For our second event, observed at around 12:07 UT on 10 September 2002, the Alfvén wave data is not as complete or accurate, but this event is interesting because we have better data for n e and n O+ and because this event has a somewhat different field line distribution. This event occurred when the C1 spacecraft was again at perigee (geocentric radius R=4.8 R E ) at MLAT=−0.2 • and MLT=11.55 (noon). The Alfvén frequencies are found as described for the 28 October 2002 event and the frequencies and uncertainties are listed in Table 3 using magnetic and electric field data from the C1 spacecraft. As was the case in the first event, the uncertainties of the frequencies listed in the second and third columns found from the measured frequency ratios were judged to be too small. Based on frequency ratios found from the C2 spacecraft (for which the wave spectra were of lower quality; otherwise we would have used these to help calculate the harmonic frequencies), we used an error model like that used for the first event to find the uncertainties of the frequency ratios listed in the next-to-last column of Table 3, and those uncertainties were used in the following calculations.
Using a three polynomial fit with respect to the square of the Alfvén crossing coordinate τ (Eq. 2), we ran our mass density inversion code to get the field line distributions shown in Fig. 7a-d (black solid curves). There is a slight peak in ρ near the magnetic equator at MLAT =0 (MLAT =0 is at the position where the geocentric radius is maximum and the magnetic field is minimum as described in Sect. 3). The height of this peak is only 23% (as compared to a factor of 2 in the afternoon local time sector at large L 6 (Takahashi and Denton, 2007)). The distribution of solutions from the MLAT (˚) Here the values of ρ are divided by the mass density found using the peak frequencies in order to better show the difference between the solutions. (b) Same curves for ρ plus or minus one standard deviation (upper and lower curves). The middle solid (dotted) curve is the solution using the peak frequencies but with a 7 polynomial term fit with respect to τ (a 5 polynomial fit with respect to z). Both of these functions are also divided by the value of ρ found using the peak frequencies using the 5 polynomial fit with respect to τ . Monte Carlo simulation indicated by the spread of the upper and lower black thin solid curves in Fig. 7b does not exclude the possibility that the real distribution is flat. The solution using three polynomial terms with respect to z is similar to that using τ and is plotted as the black dotted curves in Fig. 7a-d. This solution is for ρ expressed in amu/cm 3 . (Fig. 7e-h will be discussed in the Discussion section.)

Spatial electron density distribution for 10 September 2002 found from data from all four spacecraft
For this event, there is a strong band of upper hybrid emissions observed by all four CLUSTER spacecraft, and the electron density n e can be easily determined. Because of this, we are able to estimate the distribution of the electron density using the method of Sect. 3. Figure 8 shows the n e values and position for all four spacecraft using the definitions in Sect. 3 (see also Fig. 1). For the purposes of using the method described in Sect. 3, we used the n e data shown  in Fig. 8a but excluded the dotted portion of the green curve (C3). The reason for excluding this segment of data is that only C3 sampled the large L d and large negative MLAT values. The low n e values of C3 could be modeled by alteration of either the terms describing the L d or MLAT dependence, and which of these was affected was not well determined, leading to unnecessary variation in the solutions depending on the exact terms included in the model. Aside from C3, for which we used a start time of 12.1 UT (h), the data was used between 11.75 and 12.7 UT (times plotted in Fig. 8). These limits were chosen because of large oscillations in the C3 n e (plasma frequency) at earlier times and because C4 penetrated into the plasmasphere at later times (not shown). However, a more general treatment would allow different time limits for each spacecraft.
As described in Sect. 3, we find an optimal fit to the observed values of ln(n e ) using a superposition of particular model functions. The models we used are listed in Table 1. For the first model, fit 1, (using the L shell dependence in Eq. (8) and the parallel dependence in Eq. 14) where we have neglected the (small) MLAT dependence in the last line and the power law coefficients for the parallel and L d directions are α=1.81 and β=4.21. This distribution of density is extremely reasonable. Denton et al. (2004) state that α=2-3 is typical in the plasmatrough. Values of β equal to 4.5 (Carpenter and Anderson, 1992), 4.0 (Sheeley et al., 2001), and 3.5  have been found in statistical studies of plasmatrough density. The value β=4 results from the assumption that the flux tube content per magnetic flux is constant . The standard difference (of the natural logarithm) is χ=0.0845, leading to a multiplicative error factor exp(χ )=1.088 (approximately a 9% error in the linear density). In Fig. 9a, we plot the observed n e values (solid curves) and model values for fit 1 (dashed curves) for all four CLUS-TER spacecraft. The model curves fit the data fairly well. We noticed, however, that the model curves in Fig. 9a have lower density than the observed n e at the places where the spacecraft cross MLAT =0. These locations are indicated in Fig. 8 as vertical dashed lines for C1 (black) and C2 (red). In addition, C4 reaches MLAT =0 at the right side of the plot. Note the local increase in n e at these positions. This observation led us to explore sets of model functions with a more general field line dependence. The second set of model functions for the natural logarithm of n e (fit 2 in Table 1) includes terms proportional to z 2 and z 4 (parallel dependence (Eq. 15)) instead of the power law field line dependence. (Here, symmetry about MLAT =z =0 is assumed; see definitions in Sect. 3.) Fit 2 also has a quadratic L d dependence (constant, L d , and L 2 d terms). For fit 2, χ =0.0504 or exp(χ )=0.0504 (roughly a 5% error), a significant improvement. The model functions for fit 2 are shown in Fig. 9b (dashed curves). The field line dependence n e, (see Eq. 6) for fits 1 and 2 is shown in Fig. 10. Fit 2 yields a local peak in n e at MLAT =0 with a drop of n e of about 8% to MLAT ∼13 • before n e starts to increase again at larger values of |MLAT |. It should be kept in mind that the parallel dependence inferred here is derived from and therefore representative of the densities observed by all four spacecraft, and may not exactly correspond to the exact parallel dependence on any particular field line; we expect that it represents a typical or average field line dependence within the region sampled by the four spacecraft.
The L d dependence of n e , n e,L d , normalized to its value at L d =5 is shown in Fig. 11 for fit 1 (solid black curve) and for fit 2 (dashed blue curve). Both of these show a monotonic decrease in n e,L d with respect to L d , but the quadratic fit 2 is slightly less steep. Motivated by the better agreement of fit 2, we tried fit 3 with a cubic dependence in L d and z 2 , z 4 , and z 6 terms for the field line dependence. In this case, however, there was not a significant improvement in the fit. The value of χ was 0.0490 (versus 0.0504 for fit 2) and exp(χ ) was 1.050 (versus 1.052 for fit 2). In Fig. 11, n e,L d for fit 3 is plotted as the dashed red curve, and it is not significantly different from the blue dashed curve (fit 2). Based on this comparison, we used only a quadratic fit for L d and terms only up to z 4 for the remaining models.
The remaining models were partly motivated by the disagreement between the model and observed curves for C2 (red curves) and C4 (blue curves) in Fig. 9b. Noting that these differences occur at the right side of the plot where the spacecraft are moving to larger values of MLAT , we tried a model with an asymmetric distribution with respect to MLAT . Fit 4 included terms for z , z 2 , z 3 , and z 4 . The red solid curve in Fig. 10 Fig. 7. (a, b) Mass density ρ (black curves) and electron density n e (red curves) versus MLAT (left) and R (right). The black solid curves are all found using a three term polynomial expansion for ρ with respect to τ . The thick black solid curve is the solution based on the peak frequencies in Table 3, the middle thin black solid curve (obscured in (a, b) by the thick black solid curve) is the mean solution from the Monte Carlo simulation using the uncertainties listed in Table 3, and the upper and lower thin black solid curves are the mean solution plus or minus one standard deviation (in the log value). The black dotted curves are found using a three polynomial expansion with respect to z. The red curve is n e found from fit 5 as described in the text. ment is achieved with fit 5, which uses only even terms with respect to z (z 2 and z 4 ), but includes terms proportional to both MLT and UT (rather than MLT ≡MLT −UT for fit 4). (As described in Sect. 3, the primes for MLT and UT merely indicate that the average value has been subtracted.) In this case, χ and exp(χ ) are 0.0431 and 1.044, respectively, versus 0.0451 and 1.046, respectively, for fit 4 (Table 1).  Fig. 9. Observed values of n e (solid curves) and model values (dashed curves) versus UT for all four CLUSTER spacecraft (black, red, green, and blue color for C1, C2, C3, and C4, respectively) for (a) fit 1, (b) fit 2, and (c) fit 5 in Table 1. MLAT ' e, n Fig. 10. Parallel dependence of n e , n e, , versus the magnetic latitude MLAT for fits with n e, symmetric with respect to MLAT (black curves): fit 1 (solid black curve), fit 2 (dotted black curve), and fit 5 (dashed black curve); and for asymmetric (with respect to MLAT ) dependencies (red curves): fit 4 (solid red curve), and fit 6 (dashed red curve).
In Fig. 9c, we compare the observed and model densities for fit 5. The agreement is good in most regions; there is still some disagreement for C2 (red curves) at the largest values of UT. This remaining disagreement is either because we have not used general enough model functions, or the assumption of separability in Eq. (6)  we regard the fit shown in Fig. 9c to be very good. The L d dependence for fit 5 is the dashed green curve in Fig. 11, and it's clearly very close to the curves for fits 2 and 3 (blue and red dashed curves, respectively). The parallel dependence for fit 5 is the black dashed curve in Fig. 10. There is still a local peak in n e at MLAT =0, but the peak has slightly less amplitude (6% drop in n e from the value at MLAT =0 to the minimum value versus a 9% drop for fit 2). Fit 6 is the same as fit 5 (with terms for MLT and UT ) but with asymmetric terms in z allowed like for fit 4. That is, terms for z , z 2 , z 3 , and z 4 were all included. The χ and exp(χ ) values, 0.0430 and 1.044, respectively, were almost identical to those for fit 5 (0.0431 and 1.044, respectively), as shown in Table 1. The coefficients of the odd terms in z (z and z 3 ) were very small compared to those of the even terms, so that the field line dependence was almost the same, as can be seen by comparing the red dashed curve (fit 6) to the black dashed curve (fit 5) in Fig. 10. While the value of χ is not greatly smaller for fit 5 than for fit 4 (0.0431 for fit 5 as compared to 0.0451 for fit 4), we believe that fit 5 better represents the distribution of n e for several reasons. First of all, fit 5 has one less fit parameter than fit 4, yet achieves a better fit. Secondly, it is hard to understand physically why the distribution of density with respect to MLAT would be asymmetric, especially for 10 October, a date near the equinox, for which the ionospheric conditions should be similar at both ends of the field line. Third, fit 5 yields a more reasonable dependence for MLT. Fit 5 has a positive coefficient for the term proportional to MLT , indicating that the density increases with respect to MLT . Fit 4, on the other hand, has a negative coefficient for MLT ≡MLT − UT , indicating that at UT =0, the density decreases with respect to MLT. Seeing as the position is near local noon (Table 4), an increase in n e with respect to MLT is expected since n e generally peaks in the afternoon local time sector. Finally, when both explicit MLT and UT terms and odd terms in z are included (fit 6), the resulting fit is nearly identical to that of fit 5. Now we return to Fig. 7, where we make use of the parallel distribution of n e from fit 5. In Fig. 7a and b, the red curve is the product of the L d dependence of fit 5 evaluated at L d =4.8 (the L value of C1 at the time the Alfvén frequencies were observed) and the parallel dependence of fit 5, n 5, . In Fig. 7c and d, the same values are plotted but multiplied by 2.5 (to aid comparison to the curves for ρ). Fit 5 for n e has a local peak in n e at the magnetic equator, like ρ, but the amplitude of the peak in n e is smaller, a 6% drop in n e from MLAT =0 to the minimum value versus a 23% drop in ρ for the solution based on the peak frequencies (thick solid black curve in Fig. 7a). Using this last curve for ρ with the red curve for n e , we can derive the average ion mass M≡ρ/n e within the region of overlap in MLAT ; M is plotted in Fig. 7e and f. The value of M is largest at MLAT =0.
Some earlier work has indicated that the He+ density, which is usually small, is not nearly as sensitive to geomagnetic activity as is that of O+ (Craven et al., 1997;Krall et al., 2007). This suggests that if M is significantly greater than unity, n O+ /n e is as a first approximation equal to (M−1)/15. Combining this relation with quasineutrality for an H+/O+ plasma (1=n H+ /n e +n O+ /n e ), we can use the field line distribution of M to derive the field line distribution of H+ and O+ and these curves are plotted in Fig. 7g and h. The field line distribution of H+ is less peaked than that of n e (3% peak for H+ versus a 6% peak for n e ), while the field line distribution of O+ is more peaked than that of ρ (34% peak for O+ versus a 23% peak for ρ).

Comparison to other density measurements
The purpose of this section is to check for consistency of the measured density values, and then to compare the two events. For consistency, the average ion mass M≡ρ/n e should be greater or equal to 1 amu (corresponding to pure H+ plasma), and the directly measured heavy ions should be less than that needed to account for the mass density inferred from the Alfvén wave frequencies.

Other measurements for 28 October 2002
We were able to estimate the electron density for the C1 spacecraft at 02:48 UT on 28 October 2002 (when the Alfvén frequencies were measured) using the WHISPER instrument sounding data (Rauch et al., 2005). Natural plasma wave emissions were at a high level, saturating the receivers. Active sounder emissions are shown in Fig. 12 smoothed spectrum. The smoothing uses a common ad hoc filter function, providing in all regions a first approximation of the plasma frequency Fp. It is based on the fact that all resonance peaks (Fce, Fqs, Fp, Fuh) appear in a frequency range more or less centered near Fp. A value of Fp between 23 and 27.2 kHz leads to an electron density n e of about 8 cm −3 . Oxygen ions (O + ) can be detected using the Cluster Ion Spectrometry (CIS) instrument (Rème et al., 2001). A spectrogram of the flux of O + ions with energy from 1-40 keV is shown in the top panel of Fig. 13, and the density of O + integrated over the same energy range is shown in the bottom panel. Below 1 keV, the signal to noise ratio was too low to get a reasonable estimate. At about 02:30 UT, the measured O + density is about 0.5 cm −3 . Table 4 summarizes our information about the 28 October 2002 event. The negative D st value = −49 nT is indicative of a significant ring current population. The instantaneous K p and K p 3 (average of previous K p values weighted with exp(−(t−t Kp )/(3 days))) values (4.6 and 4.1, respectively) indicate a moderately high level of geomagnetic fluctuations. The C1 spacecraft was close to the magnetic equator (MLAT = −5.4 • ), and the inferred ρ at the spacecraft location was 37.6 amu/cm 3 (Fig. 17c and Eq. 22). The shape of the field line distribution was flat in the middle region (around MLAT=0) and very steep at large |MLAT| (close to the ionosphere), as shown in Fig. 16. Using ρ=37.6 amu/cm 3 and n e =8 cm −3 , the average ion mass M≡ρ/n e =4.7 amu. In a plasma composed of H+, He+ and O+, the average ion mass will be M = 1 + 3 n He+ n e + 15 n O+ n e .
To have a value of M=4.7 would require n O+ /n e =0.25 in an H+/O+ plasma, or n O+ =2 cm −3 . If n He+ /n e =0.20 (a typical upper limit), n O+ /n e =0.21 is needed, corresponding to n O+ =1.7 cm −3 . We observed n O+ =0.5 cm −3 with energy >1 keV, and the plot of the flux of O+ (bottom panel of Fig. 13) indicates that there is more O+ at lower energies. Thus we observed a substantial fraction of the needed O+, enough to raise M to the value M O+,obs =1.9 in Table 4, but not enough to raise it to M=4.7.

Other density measurements for 10 September 2002
At the time of this event, n e =22 cm −3 at the position of C1. From the solution based on the peak frequencies (thick solid curve) plotted in Fig. 7, ρ=63.0 amu/cm 3 there, leading to M=2.8. In this case, we were able to measure O+ down to 40 eV (Fig. 14), and find an integrated oxygen density n O+ =0.6 cm −3 . The measured O+ is by itself enough to raise M to a value of 1.4. These values are listed in Table 4. As in the previous case, there is a significant amount of flux at the lower range of energy, so it is likely that there was more oxygen present. From Eq. (25), a value of M=2.8 would require n O+ /n e =0.12 in an H+/O+ plasma, or n O+ =2.7 cm −3 (as compared to 1.7-2 for the 28 October 2002 event). With n He+ /n e =0.2 (a typical upper limit), there would need to be n O+ /n e =0.08, or n O+ =1.8 cm −3 (as compared to 1.7-2 for the 28 October 2002 event). Based on these numbers, the total amount of O+ could be similar for the two events. As can be seen in Fig. 14, there are two populations of O+, one with energy lower than 1 keV, and another with higher energy. About 0.5 cm −3 of the total O+ (0.6 cm −3 ) is in the higher energy range, and this is the same amount measured for the 28 October 2002 event (Table 4). Based on the pitch angle distribution (not shown), the higher energy population is trapped (distribution peaked around 90 • pitch angles), but the lower energy population is field aligned (distribution peaked around 0 and 180 • pitch angles).

Comparison of the two events
The density of ring current (>1 keV) O+ is about the same for both cases and the total amount of O+ could be similar, but the electron density is larger for the 10 September 2002 event (22 cm −3 as compared to 8 cm −3 for 28 October 2002), leading to a smaller value of M for the 10 September 2002 event (2.8 as compared to 4.7 for the 28 October 2002 event). The D st values for the two events are similar (Table 4), consistent with similar values for the O+ density with energy >1 keV. The K p value for the 10 September 2002 event is somewhat lower (3.6 as compared to 4.6 for the 28 October 2002 event). There is a greater difference in K p 3 (average of K p weighted with 3 day exponential decay as described in Sect. 3.3). The indices D st and K p are plotted for the two events versus time in Fig. 15, showing that K p was consistently higher during the two days leading up to the 28 October 2002 event. This probably correlates with greater convection leading to lower n e even though for both events, the CLUSTER spacecraft are in the plasmatrough (based on IM-AGE EUV images not shown). In summary, the ring current populations of the two events are similar (at least for O+), but the cold electron density is larger for the 10 September 2002 event, leading to a lower average ion mass M.

The difference in the two events
Before discussing the results for the field line distribution of density, we discuss the results presented in Sect. 6. Ideally, we would like to show that the value of mass density ρ inferred from Alfvén wave frequencies is equal to the sum of measured ion densities. Unfortunately, we cannot show that because the ion composition instrument on CLUSTER does not measure the cold (∼eV energy) heavy ions (ions heavier than H+). (Given the heavy ion densities, we can infer the proton density from the electron density n e .) We have shown that the measured O+ density n O+,obs leads to an average ion mass that is significantly greater than unity, yet still less than the value implied by ρ. and ρ are at least possibly correct. In the two events considered (Table 4), both of which were in the plasmatrough, the D st values and density of O+ with energy >1 keV were similar (see discussion in Sect. 6.2), and the total amount of O+ may also possibly be similar. The main difference between the two events seems to be that the total density was greater for 10 September 2002 than for 28 October 2002 (Table 4).
Although the field line distribution is slightly peaked for 10 September 2002, both field line distributions are relatively flat compared to field line distributions studied by  (23% peak for 10 September 2002 versus about a factor of 2 for L=6-8 in Fig. 8 of Denton et al.). Denton et al. found that the field line distribution of mass density was typically flat for L=4-5, but slightly peaked for L=5-6; the value of L in this study, 4.8, is close to the boundary between these two ranges.

Variation in solutions for the mass density for 28 October 2002
Here we consider how the solution for the field line dependence varies if we use different assumptions for the field line dependence or if we use a smaller number of harmonics. In Fig. 16, we show the solutions for ρ using different numbers of polynomial terms N poly in Eq. (1). Solutions with N poly <7 are shifted down by factors of 10 so that the solutions do not overlap, as indicated in the caption. The three solid curves represent solutions using all the observed frequencies N freq =8, but the dotted curves are solutions using N freq =N poly . (The vertical lines are there to help the viewer more accurately compare the different curves.) There are no solid curves for Fig. 16 with N poly ≤4 because when using 8 input frequencies (N freq =8) our code only converged on valid solutions when N poly was at least 5. This is probably because the field line distribution of ρ is very flat for  |MLAT| for |MLAT| 40 • , and this distribution cannot be fit with a small number of polynomials. The solid curves for N poly =5 and 6 are nearly identical. The solid curve for N poly =7 is very similar but increases less at the largest values of |MLAT|. Despite this, the solution for N poly =7 still lies within the range of solutions for the Monte Carlo simulation with N poly =5, as can be seen from Fig. 6b. While there are no solutions for N poly ≤4 using all eight input frequencies (Fig. 16), we can find solutions if we reduce N freq . The dotted curves in Fig. 16 Fig. 17. The bold solid curve is our solution with N poly =5 (τ expansion) with N freq =8. The thin solid curve is our solution using N poly =3 with N freq =3. The dashed curves are found using only the fundamental frequency (N freq =1), and assuming the power law form (Eq. 3) with α=1 (dotted), 2 (small dashes), 4 (medium dashes), and 6 (large dashes), as labelled in Fig. 17b. Using N freq =3, the three term polynomial solution does a fairly good job of modelling the field line distribution up to |MLAT|=35 • , but does not model the steep increase at large |MLAT| (at ionospheric altitudes) represented in the 5 polynomial solution. Clearly the power law form does not well describe the variation of ρ over the entire range of MLAT either. The field line distribution of ρ is very flat within |MLAT|=20 • (thick solid curve), and the best power law fit for this region is with α=1 (dotted curve). The value α=6 (large dashes) yields best agreement at large |MLAT| (close to the ionosphere). If three frequencies are used to determine the best power law fit, the inferred power law coefficient is α=2.25 and the resulting solution (not shown) is very close to the α=2 curve in Fig. 17 (found using only the fundamental frequency). As stated previously, our code does not converge on a valid solution using N freq =8 with the uncertainties in Table 2. However, if we increase the uncertainties by a factor of 10, we can find the best fitting power law solution, which has α=4.05; the resulting solution (not shown) is very close to the α = 4 solution in Fig. 17. , using CRRES data, suggested that α=2 was the best choice for L=4-5 if the power law form is used. These results also show that α=2 would work fairly well for the 28 October 2002 event up to about |MLAT|=35 • . The use of a power law solution with α = 4 better represents the values of ρ at larger |MLAT|, but as can be seen from Fig. 17, this solution leads to equatorial and topside ionosphere values of ρ that are too low, and mid-range values that are too high (comparing the α=4 medium dashed curve with the thick solid curve). Fig. 17. The bold solid curve is our solution for ρ versus MLAT (left) and geocentric radius R (right) using 5 polynomial terms to fit the 8 peak frequencies in Table 2. The thin solid curve is the solution using 3 polynomial terms to fit 3 frequencies. The four dashed curves are the solutions using only one frequency (the fundamental) assuming the power law form (Eq. 3) for α=1 (dotted), 2 (small dashes), 4 (medium dashes), and 6 (large dashes), as labelled in (b).
The bottom panels are the same as the top panels, except that the range of ρ values plotted is smaller.

Discussion of the mass density field line distribution for 28 October 2002
For the 28 October 2002 event, we were able to determine a large number of harmonic frequencies with unprecedented precision. The resulting solution for mass density was also very precise because of the large number of harmonics measured and the smaller relative uncertainties in the harmonic frequencies (Fig. 6). The net result is that the error in the inferred mass density is probably dominated by factors other than the uncertainty in frequency (e.g., magnetic field model and theoretical wave equation). As shown in Fig. 6, our solution for the field line distribution on 28 October 2002 (at L=4.8) is very flat for |MLAT| 20 • but very steeply increasing with respect to |MLAT| for |MLAT| 40 • . Because of the improved precision of the observed frequencies, we are able to see a steep increase in ρ as MLAT approaches ionospheric values ( Fig. 5a and b). And because of this increase, the Alfvén speed V A ≡B/ √ 4πρ (in CGS units) does not keep increasing as |MLAT| increases; the increase in ρ at large |MLAT| causes a decrease in V A at the largest values of |MLAT| as shown in the middle panels of Fig. 5. (Thus we observe the ionospheric Alfvén resonator (Polyakov, 1976;Lysak, 1993), a dip in the value of V A at low altitudes.) Because V A does not increase greatly as |MLAT| increases (less than a factor of 4 from Fig. 5c and d), the mass density at all values of ρ significantly Fig. 18. The solid curves are the same solutions for ρ that were plotted in Fig. 5a and b, but shown with respect to height h from the Earths surface within the southern (a) and northern (b) ionospheres. The large dashed curve is the mass density from the International Reference Ionosphere (IRI) model (Bilitza, 2001). The dotted curve is the mass density from the IRI due to O + .
affects the solutions for the frequency of the Alfvén wave harmonics. This can be seen from the values of τ (the Alfvén crossing time coordinate defined in Eq. 2) shown in Fig. 5e and f. In the WKB approximation, the Alfvén frequency ∼1/( N S ds/V A )=τ N −τ S , where the integral is evaluated from the southern ionosphere (S) to the northern ionosphere (N). The contribution to τ N −τ S in each differential change in latitude dMLAT goes like the slope of τ in Fig. 5e, and the slope of τ versus MLAT is non-negligible at every value of MLAT, showing that all regions of MLAT contribute to the Alfvén frequencies. This is the second reason for the improved precision at large values of |MLAT|; the values of ρ at large |MLAT| are having a significant effect on the frequencies. (In our previous solutions (e.g., Denton et al., 2004), ρ did not increase nearly so much at large |MLAT| so that V A became very large at large |MLAT| (due to the large value of B at low altitude) and the differential contribution to τ , dτ =ds/V A , became small.) Considering the tremendous imprecision of our previous solutions at large values of |MLAT| (Denton et al., 2001Takahashi et al., 2004;Takahashi and Denton, 2007), one might well question whether our technique can actually determine a variation in ρ of more than four orders of magnitude. While we cannot prove that this field line distribution of ρ is accurate, several pieces of evidence argue that the values of ρ may actually increase steeply as we have found. First of all, we do not find valid solutions for the field line distribution of log 10 ρ using less than five polyno-mial terms. (We also do not find a solution using the power law description with ρ∼R −α .) Secondly, polynomial fits to ρ (rather than log 10 ρ) did not yield valid solutions, indicating that we need a functional form that can represent a very steep increase in ρ at large values of |MLAT|. Third, the distribution seems to have converged with respect to N poly for N poly =5. Fourth, the values of ρ that we find at the largest values of |MLAT| compare reasonably well with ionospheric mass densities. This is shown in Fig. 18 for the southern (a) and northern (b) ionospheres. Note that our values of ρ are intermediate between the IRI values at height h 600 km and those at h 1400 km. Therefore, our estimate of mass density based on Alfvén frequencies may be probing into the ionosphere while not resolving its detailed structure.
Further work should be done to verify whether the four order of magnitude variation in ρ with respect to MLAT is justified based on our method, but it is encouraging that the ionospheric values of ρ are consistent with values from the IRI model. Such a dependence with increasing steepness away from the magnetic equator is also consistent with the field line dependence of n e inferred from active radio sounding by IMAGE RPI (Reinisch et al., 2004;Denton et al., 2006, and references therein). An important effect of the steep increase in ρ at large |MLAT| is that the large |MLAT| (low altitude) portion of the field line contributes to the global frequency of oscillation, even for the fundamental mode (Fig. 5e).

Discussion of the field line distribution of mass density for 10 September 2002
For the second event, 10 September 2002, we measured only three harmonic frequencies and the solution was not as accurate (partly because of the small number of Alfvén wave harmonics observed). In this case, we found a small local peak in mass density near the magnetic equator. This peak is similar to those we have found in earlier studies, but smaller in amplitude, 23% (as compared to a factor of 2 in the afternoon local time sector at large L 6 (Takahashi and Denton, 2007)). Because this peak is small, and because the range of solutions from the Monte Carlo simulation does not exclude a flat distribution (Fig. 7), we would be reluctant to confidently state that this peak is real based on this data alone. However, our analysis of the field line distribution of electron density for this event also indicates that the density is peaked at the magnetic equator.

Discussion of the electron density distribution for 10 September 2002
We first tried to determine the spatial distribution of n e using the gradient calculated from the measurements by the four CLUSTER spacecraft, as described by Darrouzet et al. (2006). However, tests of the method using simulated data with a known spatial distribution at the positions of the real spacecraft on 10 September 2002 showed that the calculated Ann. Geophys., 27, 705-724, 2009 www.ann-geophys.net/27/705/2009/ gradient did not accurately represent the actual gradient of the known distribution at the center of mass of the four spacecraft. Because of this, we developed a least squares fitting technique that can infer the density distribution under the key assumptions of separability and smoothness (Sect. 3). Unfortunately, the results depend somewhat on the exact model used. Nevertheless, tests using a number of models give a good indication of the terms that are necessary or not necessary. In particular, the model associated with fit 5 in Table 1 appeared to do a good job of describing the n e distribution. This fit modeled the natural logarithm of n e with a superposition of terms including a quadratic dependence in the dipole L value, a quadratic dependence with respect to the square of the sine of the magnetic latitude (z 2 and z 4 terms; see Sect. 3), and a linear dependence with respect to MLT and UT. Tests with other models showed that a cubic dependence with respect to L and odd terms with respect to MLAT were not necessary.
The model fits yielded the following results: there is a small rate of increase of the density with respect to both MLT and UT. The L d dependence, while not exactly described by a power law, is roughly consistent with L −4 d at L d ∼4.8 (see discussion immediately following Eq. 24). The field line dependence has a local peak at the magnetic equator (MLAT =0), decreases to a minimum value at a magnetic latitude |MLAT | of about 12.5 • , then increases steeply away from the magnetic equator (Fig. 10). Besides the fact that this dependence results from our model fits, it is noteworthy that we can see peaks in the raw n e data as the CLUSTER spacecraft cross MLAT =0 (Fig. 8a). To our knowledge, this is the first indication that the electron density in the plasmatrough may be locally peaked at the magnetic equator.
There is furthermore a remarkable similarity between the inferred field line distribution of n e based on the method described in Sect. 3 and that of ρ based on Alfvén wave frequencies (Fig. 7). The inferred distribution of ρ also has a local peak at the magnetic equator, decreasing to a minimum value at about 15 • . We brazenly used the inferred distribution of ρ and n e to calculate the field line distribution of H+ and O+. Considering the possible errors in both ρ and n e , it is hard to know at this point how seriously to take these distributions. Nevertheless, they are consistent with our expectations for a distribution of O+ gravitationally trapped at the magnetic equator (K. Ferriere, private communication, 2005). Given the same temperature, O+ would have a much smaller thermal speed than H+ and would be more strongly affected by the centrifugal force. Therefore, the H+ distribution might be relatively flat at MLAT =0 while the distribution of O+ is peaked there. The electrons are peaked just enough to bring about quasineutrality (Fig. 7).  plotted the log average value of n e from the CR-RES spacecraft versus MLAT (their Fig. 5). They did not emphasize the small peak near MLAT=0 because that peak was much less than what they found for ρ. Here we find the same result for 10 September 2002: Both ρ and n e show a local peak in density at the magnetic equator, but the peak is stronger for ρ. On the other hand, for the 28 October 2002 event we do not have any evidence that the density is peaked.

Conclusions
These results show some promise for using CLUSTER data to infer the distribution of both mass density and electron density near perigee. Unfortunately, a quick examination of a number of magnetic field spectrograms showed that Alfvén harmonic frequencies are not often clearly visible. What is most needed is an automated method of determining toroidal Alfvén frequencies using both magnetic and electric field data and using the polarization information (e.g., ratio of magnetic to electric power ) to identify the harmonics. Semi-automatic techniques are being used to build a database of CLUSTER n e measurements and some of these are available on the CLUSTER Active Archive (CAA).