A&A 487, 723-729 (2008)
DOI: 10.1051/0004-6361:200809950
D. Verscharen - H.-J. Fahr
Argelander Institute for Astronomy, University of Bonn, Auf dem Hügel 71, 53121 Bonn, Germany
Received 11 April 2008 / Accepted 28 May 2008
Abstract
Context. As a special case of astrophysical MHD shock waves, the solar wind termination shock is typically treated using the MHD jump conditions as they have been determined by Rankine and Hugoniot. A kinetic analysis becomes necessary for both a more detailed view of the governing processes and a deeper understanding of the plasma behaviour.
Aims. In the case of a parallel shock, only an electric field can be considered as the main process decelerating the solar wind ions. This field leads to a strong acceleration for the electrons due to the other sign of their charge and the much lower mass of the electrons than of the ions. This situation enforces a two-stream instability, which is considered to be compensated by wave-particle interactions with electrostatic plasma waves.
Methods. The kinetic approach leads to an equation in Fokker-Planck form, which can be solved by using Ito's calculus for stochastic differential equations.
Results. These two processes (electric field and wave-particle interaction) yield a decelerated subsonic solar wind on the downstream side of the termination shock, showing some new features in the ion distribution function, such as a double-hump structure and a comparatively large number of reflected ions. Within these considerations, we estimate of the spatial size of the shock region.
Key words: plasmas - solar wind - shock waves - magnetohydrodynamics (MHD)
In the inner heliosphere the solar wind represents a supersonic quasineutral plasma flow with typical bulk velocities in the range of about 400 km s^{-1} for the slow solar wind type or 800 km s^{-1} for the fast solar wind. The solar wind density and the kinetic ram pressure of the plasma flow in this region systematically decrease with increasing solar distance. In order to be able to adapt to the local interstellar medium (LISM) with its finite pressure , the solar wind plasma thus has to undergo a shock, the so-called solar wind termination shock, where it transforms into a subsonic downstream plasma flow with downstream Mach numbers (Zank 1999; Fahr 2004; Lee 1997). This termination shock has to arrange the dissipative deceleration of a collisionless plasma flow. The specific structure of astrophysical magnetohydrodynamic (i.e. MHD) shocks is strongly characterised by the orientation of the local magnetic field with respect to the shock surface normal . In this respect ``parallel'' means that the magnetic field upstream of the shock is parallel to the normal of the shock surface, i.e. , and ``perpendicular'' means that the scalar product upstream of the shock vanishes. In the ``parallel case'' there evidently is no magnetic Lorentz force component acting upon the plasma stream that can decelerate the plasma flow. This means that the bulk Lorentz force ( ) cannot contribute to the decrease in the bulk velocity of the plasma in n-direction. Thus, the magnetic field at the first order does not play a dynamic role in the treatment of parallel MHD shock waves. Only by electromagnetic or hydromagnetic waves could it play an indirect role. Thus, at first glance, only one process remains that could enforce the neccessary plasma deceleration, namely the action of the selfconsistent shock-electric field. Ions are slowed down by running against this electric potential wall associated with this field , thereby reducing the bulk momentum of the plasma. If its selfconsistent form slows down the ions, it will also, however, simultaneously lead to an acceleration of the plasma electrons during the shock transit. Of course the electron mass is much lower than the mass of a proton, which is the most abundant ion in the solar wind plasma (i.e. ).
The first result of this electric field influence is a slowed-down ion flow and, in contrast, a strongly over-shooting fast electron flow. This situation creates electric space charges and represents a typical electrostatic two-stream instability. Thus, it cannot represent the final dynamical state of the downstream plasma flow for evident reasons: the plasma would constitute a non-neutral state due to different local densities of electrons and ions. Several different plasma micro-processes are consequently put in operation by this primarily unstable situation which tend to asymptotically bring both particle species to one identical subsonic bulk velocity and to an asymptotic quasineutral plasma condition. In this respect we investigate in detail in the present paper the collisionless wave-particle interaction with consistently excited electrostatic plasma waves (also called Langmuir waves) that are dynamically coupling ions and electrons.
This type of quasiparallel MHD shock envisaged in this paper in fact often appears in astrophysical reality, for example at different regions of the solar wind termination shock where the upstream frozen-in magnetic fields either systematically or temporarily attain small tilt angles with respect to the shock normal (e.g. see Li & Zank 2006; Chalov & Fahr 1996; Zank et al. 2004; Kóta & Jokipii 2006; Schwadron & McComas 2006). In some of the termination shock areas, the periodically changing magnetic field orientation at current sheet layers temporarily leads to parallel shock conditions. In addition, the parallel shock more or less is a permanent feature at high solar latitudes. Furthermore it is worth mentioning that nearly parallel shock conditions also permanently occur at the flanks of the earth's bow shock (Treumann & Scholer 2002). In all these cases typical nonstationarities and strong wave generation have been theoretically predicted as characteristic features of such quasiparallel shocks (Krasnoselskikh et al. 2002) and also observed by in-situ CLUSTER observations (Lobzin et al. 2007) or VOYAGER observations (Burlaga et al. 2006).
It has been clear from the very beginning that the physical details of the above-mentioned dissipative plasma-wave microprocesses are not embedded in the internal physics of the Rankine-Hugoniot conservation laws of magnetohydrodynamic shock theory, since microprocesses are not specified therein. Thus, a kinetic treatment of the shock transition definitely becomes necessary for the sake of a better physical understanding of the internal shock properties. For the situation of a perpendicular shock, analytic solutions of the ion distribution function have more recently been worked out from an adequately formulated Boltzmann-Vlasov equation (Fahr & Siewert 2006); however, up to now no analytic solutions could be provided for the parallel case due to the mathematically much less accessible and treatable form of the Boltzmann-Vlasov equation governing this case (see Siewert & Fahr 2008). In this paper we solve the adequate Boltzmann-Vlasov equation for the parallel case including stochastic, entropy-generating, wave-particle interaction processes. The numerical solution of this equation is obtained using the Ito-stochastic calculus transforming the governing second-order partial differential equation into a set of stochastic linear differential equations. This procedure has already been successfully applied for different space plasma problems in the past (e.g., Chalov et al. 1995; Chalov & Fahr 1998; Dworsky & Fahr 2000). With this solution method, we then obtain the ion distribution function on the downstream side of the parallel shock in dependence on specific values for the introduced shock parameters, such as the shock thickness, the upstream solar wind density, and solar wind bulk velocity. Some remarks on the observability of the phenomena predicted by our theoretical approach are given in the conclusions of this paper.
We assume a one-dimensional shock configuration where both the plasma bulk flow vector
and the upstream local magnetic field
are
parallel to the shock normal .
The one-dimensional steady-state Boltzmann-Vlasov equation for a collisionless medium moving along the normal
space coordinate s can then be given in the form (see Fahr & Siewert 2006)
As described above, the source of deceleration of the solar wind ions is a space-charge-induced electric field . This field must be consistent with the requirements formulated by the Rankine-Hugoniot MHD shock relations.
Thus it is at first necessary to find an expression that can describe an electric field consistent with the MHD shock conditions. As described
before, the magnetic field can be omitted in these considerations. Therefore, one is left with the fairly simple Euler equation in connection
with a force term that is the result of the electric field:
(2) |
(4) |
(5) |
(6) |
The field given in the form of Eq. (7) takes care of arranging the primarily MHD-desired behaviour to decelerate the ions. The conservation of
energy then automatically leads to an expression for the electron bulk velocity
.
Thus, starting from
(8) |
The quasilinear theory of Landau damping is used to obtain the temporal change of the ion distribution function due to the interaction between particles and electrostatic plasma waves. Considering in the frame of the ions the electrons that are in resonance with the excited plasma waves
(k = wavenumber,
= electron plasma frequency,
= ion velocity in electron bulk frame), one obtains as an expression for the change of the ion distribution function the general Fokker-Planck term in the form
(11) |
(12) |
(15) |
(18) |
Now with the terms derived above the adequate Boltzmann-Vlasov equation attains the form
(21) |
(22) |
The associated Langevin equation for the general case is given by
From this,
is the increment of the Wiener process given by
(25) |
(26) |
Looking at Eq. (24), the splitting of the two acting processes becomes apparent. The effect of the decelerating electric field is described by the first term (the drift term) and the stochastic effect of the wave-particle interaction appears within the second term (i.e. the diffusion term).
To continue at this point, the specification of a bulk velocity profile U=U(s) is necessary. In our case here we represent the velocity step over the shock by a ``''-profile (compare with Lee et al. 1986) in the form
(28) |
(29) |
(30) |
Each integration of Eq. (24) corresponds to the phasespace propagation of one selected ion over the shock structure influenced by drift and diffusion (MacKinnon & Craig 1991). Thus, such integrations must be carried out with a statistically relevant sample of ions (in this case for 30 000
ions) to reach statistical significance. Counting the individual ions on the downstream side in their different velocity space intervals thus finally yields the ion distribution function
.
As the most important velocity moments of
,
the resulting downstream bulk velocity
(32) |
(33) |
From arguments of energy conservation the electron temperature can also then indirectly be achieved. For the further iteration of Eq. (24), the downstream bulk velocity in Eq. (27) now is set equal to the newly found value . Hence, one only obtains the asymptotic, ``real'' value of the resulting downstream bulk velocity after several iterations after it results from the action of the processes.
To estimate the size
of the space charge region where the electric field
dominates, an approximative solution
of Poisson's equation in the following form can be used,
(34) |
To form a rather simple estimate, the densities are derived from the particle flux conservation
Then, the above estimation leads to
(36) |
Using standard parameters (good compilations are given by Scherer et al. 2000; or can be found in the review of Zank 1999; and measurments from space probes can be found at Gazis et al. 1994) for the quiet solar periods; i.e., like , , , , the resulting ion distribution function can be calculated and plotted. The tilde sign on top of velocities indicates that these values are normalised by U_{1} for a better numerical handling; i.e., velocities are given in units of the upstream bulk velocity (e.g., ) and length scales in units of , which is set to the estimated value of . The shock extent is discretised into 10 000 steps over the range of . The initial ion velocities are taken from an upstream distribution function taken as a Maxwell-Boltzmann distribution with . The parameter b is set to , whereas is always chosen so that the ratio amounts to 4:1. This assures that the electron bulk velocity profile decreases continuously from its overshoot value to its asymptotic value . Calculations with different values for this ratio do not show any severe dependence on . The distribution function is plotted in Fig. 1.
Figure 1: Normalised ion distribution function on the downstream side of the termination shock. The underlying parameters are described in the text. The sample consists of 30 000 particles. Two humps have been developed in the distribution with one reflected and one continuing beam. The dotted curves are the fit results described in Sect. 3.4. | |
Open with DEXTER |
The real bulk velocity on the downstream side as calculated by Eq. (31) has a value of
,
which leads to a true compression ratio of
.
The ion temperature amounts to
,
the electron temperature to
.
Thus, the downstream sound velocity takes a value of
(38) |
In view of Fig. 1 one recognises that partial ion reflection obviously occurs at the termination shock. The hump on the left-hand side of the downstream ion distribution function is developed due to ion scattering at electrostatic turbulences and represents ions that have attained negative particular velocity with respect to the shock rest frame. The amount of these ions is determined by the parameters and b. For smaller shock size parameters, the fraction of these ions decreases and vice versa. Within the usual ranges the other parameters have only a very weak influence on the plasma stream properties. Even a change in the classical compression ratio does not change the downstream bulk velocity , since the turbulent interaction is so dominant that it superposes the classical velocity profile. This holds for , whereas it is usually assumed to be approximately 4.
The energy densities contained in the electric field and in the wave field can be derived from the bulk velocity behaviour. The electric field energy
leads to the change in the ion velocity and, therefore, the mean field from a global view is given by
(39) |
(40) |
(41) |
At this point the interest in the two most important parameters and b arises. A set of calculations of the downstream distribution function is carried out by us for different parameters and the moments and T_{2} (temperature both for ions and electrons) are calculated. The dependence on is shown in Fig. 2. The real value of b does not change during the iterations, thus must change because of the rescaling by units of .
Figure 2: Dependence of the moments on parameter . The ordinate for units of velocities is on the left, for the temperatures on the right. Velocities are scaled in units of U_{1}, and b has a value of . | |
Open with DEXTER |
In addition to the above-mentioned moments, the sound velocity on the downstream side of the shock is plotted into the diagram in the same scaled units as . Thus, it is possible to see which cases have no physical meaning, i.e. where the downstream bulk velocity is higher than the speed of sound. In these cases the solar wind would not be shocked down to a subsonic flow. This situation is given for and, therefore, an upper limit is obtained from the numerical simulation. This result is in very good agreement with the above approximative value of . The parameter obviously has an effect on both the deceleration and the change in the temperature.
The other important parameter for the calculation is the size of the interaction region, which is determined by b. In Fig. 3 we show the behaviour of the moments of the ion distribution function in dependence on b.
Figure 3: Dependence of the moments on parameter b. The ordinate for units of velocities is on the left, for the temperatures on the right. Velocities are scaled in units of U_{1}, is scaled as . is chosen to be . | |
Open with DEXTER |
Essentially, this parameter only has an influence on the temperature that is identically represented by the width of the distribution function. The downstream bulk velocity does not show a pronounced relationship with the size of the region where wave-particle interaction occurs. Because of the energy conservation, increasing ion temperature at constant bulk velocity means decreasing electron temperature, and the other way around. However, there is one point where isothermal conditions appear. At the ion and the electron temperatures turn out to be equal. A further transfer of thermal energy from electrons to ions during the ongoing of the propagation further downstream should be forbidden by the second law of thermodynamics (anti-entropic behaviour!). On the other side, the same argument as above holds for , namely that the downstream velocity is higher there than the sound velocity. Hence, the possible range for can be restricted to , which means that, in rescaled units, b takes a value between and . This strong restriction shows that the estimated size of is a good choice for the size of the interaction region.
A better chance to confirm the results of the kinetic shock model presented here is perhaps given through the spectral observation of energetic neutral atoms (ENAs). The interstellar boundary explorer (IBEX) satellite will soon (scheduled launch in August 2008!) present the opportunity to detect spectral intensities of ENA fluxes of shock-generated ENAs from a vantage point near the earth (McComas et al. 2005). From these measurements the ion
distribution function at the termination shock can be reconstructed (Gruntman et al. 2001). For this reason it appears to be adviced to offer an analytical fit result for the numerically obtained ion distribution function of this paper here. Therefore, the superposition of two Gaussians of the
form
Table 1: Result of the fit according to Eq. (42).
The small errors indicate the quality of the fit result. Integrating both Gaussians and summing up yields 1.01, which emphasises this result (i.e., the distribution function is normalised to 1).
The ion distribution function shows two beams, both of which are well-fitted by a Gaussian profile. One of these beams represents ions that are reflected at the shock. The number of these reflected ions are determined by the shock size parameters and b. From classical MHD theory, this detailed double-hump behaviour cannot be obtained. The MHD shock treatment only relates several global quantities to one another and does not present the kinetic details. In our model global assumptions are partially taken from MHD. On the other hand, the moments of the numerically calculated distribution function determine the MHD shock characterisation. Thus, our results are in good agreement with the magnetohydrodynamic theory; however, they show several additional details.
The predominant quantities for downstream bulk velocity and temperature are the parameters that characterise the shock extent. The others do not change the plasma properties on the downstream side appreciably, if they are chosen within the typical range. Connecting the differently gained results leads to the solar wind termination shock having a small shock thickness of about in comparison with other considerations. For example, Siewert & Fahr (2007) find a shock thickness of about for similar conditions for the parallel case (using a typical magnetic field of ). They consider Whistler waves as the process that compensates for the two-stream instability. It is, however, obvious from basic considerations in this paper that electrostatic plasma waves couple the electron and ion bulk velocities much more efficiently.
For the typically chosen parameters, the real compression ratio turns out to be r=2.5. A comparison with measured data is hard to do because of the very poor data basis. Different data from the shock passage of VOYAGER-1 leads to a derivation of the compression ratio of r=2.6_{-0.2}^{+0.4} (Stone et al. 2005), which interestingly enough lies in perfect accordance with our value here. However, this result must be handled with caution, because VOYAGER-1 is supposed to have observed in a termination shock region with a quasi-perpendicular shock situation. On the other hand, this space probe is the only instrument to now that has been able to achieve in-situ measurements of this parameter.
Our calculations in addition show a higher temperature for the electrons than for the ions, except for the situation where isothermal conditions can be reached at as described before. For most of the MHD simulations, this is not the typical result, and even multifluid simulations in contrast indicate a lower downstream temperature for the electrons (e.g., Zank 1999).
The true size of the shock layer is the crucial factor in our model; hence, one should study the possibility of a simultaneous kinetic treatment of the involved electrons. This could yield the ability to find a self-consistent plasma description, wherein the shock extent does not need to be a given parameter but a quantity that results from the calculation itself. Due to the fact that - after the again singular VOYAGER-2 passage - no more in-situ measurements can be expected in the very near future, remote sensing by satellites like IBEX represents the best capability for an observational study of our results with a high statistic significance. Within the energy bands of IBEX, the conspicuous double-hump structure should be observable.
Our model shows that the connection of classical MHD with kinetic methods is a good way to obtain detailed results for the behaviour of shocked space plasmas. The very complex closed solution of the full Boltzmann equation becomes unnecessary using this approach.
Acknowledgements
We are grateful for financial support to the DFG within the frame of the DFG-Project Fa 97/31-2. In addition, we express our gratitude to A. Dworsky for his help in applying the Ito-stochastic calculus.