A&A 368, 622-634 (2001)
DOI: 10.1051/0004-6361:20010028

Anomalous cosmic rays and the generation of energetic neutrals in the region beyond the termination shock

A. Czechowski1 - H. Fichtner2 - S. Grzedzielski1,3 - M. Hilchenbach4 - K. C. Hsieh5 - J. R. Jokipii6 -
T. Kausch7 - J. Kota6 - A. Shaw5


1 - Space Research Centre, Polish Academy of Sciences, Bartycka 18A, 00-716 Warsaw, Poland
2 - Institut für Theoretische Physik IV: Weltraum- und Astrophysik, Ruhr-Universität Bochum, 44780 Bochum, Germany
3 - Service d'Aéronomie du CNRS, BP 3, Verrières le Buisson, 91371, France
4 - Max-Planck-Institut für Aeronomie, Max-Planck-Str. 2, 37191 Katlenburg-Lindau, Germany
5 - Physics Department, University of Arizona, Tucson, AZ 85721, USA
6 - Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721, USA
7 - Institut für Astrophysik und Extraterrestrische Forschung, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany

Received 22 June 2000 / Accepted 19 December 2000

Abstract
Some characteristics of the neutral hydrogen flux of 55-80 keV detected by CELIAS/HSTOF on board SOHO (in particular, the intensity peaks during the periods when the instrument was directed towards the anti-apex of the local interstellar medium) suggest that the observed atoms may be created in the distant heliosphere, beyond the termination shock, from the anomalous cosmic ray protons which become neutralized by charge-exchange with the background gas. The theoretical models on which this conclusion was first based were very simple and missed some important features. Here we report about the first study of transport of low energy anomalous cosmic rays beyond the solar-wind termination shock employing a model of the heliosphere treating self-consistently the flows of plasma and neutral atoms and including the pressure of cosmic rays. The energetic neutral-atom (ENA) flux created by charge-exchange is calculated and compared with the observations of CELIAS/HSTOF. The viability of using ENAs as a means of imaging the structure of the termination shock, the distant heliosphere and the very local interstellar medium is discussed. We also include some results which may be interesting in connection to possible observations by INCA/Cassini and HENA/IMAGE.

Key words: solar wind - interstellar medium - cosmic rays


1 Introduction

Although the distant heliosphere (by which we mean what is also called the transition or interface region or the heliosheath) still cannot be observed directly, the models of the interaction between the solar wind and the local interstellar medium (the LISM) (Baranov & Malama 1993; Zank et al. 1996; Linde et al. 1998; Ratkiewicz et al. 1998; Pogorelov & Matsuda 1998; Tanaka & Washimi 1999) support the general picture first proposed by Parker (1963). According to it, the solar wind plasma flow decelerates to subsonic speed at the termination shock and forms a stagnation point towards the LISM apex as well as a long wake (the heliotail) towards the anti-apex. The solar plasma is separated from the interstellar plasma by a boundary surface (the heliopause) which can, however, be crossed by the neutral atoms from the LISM. In the heliosheath, their charge-exchange interaction with the plasma ions, which converts the neutrals into pick-up ions (PUIs) and the plasma ions into neutral solar wind particles, affects significantly the structure of the heliosphere and the near LISM (Baranov & Malama 1993). Recently developed models (Linde et al. 1998) include this interaction as well as the effect of the interstellar magnetic field, which may contribute significantly to the LISM pressure. The time-dependence of the solar plasma flow, in consequence of the solar cycle, has been included in the models by Pogorelov (1995) and Tanaka & Washimi (1999).

In this paper we are concerned with the energetic ion population in the distant heliosphere, beyond the solar wind termination shock, and the possibility of imaging their distribution by means of the energetic neutral atoms (ENAs) into which these ions are converted by charge-exchange with the atoms of the background gas (Fig. 1) As the ion transport reflects the features of plasma flow and the magnetic field structure, this approach could provide information about the structure of the heliosphere at distances still not accessible by other means. The ions are the PUIs carried outwards by the solar wind and (the case we concentrate on) the anomalous cosmic ray (ACR) particles accelerated at the solar wind termination shock (for a review of ACR see Jokipii 1990). The energy range considered is a few keV to 100-200 keV, because the charge-exchange rate and the ion number density are then sufficiently high to produce an ENA flux which may be observed from the Earth's orbit. In our approach we combine calculations of the distribution and modulation of the ACR spectrum based on a model of the heliosphere, that is more realistic than those used in the past, with a discussion of the ENA data. Preliminary results of this approach were presented at the 26th International Cosmic Ray Conference (Czechowski et al. 1999a-c)

  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{h2277f1.ps}\end{figure} Figure 1: A view of the heliosphere and its neighborhood. A hydrogen atom H from the LISM enters the heliosphere, is ionized at A to become a pickup ion, convected to the termination shock, accelerated, transported by combined diffusion and convection into the heliotail, neutralized at B by charge-exchange with another atom from the LISM to become an energetic neutral atom
Open with DEXTER

The first potential evidence of heliospheric ENA was obtained by CELIAS/HSTOF on SOHO (Hilchenbach et al. 1998). The flux of 55-80 keV ${\rm mass} = 1$ particles observed by CELIAS/HSTOF during quiet times in 1996 and 1997, was interpreted by Hilchenbach et al. (1998) as that of hydrogen atoms converted from ACR protons in the outer heliosphere, as anticipated by Hsieh et al. (1992). The peaking of the flux during periods, when the instrument field-of-view included the LISM anti-apex direction, agrees with the anisotropy predicted by Czechowski et al. (1995) and Czechowski & Grzedzielski (1997) on the basis of simulations of the ACR transport in the heliosheath. The data obtained in the following years are subject to much larger uncertainty because of higher ion background when approaching solar maximum. A more detailed discussion is given in Sect. 4.

While some qualitative conclusions from the simulations by Czechowski & Grzedzielski (1995, 1997) appear to have a more general validity, the models used were very simple and lacked some of the features expected to be important. As a next step towards a more realistic description of the ACR transport in the distant heliosphere, we present calculations of the ACR distribution beyond the termination shock and the associated ENA flux based on the gas dynamical 5-fluid model of the heliosphere constructed by Kausch (1998), see also Fahr et al. (2000). The model does not include the magnetic field and (as the simple models used in earlier simulations of the ACR distribution beyond the termination shock) is axisymmetric relative to the LISM apex-antiapex line. On the other hand, compared with the best available results, Kausch's model gives a reasonable approximation of the termination shock and the heliopause structure as well as the density distribution of the neutral background. The model allows us to study the adiabatic acceleration or deceleration of ACR beyond the termination shock, which was impossible in the incompressible flow models used previously. In result, we can describe the modulation of the ACR spectrum in the heliosheath region.

In Sect. 4 we discuss the observations by CELIAS/HSTOF, including some data from later periods. Although the directional dependence of the flux is consistent with theoretical calculations, the predicted ACR ENA intensity is low compared to the most recent estimation of the observed flux (we take into account the results of the recent cross-calibration of HSTOF with ACE and IMP8). We point out the signatures of the ENA/ACR signal which seem to be relatively model-independent. Finally, we present some simulations which can be interesting in view of possible observations by the CASSINI spacecraft.

2 Transport and acceleration of ACR in the outer heliosphere

2.1 Transport equation

The transport equation for the pitch-angle-averaged particle distribution $f({\vec x},p,t)$ is (Parker 1965)

\begin{displaymath}\partial_t f=\nabla \cdot {\vec \kappa} \cdot \nabla f-{\vec ...
...\frac{\partial f}{\partial \ln p}\nabla \cdot {\vec V}-\beta f
\end{displaymath} (1)

where ${\vec V}$ is the velocity of the plasma flow, ${\vec \kappa}$ the particle diffusion tensor, p the absolute value of the particle momentum, and $\beta $ the loss rate. The loss rate includes the losses due to charge-exchange ( $\beta_{\rm cx}$) with low-energy atoms of the background and the losses caused by the energy-changing processes (Coulomb collisions and ionization of the background atoms) which remove particles from the considered energy interval. We use the test particle approximation, so that the feedback effect of low energy ACR on the background flow is disregarded (the 5-fluid model of Kausch 1998, includes the cosmic ray pressure terms, but these are not derived using our solution of Eq. (1)). Although solar cycle effects introduce a natural time dependence, we have so far restricted our calculations to stationary cases.

The plasma flow velocity ${\vec V}$ and the number density of the neutral atoms (which determines the charge-exchange rate) are given by the specific model of the background heliospheric flow (Kausch's gas-dynamical model or Parker's incompressible flow). We assume axial symmetry in space (with respect to the LISM apex-antiapex line) which makes the problem 3-dimensional (the third coordinate is energy).

The components of the diffusion tensor ${\vec \kappa}$ should be chosen to be consistent with the structure of the heliospheric and interstellar magnetic field. The background flow models we use do not include the magnetic field, so that a separate argument is needed (see the following subsection).

The energy dependence of the parameters appearing in the transport equation is illustrated in Fig. 2: the loss rate term for ACR protons and the diffusion coefficient assumed inside the heliopause (see below). The cross sections are obtained using the interpolations given in Barnett et al. (1990). The charge-exchange rate on hydrogen $\beta_{{\rm cx,H}}$ is shown for the sample value of neutral hydrogen density $n_{\rm H}=0.03~{\rm cm}^{-3}$, suggested by the Kausch model (see Fig. 5 below) as a typical value between the shock and the heliopause (in the simulation $\beta_{{\rm cx,H}}$ is weighted by the local number density $n_{\rm H}$). This is lower than the LISM $n_{\rm H}$ (taken to be 0.1  ${\rm cm}^{-3}$) due to filtering caused by the interactions with plasma flow. The neutral helium density, needed for $\beta_{\rm cx,He}$, is not given by the model. In view of helium filtration being less effective than for hydrogen (see e.g. Fahr 1990), we assume that its number density between the shock and the heliopause is the same as in the LISM ( $n_{\rm He}=0.01~{\rm cm}^{-3}$).


  \begin{figure}
\par\includegraphics[width=8.6cm,clip]{h2277f2.ps}
\end{figure} Figure 2: The loss rate $\beta $ and the diffusion coefficient $\kappa _1$ (inside the heliopause) for ACR protons plotted against particle energy. Also shown are the contributions to $\beta $ from ACR protons charge-exchange on neutral hydrogen and neutral helium (dotted lines) and their sum (dashed line) assuming $n_{\rm H}=0.03$  ${\rm cm}^{-3}$ and $n_{\rm He}=0.01$  ${\rm cm}^{-3}$. Low value of $n_{\rm H}$ is due to filtering. For most of the energy interval, $\beta $ is well approximated by the charge-exchange terms.The remaining contribution is due to inelastic scattering, which shifts particles into a different energy range
Open with DEXTER

2.2 The diffusion tensor and the magnetic field

The ACR diffusion tensor may be expressed through three components: $\kappa_{\parallel}$, $\kappa_{\perp}$ and $\kappa_{\rm D}$ with two first components corresponding to parallel and perpendicular diffusion respective to magnetic field and the last term describing the drift (alternatively, one may introduce the drift term ${\vec V_D}\cdot\nabla f$ into the transport equation). The magnitudes of the parallel and perpendicular components are usually estimated from the observed ACR spectra using modulation models (e.g. le Roux et al. 1996) or from models of the turbulence in the solar wind (e.g. Zank et al. 1998).

Downstream of the termination shock the structure of the magnetic field must strongly deviate from the upstream region (where the Parker's spiral is a useful first approximation). Convection by the downstream plasma flow with a stagnation point structure makes the field essentially three-dimensional and enhanced in magnitude close to the surface of the heliopause (Nerney et al. 1991, 1993). Even more important is the effect of the solar cycle: changes in shape of the neutral sheet and field reversals with a quasi-period of about 11 years. As the plasma flow beyond the shock is slow, the effects on the magnetic field will pile up (Nerney et al. 1995). At any time the field will vary strongly in space with a typical scale of 100 AU (the heliotail region) or less in the forward (upwind) region. The expected characteristic scales of the ACR distribution are of similar order.

In result of this, contrary to the case of the inner heliosphere, one must deal with a very complicated magnetic field structure. Most of the available numerical simulations of the heliospheric flow do not include the field inside the heliopause. (It is taken into account, together with the time-dependence, in Tanaka & Washimi 1999, but their model has no interaction with neutrals. In Linde et al. 1998, the time-dependence of the inner field is not included.) As we have to restrict our model to the axisymmetric case, we are forced to make a very crude approximation. We assume that the field is disordered and make the scalar (isotropic) diffusion approximation (in the following we show that at low energy, where the loss terms are large, diffusion effects are not dominant). The drift term can be shown to be relatively unimportant in the low energy region and is not included. The diffusion coefficient between the termination shock and the heliopause we approximate by $\kappa_1=(1/3)\kappa_{\parallel}$ (corresponding to an isotropically disordered field) with $\kappa_{\parallel}$ given by the phenomenological formula of le Roux et al. (1996)

\begin{displaymath}\kappa_\parallel=\kappa_0 (v/c) F(P)\ (B_{\rm e}/B)
\end{displaymath} (2)

where $\kappa_0=3.75\ 10^{22}~{\rm cm}^2/{\rm s}$, v/c is the particle speed divided by speed of light, P is the rigidity in GV, F(P)=P for P>0.4 and F(P)=0.4 for P<0.4, B is the magnetic field and $B_{\rm e}$ the magnetic field at 1 AU. The values of $\kappa _1$ used in our simulations are shown in Fig. 2.

Beyond the heliopause the field must approach the local interstellar field, although it would be modified close to the heliopause. Neither the direction nor the magnitude of the local interstellar field is known. In addition, one cannot expect that the formulae obtained for the solar system would work for the diffusion tensor components in the LISM. Estimations used for galactic transport of cosmic rays (see Axford 1981) are derived for very high energy. When extrapolated to energies of interest to us the resulting values are very high. We can use this as an indication that the diffusion coefficients outside the heliopause are probably higher than inside (which may seem natural in view of a much larger spatial scale and possibly lower turbulence levels). Because of axial symmetry we have to keep the scalar diffusion approximation also in the LISM region. As the ACR particle density in this region is low, the error due to this assumption should not be very important for most of the results. In most of the calculations we take the diffusion coefficient outside the heliosphere to be $\kappa _2=10^2\kappa _1$ (a larger value is likely, but a change would only affect the regions where the ACR distribution is very small: see Fig. 8 in the next section).

2.3 The background flow

The model of Parker (1963) is defined by a simple analytic expression for the (axisymmetric) flow potential function $\Phi$, which defines the plasma flow velocity by $(\rho)^{1/2}{\vec V}=\nabla \Phi$ ($\rho$ being the plasma density). The flow is incompressible ( $\nabla\cdot{\vec V}=0$) because the potential is assumed to have $\Delta\Phi=0$ (the exception is the heliopause surface, at which $\rho$ can have a jump). The form of $\Phi$is such that the boundary conditions for the flow speed are satisfied, although only in the limiting case of the radius of the termination shock (taken to be a sphere) being much smaller than the distance to the heliopause. A modification of this model for the case of larger shock radius was obtained by Suess & Nerney (1990). In our application, we assume the following values of the relevant parameters (unless specified otherwise): shock radius = 79 AU, distance Sun-stagnation point = 110 AU, direction-averaged outflow speed from the shock = 100 kms-1, speed of the LISM relative to Sun = 26 kms-1, neutral hydrogen density between the shock and the heliopause = 0.1  ${\rm cm}^{-3}$, the same for helium = 0.01  ${\rm cm}^{-3}$. While Parker's or Suess & Nerney's models give at best a rough approximation to the heliospheric flow, they produce a stagnation point and a heliotail, features of particular importance for the ACR distribution. We use these models for comparison and illustration.
  \begin{figure}
\par\includegraphics[width=7.8cm,clip]{h2277f3.ps}\end{figure} Figure 3: The heliopause and the termination shock in the models of the heliospheric flow by Parker, Suess & Nerney and Kausch. For Parker's and Suess & Nerney's models we take the shock to be a sphere of radius 79 AU centred at the Sun and assume a 110 AU Sun-stagnation point distance
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8cm,clip]{h2277f4.ps}\end{figure} Figure 4: Number density of the neutral hydrogen in Kausch's model outside the termination shock. The density profiles for three directions: tail, side ( $\theta =90\hbox {$^\circ $ }$) and apex are shown as functions of heliocentric distance. Note the hydrogen wall in the apex and side profiles and the non-spherical shock geometry. The LISM value of 0.1  ${\rm cm}^{-3}$ assumed for Kausch's solution is probably too low (Lallement et al. 1991)
Open with DEXTER

Most of the calculations presented here employ a more refined model by Kausch (1998), see also Fahr et al. (2000). This is a solution of gas-dynamical equations describing the interaction of the solar wind with the interstellar medium. Five fluid components are included: the bulk plasma (in one-component description), pick-up protons, neutral hydrogen and the anomalous and galactic cosmic rays, the latter two represented by the pressure terms (determined by higher energy values than those we consider). The flow is two-dimensional with axial symmetry respective to the LISM inflow direction. Magnetic field is not considered. Although the treatment of the neutral component in the model is simplified compared to the kinetic approach by Baranov & Malama (1993) or the multi-component method of Zank et al. (1996), the general structure of the solution is not very different from Zank et al. (1996).

Figures 3 to 6 illustrate some of the relevant features of Kausch's model. In Fig. 3 the shapes of the termination shock and the heliopause are shown, together with those for the simpler models we use for comparison (in Kausch's model there is also a bow shock, not shown). The neutral hydrogen density distribution (which is important for the loss rate term as well as for the creation of ENA) is presented in Fig. 4 (note the hydrogen wall in front of the heliopause, dashed line). Since neutral helium is not considered separately in Kausch's model, we assume a constant value $n_{\rm He}=0.01$ cm-3. The shock compression ratio (which varies over the shock surface), the upstream pick-up proton density and the implied boundary values of the slope and the amplitude of f(p)are displayed in Fig. 5.

Figure 6 presents sample profiles of the divergence of the plasma flow velocity $\nabla\cdot {\vec V}$, which determines the adiabatic energy changes of the ACR. This divergence is not directly calculated by the gas-dynamical code and is less precisely known than other quantities (we have to estimate it from the solution on the output grid). While the flow divergence upstream of the shock is well approximated by 2V/r (so that it is of the order of $5\ 10^{-8}$ ${\rm s}^{-1}$ at the shock), between the shock and the heliopause it can be of either sign and typically of the order of few 10-9 ${\rm s}^{-1}$, comparable to or exceeding the loss rate term. The divergence term is thus important and constitutes an effective source (if $\nabla\cdot{\vec V}<0$) or a sink ( $\nabla\cdot{\vec V}>0$) of particles of given energy and has a part in the modulation of the spectrum.

2.4 ENA from ACR

The calculated ACR distribution can be used to determine the flux intensity of ENA from charge-exchange. It is given by the line integral of the ACR flux $j_{\rm ACR}=p^2f({\vec x},p)$ weighted by the charge-exchange cross section, the neutral density and the extinction factor ${\rm ext}({\vec x},E)$:

$\displaystyle j_{\rm ENA}(\vec{n},E)$ = $\displaystyle \int_{s_0}^{\infty} {\rm d}s
\left(n_{\rm H}({\vec x})
\sigma_{\rm cx,H}(E)\right.$  
    $\displaystyle \left.\,+\,\, n_{\rm He}({\vec x})\sigma_{\rm cx,He}(E)\right)$  
    $\displaystyle \times\,\, j_{\rm ACR}({\vec x},E){\rm ext}({\vec x},E)$ (3)

where ${\vec x}={\vec x}_0+\vec{n}(s-s_0)$ is the line of sight (in the direction given by the unit vector $\vec{n}$) from the observer's position (at ${\vec x}_0$). As our model is time-stationary, the time-delay is not considered. The extinction factor is calculated using the charge-exchange and stripping cross sections from Barnett et al. (1990). In the considered energy range these cross sections are of the order 10-16  ${\rm cm}^2$ or below, so that the extinction factor is close to 1.

The scale of the ACR flux intensity at the shock, which sets the boundary condition for our simulation, is obtained by matching with the shock spectrum of Stone et al. (1996). We use the higher of their values for the flux, corresponding to the year 1987. As the spectrum in our model is non-uniform over the shock, the matching is done for the apex direction, at energy $E_{\max}$.


  \begin{figure}
\par\includegraphics[width=7.1cm,clip]{h2277f5.ps}\end{figure} Figure 5: Boundary conditions at the shock in Kausch's 5-fluid model. The shock compression ratio r, the power q and the amplitude a (dotted line) of the assumed proton shock spectrum together with the upstream pick-up proton number density (relative units) are plotted against the position on the shock surface (labelled by the angle $\theta $ counted from the apex direction). The irregularities are due to the output grid
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.3cm,clip]{h2277f6.ps}\end{figure} Figure 6: Flow divergence outside the shock in Kausch's 5-fluid model plotted as a function of the distance from the Sun for three directions: apex, side ( $\theta =\pi /2$) and tail. Note the non-spherical geometry of the termination shock
Open with DEXTER

3 ACR distribution beyond the shock: Results from simulations

In this and the following section we present the results of our model calculations. The ACR density distributions are obtained by solving numerically the transport Eq. (1). The ENA flux is obtained from Eq. (3).

The models used to specify the parameters are described in Sect. 2. The boundary conditions are discussed in Appendix A.

By Kausch's model (A and B) we refer to the simulations in which the plasma flow and neutral hydrogen distribution are given by Kausch's solution (Figs. 3-6). Similarly, Parker's model denotes the case when the plasma flow is defined by Parker's formula with a particular choice of parameters (see Fig. 3).

The electron and neutral hydrogen densities in the LISM are taken to be $n_{\rm e}=0.1$ cm-3, $n_{\rm H}=0.1$ cm-3. For the neutral helium we assume a constant density $n_{\rm He}=0.01~{\rm cm}^{-3}$. The ACR diffusion coefficient inside the heliopause $\kappa _1$ and the loss rate $\beta $ used in the transport equation are shown in Fig. 2. The diffusion coefficient outside the heliopause $\kappa _2$ is assumed (unless explicitly specified otherwise) to be $\kappa _2=10^2\kappa _1$. The cross sections for charge-exchange are from Barnett et al. (1990). We assume that the 100 keV ACR protons flux at the termination shock (apex point) is equal to $7.0 \ 10^{-3}$  $({\rm cm^2~s~sr~keV})^{-1}$ (Stone et al. 1996). The calculated flux intensities of ACR and ACR ENA are directly proportional to this value.

We use different versions of Kausch's model to illustrate the dependence of the solutions on the model parameters. The most complete case (model A, basic model) corresponds to the boundary conditions at the shock as shown in Fig. 5 (see discussion in the Appendix). The ACR shock spectrum varies over the shock surface and is weighted by the upstream pick-up proton density.

Model B (pure power law) corresponds to ACR spectrum at the shock having a constant slope (flux $\sim $E-1.42) and constant amplitude all over the shock surface. This was also assumed for the simulations using Parker's flow.


  \begin{figure}
\par\includegraphics[width=8.2cm,clip]{h2277f7.ps}
\end{figure} Figure 7: ACR proton distribution for Kausch's model A (solid lines) and Parker's model (dashed) presented as a function of distance from the Sun for different directions: $\theta =0\hbox {$^\circ $ },\ 104\hbox {$^\circ $ },\ 133\hbox {$^\circ $ },\ 162\hbox {$^\circ $ },\ 180\hbox {$^\circ $ }$ ( $\theta =0\hbox {$^\circ $ }$is the LISM apex direction). The lines start at the termination shock (in Parker's model a sphere, in Kausch's a bullet shape). On crossing the heliopause the profile slopes change sharply due to a jump in the coefficient of diffusion: ( $\kappa _2=10^2\kappa _1$)
Open with DEXTER

Figure 7 presents, for the cases of Kausch's model A and Parker's model, the spatial distribution of the ACR protons at 63 keV (approximately the average energy of the ENA H observed by CELIAS/HSTOF on SOHO) in terms of constant-direction profiles: the density is shown as a function of the heliocentric distance for different directions labelled by the angle $\theta $ ( $\theta =0\hbox {$^\circ $ }$ is the LISM apex and $\theta=180\hbox{$^\circ$ }$ the anti-apex direction). The most visible difference is a much slower fall-off of the ACR number density down the heliotail in the basic model, caused in part by a lower loss rate (in our calculations based on Parker's model the neutral hydrogen density is set to $0.1~{\rm cm}^{-3}$: compare Fig. 4). The behaviour of the profiles for different angles near the shock is the result of non-uniform shock spectrum and non-uniform pick-up ion density at the shock in model A.

3.1 The apex-antiapex asymmetry

The main feature of the spatial ACR distribution beyond the termination shock is the apex-antiapex asymmetry, which is well-pronounced also in the ENA flux. The ACR particles concentrate in the region of the heliotail (LISM anti-apex) and consequently the ENA flux is highest from this direction. This is a quite general conclusion, not dependent on the details of the model. There are two main mechanisms contributing to this asymmetry:

1.
Convection of the ACR ions by the plasma flow with a stagnation point structure. The ions are convected from the shock surface ultimately towards the heliotail. In the apex direction, beyond the stagnation point, the outward ACR diffusion competes against the local plasma velocity so that the stationary ACR distribution will decay quickly with distance. In the heliotail direction the plasma flows outward so that convection and diffusion operate in the same direction;
2.
The second mechanism applies only if the diffusion coefficient is much larger outside than inside the heliopause ( $\kappa _2 \gg \kappa _1$), as assumed in most of our calculations. The heliosphere is then a region of relatively slow diffusion and the heliopause is, approximately, a free escape boundary (at which the ACR number density must fall to a low value). The asymmetry appears because the distance from the shock to the heliopause increases toward the heliotail. For high energy ions, with the diffusion coefficient sufficiently high to make the convection term negligible, this mechanism will still cause an asymmetry in the spatial distribution as long as the diffusion length is not too large to make the diffusion approach invalid.


  \begin{figure}
\par\includegraphics[width=8.1cm,clip]{h2277f8.ps}\end{figure} Figure 8: ACR proton distribution at 63 keV, Kausch's model B (pure power law spectrum at the shock). Four different assumptions about the outside diffusion coefficient $\kappa _2$: $\kappa _2/\kappa _1=10^3,\ 10^2,\ 10^1,\ 1$. When $\kappa _2/\kappa _1=1$ the heliosphere is not a free-escape surface, but the ACR distribution is still asymmetric due to mechanism 1 (convection, see the text)
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{h2277f9.ps}\end{figure} Figure 9: ENA 63 keV hydrogen flux from ACR charge-exchange corresponding to the ACR proton distributions of Fig. 8. The flux is plotted as a function of direction (angle $\theta =0\hbox {$^\circ $ }$is the LISM apex) for the case of an observer in the inner solar system with the line of sight not passing close to the Sun. The choice of the outer diffusion coefficient $\kappa _2$ does not much affect the ENA flux in this energy range
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=7.8cm,clip]{h2277f10.ps}\end{figure} Figure 10: ACR proton distribution at 1000 MeV, Kausch's model B (pure power law spectrum assumed at the shock) for $\kappa _2/\kappa _1=10^2,\ 1$. The values of $\theta $ are as in Fig. 8. Although at this energy, which is high for ACR, the diffusion coefficient is large ( $\kappa_1=5.55\ 10^{23}~{\rm cm}^2/{\rm s}$) the apex-antiapex asymmetry is pronounced as long as $\kappa _2 \gg \kappa _1$ (mechanism 2, see text)
Open with DEXTER

Figures 8 and 9 show the ACR density and the ENA flux for four different assumptions about the diffusion coefficients $\kappa _2$ and $\kappa _1$ (outside and inside the heliopause). The case $\kappa_2=\kappa_1$ illustrates the first mechanism, with the heliopause not being a free escape surface. Although the ACR densities look different in detail in all three cases, there is still a concentration of the ACR ions in the heliotail and the ENA flux in all those cases has its maximum in the anti-apex direction. For high energy ACR, with convection insignificant relative to diffusion, the asymmetry in the density distribution is still present due to mechanism two, assuming that the diffusion coefficient is much higher outside the heliopause (Fig. 10). The energy of 1000 MeV, high for ACR, is used to illustrate the case of large $\kappa _1$ (Eq. (2)).


  \begin{figure}
\par\includegraphics[width=7.8cm,clip]{h2277f11.ps}\end{figure} Figure 11: ACR distribution at 16 keV: complete solution compared to the no-diffusion approximation (convection + losses). The irregularities of the latter are caused by the flow divergence term. The profiles are shown for directions $\theta =0\hbox {$^\circ $ },\ 104\hbox {$^\circ $ },\ 133\hbox {$^\circ $ },\ 162\hbox {$^\circ $ },\ 180\hbox {$^\circ $ }$
Open with DEXTER

3.2 ACR transport in the heliosheath: Different contributions

We can use the simulations to estimate the relative importance of different terms in the transport equation Eq. (1). At lowest energies (up to $\sim $25 keV/nucleon in Kausch's model) the no-diffusion approximation is reasonably good (Fig. 11). In this case the ACR distribution can be estimated by integrating the ordinary differential equation $V{\rm d}f/{\rm d}s=-\beta f$describing convection and losses (the adiabatic acceleration term can be included if one assumes a fixed power law spectrum) along the plasma flow line. At higher energy (63 keV) this approximation is no longer good and diffusion must be included.

The effect of the adiabatic acceleration term is presented in Fig. 12. There is a region in near heliotail where adiabatic compression (due to $\nabla\cdot{\vec V}<0$) acts as an effective source of the ACRs and so increases the ion density in the anti-apex direction. This was already considered in an earlier work (Czechowski et al. 1995), in which an ad hoc modification of Parker's flow was used. The effect is negligible in the forward region, where large gradients of the ACR number density make the diffusion term more important (although the absolute magnitude of the divergence term can be larger there than in the heliotail region).


  \begin{figure}
\par\includegraphics[width=8cm,clip]{h2277f12.ps}\end{figure} Figure 12: ACR protons distribution in Kausch's model A for 10, 100 and 1000 keV (only apex and anti-apex directions shown). Dotted lines correspond to the case with the adiabatic acceleration ($\sim $ $\nabla\cdot {\vec V}$) term set to zero
Open with DEXTER

The ACR protons are expected to be derived mainly from pick-up protons. Models suggest that the density of the pick-up protons upstream of the shock should be larger in apex than in anti-apex direction (opposite to our results for the ACR). If one compares the ACR distribution for model A with the one obtained by assuming a uniform pick-up ion density over the shock (equal to the maximum value) one finds that the ACR density in the heliotail is in the first case diminished by a factor approximately equal to the pick-up proton apex-antiapex density ratio ($\sim $3 in Kausch's model A, see Fig. 5), but remains much larger than in the forward region. Note that in the case of helium, the pick-up ion density is supposed to be larger in the anti-apex direction due to Sun's gravitational focusing of LISM neutral helium, thus suggesting even larger asymmetry in ACR ${\rm He}^+$. Also large is the effect of the shock strength (and hence the slope of the spectrum) varying with position on the shock. In Kausch's model the shock is stronger (hard spectrum) close to the anti-apex direction, causing a ridge in the spatial ACR distribution.

3.3 Modulation of the ACR energy spectrum downstream of the shock

In the supersonic solar wind region the cosmic ray distribution is shaped by inward diffusion operating against convection by outflowing plasma, with adiabatic cooling due to the positive divergence of the flow. Taking into account also the energy range (the lower energies are confined to the vicinity of the shock) the modulation is then determined by the dependence of the diffusion coefficient on energy. On the contrary, downstream of the shock the diffusion and convection are both directed outward. If one is interested in the lower energy range (<102 keV) the charge exchange rate is high ( $10^{-9}\ {\rm s}^{-1}$) and strongly energy-dependent: this is the main factor in the modulation of the low-energy ACR spectrum in this region. The energy dependence of the diffusion coefficient can be seen to be less important. The high-energy (>102 keV) part of the spectrum downstream of the shock is only weakly modulated (the loss term $\beta $ is small at high energy).


  \begin{figure}
\par\includegraphics[width=8cm,clip]{h2277f13.ps}\end{figure} Figure 13: Evolution of the ACR energy spectrum in Kausch's model B (apex direction). The spectra at distances (top to bottom) of 89 (shock), 188, 305, 551 and 990 AU are shown. The dotted lines correspond to disregarding the $\sim $ $\nabla\cdot {\vec V}$ term
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.1cm,clip]{h2277f14.ps}\end{figure} Figure 14: Evolution of the ACR energy spectrum in Kausch's model B (anti-apex direction). The spectra at distances (top to bottom) of 187 (shock), 275, 380, 600 and 990 AU are shown. The dotted lines correspond to disregarding the $\sim $ $\nabla\cdot {\vec V}$ term
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.1cm,clip]{h2277f15.ps}\end{figure} Figure 15: Evolution of the ACR energy spectrum in Parker's model (apex direction). The spectra at distances (top to bottom) of 79 (shock), 181, 300, 552 and 1000 AU are shown
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=8.1cm,clip]{h2277f16.ps}\end{figure} Figure 16: Evolution of the ACR energy spectrum in Parker's model (anti-apex direction). The spectra at distances (top to bottom) of 79 (shock), 181, 300, 552 and 1000 AU are shown
Open with DEXTER

Another source of modulation is the adiabatic acceleration (deceleration) due to the divergence of the plasma flow. Although in Kausch's model (A and B) the magnitude of the divergence term is comparable with (or greater than) the charge-exchange rate $\beta $, the effect of the divergence term on modulation is small. Generally, we find that the divergence term tends to leave the pure power law spectrum unchanged. This is easy to see if the divergence term and convection are dominant in the transport equation (for a power law the divergence term is equivalent to the energy-independent loss term, which is not causing any modulation).

Figures 13 to 16 show, for Parker's and Kausch's models, the calculated energy spectra at different distances from the Sun, both in the apex and anti-apex (heliotail) direction. The strong modulation seen in the case of Parker's model is because the hydrogen number density inside the heliopause was (for purpose of comparison) assumed to be much higher (0.1  ${\rm cm}^{-3}$) than the value which follows from Kausch's simulation (about 0.02  ${\rm cm}^{-3}$ in the heliotail). The results obtained by disregarding vthe divergence term are shown by the dotted lines in Figs. 13 and 14. The modulation is predominantly due to the energy dependence of the loss rate $\beta $ (or, charge-exchange cross section). Kausch's model assumption of low outside $n_{\rm H}$ makes the modulation relatively weak in this case.

The relatively weak effect of the divergence term is important for future modelling: it implies that a reasonable description can be obtained by assuming that the modulation is solely due to the energy dependence of the loss term. As the divergence term is the only important one in which the energy derivative appears, one could in effect decouple the simulations for different energies (the divergence term can be estimated by assuming the form of the energy spectrum). It is then possible to do useful 3-D simulations without having to introduce the energy variable as an extra dimension.


  \begin{figure}
\par\includegraphics[width=8.1cm,clip]{h2277f17.ps}\end{figure} Figure 17: ENA H energy spectrum and the ACR proton shock spectrum for Kausch's model A (apex and anti-apex direction). The low energy ACR flux is higher in apex than in anti-apex direction due to larger pick-up proton density there
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=17cm,clip]{h2277f18.ps}\end{figure} Figure 18: ENA flux: the data points show the flux (units (cm2 s sr keV)-1) of mass = 1 particles in low energy channels ($\sim $58-88 keV) observed during quiet times (low ion flux) by CELIAS/HSTOF instrument. The results are the revised ones, taking the cross-calibration with ACE and IMP8 into account (Hilchenbach et al. 2000). Superimposed curves show the calculated flux (scaled up by a factor of 10) in the field of view of the instrument, of 63 keV hydrogen ENAs from ACR transcharging for different models of the heliosphere: Kausch's model A (solid line), model B (dotted line) and Parker's model (dashed line)
Open with DEXTER

3.4 ENA energy spectrum

Figure 17 shows the calculated ENA H energy spectrum for Kausch's model A together with the ACR spectrum assumed at the shock. The shape of the ENA spectrum is steeper than the ACR shock spectrum due to energy dependence of the charge-exchange rate (or cross section).

4 ENA: Observations against models

Most of the ENA data points from CELIAS/HSTOF published so far correspond to two first years of observations: for the rest of the period the published data are fragmentary. The main reason is that in 1996 and 1997 the solar activity was low, with comparatively many quiet periods of low ion flux. The ENA hydrogen events were identified as those ${\rm mass} = 1$events in low-energy channels ($\sim $55-80 keV) that occurred during quiet periods (after subtracting the accidental coincidence background). The low-energy channels are the ones for which the probability of the ion transmission through the deflecting fields is the lowest.

Figure 18 presents the flux of mass = 1 particles corresponding to the quiet-time events as a function of time. As the CELIAS/HSTOF line-of-sight has a fixed direction relative to Sun, over a year it scans all directions in the ecliptic plane. The anti-apex direction corresponds to DOY 194. In both 1996 and 1997 there are peaks of the flux intensity about that time. While in 1998 the contact with the spacecraft was lost for a time period which included the DOY 194, there are some data obtained during 1999 (first presented at IAGA 1999, Hilchenbach et al. 1999). In this case most of the data for low energy channels come from active periods. They set an upper limit rather than an estimate of the neutral flux. The data points shown in Fig. 18 are for quiet periods. Only one short quiet period was found near the anti-apex direction. The data for the corresponding period in year 2000 are absent because of the big solar flare.

The data shown in Fig. 18 differ from those published earlier (Hilchenbach et al. 1998). The reason is the cross-calibration of the CELIAS/HSTOF data with those from ACE and IMP8 which was recently carried out (Hilchenbach et al. 2000). In result, the previous estimation of the efficiency of HSTOF and in consequence of the flux level is now in doubt. The presently suggested value of the flux is higher by a factor of 10. Also affected is the relative magnitude of the 1996 and 1997 peaks, which now seem to be the same within the experimental errors (in the previous estimations, the 1996 peak was higher by a factor of 1.5-2).

The curves show the results of calculations (scaled up by a factor of 10) for Kausch's model (A and B) and for Parker's model. To fit the present HSTOF data one would have to take the ACR flux at the shock to be a factor of 10 higher than the value assumed in our simulations.

5 Conclusions and outlook

In the calculations reported here, the ACR distribution beyond the termination shock was for the first time treated using a reasonably realistic model of the heliosphere, which incorporates the very important interaction with neutral background. In particular, the calculations included a non-spherical termination shock, with parameters varying over its surface, and the effect of adiabatic energy changes in a compressive ( $\nabla \cdot {\vec V} \neq \ 0$) plasma flow. The prediction of an anisotropy of the ENA flux (highest flux from the anti-apex) was confirmed. Despite the uncertainty due to large experimental errors, this prediction is also consistent with the observations by HSTOF.

The change in the HSTOF data following the recent cross-calibration affects the comparison with the ACR ENA model in two ways. The difference in intensity between the 1996 and 1997 peaks, which was difficult to explain in the ACR ENA model (Hilchenbach et al. 1998; Czechowski et al. 1999c), is now reduced or absent. On the other hand, the predicted ENA flux intensity, which seemed previously to agree well with the data, is now too low. To restore the agreement, the parameter setting the flux intensity scale in the model (the flux of ACR protons of reference energy at the apex point of the termination shock) would have to be increased.

The ACR flux at the termination shock needed to explain the present HSTOF data would have to be of the order of 0.05-0.1 (cm2 s sr keV)-1 at $\sim $60 keV, a factor of 10 higher than our estimation based on extrapolating the power law spectrum of Stone et al. (1996) (taking into account the presence of another component, the CIR ENA, may reduce the discrepancy: see the discussion further down). The high proton flux at the shock in this energy region, leading to the ENA flux consistent with the present HSTOF data interpretation, was recently derived from the models which include pre-acceleration of the PUI by solar wind turbulence in the region upstream of the termination shock (Chalov et al. 1995; Fichtner et al. 1996; Czechowski et al. 1999b; Fahr & Lay 2000). The ENA flux in the last paper (Fahr & Lay 2000) is too high by a factor of 10 due to a numerical error (H. Fahr, private communication).


  \begin{figure}
\par\includegraphics[width=7.3cm,clip]{h2277f19.ps}\end{figure} Figure 19: ENA 10 keV hydrogen flux intensity map for Kausch's model. Heliocentric longitude/latitude plot. The anti-apex is at the longitude $\approx $74 $\hbox {$^\circ $ }$, the latitude = $-7\hbox {$^\circ $ }$. The elongation of the image is a result of the difference between the vertical and horizontal scale
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=7.3cm,clip]{h2277f20.ps}\end{figure} Figure 20: ENA 63 keV hydrogen flux intensity map for Kausch's model. Heliocentric longitude/latitude plot. The anti-apex is at the longitude $\approx $74 $\hbox {$^\circ $ }$, the latitude = $-7\hbox {$^\circ $ }$. The elongation of the image is a result of the difference between the vertical and horizontal scale
Open with DEXTER

The main problem in interpreting the observations from the point of view of the ACR-ENA model is the uncertainty in the predicted flux intensity scale. The low energy ACR flux at the shock can only be estimated by de-modulation of the observed spectra and extrapolation to low energy. If the observations in the region close to the post-shock pick-up ion energy (possibly about 10 keV/nucleon, which is at the lower limit of the energy range of INCA/Cassini and HENA/IMAGE) would be available, the theoretical estimation of the ion flux intensity would be easier and the conclusions drawn from the observations could be stronger.

Although the experimental errors are too large to make a clear identification, the shift of the 1996 ENA peak from the anti-apex direction, also visible in the 1997 peak and possibly repeated in 1999 (although there is only a single quiet period data point available for the 1999 peak, the shift seems consistent with the upper limits on the neutral flux derived from the observations during high ion flux periods) could be an indication of a shift in the heliotail direction towards larger ecliptic longitude. Such a shift may be produced by the interstellar magnetic field inclined with respect to the apex-antiapex axis (Fahr et al. 1988; Ratkiewicz et al. 1998; Czechowski & Grzedzielski 1998). The data are insufficient to make anything more than an order-of-magnitude estimation of the shift ( $10\hbox{$^\circ$ }{-}15{\hbox{$^\circ$ }}$). Unfortunately, the only available models of the heliosphere able to deal with the interstellar field of arbitrary orientation (Ratkiewicz et al. 1998; Pogorelov & Matsuda 1998) include no neutral component, which makes any estimation of the interstellar field based on the observed shift uncertain.

There is a possibility (Kota et al. 2000) that the anti-apex peak of the CELIAS/HSTOF ENA signal is related to the neutral helium density enhancement in the anti-apex direction due to gravitational focusing by the Sun (this would also give a natural explanation of the shift in the position of the peak). The observed ENA would then originate not from the ACR ions but from the ions accelerated in interplanetary space, presumably by the co-rotating interaction regions (CIRs). The ACR ENA from the heliotail would then have to be looked for at lower energy, where the charge-exchange between protons and helium is suppressed. In any case, the ENA of the CIR origin should be considered as an additional contribution, besides the ACR ENA, to the anti-apex flux intensity peak.

In addition to CELIAS/HSTOF, two other instruments with capability to observe the heliospheric ENA are in operation: INCA on Cassini (Krimigis et al. 2000) and HENA on the IMAGE satellite (Mitchell 2000). The observations by INCA and IMAGE need not be restricted to the vicinity of the ecliptic but can cover other sectors of the sky, depending on the spacecraft orientation, so that in principle a map of the ENA flux intensity for all the sky could be obtained. In Figs. 19 and 20 our results for the ENA flux are presented as density plots covering all the sky. Our model being axially-symmetric, we cannot examine the effect of a latitudinal dependence of the solar wind parameters on the ACR and ENA fluxes. The heliotail should be visible as a peak in the ENA flux, well-defined particularly at higher energy. At low energy this peak is wider. Possibility of using the ENA for observing the effect of the outer magnetic field on the shape of the heliotail was considered in Czechowski & Grzedzielski (1998).

Acknowledgements
A.C. acknowledges support by KBN grant 2 P03C 004 14. The work at University of Arizona was supported in part by NSF grant ATM9727080 and NASA grant NAG57966. A.C. and H.F. are grateful for financial support within the framework of a Polish-German cooperation (project 436 POL 113/80/0).

Appendix A: Boundary conditions

The region in space in which Eq. (1) is solved is between the termination shock and the outer boundary (a sphere of radius 990 AU for the Kausch model). The domain of integration is the direct product of the region in space with the energy interval $(E_{\min},E_{\max})$ (for most calculations (10,103) keV). The boundary conditions for the distribution $f({\vec x},p)$ are specified at the boundaries in space (for the whole energy interval) and at the limits in energy (for the whole region in space).

At the shock surface we specify the function $f({\vec x},p)$. This requires an assumption about the ACR shock spectrum. The acceleration by a quasi-perpendicular shock involves a drift along the shock surface (Jokipii 1982, 1987; Steenkamp & Moraal 1993), but here we are interested in the low-energy particles. For those, the diffusion process is quite slow, and we may approximate the accelerated spectrum with that obtained from simple planar-shock theory, taking local values of the shock parameters. That this is reasonable is readily seen by noting that the relevant diffusion scale $\kappa /V$ is much less than the macroscopic length scales, which are of the order of tens of AU For a locally planar shock, it is a simple matter (see, e.g., Drury 1983) to show that the accelerated particles have the spectrum $f(p) \propto (1/\Delta V) p^{-q}$ where $\Delta V$ is the (normal) change of the velocity across the shock, and where q = 3 r/(r - 1) (with r the shock compression). This should be weighted by the factor representing the source of the ACR ions, which can be taken as given by the number density $n_{\rm PUI}$of the corresponding PUIs just upstream of the shock. We use

\begin{displaymath}f({\vec x},p)\vert _{\rm shock}={\rm const}\,\, n_{\rm PUI} \frac{1}{\Delta
V}(E/E_0)^{-q/2}.
\end{displaymath} (A.1)

For E0, the value at which the power law spectra converge, we use 3.5 keV (characteristic pick-up proton energy; in some of our earlier calculations we used 1 keV).

Most of our calculations assume one of the following cases: a fixed power law in energy and no variation of $f({\vec x},p)$ over the shock; a power law dependent on the local shock strength and an amplitude dependent on the flow velocity change across the shock (with the factor of $n_{\rm PUI}$ included or not). The scale of f at the shock as well as the fixed power law was obtained using an extrapolation of the results of Stone et al. (1996).

The boundary condition at large distance should ensure $f_{r\to \infty}\to 0$. As the background flow solution which we use is restricted to a sphere of a finite radius (990 AU), we require that the radial derivatives of f match with the approximate analytical solution (for which $f\to 0$). This procedure was checked against other choices.

The boundary conditions in energy require a comment. We specify the slope ${\rm d}\ln f/{\rm d}\ln p$ of the energy spectrum at $E_{\max}$ and $E_{\min}$, for the whole region in space. A fully consistent choice would require knowledge of the solution for the modulated spectrum. At the high energy limit, $E_{\max}$, the ACR number density is low and an uncertainty in the boundary condition for the spectrum is not likely to affect the solution at lower energy. This is not true for the low energy limit. We found, however, that the modulation of the energy spectrum at the low energy limit is weak, because the loss rate $\beta $ at $E\sim E_{\min}$ is relatively flat as a function of energy. In result, choosing a fixed power law spectrum (for example, an average over the shock surface) at minimum energy will not induce large errors in the solution at intermediate energy, where most of the modulation occurs. We have checked the effect of boundary conditions by comparing the solutions for different choices of the boundary values of energy ((10,103) keV and (1,106) keV) and found a good agreement in the intermediate region.

References

 


Copyright ESO 2001