A&A 429, 1069-1080 (2005)
DOI: 10.1051/0004-6361:20041348
V. Izmodenov ^{1,} ^{2} - Y. Malama ^{2} - M. S. Ruderman ^{3}
1 - Lomonosov Moscow State University, Department of
Aeromechanics and Gas Dynamics, Faculty of Mechanics and
Mathematics, and Institute of Mechanics, Moscow 119899, Russia
2 - Institute for Problems in Mechanics,
Russian Academy of Sciences, Prospect Vernadskogo 101-1, Moscow
117526, Russia
3 - Department of Applied Mathematics, University
of Sheffield, Hicks Building, Hounsfield Road, Sheffield S3 7RH,
UK
Received 25 May 2004 / Accepted 3 September 2004
Abstract
We present results of a new time-dependent kinetic model
of the H atom penetration through the solar wind - interstellar
medium interaction region. A kinetic 6D (time, two dimensions in
space, and three dimensions in velocity-space) equation for
interstellar H atoms was solved self-consistently with
time-dependent Euler equations for the solar wind and interstellar
charged components. We study the response of the interaction
region to 11-year solar cycle variations of the solar wind dynamic
pressure. It is shown that the termination shock location varies
within 7 AU, the heliopause variation is 4 AU, and the
bow shock variation is negligible. At large heliocentric
distances, the solar cycle induces 10-12% fluctuations in the
number density of both primary and secondary interstellar H atoms
and atoms created in the inner heliosheath. We underline the
kinetic behavior of the fluctuations of the H atom populations.
Closer to the Sun the fluctuations increase up to 30-35% at 5 AU
due to solar cycle variation of the charge exchange rate. Solar
cycle variations of interstellar H atoms in the heliospheric
interface and within the heliosphere may have major importance for
the interpretation of H atom observations inside the heliosphere.
Key words: Sun: solar wind - interplanetary medium - ISM : atoms - Sun: activity - hydrodynamics
More than 30 years (three solar cycles) of solar wind observations show that its momentum flux varies by factor of 2 from solar maximum to solar minimum (Gazis 1996; Richardson 1997). It was shown theoretically that such variations of the solar wind momentum flux strongly influence the structure of the heliospheric interface - the region of the solar wind interaction with the local interstellar medium (e.g., Karmesin et al. 1995; Wang & Belcher 1998, 1999; Baranov & Zaitsev 1998; Zank 1999; Zaitsev & Izmodenov 2001; Scherer & Fahr 2003a,b; Zank & Mueller 2003; Izmodenov & Malama 2004a,b).
Most global models studying the solar cycle effects ignore the interstellar H atom component or took this component into account by using simplified fluid (Scherer & Fahr 2003a,b) or multi-fluid (Zank & Mueller 2003) approximations. These simplifications were done because it is difficult to solve 6D (time, two dimensions in space, and three dimensions in velocity-space) kinetic equation for the interstellar H atom component. A kinetic description of interstellar atoms is necessary due to their large mean free path comparable with the size of the heliospheric interface. Recently, it has been shown explicitly by Izmodenov (2001), Izmodenov et al. (2001) that the velocity distribution function of interstellar atoms is not Maxwellian. Fluid or multi-fluid approaches assume the velocity distribution of H atoms to be Maxwellian or the sum of several Maxwellian distributions. This assumption introduces additional uncertainties in the model. (To estimate these uncertainties one needs to compare a multi-fluid model with a kinetic model for each specific set of model parameters.) Consequently, the multi-fluid approaches adopted for H atoms in the heliospheric interface may result in misleading interpretations of observational data. At the same time, the momentum flux variations of the solar wind may have a significant effect on the distribution of interstellar H atoms in the heliosphere due to coupling of the atom and plasma components by charge exchange. Therefore, it is necessary to study solar cycle effects by solving 6D kinetic equation self-consistently with fluid plasma equations. In this paper we present the results of such a model.
The qualitative pattern of the heliospheric interface is shown in Fig. 1. The solar wind plasma decelerates and turns to the tail at the termination shock (TS), while the interstellar plasma decelerates and turns outward of the axis of symmetry at the bow shock (BS). The heliopause (HP), which is a contact discontinuity, separates these two plasma flows. Hydrogen atoms newly created by charge exchange have the velocities of their ion partners in the charge exchange collisions. Therefore, the parameters of these new atoms depend on the local plasma properties. It is convenient to distinguish four different populations of H atoms: 1) atoms created in the supersonic solar wind; 2) atoms originating in the heliosheath and known as heliospheric ENAs (Gruntman et al. 2001); 3) atoms created in the disturbed interstellar wind; 4) original (or primary) interstellar atoms.
In this work we advance the heliospheric interface model developed by the Moscow group (e.g. Baranov & Malama 1993; Izmodenov et al. 1999; Alexashov et al. 2000, 2004a,b; Myasnikov et al. 2000; Zaitsev & Izmodenov 2001; Izmodenov & Alexashov 2003; Izmodenov et al. 2003a,b; Izmodenov et al. 2004; for reviews see Izmodenov 2003, 2004) by introducing 11 year sinusoidal variations of the solar wind in the model. The goal of this paper is to explore theoretical aspects of non-stationary interaction of the solar wind and the local interstellar cloud. The results presented here cannot be directly applied to interpretation of observational data.
Figure 1: Structure of the heliospheric interface, the region of the interaction of the solar wind and the Local Interstellar Cloud. | |
Open with DEXTER |
We consider all the plasma components (electrons, protons, pickup ions, solar wind alpha particles and interstellar He ions) as one fluid with total density and bulk velocity . It is assumed that all ionized components have the same temperature T. Although this assumption cannot be made in the case of the solar wind, the one-fluid model is based on mass, momentum and energy conservation laws and predicts the plasma bulk velocity and locations of the shocks very well.
The plasma is quasi-neutral; i.e., for interstellar plasma, and for the solar wind. The interstellar and interplanetary magnetic fields are ignored in the paper. Protons interact with H atoms by charge exchange. While the interaction of interstellar H atoms with protons by charge exchange is important, for helium ions the process of charge exchange is negligible because of small cross sections for the charge exchange of helium atoms. Hydrodynamic equations for the charged component are solved self-consistently with the kinetic equation for the H atom component.
The governing equations for the charged component are the
following:
The right hand parts, q_{1},
and q_{3}, are
sources of mass, momentum and energy due to charge exchange of H atoms and protons, photoionization and electron impact ionization
(Baranov &
Malama 1993):
The velocity distribution of H atoms
is calculated from the linear kinetic
equation:
The boundary conditions for Euler equations are the following. We
assume that the solar wind is spherically symmetric at the Earth
orbit. It makes our model axissymmetric. We also assume that, at
the Earth orbit, the solar wind number density oscillates
harmonically, while the bulk velocity and temperature remain
constant:
The following solar wind parameters averaged over a few solar cycles were used: cm^{-3}, km s^{-1}. We also used the following parameters of the interstellar gas in the unperturbed interstellar medium at the outer boundary: km s^{-1}, T K, cm^{-3}, cm^{-3}. These particular values of the interstellar velocity and temperature were chosen on the basis of the recent observations of the interstellar He atoms by GAS/Ulysses (Witte et al. 1996; Witte 2004; Gloeckler et al. 2004). The choice of and is based on our analysis of the Ulysses pickup ion measurements (see, e.g., Izmodenov et al. 2003a,b, 2004).
Now we explain why we assumed "idealized'' harmonic perturbations at the Earth orbit rather than the realistic perturbations of the solar wind parameters obtained from observations. The statistical Monte-Carlo method used to obtain periodic solutions of the kinetic equation requires us to fix the time-period. The non-linear nature of the system may lead to the interaction of the external 11-year fluctuations with the internal oscillations of the heliospheric interface. As a result, oscillations with periods different from 11 years may appear in the self-consistent solution of the governing Eqs. (1)-(3) with the boundary conditions (4). One of the main objectives of our study was to verify if such oscillations do appear. To do this we increased the time-period of our Monte-Carlo calculations by 6 times. If oscillations with periods different from 11 years are present in our solution we would be able to determine these periods. If we use the realistic perturbations of the solar wind parameters that are, in general, not periodic, as the boundary conditions, it would be rather difficult to detect the presence of oscillations with periods different from 11 years.
To solve the Euler equations self-consistently with the kinetic equation we used as iterative procedure as suggested by Baranov et al. (1991) for the stationary model. In the first step of this iterative procedure the Euler equations with the constant source terms q_{1}, and q_{3} were solved with the use of the boundary conditions (4). The source terms were taken from the stationary solution with the average solar wind parameters. We performed calculations over 300 solar cycles. As a result, we got the distribution of the plasma parameters. We analyzed this distribution and found that there is the 11 year periodicity only.
In the second step we solved the kinetic equation by a Monte Carlo method with splitting of trajectories (Malama 1991). To increase the statistical efficiency of the method, we assumed periodicity with the time period = 66 years. To minimize statistical errors we averaged the statistical results over year. When doing so we used a distribution of the plasma parameters for the last 66 years obtained in the first step. As a result, we obtained the periodic (66 year) q_{1}, and q_{3} source terms. In the third step we solved the Euler equations with the boundary conditions (4) and the periodic source terms obtained in the second step. Again, we performed gas dynamic calculations over 300 solar cycles. Analysis of the plasma distributions shows the 11 year periodicity only. Then we solved the kinetic equation by the Monte Carlo method with the distribution of the plasma parameters for the last 66 years obtained in the third step. We continued this process of iteration until the results of two subsequent iterations were practically the same.
Our method allows us to obtain the self-consistent solution of the system of Euler Eqs. (1) and the 6D kinetic Eq. (3) with the boundary conditions (4). Since the uniqueness of the solution for this system is not proven, we cannot exclude that other solutions of this system of equations may exist. These solutions may have periods different from 11 years. The numerical method used to solve Euler equations does not need any restricting assumptions. Remarkably, as it is reported in the next section, our numerical solution does not contain oscillations with periods different from 11 years.
Figure 2: Time variations of the heliocentric distances to the termination shock, bow shock and the heliopause in the upwind direction, and to the termination shock in the downwind direction. Bottom plot shows variations of the solar wind momentum flux, , with time. Dashed curves correspond to the solution of the problem with the "broken'' first of six solar cycles, when we increased the solar wind ram pressure by factor of 1.5 during the first 11 years of our 66-year periodic calculations. | |
Open with DEXTER |
In addition to the "ideal'' solar cycle calculations we carried out some complimentary calculations. We increased the solar wind ram pressure by a factor of 1.5 during the first 11 years of our 66-year period (Fig. 2, bottom plot). This study was inspired by the fact that observations of the solar wind are restricted to the recent years. We wanted to understand how the solar wind conditions in the past, when they were not observed, influence the heliospheric interface and observational quantities at present and in the future.
The variations of the heliocentric distances to the termination shock, heliopause and the bow shock are shown in Fig. 2. The discontinuities vary with a 11-year time-period under the action of 11-year fluctuations of the solar wind dynamic pressure at the inner boundary of our computational grid. The termination shock oscillates from its minimal distance of 93 AU, which is reached in the last (11th) year of the "ideal'' solar cycle, to its maximum distance of 107 AU, which is reached during the fourth year of the cycle. Fluctuations of the TS are bigger in the downwind direction than in the upwind direction. By the upwind direction we mean the direction that is opposite to the direction of the Sun-LIC relative motion. The downwind direction coincides with the Sun-LIC velocity vector. The variation of the TS in the downwind direction is 25 AU from its minimum value of AU during the third year of the cycle to its maximum value of 188 AU in the 9th year of the cycle. In the upwind direction the most distant position of the termination shock is reached 1.5-year after the maximum of the solar wind dynamic pressure at 1 AU. The variations of the solar wind dynamic pressure are shown at the bottom of Fig. 2 for convenience. The phase of downwind fluctuations of the TS is shifted by years compared the phase of the upwind fluctuations. It is interesting to compare our results with the results obtained by Zank & Mueller (2003) and Scherer & Fahr (2003a,b) on the basis of (multi-) fluid descriptions of the interstellar H atoms (Table 1).
Table 1: Results of calculations.
The strength of the TS has important consequences for spectra of anomalous cosmic rays (ACRs) because the velocity jump at the TS is related to the spectral index of ACRs that determines the variation of intensity of the cosmic rays jwith energy E: . We have computed the variation of the velocity jump at the TS. In the upwind direction the jump of the plasma velocity at the TS, or, in other words, the strength of the TS, varies from its minimum value of 2.92 to its maximum value of 3.09. This corresponds to the variation of from 1.28 to 1.22. The strength of the TS varies from 2.92 to 3.17 in the downwind region. This strength variation is translated into the variation of from 1.28 to 1.19.
The heliopause fluctuates with smaller amplitude as compared to the termination shock. It varies from 169 AU, which is reached during 4th year of the solar cycle, to 173 AU reached in the 9th year of the solar cycle. The distance to the heliopause averaged over the solar cycle is 171 AU. This coincides with the stationary solution. The solar-cycle induced fluctuations of the BS is less than 0.1 AU in the upwind direction. The fluctuations are not visible in Fig. 2. The distance to the BS averaged over the solar cycle is 308 AU, while this distance is 311 AU in the case of the stationary solution.
Dashed curves in Fig. 2 correspond to solution of the problem with the "broken'' first of six solar cycles, when we increased the solar wind ram pressure by a factor of 1.5 during the first 11 years of our 66-year periodic calculations (bottom plot in Fig. 2). The termination shock in the upwind direction "feels'' the increase of the solar wind dynamic pressure during approximately 4-5 years after the increase ended. In the downwind direction the "feeling'' is somewhat longer and lasts another solar cycle. The second maximum of the TS both in the upwind direction (dashed curve on the top plot in Fig. 2) and in the downwind direction is closer to the Sun compared to the subsequent maxima.
The post-reaction of the heliopause to the 50% increase of the solar wind dynamic pressure is much longer compared to the reaction of the termination shock. The heliopause does not return to its periodic fluctuations even at the end of the 66-year time-period. One can see from the figures that the heliocentric distances to the termination shock, heliopause and bow shock are always larger in the case of a "broken'' solar cycle compared to our "regular'' solar cycles. This is related to the fact that the solar wind ram pressure averaged over 66 years is 8% greater compared to the "regular'' cycle calculations. The effect is the most pronounced for the bow shock. It appears that 66 years are not enough for the BS to relax to its "regular''-cycle position. As a result, the BS is AU away for the "broken''-cycle calculations compared to the "regular'' cycle.
Figure 3: Interstellar plasma number density, velocity, pressure and temperature as functions of the heliocentric distance for two different moments of time: t_{1} = 1 year (curves 1), t_{2}= 6 year (curves 2). Stationary solution (curves 3) and averaged over 11 years time-dependent solution (dots) are shown. | |
Open with DEXTER |
Plasma parameters undergo 11-year fluctuations in the entire computation region. However, the wave-length of the plasma fluctuations in the solar wind is apparently larger compared to the distances to the TS and HP. This means that time snap-shots of the distributions of plasma parameters (density, velocity and temperature) are not qualitatively different from stationary solutions. The situation is different in the outer heliosheath, which is the region between the HP and BS. 11-year periodic motion of the heliopause produces a number of additional weak shocks and rarefaction waves (Baranov & Zaitsev 1995). The amplitudes of these shocks and rarefaction waves decrease while they propagate away from the Sun due to the increase of their surface areas, interaction between the shocks and rarefaction waves, and the dissipative attenuation of the shocks. To resolve the wave structure we increased the resolution of our computational grid by three times in the region. We have checked also that an additional increase of the resolution of our computational grid does not change the results. Figure 3 presents distributions of plasma density, velocity, pressure and temperature as functions of the heliocentric distance in the upwind direction at two different moments (t_{1} = 1 year (curves 1) and t_{2} = 6 year (curves 2). It is seen that the characteristic wavelength in the region is AU. Long-scale waves are also seen in plasma distributions in the post-shocked plasma of the downwind region (Fig. 3, left column). Amplitudes of the waves are much less than in the upwind direction and the wavelength is AU.
Figure 4: Number densities ( top row), bulk velocities ( middle row) and kinetic temperatures ( bottom row) of primary and secondary interstellar atom populations ( left column) and atoms created in the supersonic solar wind and inner heliosheath ( right column) in the upwind direction as functions of the heliocentric distance. Dots, which represent the stationary solution, are practically coincident with solid curves, which represent the 11 year averaged time-dependent solution. | |
Open with DEXTER |
Figure 3 shows a comparison of the 11 year average distributions of interstellar plasma parameters (dots) with those obtained from stationary solution. The stationary calculations were performed with exactly the same inner and outer boundary conditions as used in the time-dependent calculations. At the Earth's orbit we assume 11 year average values of the solar wind density. It is seen that the two distributions practically coincide. The congruence with the stationary solution is additional evidence of sufficient resolution of our computational grid and the lack of significant numerical dissipation in our numerical calculations. Our results contradict the conclusion made by Zank & Müller (2003) that "the shocks provide additional heating in the heliotail and outer heliosheath''. According to our results the heating is very small and it is not noticeable in our calculations. However to draw conclusion on the plasma heating in the heliosheath, Zank & Müller (2003) had compared their time-dependent results with a stationary model that assumed smaller solar wind dynamic pressure compared with the 11-year averaged value. Therefore, the observed heating could be 1) due to the shock heating; 2) due to different boundary conditions. Additional multi-fluid study is needed to distinguish between these two mechanisms.
To better understand why the variation of the solar wind parameters almost does not disturb the bow shock, we studied the propagation of perturbations in the outer heliosheath analytically. We studied the propagation of perturbations only near the symmetry axis, and assumed that the wavelength is small in comparison to the characteristic scale of inhomogeneity. The latter assumption enabled us to use the WKB approximation. In addition, we neglected the interaction of the plasma perturbations with the H atoms. Then, using the reductive perturbation method, we derived the governing equation for the plasma perturbations. This equation is a generalization of the nonlinear equation used in nonlinear acoustics for the description of sound waves. We solved this equation assuming the boundary conditions at the heliopause corresponding to the harmonic oscillation of the heliopause with period of 11 years and amplitude of 2 AU. Our main result is that, due to the nonlinear steepening, the shock forms in the wave profile at about the middle of the distance between the heliopause and the bow shock. The wave energy dissipation in this shock causes strong attenuation of the perturbations in their way from the heliopause to the bow shock. As a result, the wave amplitude at the bow shock is about 3 times smaller than that predicted by the linear theory. Since the wave energy flux is proportional to the amplitude squared, this implies that almost 90% of the wave energy is dissipated in the shocks. On the basis of this result we conclude that the main reason why the solar cycle variation almost does not disturb the bow shock is that the perturbations propagating from the heliopause to the bow shock are strongly attenuated due to dissipation in the shocks.
Comparison with numerical results reveals that the attenuation of the perturbations obtained in the numerical simulation is even stronger than that predicted by the analytical solution. The most probable cause of this difference is that the interaction between the plasma perturbations and the H atoms provides an additional wave dissipation.
Using the analytical solution describing the wave propagation in the outer heliosheath, we estimated the rate of the plasma heating due to wave dissipation. We found that the mean temperature of the plasma in the outer heliosheath can be increased by about 280 K during one solar cycle. The plasma heating due to the wave dissipation is compensated by the energy loss due to convective plasma motion and due to the interaction between the plasma and the H atoms as seen from the results of our numerical calculations.
The main advantage of our model compared to previously published multi-fluid models (Scherer & Fahr 2003; Zank & Mueller 2003) is a rigorous kinetic description of the interstellar H atoms. Charge exchange significantly disturbs the interstellar atom flow penetrating the heliospheric interface. The atoms newly created by charge exchange have velocities of their ion partners in the charge exchange collisions. Therefore, the velocity distribution of these new atoms depends on the local plasma properties at the place of their origin. As it was discussed in the introduction, it is convenient to distinguish four different populations of H atoms depending on the region in the heliospheric interface where the atoms originate. Figure 4 compares the distributions of the populations of H atoms obtained by the stationary model (dots) with the time-dependent solution averaged over 11 years. For the plasma component there is no noticeable difference between these two distributions. Although we present only distributions in the upwind direction, our conclusion remains valid for all the computational domain. The stationary distributions of the H atom parameters for directions different from upwind could be found in our earlier papers (see, e.g., Izmodenov 2000; Izmodenov et al. 2001).
Figure 5: Time-variation of the number densities of primary and secondary interstellar atoms ( top panels), H atoms created in the inner heliosheath and H atoms created in the supersonic solar wind ( bottom panels) as functions of heliocentric distance for two different moments in the solar cycle. | |
Open with DEXTER |
To evaluate time-dependent features in the distribution of H atoms in the heliospheric interface we plot the number densities of the four populations of H atoms normalized to the densities obtained in the stationary solution. By doing this we suppressed spatial gradients of the densities, which are apparently larger than the time-variations of the densities. Figure 5 shows the normalized densities for two different years of the solar cycle. Solid curves correspond to t_{1} = 1 year and dashed curves to t_{2} = 6 year. It is seen that the variation of the density is within 5% of its mean value for the primary and secondary interstellar populations, and for the atoms created in the inner heliosheath. Closer to the Sun, for distances less than 10 AU, the amplitude of the fluctuations increases up to 15%. The variation of the number density of H atoms created in the supersonic solar wind is 30% about its mean value.
Figure 6 shows the time-variation of the number densities, bulk velocities and kinetic temperatures of three populations of H atoms at different heliocentric distances in the upwind direction. All parameters are normalized to their initial values at t= 0. Clear 11-year periodicity is seen for the number densities of the atoms. Deviation from the exact 11-year periodicity is related to the errors of our statistical calculations, which are 2-3%. Less than 10% variation (from maximum to minimum) is seen for number densities of all populations at distances greater than 10 AU. At 5 AU the variations are of the order of 30%. Variations of the bulk velocity and kinetic temperature are negligibly small for both primary and secondary interstellar populations. However, the bulk velocity and kinetic temperature of atoms created in the inner heliosheath vary with the solar cycle by 10-12%. This is related to the fact that most of the H atoms of the latter population are created in the vicinity of the heliopause (Fig. 4) and they reflect long wavelength plasma variations in this region. The correlation of parameters of the H atom population created in the inner heliosheath with the plasma parameters in the vicinity of the heliopause is illustrated in Fig. 7.
Figure 6: The number densities ( second column from the left), bulk velocities (second column from the right) and kinetic temperatures ( right column) of the primary (solid curves) and secondary (dashed curves) interstellar atom populations and the atoms created in the inner heliosheath (curves with diamonds) at different heliocentric distances in the upwind direction as functions of time. For comparison, the number density of the plasma is shown ( left column). All parameters are non-dimensionlized to their values at t=0. | |
Open with DEXTER |
Figure 7: Comparison of variations of density, velocity and kinetic temperature of H atoms created in the inner heliosheath (left column, lines with triangles) with those plasma parameters (right column). The variations are shown at R = 160 AU in the vicinity of the heliopause in the upwind direction for "broken'' solar cycle calculations, where we increase the solar wind ram pressure by factor of 1.5 during the first 11 years of chosen 66-year time-period. | |
Open with DEXTER |
It is important to note that number densities of all three components of H atoms fluctuate in the same phase. Such coherent behavior of fluctuations remains in the entire supersonic solar wind region (R<90 AU) for the three populations of H atoms and in the inner heliosheath for the primary and secondary atoms. The reason for such coherent behavior of the variations of H atom densities becomes evident when we compare them with the plasma density variations (Fig. 6, left column). The two quantities vary almost in anti-phase. Apparently, such a correlation is only possible when temporal variations of the H atom densities are caused by variation of the local loss of the neutrals due to charge exchange and ionization processes. The local fluctuations are not transported over large distances because the velocities of individual atoms are chaotic and their mean free path is large.
However, coherent fluctuations of different populations of H atoms disappear in the regions where the populations originate, and the process of creation dominates the losses. Indeed, in the inner heliosheath (for example, at 160 AU in the upwind direction as shown in Fig. 6) fluctuations of number density of H atoms created in this region are shifted with respect to the coherent fluctuations of the primary and secondary interstellar atom populations, and are in phase with the variation of the proton number density near the heliopause. Variations of the secondary interstellar atom populations are in anti-phase with variations of primaries in the outer heliosheath (see, R=190 AU in Fig. 6) and almost in phase with plasma fluctuations in the region. Again, the creation processes are dominant in the outer heliosheath for the population of the secondary interstellar atoms.
Finally, it is important to note that the behavior of the H atom populations in the heliospheric interface has kinetic nature. Variations of the atom parameters are determined by the loss and creation processes rather than by the convection and pressure gradient terms as it would be in the fluid description. The fluid description is valid if the Knudsen number , where l and L are the mean free path of the particles and the characteristic spatial scale of the problem, respectively. For the stationary problem the distance between the HP and BS, which is approximately 100 AU, can be chosen as L. The mean free path of H atoms in the region is AU. Therefore, . The results obtained on the basis of the kinetic and fluid descriptions where compared by Baranov et al. (1998) and Izmodenov et al. (2001). This comparison has shown explicitly that the velocity distribution function of H atoms is non-Maxwellian everywhere in the interface. For the time-dependent problem considered in this paper the characteristic size, L, is determined as a half of the wavelength of plasma fluctuations. In the region between the HP and BS AU as follows from Fig. 3. Therefore, and a fluid description is even less appropriate than for the stationary model. The fact that the fluid description is inappropriate for the atom motion is the most probable cause of the big discrepancy between our results and the results obtained by Zank & Mueller (2003) and Scherer & Fahr (2003) who used the multi-fluid and fluid approaches, respectively. It is interesting to note the qualitative difference between the Zank & Mueller (2003) and Scherer & Fahr (2003) results. The reason for this discrepancy could be once again the different description of H atoms used in these two papers. Scherer & Fahr (2003) used a one-fluid description for H atoms, while Zank & Mueller (2003) used a three-fluid description for H atoms in the interface.
We developed a non-stationary self-consistent model of the heliospheric interface and used it to explore the solar cycle variations of the interface. We obtained the periodic solution of the system of Euler Eqs. (1) for plasma, and the kinetic Eq. (2) for interstellar H atoms with the periodic boundary conditions (4) for the solar wind at the Earth's orbit. The period of the solution is 11 years.
Our basic results for the plasma component confirm the results obtained previously:
Figure 8: Comparison of the time-dependent self-consistent model results (solid curves) with the solution of the time-dependent Euler equation with the source terms (3) taken from the corresponding stationary solution (connected dotted). The figure presents the distribution of plasma number density in the upwind direction for two different moments of the solar cycle. The difference between two models is a few percent. | |
Open with DEXTER |
Acknowledgements
We thank Prof. A. I. Khisamutdinov for a fruitful discussion regarding periodicity of the global heliospheric interface structure and Monte Carlo methods. The calculations were performed by using the supercomputer of the Russian Academy of Sciences. This work was supported in part by INTAS Award 2001-0270, RFBR grants 04-02-16559, 04-01-00594, RFBR-NSFC grant 03-01-39004, NASA grant NNG04GB80G and International Space Science Institute in Bern in the frame of the ISSI team "Physics and Gas Dynamics of the heliotail'' "Determination of H atom parameters of the LIC from within the heliosphere (PI - E. Moebius)''.
We also assume that the wave motion is in the radial direction and depends on x only. Then the wave motion is governed by the system of gasdynamic equations with the right-hand sides describing the interaction between the plasma and the neutral gas.
We consider the heliopause oscillation as an external driving
force and put the boundary condition
An important quantity for our analysis is the sound speed . In what follows we take .
Let us now derive an approximate equation governing the propagation of nonlinear waves driven by the heliopause oscillation in the outer heliosheath. The typical sound speed in the outer heliosheath is 20 km s^{-1}. The wave period is 11 years s. This gives the typical wavelength of about AU. Since the outer heliosheath size along the symmetry axis is about 150 AU, we conclude that the characteristic wavelength is smaller that the characteristic scale of inhomogeneity. This observation enables us to use the WKB approximation.
The propagation of nonlinear sound waves in homogeneous as well as in inhomogeneous media has been extensively studied in nonlinear acoustics (e.g. Rudenko & Soluyan 1977) and in applications to the solar atmosphere (see, e.g., Narain & Ulmschneider 1990, and references therein). In particular, the equation governing the wave propagation has been derived under various assumptions about background quantities. However, for the reader's convenience, we will briefly outline the derivation of the governing equation in what follows.
Our aim is to derive an equation that describes the combined
effect of the nonlinearity and inhomogeneity. In accordance with
this we assume that the small parameter characterizing the effect
of inhomogeneity is of the same order as the small parameter
characterizing the effect of nonlinearity, which is .
In other words, we assume that the ratio of the characteristic
scales of perturbation and inhomogeneity is .
To show
this explicitly we introduce the stretching variable
and assume that ,
u_{0} and p_{0} are
functions of .
In line with the WKB method we also
introduce the running variable
,
where
is the
eiconal. In the new variables we rewrite the system of governing
equations as
In the zero order approximation we collect terms of order unity in
Eqs. (A.2)-(A.4). As a result we obtain the
equations for the background quantities ,
u_{0} and
p_{0}. The solution of these equations is shown in
Fig. A.2, where the background density, ,
velocity, u_{0}, and temperature, T_{0}, are shown as
functions of x. The background pressure p_{0} is calculated
using the relation
,
where
is
the Boltzmann constant and
the proton mass.
In the first order approximation we collect terms of order
in Eqs. (A.2)-(A.4). Since we
have assumed that perturbations do not affect q_{1}, q_{4} and
q_{5}, these quantities do not contribute in the equations of
any order approximation starting from the first one. As a result,
we obtain that the equations of the first order approximation
constitute a system of homogeneous algebraic equations with
respect to
,
and
.
It has a non-trivial solution only
if its determinant is zero. This condition results in the equation
determining the eiconal
Let us return to the original independent variables, use the
approximation
,
and make
the variable substitution
Equation (A.19) determines U as an implicit function of t and x. In Fig. A.3 the evolution of the wave
shape with the distance x from the heliopause is shown. The
quantity U is displayed as a function of
for three different values of x. The upper panel corresponds
to x = 0, and U is a sinusoidal function of .
When
x increases, the profile starts to steepen due to the action of
nonlinearity. Eventually, this steepening results in a gradient
catastrophe, which is the appearance of infinite gradient at one
particular point of the wave profile as shown in the middle panel
of Fig. A.3.
Figure A.3: Evolution of the shape of the wave profile with the distance from the heliopause. is shown as a function of . The upper, middle and lower panes correspond to x = 0, AU, and x = 70 AU. In the lower panel the dashed line shows the unphysical branch of the multivalues function determined by Eq. (A.19). The vertical solid line is the shock. |
In Fig. A.4 the dependence of
on x is shown
for
,
where
AU is the
distance between the heliopause and the bow shock.
Figure A.4: The dependence of the velocity jump at the shock on the distance xfrom the heliopause. The velocity is given in km s^{-1}, and the distance x in AU. |
Let us now consider the plasma heating due to the wave energy dissipation. The standard approach to solving this problem is as follows. There are different timescales in the system. The first one is the wave period. The second is the characteristic time of variation of the background quantities caused by the additional momentum input due to the wave pressure and the energy input due to wave dissipation. To take the two different timescales into account we introduce two different times, "fast'' and "slow''. Then we average governing equations with respect to the fast time. As a result we obtain the system of equations describing the slow evolution of the background quantities. This system is very complicated and can be solved only numerically. So, instead of using this approach, we will give only an estimate of plasma heating rate.
Let us consider the domain between the heliopause and the bow
shock in the form of the frustum of a cone with the smaller base
of radius
,
larger base of radius
,
and height x_{m}, where
.
Now we estimate
the rate of the energy increase in this domain due to wave
damping. Up to now we have used the ideal description of the
plasma motion. However, to calculate the rate of wave energy
dissipation we need to take some dissipation mechanism/mechanisms
into account. Since the dissipative coefficients are small,
dissipation occurs only in thin layers that reduce to shocks in
the limit of the ideal description. As a result the rate of wave
energy dissipation does not depend on what particular dissipation
mechanism/mechanisms we take into account. Having this in mind we
choose viscosity. The rate of energy dissipation per unit volume
due to viscosity is
given by