Modeled Interaction of Comet 67P/Churyumov-
Gerasimenko with the Solar Wind Inside 2 AU
M. Rubin
1
•
T. I. Gombosi
2
•
K. C. Hansen
2
•
W.-H. Ip
3,4
•
M. D. Kartalev
5
•
C. Koenders
6
•
G. To
´
th
2
Received: 14 May 2015 / Accepted: 10 August 2015 / Published online: 18 August 2015
Springer Science+Business Media Dordrecht 2015
Abstract Periodic comets move around the Sun on elliptical orbits. As such comet 67P/
Churyumov-Gerasimenko (hereafter 67P) spends a portion of time in the inner solar system
where it is exposed to increased solar insolation. Therefore given the change in heliocentric
distance, in case of 67P from aphelion at 5.68 AU to perihelion at *1.24 AU, the comet’s
activity—the production of neutral gas and dust—undergoes significant variations. As a
consequence, during the inbound portion, the mass loading of the solar wind increases and
extends to larger spatial scales. This paper investigates how this interaction changes the
character of the plasma environment of the comet by means of multifluid MHD simula-
tions. The multifluid MHD model is capable of separating the dynamics of the solar wind
ions and the pick-up ions created through photoionization and electron impact ionization in
the coma of the comet. We show how two of the major boundaries, the bow shock and the
diamagnetic cavity, form and develop as the comet moves through the inner solar system.
Likewise for 67P, although most likely shifted back in time with respect to perihelion
passage, this process is reversed on the outbound portion of the orbit. The presented model
herein is able to reproduce some of the key features previously only accessible to particle-
based models that take full account of the ions’ gyration. The results shown herein are in
decent agreement to these hybrid-type kinetic simulations.
& M. Rubin
martin.rubin@space.unibe.ch
1
Physikalisches Institut, University of Bern, Sidlerstrasse 5, 3012 Bern, Switzerland
2
Atmospheric, Oceanic and Space Sciences, University of Michigan, 2455 Hayward Street,
Ann Arbor 48105, MI, USA
3
Graduate Institute of Astronomy, National Central University, Chung-Li 32054, Taiwan
4
Space Science Institute, Macau University of Science and Technology, Macao, China
5
Geospace Modeling and Forecasting Center, Institute of Mechanics, Bulgarian Academy of
Sciences, 1113 Sofia, Bulgaria
6
Institut fu
¨
r Geophysik und extraterrestrische Physik, Technische Universita
¨
t Braunschweig,
Mendelssohnstr. 3, 38106 Brunswick, Germany
123
Earth Moon Planets (2015) 116:141–157
DOI 10.1007/s11038-015-9476-8
Keywords Comet 67P/Churyumov-Gerasimenko Comet–solar wind interaction
Comet Plasma MHD Multifluid MHD Numerical simulation Rosetta
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;Re
`
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 (Ha
¨
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 dis-
tance 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 for-
mation 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
142 M. Rubin et al.
123
(Cravens 1989; Gombosi et al. 1996), which has been observed by the Giotto magne-
tometer at comet 1P/Halley (Neubauer 1987). The outflowing plasma is cold and super-
sonic 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 indi-
vidual ion species. Furthermore the electron pressure is calculated independently to derive
a self-consistent electron temperature (IV). The magnetic field is obtained from the com-
bination of Faraday’s law of induction (V) and the generalized Ohm’s law (VI).
I :
oq
s
ot
þr q
s
u
s
ðÞ¼
dq
s
dt
II :
oq
s
u
s
ot
þr q
s
u
s
u
s
þ Ip
s
ðÞZ
s
e
q
s
m
s
E þ u
s
BðÞ¼
dq
s
u
s
dt
III :
op
s
ot
þ u
s
rðÞp
s
þ cp
s
ru
s
ðÞ¼
dp
s
dt
IV :
op
e
ot
þ u
e
rðÞp
e
þ cp
e
ru
e
ðÞ¼
dp
e
dt
V :
oB
ot
¼rE
VI : E ¼u
e
B
1
n
e
e
rp
e
ð1Þ
with q
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 c the adiabatic index for which we assumed 5/3, i.e. the ratio of the
specific heats c ¼ f þ2ðÞ=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, H
2
O
?
,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
¼
P
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
Modeled Interaction of Comet 67P/Churyumov-Gerasimenko with the… 143
123
Table 1 The boundary conditions for the 6 cases between 1.3 and 1.8 AU are scaled from nominal values at 1 AU (n
SW
= 10 cm
-3
, |B| = 7 nT, f
i
= 1.0 9 10
-6
s
-1
)
according to Hansen et al. (2007)
Distance from the Sun (AU) 1.3 1.4 1.5 1.6 1.7 1.8
Neutral gas production rate (s
-1
) 5.0 9 10
27
4.0 9 10
27
3.0 9 10
27
2.2 9 10
26
1.7 9 10
27
1.2 9 10
27
Photoionization frequency (s
-1
) 6.0 9 10
-7
5.1 9 10
-7
4.5 9 10
-7
3.9 9 10
-7
3.5 9 10
-7
3.1 9 10
-7
Solar wind density (cm
-3
) 6.00 5.13 4.48 3.94 3.49 3.11
Solar wind electron temp. (K) 1.32 9 10
5
1.28 9 10
5
1.23 9 10
5
1.20 9 10
5
1.16 9 10
5
1.13 9 10
5
Solar wind ion temp. (K) 7.92 9 10
4
7.68 9 10
4
7.38 9 10
4
7.20 9 10
4
6.96 9 10
4
6.78 9 10
4
IMF B
x
component (nT) -2.97 -2.60 -2.27 -2.00 -1.76 -1.56
IMF B
y
component (nT) -3.84 -3.61 -3.38 -3.17 -2.97 -2.79
For all cases the solar wind velocity, u
SW
= 400 km/s, the neutral water outflow velocity, u
n
= 800 m/s, the neutral gas temperature, T
n
= 50 K, and the mean cometary
pick-up ion mass, m = 17 amu, are assumed constant
144 M. Rubin et al.
123
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 10
15
m
3
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 (Ha
¨
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
¼
Q
n
4pu
n
rjj
2
e
rjjt
io
u
n
ð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 t
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.
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 (To
´
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 9 10
6
to
16 9 10
6
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 9 10
6
to 8 9 10
6
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
Modeled Interaction of Comet 67P/Churyumov-Gerasimenko with the… 145
123
