EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 517 (July 2010) A&A, 517 (2010) A35 Full HTML
Free Access
Issue
A&A
Volume 517, July 2010
Article Number A35
Number of page(s) 16
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/200912376
Published online 29 July 2010
A&A 517, A35 (2010)

Numerical simulation of viscous-like flow in and around the plasma tail of a comet

M. Reyes-Ruiz1 - H. Pérez-de-Tejada2 - H. Aceves1 - R. Vázquez1

1 - Instituto de Astronomía, Universidad Nacional Autónoma de Mexico, Apdo Postal 877, Ensenada 22800, B.C., Mexico
2 - Instituto de Geofísica, Universidad Nacional Autónoma de Mexico, Ciudad Universitaria, C.P. 04510, Mexico D.F.

Received 22 April 2009 / Accepted 2 March 2010

Abstract
Aims. We model the interaction of the solar wind with the plasma tail of a comet using numerical simulations, taking into account the effects of viscous-like forces.
Methods. We developed a 2D hydrodynamical, two species, finite difference code to solve the time-dependent continuity, momentum, and energy conservation equations, and model the interaction of the solar wind with a cometary plasma tail. We compute the evolution of the plasma of cometary origin in the tail as well as the properties of the shocked solar wind plasma around it, as it transfers momentum on its passage by the tail. Velocity, density and temperature profiles across the tail are obtained. Several models with different flow parameters are considered to study the relative importance of viscous-like effects and the coupling between species on the flow dynamics.
Results. Assuming a Mach number equal to 2 for the incident solar wind as it flows past the comet's nucleus, the flow exhibits three transitions with location and properties depending on the Reynolds number of each species and on the ratio of the timescale for inter-species coupling to the crossing time of the free-flowing solar wind. By comparing our results with the measurements taken $in \ situ$ by the Giotto spacecraft during its flyby of comet Halley, we constrain the flow parameters for both plasmas.
Conclusions. In the context of our approximations, we find that our model is qualitatively consistent with the $in \ situ$ measurements as long as the Reynolds number of both solar wind protons and cometary H2O+ ions is low, less than 100, suggesting that viscous-like momentum transport processes may play an important role in the interaction of the solar wind and the plasma environment of comets.

Key words: comets: general - solar wind - hydrodynamics - comets: individual: Halley - methods: numerical - comets: individual: Giacobinni-Zinner

1 Introduction

The interaction of the solar wind with the plasma environment of comets as they approach the Sun, has been investigated since the early days of space physics (see Biermann 1951; Alfvén 1957; reviews by Cravens & Gombosi 2004, and Ip 2004). The basic elements of the interaction were developed following the work of Biermann (1951), who proposed that the interaction between the solar wind and the comet's plasma is responsible for the observed aberration angle of plasma tails with respect to the Sun-comet radius vector. Due to the inefficiency of Coulomb collisional processes in the coupling of the solar wind with cometary plasmas, Alfven (1957) proposed that the interplanetary magnetic field (IMF) is a fundamental ingredient in the interaction between the solar wind and comets, responsible for channelling the cometary ions as it drapes into a magnetic tail. Biermann et al. (1967) suggested that as cometary ions are created and incorporated (picked-up) into the solar wind, the loading of the flow with this additional mass results in a modification of the flow properties as the solar wind approaches a comet, an idea further developed by Wallis (1973) (for a review see Szego et al. 2000). The IMF and mass loading are thus the main dynamical agents generally considered when developing models for the interaction of the solar wind with cometary ionospheres, as well as with other solar system bodies with an ionosphere and without a strong intrinsic magnetic field.

However, as pointed out by Perez-de-Tejada et al. (1980) and Perez-de-Tejada (1989), several features of the flow dynamics in the cometosheath and plasma tail of comets can be attributed to the action of viscous-like forces as the solar wind interacts with cometary plasma. These interaction processes are believed to be similar to those known occur in other solar system bodies that have an ionosphere and no significant intrinsic magnetic field, particularly Venus and Mars (for a review see Perez-de-Tejada 1995; Perez-de-Tejada et al. 2009 and references therein). In situ measurements indicate that, as in Venus and Mars, the solar wind flow in the ionosheath of comet Halley exhibits an intermediate transition, also called the ``mystery transition'', located approximately half-way between the bow shock and the cometopause (Johnstone et al. 1986; Goldstein et al. 1986; Reme 1991; Perez-de-Tejada 1989, and references therein). Below this transition, as we approach the cometopause, the antisunward velocity of the shocked solar wind decreases in a manner consistent with a viscous boundary layer (Perez-de-Tejada 1989). The presence of viscous-like processes is also supported by the increase in the gas temperature, and the decrease in density, as we move from the intermediate transition to the cometopause. Taking the distance between the intermediate transition and the cometopause to represent the thickness of a viscous boundary layer, which depends on the effective Reynolds number of the flow ( $R_{\rm eff}$), Perez-de-Tejada (1989) estimated that $R_{\rm eff} \approx 300$ for the solar wind flow in the cometosheath is necessary to reproduce the flow properties measured in situ by the Giotto spacecraft on its flyby of comet Halley.

Additional support for the importance of viscous-like effects in the dynamics of the flow in the cometosheath and tail regions, is provided by comparing the magnitude of the terms corresponding to momentum transport caused by viscous-like forces with the ${\bf J} \times {\bf B}$ forces in the momentum conservation equation. Perez-de-Tejada (1999, 2000) argued that downstream of the terminator in the ionosheath of Venus, a scenario analogous to that considered in this paper, the fact that the flow is superalfvenic, as found from the in situ measurements of the Mariner 5 and Venera 10 spacecraft, suggests that viscous-like forces may dominate over ${\bf J} \times {\bf B}$ forces in the flow dynamics in the boundary layer formed in the interaction between the solar wind and ionospheric plasma. If the flow were characterized by a low effective Reynolds number, $R_{\rm eff}$, this layer would extend over a significant portion of the ionosheath of the planet.

For comets, in situ measurements obtained during the passage of the ICE spacecraft through the tail of comet Giacobinni-Zinner (Bame et al. 1986; Slavin et al. 1986; Meyer-Vernet 1986; Reme 1991) indicated that along the inbound trajectory (which lies slightly tailward of the comet nucleus) the magnetic field in the so-called transition and sheath regions, is approximately 10 nT, the number density is approximately 10 cm-3 and the tailward flow velocity varies from $\sim$400 km s-1 (near the bow shock) to 100 km s-1. According to Perez-de-Tejada (1999), the ratio of viscous-like to magnetic forces is essentially the square of the alfvenic Mach number, $M_{\rm A}^2 = V^2 / (B^2/8 \pi \rho)$. From the data of the ICE spacecraft cited above, we find that $M_{\rm A}^2$ranges between 4 and 40 across the cometosheath and hence, viscous-like stresses may dominate over ${\bf J} \times {\bf B}$ forces by a similar amount, or more, throughout the cometosheath region tailward of the nucleus. In the vicinity of the plasma tail, the measurements of the ICE spacecraft (Bame et al. 1986; Slavin et al. 1986) indicate that the midplane density, dominated by cometary ions, reaches values of 200 cm-3 at the point where the magnetic field has its maximum strength of 50 nT. With flow speeds of approximately 20 km s-1, the square of the alfvenic Mach number reaches a minimum value of 2-3 so that, even in the plasma tail, viscous-like forces are, at least, as important as ${\bf J} \times {\bf B}$ forces following the arguments of Perez-de-Tejada (2009).

That $M_{\rm A}^2 \gg 1$ in the cometosheath means that the magnetic energy density is much smaller than the kinetic energy density associated with the inertia of the plasma. This implies that ${\bf J} \times {\bf B}$ forces are not the dominant dynamical factor responsible for the large-scale properties of the flow in the region. One can argue that the formation of a magnetic tail is an indication that in the cometosheath, the large-scale magnetic field does not dominate the dynamics, but is merely carried around by the superalfvenic flow. If the dynamics were controlled by the magnetic forces, field lines would not bend onto a magnetic tail and the direction of the ion tail would not be essentially in the direction of the local solar wind velocity. We believe that the magnetic field does play a crucial role in the momentum transfer between the solar wind and the cometary plasma, but it is the small-scale, ``turbulent'' magnetic field component that mediates the microscopic interaction between charged particles, leading to the transfer of momentum that we model as an effectively viscous process.

In this paper, we present results of 2D hydrodynamical, numerical simulations of the flow of solar wind and cometary H2O+ ions in the tail and tailward cometosheath of a comet. This is our first attempt to model the interaction of the solar wind with the plasma environment of a comet taking into account viscous-like forces. We review the estimation of the effective Reynolds number of Perez-de-Tejada (1989), by comparing in situ measurements of comet Halley with results from numerical simulations of the viscous-like, compressible flow of the solar wind over a dense, cold, and slow velocity gas representing the plasma tail of a comet. We also study the relative importance of viscous-like forces and the coupling between the fast moving protons of the solar wind and the slow H2O+ ions in the tail by comparing models with different values of the effective Reynolds number, the parameter controlling viscous-like effects, and the effective coupling timescale between both species.

The paper is organized as follows. In Sect. 2 we briefly review our ideas about the origin of the viscous-like forces believed to operate in the solar wind-comet plasma interaction. A description of the problem that we attempt solve and the methodology used, the basic equations, approximations, and main parameters, is given in Sect. 3. Section 4 presents results of a series of simulations with different model parameters. A comparison of our results with in situ measurements at comet Halley is discussed in Sect. 5. Finally, in Sect. 6 we summarize our main results and present our conclusions.

2 On the origin of ``viscosity''

The precise origin of the viscous-like momentum transfer processes invoked in the viscous flow interpretation of the intermediate transition, in both the ionosheath of comet Halley and other ionospheric obstacles to the solar wind, remains unclear. Typical properties of solar wind and cometosheath plasma result in a ``normal'' or ``classical'' viscosity, as it appears in the Navier-Stokes equations when derived from Boltzmann's equation, which can be considered negligible in the dynamics of the flow. Using, for example, properties of the shocked solar wind in the vicinity of the plasma tail measured at comet Giacobini-Zinner, ni = 10 cm-3, |B| = 10 nT and $T = 3 \times 10^5$ K (Bame et al. 1986; Slavin et al. 1986), one can calculate the ``classical'' viscosity coefficient resulting from particle interactions according to Spitzer (1962, Eq. (5)-(55)) to be $\mu \sim 10^{-17}$ g cm-1 s-1 . This extremely low value for the ``normal'' viscosity coefficient, reflects the inefficiency of momentum transport across field lines via collisions in a plasma threaded by a strong, uniform magnetic field.

The relative importance of viscous effects on the flow dynamics can be estimated from the Reynolds number, Re, characterizing the flow. For a viscosity coefficient $\mu$ and gas density $\rho$, $Re = \rho V L / \mu$ measures the ratio of the inertial terms to the viscous terms in the momentum conservation equation. For typical values of the solar wind velocity and mass density in the cometosheath around the tail of comet Giacobini-Zinner, V = 200 km s-1 and $\rho = 1.67 \times 10^{-23}$ gm cm-3, respectively (Bame et al. 1986), and adopting a characteristic lengthscale of 105 km for the variation in the flow velocity (roughly the thickness of the sheath region), we find that the corresponding Reynolds number for the flow, based on the normal viscosity coefficient estimated above ( $\mu \sim 10^{-17}$ g cm-1 s-1), is $Re \approx 10^{11}$. This indicates that viscous forces resulting from the collisions between particles in this environment are negligible in the flow dynamics. If viscous-like processes such as those invoked in this study were to play a role in the dynamics of the flow, they would have to arise from non-collisional processes. In the following paragraphs, we discuss some possible mechanisms responsible for the origin of this effective or anomalous viscosity.

As in Venus and Mars, strong fluctuations in the magnetic field have been measured in the ionosheath of comets Halley and Giacobinni-Zinner (Baker et al. 1986; Scarf et al. 1986; Klimov et al. 1986; Tsurutani & Smith 1986) which possibly affect the transport properties in these plasmas. Cravens et al. (1980) derived an expression for the effective heat conductivity in a plasma threaded by a fluctuating magnetic field, on the basis of a quasi-linear theory of particle transport. Assuming that the Prandtl number of the flow is near unity, Perez-de-Tejada (2005) calculated the effective viscosity coefficient for the ionosheath of Venus based on these results. If we use the same procedure to calculate the effective viscosity coefficient for the solar wind around the tail of a comet (with the conditions measured at comet Giacobini-Zinner, Bame et al. 1986; Slavin et al. 1986), we find $\mu \sim 10^{-11}$ g cm-1 s-1. Even this value for the viscosity coefficient, six orders of magnitude greater than the ``normal'' coefficient described above, still yields viscous forces that are dynamically unimportant. The effective Reynolds number, based on this effective viscosity coefficient and the flow properties described above is Re > 105, indicating that ``viscosity'' resulting from magnetic field fluctuations is negligible, at least as estimated by Cravens et al. (1980). Assuming that the Prandtl number is not very different from unity, as argued to be true by Perez-de-Tejada (2005), we can also neglect the heat conduction caused by these processes.

Another possibility is that, as occurs in many fluid dynamical applications, turbulence, such as that detected in comet Giacobinni-Zinner (Tsurutani & Smith 1986), is characterized (sometimes even defined) by a dramatic increase in the efficiency of transport processes in the flow. In the case of momentum transport this leads to the concept of turbulent or eddy viscosity (Lesieur 1990). Notwithstanding the significant uncertainties in the precise formulation of turbulent ``viscous'' stresses in terms of the large-scale flow properties, particularly in MHD turbulence, the likely importance of eddy viscosity in the flow within the cometosheath is expected because of the large value of the Reynolds number as estimated above.

Finally, as discussed by Shapiro et al. (1995) and Dobe et al. (1999, and references therein), conditions in the ionosheath of Venus and Mars favour the development of plasma instabilities leading to effective wave-particle interactions. If this mechanism were also to operate in cometosheaths, it may lead, as in these planets, to increased coupling between the solar wind and cometary plasma in a viscous-like manner, as suggested by Perez-de-Tejada (1989).

The study of the precise mechanism(s) leading to efficient momentum transfer between solar wind and cometary plasma, as suggested by the present study, is beyond the scope of this work. However,we consider that a detailed study of the importance of viscous-like effects on the flow dynamics in solar wind-comet interactions, in the context of a phenomenological model, is justified by the preceding arguments and may serve to constrain their relative importance of the different flow parameters. It is the purpose of this paper to begin these investigations.

3 Model description

We model the interaction of the solar wind with the plasma tail of a comet using a 2D hydrodynamic, finite difference code for two species (a and b), that is an extension of the single species version presented in Reyes-Ruiz et al. (2008). We include in the dynamical equations a coupling term between both of the species, solar wind protons and cometary ions, which we assume to be H2O+ ions. This term allows the solar wind flow to become mass loaded with cometary ions as they diffuse upwards from the tail, and cometary ions to be accelerated by the fast, streaming solar wind. The coupling term is taken from the work of Szego et al. (2000), who describe the treatment of mass-loaded plasmas. However, to isolate the effects of the viscous-like forces, we do not consider the ongoing creation of new ions in the flow, by photoionization or any other mechanism, as is usually performed in mass loading studies. Since we model only the flow in and around the tail of the comet, the only source of additional ions in our problem is by means of the boundary condition at the left-hand edge of our simulation box (Sect. 3.3). It is clear that the 2D character of our simulations is an approximation of the true problem and may not allow us to study some processes that may be essential to the dynamical evolution of the flow. We make this approximation by considering that this is the first approach to the problem in which viscous-like forces are taken into account. We also neglect the effect of the IMF entrained in the solar wind flow, and leave to future work the study of the dynamical effects of ${\bf J} \times {\bf B}$ forces, although we do not expect these to be dominant in the region (see our arguments in the Introduction section).

Since we are interested in the gas dynamics in the tail, we focus on the region behind the coma, starting from a few times 104 km behind the comet's nucleus and extending downstream to a few times 105 km as illustrated in Fig. 1.

\begin{figure} \par\includegraphics[width=7.8cm,clip]{12376fig1.ps} \end{figure} Figure 1:

Illustration of the computational domain we use in our simulations. The box provides an approximate scale of the simulated region. Image of comet Halley taken the day of the encounter with Giotto by Miller, University of Michigan/CTIO (Brandt et al. 1992).

Open with DEXTER

3.1 Basic equations

The present code solves the equations of mass, momentum, and energy conservation, including terms representing the viscous-like effects and interspecies coupling due to turbulence and/or wave-particle interactions. In cartesian coordinates and in conservative form, for species a, these can be written as

\begin{displaymath} \frac{\partial \vec{U}^a}{\partial t} + \frac{\partial \vec{... ...al x} + \frac{\partial \vec{F}^a}{\partial y} = \vec{S}^{ab} , \end{displaymath} (1)

where

\begin{displaymath}\vec{U}^a = \left( \begin{array}{c} \rho^a \\ \rho^a V_x^a \\ \rho^a V_y^a \\ E_t^a \end{array} \right), \end{displaymath} (2)

\begin{displaymath}\vec{E}^a = \left( \begin{array}{c} \rho^a V_x^a \\ \rho^a... ...+ V_y^a T_{xy}^a) + k_5^a \dot{q}_x^a \end{array} \right) , \end{displaymath} (3)

\begin{displaymath}\vec{F}^a = \left( \begin{array}{c} \rho^a V_y^a \\ \rho^a... ...a + V_y^a T_{yy}^a) + k_5^a \dot{q}_y^a \end{array} \right), \end{displaymath} (4)

and the inter-species coupling term

\begin{displaymath}\vec{S}^{ab} = \left( \begin{array}{c} O \\ \rho^a \nu_{ab... ... \left[ {\bf V}^b - {\bf V}^a \right]^2 \end{array} \right). \end{displaymath} (5)

In the preceding equations, $\rho^a$ is the mass density of gas a, Vxa and Vya are its velocity components, Ta is its temperature, and Eta is the total energy density of species a defined by

\begin{displaymath}E_t^a = \rho^a e^a + \frac{1}{2} \rho^a (V^a)^2 , \end{displaymath} (6)

with ea being the internal energy per unit mass. In Eqs. (3) and (4), the coefficients kia (i=1,5) are the following combinations of dimensionless numbers and the adiabatic index for the gas, $\gamma^a$:

\begin{displaymath}k_1^a = \frac{1}{\gamma^a M_o^2}, \end{displaymath} (7)

\begin{displaymath}k_2^a = \frac{1}{R_{\rm eff}^a}, \end{displaymath} (8)

\begin{displaymath}k_3^a = (\gamma^a - 1), \end{displaymath} (9)

\begin{displaymath}k_4^a = \frac{\gamma^a (\gamma^a-1) M_o^2}{R_{\rm eff}^a}, \end{displaymath} (10)

\begin{displaymath}k_5^a = \frac{\gamma^a}{R_{\rm eff}^a Pr^a}, \end{displaymath} (11)

where the Mach number (Mo), the Reynolds number ( $R_{\rm eff}^a$), and Prandtl number (Pra) for the flow of gas a, are defined, respectively, as

\begin{displaymath}M_o = \frac{V_o}{C_{so}}, \end{displaymath} (12)

\begin{displaymath}R_{\rm eff}^a = \frac{\rho_o V_o L}{\mu_o^a}, \end{displaymath} (13)

\begin{displaymath}Pr^a = \frac{\mu_o^a c_p^a}{\kappa_o^a}\cdot \end{displaymath} (14)

The quantities with subindex o are those used for the normalization of the flow variables and parameters, the reference sound speed is defined as $C_{so} = \sqrt{\gamma^a P_o /\rho_o}$, cpa is the specific heat at constant pressure for gas a, and L is the normalization for the spatial coordinates. For simplicity, we assume that the flow parameters, $\mu$ and $\kappa$, are uniform and that $\mu^a = \mu_o^a$, $\kappa^a = \kappa_o^a$, $\mu^b = \mu_o^b$ and $\kappa^b = \kappa_o^b$.

We note that this form of the interspecies coupling term, although widely used in multispecies gas modelling in various astrophysical scenarios (e.g. Schunk & Nagy 1980; Draine 1986; Cravens 1991; Falle 2003; Van Loo et al. 2009; Szego et al. 2000, and references therein) can be derived strictly from the Boltzmann collision integral only for the case corresponding to Maxwell molecules (see for example Gombosi 1994). We use it to represent a similarly simple, alternative expression for charged particles, and it must be considered an approximation of uncertain validity in our case. Schunk (1977) discussed the modifications to these expressions for the interspecies coupling of electrically charged molecules and in future contributions we shall explore the effect of these modifications. In the present calculations, we chose this approach to modelling a multispecies flow, which follows the dynamics of each species separately, instead of an approach following a single fluid composed of many different species, to clearly distinguish the widely different properties ( $\rho, \ {\bf V}, \ T$, etc.) of solar wind and cometary ions.

The coupling between species represented by the term $\vec{S}^{ab}$ in Eq. (1) is taken from the work of Szego et al. (2000), and has the form of the traditional coupling resulting from binary collisions. The term $\nu_{ab}$ contained in $\vec{S}^{ab}$, reflects the effective result of all processes able to transfer momentum and energy from one species to another. We note that in the adimensional form of the equations $\nu_{ab}$ is actually $t_o \nu_{ab}$, which can be considered to be the ratio of the flow crossing time, to = L/Vo, to the inter-species coupling timescale, $1/\nu_{ab}$. To preserve the symmetry between the coupling terms for both species, guaranteed by the identity $\rho^a \nu_{ab} = \rho^b \nu_{ba}$, we scale $\nu_{ab}$ as $\rho^b$ and $\nu_{ba}$ as $\rho^a$ with a single proportionality constant, $\nu _o$, which we assume to be uniform and constant. In our present code, $\nu _o$ enters as a parameter that can be varied to compare the importance of inter-species coupling to viscous-like forces.

The terms Txxa, Txya and Tyya in Eqs. (3) and (4) represent the components of the viscous-like stress tensor given by

\begin{displaymath}T_{xx}^a = \frac{4}{3} \frac{\partial V_x^a}{\partial x} - \frac{2}{3} \frac{\partial V_y^a}{\partial y}, \end{displaymath} (15)

\begin{displaymath}T_{xy}^a = \frac{\partial V_y^a}{\partial x} + \frac{\partial V_x^a}{\partial y} , \end{displaymath} (16)

and

\begin{displaymath}T_{yy}^a = - \frac{2}{3} \frac{\partial V_x^a}{\partial x} + \frac{4}{3} \frac{\partial V_y^a}{\partial y} \cdot \end{displaymath} (17)

As in multiple fluid dynamics applications (Lesieur 1990), we use the Boussinesq hypothesis in writing the Reynolds stress tensor, we adopt a ``standard'' form for the relation between the viscous-like stress tensor and the large-scale flow velocity, using an effective viscosity coefficient that encapsulates turbulent viscosity as well as the possible effect of wave-particle interactions or any other plasma instabilities leading to an increased coupling between ions in these collisionless plasmas.

In Eqs. (3) and (4), $\dot{q}_x^a$ and $\dot{q}_y^a$ are the components of the effective heat flux vector for species a(based on the Boussinesq hypothesis)

\begin{displaymath}\dot{q}_x^a = - \frac{\partial T^a}{\partial x} , \end{displaymath} (18)

and

\begin{displaymath}\dot{q}_y^a = - \frac{\partial T^a}{\partial y} \cdot \end{displaymath} (19)

Furthermore, we assume throughout this work that both gases are ideal so that

\begin{displaymath}e^a = \frac{1}{\gamma^a -1} \frac{p^a}{\rho^a}, \end{displaymath} (20)

for the equation of state, $p^a = \rho^a R T^a$. We assume that both the solar wind plasma and the cometary plasma, in the tail region, are characterized by an adiabatic index, $\gamma^a = \gamma^b = 5/3$. In Sect. 4, we present some results for $\gamma ^b = 1.25$, and discuss the effects of changing this property of the cometary plasma.

An analogous set of equations and definitions are written for species bby interchanging the superscripts a and b in Eqs. (1)-(20), with $\nu_{ab} = \nu_{ba}$. All equations are coupled by the source term $\vec{S}^{ab}$ in Eq. (1), and solved simultaneously.

3.2 Numerical code

The set of equations described above is discretized in space using 2nd order finite differences, and is advanced in time using an explicit, 2nd order MacCormack scheme (Anderson 1995). The implementation of the scheme is an extension of that described in Reyes-Ruiz et al. (2008), but now with the additional source term $\vec{S}$ in the equations of motion. In MacCormack's scheme, the solution is advanced over one timestep by a sequence of intermediate steps, called the predictor and corrector steps. In the predictor step, an intermediate solution ($U^{\ast}$) is calculated from the values of the physical variables, $\vec{U}_{i,j}^t$, at a given time, t, and position, (xi, yj), according to

\begin{displaymath}\vec{U}_{i,j}^{\ast} = \vec{U}_{i,j}^{t} - c_1 \left[ \vec{... ...} - \vec{F}_{i,j}^{t} \right] + \Delta t \ \vec{S}_{i,j}^t , \end{displaymath} (21)

where $c_1 = \Delta t / \Delta x$, $c_2 = \Delta t / \Delta y$and $\vec{E}^{t}$, $\vec{F}^{t}$, and $\vec{S}^t$are evaluated with $\vec{U}^{t}$ according to Eqs. (3)-(5). This predicted solution is then corrected to obtain the solution at the subsequent time, $t+\Delta t$, using
           $\displaystyle \vec{U}_{i,j}^{t + \Delta t} = \frac{1}{2} \left[ \vec{U}_{i,j}^{t} \right.$ + $\displaystyle \vec{U}_{i,j}^{\ast} - c_1 \left( \vec{E}_{i,j}^{\ast} - \vec{E}_{i-1,j}^{\ast} \right)$  
  - $\displaystyle c_2 \left( \vec{F}_{i,j}^{\ast} - \vec{F}_{i,j-1}^{\ast} \right) + \left. \Delta t \ \vec{S}_{i,j}^{\ast}, \right]$ (22)

where $\vec{E}^{\ast}$, $\vec{F}^{\ast}$, and $\vec{S}^{\ast}$ are computed from $\vec{U}^{\ast}$ using Eqs. (3)-(5). Additional details of the implementation of MacCormack's scheme are given in Reyes-Ruiz et al. (2008).

A final upgrade to our previous code is the ability to handle some types of non-uniform, cartesian grids. For the simulations performed in this work, the grid is defined by a series of (xi, yj) coordinates for which the spacing is arbitrary. In our simulations, the xi points are geometrically distributed from $x_{\rm min}$ to $x_{\max}$ with nx elements. The yjpoints are equispaced from the initial location of the tail (i.e. from y = 0 to y = 1 at 30 gridpoints) and geometrically distributed from y = 1 to $y = y_{\max}$. In both series, the common ratio is 1.02. The 2nd order approximation of the x-derivative of a function f at xi can be easily obtained from the Taylor series expansion of the function at xi-1 and xi+1, and is given by

\begin{displaymath}\left(\frac{df}{dx}\right)_i = \frac{\Delta x_{i-1}^2 f_{i+1... ...lta x_{i-1} \Delta x_{i}^2 - \Delta x_{i} \Delta x_{i-1}^2]} , \end{displaymath}

where $\Delta x_{i} = x_{i+1} - x_{i}$. An analogous expression exists for the y-derivative. This grid allows a higher resolution in the vicinity of the region of strong interaction, while placing the $y = y_{\max}$ boundary sufficiently far away to avoid numerical artefacts in our results.

3.3 Initial and boundary conditions

The solution for the flow is evolved from the following initial conditions. A dense, cold, slow moving plasma representing the tail is located between y=0 and y=1. Both H2O+ and H+ ions are present in the tail, but protons are far less abundant than H2O+. Between y=1.5 and $y=y_{\rm max}$, the gas has the properties of a shocked, hot, fast-moving solar wind that contains both H2O+ and H+ ions, its number density of protons being 50 times higher than H2O+. In all the calculations presented here, we adopted a value Mo = 2 for the Mach number of the shocked solar wind incident on our computational domain. This assumption is made based on the results of Spreiter & Stahara (1980), who computed the gas dynamics of the flow of the shocked solar wind in the ionosheath of Venus. Spreiter & Stahara (1980) found that the flow is characterized by M = 2, as the solar wind crosses the terminator of the planet (the line separating the day and night sides) and heads tailwards. In comets, we take the terminator to coincide approximately with the location of the nucleus. Between y=1.0 and y=1.5, there is transition region where the flow properties change smoothly in an exponential manner from those in the tail to those in the solar wind. The initial density of each species is taken to be

\begin{displaymath}\rho^{a}(t=0) = \left\{ \begin{array}{rl} 0.025 \rho_{\rm ta... ...5mm] \rho_{\rm sw} &\mbox{ if $y>1.5$ } \end{array} \right. , \end{displaymath} (23)

\begin{displaymath}\rho^{b}(t=0) = \left\{ \begin{array}{rl} \rho_{\rm tail} &\... ...0.32 \rho_{\rm sw} &\mbox{ if $y>1.5$ } \end{array} \right. . \end{displaymath} (24)

We assume that both species move initially with the same velocity

\begin{displaymath}V_x^{a,b}(t=0) = \left\{ \begin{array}{rl} V_{\rm tail} &\mb... ....5mm] V_{\rm sw} &\mbox{ if $y>1.5$ } \end{array} \right. , \end{displaymath} (25)

Vya,b(t=0) = 0. (26)

In normalized quantities, $V_{\rm sw} = 1$ and $V_{\rm tail}= 0.01$. For the results shown here, we use, in normalized variables, $\rho_{\rm sw} = 1$ and $\rho_{\rm tail}= 400$. The local temperature of both species is assumed to be the same inside the tail, where $T_{\rm tail}^a = T_{\rm tail}^b = T_{\rm tail}$ and $T_{\rm sw} = 100 T_{\rm tail}$, with $T_{\rm sw} = 1$ being in normalized units. Outside the tail, for y > 1.5, cometary ions are injected with a temperature an order of magnitude lower than the streaming solar wind protons, $T_{\rm sw}^a = 10 T_{\rm sw}^b = T_{\rm sw}$. These choices of temperatures and densities are made to ensure an initial pressure balance between the H2O+ plasma (species b) inside the tail, and the proton plasma (solar wind, species a) outside. For $p^a = \rho^a T^a$ and $p^b = ( m_p^a / m_p^b) \rho^b T^b$, where mpa and mpb are the particle masses of species a and b, respectively, we find that our choice of initial conditions is characterized by a pressure in-balance among each species. In future studies, we propose to determine whether the rapid movement of cometary ions resulting from this initial condition is prevented by the wrapped-around IMF over the comet's tail. Although significantly different from the flow properties at later times in the simulations, for all cases we have studied these initial conditions do not give rise to any long-lasting instability in the flow, so the final state does not depend on their precise form or value.

The boundary conditions are chosen to be consistent with the initial conditions. At the left boundary, $x=x_{\min}$, the flow density and velocity follow exactly those given by the initial conditions in Eqs. (23)-(26). Since the inflow to the comet's tail (y<1) is subsonic, we allow the inflow pressure to float freely as a linear extrapolation of the active mesh values (e.g. Anderson 1995). The right-hand boundary, $x=x_{\max}$, corresponds to the commonly used outflow conditions for supersonic flows, namely the derivatives being zero for all flow variables. We also performed simulations with an outer boundary condition obtained by linearly extrapolating the flow variables, resulting only in minor differences in the last gridpoints before the $x=x_{\max}$ boundary.

4 Results

We performed a series of simulations with different set of parameters $R_{\rm eff}^{a,b}$ and $\nu _o$ to determine the effect of viscous forces and inter-species coupling in the flow dynamics. In all cases considered, the flow evolves from the prescribed initial condition (Eqs. (23)-(26)), passing through a fast transient phase, during which a considerable portion of the mass originally in the tail is eroded by the solar wind exiting our simulation domain. The relevance of this transient phase, which lasts a few tens of solar wind crossing times ( to = L/Vo), relative to the observed features that form in the ion tail, will be analysed in a future publication. In this work, we concentrate on the quiescent stage of evolution because its longer timescale implies that it is more likely to be observed. In all cases, we present results for the flow velocity, density and temperature after a time long enough for the quasi-steady state to have been reached. All results are presented in terms of the normalized quantities defined in the previous section. The appropriate values of L, Vo, $\rho_o$, and To can be chosen for a particular application, as exemplified in Sect. 4 for comet Halley.

To determine the appropriate effective Reynolds number for each species, we consider the following. According to Perez-de-Tejada (1989), the geometry of the flow, measured in situ by the Giotto spacecraft in its fly-by comet Halley in March 1986, implies an effective Reynolds number of about 350 for the shocked solar wind flow above the cometopause along the spacecraft trajectory. For a similar region in the ionosheath of Venus, Perez-de-Tejada (1999) and Reyes-Ruiz et al. (2008), in contrast, estimate a value of the Reynolds number an order of magnitude lower ( $R_{\rm eff} = 20$); this is based on a comparison of in situ measurements (by the Venera 10 and Mariner 5 spacecraft) at Venus with the flow properties derived from a numerical simulation of the viscous-like solar wind-ionosphere interaction. To assess the estimation of Perez-de-Tejada (1989) we conducted simulations with 3 different values of the Reynolds number: (i) A high value, $R_{\rm eff} = 100$, similar to that estimated by Perez-de-Tejada (1989) for comet Halley; (ii) an intermediate value, $R_{\rm eff} = 30$, comparable to the value estimated by Reyes-Ruiz et al. (2008) for the solar wind flow in the ionosheath of Venus; and (iii) a low value, $R_{\rm eff} = 10$, used to verify the trend detected in the results as $R_{\rm eff}$ is decreased.

In most cases, the value of the effective Reynolds number for both species is assumed to be the same. In our view, the lack of knowledge of the precise mechanisms causing the effective viscosity in these plasmas justifies this assumption. However, we did analyse a case with different values of the effective Reynolds number for each species which we describe in the Discussion section.

The value of $\nu _o$ is also varied to explore the relative importance of inter-species coupling, versus viscous forces, which are proportional to $R_{\rm eff}^{a,b}$. We show results for 3 different cases: a strong coupling case characterized by $\nu_o = 1$, which can be interpreted as having a timescale for inter-species coupling equal to the solar wind crossing timescale, to = L/Vo; a medium coupling case, in which the coupling timescale is an order of magnitude greater than the crossing timescale, $\nu _o = 0.1$; and a weak coupling case, for which $\nu _o = 0.01$, so that inter-species coupling effects are much smaller than other dynamical effects.

To compare the state of the flow at the same time in its evolution for all cases, starting from the same initial condition, we have chosen, arbitrarily, to show results at t = 1234, with time units in multiples of the solar wind crossing time. The number of timesteps required to reach this time depends on the model parameters, but in most cases is less than 200 000.

\begin{figure} \par\includegraphics[width=8.5cm,clip]{12376fig2.ps} \end{figure} Figure 2:

Density contours (shades of gray) and flow geometry (velocity vectors) for case 1 ( Rea,b = 30, $\nu _o = 0.1$) after 1234 simulation time units. The top panel shows the configuration for the proton plasma (species a) and the right side panel shows the ``equilibrium'' configuration for cometary H2O+ ions. Density and velocity are in normalized units.

Open with DEXTER

\begin{figure} \par\includegraphics[width=16cm,clip]{12376fig3.ps} \end{figure} Figure 3:

Vertical profiles of the flow properties for case 1 at three different positions; x = 2 ( left column of panels), x=5 ( middle column), and x=8 ( right column). In all cases, grey lines indicate the properties of the proton plasma (species a) and black lines denote the properties of the H2O+ plasma (species b). The top row shows the x component of velocity, Vxa,b, the middle row shows the temperature, Ta,b, and the bottom row shows the mass density, $\rho ^{a,b}$. All quantities are in normalized units.

Open with DEXTER

\begin{figure} \par\includegraphics[width=16.5cm,clip]{12376fig4.ps} \end{figure} Figure 4:

Same as in Fig. 3 but for simulation case 2 ( $R_{\rm eff}^{a,b} = 30$, $\nu _o = 1.0$).

Open with DEXTER

4.1 Effect of inter-species coupling

Our fiducial model, case 1, is characterized by model parameters $R_{\rm eff}^{a,b} = 30$ and $\nu _o = 0.1$. In Fig. 2, we show density contours and the flow velocity for each species. The tail initially contains a uniform density, $\rho^a = 10$ (protons) and $\rho^b = 400$ (H2O+ ions) for y < 1, and after t = 1234 to, a significant portion of the tail has been eroded by the effect of viscous forces and inter-species coupling. A shock wave is evident in the deflection of the flow velocity from the initial uniform distribution imposed by the boundary condition at $x = x_{\rm min}$. Also noticeable is the strong velocity gradient around y = 2, which corresponds to the viscous boundary layer. Both effects are also shown in Fig. 3, where vertical profiles of the x component of velocity, Vx, temperature, T, and mass density, $\rho$, are shown for three different x-positions, $x = 2, \ 5$, and 8; in the left, middle, and right columns of each figure, respectively.

In Fig. 3, the shock front and the boundary layer can be identified at all 3 positions, but they are well separated only for x = 5 and x = 8, shown in the middle and right hand columns, respectively. For a given Vx profile, the shock front corresponds to the uppermost decrease from the uniform velocity (Vx = 1) in the free-flowing solar wind. In the middle panel, corresponding to x = 5, this transition is located approximately at y = 5. A second transition, located approximately at y = 2.5 for x = 5, marks the top of the viscous boundary layer, below which the velocity drops sharply to the very low flow velocities in the middle of the tail. The shock front can also be seen as an increase in both temperature and density in the corresponding panels for each position. The temperature increase and density decrease that are characteristic of viscous boundary layers and have been found in previous studies of viscous flow over a flat plate (e.g. Reyes-Ruiz et al. 2008), is also observed in other cases modelled here. This clearly indicates that the region around y = 2 (at x = 5) is indeed a viscous boundary layer.

\begin{figure} \par\includegraphics[width=16.5cm,clip]{12376fig5.ps} \end{figure} Figure 5:

Same as in Fig. 3 but for simulation case 3 ( $R_{\rm eff}^{a,b} = 30$, $\nu _o = 0.01$).

Open with DEXTER

For case 2, we use the same Reynolds number, $R_{\rm eff}^{a,b} = 30$, as in case 1, but increase the importance of inter-species coupling by using $\nu _o = 1.0$. A figure showing the general flow geometry is not shown since no appreciable differences are found with case 1 (shown in Fig. 2). However, the vertical profiles of flow properties, shown in Fig. 4, clearly illustrate the effect of a much stronger inter-species coupling used in this model. Since, as both species are more tightly coupled, their velocity and temperature distributions tend to be almost identical. The density distribution conforms to the respective boundary conditions of each species; since these conditions are different, there is no reason why both densities should tend to equalize and they do not. Figure 4 shows that the shock front and the boundary layer are not well separated at the rightmost position shown, x = 2. From the shock front height and boundary layer thickness shown in the middle and right columns of Fig. 4, we see that both are proportional to the inter-species coupling (see Sect. 4). The shock front height at x = 5, for example, changes from y = 4.7 for case 1, to y = 5.5 in this case, while the thickness of the boundary layer goes from y = 2.8 to y = 3.2 as we increase the inter-species coupling parameter from 0.1 to 1.

In Fig. 5 we show the results for case 3 characterized by a very weak inter-species coupling, $\nu _o = 0.01$. The general flow geometry (not shown) is very similar to that in Fig. 2. A comparison of Fig. 5 (weak coupling) with Figs. 3 and 4 (medium and strong coupling, respectively), clearly shows that in case 3 the dynamics of both species is essentially uncoupled. The location of the shock front and the top of the boundary layer are different for each species. For example, at x = 5, the shock front and boundary layer can be clearly distinguished only for the cometary H2O+ ions. For the H2O+ ions, the shock front is located at approximately y = 5 and the top of the boundary layer is at y = 2.5, while for protons the shock front and the top of the boundary layer are both located around y = 2.

\begin{figure} \par\includegraphics[width=8.5cm,clip]{12376fig6.ps} \end{figure} Figure 6:

Same as in Fig. 2 but for simulation case 4 ( $R_{\rm eff}^{a,b} = 10$, $\nu _o = 0.1$).

Open with DEXTER

4.2 Effect of viscous-like forces

To analyse the effect of the viscous-like momentum transport between the solar wind and material in the comet's plasma tail, we compare 3 simulations with the same inter-species coupling parameter, $\nu _o = 0.1$, but different values of the effective Reynolds number. Figures 6 and 7 shown the resulting flow geometry and vertical profiles, respectively, for our case 4, characterized by a higher viscosity corresponding to a lower effective Reynolds number, $R_{\rm eff} = 10$, than case 1. Comparing the global geometry of the flow in this case (Fig. 6) with that in a case with a higher Reynolds number, $R_{\rm eff}^{a,b} = 30$ (Fig. 2), we see that after 1234 crossing times, the erosion of the tail is much greater in this high viscosity case for both species. This result is expected as well as the increase in the thickness of the boundary layer as we decrease the effective Reynolds number. This can be clearly seen when comparing the vertical profiles of the flow properties shown in Fig. 3 (medium viscosity) and Fig. 7 (high viscosity). For example, as shown in Fig. 7 for x = 5, the top of the boundary layer increases from y = 2.8 for Rea,b = 30 to approximately y = 3.7 for $R_{\rm eff}^{a,b} = 10$. The increased thickness of the boundary layer as we decrease $R_{\rm eff}^{a,b}$, effectively represents a more blunt obstacle to the solar wind flow. Hence, the height of the boundary layer also increases as we decrease Rea,b. This is also shown in Fig. 7 where, for example at x = 5, the height of the shock front is located at approximately y = 7, about 2 scale units higher than the shock front location for the model with lower effective viscosity (Fig. 3).

\begin{figure} \par\includegraphics[width=16.5cm,clip]{12376fig7.ps} \end{figure} Figure 7:

Same as in Fig. 3 but for simulation case 4 ( $R_{\rm eff}^{a,b} = 10$, $\nu _o = 0.1$).

Open with DEXTER

The trend noted above between high ( $R_{\rm eff}^{a,b} = 10$) and medium effective viscosity ( $R_{\rm eff}^{a,b} = 30$) is confirmed by comparing with results for an even lower viscosity, such as case 5, which corresponds to a model with Rea,b = 100, shown in Figs. 8 and 9. As expected, a decreased viscosity leads to significantly less erosion of the tail than in cases 1 and 4 (medium and high viscosity respectively) as shown in Fig. 8. As discussed above and shown in Fig. 9, the top of the boundary layer also decreases as we increase the Reynolds number, and the location of the shock front also consequently decreases. For example, in the profiles corresponding to x = 5 in Fig. 9, we find that the top of the boundary layer decreases from y = 2.8 for $R_{\rm eff}^{a,b} = 30$, to y = 2.2 for $R_{\rm eff}^{a,b} = 100$. The location of the shock front goes from y = 4.8 for $R_{\rm eff}^{a,b} = 30$ to approximately y = 3.7 for $R_{\rm eff}^{a,b} = 100$.

5 Discussion

Owing to the uncertainty in establishing the physical mechanisms causing the effective viscosity, we have assumed that the effective Reynolds number of both species are identical in the calculations presented above. However, we have also performed simulations in which the effective Reynolds number of each species differ and find that for a medium value of the inter-species coupling, $\nu _o = 0.1$, the results are almost identical to those with a single value for the effective Reynolds number of both species (equal to the effective Reynolds number of species b). For example, for a case with $R_{\rm eff}^a = 100$, $R_{\rm eff}^b = 10$, and $\nu _o = 0.1$, the vertical profile of flow properties for the H2O+ ions at all x locations is almost identical to that shown for case 4 in Fig. 7, which is characterized by $R_{\rm eff}^a = R_{\rm eff}^b = 10$and $\nu _o = 0.1$. The vertical profile of flow properties for the protons, while not identical, is still very similar to case 4. This suggests that viscous stresses, particularly in the species that dominates the mass of the problem, are the dominant factor in the flow dynamics in the context of our model.

As mentioned in Sect. 3, we have assumed in the results presented above that the adiabatic index of both species is the same, i.e. $\gamma^a = \gamma^b = 1.67$. While this value of $\gamma$ can be safely assumed for the solar wind plasma (assuming thermal equilibrium for the species), it is not so clearly valid for the H2O+ plasma in which the excitation of rotational and vibrational degrees of freedom may lead to a lower value of $\gamma$ (again assuming thermal equilibrium for the species). To illustrate the effect of a different, lower value of the adiabatic index for cometary plasma, we have also conducted simulations with a value $\gamma ^b = 1.25$ for the adiabatic index of the H2O+ plasma. This value corresponds to a gas composed of triatomic molecules in thermal equilibrium at a high enough temperature for all molecular degrees of freedom to be excited. The results for this case, $\gamma^a = 1.67$ and $\gamma ^b = 1.25$, with the same effective viscosity and interspecies coupling parameters as case 1 ( $R_{\rm eff}^{ab} = 30$ and $\nu _o = 0.1$) are shown in Fig. 10, which shows the vertical profiles of Vx, T, and $\rho$ for both species in both cases.

When comparing Fig. 10 ( $\gamma ^b = 1.25$) and Fig. 3 ( $\gamma^b = 1.67$), it is clearly evident that if the cometary plasma is characterized by a lower value of the adiabatic index, the heating of the H2O+ plasma is significantly reduced in the boundary layer, since part of the dissipated energy goes to the excitation of the additional degrees of freedom corresponding to the lower value of $\gamma$. This leads to a smaller amount of plasma expansion in the region and a thinner velocity boundary layer. The height of the shock front is consequently reduced. In future contributions, we shall consider the appropriate value of $\gamma$ for the cometary plasma.

5.1 Comparison with in situ measurements

We have found that a comparison of our results with the in situ measurements performed by the Giotto spacecraft, as it flew by comet Halley in March of 1986, is not straightforward. The simplified geometry we use in our simulations to study the interaction in the tail region exclusively, precludes a direct comparison. Nevertheless, some insight into the implications of our results can be obtained from a simplified comparison.

Once a particular application scenario has been chosen, we can establish the values for the characteristic length, velocity, density, and temperature used in the adimensionalization of the equations of motion (Sect. 3). For comet Halley, we interpret the in situ measurements reported in Goldstein et al. (1986), Johnstone et al. (1986), and Perez-de-Tejada (1989) by adopting the parameters L = 150 000 km, Vo = 250 km s-1, $\rho_o = 1.67 \times 10^{-23}$ gm/cm3, and $T_o = 2.5 \times 10^5$ K.

According to Johnstone et al. (1986), the Giotto spacecraft observed 3 distinct transitions in the plasma properties on its inbound trajectory towards comet Halley's nuclear region: (1) the outermost transition occurs about 900 000 km from the point of closest approach and can be identified with the bow shock crossing; (2) the cometopause, where a strong increase in the density of cometary ions is observed, can be located at around 150 000 km from closest approach (Perez-de-Tejada 1989); (3) approximately midway between the shock location and the cometopause, at about 400 000 km from closest approach, the so-called intermediate transition delineates the top of the viscous boundary layer according to the viscous flow interpretation of the solar-wind-comet interaction given by Perez-de-Tejada (1989). Pending a more detailed comparison of the Giotto measurements with the results of our simulations, which should take into account the full geometry of the problem, we identify the cometopause detected in the measurements with the region of very strong H2O+ density increase in our simulations, located at y = 1.0 (approximately) in our normalized units. Based on this assumption, in Fig. 11 we compare the thickness of the boundary layer and the height of the shock front evaluated from our simulation results at x = 5, for models with different effective Reynolds number ( $R_{\rm eff}^{a,b}$) and inter-species coupling parameter ($\nu _o$). As already seen, both the thickness of the boundary layer and the height of the shock front decrease with increasing Reynolds number so that, almost irrespective of the value of $\nu _o$, a low value of $R_{\rm eff}^{a,b}$ is required to explain the measured transition locations.

\begin{figure} \par\includegraphics[width=8.5cm,clip]{12376fig8.ps} \end{figure} Figure 8:

Same as in Fig. 2 but for simulation case 5 ( $R_{\rm eff}^{a,b} = 100$, $\nu _o = 0.1$).

Open with DEXTER

In Fig. 11, the dependence of the transition locations on the value of the inter-species coupling parameter is clearly evident. In simulations with a strong inter-species coupling, the solar wind ions are able to transfer momentum to cometary ions more efficiently giving rise to a thicker boundary layer and higher shock front. The opposite is true when both species are weakly coupled ( $\nu _o = 0.01$). In this case, solar wind ions flow by cometary plasma interacting very weakly. A smaller amount of momentum is transferred between the solar wind and cometary plasma in a situation reminiscent of a high Reynolds number case. Our analysis of scale-heights is not only based on the properties of the velocity profiles in our simulations. As pointed out by Perez-de-Tejada (1989), there are corresponding changes in the density and temperature of the gas as it enters a boundary layer. The heating and expansion characteristics of the viscous boundary layers can also be discerned our results, particularly for cases with a low Reynolds number and high inter-species coupling parameter.

\begin{figure} \par\includegraphics[width=16.5cm,clip]{12376fig9.ps} \end{figure} Figure 9:

Same as in Fig. 3 but for simulation case 5 ( $R_{\rm eff}^{a,b} = 100$, $\nu _o = 0.1$).

Open with DEXTER

\begin{figure} \par\includegraphics[width=16.5cm,clip]{12376fig10.ps} \end{figure} Figure 10:

Comparison of the vertical profiles of flow properties for a case characterized by the same value of $R_{\rm eff}$ and $\nu _o$ as case 1, but with $\gamma ^b = 1.25$. Profiles at x = 2 ( left column of panels), x=5 ( middle column) and x=8 ( right column) are shown. Gray lines indicate the properties of the proton plasma and black lines denote the properties of the H2O+ plasma. All quantities are in normalized units.

Open with DEXTER

We note that similar properties of the flow were also measured by the ICE spacecraft in its flyby through comet Giacobinni-Zinner as discussed in Ip (2004). A comparison of the location of the transition from the sheath region to the so-called transition region and the bow shock location as estimated by Reme (1991), corrected for the different height of the cometopause, yields very similar relative positions to those shown by the dotted lines in Fig. 11 for the transitions in comet Halley. In future work, we will address the differences between the measured flow properties of comets Halley and Giacobinni-Zinner.

5.2 Implications for 3D geometry

We emphasize that the geometry presented in this paper corresponds to a 2D model. In Venus, as discussed by Perez-de-Tejada (1995), the viscous-like interaction between solar wind and ionospheric plasmas occurs preferentially over the magnetic poles of the planet (defined in terms of the incident IMF), where the pile-up of magnetic field lines is less than around equatorial latitudes. According to Perez-de-Tejada (1995), up to about 80$^{\rm o}$ SZA, the piled-up magnetic field over the dayside ionosphere and along the flanks, inhibits to some degree a direct, viscous-like interaction between solar wind and ionospheric plasmas.

If we apply these ideas to the solar wind-comet interaction, this implies that the flow properties we have computed here, correspond more closely to locations over and downstream from the magnetic poles of the comet. For different locations along the tail, the piled-up magnetic field may prevent an efficient viscous-like dragging of ionospheric material, ${\bf J} \times {\bf B}$ forces may be more important, and the flow dynamics may be more accurately modelled in terms of an MHD model, such as those of Wegmann (2002) and Jia et al. (2007). As the IMF is constantly changing direction on a wide range of amplitudes and timescales, the region of viscous-like interaction between the solar wind and cometary plasma, changes with time. Given that the typical IMF orientation is approximately in the ecliptic plane, one should expect the flow within +/-20$^{\rm o}$, measured in the y-direction (as typically defined) from the magnetic poles of the comet, to be more accurately described by our model.

6 Conclusions

We have presented results for the numerical simulation of the interaction between the solar wind and the plasma in the tail of a comet, taking into account the effect of viscous-like stresses previously argued to be important by Perez-de-Tejada et al. (1980). To our knowledge, this is the first time that viscous-like effects have been incorporated into these studies. Our results indicate that there are of 3 distinct transitions in the flow properties: an outermost transition corresponding to the shock front, and innermost transition related to the cometopause and an intermediate transition, which we can identify with the height of the boundary layer characterized by a steep decline in the anti-sunward flow velocity, and the onset of plasma heating and expansion due to viscous-like dissipation. The locations of these transitions depend on the flow parameters, namely the effective Reynolds number of the flow for each species, $R_{\rm eff}^{a,b}$, and the inter-species coupling parameter, $\nu _o$.

By comparing the flow properties from our numerical simulations to the locations of the shock front and intermediate transition measured by the Giotto spacecraft as it approached the nucleus of comet Halley, we have found that, almost irrespective of the strength of the inter-species coupling, $\nu _o$, a low value of the effective Reynolds number, of approximately $R_{\rm eff}^{a,b} \lesssim 20$ for both species, is required to reproduce the measured transition locations. In the context of our model, this implies that the measured flow properties cannot be explained if one does not take into account the viscous-like forces in the interaction of the solar wind and the plasma tail of a comet. Although the conclusions drawn from this study are strictly applicable only to comet Halley and solar wind conditions at the time the in situ measurements were taken, one may speculate that viscous-like processes may be important in general to the solar wind-comet interaction.

This has been the first attempt to include viscous-like forces in the numerical simulation of the interaction of the solar wind with a comet's plasma environment. We note that there are many pending issues still to be addressed that could have potentially important consequences for the details of the solutions obtained for our simplified treatment. First and foremost, the precise forms that we have used for the viscous-like stress and effective interspecies coupling, may need to be revised. As we have argued in the Introduction, plasma properties imply that ``normal'' viscosity is negligible in the region under consideration. Hence, we are invoking an effective viscosity presumably resulting from plasma turbulence and/or wave-particle interactions. However, the precise form of the terms corresponding to viscous-like momentum transfer in the equations of motion (Bousinessq hypothesis) has not been formally determined. As we have also discussed in Sect. 3, the interspecies coupling terms we have used cannot be strictly derived for a plasma such as the one we have modelled. In future contributions, we will analyse the effect of alternative approaches to studying this multispecies system, such as the single fluid approach widely used in laboratory plasma physics.

\begin{figure} \par\includegraphics[width=8cm,clip]{12376fig11.ps} \end{figure} Figure 11:

Height of the shock front, $H_{\rm sh}$ (grey lines with squares) and thickness of the boundary layer, $\delta $(black lines with triangles), as a function of effective Reynolds number, for a set of models with different value for the inter-species coupling parameter, $\nu _o$. Values for these scale-heights correspond to x = 5 in our model. The dotted lines indicate the height of the shock front (grey) and the location of the intermediate transition (black) during the inbound portion of Giotto's flyby through the tail of comet Halley.

Open with DEXTER

There are additional important effects that should also be considered. These include geometrical effects caused by the curvature of the ionosphere, which need to be studied to ensure a more direct, quantitative comparison between in situ measurements by the Giotto spacecraft and the results of simulations. We should study more carefully the interaction of the charged species with neutral gas ejected from the comet, which, especially in the vicinity of the nucleus, is the most abundant species. We should also investigate the effect of the magnetic field on the flow (particularly in the dayside and around the midplane of the near-tail region), and in particular, 3D effects and incoming flow time dependence. Assessing the relevance of these factors was beyond the scope of the present study. This is the subject of work currently in progress and will be reported in future contributions.

Acknowledgements
The authors acknowledge the many useful comments made by the anonymous referee which have lead to its considerable improvement. M.R.R. acknowledges the support of research grant IN109409 of DGAPA-UNAM. H.A. acknowledges the support of research grant IN121406 of DGAPA-UNAM and CONACYT-México grant No. 508807.

References

  1. Alfven, H. 1957, Tellus, 9, 92 [NASA ADS] [CrossRef] [Google Scholar]
  2. Anderson, J. D. 1995, Computational fluid dynamics: The basics with applications, (New York: McGraw-Hill) [Google Scholar]
  3. Baker, D. N., Feldman, W. C., Gary, S. P., et al. 1986, Geophys. Res. Lett., 13, 271 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bame, S. J., Anderson, R. C., Asbridge, J. R., et al. 1986, Science, 232, 356 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  5. Biermann, L. 1951, Z. Astrophys., 29, 274 [NASA ADS] [Google Scholar]
  6. Biermann, L., et al. 1967, Sol. Phys., 1, 264 [Google Scholar]
  7. Brandt, J. C., Rahe, J., Niedner, M. B., et al. 1992, The International Halley Watch Atlas of Large-Scale Phenomena, The University of Colorado, Boulder [Google Scholar]
  8. Cravens, T. E. 1980, J. Geophys. Res., 85, A7778 [NASA ADS] [CrossRef] [Google Scholar]
  9. Cravens, T. E. 1991, in Comets in the Post-Halley Era, ed. R.L. Newburn Jr., et al., 1211 (The Netherlands: Kluwer Acad. Pub.) [Google Scholar]
  10. Cravens, T. E., & Gombosi, T. I. 2004, Adv. Space Res., 33, 1968 [NASA ADS] [CrossRef] [Google Scholar]
  11. Dobe, Z., Quest, K. B., Shapiro, V. D., et al. 1999, Phys. Rev. Lett., 83, 260 [NASA ADS] [CrossRef] [Google Scholar]
  12. Draine, B. T. 1986, MNRAS, 220, 133 [NASA ADS] [Google Scholar]
  13. Falle, S. A. E. G. 2003, MNRAS, 344, 1210 [NASA ADS] [CrossRef] [Google Scholar]
  14. Goldstein et al. 1986, in Proc. 20th ESLAB Symposium/on the Exploration of Halley's Comet, ESA SP-250 [Google Scholar]
  15. Gombosi, T. I. 1994, Gaskinetic Theory, Cambridge University Press (UK: Cambridge) [Google Scholar]
  16. Ip, W.-H. 2004, Comets II, ed. M. C. Festou, H. U. Keller, & H. A. Weaver (Tucson: University of Arizona Press), 605 [Google Scholar]
  17. Jia, Y. D. 2007, J. Geophys Res., 112, A05223 [CrossRef] [Google Scholar]
  18. Johnstone, A., Coates, A., Kellock, S., et al. 1986, Nature, 321, 344 [NASA ADS] [CrossRef] [Google Scholar]
  19. Klimov, S., Savin, S., Aleksevich, Ya., et al. 1986, Nature, 321, 292 [NASA ADS] [CrossRef] [Google Scholar]
  20. Lesieur, M. 1990, Turbulence in fluids, 2nd ed., (The Netherlands: Kluwer, Dordrecht) [Google Scholar]
  21. Meyer-Vernet, N. 1986, Science, 232, 370 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  22. Perez-de-Tejada, H. 1989, J. Geophys. Res., 94, A10131 [NASA ADS] [CrossRef] [Google Scholar]
  23. Perez-de-Tejada, H. 1995, Space Sci. Rev., 72,655 [NASA ADS] [CrossRef] [Google Scholar]
  24. Perez-de-Tejada, H. 1999, ApJ, 525, L65 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  25. Perez-de-Tejada, H. 2000, Planet. Space Sci., 48, 871 [NASA ADS] [CrossRef] [Google Scholar]
  26. Perez-de-Tejada, H. 2005, ApJ, 618, L145 [NASA ADS] [CrossRef] [Google Scholar]
  27. Perez-de-Tejada, H., Orozco, A., Dryer, M. 1980, Astrophys. Space Sci., 68, 233 [NASA ADS] [CrossRef] [Google Scholar]
  28. Perez-de-Tejada, H., Lundin, R., Durand-Manterola, H., et al. 2009, J. Geophys. Res., 114, A02106 [CrossRef] [Google Scholar]
  29. Reme, H. 1991, in Cometary Plasma Processes, ed. A. D. Johnstone, 86, AGU monograph 61, AGU, Washington, D.C. [Google Scholar]
  30. Reyes-Ruiz, M., Díaz-Méndez, E., Perez-de-Tejada, H., et al. 2008, A&A, 489, 1319 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Scarf, F. L., Ferdinand, V., Coroniti, V., et al. 1986, Science, 232, 377 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  32. Schunk, R. W. 1977, Rev. Geophys. Space Phys., 15, 429 [NASA ADS] [CrossRef] [Google Scholar]
  33. Schunk, R. W., & Nagy, A. F. 1980, Rev. Geophys., 18, 813 [NASA ADS] [CrossRef] [Google Scholar]
  34. Shapiro, V. D. et al. 1995, J. Geophys Res., 100, 21289 [NASA ADS] [CrossRef] [Google Scholar]
  35. Slavin, J. A., Smith, E. J., Tsurutani, B. T., et al. 1986, Geophys. Res. Lett., 13, 283 [NASA ADS] [CrossRef] [Google Scholar]
  36. Spitzer, L. 1962, Physics of Fully Ionized Gases, 2nd edn, (New York: Interscience Publishers) [Google Scholar]
  37. Spreiter, J., & Stahara, S. 1980, J. Geophys Res., 85, 7715 [Google Scholar]
  38. Szego, K., Glassmeier, K.-H, Bingham, R., et al. 2000, Space Sci. Rev., 94, 429 [NASA ADS] [CrossRef] [Google Scholar]
  39. Tsurutani, B. T., & Smith, E. J. 1986, Geophys. Res. Lett., 13, 263 [NASA ADS] [CrossRef] [Google Scholar]
  40. Van Loo, S. 2009, MNRAS, 395, 319 [NASA ADS] [CrossRef] [Google Scholar]
  41. Wallis, M. K. 1973, Planet. Space. Sci., 21, 1647 [Google Scholar]
  42. Wegmann, R. 2002, A&A, 389, 1039 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Figures

  \begin{figure} \par\includegraphics[width=7.8cm,clip]{12376fig1.ps} \end{figure} Figure 1:

Illustration of the computational domain we use in our simulations. The box provides an approximate scale of the simulated region. Image of comet Halley taken the day of the encounter with Giotto by Miller, University of Michigan/CTIO (Brandt et al. 1992).

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=8.5cm,clip]{12376fig2.ps} \end{figure} Figure 2:

Density contours (shades of gray) and flow geometry (velocity vectors) for case 1 ( Rea,b = 30, $\nu _o = 0.1$) after 1234 simulation time units. The top panel shows the configuration for the proton plasma (species a) and the right side panel shows the ``equilibrium'' configuration for cometary H2O+ ions. Density and velocity are in normalized units.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=16cm,clip]{12376fig3.ps} \end{figure} Figure 3:

Vertical profiles of the flow properties for case 1 at three different positions; x = 2 ( left column of panels), x=5 ( middle column), and x=8 ( right column). In all cases, grey lines indicate the properties of the proton plasma (species a) and black lines denote the properties of the H2O+ plasma (species b). The top row shows the x component of velocity, Vxa,b, the middle row shows the temperature, Ta,b, and the bottom row shows the mass density, $\rho ^{a,b}$. All quantities are in normalized units.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=16.5cm,clip]{12376fig4.ps} \end{figure} Figure 4:

Same as in Fig. 3 but for simulation case 2 ( $R_{\rm eff}^{a,b} = 30$, $\nu _o = 1.0$).

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=16.5cm,clip]{12376fig5.ps} \end{figure} Figure 5:

Same as in Fig. 3 but for simulation case 3 ( $R_{\rm eff}^{a,b} = 30$, $\nu _o = 0.01$).

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=8.5cm,clip]{12376fig6.ps} \end{figure} Figure 6:

Same as in Fig. 2 but for simulation case 4 ( $R_{\rm eff}^{a,b} = 10$, $\nu _o = 0.1$).

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=16.5cm,clip]{12376fig7.ps} \end{figure} Figure 7:

Same as in Fig. 3 but for simulation case 4 ( $R_{\rm eff}^{a,b} = 10$, $\nu _o = 0.1$).

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=8.5cm,clip]{12376fig8.ps} \end{figure} Figure 8:

Same as in Fig. 2 but for simulation case 5 ( $R_{\rm eff}^{a,b} = 100$, $\nu _o = 0.1$).

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=16.5cm,clip]{12376fig9.ps} \end{figure} Figure 9:

Same as in Fig. 3 but for simulation case 5 ( $R_{\rm eff}^{a,b} = 100$, $\nu _o = 0.1$).

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=16.5cm,clip]{12376fig10.ps} \end{figure} Figure 10:

Comparison of the vertical profiles of flow properties for a case characterized by the same value of $R_{\rm eff}$ and $\nu _o$ as case 1, but with $\gamma ^b = 1.25$. Profiles at x = 2 ( left column of panels), x=5 ( middle column) and x=8 ( right column) are shown. Gray lines indicate the properties of the proton plasma and black lines denote the properties of the H2O+ plasma. All quantities are in normalized units.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=8cm,clip]{12376fig11.ps} \end{figure} Figure 11:

Height of the shock front, $H_{\rm sh}$ (grey lines with squares) and thickness of the boundary layer, $\delta $(black lines with triangles), as a function of effective Reynolds number, for a set of models with different value for the inter-species coupling parameter, $\nu _o$. Values for these scale-heights correspond to x = 5 in our model. The dotted lines indicate the height of the shock front (grey) and the location of the intermediate transition (black) during the inbound portion of Giotto's flyby through the tail of comet Halley.

Open with DEXTER
In the text

Copyright ESO 2010

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.