Origin of some anisotropic tailward flows in the plasma sheet

We use a test particle model to explore anisotropy and fast flows in the central plasma sheet (CPS) that are a consequence of plasma sheet boundary layer (PSBL) ion beam dynamics. Ion distributions and flows (velocity moments) in the CPS and equatorial current sheet (CS) are compared and we find that mirroring of initially earthward beams from the PSBL, and their subsequent convection to the CS region, results in strong anisotropy throughout the CPS. At higher latitudes, velocity moments are field-aligned and feature earthward flow. Deeper in the CPS, velocity moments yield flows in the anti-earthward direction. There is no clear distinction between the PSBL and CPS, since velocity distributions with large streaming components occur throughout the model CPS, but in the CS region they are anisotropic and nongyrotropic. In the CS region velocity moments can feature anti-earthward cross field flows. These tailward flows (> 400 km/s) are observed in the CS region betweenX = −20 to −30 RE due to nonadiabatic effects. Model results suggest that fast tailward plasma flows can be obtained without necessarily appealing to magnetotail processes associated with dynamic geomagnetic activity.


Introduction
While, to the first order, the plasma sheet (PS) comprises a warm ion population with bulk flows primarily dominated by large scale E × B convection, careful inspection of distribution functions reveals many anisotropic flows encompassing a wide range of time scales.Early studies by Huang andFrank (1986, 1987), based on distribution functions from data acquired over 2 min by the ISEE 1 satellite between −10 Correspondence to: J. A. Wanliss (james.wanliss@erau.edu) to −20 R E , suggested flows were relatively slow (∼50 -100 km/s) and predominantly earthward.In another study, Paterson et al. (1998) surveyed the PS between 10 and 50 R E and found that flows perpendicular to the magnetic field were slow, and on the whole consistent with plasma convection.Another recent study by Borovsky et al. (1997) found that at low velocities (< 250 km/s) PS turbulence appears to be isotropic.
Unlike the slower flows, fast flows above 400 km/s appear more structured.Cattell and Elphic (1987) noted that ISEE 2 detectors were able to find much faster flows in the central plasma sheet (CPS) and subsequently, using distribution functions from AMPTE data acquired over 5 s, Baumjohann et al. (1989) reported that almost all fast flows inside the apogee of AMPTE IRM (viz.∼19 R E ) were earthward.Angelopoulos et al. (1992) were able to detect short-lived bursts of ion flow (> 400 km/s) which they termed bursty bulk flows (BBFs).Fast ion flows in the CPS have been studied extensively in the past two decades (e.g.Angelopoulos, 1996, and references therein).In the PSBL the flows are nearly field-aligned, whereas in the CPS and near the current sheet (CS) region about 70% of the fast flows are predominantly perpendicular to the magnetic field direction (Baker et al., 1996).Thus, it is believed that in the CPS much of the fast plasma transport is across field lines rather than along them.Though fast flows in the CPS tend to be bursty with durations of no more than a few minutes, it is claimed that they are responsible for most of the tail plasma and magnetic flux transport (Angelopoulos et al., 1992(Angelopoulos et al., , 1994)).There is, however, disagreement in the field regarding this point (see Paterson et al., 1998Paterson et al., , 1999;;Angelopoulos et al., 1999).Several studies do associate some bursty flows with neutral lines and substorms (e.g.Nagai et al., 1998;Baker et al., 1999), but a one-to-one correlation between fast bursty flows and the substorm phase is not clear (Angelopoulos et al., 1992).Nakamura et al. (1991Nakamura et al. ( , 1992Nakamura et al. ( , 1994) ) have studied the characteristics of the plasma flows in the CS region, where the magnetic field is weak.They found that high speed flow events are usually restricted to the region tailward of X = Fig. 1.CANOPUS magnetometer stations on a geographical coordinate system grid (cf.Rostoker et al., 1995).
-25 R E , and are predominantly sunward or tailward.Recently, Lyons et al. (1999Lyons et al. ( , 2000) ) found evidence that the flows within the CPS and those within the PSBL may not be distinct phenomena.They found that when they are observed, significant flows appear to exist throughout the entire height of the tail plasma sheet, and are associated with structured currents in the tail.Thus, bursty flows could be associated with some flapping motion of the magnetotail in the Z-direction.This is different from the standard picture, and means that bursty flows in the current sheet region might be present regardless of the presence or absence of an acceleration mechanism, such as a neutral line.In this paper we explore the anisotropy and fast tailward flows in the CPS that can arise from the dynamics of ion beams from the PSBL.In particular, we compare the ion distributions and flows (velocity moments) in the CPS and equatorial CS.

Observational motivation
To consider the origin of some fast tailward CPS flows, we begin by considering an example of a fast tailward flow in the current sheet region that appears unrelated to any significant magnetic activity.The data comes from the CANO-PUS array of magnetometers (Rostoker et al., 1995) and the GEOTAIL LEP instrument (Mukai et al., 1994), and the interval of interest is 04:00-05:00 UT on 4 February 1998.The CANOPUS magnetometer station locations are indicated in Fig. 1 by solid diamonds.The X-component magnetometer (roughly geomagnetic north-south) traces around the interval of interest are shown in Fig. 2.There is little way of magnetic activity prior to 05:30 UT, and the interval 04:00-05:30 UT is essentially characterized by a steady-state magnetic field, where almost no ground activity is evident.The IMAGE network also shows quiet magnetic fields during this interval.The WIND spacecraft was in the solar wind at about 230 R E from the Earth, and observed a southward IMF av- eraging about −2 nT, but slowly decreasing towards the end of the interval of interest.Shortly after 05:30 UT, there is a small poleward border intensification which is centred north of RANK.However, this is well after our interval of interest between 04:00-05:00 UT.We associate the RANK border intensification with a rapid change in the north-south IMF component, from −6 to 0 nT in less than a minute.The poleward border intensification at about 05:40 UT provides evidence of release of magnetic energy via a small dipolarisation, and indicates that there is a stretched field geometry which recovers after 06:00 UT.
GEOTAIL's orbit is at (X, Y, Z) GSM = (-30.4,2.3,-1.8)R E between 04:00-05:00 UT. Figure 1 shows the GEOTAIL footpoint as a pentagon, mapped to the ionosphere using the Tsyganenko (1987) model.GEOTAIL's footpoint is in the vicinity of the magnetometer stations.Figure 3 shows the (a) B x , (b) B y , and (c) B z components of the GEOTAIL magnetic field in GSM coordinates, and the (d) V x component of the velocity moment (dotted).The perpendicular component of V x is drawn as the solid line.For much of the interval GEOTAIL is crisscrossing the neutral sheet, defined by B x = 0. We draw attention to the interesting bipolar flow between 04:30-04:40 UT.At 04:32 UT a strong earthward flow is observed, which has a perpendicular (to the magnetic field) component of more than 800 km/s.This sort of flow would typically be classified as a bursty bulk flow, and possibly attributed to a neutral line or point; the perpendicular flow is observed close to the neutral sheet, and is associated with a large positive B z .Neither the magnetometers nor meridian scanning photometers (data not shown) of the CANOPUS array detect any of the usual evidence for reconnection, and the interval between 04:30-04:40 UT gives evidence only of a steady state.What is more intriguing than the earthward flow is the fast perpendicular (to the magnetic field) tailward flow between 04:36-04:40 UT.This flow is again associated with the neutral sheet region, but instead of a negative B z , as one might expect for a neutral line associated flow, B z ∼ +5 nT.Thus, the tailward perpendicular flow is opposite to the expected convection direction.A fluids-based intuition is not helpful in explaining this flow.It is precisely this kind of tailward flow which we attempt to explain in this paper.Our hypothesis is that, within a simple model framework, this type of tailward flow observation may be explained as being the natural result of PSBL dynamics in a stretched magnetotail.

Modelling technique
With the exception of disturbed intervals, when the geomagnetic topology undergoes major reconfiguration, it is usually assumed that to the first order the CPS is isotropic.However, highly anisotropic plasmas have been observed in the near-Earth plasma sheet (Frank et al., 1996;Chen et al., 2000).Chen et al. (2000) used the WIND satellite to demonstrate that the plasma distribution in the tail CS is multicomponent especially, but not exclusively, during geomagnetically active periods.Such multicomponent distribution functions, when convolved, can result in large velocity moments.Whereas magnetohydrodynamic (MHD) simulations have been very successful in modelling global magnetospheric structure in response to the solar wind (Raeder et al., 1995;Fedder et al., 1995a, b;Frank et al., 1995;Walker and Ogino, 1996;Wiechen, 1999), they are limited because the MHD description deals only with the moments of the distribution function, and cannot explain the non-Maxwellian features, such as beams and non-gyrotropic distribution functions found in the PS and PSBL.In describing the magnetotail as a convecting fluid, MHD provides a good idea of the global flows, but fails on the intermediate and small scale, and does not include nonadiabatic particle effects.Non-gyrotropic and non-Maxwellian distribution functions in the PSBL and near reconnection regions have been modelled using kinetic simulations (Krauss-Varban and Omidi, 1995;Hoshino et al., 2001).However, even full kinetic (or hybrid) simulations are limited at present by computational power, so that studies consider only small-scale and limited regions of the magnetotail.The noninteracting test-particle approach has been used successfully to model numerous regions of the magnetosphere, for example, the magnetopause (Speiser, 1981), the polar region and plasma sheet (Horwitz, 1986a, b;Horwitz et al., 1989), the current sheet region (Lyons and Speiser, 1982) and neutral lines (Speiser and Martin, 1992;Martin et al., 1994;Smets et al., 1996Smets et al., , 1998)), as well as the near-Earth regions (Delcourt et al., 1990(Delcourt et al., , 1991(Delcourt et al., , 1996)).A great strength of the noninteracting test-particle approach is that it is computationally efficient compared to the Vlasov and MHD methods; thus, even full magnetospheric models can be used with inclusion of nonadiabatic particle effects (Ashour-Abdalla et al., 1991a,b, 1993, 1994, 1995, 2000).Not withstanding its limitations, with non-selfconsistency being the most obvious, the method has been very successful in modelling different magnetospheric regions, where the collisionless nature of the plasma, accompanied by the apparent unimportance of collective effects, results in the method being appropriate.
When large velocity moments are observed in the CS region, they are often interpreted in terms of a neutral line mechanism.However, theoretical considerations indicate that the attribution of fast bursty flows to a neutral line is not inevitable or necessary (Liu, 2001).Furthermore, our observational example in the previous section shows that there are also cases that cannot easily be explained by MHD-like mechanisms, in particular fast tailward CS flows associated with positive B z .Since ion distribution functions can be complex, a kinetic description is necessary to attempt to explain such observations.In this paper we use test-particle simulations to model complex three-dimensional velocity distributions in the magnetotail.We will show that a simple model of the stretched magnetotail can explain some fast tailward flow measurements in the CPS, and how anisotropy in the CPS can arise as a natural consequence of ion beams in the PSBL.

Magnetic field
The magnetotail model used here is characterized by its highly stretched tail configuration.It was used by Wanliss et al. (2000a) to model the substorm growth phase, during which the midnight sector magnetic field becomes more stretched than at other times.We chose to use this simple model because it is easier to attain a stretched configuration than for the most commonly used statistical models.The Tsyganenko (Tsyganenko and Usmanov, 1982;Tsyganenko, 1987Tsyganenko, , 1989Tsyganenko, , 1996) ) models involve a crosstail current based on a statistical treatment of a large number of observations that represent an average over all configurations of the mag- netosphere.There is good reason, therefore, to expect that these models do not adequately fit substorm growth phase and other times when the magnetotail is very stretched and extremely thin current sheets exist.Modifications to these models, determined on an event by event basis, bear out this conviction (Pulkkinen et al., 1991(Pulkkinen et al., , 1992(Pulkkinen et al., , 1994(Pulkkinen et al., , 1998;;Lu et al., 1999).
Our model comprises only three components; viz. a dipole field, an equilibrium tail component (Zwingmann, 1983), and an azimuthally and radially confined weak magnetic field region (WFR).All three modular components are used to build up the final magnetic field model configuration.The WFR creates a region of minimum-B in the near-Earth plasma sheet.Figure 1 of Wanliss et al. (2000a) (not shown here) gives a profile of our model magnetic field in the current sheet, compared to observations, and shows the region of minimum-B.The WFR is located around 11 R E to accord with observations (Kaufmann, 1987;Baker and McPherron, 1990;Iijima et al., 1993;Baker et al., 1993;Sergeev et al., 1993;Nakai et al., 1997).The tail component is the Zwingmann (1983) two-dimensional equilibrium, which is a daughter of the earlier Harris equilibrium (Harris, 1962).For this component, the most important parameters are the lobe field and the crosstail current thickness.The current sheet thickness in the model increases with distance from the Earth.Figure 4 shows the model lobe field compared to satellite observations (Slavin et al., 1985;Nakai et al., 1991).The coordinate system is centred on the Earth, with the X-axis positive, sunward long the Sun-Earth line, the Z-axis antiparallel to the Earth's dipole moment and the Y -axis perpendicular to the other two axes.The reader is referred to the paper by Wanliss et al. (2000a) for further mathematical details of the magnetic field model.The most important variable in the model is the crosstail current sheet thickness, which acts as a measure of stretching in the magnetotail.Fig- ure 5 shows a meridional comparison of dipole field lines (dot-dash) to the model's magnetic field lines for two different current sheet thicknesses (indicated by dashed and solid lines, respectively).Based on the analysis of Wanliss et al. (2000a), who found that at the end of the growth phase the current sheet can remain extremely thin for tens of minutes, we chose to use a model thickness of 0.1 R E , as measured at 10 R E .This value compares well with observations of extremely thin growth phase current sheets (Sergeev et al., 1993;Sanny et al., 1994;Kubyshkina et al., 1999).Thus, we are considering only situations that arise when the magnetic field is stretched more than average.

Electric field
Instead of modelling the electric field as a constant, it is allowed to vary with position, and is directed primarily in the Y -direction.Many previous studies of test-particle motion have tended to use a model featuring a constant electric field, directed in the Y -direction (Onsager et al., 1991;Ashour-Abdalla et al., 1991a,b, 1996;Larson and Kauffman, 1996), although there are some models that have used non-constant electric fields (Cladis, 1986;Sauvaud and Delcourt, 1987;Joyce et al., 1995).We begin with a dawn-dusk directed electric field that varies spatially in proportion to the square root of the magnetic field, corresponding approximately to equipotential field lines, with the magnetic field scaling with the cross-sectional area of a flux tube (cf.Cladis, 1986;Onsager et al., 1996), viz.
where α is a constant, related to the cross polar cap potential, and B is the magnetic field magnitude.

Cross polar cap potential
To select a reasonable value for α requires knowledge of the cross polar cap potential difference.This potential can be mapped along open field lines, and assuming a crosstail magnetosphere width of 40 R E , the appropriate large-scale crosstail electric field can be approximated by (2) Kan (1988) demonstrated that the polar cap potential drop needed to exceed a minimum level, 70 kV in their model, prior to the substorm expansive phase, when the magnetic topology is highly stretched.Thus, is the minimum allowable potential drop in the new model, which corresponds to E y ∼ 0.3 mV/m in the distant tail, where B z ∼ 1 nT, consistent with other studies (Alexeev et al., 1993).The observational study of Weimer et al. (1992) is of special interest, since they studied variations of polar cap potential, measured during different stages of 64 isolated substorms, and grouped according to AE index.Their Figs. 3 and 7 show the polar cap potential measured as a function of time relative to the time from substorm expansive phase onset.For the small substorms (peak AE < 600 nT) the polar cap potential is in the range of 40 to 80 kV during the late growth phase, shortly before expansive phase onset.Intermediate substorms (600 nT < peak AE < 1000 nT) showed that in some cases potentials over 70 kV were measured more than 2 h before onset, and potentials over 110 kV were measured at 1 h before onset.The large substorms (peak AE > 1000 nT) have cross polar cap potentials from 90 kV to over 150 kV during the growth phase (i.e.α = 11.4 to> 18.6(V /s) 1/2 ).These data indicate that a cross polar potential drop of 60-150 kV (α = 8 − 20(V /s) 1/2 ) is reasonable during the latter stages of a substorm growth phase.We will use a value of α = 16(V /s) 1/2 .The recent study of Miyashita et al. (1999) serves to further confirm that a value α=16 (V /s) 1/2 is reasonable for the stretched configurations of the late substorm growth phase.Miyashita et al. (1999) selected 263 substorm events and statistically investigated, amongst other parameters, the time variations of the electric field in the region −7.5 ≥ X ≥ −52.5RE and −15 ≤ Y ≤ 15R E in GSM coordinates.The electric field was calculated using the frozen-in relation, E = −V × B. They found that, at 10 min before substorm expansive phase onset, and for X ≈ −10R E and X ≈ −28R E , the electric field E y values are large, with an average of about 0.7 mV/m.In the intermediate region, around X ≈ −18 R E , the electric field E y is about 0.4 mV/m.
Since we are also interested in the particle dynamics and distribution functions in the CS region, it is important that the electric field in this region is curl free.To ensure this, E x and E z components are included, given by where it is assumed that E y is slowly varying in the Ydirection.With the full electric field expressed by the terms thus described, one finds that in general E x , E z E y .The validity of these assumptions are considered in Appendix A.

Distribution function mapping technique
The mapping technique we use to calculate the ion distribution functions utilizes Liouville's theorem.To calculate the The upper distribution function shows contours of constant phase space density.Numerous pseudo-particles are followed from unique initial locations in velocity space (lower distribution) to the source (upper distribution).Since the velocity distribution is known at the source, Liouville's theorem allows each initial pseudoparticle to be tagged with a known phase space density.distribution function at a certain location r in the magnetosphere, a three-dimensional array of initial velocities is specified in directions parallel to the magnetic field, V B , in the electric field direction, V E , and in the E × B or convection direction, V C .A velocity grid is chosen which divides each velocity direction into 31 bins, from −2500 to +2500 km/s.Every pseudo-particle has a unique velocity condition with coordinates (r, v i ), where i is one of the 31 3 initial velocity components chosen.Each particle represents a system of the phase space ensemble, which moves along an orbit in phase space, carrying a piece of probability along with it.Each phase particle is used to derive a reduced distribution function f (r, v i ) that depends only on the phase space coordinates of one pseudo-particle (Wanliss et al., 2000b).
Thus, the calculation begins with a single ion (protons are used) placed at its initial position in phase space (r, v i ).This ith ion will have a value f i assigned to its particular velocity, v i .However, before f i can be known, the trajectory must be integrated backward in time until the ion reaches an assumed source region.The integration is achieved via a fourth-order adaptive Runge-Kutta scheme that solves the full Lorentz equation for protons in the given electric and magnetic fields.Once the integration is complete the velocity v s , and position r s , at the source is known.Thus, the phase space density at the source, f s (r s , v s ), may be computed according to some assumed dependence on (r s , v s ).For example, the source distribution may have a Maxwellian form.Liouville's theorem implies that the phase space density attached to r and velocity v i is Otherwise, if the ion never reaches the source region, the phase space density for that initial condition is assigned the value zero: Zeros in phase space correspond to phase space holes in the velocity distribution functions.Figure 6a (left-hand side) shows how a particle, with an initial position and velocity of (r, v i ), moves through space to the source boundary, where it now has position and velocity (r s , v s ).After this integration, the next particle of the ensemble of 31 3 particles is placed at the same starting position r, with a different initial velocity, v i+1 , and integrated to the model boundary, and so on, until all the particles in the initial ensemble have reached the boundary.Even though the starting position r is the same for all 31 3 particles, each has a unique phase space density assigned to it because each v i is unique.Figure 6b (righthand side) shows how two pseudo-particles move backwards in time from their known position in velocity space to a new position in velocity space, but now at the source boundary.
Once the particle reaches the source, where the distribution function is known, Liouville's theorem is invoked to yield the phase space density at the initial location in velocity space.
After the phase space density has been computed for all the initial conditions, i, the three-dimensional velocity space at position r is completely filled.An array of phase space densities in velocity space exists at the starting position r, and the full three-dimensional velocity distribution function at the position r can be plotted.Their source distributions were assumed to be Maxwellian.
We propose to use a boundary condition that does not assume that the distributions in the CS are Maxwellian.Our boundary condition is dependent on ion beams in the PSBL, which appear to be ubiquitous there (Lui et al., 1983;Eastman et al., 1984).Therefore, instead of specifying the boundary condition in the model CS, and then finding the distributions in the PSBL, what is proposed here is to specify the boundary condition at higher latitudes, including the PSBL, and then to model the ion distributions in the CPS, especially around the CS region.
We constrain the model boundary condition by satellite data and previous theoretical studies.The location of the model boundary is specified in space as the last closed magnetic field line in the noon-midnight meridian that crosses the CS at X=-100 R E , which is shown as the solid line in the schematic of Fig. 7. Eastman et al. (1998) used the GEO-TAIL data set to demonstrate that the probability of observing the plasma mantle increases with distance downtail, from 18% near the Earth, to about 50% near 205 R E .On the other hand, the number of intervals during which the plasma sheet is detected, including the CPS and PSBL, decrease with increasing tailward distance.Pure lobe occurrences also decrease in observation probability with distance downtail.
To take these observations into account, the distribution functions along the model boundary are PSBL-and lobelike near the Earth, but become increasingly mantle-like with downtail distance.This is schematically represented in Fig. 7 by the model distribution functions along the boundary.The ion distributions along the model boundary comprise two components: a beam component featuring motion towards the Earth and parallel to the magnetic field, and a core component that has a small field-aligned drift.In Fig. 8 the distributions, shown at different X-locations along the boundary, are shown as isocontour plots of phase space density in the V − V ⊥ plane, where the V direction is towards the earth and parallel to the magnetic field, and V ⊥ is perpendicular to the magnetic field.The downtail distance along the boundary is shown in each of the plots.
In a model of the PSBL, Ashour-Abdalla et al. (1993) demonstrated that the PSBL ion density decreases with distance downtail.Therefore, in our model, as the distance from the Earth increases, the density of the beam term decreases.This is also consistent with the fact that the probability of observing the PSBL decreases with downtail distance (Eastman et al., 1998), assuming that the beam density drops below the detector threshold.Ashour-Abdalla et al. (1991a, 1993) demonstrated that small-scale "beamlets" can develop in the beam component.However, these small-scale features are likely to be quickly smeared by wave-particle interactions, to form a single solid, homogeneous ion beam (Schriver and Ashour-Abdalla, 1991).Accordingly, only a single solid bi-Maxwellian distribution is used for the beam component, with T ⊥ /T || = 5, which is consistent with observations (DeCoster and Frank, 1979;Eastman et al., 1984).The drift speed of the model beam component varies along the boundary, fastest near the earth, and slowest in the distant tail.It is given a maximum of V b ≈ 1200 km/s, at X = −10 R E , which is a reasonable value (Eastman et al., 1984), although faster PSBL beam velocities have been detected.The beam parameters chosen have V bth|| ∼ 0.08V b (Eastman et al., 1984), where V bth|| is the thermal velocity in the magnetic field direction, and V b is the beam speed.
The various parameters for the beam component are shown as the dashed curves in Fig. 9, and the appropriate parameters for the core term are shown as the solid curves.All the variables are allowed to vary linearly between the values chosen at X = −10 and −100 R E .Near the earth the core term is cold, and represents lobe ions (T ⊥ = 25 eV ).The temperature ratio is T ⊥ /T || = 2, which is a reasonable choice (Rosenbauer et al., 1975).Also, the density is very low in the near-earth (n = 0.02 cm −3 ).It increases linearly with distance downtail in order to simulate a transition between pure-lobe to more mantle-like plasma.At X = −100 R E the core term has n = 0.11 cm −3 and T ⊥ = 560 eV.The complete model boundary condition near the earth is typical for the outer edge of the PSBL.Since the variables vary linearly downtail, a region in the midtail exists where the core term of the distributions is intermediate between mantle-like and lobe-like.In the distant tail, the beam term is negligible, and the complete boundary condition resembles a mantle-like ion distribution.
Although the boundary condition does not include sources other than the PSBL-lobe-mantle sources, the use of this approach does not mean that other valid possibilities are not recognised.For example, there is little doubt that part of the plasma sheet population comes from the ionosphere (e.g.O + ).However, it is assumed here that the different sources act independently, and the focus is on the population of the tail that comes from the sources specified.Ionospheric sources are most important in the near-tail, while sources from the plasma mantle are more appropriate further downtail (Ashour-Abdalla et al., 1993).

CPS region
In Fig. 10 typical model PSBL ion distributions are shown at various distances above the neutral sheet (NS), defined as the region in the magnetotail where B x = 0.The velocities are expressed in terms of the direction of the magnetic field (V B ), electric field (V E ), and convection (V C ).Each distribution function is calculated at X = −20 R E in the noon-midnight meridian at different Z locations above the NS.Figures 10a-d correspond to locations at Z = 3.2, 2.25, 0.75, 0.3 R E , respectively.The left column shows the distributions in two-dimensional cuts in the V B − V C plane, and the right column shows one-dimensional cuts along the V B direction.The two-dimensional distributions show contours of constant phase space density, with contours every two decades from −23 to −10 s 3 /m 6 .All these distributions are calculated quite far from the NS, so the magnetic field is almost entirely in the X-direction.Thus, V B corresponds approximately to the earthward direction, and V C corresponds to the −V z direction.
Of course, immediately equatorward of the model boundary the model BC is recovered, viz. a cold core plus an earthward streaming beam with drift velocity of ≈1200 km/s.This is typical for the outer edge of the PSBL. Figure 10a illustrates the distributions deeper in the boundary layer, at Z = 3.2 R E .Here, the earthward beam is still present, but its distribution function has begun to assume a "kidney bean" shape.At Z = 2.25 R E both earthward and tailward stream- ing ions are observed, with the tailward ions at slightly higher field-aligned speeds than the earthward streaming ions.The perpendicular temperature of the tailward beam is slightly smaller than that for the earthward beams.Also, at this location, both beam distributions appear somewhat "kidney bean" shaped, due to their velocity space mapping into a region of increasing magnetic field strength, as described by Onsager et al. (1991).The ion core has a higher phase space density than the beams, and the earthward beam (+V B ) has a higher phase space density than the tailward beam.
At Z = 0.75 R E the earthward beam is still present, but has a lower field-aligned drift velocity, and is colder in the perpendicular direction.At this point the tailward beam has a larger perpendicular temperature than the earthward beam.Its phase space density continues to increase towards lower latitudes.Still closer to the NS, at Z = 0.3 R E , the earthward beam has a lower phase space density than the tailward beam.
In summary, the distributions retain all the characteristics of the PSBL, with a core component, and with the beam components having T ⊥ /T > 1.Initially, only earthward beams are observed; slightly deeper in the model PSBL, a counterstreaming beam is observed.The counterstreaming beam has a consistently larger field-aligned drift velocity than its earthward counterpart.When it is first observed, the counterstreaming beam has a perpendicular temperature that is cooler than that of the earthward beam.However, closer to the NS its perpendicular temperature is larger than that of the earthward beam.Finally, only the tailward beam is clearly observed, as the phase space density of the earthward beam becomes increasingly smaller.Figure 11 shows different two-dimensional cuts of the full velocity distribution functions calculated at Z=3.5, 0.75 R E .The top row shows the cut in the V B −V C plane for the two different Z-positions.The slight drift velocity in the V C direction is evident in this cut.However, the second row indicates that there is almost no drift in the V E direction.The third row shows the V C −V E cuts.The distribution function is isotropic in this plane, and only the slow convective drift is evident.As can be seen from the preceeding figures, the velocity distributions in the CPS can have complicated forms, including asymmetry and counterstreaming.This is contrary to the view that the CPS is a region of warm isotropic ions, but is similar to previous simulation results (Ashour-Abdalla et al., 1993), and consistent with observations of counterstreaming near the NS (Eastman et al., 1984;Nakamura et al., 1991).It demonstrates that PSBL-like distribution functions can occur deep in the CPS.
Several initial positions in velocity space were selected in order to determine the behaviour of particles comprising the core and beam components of the distribution.Figure 12 shows examples of representative proton orbits for different initial positions in velocity space.The upper figure shows the model's two-dimensional ion velocity distribution in the The final position in velocity space of particles comprising part of the earthward moving beam are labelled with crosses and A1, A2.For the core component the labels are circles and B1, B2.For the counterstreaming component the labels are squares and C1, C2, C3.The lower figure shows trajectory traces corresponding to the initial conditions chosen in the upper figure.All of the trajectories converge to a single spatial location, but with different positions in velocity space.The protons from the core component map to the core term in the model BC (see Fig. 13).Therefore, they retain the essential features of the core part of the distribution function at the model boundary.The particles forming the earthward beam come from further downtail and have arrived at the observation location after convection under the model electric field.They map to the earthward beam BC, and have a large field-aligned component of velocity throughout their motion, from source to observation point.Particles comprising the counterstreaming beam also map to the earthward beam part of the source distribution, but mirror in the dipolar magnetic field near the Earth, whereafter they continue to convect towards the midplane until they reach the observation point.

CPS region moment calculations
Figure 13 shows moments calculated from the full threedimensional velocity distribution functions computed at various Z-values at X = −20 R E .The density slowly increases, from the model boundary (at Z = 4.4 R E ) to about Z = 3.5 R E .This increase in density is due to the slow heating of the plasma while there is only a small decrease in the maximum phase space density of the different components.The distribution function of the earthward beam becomes somewhat "kidney bean" shaped, contributing to this effect.At lower latitudes the observed earthward beam has lower phase space density and velocity, and the temperature becomes cooler (cf.Fig. 10).Thus, between about Z = 3.5 − 2.5 R E , where the earthward beam term is decreasing rapidly in phase space density, the overall particle density decreases.At Z ≈ 2.5 R E the counterstreaming beam is easily identifiable in the distribution function.At At this location V B ≈ V z and V C ≈ V x .Several initial positions in velocity space are selected in order to demonstrate the representative behaviour of particles comprising the different components of the distribution function.Particles from the convective phase space density enhancements (CPSDE), are labelled 1-4.The CPSDE at negative/positive V B is formed from mirrored ion beams that initially come from the Northern/Southern Hemisphere.The phase space density of the negative V B CPSDE is slightly higher than its +V B counterstreaming counterpart.The lower figure shows the trajectory traces corresponding to the initial conditions chosen in the upper figure.Arrows show the direction of particle motion.
first the counterstreaming beam is very weak, and contributes little to the overall particle density, but as one moves to lower latitudes the temperature and phase space density of the tailward beam increase, which results in a slight increase in total particle density.The average density through the plasma sheet is about 0.1 cm −3 .
The velocity moments are shown in the lower plot of Fig. 13.The velocities are plotted in terms of the magnetic field (V B ), electric field (V E ), and convection (V C ) directions.At large Z-values V B ≈ 900 km/s and decreases steadily towards lower latitudes.This steady decrease in V B corresponds to the slow decrease of the drift velocity of the earthward beam at lower latitudes.Below Z ≈ 3.2 R E the decrease is faster, due to the appearance of the tailward beam.At Z = 2.0 R E the tailward beam is well developed and results in the overall bulk V B < 300 km/s.Between 0.5-2 R E only a very small magnitude flow is observed in the V B moment, even though there are very fast counterstreaming beams in the distribution functions.Soon the tailward beam dominates the distribution function and results in negative V B .The model predicts that PSBL-like signatures can be observed deep in the CPS.Thus, qualitatively, there is very little difference between the PSBL and CPS ion distributions.The model V C and V E are very constant.The V C component is almost identical to the predicted convection velocity from (E × B)/B 2 (∼ 100 km/s), and the V E component is zero.

CS region
CS region distribution functions were all calculated at Z = 0.0001 R E .Figure 14 shows two-dimensional cuts of the model distribution functions in the V B − V C plane.Each distribution function is calculated in the CS, for the noonmidnight meridian, and at different X-positions, which are indicated in each of the sub-plots.At these CS region locations, V B ≈ V z and V C ≈ V x .Near the Earth, in general, the CS distribution functions comprise a quasi-isotropic core term, as well as phase space density enhancements in the −V C direction (tailwards).Beyond approximately X = −35 R E , the distribution function comprises only a core term.The core term has T /T ⊥ > 1, and the ratio increases with downtail distance.The phase space density enhancements (PSDE) are largest between about 20 and 30 R E downtail.The sequence of distribution functions in Fig. 14 indicates that the presence of the PSDE may perhaps be responsible for tailward flows.
Figure 15 shows a (top) model two-dimensional ion velocity distribution in the V B − V C plane.It is calculated at (X, Y, Z) = (−24, 0, 0.0001) R E .Several initial positions in velocity space are selected, in order to demonstrate the representative behaviour of particles comprising the different components of the distribution function (crosses, circles, squares).Particles from the PSDE regions are labelled 1-4.The PSDE at negative/positive V B is formed from mirrored ion beams that initially come from the Northern/Southern Hemisphere.The phase space density of the negative V B PSDE is slightly higher than its +V B counterstreaming counterpart.The lower figure shows the trajectory traces corresponding to the initial conditions chosen in the upper figure .The magnetic moment for one of these particles, with velocity at the observation location of (V B , V C , V E ) = (500, −1060, 0) km/s is shown, as the dashed-dot line in Fig. 16.The solid line, corresponding to the left ordinate, shows the trajectory trace in the X − Z plane.The magnetic moment is well conserved until the particle reaches Z ≈ −0.02 R E .It becomes unmagnetised at this point and the conservation of the magnetic moment is subsequently vi-Fig.16.The trajectory of a particle in the X − Z plane is shown (left ordinate, solid line) as a reference against which to compare the variation of the magnetic moment (µ) which is the dashed-dot line and right ordinate.The magnetic moment is normalised.The particle comes from the model boundary in the Southern Hemisphere, and arrives at the observation location at (X, Z) = (−24, 0.0001) R E with a velocity (V B , V C , V E ) = (500, −1060, 0) km/s.olated as the particle moves in the thin current sheet to the observation location.
Figure 17 shows two-dimensional cuts for the distribution calculated in the CS region at X = −22 R E .The top row shows cuts in the V B − V C plane at V E = 0 km/s (left) and 200 km/s (right).The middle row shows cuts in the V B − V E plane at V C = 0 km/s (left) and 200 km/s (right), and the bottom row shows cuts in the V C − V E plane at V B = 0 km/s (left) and 200 km/s (right).The middle row, with cuts in the V B − V E plane, shows that the PSDE components extend to the positive V E plane.Thus, the ions comprising the PSDE begin to gain significant momentum in the V E direction once demagnetisation has occurred in the CS region.

CS region moment calculations
The density and velocity moments are calculated for the CS, and shown in Fig. 18.The upper plot shows the density plotted as a function of the downtail position in the CS region.In the distant tail the density is slightly lower than 0.1 cm −3 , and decreases slowly until about X = −30 R E .The decrease in density towards the Earth occurs in the region where only the core term exists (cf.Fig. 14).As was shown in the previous section, the core term at the observation location maps to the core component of the distribution function at the model boundary.Since the model boundary has the core term becoming warmer and more mantle-like with distance from the Earth, it follows that the density would slowly decrease closer to the Earth.At X = −30 R E the PSDE components, which come from the mirrored beams, are well developed and clearly seen in the distribution functions (Fig. 14).It is these components that contribute to the increasing ion density in the near-Earth region.
The velocity moments are even more interesting.The component along the magnetic field direction, V B , is typically very small; in the distant tail it is slightly negative, in the midtail it becomes slightly positive, and in the near tail, closer than about 30 R E , is becomes quite variable but with a negative average.This is due to the PSDE terms; the phase space density of the negative V B term is slightly higher than its +V B counterpart, resulting in an overall neg- ative V B drift.The V E (≈ V y ) term is positive for all distances from the Earth.This is due in part to particle trapping in the crosstail current sheet, and the subsequent meandering motion in the electric field direction.Strong duskward flows occur between 20-30 R E due to the gyration of the meandering ions from the beam components originating in the model boundary.Ashour-Abdalla et al. (1995) reported similar qualitative behaviour in the CS.Of greatest interest is the V C term (earthward/tailward flows), which is fairly constant downtail of about 30 R E .However, once the PSDE components appear in the distribution functions, they contribute to negative V C , viz.tailward flows.Similar non-MHD tailward flows were reported in the modelling work by Ashour-Abdalla et al. (1995).In our model, at about X = −22 R E , the velocity is V C ∼ −570 km/s.Still closer to the Earth the phase space density of the PSDE becomes lower and V C begins to increase to smaller values.
Figure 19 shows the earthward velocity moment plotted for three different current sheet half-thicknesses.Thicker current sheets correspond to more dipolar topologies, and thinner sheets to a stretched configuration (Wanliss et al., 2000a).L z = 1 corresponds to the magnetic field model used thus far.Model current sheet thickness at X = −10 R E equals 0.1 × L z (R E ).Tailward flows are still present for a rather thick current sheet (L z = 11), but are absent by the time L z = 19.For the L z = 11 case the region over which tailward flows occur is narrower, and the minimum is also smaller than for the L z = 1 case.This means that although tailward flows are not as pronounced for thicker current sheets, one could expect to observe such flows during much of the course of the substorm growth phase, during which period the current sheet thins and the magnetotail topology stretches (Wanliss et al., 2000a).

Summary and conclusions
To the first order, ion distributions in the CPS are often considered to be isotropic during quiet times.Strong anisotropy and fast cross-magnetic field flows are usually attributed to very dynamic processes, for instance, during the substorm expansive phase (Baker et al., 1999).However, some theoretical considerations indicate that the attribution of fast flows to processes associated with expansive phase onset is neither inevitable or necessary (Liu, 2001).
Recent studies have demonstrated that the ion plasma distribution in the CPS and CS regions are multicomponent (Frank et al., 1996;Chen et al., 2000), and not only during geomagnetically dynamic periods.When convolved, such anisotropic distributions can appear as fast flows with large velocity moments.Ion flows in the PSBL are nearly fieldaligned, but in the CPS and near the CS region about 70% of the fast flows are predominantly perpendicular to the magnetic field direction (Baker et al., 1996).Thus, it is believed that in the CPS much of the fast plasma transport is across field lines rather than along them.Though fast flows in the CPS tend to be bursty with durations of no more than a few minutes, it is claimed that they are responsible for most of the tail plasma and magnetic flux transport (Angelopoulos et al., 1992(Angelopoulos et al., , 1994)).
Our study has explored some of the anisotropy related to fast flows in the PSBL and tailward flows in the CPS; we have studied the causes of fast earthward CPS flows.In particular, we investigated the possibility that some tailward cross magnetic field flows are the result of normal PSBL dynamics, and that they should not automatically be associated with neutral lines and expansive phase onsets.The procedure followed was based on the numerical modelling of ion distributions in a stretched magnetotail geometry.This stretched geometry was designed to simulate the magnetic topology of the near-Earth during the latter stages of the substorm growth phase.In essence, the investigation considered the consequences of ion beams in the PSBL in terms of possible flows at lower latitudes.We found that mirroring of initially earthward beams from the PSBL, and their subsequent convection to the CS region, resulted in strong anisotropy throughout the height of the CPS; the distributions contained complicated forms, including asymmetries and counterstreaming.At higher latitudes the velocity moments yielded flows that were field-aligned and earthward.Deeper in the CPS velocity moments of the distribution functions yielded fast flows in the anti-earthward direction.There was no clear distinction between the PSBL and CPS, since velocity distribution functions with large streaming components occured throughout the model CPS.This is consistent with observational studies (Lyons et al., 1999(Lyons et al., , 2000)).In addition, the distributions were primarily gyrotropic throughout the CPS.Ion distributions computed in the CS region were anisotropic and no longer gyrotropic, thus providing a clear distinction between the CS and CPS regions.In the CS region the velocity moments were primarily anti-earthward, and these flows were across the magnetic field.Tailward flows (> 400 km/s) were observed in the CS region at X = −20 to −30 R E .
We found that in order to generate fast tailward flows in the CS region, a sufficiently stretched magnetic topology was required.Previous work (Wanliss et al., 2000a) indicates that such topologies are often observed during the latter stages of the substorm growth phase.Thus, our model predicts that fast tailward CS region flows are likely to be observed during the late substorm growth phase, or at other times when the magnetotail is highly stretched.This suggests that fast tailward plasma flows can be obtained without necessarily appealing to magnetotail processes which are associated with dynamic geomagnetic activity such as substorm expansive phases.Although we have not considered flapping in our model, we note that any flapping of the magnetotail would complicate matters futher and cause these flows to appear bursty.At expansive phase onset dynamic magnetotail behaviour occurs and rapid flapping or thinning and thickening of the CPS can occur.If pre-existing tailward flows, generated through PSBL dynamics, encountered such variable fields, then moment calculations would be more complex to interpret, since the magnetic field observed by an instrument might change rapidly enough to cause the distinction between directions parallel and perpendicular to the magnetic field to be unclear in practice.Thus, some measured cross field flows could be associated with the flapping motion of the magnetotail, as suggested by Lyons et al. (1999Lyons et al. ( , 2000)).It is likely that these difficulties cannot be resolved through an examination of moments alone, but require a study of the velocity distribution functions as well, as our results indicate.
Despite the limitations of the test-particle model -selfconsistency being the major one -it has the advantage that it Fig. 2. X-component traces from the Churchill line of CANOPUS magnetometers for 4 February 1998.

Fig. 3 .
Fig. 3. GEOTAIL spacecraft 4 February 1998 observations of (a) B x , (b) B y , (c) B z magnetic field traces in GSM coordinates.(d) The V x component of velocity (dotted) and the V x⊥ cross magnetic field component (solid).

Fig. 5 .
Fig. 5. Schematic showing the dipole magnetic field lines (dotdashed) compared to our stretched magnetotail for two different current sheet thicknesses.

Fig. 6 .
Fig. 6.(a) The diagram illustrates how pseudo-particles, with initial position and velocity (r, v i ) move through space to a source region.(b)The upper distribution function shows contours of constant phase space density.Numerous pseudo-particles are followed from unique initial locations in velocity space (lower distribution) to the source (upper distribution).Since the velocity distribution is known at the source, Liouville's theorem allows each initial pseudoparticle to be tagged with a known phase space density.

Fig. 7 .
Fig. 7.A schematic diagram showing the model boundary, which is the magnetic field line that crosses the NS at X = −100 R E (solid line).The probability of observing mantle/PSBL plasma increases/decreases with distance downtail.Model velocity distribution functions at various positions along the boundary are shown as line drawings with the ordinate corresponding to phase space density and the abscissa to the earthward velocity, V x .
Onsager et al. (1991) demonstrated that a two-dimensional magnetic field model, based on quasi-steady reconnection in the distant magnetotail, is able to reproduce all the observed features in the observed electron and ion distribution functions in the PSBL.Their technique of mapping distribution functions is essentially the same as is used in this work, although the magnetic field and electric field models are different, and we use the full Lorentz equation of motion.They used a plasma sheet, as well as a lobe source boundary condition, to calculate model distribution functions in the PSBL.

Fig. 8 .
Fig. 8. Two-dimensional cuts in the V − V ⊥ plane of the model boundary condition (BC) distribution functions at different downtail positions for positive Z.The V direction is earthwards.In the near-Earth region the distributions are representative of the outer edge of the PSBL.In the midtail region the distributions are representative of an intermediate region including lobe-mantle plasma.The distant tail is representative of an almost pure mantle BC.The contours of constant phase space density are plotted every two decades from −23 to −10 s 3 /m 6 .

Fig. 9 .
Fig. 9. Model variables against downtail position for (top) the perpendicular temperature (T ⊥ ), (middle) number density (n), and (bottom) velocity (V b ) are shown for the core term (solid curves) and the beam term (dashed curves).

Fig. 10 .
Fig. 10.Model ion velocity distributions shown as (left column) two-dimensional distributions in the V B -V C plane, and (right column) cuts in the V B direction (V C =0) at four locations in the noonmidnight meridian, (a) Z=3.2,(b) 2.25, (c) 0.75, and (d) 0.3 R E , all at X=-20 R E .

Fig. 11 .
Fig. 11.Two dimensional cuts through model distribution functions calculated at X=-20 R E and (left column) Z=3.5 R E and (right column) Z=0.75 R E .The velocities are expressed in terms of the direction of the magnetic field (V B ), electric field (V E ), and convection (V C ).The top row shows the cut in the V B − V C plane, the middle row shows the cut in the V B −V E plane, and the bottom row is the cut in the V C − V E plane.

Fig. 12 .
Fig. 12. (Top) Model two-dimensional ion velocity distribution in the V B − V C plane, calculated at (X, Y, Z) = (−20, 0, 1) R E .Several initial positions in velocity space are selected in order to demonstrate the representative behaviour of particles comprising the core and beam components of the distribution.The particles comprising part of the earthward moving beam are labelled with crosses and points A1, A2.The sample particle orbits from the core component are labelled as circles and B1, B2.For the counterstreaming component the labels are squares and C1, C2, C3.The lower figure shows the trajectory traces corresponding to the initial conditions chosen in the upper figure.Arrows show the direction of particle motion.

Fig. 13 .
Fig. 13.Model (top) number density and (bottom) velocity moments calculated at various Z-distances from the NS at X=-20 R E .The velocities are expressed in terms of the directions of the magnetic field (V B ), electric field (V E ) and convection (V C ).

Fig. 14 .
Fig. 14.Model ion velocity distributions shown as two-dimensional distributions in the V B -V C plane.The distribution functions are plotted for different downtail positions, from X = −20 R E to −75 R E .The X-positions are indicated in each subplot.

Fig. 15 .
Fig. 15.(Top) Model two-dimensional ion velocity distribution in the V B −V C plane, calculated at (X, Y, Z) = (−24, 0, 0.0001) R E .At this location V B ≈ V z and V C ≈ V x .Several initial positions in velocity space are selected in order to demonstrate the representative behaviour of particles comprising the different components of the distribution function.Particles from the convective phase space density enhancements (CPSDE), are labelled 1-4.The CPSDE at negative/positive V B is formed from mirrored ion beams that initially come from the Northern/Southern Hemisphere.The phase space density of the negative V B CPSDE is slightly higher than its

Fig. 17 .
Fig. 17. cuts of the model velocity distribution calculated at (X, Z) = (−22, 0.0001) R E .The velocities are expressed in terms of the directions of the magnetic field (V B ), electric field (V E ) and convection velocity(V C ).The top row shows cuts in the V B − V C plane at V E = 0 km/s (left) and 200 km/s (right).The middle row shows cuts in the V B − V E plane at V C = 0 km/s (left) and 200 km/s (right).The bottom row shows cuts in the V C − V E plane at V B = 0 km/s (left) and 200 km/s (right).At this locationV C ≈ V x , V B ≈ V z , V E ≈ V y .

Fig. 18 .
Fig. 18.Model (top) density and (bottom) velocity moments calculated in the magnetotail CS region (Z = 0.0001 R E ) in the noon-midnight meridian.The velocities are expressed in terms of the directions of the magnetic field (V B ), electric field (V E ) and convection velocity (V C ).At these locations in the CS regionV C ≈ V x , V B ≈ V z , V E ≈ V y .

Fig. 19 .
Fig. 19.Model velocity moments in the V C (≈ V x ) direction calculated in the magnetotail CS region (Z = 0.0001 R E ).The moments are shown for three different current sheet half-thickness parameters (L z = 1, 11, 19).L z has been normalised to L z = 1 ≡ 0.05 R E .