A&A 429, 1069-1080 (2005)
DOI: 10.1051/0004-6361:20041348

Solar cycle influence on the interaction of the solar wind with Local Interstellar Cloud

V. Izmodenov ^{1, 2} - Y. Malama ² - M. S. Ruderman ³

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 $\pm$7 AU, the heliopause variation is $\sim$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

1 Introduction

More than 30 years (three solar cycles) of solar wind observations show that its momentum flux varies by factor of $\sim$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.

2 Model

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 $\rho$ and bulk velocity $\vec{v}$. 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., $n_{\rm e} = n_{\rm p}+n_{\rm He^+}$ for interstellar plasma, and $n_{\rm e}= n_{\rm p}+2 n_{\rm He^{++}}$ 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:

 

$$\frac{\partial \rho }{\partial t}+ \nabla\cdot(\rho \vec{v})=q_1,$$

$$\frac{\partial(\rho \vec{v})}{\partial t}+ \nabla\cdot(\rho \vec{v}\vec{v}+ p\/\/\/\/\mathsf{I}) = \vec{q}_2,$$ (1)

$$\frac{\partial E}{\partial t} + \nabla\cdot[(E+p)\vec{v}] = q_3 .$$

Here $\rho$ is the total density of the ionized component, pis the total pressure of the ionized component, $E=\rho (\varepsilon+ \vec{v}^2/2)$ is the total energy per unit volume, $\varepsilon = p[(\gamma-1)\rho]^{-1}$ is the specific internal energy, and $\mathsf{I}$ is the unit tensor. The temperature of the plasma is determined from the equation of state $p=2(n_{\rm p} +n_{\rm He^+}) k T$ for the interstellar plasma and $p=(2n_{\rm p} +3n_{\rm He^{++}}) k T$ for the solar wind, where k is Boltzman's constant; $n_{\rm p}$, $n_{\rm He^+}$ and  $n_{\rm He^{++}}$ are the proton, interstellar He ion and solar wind alpha particle number densities. In addition to Eqs. (1) we solve the continuity equation for He⁺ in the interstellar medium and for $\alpha$-particles in the solar wind. Then the proton number density is calculated as $n_{\rm p} = (\rho-m_{\rm He}n_{\rm He})/m_{\rm p}$, where $n_{\rm He}$ denotes the He⁺ number density in the interstellar medium, and the He⁺⁺ number density in the solar wind.

The right hand parts, q₁, $\vec{q}_2$ and q₃, are sources of mass, momentum and energy due to charge exchange of H atoms and protons, photoionization and electron impact ionization (Baranov & Malama 1993):

 

$$m_{\rm H} n_{\rm H} (\nu_{\rm ph}+\nu_{\rm imp}),$$                               q₁ =

$$\vec{q}_{2} m_{\rm H} \int\!\!\int u \sigma^{\rm HP}_{\rm ex}(u) (\vec{w}_{\rm H} - \vec{w}_{\rm p}) f_{\rm H}(\vec{w}_{\rm H})$$ =

$$\times f_{\rm p}(\vec{w}_{\rm p})~ {\rm d}\vec{w}_{\rm H}~ {\rm d} \vec{w}_{\rm p}$$

$$+ m_{\rm H} \int (\nu_{\rm ph}+\nu_{\rm imp}) \vec{w}_{\rm H} f_{\rm H}(\vec{w}_{\rm H})~ {\rm d}\vec{w}_{\rm H},$$ (2)

$$m_{\rm H} \int \left( \nu_{\rm ph} + \nu_{\rm imp} \right) \frac{\vec{w}_{\rm H}^2}{2} f_{\rm H}(\vec{w}_{\rm H})~ {\rm d}\vec{w}_{\rm H}$$ q₃ =

$$+ m_{\rm H} \int\!\! \int u \sigma^{\rm HP}_{\rm ex}(u) \frac{\vec{w}_{\rm H}^2- \vec{w}_{\rm p}^2}{2}$$

$$\times f_{\rm H}(\vec{w}_{\rm H}) f_{\rm p}(\vec{w}_{\rm p})~ {\rm d}\vec{w}_{\rm p}~ {\rm d}\vec{w}_{\rm H} .$$

Here, $u=\vert\vec{w}_{\rm p} - \vec{w}_{\rm H}\vert$ is the relative atom-proton velocity, $\sigma^{\rm HP}_{\rm ex}(u)$ is the charge exchange cross section of an H atom with a proton; $\nu_{\rm ph}$ is the photoionization rate; $m_{\rm H}$ is the atomic mass; $\nu_{\rm imp}$ is the electron impact ionization rate. Numerical values for $\sigma^{\rm HP}_{\rm ex}$, $\nu_{\rm ph}$ and $\nu_{\rm imp}$are given in Izmodenov et al. (1999) and Izmodenov et al. (2004). The function $f_{\rm p}(\vec{w}_{\rm p}) = n_{\rm p} (\!\sqrt{\pi}c_{\rm p})^{-3} \exp(-(\vec{w}_{\rm p} - \vec{v})^2/c_{\rm p}^2)$ is the local Maxwellian velocity distribution function of protons, with the parameters ($n_{\rm p}$, $\vec{v}$, $c_{\rm p} = (2 k T/m_{\rm H})^{1/2}$) determined by the solution of the Euler Eqs. (1). The function $f_{\rm H}(\vec{w}_{\rm H})$is the velocity distribution function of H atoms.

The velocity distribution of H atoms $f_{\rm H}(\vec{r}, \vec{w}_{\rm H}, t)$ is calculated from the linear kinetic equation:

 

$$\frac{\partial f_{\rm H}}{\partial t} + \vec{w}_{\rm H} \cdot \fr... ...{\vec{F}}{m_{\rm H}} \cdot \frac{\partial f_{\rm H}}{\partial \vec{w}_{\rm H} } -f_{\rm H} \int \vert\vec{w}_{\rm H} - \vec{w}_{\rm p} \vert \sigma^{\rm HP}_{\rm ex} f_{\rm p} (\vec{r}, \vec{w}_{\rm p})~ {\rm d}\vec{w}_{\rm p}$$ = (3)

$$+ f_{\rm p} (\vec{r}, \vec{w}_{\rm H}) \int \vert \vec{w}_{\rm H}... ... HP}_{\rm ex} f_{\rm H} (\vec{r}, \vec{w}_{\rm H}^* )~ {\rm d}\vec{w}_{\rm H}^*$$

$$- ( \nu_{\rm ph} + \nu_{\rm imp} ) f_{\rm H} ( \vec{r}, \vec{w}_{\rm H} ) .$$

Here $\vec{w}_{\rm p}$ and $\vec{w}_{\rm H}$ are the individual proton and H atom velocities, respectively, and $\vec{F}$ is the sum of the solar gravitational force and the solar radiation pressure force. The plasma and neutral components interact mainly by charge exchange. Photoionization, solar gravity and radiation pressure taken into account in Eq. (3) are only important at small heliocentric distances. Electron impact ionization may be important in the inner heliosheath, which is the region between the termination shock and the heliopause. Note that photoionization rate and solar radiation pressure may vary with the solar cycle. At present we do not include solar cycle variations of these parameters in our model, but we plan to do so in the future.

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:

 

$$n_{{\rm p}E} n_{{\rm p}E,0}(1+\delta_n\sin\omega t),$$ =

$$v_{{\rm p}E,0},$$ v_E = (4)

$$T_{E0} = \mbox{const}.$$ T_E =

For the solar wind disturbances determined by Eq. (4), the ratio of the maximum to minimum momentum flux is equal to $\Delta=(1+\delta_n)/(1-\delta_n)$. Following Baranov & Zaitsev (1998) we made calculations for $\delta_n = 1/3$, so that $\Delta = 2$. As it was pointed out by Zank & Mueller (2003), the solar cycle effects in the heliospheric interface remain the same when the variation of the solar wind dynamic pressure is caused by the solar wind velocity variation. In this paper we consider variations of the solar wind density which is sufficient for the purposes of this paper. The effects of the realistic solar cycle on the TS variation were studied by Izmodenov et al. (2003a).

The following solar wind parameters averaged over a few solar cycles were used: $n_{{\rm p}E,0} = 8$ cm^-3, $V_{{\rm p}E,0}= 445$ km s^-1. We also used the following parameters of the interstellar gas in the unperturbed interstellar medium at the outer boundary: $V_{\rm LIC} = 26.4$ km s^-1, T $_{\rm LIC} = 6500$ K, $n_{\rm H,LIC} = 0.18$ cm^-3, $n_{\rm p,LIC} = 0.06$ 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 $n_{\rm H,LIC}$ and $n_{\rm p,LIC}$ 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₁, $\vec{q}_2$ and q₃ 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 $t_{\rm period}$ = 66 years. To minimize statistical errors we averaged the statistical results over $t_{\rm mc}= 1$ 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₁, $\vec{q}_2$ and q₃ 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, $\rho _E V^2_E$, 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 ${\sim} 30$ 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.

3 Results

3.1 Plasma component

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 $\sim$93 AU, which is reached in the last (11th) year of the "ideal'' solar cycle, to its maximum distance of $\sim$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 $\sim$25 AU from its minimum value of ${\sim} 163$ AU during the third year of the cycle to its maximum value of $\sim$188 AU in the 9th year of the cycle. In the upwind direction the most distant position of the termination shock is reached $\sim$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 ${\sim} 3.5$ 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 $\beta$that determines the variation of intensity of the cosmic rays jwith energy E: $j \sim E^{\beta}$. 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 $\beta$ 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 $\beta$ 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 $\sim$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 $\sim$308 AU, while this distance is $\sim$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 ${\sim} 10$ 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: t1 = 1 year (curves 1), t2= 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 year (curves 1) and t₂ = 6 year (curves 2). It is seen that the characteristic wavelength in the region is ${\sim} 40$ 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 ${\sim} 200$ 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.

3.2 H atoms

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 year and dashed curves to t₂ = 6 year. It is seen that the variation of the density is within $\pm$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 $\pm$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 $\sim$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 ${\rm Kn} = l/L \ll 1$, 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 ${\sim} 50$ AU. Therefore, $\rm Kn_{stationary} \approx 0.5$. 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 $L \approx 20$AU as follows from Fig. 3. Therefore, $\rm Kn_{time}\approx 2$ 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.

4 Summary

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:

1.

The solar cycle variation of the TS location is $\pm$7 AU about its mean value.

2.

The heliopause varies by $\pm 2$ AU about its mean value.

3.

The variation of the bow shock location is negligible.

4.

There is a sequence of additional weak shocks and rarefaction waves in the region between the heliopause and the bow shock. The additional heat of the plasma in the outer heliosheath induced by the shock waves is small and it is not observable in our calculations.

5.

Our numerical results in the region between the HP and BS are confirmed by an analytical solution based on the WKB approximation.

For the interstellar H atom component we obtain the following new results:

1.

Variation of the number density of the H atoms in the outer heliosphere is within 10%. The variation increases at 5 AU up to 30% due to strong ionization processes in the vicinity of the Sun.

2.

The variations of the number densities of three populations of H atoms - primary and secondary interstellar atoms, and atoms created in the inner heliosheath - are coherent in the entire supersonic solar wind region and determined by loss due to charge exchange. The coherent behavior of fluctuations disappear in the regions where the production process is dominant.

3.

There is no significant variation of the temperature and bulk velocity of the primary and secondary interstellar H atoms with the solar cycle. However, the bulk velocity and kinetic temperature of atoms created in the inner heliosheath vary with the solar cycle by 10-12%. It is shown that this variation reflects the plasma properties at the heliopause.

4.

There is a qualitative difference between our results and the results obtained by using the fluid or multi-fluid description for the interstellar H atoms. It was shown that the multi-fluid description is less appropriate for the time-dependent case than for the stationary case because the Knudsen number is larger for the time-dependent problem.

In addition, we have performed specific calculation of the time-dependent Euler equation with the source terms (3) taken from the corresponding stationary solution. The difference in plasma distribution with the self-consistent model is a few percent (Fig. 8). Therefore, we suggest that time-dependent multi-fluid models may produce results that are closer to the kinetic time-dependent model in the case when the source terms are taken from the stationary solution.

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)''.

References

• Alexashov, D. B., Baranov, V. B., Barsky, E. V., & Myasnikov, A. V. 2000, Astron. Lett., 26, 743 [NASA ADS] [CrossRef] (In the text)
• Alexashov, D. B., Chalov, S. V., Myasnikov, A. V., Izmodenov, V. V., & Kallenbach, R. 2004a, A&A, 420, 729 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Alexashov, D. B., Izmodenov, V. D. B., & Grzedzielski, S. 2004b, Adv. Spca Res., 34, 1, 109
• Baranov, V. B., & Malama, Yu. G. 1993, J. Geophys. Res., 98, 15157 [NASA ADS] (In the text)
• Baranov, V. B., & Malama, Yu. G. 1996, Space Sci. Rev., 78, 305 [NASA ADS] [CrossRef]
• Baranov, V. B., & Zaitsev, N. A. 1998, Geophys. Res. Lett., 25, 4051 [NASA ADS] [CrossRef] (In the text)
• Baranov, V. B., Lebedev, M. G., & Malama, Yu. G. 1991, ApJ, 375, 347 [NASA ADS] [CrossRef] (In the text)
• Baranov, V. B., Izmodenov, V., & Malama, Yu. G. 1998, J. Geophys. Res., 103, A5, 9575 (In the text)
• Gazis, P. R. 1996, Rev. Geophys., 34, 379 [NASA ADS] [CrossRef] (In the text)
• Gloeckler, G., Moebius, E., & Bzowski, M. 2004, A&A, in press (In the text)
• Gruntman, M., Roelof, E., Mitchell, D., et al. 2001, J. Geophys. Res., 106, 15767 [NASA ADS] [CrossRef] (In the text)
• Izmodenov, V. V. 2000, Astrophys. Space Sci., 274, issue 1/2, 55 (In the text)
• Izmodenov, V. V. 2001, Space Sci. Rev., 97(1/4), 385 (In the text)
• Izmodenov, V. V. 2004, in The Sun and the Heliosphere as an Integrated System, ed. G. Poletto, & S. T. Suess (Kluwer), in press (In the text)
• Izmodenov, V. V. 2003, Proceedings of the Interstellar enviroment of the heliosphere, COSPAR Colloqium in Honour of Stanislaw Grzedzielski, ed. D. Breitschwerdt, & G. Haerendel, MPE Report, 285, 113 (In the text)
• Izmodenov, V. V., & Alexashov, D. 2003, Astron. Lett., 29, 58 [NASA ADS] [CrossRef] (In the text)
• Izmodenov, V. V., & Malama, Y. G. 2004a, Adv. Space Res., 34(1), 74 (In the text)
• Izmodenov, V. V., & Malama, Y. G. 2004b, Proc. 3rd IGPP Conf. on Astrophysics - Physics of the Outer Heliosphere, Riverside, California, 8-13 February 2004, AIP, Conf. Proc., 719, 47
• Izmodenov, V., Gruntman, M., & Malama, Yu. G. 2001, J. Geophys. Res., 106, 10681 [NASA ADS] [CrossRef] (In the text)
• Izmodenov, V. V., Geiss, J., Lallement, R., et al. 1999, J. Geophys. Res., 104, 4731 [NASA ADS] [CrossRef] (In the text)
• Izmodenov, V., Gloeckler, G., & Malama, Yu. G. 2003a, Geophys. Res. Lett., 30, 1351 [NASA ADS] (In the text)
• Izmodenov, V. V., Malama, Yu. G., Gloeckler, G., & Geiss, J. 2003b, ApJ, 954, L59 [NASA ADS]
• Izmodenov, V., Malama, Yu. G., Gloeckler, G., & Geiss, J. 2004, A&A, 414, L29 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Karmesin, S. R., Liewer, P. C., & Brackbill, J. U. 1995, Geophys. Res. Lett., 22, 1153 [NASA ADS] [CrossRef] (In the text)
• Malama, Yu. G. 1991, Astrophys. Space Sci., 176, 21 [NASA ADS] [CrossRef] (In the text)
• Myasnikov, A., Alexashov, D., Izmodenov, V., & Chalov, S. 2000, J. Geophys. Res., 105, 5167 [NASA ADS] [CrossRef] (In the text)
• Narain, U., & Ulmschneider, P. 1990, Space. Sci. Rev., 54, 377 [NASA ADS] (In the text)
• Priest, E. R. 1982, Solar Magnetohydrodynamics (Dordrecht: Reidel)
• Richardson, J. D. 1997, Geophys. Res. Lett., 24, 2889 [NASA ADS] [CrossRef] (In the text)
• Rudenko, O. V., & Soluyan, S. I. 1977, Theoretical Foundations of Nonlinear Acoustics (New York, London: Plenum Publishing Corporation, Consultant Bureau) (In the text)
• Scherer, K., & Fahr, H. J. 2003a, Geophys. Res. Lett., 30 (In the text)
• Scherer, K., & Fahr, H. J. 2003b, Annales Geophys., 21, 1303
• Wang, C., & Belcher, J. W. 1998, J. Geophys. Res., 103, 247 [NASA ADS] [CrossRef] (In the text)
• Wang, C., & Belcher, J. W. 1999, J. Geophys. Res., 104, 549 [NASA ADS] [CrossRef] (In the text)
• Witte, M., Banaszkiewicz, M., & Rosenbauer, H. 1996, Space Sci. Rev., 78, 289 [NASA ADS] [CrossRef] (In the text)
• Witte, M. 2004, A&A, 426, 835 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Zaitsev, N. A., & Izmodenov, V. V. 2001, in The Outer Heliosphere: The Next Frontiers, ed. K. Scherer, H. Fichtner, H. J. Fahr, & E. Marsch, 65 (In the text)
• Zank, G. P. 1999, in Solar Wind 9, ed. S. Habbal, R. Esser, J. Hollweg, & P. Isenberg, AIP, 783 (In the text)
• Zank, G., & Mueller, H. 2003, J. Geophys. Res., 108, A6, 1240 (In the text)

 

  
Online Material

Appendix A

In this appendix we study the propagation of perturbations caused by the solar cycle variation in the outer heliosheath. We restrict our analysis to the vicinity of the symmetry axis. In this vicinity the heliopause can be approximated by a spherical surface with the radius R₀ equal to the curvature radius of the heliopause at the symmetry axis, and with center on the symmetry axis at a distance R₀ from the heliopause in the solar direction. We introduce spherical coordinates with center at the center of spherical surface and the radial coordinate R = R₀ + x (see Fig. A.1). In what follows we neglect the dependence of the background quantities on the angle coordinates and assume that they are functions of x only. In addition, we neglect the velocity components in the angular directions. Hence, in the unperturbed state, $\rho = \rho_0(x)$, u = u₀(x) and p = p₀(x), where u is the velocity in the radial direction, and the two other components of the velocity are zero.

Figure A.1: The sketch of spherical coordinate system. The heliopause and bow shock are shown by the solid lines. The spherical surface approximating the heliopause in the vicinity of the symmetry axis is shown by the dashed line.

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

$$u = \epsilon\hat{u}\sin~(\omega t)$$ (A.1)

at the heliopause (x = 0). Here $\hat{u}$ is a characteristic speed near the heliopause and $\epsilon$ is a dimensionless parameter, $\epsilon \ll 1$.

An important quantity for our analysis is the sound speed $c_S = (\gamma p_0/\rho_0)^{1/2}$. In what follows we take $\hat{u} = c_S(0)$.

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  $\approx 3.5\times 10^8$ s. This gives the typical wavelength of about $7\times 10^9\;\mbox{km} \approx 50$ 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 $\epsilon$. In other words, we assume that the ratio of the characteristic scales of perturbation and inhomogeneity is $\epsilon$. To show this explicitly we introduce the stretching variable $\sigma = \epsilon x$ and assume that $\rho_0$, u₀ and p₀ are functions of $\sigma$. In line with the WKB method we also introduce the running variable $\tau = t - \epsilon^{-1}\varphi(\sigma)$, where $\varphi(\sigma)$ is the eiconal. In the new variables we rewrite the system of governing equations as

$$\frac{\partial\rho}{\partial\tau} - \frac{{\rm d}\varphi}{{\... ...n}{\zeta^2}\frac{\partial(u\rho\zeta^2)}{\partial\sigma} =q_1,$$ (A.2)

$$\frac{\partial u}{\partial\tau} - u\frac{{\rm d}\varphi}{{\r... ...\frac{\epsilon}{\rho}\frac{\partial p}{\partial\sigma} = q_4 ,$$ (A.3)

$$\frac{\partial}{\partial\tau}\left(\frac{p}{\rho^\gamma}\righ... ...l}{\partial\sigma} \left(\frac{p}{\rho^\gamma}\right) = q_5 ,$$ (A.4)

where $\zeta = \epsilon R = \epsilon R_0 + \sigma$, and q₄ and q₅ can be expressed in terms of q₁, $\vec{q}_2$ and q₃. Recall that the right-hand sides of these equations describe the interaction between the plasma and the neutral gas. In what follows we neglect the effect of the neutral gas on the perturbations and assume that q₁, q₄ and q₅ do not change when the perturbations are introduced. We are looking for the solution to the system of Eqs. (A.2)-(A.4) in the form of expansions in power series with respect to $\epsilon$,

$$f = f_0 + \epsilon f_1 + \epsilon^2 f_2 + \dots ,$$ (A.5)

where f represents any of the quantities $\rho$, u and p.

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 $\rho_0$, u₀ and p₀. The solution of these equations is shown in Fig. A.2, where the background density, $\rho_0$, velocity, u₀, and temperature, T₀, are shown as functions of x. The background pressure p₀ is calculated using the relation $p_0 = (2k_{\rm B}/m_{\rm p})\rho_0 T_0$, where $k_{\rm B}$ is the Boltzmann constant and $m_{\rm p}$ the proton mass.

Figure A.2: The dependences of background quantities in the outer heliosheath on the distance x from the heliopause at the symmetry axis. The upper, middle and lower panels correspond to the density, velocity and temperature. The density is given in $\rho _{\rm p,\infty }$, the velocity in $v_{\rm p,\infty }$, and the temperature in 104 K. The distance x from the heliopause is given in AU.

In the first order approximation we collect terms of order $\epsilon$ in Eqs. (A.2)-(A.4). Since we have assumed that perturbations do not affect q₁, q₄ and q₅, 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 $\partial\rho_1/\partial\tau$, $\partial u_1/\partial\tau$ and $\partial p_1/\partial\tau$. It has a non-trivial solution only if its determinant is zero. This condition results in the equation determining the eiconal

$$\left(1 - u_0\frac{{\rm d}\varphi}{{\rm d}\sigma}\right)^2 = c_S^2\left(\frac{{\rm d}\varphi}{{\rm d}\sigma}\right)^2 \cdot$$ (A.6)

In what follows we only consider waves propagating in the positive x-direction. Then it follows from Eq. (A.6) that the eiconal is given by

$$\frac{{\rm d}\varphi}{{\rm d}\sigma} = \frac1{c_S + u_0} \cdot$$ (A.7)

In addition, it follows that $\rho_1$ and p₁ can be expressed in terms of u₁ as

$$\rho_1 = \frac{\rho_0}{c_S}~u_1 , \quad p_1 = \rho_0 c_S u_1 .$$ (A.8)

In the second order approximation we collect terms of order  $\epsilon^2$ in Eqs. (A.2)-(A.4). This results in the system of equations that, with the aid of Eqs. (A.7) and (A.8), can be written as

 

$$c_S\frac{\partial\rho_2}{\partial\tau} - \rho_0\frac{\partial u_2... ...rtial}{\partial\sigma}~ \left[\frac{\rho_0\zeta^2(c_S + u_0)}{c_S}~u_1\right] ,$$ (A.9)

 

$$c_S\frac{\partial u_2}{\partial\tau} - \frac1{\rho_0}\frac{\parti... ...\rho_0 c_S)}{\partial\sigma} - \frac{\partial u_0}{\partial\sigma}\right] u_1 ,$$ (A.10)

 

$$\frac{\partial p_2}{\partial\tau} - c_S^2\frac{\partial\rho_2}{\p... ..._S}~ \frac{\partial}{\partial\sigma}\left(\frac{p_0}{\rho_0^\gamma}\right)\cdot$$ (A.11)

This is a linear system of inhomogeneous algebraic equations with respect to $\partial\rho_2/\partial\tau$, $\partial u_2/\partial\tau$ and $\partial p_2/\partial\tau$. The determinant of this system is zero. This means that it has a non-trivial solution only if its right-hand side satisfies the compatibility condition. To obtain this condition we eliminate all the variables of the second order approximation from Eqs. (A.9)-(A.11). As a result we arrive at

$$\frac{\partial u_1}{\partial\sigma} - \lambda u_1\frac{\partial u_1}{\partial\tau} + \mu u_1 = 0,$$ (A.12)

where the coefficient functions $\lambda(\sigma)$ and $\mu(\sigma)$ are given by

$$\lambda = \frac{\gamma + 1}{2(c_S + u_0)^2} ,$$ (A.13)

$$\mu = \frac{c_S}{2\rho_0\zeta^2(c_S + u_0)^2} \frac{\partial... ...ma}\left[ \frac{\rho_0\zeta^2(c_S + u_0)^2}{c_S}\right] \cdot$$ (A.14)

When the nonlinear term in Eq. (A.12) is neglected, it reduces to the equation of conservation of the wave action flux in a radially expanding tube with the cross-section proportional to $\zeta^2$, which is $\rho_0\zeta^2(c_S + u_0)^2 u_1^2/c_S = \mbox{const}$.

Let us return to the original independent variables, use the approximation $u' \equiv u-u_0 \approx \epsilon u_1$, and make the variable substitution

 

u' = Us(0)/s(x) , (A.15)

where

$$u' = Us(0)/s(x) , \quad s = R\rho_0^{1/2}c_S^{-1/2}(c_S + u_0) .$$ (A.16)

Then Eq. (A.12) transforms to

$$\frac{\partial U}{\partial x} + \frac1{c_S + u_0}\frac{\part... ...U}{\partial t} - \Lambda U\frac{\partial U}{\partial t} = 0 ,$$ (A.17)

where $\Lambda = \lambda s(0)/s(x)$. In order to solve Eq. (A.17) we consider x and U as independent variables, and t as a dependent variable, i.e. t = t(x,U). Then Eq. (A.17) transforms to

$$\frac{\partial t}{\partial x} = \frac1{c_S + u_0} - \Lambda U .$$ (A.18)

Integrating this equation we obtain

$$t = \varphi(x) - U\psi(x) + F(U) ,$$ (A.19)

where

$$\varphi(x) = \int_0^x \frac{{\rm d}\tilde{x}}{c_S(\tilde{x}) ... ...\quad \psi(x) = \int_0^x \Lambda(\tilde{x})~{\rm d}\tilde{x} ,$$ (A.20)

and the function F(U) is determined by the boundary condition (A.1). The right-hand side of Eq. (A.1) is a non-monotonic function of t. This implies that F(U) is a multi-valued function. The periodic driving described by Eq. (A.1) causes periodic motion in the region x > 0. This observation enables us to consider only one wave period corresponding to $-\pi \leq \omega t \leq \pi$ at x = 0. Then F(U) is given by

$$F(U) = \left\{\begin{array}{ll} \displaystyle \frac1\omega\l... ...antom{-}\frac\pi2 \leq \omega t \leq \pi . \end{array} \right.$$ (A.21)

The inequalities for $\omega t$ have to be satisfied at x = 0.

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 $\tau = t - \varphi(x)$for three different values of x. The upper panel corresponds to x = 0, and U is a sinusoidal function of $\tau = t$. 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. $U/\epsilon\hat{u}$ is shown as a function of $\omega \tau = \omega [t - \varphi (x)]$. The upper, middle and lower panes correspond to x = 0, $x = x_{\rm c} \approx 44.5$ AU, and x = 70 AU. In the lower panel the dashed line shows the unphysical branch of the multivalues function $U(\tau )$determined by Eq. (A.19). The vertical solid line is the shock.

In order to determine the spatial position $x_{\rm c}$ and the point of the wave profile where the gradient catastrophe occurs, we note that at this point $\partial\tau/\partial U = \partial^2\tau/\partial U^2 = 0$. Using Eqs. (A.19) and (A.21) we rewrite these conditions as

$$\omega\hat{u}\psi(x_{\rm c}) = \left(\epsilon^2 - U^2\right)^{-1/2} , \quad \omega U\left(\epsilon^2 - U^2\right)^{-3/2} = 0 .$$ (A.22)

It follows from these equations that the gradient catastrophe occurs at the point of the wave profile where U = 0, and at the spatial position $x_{\rm c}$ determined by the equation

$$\psi(x_{\rm c}) = \frac1{\epsilon\omega\hat{u}} \cdot$$ (A.23)

If we formally use Eq. (A.19) for $x > x_{\rm c}$, then we obtain U as a multi-valued function of $\tau$ as shown in the lower panel of Fig. A.3. Such a multi-valued solution does not make physical sense. To obtain a physically meaningful solution we have to allow a discontinuity at a particular position of the wave profile. This discontinuity corresponds to a shock. It follows from the Rankine-Hugoniot relations at shocks that the velocity perturbation at the two sides of the shock have the same magnitude and opposite signs (see, e.g., Rudenko & Soluyan 1977). It immediately follows from this condition and Eq. (A.19) that the shock position in the wave profile is $\tau = 0$, and the shock intensity $U_{\rm s}$ is determined by the equation and inequality

$$U_{\rm s} = \epsilon\hat{u}\sin~(\omega U_{\rm s}\psi) , \quad 0 < \omega U_{\rm s}\psi \leq \pi/2 .$$ (A.24)

We have used Eqs. (A.23) and (A.24) to calculate $x_{\rm c}$ and the dependence of the velocity jump at the shock, $\Delta u = 2U_{\rm s} s(0)/s(x)$. Using the numerically obtained value for the amplitude of the heliopause oscillation, 2 AU, and the fact that the oscillation period is equal to the solar cycle period, we immediately obtain that $\epsilon\hat{u} \approx 5.4$ km s^-1. Since $\hat{u} = c_{\rm S}(0) \approx 40$ km s^-1, we obtain $\epsilon \approx 0.135$, so that the use of $\epsilon$ as a small parameter is justified. Using $\epsilon\hat{u} \approx 5.4$ km s^-1and numerically calculated functions $\rho_0(x)$, u₀(x) and $c_{\rm S}(x)$, we found that $x_{\rm c} \approx 44.5$ AU. Since the distance between the heliopause and the bow shock is about 140 AU, we conclude that the gradient catastrophe occurs inside the outer heliosheath.

In Fig. A.4 the dependence of $\Delta u$ on x is shown for $x_{\rm c} \leq x \leq x_m$, where $x_m \approx 140$ 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.

Figure A.5 displays snapshots of the wave profile at three different moments of time. We can see from this figure that the wave amplitude at the bow shock is approximately equal to $\Delta u/2 \approx 2$ km s^-1. If we neglect nonlinearity and use conservation of the wave action, then we obtain that the wave amplitude at the bow shock would be about 10.6 km s^-1, which is almost two times larger than $\epsilon\hat{u}$. Hence, we conclude that the wave amplitude decreases by more than five times due to dissipation at the shock. Since the wave energy flux is proportional to the wave amplitude, it follows that dissipation at the shock reduces the wave energy flux by one and a half orders of magnitude in comparison with the value given by the linear theory.

Figure A.5: The snapshots of the wave profile at $\omega t = n\pi$, $\omega t = (n+1/2)\pi$ and $\omega t = (n+1)\pi$, where n is any integer number. The dotted vertical lines indicate the bow shock position. The velocity is given in km s-1 and the distance 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 $\epsilon R_0$, larger base of radius $\epsilon(R_0 + x+m)$, and height x_m, where $\epsilon \ll 1$. 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

$$H_v = \frac{4\rho\nu}3\left(\frac{\partial u'}{\partial x}\right)^2 ,$$ (A.25)

where $\nu$ is the kinematic viscosity. For waves with small amplitude we can substitute $\rho_0$ for $\rho$ in this expression. Since viscosity is very weak, it is only important in the region of large gradients, i.e. in the vicinity of the shock position in the ideal description. This observation enables us to use the approximate description of the wave profile (see, e.g., Rudenko & Soluyan 1977). For $x < x_{\rm c}$ we use the ideal description. For $x > x_{\rm c}$ we use the ideal description for all moments of time except a small time interval embracing the moment of time when the shock passes this particular spatial position. This interval is determined by the inequality $\vert\omega\tau\vert \lesssim \delta = \nu\omega/\hat{u}^2 \ll 1$ (recall that $\tau = t - \varphi(x)$). The solution of viscous hydrodynamic equations inside this interval describes the structure of the shock. It is given by (e.g., Rudenko & Soluyan 1977)

$$u' = u'_v \equiv \frac{\Delta u}2\tanh\frac{\kappa(x)\omega\tau}{\delta} ,$$ (A.26)

where $\kappa(x) = (\gamma+1)c_S\Delta u/4\hat{u}^2$. Using the relation $\partial u'_v/\partial x \approx (1/c_S)\partial u'_v/\partial\tau$ and Eq. (A.26), and taking into account that dissipation occurs only in the thin layer $\vert\omega\tau\vert \lesssim \delta$, we obtain that the rate of energy dissipation per unit volume at the position x > x₀averaged over the wave period is

 

$$\langle H_v \rangle \approx \frac{4\rho_0\nu}{3c_S^2(2\pi/\omega)} \int_{-\delta/\omega}^{\delta/\omega} \left(\frac{\partial u'_v}{\partial\tau}\right)^2{\rm d}\tau$$ (A.27)

$$\approx \frac{(\gamma+1)\omega\rho_0(\Delta u)^3}{24\pi c_S} \int_{-\inft... ...d}\xi}{\cosh^4\xi} = \frac{(\gamma+1)\omega\rho_0(\Delta u)^3}{18\pi c_S} \cdot$$

When deriving this expression we have used the fact that $1/\cosh\xi$ decreases very fast when $\vert\xi\vert\to\infty$ and substituted the lower and upper limits of integration by $-\infty$and $\infty$ respectively. To obtaine the rate of energy dissipation in volume V averaged over the period, L_v, we have to integrate $\langle H_v \rangle$ over V. Taking into account that $\langle H_v \rangle = 0$ for $x < x_{\rm c}$ and using the numerical results, we obtain

 

$$L_v = \pi\epsilon^2\int_{x_{\rm c}}^{x_m}\langle H_v \rangle (R_0... ...0^{38}\epsilon^2\omega\rho_{{\rm p},\infty}v_{{\rm p},\infty}^2\; ~\mbox{J/s} ,$$ (A.28)

where $\omega$, $\rho _{\rm p,\infty }$ and $v_{\rm p,\infty }$ has to be measured in the SI units, i.e. in s^-1, kg/m³ and m/s respectively. The total mass of the plasma in the considered domain is $M \approx 10^{41}\pi\epsilon^2\rho_{{\rm p},\infty}\; \mbox{kg}$. We take c_v M as an estimate for the thermal capacity of the plasma in the domain, where c_v is specific heat at constant density. Since $c_v = 2k_{\rm B}[(\gamma-1)m_{\rm p}]^{-1} \approx 2.5\times 10^4\;~ \mbox{J}~\mbox{K}^{-1}~\mbox{kg}^{-1}$, we obtain $c_v M \approx 2.5\times 10^{45}\pi\epsilon^2\rho_{\rm p,\infty}\;~\mbox{J/K}$. Then, if we neglect all energy losses from the domain, then its mean temperature will increase by

$$(2\pi/\omega)L_v/c_v M \approx 2.2\times10^{-7}v_{\rm p,\infty}^2\;\mbox{K} \approx 150\;\mbox{K}.$$ (A.29)

This quantity gives us only an upper boundary of the temperature increase. In reality there is thermal energy loss due to interaction with neutral hydrogen. Hence, the real temperature increase will be smaller than 150 K per one solar cycle. We can also expect that the temperature growth will saturate after a few cycles.