Interaction of Low-Activity Comets with the Solar Wind

Abstract

The solar wind flow around a cometary atmosphere is calculated using the three-dimensional magnetohydrodynamic model developed by the authors. Emphasis is placed on the case of low-activity comets in which some special features, both quantitative and qualitative, are inherent. The behavior of the flowfield and the magnetic field disturbed by the cometary outflow is analyzed. Some similarity laws that govern the pattern of the interaction between the solar wind and a cometary atmosphere are derived on the basis of the calculated results.

1 Introduction

For about half a century the theoretical investigation of the solar wind flow past cometary comae has been performed on the basis of the model developed by Biermann et al. (1967). This model regards the interaction between the cometary-origin plasma flowing out from the comet surface and the oncoming solar wind as the head-on collision of two counter-streaming supersonic plasma flows with the formation of two shock waves, across which the two flows are decelerated, and a contact discontinuity surface separating the decelerated flows. Originally, it was a two-dimensional (axisymmetric) gasdynamic model, whose constitutive equations included source terms on their right sides that described photoionization of the neutral component of the cometary outflow and the ingress of mass, momentum, and energy into the plasma flow (the so-called mass loading effect). The first self-consistent numerical implementation of this model was given by Baranov and Lebedev (1986).

In the later studies (Baranov and Lebedev 1988, 1993) the parametric investigation was performed and the results of the theoretical predictions were compared with the data measured onboard the Vega, Giotto, and Suisei spacecraft during their mission to comet Halley in March 1986. The agreement between the calculated and experimental plasma parameter (density and velocity) profiles along the spacecraft trajectories turned out to be good throughout almost the entire region ahead of the contact discontinuity surface and somewhat worse in the inner shock layer formed by the cometary outflow. This disagreement was attributed, in particular, to the fact that the model employed did not take account for the interplanetary magnetic field (IMF) and its effect on the plasma flow under consideration.

The further development of the model of the solar wind/cometary coma interaction followed different lines of research. First and foremost, the IMF influence on the interaction pattern and, conversely, the IMF disturbances induced by the cometary ionosphere were investigated within the framework of both two-dimensional (Gombosi et al. 1993; Murawski et al. 1998) and three-dimensional (Gombosi et al. 1996) models. The magnetohydrodynamic (MHD) simulation made it possible, in particular, to reproduce in the calculations the so-called magnetic pile-up region and magnetic cavity surrounding comet Halley, which were detected in the measurements onboard the Giotto spacecraft. Another advancement of the model consisted in taking account for some physical and chemical reactions proceeding in the flowfield, other than the predominant photoionization process taken in consideration in the original model of Biermann et al. (1967). Thus, even in the earlier study by Baranov and Lebedev (1986) the effect of charge exchange between the neutral and charged components was taken into account. Some other processes, such as electron impact ionization, ion-neutral friction, ion-electron recombination, were considered and shown to have an influence on the interaction pattern, particularly, in the near-comet region between the nucleus and the contact discontinuity surface (see, for example, Gombosi et al. 1996; Keremidarska et al. 2011; Kartalev et al. 2012).

An important question concerning the applicability of single-fluid gasdynamic models in simulating the interpenetration of different species is raised in the title of the article by Kartalev et al. (2012). The affirmative answer to this question can be given only within certain limits. In fact, the original Biermann model does not distinguish between the solar-wind protons and the picked-up heavy ions of cometary origin treating both as hypothetical averaged nucleons. For this reason, in recent times one more line of investigation has included simulations on the basis of multifluid models developed, for example, in Benna and Mahaffy (2007) and Rubin et al. (2014).

Nevertheless, it should be noted that, while the multifluid models are usually centered on the behavior of different species in the proximities of comets and the associated chemistry, the single-fluid model is capable of predicting the global interaction pattern. Moreover, as shown by Baranov and Lebedev (2014), it is possible to separately calculate the distributions of the protons and the heavy cometary ions after the calculations within the framework of the single-fluid model have been completed.

A new MHD model of the interaction between the solar wind and cometary ionospheres developed by the authors in the last few years is based precisely on the single-fluid approach. The detailed physical and mathematical description and validation of the model and its numerical implementation can be found in Alexashov et al. (2015) and Baranov et al. (2015). The numerical solutions were obtained using the second-order Godunov method generalized by Miyoshi and Kusano (2005) to include the case of MHD flows. Like in our earlier studies (Baranov and Lebedev 1986, 1988, 1993) and Kartalev’s studies (Keremidarska et al. 2011; Kartalev et al. 2012), the shock-fitting technique was employed in numerically solving the system of governing equations. The advantage of this approach over the shock-capturing technique used, for example, by Gombosi et al. (1996), Benna and Mahaffy (2007) and Rubin et al. (2014) is that it makes it possible to accurately localize the positions of strong discontinuity surfaces (shock waves and contact discontinuities) in the flowfield. Moreover, it allows one to perform calculations on much sparser difference grids, thus reducing computer expenditures. The results of a numerical investigation of the disturbed magnetic fields near spacecraft-explored comets obtained on the basis of this model are presented in our recent paper (Baranov et al. 2015).

One more important question arises as to the range of applicability of continuum (gasdynamic) models to the investigation of such fairly rarefied streams as the solar wind and, particularly, the cometary outflow. That it performs well in the case of medium-activity comets, as comet Halley at perihelion, was brilliantly confirmed by many comparisons of the theoretical predictions for that comet with the spaceboard measurement data (see, for example Baranov and Lebedev 1988, 1993; Gombosi et al. 1993, 1996). However, this is not all so clear for lower-activity comets, the mass production of a comet being, in addition, a function of the heliocentric distance. For this reason, at present attempts are being made to develop kinetic or hybrid models which would be more appropriate in studying the interaction of the solar wind with low-gas-production comets, including the case of comet Churyumov–Gerasimenko, the target comet of the Rosetta mission (Wiele et al. 2011; Rubin et al. 2014; and others).

In this study, we will first revisit the case of comet Grigg–Skjellerup basing upon our new MHD model and invoking also the results of some other studies (Sect. 2). That the theoretical predictions for this comet are a close fit to the data of the measurements onboard the Giotto spacecraft gives a good reason to believe that the model employed shall also provide adequate results for the case of comet Churyumov–Gerasimenko which belongs to the same Jovian family as comet Grigg–Skjellerup. The predictions for the former comet are given in Sect. 3. Finally, we will derive certain similarity laws that govern the pattern of the interaction between the solar wind and cometary comae for low-activity comets, when the continuum model remains still adequate.

2 Revisiting the Case of Comet Grigg–Skjellerup

Comet P/Grigg–Skjellerup (26P/GS) was the third (after P/Giacobini–Zinner and P/Halley) spacecraft-explored comet. The Giotto spacecraft encounter with this comet took place in July 1992, six years after the comet Halley flyby, at a heliocentric distance of precisely 1 a.e. The first results of the onboard measurements were presented in Johnstone et al. (1993) (plasma parameter profiles) and Neubauer et al. (1993) (magnetic field behavior).

A special feature of this event was the low activity of 26P/GS, its gas production at perihelion being two orders less than that of comet Halley ($7 \times 10^{27}$ vs. $5.5 \times 10^{29}$ s⁻¹). This suggested that the interaction between the solar wind and the 26P/GS would be kinetic in nature, because of large gyroradii of newborn cometary-origin ions compared with the interaction region dimensions (Johnstone et al. 1993; Neubauer et al. 1993). However, contrary to the expectations, the qualitative pattern of the interaction turned out to be the same as in the case of the larger Halley comet. The later numerical simulations and their comparisons with the measured data confirmed this fact.

It should be noted that the Giotto–26P/GS encounter occurred under the so-called magnetic cloud conditions, when the IMF strength was unusually high amounting to about 16 nT and, moreover, the magnetic field was almost perpendicular to the solar wind velocity. Under these conditions the cometary ion gyroradius was, according to the estimations made by Flammer and Mendis (1993), about 660 km. Thus, it was one and a half order smaller than the cometocentric bow shock stand-off distance during the Giotto-26P/GS encounter (about 10⁴ km). Of course, this ratio of 1/15 is rather large compared with that in the case of comet Halley which was about 1/250 (under standard IMF conditions the gyroradius was about 2000 km and the bow shock stand-off distance was about 500,000 km). However, it may turn out that it would suffice in order for the hydrodynamic model perform well in this case.

The first simulations for the 26P/GS case were carried out as early as in 1993 (Flammer and Mendis 1993; Schmidt et al. 1993). In those studies the interaction pattern was roughly estimated on the basis of approximate, one-dimensional models.

Lebedev (2000) calculated the solar wind-26P/GS coma interaction using the same axisymmetric, purely gasdynamic model and the same computational code as in Baranov and Lebedev (1986), Baranov and Lebedev (1988) and Baranov and Lebedev (1993). The calculated profiles of the plasma velocity along the Giotto trajectory and the crossings of the bow shock ahead of the comet turned out to be in reasonable agreement with the results of the onboard measurements.

Benna and Mahaffy (2006) considered the same problem within the framework of their multifluid MHD model. Many chemical reactions were taken into account. The calculated ion density, velocity, and field strength distributions along the Giotto trajectory were in good agreement with the recorded observations. The model was also extrapolated to the case of comet Churyumov–Gerasimenko at perihelion in August 2015.

In this study, we revisit the case of 26P/GS employing our new MHD model developed in Alexashov et al. (2015) and Baranov et al. (2015). For the sake of comparison, we take the same values of the relevant parameters of the solar wind and the cometary outflow as in Benna and Mahaffy (2006), namely, the solar wind velocity $V_\infty =380$ km s⁻¹, the solar wind number density $n_{\infty }= 9$ cm⁻³, the cometary outflow velocity $W=0.6$ km s⁻¹, the comet gas production rate $Q=7 \times 27$ s⁻¹, the molecular weight of the cometary-origin particles $\mu =18$ amu, and the characteristic photoionization time $\tau =1.5 \times 10^6$ s. Of the physical and chemical processes other than photoionization only the charge exchange between neutral particles and ions is taken into account in the inner shock layer, the charge-exchange cross-section taken to be $Q_{\mathrm{ch}}=10^{-15}$ cm⁻². In accordance with the above-said, the IMF strength is taken to be 16 nT and the angle between the undisturbed solar wind and IMF vectors is 90°.

The geometric flow pattern is presented in Fig. 1 on both large (a) and small (b) scales. The Giotto trajectory is plotted on the right panel, the closest approach of the spacecraft to the comet being about 200 km. In the case of the strong IMF under consideration the magnetic field has a more pronounced effect on the bow shock (BS) position and shape than in the case of the comet Halley flyby (Baranov et al. 2015) and, as shown below, in the case of comet Churyumov–Gerasimenko under the usual IMF. The magnetic field enlarges the outer shock layer pushing the BS off the comet and drawing the contact discontinuity CD and the inner shock IS nearer to the comet. The calculated crossing of the BS by the Giotto trajectory is in very good agreement with the recorded value at the inbound leg and gives a reasonable fit to the measured result at the outbound leg. We might observe in passing that the measured BS shape is slightly asymmetric in the plane containing the comet–Sun axis and the Giotto path.

The IMF effect on the inner coma structure is all the more strong. In the presence of the magnetic field the CD and IS surfaces are located very close to the nucleus on the sunward side of the comet, the layer between them being very thin. On the night side the presence of the magnetic field makes the IS structure less smeared, namely, the triple point on the IS surface, from which a reflected shock wave and a contact surface proceed, is more pronounced than in the no-field case. In the comet tail the contact discontinuity CD forms a slowly converging tube within which a nonmagnetized fluid flows (since the external magnetic field cannot penetrate through the contact discontinuity surface).

The results of the numerical modeling confirm the conclusion made by Neubauer et al. (1993) in analyzing the measured data that the Giotto spacecraft had not crossed the boundaries of the magnetic cavity ahead of the comet, as distinct from the comet Halley flyby. The same result was obtained in the earlier simulations (Lebedev 2000; Benna and Mahaffy 2006).

Fig. 1

Discontinuity surfaces in the solar-wind flow around the 26P/GS ionosphere in the presence and absence of the IMF (solid and dashed curves, respectively). The comet nucleus is at origin, BS bow shock, CD contact discontinuity surface, and IS inner shock. The detected (Johnstone et al. 1993) and calculated crossings of the bow shock at the Giotto trajectory are marked by white and colored circles, respectively. The numbers at them denote the corresponding distances from the closest approach measured in the units of km

In Fig. 2 we have plotted the distributions of the ion number density and velocity along the Giotto trajectory. On average, the calculated and measured (Johnstone et al. 1993) profiles are in satisfactory agreement. A more or less considerable discrepancy can be noticed in the vicinity of the closest approach, where in the calculations the ion velocity almost vanishes, while its measured value remains fairly perceptible (about 80 km s⁻¹). It should be noted that this situation is repeated in all the three simulations considered, despite the fact that they are based on different physical models and use different numerical methods. In fact, these are a single-fluid, no-MHD model in Lebedev (2000), a multifluid MHD model in Benna and Mahaffy (2006), and a single-fluid MHD model in this study. Nevertheless, they give the same near-vanishing velocity at the comet–Sun line, in contrast with the observed finite value of the velocity.

However, it may be noted that the multifluid model by Benna and Mahaffy (2006) does not dynamically separate the solar-wind protons and the cometary ions. In this sense, all the three models are single-fluid in nature and this might be the reason for the mismatch with the measured velocity.

We also note that it is the solar wind proton distribution that is presented in Johnstone et al. (1993). For this reason, the agreement with the calculated ‘nucleon’ distribution is good only outside the heavy ion mantle, where the solar wind is strongly contaminated by the cometary ions, that is, outside the closest approach vicinity.

Fig. 2

Calculated (dash) and measured (solid, Johnstone et al. 1993) profiles of the ion number density (a) and the velocity (b) along the Giotto trajectory during the 26P/GS flyby

Fig. 3

Calculated (dash) and measured (solid, Neubauer et al. 1993) profiles of the magnetic field strength along the Giotto trajectory during the 26P/GS flyby

The calculated and measured (Neubauer et al. 1993) magnetic field strength profiles are compared in Fig. 3. Again the agreement is good, except for the proximity of the closest approach. Here, the measured value amounts to about 88 nT, while the calculations give only 67 nT [the calculations of Benna and Mahaffy (2006) give the value of 76 nT]. Of course, it should be borne in mind that the actual magnetic field disturbed by the cometary outflow is highly fluctuative in nature and it may fluctuate not only in space, as in Fig. 3, but in time as well. However, the discrepancy in the velocity and field strength distributions near the closest approach remains is as yet poorly understood.

We note that in the spaceboard experiment the jumpwise behavior of the velocity and magnetic field was observed only at the outbound leg of the Giotto trajectory (cf. Figs. 2, 3), while at the inbound leg the spacecraft entry into the shock layer was characterized rather by the fluctuation enhancement than by a step in any parameter distribution at some point of the trajectory. Many reasons for this special feature of the disturbed flow and the magnetic field, together with the existing discrepancies in the calculated and measured parameters at the closest approach, can be suggested. One of these reasons can be the nonuniformity and the time variability of the cometary outflow which must manifest themselves particularly clearly at small distances from the nucleus.

However, except from these details, the single-fluid MHD model used in this study turns out to be fairly adequate in describing the interaction between the solar wind and the cometary atmosphere for the low-activity Grigg–Skjellerup comet.

We note in conclusion that the calculations described above were carried out in the computation domain measuring $3\times 10^5 \ \mathrm{km} \ \times 3\times 10^5 \ \mathrm{km} \ \times 3\times 10^5 \ \mathrm{km}$ and including about 279,600 computation cells. Compared with 4.4 millions cells spreading over 19 levels of refinement involved in the calculations of Benna and Mahaffy (2006), the number 279,600 seems rather modest. This saving in computer resources was achieved thanks to the use of the shock-fitting technique.

3 Results of the Calculations for Comet Churyumov–Gerasimenko

The expected values of the solar wind and the cometary outflow parameters for the case of comet 67P/Churyumov–Gerasimenko are presented in Hansen et al. (2007) for four heliocentric distances ranging from 3.25 to 1.3 au, the latter value corresponding to the perihelion. These reference cases are summarized in Table 1.

Table 1 Expected parameters of the solar wind/cometary coma interactions for 67P/CG

Here, R is the heliocentric distance, μ is the molecular weight of the cometary ions, and λ is the characteristic photoionization length. The other designations are the same as in Sect. 2.

Since in the case of comet Grigg–Skjellerup the hydrodynamic model of the solar wind/cometary coma interaction turned out to be adequate beyond its apparent range of applicability, this offers hope that it would also be fruitful in describing the flow and magnetic field behavior in the case of comet 67P/CG, the target comet of the Rosetta mission, that belongs to the same Jovian family, as comet 26P/GS. This seems reasonable, at least, for the cases d and c in Table 1 in which the comet escorted by the Rosetta spacecraft would be at distances of 1.3 and 2.0 AU from the Sun with the gas production rates from the comet surface of $5 \times 10^{27}$ and $8 \times 10^{26}$ s⁻¹, whereas in the cases b and a (2.7 and 3.25 AU and $8 \times 10^{25}$ and $1 \times 10^{24}$ s⁻¹, respectively) it is agreed that the kinetic model would be more appropriate. Moreover, the early results of the Rosetta mission show a highly nonuniform, spatially and temporarily varying neutral gas atmosphere around the comet (Hässig et al. 2015).

Below we present the results of the calculations within the framework of our MHD model (Alexashov et al. 2015; Baranov et al. 2015) for all the four cases a to d taken from Hansen et al. (2007) (see Table 1) supposing that the onboard measurements will decide what model more closely parallels the actuality in each of the cases listed.

We note that in Hansen et al. (2007) it is supposed that the joint 67P/CG–Rosetta flight would occur under normal IMF conditions, its strength B varying from 1.6 to 4.9 nT and the angle that the IMF vector makes with the distant solar wind varying from 76° to 52° with passage from point a to d. Thus, the IMF effect on the solar wind/comet interaction would be considerably weaker than in the case of 67P/GS.

Figure 4 presents the sections of the three-dimensional grid, on which the computations for the case d were performed, by the xz and xy planes. Here, the x axis coincides with the comet–Sun line, while the IMF vector $\mathbf{B}_\infty$ vector in the undisturbed solar wind lies in the xz plane. Figure 4a gives an idea of the dimensions and the shape of the computation domain which included $160 \times 40 \times 40$ cells (the first number relates to the radial direction and the two last numbers to two angular directions). The tail region boundary is approximately at x = $-$4000 km, while the boundary of the undisturbed solar wind is at the cometocentric distance x = 150,000 km. It was necessary to take so distant entry boundary of the computation domain because of very slow decay of the picked-up ion effect on the preshocked solar wind in the case of the low-activity comet.

The radial size of a computation cell varied from about 10,000 km in the distant solar wind (at $R \approx 150{,}00$ km) to about 50 km in the vicinity of BS ($R \approx 3500$ km) and then diminished to 0.2–0.5 km in the vicinity of CD and in the inner shock layer between IS and CD. The angular sizes of the cells were 4.5° in the meridional direction and 9° in the azimuthal direction.

In this figure, the cross-section of the bow shock BS is visible; it intersects the x axis at a distance of about 3500 km from the nucleus. This is in very close agreement with the data presented in Hansen et al. (2007), where the BS stand-off distance of 3600 km was obtained using the MHD BATSRUS model described in detail in Gombosi et al. (1996). At the same time, the hybrid simulation reported in Koenders et al. (2013) gives the value of only 2000 km.

Fig. 4

Computation domain and difference grid used in the calculations of the solar wind-67P/CG interaction (case d in Table 1) on small (a) and large (b) scales in the xz plane and in the yz plane (c)

In Fig. 4b the grid structure is presented on the larger scale, as the comet nucleus is approached; here, the inner boundary of the computation domain and the cross-sections of the inner shock IS and the contact discontinuity surface CD can be visible, including IS and CD in a part of the tail region. The inner boundary of the computation domain is a sphere (white circle in Fig. 4b) on which the boundary condition of the undisturbed, spherically-symmetric, hypersonic cometary outflow is imposed.

Grids with different structures were used in computing flows far away and near the comet (Fig. 4a, b); they were matched with grid refinement which looks as a shaded sector perpendicular to the x axis in Fig. 4a. Finally, the azimuthal structure of the grid is shown in Fig. 4c in the region between the inner boundary (black circle) and the IS front in the xz (terminator) plane.

It should be noted that, due to the magnetic field influence, the picture plotted in Fig. 4 is slightly asymmetric, although this asymmetry is hardly visible in the plot. Contrariwise, the hybrid simulation for this case (Koenders et al. 2013; Rubin et al. 2014) gives a fairly asymmetric flow pattern. This, together with considerably smaller dimensions of the disturbed flow region, is the main difference of the solution obtained within the framework of the hybrid model from the results of our MHD simulations¹.

The IMF effect is better visible in Fig. 5 in which the shock-wave patterns of the solar wind-67P/CG interaction are compared for case d (Hansen et al. 2007) in the absence and the presence of the magnetic field. The IMF effect on the BS position and shape is marginal (within a few percents) but IS and CD are noticeably influenced by the about 5 nT IMF. Thus, at the comet–Sun line the cometocentric distances of these surfaces in the presence of the magnetic field are half as large as in the case of its absence. As in the case of 26P/GS, the magnetic field makes the IS shape more angulate (in the absence of the field it transforms into an oval).

Fig. 5

Shock-wave patterns of the solar wind-67P/CG interaction (case d in Table 1) on the large (a) and small (b) scales. Solid and dashed curves correspond to the presence (4.9 nT) and the absence of the IMF

The fields of the density and the magnetic field strength are plotted in Fig. 6a, b, together with the streamlines and the magnetic field lines in the proximity of the comet. The ion and magnetic field pile-up regions can be discerned. Thus, on the interval of the comet–Sun line from x = 100 km to the contact surface CD the plasma density increases by about two orders, while the field strength increases approximately twofold. It should be noted that the density increase is monotonic up to the contact surface, whereas the magnetic field profile has a maximum near CD, whereupon the field strength diminishes somewhat. The density increases jumpwise across CD; on the night side of the comet there is formed an extended high-density region inside the contact surface; it is associated with the magnetic cavity within the same surface (black region in Fig. 5b). Of course, it should be borne in mind that this flow pattern can exist only in the case of uniform outflow from the comet surface; the actual pattern will apparently be more complicated.

Fig. 6

Density field and streamlines (a) and magnetic field structure (b) near comet 67P/CG in interaction with the solar wind (case d in Table 1)

The radial density distributions along the comet–Sun axis (0°), the terminator (90°), and in the antisunward direction (180°) are plotted in Fig. 7a. Since on two terminator sides (±90°) the density distributions differ only slightly (low flow asymmetry, as noted above), this distribution is presented only for one side. In this figure, the formation of the so-called heavy ion mantle is clearly visible: within the outer shock layer in between BS and CD the density increases by four orders and more. At the night side of the comet after the passage of the closed CD the density varies only slightly.

The distributions of two components of the magnetic field along the comet–Sun line and the terminator (the third component is zero) are presented in Fig. 7b. An increase in the magnetic field strength in the vicinity of the contact surface is due to the z-component perpendicular to the undisturbed solar wind, especially, at the comet–Sun line. Since the magnetic field cannot penetrate through the tangential discontinuity which is not closed at the night side of the comet, this leads to the formation of an extended narrow region aligned with the negative x direction, from which the magnetic field is absent (cf. Fig. 6b).

Fig. 7

Radial distributions of the density (a) and the magnetic field components (b) for different values of the angle $\varphi$ measured from the comet–Sun line in the xz plane (case d of the solar wind-67P/CG interaction in Table 1)

Finally, in Fig. 8a–c, we have plotted the profiles of the density, the velocity, and the field strength along the comet–Sun line for all the four cases presented in Hansen et al. (2007) and Table 1. Clearly that the plasma parameter distributions are similar and can be reduced into the same curves when using the properly chosen scale length. This will be done in Sect. 4. A slight deviation from the similarity, which can hardly be noticed with an unaided eye, is attributable to the magnetic field effect, whose strength and, particularly, direction vary from case a to case d.

It is interesting to compare the calculated dimensions of the solar wind-67P/CG interaction region (see the caption in Fig. 8a) with the characteristic sizes of the comet (3–4 km). Clearly that in case a the entire interaction region is located within the comet surface. Thus, in this case the MHD model used is quite inapplicable. Moreover, in this case any model intended to simulate the solar wind flow past the comet must take account for the actual shape and dimensions of the nucleus and the distribution of the outflow sources over the comet surface.

In cases b and, particularly, c the bow shock is fairly distant from the nucleus but the inner shock layer lies either within the comet or on its surface. In these cases, it may be supposed that the MHD solution fairly adequately describes the interaction pattern far from the nucleus, since, as shown already in the earlier studies (Brosowski and Wegmann 1973; Baranov et al. 1986), this pattern is low-sensitive to the particular details of the cometary outflow.

Finally, in case d the entire disturbed region between the inner and bow shocks IS and BS is located at a considerable distance from the nucleus, so that there is every chance of successful simulation of this case on the basis of our MHD model, as in the case of comet Grigg–Skjellerup.

Fig. 8

Distributions of the density (a), velocity (b), and magnetic field components (c) along the comet–Sun line for cases a to d of the solar wind-67P/CG interaction in Table 1)

4 Similarity of the Solutions for Low-Activity Comets

We will now dwell on the selection of certain scale lengths which could help to bring the solutions for the problem of the solar wind/cometary coma interaction into the self-similar form.

The control parameters of the problem are the distant solar wind velocity $V_\infty$ and density $\rho _{\infty }$ (or number density $n_\infty$), the gas production rate from the comet surface Q and the cometary outflow velocity W, the molecular weight of the cometary heavy ions $\mu$, and the photoionization time $\tau$. Since the distant solar wind is hypersonic, its temperature $T_\infty$ (or pressure $p_\infty$) is irrelevant. If the magnetic field is taken into account, then the above set of parameters should be supplemented with the IMF strength $B_\infty$ and the angle $\theta _\infty$ the IMF makes with the undisturbed solar wind velocity.

Various combinations of these parameters can lead to different quantities having the dimension of the length. Of these we will choose the following length scale

$$L_\zeta =\frac{Q\mu }{4\pi n_\infty V_\infty W\tau }.$$

(1)

Dividing this expression by $W\tau$ we obtain the following dimensionless control parameter

$$\zeta =\frac{Q\mu }{4\pi n_\infty V_\infty W^2 \tau ^2}.$$

(2)

Obviously, this parameter controls the rate of the solar wind mass-loading by cometary ions via the gas outflow from the comet surface (parameter Q) and their subsequent photoionization (parameter $\tau$). Together with the velocity ratio $\alpha =W/V_\infty$, it determines the solution of the problem, at least, in the absence of a magnetic field.

Throughout this Section we will use the results of our parametric calculations performed on the basis of our model (Sect. 1) in the absence of the magnetic field. These calculations were carried out on a wide range of the control parameter $\zeta$ from the conditions corresponding to low-active comets to highly active ones and for velocity ratios from 1/600 to 1/200 which apparently covers the possible range of the solar wind and cometary outflow velocities.

In Fig. 9 we have plotted the dependence of the dimensionless (normalized by $L_\zeta$) cometocentric bow shock stand-off distance ${\varDelta }_0$ as a function of the parameter $\zeta$. Clearly, on the low-$\zeta$ range ($\zeta < 5 \times 10^{-3}$) the dimensionless solution for the bow shock no longer depends on the control parameter $\zeta$ and can be written as follows:

$${\varDelta }_0/L_\zeta =2.3$$

(3)

The same can be said about the $\zeta$-dependences of the cometocentric distances of the contact discontinuity CD (${\varDelta }_{0,CD}$) and the inner shock IS (${\varDelta }_{0,IS}$): these quantities normalized by $L_\zeta$ remain constant on the above-mentioned $\zeta$ range. However, whereas the ${\varDelta }_0(\zeta )$ function is almost independent of the second control parameter $\alpha =W/V_\infty$, the analogous dependences for CD and IS can considerably vary depending on $\alpha$. Thus, the variation of $\alpha$ on the range from 1/600 to 1/200 leads to about twofold variation in ${\varDelta }_{0,CD}$ and ${\varDelta }_{0,IS}$. However, the absolute values of the two latter distances are very small compared with ${\varDelta }_0$ (not >5 % of ${\varDelta }_0$), so that their variation affects only a small portion of the vast shock layer between CD and BS.

Fig. 9

Dependence of the dimensionless bow-shock stand-off distance and the Mach number ahead of the bow shock at the comet–Sun line on the similarity parameter $\zeta$. CG,c, CG,d, GS, H, and HB correspond to the conditions of the solar wind interaction with comets Churyumov–Gerasimenko (cases c and d in Table 1), Grigg–Skjellerup, Halley, and Hale–Bopp

In Fig. 9 the dependence of the Mach number $M_0$ directly ahead of the bow shock at the axis of symmetry is also presented. For low-active comets $M_0 \approx 1.55$.

As noted in Sect. 3, the results of the calculations indicate the similarity of the flow parameter profiles in the outer shock layer calculated for different solar wind conditions (Fig. 8a), when these parameters are referred to the corresponding values in the distant solar wind. This is confirmed by the results of many calculations made by us in this and other studies. Thus, for low-activity comets in dimensionless variables there exists the unique limiting solution of the problem of the interaction between the solar wind and the cometary atmosphere. The presence of a moderate magnetic field has only a slight effect on the interaction pattern in most of the outer shock layer (cf. Fig. 5b), although the inner region of interaction can be influenced by the magnetic field (Fig. 5a), as well as by different physical and chemical processes, other than photoionization, occurring in the interaction zone.

The self-similar limiting solution is presented along the comet–Sun line in Fig. 10. Here, the profiles of the velocity and the nucleon, solar-wind proton, and heavy ion densities are plotted. The two last distributions are obtained within the framework of the approach developed in Baranov and Lebedev (2014). The slow decay of the disturbances induced in the preshocked solar wind by the mass-loading is visible: at a cometocentric distance twice as large as the bow shock stand-off distance the solar wind velocity has already lost about one quarter of its undisturbed value.

Fig. 10

Profiles of the ion velocity (a) and the densities (b) of the nucleons (solid), protons (dash), and heavy ions (dash-and-dot) along the comet–Sun line in the dimensionless variables in the limiting self-similar case of low-activity comets. CD and BS correspond to the contact discontinuity and the bow shock, while IB _c and IB _m relate to the ion composition boundaries with respect to the concentration and the mass

Ahead of the bow shock and after the passage across it the main fraction of charged particles in the mass-loaded solar-wind flow is provided by protons. In certain studies the notion of the ion composition boundary (ICB) was introduced. It is defined as the locus of the points in the shock layer, where the solar-wind proton density becomes equal to that of the cometary ions. The ICB position can be determined from the plots presented in Fig. 10b. The thickness of the layer between ICB and BS at the axis amounts to about 0.4 of that of the entire shock layer, being slightly greater at the terminator. The above-said pertains to the mass densities. Obviously, when the number densities (concentrations) are considered, the ICB must be located nearer to CD; thus, at the comet–Sun axis the thickness ratio will be 0.08 rather than 0.4. Here, the so-called ion pile-up region, or heavy ion mantle, begins, where a small amount of protons is dwarfed to insignificance in the middle of an enormous crowd of the heavy cometary ions.

It should be noted that, of course, the limiting interaction pattern takes place not only in the near-axial region but also throughout the entire outer shock layer from the comet–Sun axis to the terminator and even in the tail region. In the sunward region the shape of the bow shock can be satisfactorily approximated by the paraboloid of revolution

$$y=k\sqrt{{\varDelta }_0({\varDelta }_0 - x)}$$

(4)

where the coefficient $k \approx 2.4$.

We will now return to dimensional variables. From Eqs. (1) and (3) we obtain the following expression for the bow shock stand-off distance for low-active comets

$${\varDelta }_0=2.3\frac{Q\mu }{4\pi n_\infty V_\infty W\tau }$$

(5)

In particular, this means that on the range considered the stand-off distance linearly depends on the cometary gas production rate Q. The same is true for the other characteristic distances ${\varDelta }_{0,CD}$ and ${\varDelta }_{0,IS}$. These results are also confirmed by the calculated results presented in Hansen et al. (2007) (see Table V in that paper).

It is interesting to compare the result given by Eq. (5) with that of the one-dimensional analysis made in the pioneering work of Biermann et al. (1967) and recently revisited by Koenders et al. (2013). We will repeat this analysis in terms of the Mach number rather than the mass flow rate, as was done in the above-mentioned studies, since this approach seems more obvious.

The equations governing the one-dimensional flow in the preshocked solar wind read as follows:

$$\begin{aligned}&\frac{d(\rho V)}{dx} = Q_m(x) = \frac{Q\mu m_p}{4\pi Wx^2\tau }\exp \left( - \frac{x}{W\tau }\right) , \\ & \quad \frac{d(\rho V^2 + p)}{dx}=Q_i \approx 0, \\ & \quad \frac{d}{dx}\left( \frac{1}{2}\rho V^3 \ + \ \frac{\gamma }{\gamma - 1}pV \right) =Q_e \approx 0. \end{aligned}$$

(6)

Here, the x axis is directed toward the Sun, with the origin at the comet nucleus, $\gamma$ is the specific heat ratio, and $m_p$ is the proton mass. The source terms $Q_m$, $Q_i$, and $Q_e$ on the right sides of the equations describe the mass, momentum, and energy ingress into the solar wind owing to the loading by newborn cometary ions. The expression for $Q_m$ is the solution of the continuity equation for the cometary neutrals with account for their loss by photoionization.

The flow one-dimensionality is a good approximation for the preshocked solar wind, since both spaceboard measurements and numerical simulations show that this flow is near-plane-parallel (see, for example, the comparison of the measured and calculated velocity profiles along the Suisei spacecraft trajectory in Fig. 4 in Baranov and Lebedev 1993). Setting the momentum and energy sources equal to zero seems more questionable but it suits for the first approximation. Integrating Eq. (6) we obtain

$$\begin{aligned}&\rho V = \rho _{\infty } V_{\infty } + \int \limits _\infty ^x Q_m(x) \, dx, \\& \quad \rho V^2 + p = \rho _{\infty } V_{\infty }^2, \\& \quad \rho V^3 + \frac{2}{\gamma - 1}pV = \rho _{\infty } V_{\infty }^3. \end{aligned}$$

(7)

Here, it is taken into account that in the hypersonic distant solar wind the pressure $p_\infty =0$.

We will introduce the Mach number $M=V/c=\sqrt{\rho V^2/\gamma p}$ and express the flow parameters in terms of M using two last equations (7). Then after some algebra we arrive at the equation determining the implicit dependence of the Mach number on the cometocentric distance x

$${\varPhi }(M(x))=\frac{\rho V}{\rho _{\infty } V_{\infty }}= \frac{\gamma ^2}{\gamma -1}\frac{\left[ 2+(\gamma -1)M^2\right] M^2}{\left( 1+\gamma M^2\right) ^2} = 1 + \frac{1}{\rho _{\infty } V_{\infty }} \int \limits _\infty ^x Q_m(x) \, dx$$

(8)

However, this solution does not admit the passage through the sonic point ($M=1$), where, as can readily be seen, the derivative $d{\varPhi }/dM=0$ and, accordingly, $dM/dx=\infty$. This point is a limiting point of the solution or, to be more precise, it is the vertex of the limiting line (surface) at the axis of symmetry. The properties of these lines were considered in detail in Mises (1958), where it was asserted that the appearance of a limiting line in the stationary solution can be considered as the mathematical criterion of the physical breakdown of the steady flow under consideration (see Chap. 5, Sect. 4 in that book). Actually, the solution with a limiting line cannot be realized and in the actual flow this line must be preceded by a shock wave. The shape and position of this shock can be determined by solving the boundary value problem of the solar wind/cometary coma interaction in the class of weak solutions, that is, admitting strong discontinuities.

Returning to the limiting line, we can easily determine its position $x_\ell$ at the axis by setting $M=1$ in Eq. (8). Moreover, in the limiting case of low-active comets we can set the exponential in the first equation (6) equal to unity and thus derive the solution in the quadrature form, as follows:

$$x_\ell =\left( \gamma ^2 - 1\right) \frac{Q\mu }{4\pi n_\infty V_\infty W\tau }$$

(9)

This is precisely the expression obtained in Biermann et al. (1967). It should be borne in mind that this expression determines the location of the purely mathematical limiting line rather than that of the physical shock wave. Since Biermann’s days many attempts have been made to relate these quantities, that is, to establish the relation of the form ${\varDelta }_0=f(x_\ell )$. However, all these attempts should be recognized as unsuccessful and the reason for this is obvious: the bow shock wave cannot be determined without solving the boundary value problem invoking the processes occurring throughout the entire interaction region. It might only be expected that the expressions for the two quantities would be of the same functional form, in which case Eq. (9) could be used for calculating ${\varDelta }_0$ by an appropriate choice of a coefficient ahead of this expression.

In fact, it can easily be seen that the functional forms of Eqs. (5) and (9) are the same and that Eq. (5) obtained on the basis of the calculated results gives the value of the required coefficient.

We will now consider the question of the range of applicability of the limiting solutions derived. As can be seen in Fig. 9, the ${\varDelta }_0 (\zeta )$ curve starts to deflect from its asymptotic value of 2.3 at $\zeta \approx 5 \times 10^{-3}$. This value can be taken for the right end of the applicability interval. As for the left end, it is determined by the applicability of the MHD model itself, although, as shown in Sect. 3, purely mathematical self-similar solutions of our problem exist for vanishing values of $\zeta$.

In Fig. 9 we have marked certain positions on the $\zeta$-axis that indicate the conditions characteristic of several comets, namely, Churyumov–Gerasimenko in cases c and d, according to Hansen et al. (2007) and Table 1, Grigg–Skjellerup, Halley, and Hale–Bopp; for comet Halley the dashed interval on this axis corresponds to the week in March 1986 between the comet Halley flyby by Suisei and the two Vegas and the Halley/Giotto encounter. Clearly that the right end of the applicability range corresponds to the 26P/GS conditions.

With further increase in $\zeta$ the dimensionless bow shock stand-off distance diminishes, while the dimensional distance continues to increase, though at a slower rate. It is interesting to note that the inflection point on the ${\varDelta }_0(\zeta )$ curve corresponds to the comet Halley conditions.

The control parameter $\zeta$ and the corresponding scale length $L_\zeta$ become unfit for analyzing the outflows of medium- and high-gas-production comets and their interaction with the solar wind. In particular, formulas of the type of Eqs. (5) or (9) cannot be used for determining the bow shock stand-off distance, which was first noticed by Galeev et al. (1985).

In this case, another scale length $L_\beta$ becomes more appropriate, namely

$$L_\beta =\left( \frac{Q\mu }{4\pi n_\infty V_\infty }\right) ^{1/2} = \zeta ^{-1/2} L_\zeta .$$

(10)

In this case, it is more convenient to take the control parameter of the problem in the form:

$$\beta =\left( \frac{Q\mu }{4\pi n_\infty V_\infty W^2 \tau ^2}\right) ^{1/2}= \zeta ^{1/2}.$$

(11)

Parameters $\beta$ and $L_\beta$ were used in our earlier studies concerned with comet Halley (Baranov and Lebedev 1988, 1993).

In Fig. 11 we have plotted the dependence of the dimensionless (normalized by $L_\beta$) cometocentric bow shock stand-off distance ${\varDelta }_0$ as a function of the parameter $\beta$. In these variables, the range of low-activity comets takes up a minimum space near the origin. Here, the ratio ${\varDelta }_0/L_\beta$ increases as the linear function of the control parameter $\beta$. This increase stops short of comet Halley, so that on the whole range between comets Halley and Hale–Bopp this ratio varies only slightly (within about 10 %).

Fig. 11

Dependence of the dimensionless bow-shock stand-off distance at the comet–Sun line on the similarity parameter $\beta$. GS, H, and HB correspond to the conditions of the solar wind interaction with comets Grigg–Skjellerup, Halley, and Hale–Bopp

Then taking the mean value of ${\varDelta }_0/L_\beta =0.62$ on this interval and recalling the definition of the scale length $L_\beta$ (Eq. 2) we arrive at the following dependence of the dimensional bow shock stand-off distance on the relevant parameters of the problem for the range of medium- and high-activity comets

$${\varDelta }_0=0.62\left[ \frac{Q\mu }{4\pi n_\infty V_\infty }\right] ^{1/2}.$$

(12)

The comparison of Eqs. (5) and (12) shows that on the medium activity range the dependence of ${\varDelta }_0$ on the cometary outflow rate Q is weaker than in the case of low-gas-production comets ($\sim$ $Q^{1/2}$ rather than $\sim$ $Q$). Moreover, in the former case ${\varDelta }_0$ turns out to be independent of the photoionization scale time $\tau$. These facts make in this case the comparison of the measured and theoretical data more favorable, since the relevant parameters of the solar wind/cometary outflow interaction are known with a poor accuracy (especially, the photoionization time $\tau$).

For high-activity comets even the dimensional bow shock stand-off distance starts to diminish. As for the cometocentric distances of CD and IS, these continue to grow. As a result, the outer shock layer starts to shrink and the interaction between two counter-streaming flows becomes increasingly more gasdynamic in nature. This is due to the fact that in the case of high-activity comets an increasingly greater portion of cometary neutrals has a time to ionize in the proximity of the nucleus, thus introducing less disturbances in the solar wind flow. However, this issue is beyond the scope of this study devoted to low-activity comets.

5 Summary

The three-dimensional magnetohydrodynamic model developed by the authors is applied to calculate the solar wind flow around low-activity comets. Despite the fact that this model seems inadequate in describing the processes occurring in this case, the comparison of the theoretical results with the measurement data obtained onboard the Giotto spacecraft during its encounter with comet Grigg–Skjellerup shows that in this case the model performs well. This gives promise that the theoretical predictions concerning the flowfield and magnetic field behavior during the encounter of the Rosetta spacecraft with comet Churyumov–Gerasimenko obtained on the basis of our model will be useful in interpreting the data to be obtained by the spaceboard experiments. Finally, some similarity laws that govern the pattern of the solar wind/cometary atmosphere interaction for low-activity comets are established on the basis of the calculated results. In the dimensionless variables introduced this interaction is described by a unique solution on the entire low-activity range.