1 Introduction

Cometary comae and their interaction with the solar wind have been the target of various studies in the past. Data sent back to Earth from missions to comet, such as the Giotto mission to comet 1P/Halley (Reinhard 1986), revealed a wealth of features in the plasma interaction around an active comet (Johnstone et al. 1986; Rème et al. 1986). Multiple plasma boundaries could be detected (Cravens 1989), ranging from the bow shock (Coates et al. 1997) at large cometocentric distances to the diamagnetic cavity, a magnetic field-free zone close to the nucleus, as measured by the magnetometer (Neubauer et al. 1986; Neubauer 1987). Further boundaries have been identified and reviewed [e.g. Coates and Jones (2009)]. The Ion Mass Spectrometer (Balsiger et al. 1986) observed a pile-up of ions associated with an increase in the electron temperature, which shuts off ion–electron recombination (Häberli et al. 1995). As comets cover a range of heliocentric distances, their interaction with the Sun undergoes significant variations. The activity of the comet increases with decreasing heliocentric distance due to solar insolation. Furthermore the interplanetary magnetic field increases as the comet approaches the Sun. Hansen et al. (2007) investigated the plasma environment of comet 67P, the target comet of ESA’s Rosetta mission, using a suite of models including single fluid MHD simulations and hybrid-kinetic plasma simulations. Koenders et al. (2013) investigated the stand-off distance of the cometary bow shock for typical production rates and solar wind conditions of the same comet and later on also studied the close nucleus interaction of 67P after formation of the diamagnetic cavity (Koenders et al. 2015).

In Rubin et al. (2014b) we investigated the solar wind neutral gas interaction of a low activity comet at 2.7 AU from the Sun. This peculiar environment is governed by large gyroradii of the cometary heavy ions. Hybrid kinetic simulations, which are able to reproduce gyro radius effects, are ideal for this situation since they capture the particle nature by modeling the interaction with a set of representative model particles that experience the same forces as real ions in the coma (Lipatov et al. 1997). We have shown however, that multifluid magnetohydrodynamics (MHD) is also able to reproduce some of the basic features of this interaction, e.g. the large scale gyration, at a much smaller computational cost compared to kinetic simulations. Corresponding observations were made by the AMPTE mission by releasing barium (Coates et al. 1988) and lithium (Coates et al. 1986; Moebius et al. 1986) ions in the solar wind with subsequent detection downstream by the United Kingdom Satellite.

Here in this work, we pursue our multifluid MHD simulations for comet 67P along its orbit, i.e. from 1.8 AU to perihelion at approximately 1.25 AU. The plasma environment is subject to change during this time period: the Mach cone that formed at large heliocentric distances steepens into a bow shock. Furthermore, the bow shock stand-off distance upstream of the comet increases as the comet becomes more active (Hansen et al. 2007; Koenders et al. 2013, 2015).

Close to the nucleus, a diamagnetic cavity forms—inside the region where outflowing neutral gas cools and drags the cometary plasma due to abundant charge-exchange and collisional interactions. Inside the cavity the magnetic field drops to very low values (Cravens 1989; Gombosi et al. 1996), which has been observed by the Giotto magnetometer at comet 1P/Halley (Neubauer 1987). The outflowing plasma is cold and supersonic and at the inner shock the plasma transitions from supersonic to subsonic speeds. This forms a natural boundary that prevents the magnetized solar wind from penetrating this innermost region, which therefore remains devoid of any significant magnetic field assuming that comets themselves do not carry any remnant large-scale magnetization as shown by Auster et al. (2015) for comet 67P. This inner shock has been observed by the Giotto ion mass spectrometer (Balsiger et al. 1986; Rubin et al. 2009). Like the bow shock, the size of the diamagnetic cavity depends on the mass loading of the comet.

This work together with our previous simulations in Rubin et al. (2014b) give an overview of the expected plasma environment of comet 67P throughout the Rosetta comet escort phase that started in August 2014 and will nominally end in December 2015.

2 Multifluid MHD Model

The baseline MHD equations and the associated source terms containing the involved physical processes are the same as in Rubin et al. (2014b). The set of Eqs. (1) contains the continuity (I), momentum (II), and pressure equations (III) that are solved for each individual ion species. Furthermore the electron pressure is calculated independently to derive a self-consistent electron temperature (IV). The magnetic field is obtained from the combination of Faraday’s law of induction (V) and the generalized Ohm’s law (VI).

$$\begin{aligned} & I:\quad \frac{{\partial \rho_{s} }}{\partial t} + \nabla \cdot \left( {\rho_{s} \varvec{u}_{s} } \right) = \frac{{\delta \rho_{s} }}{\delta t} \\ & II:\quad \frac{{\partial \rho_{s} \varvec{u}_{s} }}{\partial t} + \nabla \cdot \left( {\rho_{s} \varvec{u}_{s} \varvec{u}_{s} + Ip_{s} } \right) - Z_{s} e\frac{{\rho_{s} }}{{m_{s} }}\left( {\varvec{E} + \varvec{u}_{s} \times \varvec{B}} \right) = \frac{{\delta \rho_{s} \varvec{u}_{s} }}{\delta t} \\ & III:\quad \frac{{\partial p_{s} }}{\partial t} + \left( {\varvec{u}_{s} \cdot \nabla } \right)p_{s} + \gamma p_{s} \left( {\nabla \cdot \varvec{u}_{s} } \right) = \frac{{\delta p_{s} }}{\delta t} \\ & IV:\quad \frac{{\partial p_{e} }}{\partial t} + \left( {\varvec{u}_{e} \cdot \nabla } \right)p_{e} + \gamma p_{e} \left( {\nabla \cdot \varvec{u}_{e} } \right) = \frac{{\delta p_{e} }}{\delta t} \\ & V:\quad \frac{{\partial \varvec{B}}}{\partial t} = - \nabla \times \varvec{E} \\ & VI:\quad \varvec{E} = - \varvec{u}_{e} \times \varvec{B} - \frac{1}{{n_{e} e}}\nabla p_{e} \\ \end{aligned}$$
(1)

with ρ s , u s , p s , and Z s the mass density, velocity, pressure, and charge state of ion species s, respectively, and correspondingly with index e for the electrons. I is the identity matrix, e the unit charge, and γ the adiabatic index for which we assumed 5/3, i.e. the ratio of the specific heats \(\gamma = \left( {f + 2} \right)/f\) with f = 3 considering only motional degrees of freedom. This is of course only an approximation for the mix of solar wind protons, alpha particles, and the corresponding photoionization and dissociation products of water, H2O+, OH+, O+, and H+. The model assumes charge neutrality, i.e. the electron number density is derived from the number density and charge states of the different ions \(n_{e} = \mathop \sum_{s = ions} Z_{s} n_{s}\).

The right hand sides of I–IV in Eq. 1 contain the source and sink terms that treat sources and losses in the mass, momentum, and pressure equations due to photoionization, electron impact ionization, charge exchange, elastic as well as inelastic collisions, and ion–electron recombination. The details of these terms can be found in Rubin et al. (2014b) including the rates for the collisional interactions and ion–electron recombination. Most important for the formation of a diamagnetic cavity boundary are charge-exchange reactions between ions and the neutral gas that emanates from the nucleus. We therefore list here the ion-neutral charge exchange rate used in our model as \(k_{ns} = 1.7 \times 10^{ - 15}\,{\text{m}}^{3} \,{\text{s}}^{ - 1}\) taken from Gombosi et al. (1996).

The neutral gas distribution is derived using a Haser-type description (Haser 1957). Remote sensing and in situ instruments on the Rosetta spacecraft have shown that the neutral gas environment is not spherically symmetric (Hässig et al. 2015), which is in line with many earlier findings at other comets, for instance by the EPOXI spacecraft at comet Hartley 2 (A’Hearn et al. 2011) or during the Deep Space 1 flyby at comet 19P/Borrelly (Young et al. 2004). 67P exhibits variations in the neutral gas production with longitude and latitude above the nucleus, solar phase angle, and the position of the Sun with respect to the comet because the illumination conditions, i.e. the Sun’s elevation angle in the comet fixed frame, will change on the way to perihelion and the distribution of the neutral gas is constantly changing with the rotation of the comet (Bieler et al. 2015). Nevertheless, for the purpose of this investigation we use uniform outgassing to reduce the number of free parameters in our simulation runs and to facilitate the comparison between the simulations at different heliocentric distances. According to the Haser model the density of the neutral gas is

$$n_{h} = \frac{{Q_{n} }}{{4\pi u_{n} \left| \varvec{r} \right|^{2} }}e^{{\left( { - \frac{{\left| \varvec{r} \right| \cdot \upsilon_{io} }}{{u_{n} }}} \right)}}$$
(2)

Here Q n is the neutral gas production rate, r is the cometocentric distance, u n is the neutral gas outflow velocity, which we keep constant at 800 ms−1, and \(\upsilon_{io}\) the ionization rate.

We also keep the solar wind boundary conditions fixed for each individual heliocentric distance. The set of boundary conditions are listed in Table 1. These values have been derived according to the work by Hansen et al. (2007) and also reflect the boundary conditions of our earlier paper (Rubin et al. 2014b) with the exception that here we use a Parker spiral with a non-vanishing magnetic field component, B x , along the solar wind flow direction. As the comet approaches the Sun, the magnetic field as well as the density and temperature of the solar wind increases. Furthermore, both the comet’s activity increases and the photoionization scale length decreases leading to enhanced mass loading of the solar wind.

Table 1 The boundary conditions for the 6 cases between 1.3 and 1.8 AU are scaled from nominal values at 1 AU (nSW = 10 cm−3, |B| = 7 nT, fi = 1.0 × 10−6 s−1) according to Hansen et al. (2007)
Full size table

The model used is the Block Adaptive-Tree Solar wind Roe-type Upwind Scheme (BATS-R-US) code (Powell et al. 1999), a component of the Space Weather Modeling Framework (Tóth et al. 2005, 2012). The multifluid MHD equations are solved on a block adaptive mesh that allows resolving the basic structures of the comet–solar wind interaction, from the extended mass-loading region spanning several million kilometers down to the size of the nucleus of just a few kilometers.

The results are presented in coordinates centered at the comet such that the undisturbed solar wind magnetic field is located in the x–y plane, the solar wind flow direction is along the negative x-axis, and the convectional electric field points in the z-direction. The same computational mesh was used for all simulations. The x-axis extends from −16 × 106 to 16 × 106 km. The y-axis is in the ecliptic plane and the z-axis completes the right-handed system and extends, as the y-axis, from −8 × 106 to 8 × 106 km. The cell sizes range from 200 m in the near nucleus region to several 10,000 km far away from the comet. The comet itself is represented by a sphere with 2 km radius, the corresponding boundary conditions are chosen such that plasma outflow at the nucleus is artificially set to zero: the nucleus is not considered to be a major source of plasma. However, inflow of plasma is possible e.g. in case no diamagnetic cavity has formed and the solar wind has direct access to the surface of the nucleus.

The approximately 3 million cells are split in 5400 computational blocks with 512 cells each. Each block, independent of the physical size it represents in the simulation domain, requires more or less equal computational resources. Therefore the blocks can be efficiently distributed on the different processors for parallel execution. We used 192 processors with hyperthreading for a total of 384 parallel threads. Each block is advanced at its optimum time-step based on the local numerical stability condition. The model is also able to run fully time-dependent simulations as done in Rubin et al. (2012). In this investigation, however, the simulation results represent steady-state snapshots for the different heliocentric distances. The results presented in this paper were obtained with the local Rusanov or TVD Lax-Friedrich scheme (Rusanov 1961) with the minmod limiter. The point-implicit scheme involves a linear system of equations to be solved independently in each grid cell to avoid problems with stiff source terms. This matrix contains the partial derivatives of the source terms and is obtained numerically. As a result the numerical time-step of the simulation can be increased albeit at the cost of additional numerical computations.

The BATS-R-US MHD comet model used here, in both single-fluid and multifluid versions, has previously been validated by comparison to measurements obtained by spacecraft such as the European Space Agency’s (ESA) Giotto probe which flew by comet 1P/Halley in 1986 (Gombosi et al. 1996; Rubin et al. 2009, 2014a) or observations at comet 19P/Borrelly by NASA’s Deep Space 1 spacecraft (Jia et al. 2008).

3 Results and Discussion

A total of six simulations have been carried out with each run separated by a heliocentric distance of 0.1 AU (see Table 1 for details). Figure 1 shows an overview of the mass density of the solar wind ions, from the largest heliocentric distance of 1.8 AU (top two panels) to the smallest distance, 1.3 AU, representing the conditions at perihelion (bottom two panels). Since protons dominate the solar wind the mass density in amu/cm3 corresponds roughly to the number density of the solar wind particles when assuming 1 amu per proton. The left panels indicate the x–z plane, i.e. the plane containing the convective electric field. All overview panels share the same size range and can thus be compared directly. Also the color scale is the same for all the panels. It shows for instance the increasing upstream solar wind proton density for smaller heliocentric distances used as input in our model (Table 1). As in Rubin et al. (2014b), cometary heavy ions, dominated by water group ions, start gyrating due to the convective electric field while the solar wind ions are deflected in the opposite direction. The corresponding plasma streamlines of the solar wind are indicated in a subset of the plates. Closer to the Sun the magnetic field is stronger leading to smaller gyro radii. Furthermore, the increased mass loading leads to deceleration of the solar wind and thus follows a further decrease in the local gyro radius. The right panels show the x–y plane, i.e. the plane containing the undisturbed upstream magnetic field.

Fig. 1
figure 1

Cometary solar wind mass density for the six different heliocentric distances from 1.8 to 1.3 AU top to bottom. The left column shows the plane aligned such that the convectional electric field points in the minus z-direction. The column on the right shows the ecliptic plane containing the undisturbed upstream magnetic field. The comet is located at the origin of the coordinate system and the positive x-direction points towards the Sun. The plates show how the Mach cone evolves into a bow shock as the comet approaches the Sun. A subset of plates shows projections of the solar wind streamlines

Full size image

In the left panels the formation of the bow shock can be observed. At 1.8 AU the interaction is still dominated by the Mach cone pointing towards the positive z-axis. As the activity increases a bow shock first forms around 1.7 AU and its stand-off distance increases as the comet becomes more active [cf. Koenders et al. (2013)]. The bow shock is a consequence of the mass addition by photon and electron-impact ionization of neutrals in the extended gas coma of the comet (Cravens 1989). The panels on the right show the asymmetry introduced by the non-vanishing B-field component (B x ) along the solar wind flow direction. As discussed above, the newly born ions are accelerated by the motional electric field. The trajectory of these ions is cycloidal, resulting from the superposition of the gyration and the E × B drift. The corresponding velocity space distribution is a ring beam distribution, where the gyration speed of the ring depends on the bulk plasma velocity and the angle between the solar wind velocity and the magnetic field vectors (Coates et al. 1989; Neugebauer et al. 1989; Coates and Jones 2009; Gombosi 2015). Since we assume a Parker spiral in our boundary conditions, this angle changes with heliocentric distance. From the model results we derived the bow shock stand-off distances presented in Table 2. For the conditions listed in Table 1, our model shows the formation of the bow shock at around 1.7 AU. Subsequently with increasing activity at perihelion the bow shock stand-off distance upstream of the comet nucleus reaches 3000 km. Koenders et al. (2013) derived a somewhat smaller bow shock stand-off distance of 2000 km from their Hybrid simulations while Hansen et al. (2007) obtained 3600 km with their single-fluid MHD model At 1.8 AU as well as farther away from the Sun the interaction is dominated by the formation of the above-discussed Mach cone.

Table 2 Upstream standoff distances of the bow shock and the diamagnetic cavity boundary in the model
Full size table

Figure 2 shows the corresponding modeled mass density of the cometary heavy ions, again from 1.8 AU (top panels) to perihelion (bottom two panels). We assume an average pick-up ion mass of 17 amu (Table 1). Therefore the number density of cometary pick-up ions can be estimated by dividing the mass density by 17 amu. As the comet approaches the Sun the solar insolation on the nucleus increases and leads to enhanced activity of the comet. Combined with the higher neutral gas density in the coma, the photoionization scale length also decreases as the comet approaches the Sun (we use a 1/r 2 dependence in our inputs in Table 1). The size of the interaction region clearly increases from higher to lower heliocentric distances. As in Fig. 1, the dominant asymmetries can be found in the x–z plane, the plane containing the convectional electric field. In order that momentum is conserved, the cometary heavy ions are deflected opposite the solar wind ions as shown by the streamlines. The x–y plane which contains the undisturbed upstream magnetic field is much more symmetric throughout the modeled heliocentric distances. As discussed before, solar wind protons and cometary heavy ions are coupled to each other through the electric and magnetic fields and therefore the solar wind ion bulk flow velocity is affected in regions of increased mass loading, i.e. close to the nucleus, and hence deflected in the opposite direction from the pick-up ions due to conservation of momentum and energy. As the comet becomes more active, the asymmetry becomes less and less pronounced. This is a consequence of the smaller local gyro radii due to increasing magnetic field strength and lower plasma velocities caused by higher mass loading.

Fig. 2
figure 2

Overview of the cometary heavy ion mass density for the six different heliocentric distances from 1.8 to 1.3 AU top to bottom. The left column shows the plane aligned such that the convectional electric field points in the minus z-direction. The column on the right shows the ecliptic plane containing the undisturbed upstream magnetic field. The comet is located at the origin of the coordinate system and the positive x-direction points towards the Sun. The plates show the interaction region of the comet with the solar wind growing in size as the comet’s activity increases. A subset of plates shows projections of the cometary plasma streamlines

Full size image

Figure 3 shows the overview of the magnetosonic Mach number \(\left( {{{{\text{here}}\;M_{SW} = v_{ion} }/ {\sqrt {v_{A}^{2} + v_{S}^{2} } }}} \right)\) where we use the mass averaged velocity of the ions and the Alfvén \(v_{A} = B/\sqrt {\mu_{0} \rho }\) and sound \(v_{S} = \sqrt {\gamma p/\rho }\) speeds with ρ and p the total mass density and pressure, respectively, such that the results should be comparable to single fluid MHD results. The trans-magnetosonic points are indicated by the contour line at M SW  = 1. Also indicated for the perihelion case is the approximate location of the bow shock. The figure shows how the solar wind is gradually slowed down already way upstream of the bow shock; it is therefore a weak shock.

Fig. 3
figure 3

Overview of the magnetosonic Mach number for the six different heliocentric distances between 1.8 and 1.3 AU, top to bottom, respectively. The left column shows the plane aligned such that the convectional electric field points in the minus z-direction. The column on the right presents the ecliptic plane containing the undisturbed upstream magnetic field. Furthermore a contour line at MSW = 1 indicates the trans-magnetosonic point and the approximate location of the bow shock for the perihelion case is indicated by the red dash-dotted lines. The undisturbed solar wind flow direction would be along the negative x-direction

Full size image

Coates et al. (1997) analyzed data obtained during the Giotto flyby at comet 26P/Grigg-Skjellerup (hereafter 26P) and found M SW  = 1.7 close to the shock location. Giotto crossed the bow shock at the flanks and in our simulation, depending on the position, M SW ranges from 1.5 to 2.

Eventually Giotto passed 26P some 100–200 km behind the nucleus (Neubauer et al. 1993). The smallest value derived by Coates et al. (1997) compares quite well to the M SW  = 0.75 from our model on the night side, 200 km behind 67P. Of course this comparison is limited, the fly-by at 26P occurred closer to the Sun at 1.01 AU during stronger magnetic field conditions (Neubauer et al. 1993). The gas production rates on the other hand are expected to be similar but in combination this results in a larger stand-off distance of ~10,000 km at 26P. A more detailed comparison and discussion can be found in e.g. Benna and Mahaffy (2006).

Figure 4 presents a close-up of the solar wind density, i.e. within 200 km from the nucleus. The close-ups show all the same distances and color scales. First, at heliocentric distances farther than 1.7 AU from the Sun a strong asymmetry can be observed in the x–z plane. The solar wind is deflected towards the positive z direction as shown by the Mach cone in Fig. 1. This leads also to the deflection of the cometary heavy pick-up ion tail in the positive z direction as shown in Fig. 2. With the increase in cometary activity, more and more solar wind protons charge-exchange with the neutral gas: it could be shown that this effect is already noticeable at heliocentric distances well beyond 2 AU (Nilsson et al. 2015). Finally, with the further increasing cometary activity a diamagnetic cavity forms: solar wind ions are unable to penetrate the innermost region around the nucleus as the plasma dynamics become governed by the outflowing neutral gas through charge-exchange reactions. This is clearly shown by the absence of solar wind particles inside the diamagnetic cavity in our model. The orientation of the cavity varies too: once formed, roughly inside 1.7 AU, the cavity together with the solar wind are deflected towards the positive z direction. But similarly to the cometary ion distribution shown in Fig. 2, this deflection becomes less important as the activity increases. Table 2 lists the deflection of the diamagnetic cavity with respect to the undisturbed solar wind inflow direction. In the x–y plane some asymmetries can be observed, but in general, however, they are less prominent. According to our model, at perihelion the extent of the diamagnetic cavity ranges from the surface up to roughly 30 km in the sunward direction. Koenders et al. (2015) report a distance of approximately 25 km obtained with their Hybrid model and the single fluid MHD model presented by Hansen et al. (2007) yields 38 km.

Fig. 4
figure 4

Close-up of the solar wind plasma mass density for the six different heliocentric distances from 1.8 to 1.3 AU top to bottom. The left column shows the plane aligned such that the convectional electric field points in the minus z-direction. The column on the right shows the ecliptic plane containing the undisturbed upstream magnetic field. A white sphere of 2 km radius marks the physical nucleus of the comet. The undisturbed solar wind flow direction would be along the negative x-direction. The plates show how the diamagnetic cavity first forms and then increases in size while the tilt with respect to the undisturbed solar wind flow direction decreases

Full size image

Figure 5 shows the corresponding distribution of the cometary heavy ions. The ion tail forms in the direction of the deflected solar wind, i.e. towards the positive z-direction. Therefore the cometary heavy ions are also deflected in the same direction. The gyration of the cometary heavy ions starts in the opposite way, i.e. the negative z direction as shown in Fig. 2. Due to the slowdown of the solar wind and the increasing magnetic field, however, the gyration radius becomes smaller with decreasing heliocentric distance and does not govern the whole interaction such as in our previous study by Rubin et al. (2014b). Furthermore, the higher the local neutral gas density, the higher the probability for collisional interaction with the abundant neutrals. Inside the diamagnetic cavity the outflowing neutrals govern the ion motion. Even outside, the cometary pick-up ions are frequently interrupted in their gyration due to collisional interaction with the neutral gas coma [see e.g. Rubin et al. (2012)]. Starting from roughly inside 1.7 AU, the formation of the inner shock separating the magnetized and non-magnetized regions can be seen. Table 2 lists both the size as well as the deflection of the cavity from the anti-sunward direction. This inner shock separates the cold supersonic ions from the much hotter subsonic ions and marks the edge of the diamagnetic cavity. Inside the dynamics of the plasma is expanding together with the bulk of the neutrals due to abundant charge-exchange reactions with the outflowing gas coma as shown by the streamlines of the cometary heavy ions added to one of the plates for reference. The supersonic nature of the outflow prevents the penetration of solar wind magnetic field and correspondingly the magnetized solar wind plasma (cf. Fig. 4). Just outside the inner shock, roughly 10 km according to our model of 67P, the outflowing plasma merges with mass loaded solar wind from further upstream. This is the so-called contact surface. Outside of the inner shock the electron temperature is still low and therefore ion electron-recombination is an important process, more than the slow tailward transport of the plasma. In this recombination layer a sizeable portion of the ions and electrons are lost through dissociative ion–electron recombination (Gombosi et al. 1996). As mentioned earlier, our model also considers the physical nucleus that shuts off photoionization in the shadow of the comet. This is visible along the negative x-direction inside the cavity: the radially outflowing neutrals, which govern the ion motion, prevent the plasma from filling in the void that gets created behind the nucleus. Potentially, behind the nucleus still a minimum level of mass loading occurs due to impact ionization by diffused photoelectrons created elsewhere in the cavity.

Fig. 5
figure 5

Close-up of the cometary heavy ion mass density for the six different heliocentric distances from 1.8 to 1.3 AU top to bottom. The left column shows the plane aligned such that the convectional electric field points in the minus z-direction. The column on the right shows the ecliptic plane containing the undisturbed upstream magnetic field. A white sphere of 2 km radius marks the physical nucleus of the comet. The undisturbed solar wind flow direction would be along the negative x-direction. The plates show the formation of the diamagnetic cavity and the inner shock and its alignment with respect to the Sun-comet line. One plate shows the streamlines of cometary pick-up ions inside and outside the diamagnetic cavity for reference

Full size image

Figure 6 shows the magnetic field strength for both planes in the close-up range of ±200 km around the nucleus. The magnetic field lines are also shown in the x–y plane, which contains the undisturbed upstream magnetic field. Very prominent is the formation of the diamagnetic cavity within roughly 1.7 AU from the Sun. Inside the cavity, the magnetic field drops to very low values as the solar wind has no access to this innermost region. As the comet becomes more active the cavity increases in size. The solar wind magnetic field lines, however, are draped around the cometary obstacle and compressed on the sunward side of the cavity and therefore show a corresponding increase in the field strength. The rotation of the diamagnetic cavity discussed earlier is also well visible in this figure of the magnetic field strength. Similar results have been obtained by Koenders et al. (2015) using a hybrid kinetic plasma model.

Fig. 6
figure 6

Close-up of the magnetic field strength for the six different heliocentric distances between 1.8 and 1.3 AU, top to bottom, respectively. The left column shows the plane aligned such that the convectional electric field points in the minus z-direction. The column on the right presents the ecliptic plane containing the undisturbed upstream magnetic field. Furthermore projected field lines show the draping of the magnetic field around the cometary obstacle. A white sphere of 2 km radius marks the physical nucleus of the comet. The undisturbed solar wind flow direction would be along the negative x-direction. The plates show the formation of the diamagnetic cavity and its alignment with respect to the Sun-comet line

Full size image

4 Summary and Conclusions

We have shown how the plasma interaction of a comet of the type 67P changes in appearance as the comet approaches the Sun. At distances around 1.8 AU and farther away from the Sun the interaction is dominated by a Mach cone while closer to the Sun the increased mass-loading of the solar wind leads to the formation of a bow shock wave which increases in size (Koenders et al. 2013). In addition, the near-nucleus region shows a significant change in the interaction: at large heliocentric distances the solar wind has direct access to the comet’s surface while for higher neutral gas production rates a diamagnetic cavity forms, which prevents the magnetized solar wind from reaching the innermost region of the comet. The multifluid MHD model offers several improvements over conventional single-fluid MHD models. Some of the kinetic effects that were previously only accessible to computationally much more expensive hybrid models can be reproduced (Rubin et al. 2014b). This includes the formation of the Mach cone as well as the rotation of the diamagnetic cavity with respect to the undisturbed solar wind inflow direction due to the momentum exchange between the solar wind and pick-up plasma populations as shown by hybrid-type simulations (Koenders et al. 2015). In our earlier work the large gyro radii of the pick-up ions dominated the interaction of comet 67P (Rubin et al. 2014b). Gradually, as the comet approaches the Sun, the activity increases and so does the mass loading. Conservation of energy and momentum leads to a slow down of the solar wind and in combination with the increasing magnetic field strength the gyration radii become smaller and we arrive at the situation at a distance of 1.8 AU presented in this paper.

In summary we have now a detailed description of the varying plasma interaction of comet 67P around the Sun, from large distances where the interaction is dominated by the large gyro radii to distances much closer to the Sun where mass-loading processes influence the solar wind flow direction and velocity to the point when eventually a diamagnetic cavity devoid of any magnetized solar wind forms.