© ESO 2018

1. A global model

The plasma interaction between the solar wind and a cometary atmosphere (coma) offers a unique situation in the solar system. The absence of an intrinsic magnetic field, the typical small size of the nucleus, and its negligible gravity combined with its highly elliptical orbit result in an ever changing interaction, in which the coma continuously and completely escapes the comet, dragged away by the magnetised stream of solar particles. These properties also result in one of the largest obstacles to the solar wind in the solar system. At comet Halley, the first cometary ions were detected 7.8 million kilometres away from the nucleus by the Giotto probe (Johnstone et al. 1986), a distance comparable to the day-side extent of the Jovian magnetosphere.

A major advance to comprehend this interaction was proposed by Alfvén (1957), who emphasised the role of the solar wind magnetic field in the formation of the cometary tails. Biermann et al. (1967) proposed a model of the day-side of a comet atmosphere in a hydrodynamical description of the interaction. These and all the previous efforts were tackling the features of the comet’s head and tails that were visible from Earth, naturally directing the scientific interest towards strongly outgassing comets close to their perihelion. In the next two decades, space probes were leaving Earth targeting such active comets, and at Giacobini-Zinner and Halley, what had previously been invisible became visible: a whole set of plasma structures came within the reach of scientists (Grewing et al. 1988; Cowley 1987; Gombosi 2015).

Between 2014 and 2016, the Rosetta spacecraft cohabited for more than two years with its host body, comet 67P/Churyumov-Gerasimenko (67P), enabling for the first time observations at large heliocentric distances (>3 au). Scientists were given the opportunity to witness the early interaction between a young tenuous coma and the solar wind, far away from the Sun. The nature of such an interaction was entirely new. Indeed, whereas at previously visited comets the interaction region was much larger than the scale of the ion gyromotion, resulting in what one could call a “fluid comet”, for which the classical fluid treatment of the plasmas applies, at 67P and at large heliocentric distances, the ion transit timescale through the coma is shorter than its gyroperiod, resulting in a “kinetic comet” for which no analytical approach is available so far.

Using in situ measurements, the evolution of this interaction was followed carefully (Behar et al. 2016a; Nilsson et al. 2017) resulting in some surprising findings. Initially barely disturbed, the solar wind started displaying a peculiar behaviour as the nucleus was getting closer to the Sun. Its flow slowly diverged from the Sun-comet direction, to eventually be seen flowing almost back towards the Sun at speeds of hundreds of kilometres per second. Eventually, the flow vanished from the in situ measurements: a void of solar ions was formed around the nucleus, while no severe deceleration was observed (Behar et al. 2017). The same interaction – for large to intermediate heliocentric distances – was also tackled by several simulation efforts, using hybrid particle-in-cell models (Bagdonat & Motschmann 2002; Hansen et al. 2007; Koenders et al. 2016a, b; Behar et al. 2016a), as well as a fully kinetic model (Deca et al. 2017). At heliocentric distances down to less than 2 au, simulations result in a highly asymmetric plasma environment. In particular, the solar wind presents a structure of high ion density only seen in the hemisphere of the coma opposite to the direction of the upstream electric field, one of the typical signatures of such a kinetic comet. This structure is interpreted as a Mach cone by Bagdonat & Motschmann (2002), a term adopted in several of the cited simulation studies. However until now, this asymmetric structure was only found in numerical models with intricate physics, and generally lacks the physical interpretation which would elucidate the experimental results exposed above.

Our goal is to understand and analyse the mechanism which transfers momentum between the solar wind and the coma, leading to such a deflection of a barely decelerated flow, in the absence of collisions. Additionally, we aim at providing novel insight into the nature of this asymmetric solar wind structure, which may be considered as the very seed of a cometary magnetosphere.

The present attempt to model the interaction puts the emphasis on the role of the motional electric field by considering the parameter space region, in which currents orthogonal to the magnetic field and pressure gradients can be neglected. Under these conditions, any noticeable disturbances of the flow necessarily result from the integrated interaction with and through the smooth and extended obstacle. This is in contrast to the situation at more classical and massive solar system bodies with intrinsic magnetic fields, dense and limited atmospheres, or conductive cores inducing a magnetic feedback to the solar wind. There, contrasted plasma structures and boundaries are formed, such as bow shocks, magnetopauses, induced magnetospheric boundaries or ionopauses. The obstacle to the solar flow is therefore compact and localised. Similar plasma boundaries may also appear at comets. For example, weak bow shocks were observed at comets Halley, Giacobini-Zinner, Grigg-Skjellerup and Borrelly (Coates 2009). However, such boundaries are only formed close to the Sun, and even in that case, the neutral atmosphere extends further out (a weak bow shock was observed about a million kilometres away from comet Halley’s nucleus, after the detection of the first cometary ions Johnstone et al. 1986), and mass-loading (the addition of new-born cometary ions to the solar wind) takes place, whether boundaries are formed closer to the nucleus or not. Therefore the present model should be representative and relevant for the region beyond the potential bow shock which forms when a comet gets closer to the Sun.

To infer the global behaviour of a system, whenever possible, an analytical model may overcome intrinsic limitations of in situ data (one-point measurement, instrumental errors and limitations) and simulation data (simulation of only a finite region of space, intricate physics, numerical limitations). While doing so, it allows to encapsulate one or a few of the driving mechanisms of a system in a reduced form, though at the cost of realism. In the present series of articles, the synergy between these three approaches – experimental data, numerical simulations, and theoretical models – is explored. This article focuses on the physical model and provides an expression for the force experienced by single solar wind protons, through the extended coma. The corresponding dynamics is thoroughly solved by Saillenfest et al. (2018), a solution widely used in the present work. The semi-analytical model we propose is computationally very cheap, and allows for a systematic comparison to each and every in situ data point. This extended comparison, together with the comparisons to numerical simulations, follows in subsequent articles of the series.

The model developed in the following sections requires several sub-models. One is a description of the ionised coma and its density distribution, following the same need for simplicity-to-relevancy ratio. The second is a description of the electric field and the magnetic field, which piles up due to the local decrease in the average velocity of the electrons, as slow new-born cometary ions are added to the flow. The motional electric field is completely dependent on the motion of the particles, which itself depends on the electric and magnetic fields. This inter-dependency is tackled in the following section as a generalised gyromotion, and results in a three-dimensional (3D) model of its own.

2. Generalised gyromotion

In this section, we derive the general dynamics of the interaction between two collisionless beams of plasma that are only subject to the Lorentz force. The subscripts sw and com are used for parameters and values of the solar wind and the cometary ion populations, respectively. The characteristic length, time, and velocity of the system are ℓ, t and u. E and B are the electric and the magnetic field, u_i is the average velocity of all charges carried by ions, n is the plasma number density, e is the elementary charge, j is the electric current and P_e is the electron pressure.

Our starting point is the simplified Ohm’s law, in which the electron inertial and the resistivity/collisional terms are neglected (see, e.g. Valentini et al. 2007). The plasma is weakly magnetised and quasi-neutral. The system is considered to be at steady state: ∂_t ⋅ ≡0.

$$\begin{array}{c}\hfill \mathit{E}=-\underline{{\mathit{u}}_{\mathrm{i}}}\times \mathit{B}+\frac{1}{\mathit{ne}}\mathit{j}\times \mathit{B}-\frac{1}{\mathit{ne}}\mathbf{\nabla }{\mathit{P}}_{\mathrm{e}}.\end{array}$$(1)

The total electric field exhibits three distinct components: the motional electric field, the Hall term, and the pressure gradient term. Analysing orders of magnitude in these three terms, as well as in the Ampère’s law, the Faraday’s law, and the Lorentz force, the following orderings can be found.

$$\begin{array}{c}\hfill \begin{array}{c}\ell \gg d_{\mathrm{i}}\Rightarrow |\underline{{\mathit{u}}_{\mathrm{i}}}\times \mathit{B}|\gg 1/ne|\mathit{j}\times \mathit{B}|\hfill \\ {\ell }_{\mathrm{P}}\gg r_{\mathrm{ge}}{v}_{\mathrm{the}}/u\Rightarrow |\underline{{\mathit{u}}_{\mathrm{i}}}\times \mathit{B}|\gg 1/ne|\mathbf{\nabla }{\mathit{P}}_{\mathrm{e}}|\hfill \end{array}.\end{array}$$(2)

In these expressions, d_i is the ion inertial length, and ℓ_P is the characteristic length of the pressure gradient. r_ge  =  m_eu_e/(eB) is the electron gyroradius and ${\mathit{v}}_{\mathrm{the}}=\sqrt{2k{T}_{\mathrm{e}}/m_{\mathrm{e}}}$ the electron thermal speed. At these scales (Eq. (2)), the electric field is reduced to

$$\begin{array}{c}\hfill \mathit{E}=-\underline{{\mathit{u}}_{\mathrm{i}}}\times \mathit{B}.\end{array}$$(3)

This directly implies that the currents perpendicular to the magnetic field are negligible for scales ℓ  ≫  d_i. The average velocity u_i in the case of our two beams can be reduced to only two terms:

$$\begin{array}{c}\hfill \begin{array}{cc}\underline{{\mathit{u}}_{\mathrm{i}}}\hfill & ={\xi }_{\mathrm{sw}}\underline{{\mathit{u}}_{\mathrm{sw}}}+{\xi }_{\mathrm{com}}\underline{{\mathit{u}}_{\mathrm{com}}}\hfill \\ \hfill & \\ {\xi }_{\mathrm{sw}}\hfill & ={\displaystyle \frac{n_{\mathrm{sw}}{q}_{\mathrm{sw}}}{n_{\mathrm{sw}}{q}_{\mathrm{sw}}+n_{\mathrm{com}}{q}_{\mathrm{com}}}};{\xi }_{\mathrm{com}}={\displaystyle \frac{n_{\mathrm{com}}{q}_{\mathrm{com}}}{n_{\mathrm{sw}}{q}_{\mathrm{sw}}+n_{\mathrm{com}}{q}_{\mathrm{com}}}}·\hfill \end{array}\end{array}$$(4)

In the absence of any force other than the Lorentz force, the dynamics of a single particle in either of the two populations is described by the following system of ordinary differential equations:

$$\begin{array}{c}\hfill \begin{array}{ccc}{\dot{\mathit{u}}}_{\mathrm{sw}}\hfill & =\hfill & q_{\mathrm{sw}}/m_{\mathrm{sw}}(\mathit{E}+{\mathit{u}}_{\mathrm{sw}}\times \mathit{B}),\hfill \\ {\dot{\mathit{u}}}_{\mathrm{com}}\hfill & =\hfill & q_{\mathrm{com}}/m_{\mathrm{com}}(\mathit{E}+{\mathit{u}}_{\mathrm{com}}\times \mathit{B}).\hfill \end{array}\end{array}$$(5)

Considering two initially perfect beams in velocity space, we have u_com = u_com and u_sw = u_sw: all particles of the population experience the same acceleration at the same time. The temperature is not defined.

Without loss of generality, one can choose a frame in which the magnetic field is directed along the y-axis. With Eqs. (3)–(5), and q_sw  =  q_com  =  q:

$$\begin{array}{c}\hfill \begin{array}{ccc}{\dot{\mathit{u}}}_{\mathrm{sw}}\hfill & =\hfill & {\displaystyle \frac{q{\xi }_{\mathrm{com}}B}{m_{\mathrm{sw}}}}({\mathit{u}}_{\mathrm{sw}}-{\mathit{u}}_{\mathrm{com}})\times \widehat{\mathit{y}},\hfill \\ \hfill \\ {\dot{\mathit{u}}}_{\mathrm{com}}\hfill & =\hfill & -{\displaystyle \frac{q{\xi }_{\mathrm{sw}}B}{m_{\mathrm{com}}}}({\mathit{u}}_{\mathrm{sw}}-{\mathit{u}}_{\mathrm{com}})\times \widehat{\mathit{y}}.\hfill \end{array}\end{array}$$(6)

A first important result is that there can be no particle acceleration along the magnetic field.

A second noteworthy result is that the velocity of the centre of mass, defined as

$$\begin{array}{c}\hfill \underline{{\mathit{v}}_{\mathrm{i}}}=\frac{n_{\mathrm{sw}}{m}_{\mathrm{sw}}}{n_{\mathrm{sw}}{m}_{\mathrm{sw}}+n_{\mathrm{com}}{m}_{\mathrm{com}}}{\mathit{u}}_{\mathrm{sw}}+\frac{n_{\mathrm{com}}{m}_{\mathrm{com}}}{n_{\mathrm{sw}}{m}_{\mathrm{sw}}+n_{\mathrm{com}}{m}_{\mathrm{com}}}{\mathit{u}}_{\mathrm{com}}\end{array}$$(7)

is conserved through time, $\underline{{\dot{\mathit{v}}}_{\mathrm{i}}}=\mathbf{0}$. This holds over spatial scales shorter than ℓ and ℓ_P, and if no mass is added. The equations of motion have the general form $\dot{\mathit{u}}=\omega \mathit{u}\times \widehat{\mathit{y}}$. In velocity space, the two beams move along circles perpendicular to the magnetic field (no acceleration along the magnetic field), as shown in Fig. 1.

If the two beams have initially the same parallel velocity (velocity component along the magnetic field), the two circles are in one and the same plane, centred on v_i, independent of any change of inertial frame. This can easily be seen in Fig. 1. Then in the frame in which v_i = 0, the two populations describe circles in velocity and physical space. The generalised gyrofrequency and gyroradii are then:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathcal{R}}_{\mathrm{sw}}=|{\mathit{u}}_{\mathrm{sw}}|/\omega \hfill \\ \hfill & {\mathcal{R}}_{\mathrm{com}}=|{\mathit{u}}_{\mathrm{com}}|/\omega \hfill \\ \hfill & \\ \hfill & \omega =eB{\displaystyle \frac{n_{\mathrm{sw}}{m}_{\mathrm{sw}}+n_{\mathrm{com}}{m}_{\mathrm{com}}}{(n_{\mathrm{sw}}+n_{\mathrm{com}})m_{\mathrm{sw}}{m}_{\mathrm{com}}}}·\hfill \end{array}\end{array}$$(8)

Still considering the case where both populations have the same initial parallel velocity, one can also always choose an inertial Cartesian frame in which u_sw is along the x-axis, and u_com = 0. This frame is referred to as the comet frame, and is used in Fig. 2. Both beams describe circles in velocity space with the same angular speed and the same centre v_i. The corresponding motion of the ions in physical space is a trochoid, the most general two-dimensional (2D) gyromotion. In the comet frame, the particles belonging to the cometary population (com) describe a more classical cycloid, as they periodically reach a velocity equal to 0. In Fig. 2, one can see that as the density ratio increases, |v_i| becomes closer to the origin. This can be interpreted as the corresponding slowing down of the plasma fluid for spatial scales much larger than the generalised gyroradius.

We note that if neither the ions nor the electrons have a velocity component parallel to the magnetic field, and they will not gain such a component during the interaction. We therefore obtain u_i = u_e. More generally, this equality is only verified if at any point in time, electrons and ions have the same parallel velocity, independent of the reference frame.

The dynamics depends greatly on the density ratio. If n_sw ≫ n_com, then v_i ∼ u_sw (top-part in Fig. 2), and the seldom cometary ions behave as test particles in the almost undisturbed flow of the population sw. The cycloid has then a radius equal to the cometary ion Larmor radius. As the density ratio n_com/n_sw increases, particles of the cometary population still describe a cycloid, though the corresponding radius decreases. When the density ratio is equal to the inverse of the mass ratio, both populations move along cycloids of equal radius, as seen in the middle panel of Fig. 2.

If the beams do not have the same initial parallel velocity, independent of the choice of frame, at least one of the populations will drift along the y-axis at a constant speed. This is the case in the comet frame, when accounting for the Parker spiral angle: both populations present an additional drift perpendicular to the x-axis, with opposite directions. Flipping the sense of the magnetic field does not change the direction of the drift, and the average Parker spiral configuration of the interplanetary magnetic field induces a dawn-dusk asymmetry in the ion dynamics around the comet. This topic is tackled in Behar et al. (2018b), based on statistics over the entire mission, and is based on this generalised gyromotion.

This generalised gyromotion is a collisionless 3D description of such a beam-beam interaction, much broader than the following 2D application for a comet.

3. Generalised gyromotion in a cometary atmosphere

We will now see to what extent this generalised gyromotion can describe the dynamics of the solar wind during its interaction with a comet. In order to resolve the equation of motion for the solar wind protons in Eq. (6), we need to express three main parameters, namely the cometary ion density and velocity, and the magnetic field.

3.1. Cometary ion density

The spatial distribution of the cometary ions is a major ingredient of the model, as it defines what obstacle is presented to the solar wind. The cometary atmosphere is assumed to have a spherical symmetry. For this exercise, the size of the nucleus is negligible, and so is its mass: the neutral elements, produced at a rate Q, are expanding radially in all directions with a constant speed u₀. We assume these particles to be water molecules, H₂O. They are ionised or photo-dissociated at a rate ν_d. By writing the equation of continuity with source terms on the cometary neutral density n₀, we obtain, in this spherical symmetry,

$$\begin{array}{c}\hfill \frac{1}{r^2}\frac{\mathrm{d}(r^2{n}_0{u}_0)}{\mathrm{d}r}=-{\nu }_{\mathrm{d}}{n}_0,\end{array}$$(9)

with the following solution established and used by Haser (1957):

$$\begin{array}{c}\hfill n_0(r)=\frac{Q}{4\pi u_0{r}^2}\times {\mathrm{e}}^{-r/R_{\mathrm{d}}};R_{\mathrm{d}}=u_0/{\nu }_{\mathrm{d}}.\end{array}$$(10)

The cometary ions are created by ionisation of the neutral particles with a rate ν_i. They have the initial radial velocity u₀ but will immediately be accelerated by the local electric and magnetic fields, to eventually escape the region of the denser coma. We separate the ionised coma into two different cometary ion components: the new-born cometary ions first, which are the main obstacle to the solar wind, and second the accelerated (or pick-up) cometary ions. The dynamics of the first population is assumed to be trivial: the new-born cometary ions move radially away from the nucleus with the same speed as the neutral molecules. The dynamics of the second, however, is much more complex, driven mostly by the mass-loading mechanism, meaning that the pick-up cometary ions leave the system quicker than they would have ballistically. New born ions become pick-up ions at a rate of ν_ml. Accordingly, a destruction term appears in the continuity equation of the new-born cometary ions:

$$\begin{array}{c}\hfill \frac{1}{r^2}\frac{\mathrm{d}(r^2{n}_{\mathrm{com}}{u}_0)}{\mathrm{d}r}={\nu }_{\mathrm{i}}{n}_0-{\nu }_{\mathrm{ml}}{n}_{\mathrm{com}}.\end{array}$$(11)

With R_i = u₀/ν_i and R_ml = u₀/ν_ml

$$\begin{array}{c}\hfill n_{\mathrm{com}}(r)=\frac{1}{R_{\mathrm{i}}}\frac{R_{\mathrm{d}}{R}_{\mathrm{ml}}}{R_{\mathrm{d}}-R_{\mathrm{ml}}}(1-{\mathrm{e}}^{-r(\frac{1}{R_{\mathrm{ml}}}-\frac{1}{R_{\mathrm{d}}})})·n_0(r).\end{array}$$(12)

Three characteristic radii are found in the density profile. In Sect. 3.5, we see that R_ml  ≪  R_d <  R_i. For r <  R_ml, that is, before the new born ions are accelerated and neglected, n_com(r <  R_ml)∝1/r, a result observed by the Rosetta mission at comet 67P in the first ∼200 km from the nucleus, and discussed by Edberg et al. (2015). At the larger scales that we are interested in, R_ml ≪ ℓ ≪ R_d, the neutral and the ion densities are proportional, while the exponential term is still negligible. Subsequently, we obtain

$$\begin{array}{c}\hfill n_{\mathrm{com}}(R_{\mathrm{ml}}\ll r\ll R_{\mathrm{d}})=\frac{{\nu }_{\mathrm{i}}}{{\nu }_{\mathrm{ml}}}\frac{Q}{4\pi u_0}·\frac{1}{r^2}[{\mathrm{m}}^{-3}]·\end{array}$$(13)

In this description, we have a steady creation and disappearance of the slow, new-born cometary ions that are constituting the bulk mass of the ionised coma, which interacts electromagnetically with the solar wind.

3.2. Two-dimensional magnetic pile-up

Another important term to model in the equation of motion of the protons is the magnetic field B within the coma. We first need to express u_i, the total velocity of the ion fluid, in which the B-field is considered to remain frozen-in:

$$\begin{array}{c}\hfill \underline{{\mathit{u}}_{\mathrm{i}}}=\underline{{\mathit{u}}_{\mathrm{e}}}=\frac{n_{\mathrm{sw}}}{n_{\mathrm{sw}}+n_{\mathrm{com}}}{\mathit{u}}_{\mathrm{sw}}+\frac{n_{\mathrm{com}}}{n_{\mathrm{sw}}+n_{\mathrm{com}}}{\mathit{u}}_{\mathrm{com}}·\end{array}$$(14)

Our goal is to solve for u_sw, therefore necessarily some more assumptions have to be made in order to simplify the total ion velocity and remove degrees of freedom in the system.

We first define the Comet-Sun-Electric field reference frame, illustrated in Fig. 3, as follows: the upstream solar wind magnetic field is directed purely along the y-axis, with an amplitude B_∞. The corresponding electric field according to Eq. (3) is along the z-axis with an amplitude E_∞ = u_∞B_∞. In this precise frame (comet-centred CSE), u_com ≪ u_sw, and it can be shown with the help of Eq. (13) that in most of the interaction region, u_i ∼ ξ_swu_sw. The latter is in agreement with fully kinetic simulations of the interaction (e.g. Fig. 2 of Deca et al. 2017). In the generalised gyromotion, the asymmetry in the total ion velocity appears only because of the different masses of the two populations, as for identical masses u_i = v_i (see Sect. 2), that is, both not accelerated and remaining along the x-axis only.

As B_∞ is along the y-axis (i.e. no Parker spiral angle), u_∞ is along the x-axis and as the cometary outflow is spherical, the plane y  =  0 is a plane of symmetry of the system. Therefore, within y  =  0, neither v_i nor u_i can have a component along $\widehat{\mathit{y}}$ (corresponding to the fact that no particle acceleration happens along B). One more simplification is needed to be able to express the magnetic field. We assume that the total ion velocity remains along the x-axis and follows:

$$\begin{array}{c}\hfill {\mathit{u}}_{\mathrm{i}}=-{\xi }_{\mathrm{sw}}{u}_{\infty }\widehat{\mathit{x}},\end{array}$$(15)

with ξ_sw >  0 and u_∞ >  0. From Eq. (15) of the total ion velocity and Eq. (3) of the electric field, one finds that E  =  (0,−u_iB_z,u_iB_y) with u_i  =  |u_i|> 0. Additionally, the steady state Faraday’s law ∇  ×  E  =  0 states that neither E_y nor E_z vary along x. Finally, in the plane of interest (y  =  0):

$$\begin{array}{c}\hfill \mathit{B}=\frac{B_{\infty }}{{\xi }_{\mathrm{sw}}}\widehat{\mathit{y}}.\end{array}$$(16)

The magnetic field frozen in the ion fluid should – in the absence of a Hall term – depend on the motion of both the solar wind and cometary ions. Since the advection of the cometary ions is not solved, the total ion velocity cannot be consistently derived.

3.3. Solar wind proton dynamics

Considering the new-born cometary ion population, one can assume in the cometary frame that u_com  ≪  u_sw. Therefore, using Eqs. (6), (15) and (16), we have

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{\mathit{u}}}_{\mathrm{sw}}={\displaystyle \frac{e{\xi }_{\mathrm{com}}B}{m_{\mathrm{sw}}}}·{\mathit{u}}_{\mathrm{sw}}\times \widehat{\mathit{y}}\hfill \\ ={\displaystyle \frac{e{\nu }_{\mathrm{i}}Q{B}_{\infty }}{4\pi m_{\mathrm{sw}}{n}_{\mathrm{sw}}{\nu }_{\mathrm{ml}}{u}_0}}·{\displaystyle \frac{1}{r^2}}·{\mathit{u}}_{\mathrm{sw}}\times \widehat{\mathit{y}}[{\mathrm{m}s}^{-1}].\hfill \end{array}\end{array}$$(17)

The force experienced by an individual solar wind proton is therefore of the form

$$\begin{array}{c}\hfill \mathit{F}=\frac{m_{\mathrm{sw}}\eta }{r^2}{\mathit{u}}_{\mathrm{sw}}\times \widehat{\mathit{y}}.\end{array}$$(18)

The force is always perpendicular to the proton velocity, with a strength proportional to the inverse of the square distance 1/r². The equation of motion for protons is then

$$\begin{array}{c}\hfill \ddot{\mathit{r}}=\frac{\eta }{r^2}\dot{\mathit{r}}\times \widehat{\mathit{y}};\eta =\frac{e{\nu }_{\mathrm{i}}Q{B}_{\infty }}{4\pi {\nu }_{\mathrm{ml}}{n}_{\mathrm{sw}}{m}_{\mathrm{sw}}{u}_0}[{\mathrm{m}}^2{\mathrm{s}}^{-1}].\end{array}$$(19)

In this description, the solar wind protons do not lose energy and are only gyrating, with a gyroradius function of their distance to the nucleus only. This can also be seen as the motion of charged particles in an effective magnetic field always perpendicular to the plane of the motion, with an amplitude proportional to 1/r². This is the core of the model, the reduced form of the solar wind proton interaction with a coma.

The dynamical system defined by Eq. (19) for solar wind protons is integrable. Its solutions are thoroughly studied in Saillenfest et al. (2018); here, we recall their main features. Let us introduce the polar coordinates (r, θ) in the (x, z) plane of motion. The dynamical equations rewrite¹:

$$\begin{array}{cc}& \ddot{r}-r{\dot{\theta }}^2=-\frac{\eta }{r}\dot{\theta }\hfill \\ \hfill & r\ddot{\theta }+2\dot{r}\dot{\theta }=\frac{\eta }{r^2}\dot{r},\hfill \end{array}$$(20)

where the dot means derivative with respect to time t. These coupled differential equations imply the conservation of kinetic energy E and a generalised angular momentum C that can be expressed as two characteristic radii:

$$\begin{array}{cc}& r_E=\frac{|\eta |}{u}\text{where}\hspace{1em}u=\sqrt{{\dot{r}}^2+r^2{\dot{\theta }}^2}\hfill \\ \hfill & r_C=rexp(1-r^2\dot{\theta }/\eta )\hfill \end{array}$$(21)

Their respective values fully determine the trajectory, with a bifurcation occurring at r_C = r_E. As for 1/r³ magnetic fields (Graef & Kusaka 1938), the solution (θ, t) can be written as a function of r defined by an integral.

The solar wind particles can be considered as originating from infinity on initially parallel trajectories. Since they all have the same conserved velocity u  =  u_sw, the characteristic radius r_E acts only as a scaling parameter (whereas the particles of the solar wind span all the possible values of r_C). With this setting, Saillenfest et al. (2018) show that around the nucleus at the origin, a circular cavity totally free of particles is created with radius² r_cav  ≈  0.28 r_E.

The resulting trajectories are shown in Fig. 4. A portion of the incoming flux of particles is temporarily focussed along a very specific curve, defined as the crossing points of infinitely close neighbour trajectories. By analogy to light rays, we call it a “caustic”, resulting in an overdensity of particles. This caustic has a well-defined shape, which can be expressed as the root of the variational vector. It is plotted in Fig. 4 for different values of the scaling parameter r_E, which rescales according to the comet activity and the heliocentric distance. A similar overdensity curve can be observed in other contexts, such as a flux of charged particles in the equatorial plane of a magnetic dipole (Störmer 1930; Shaikhislamov et al. 2015), or the deflection of solar particles around the thin atmosphere of Pluto (McComas et al. 2008). This could indicate that analogous processes are at play. Further discussions regarding a general 1/rⁿ law are given by Saillenfest et al. (2018).

We note that the notion of “impact parameter”, z_∞  =  z(x  →  ∞),  has no clear meaning for a 1/r² effective magnetic field (it is infinite for every particle). We should therefore express the problem in another way: we simply deal here with a far-enough starting distance for the particles, such that we can safely assume that their trajectories are parallel.

3.4. Third dimension

Outside of the plane (y = 0), the magnetic field draping introduces an angle between the magnetic field and the normal to the Comet-Sun line. As previously mentioned, such an angle would result in an additional acceleration of the solar wind along the y-axis, and the dynamics out of the (y = 0)-plane would not be planar. Just as its pile-up, the draping of the magnetic field is also given by the bulk velocity of the ions, assuming the field is frozen in the ion flow. Therefore a generalisation of the model to the third spatial dimension would result in the same physics as this 2D approach in the (y = 0)-plane.

3.5. Parameters and scales

Table 1 gives the physical quantities in the factor η together with their evolution with the heliocentric distance R. At comet 67P, the creation rate of neutral volatiles Q was found to be asymmetric around perihelion, with a higher activity after perihelion. A multi-instrument analysis can be found in Hansen et al. (2016), resulting in an empirical analytical fit given in Table 1. When necessary, we use the notation Q_i for the pre-perihelion in-bound leg of the Rosetta mission, and Q_o for the post-perihelion out-bound leg. The drawback of this empirical model is a discontinuity in the value of Q at perihelion. The value of the destruction rate ν_d and the ionisation rate ν_i and their dependence on the heliocentric distance is taken from Crovisier (1989). The magnetic field function of the heliocentric distance can be found in Cravens (2004).

The rate ν_ml at which new-born cometary ions are considered to turn into accelerated pick-up ions is arbitrary, but plays an important role in various parameters and scales of the model. It should be much larger than ν_i and u₀/ℓ, with ℓ being the characteristic length scale of the system. A proxy of its value can be obtained from in situ data. In a mission overview, using the imaging spectrometer RPC-ICA onboard Rosetta, Nilsson et al. (2017) show that a large majority of cometary ions observed at the spacecraft are rather slow, typically below a few tenths of an electron Volt. Detaching from this cold population, accelerated cometary ions are observed with fluxes two to three orders of magnitude lower. In the case study of Behar et al. (2016b), which used data from the RPC-ICA instrument as well, most of the values needed for estimating ν_ml can be found, with the exception of the proton and the cold cometary ion densities. Both densities are given here as the integrated plasma moment of order zero, on the same day. The data were taken on 2014-11-28, a day representative of the general cometary ion spectrum during the mission. An integrated spectrum (differential flux function of the energy per charge of the ions) is given in Fig. 5. The probe was 2.88 au from the Sun, 30 km away from the nucleus on a terminator orbit.

The cold ions reach energies up to 45 eV (22 km s⁻¹), and have an average density of 16.7 cm⁻³. The spacecraft potential, observed to be negative, on average, and often below −10 eV (Odelstad et al. 2017), accelerates the positive ions towards the spacecraft, so that the ions collected by the instrument appear more energetic that they actually are in the cometary plasma. The peak in the flux of the cometary new-born ions at about 20 eV corresponds to much lower energies, closer to 0 eV: the population modelled by cometary ions at rest in the analytical model, disappearing with rate ν_ml. Unfortunately, no spacecraft potential measurement is available on that day. The solar wind protons have an average speed of 376 km s⁻¹ (734 eV), with an average density of 0.14 cm⁻³. One obtains a total ion fluid velocity of ξ_swu_sw  ∼  3.13 km s⁻¹. With an average magnetic field of about 15 nT as measured by RPC-MAG that day, the motional electric field is about 0.05 mV m⁻¹. Assuming a constant acceleration, we get a very coarse duration of 77 s for the cometary ions to reach 45 eV. In summary, after about 77 seconds, slow new-born cometary ions surpass the kinetic energy of 45 eV (22 km s⁻¹) and become “accelerated” cometary ions, or pick-up cometary ions that have a density of 0.28 cm⁻³. The corresponding rate ν_ml is 0.011 s⁻¹ (corresponding at this heliocentric distance to a value of r_E  =  77.5 km). Moreover, on that precise day, it is fair to neglect the accelerated cometary ions, 60 times less dense than the cold ions. For the rest of the study, we consider ν_ml to be on the order of 0.01 s⁻¹ in magnitude: new-born cometary ions are neglected after being accelerated for 100 s.

3.6. Particle-field feedback

The steady-state assumption together with the scales that are considered in the problem result in extremely simplified Faraday’s and Ampère’s laws. Close to the caustic, two beams are seen, one incident and one emerging from it. These two beams of identical speed but different direction have the effect of decreasing the local bulk speed. The magnetic field should also be affected, slightly increasing along this structure. In turn, particle trajectories will be corrected by this magnetic feedback. Particles and fields will affect each other until steady state is reached. The model however cannot go further than the third step – the bulk speed decrease – in the sequence shown in Fig. 6.

Yet another source of magnetic field pile-up is missing: pick-up cometary ions that gained energy from the interaction are neglected. Even though these ions are present everywhere, they will more significantly increase the total cometary ion density in the +z-hemisphere of the coma. The magnetic field will pile-up slightly more in this hemisphere and the proton deflection is expected to be somewhat higher for positive z values, close to the nucleus. An example of such effects can be found in the data and simulation analyses of Koenders et al. (2016a).

3.7. Consistency of the model and summary

We gather the main assumptions and orderings we have been working with below. We remind that ℓ is the characteristic length of this steady-state system, u is its characteristic velocity, and ℓ_P is the characteristic length of the electron pressure gradient.

• i

  ℓ ≫ d_i

• ii

  ℓ_P ≫ r_gev_the/u

• iii

  ∂_t ⋅ ≡0

• iv

  u_com = u_i and u_sw = u_sw

• v

  B ⊥ u_sw

• vi

  R_ml ≪ ℓ ≪ R_d <  R_i

• vii

  n_com ∼ n_new-born

• viii

  u_com ≪ u_sw

• ix

  $\underline{{\mathit{u}}_{\mathrm{i}}}={\xi }_{\mathrm{sw}}{u}_{\infty }\widehat{\mathit{x}}$

Since the phenomenon we model is the deflection of the solar wind, the characteristic length ℓ can be considered as the distance over which a particle is deflected by some amount. In other words, the characteristic length is a fraction of the radius of curvature of the trajectories, which in our case has the simple and convenient expression ρ  =   − u_∞r²/η. We choose the definition ℓ  ≡  ρ, corresponding to the distance over which a particle is deflected by 1 radian. The smallest radius being modelled is the radius of the cavity, r_cav  ∼  0.28 |η|/u_∞, and therefore the assumptions are not verified closer to the nucleus.

We now review these assumptions.

i ℓ  ≫  d_i – We compared the values of ℓ and $d_{\mathrm{i}}=c/e\sqrt{{ϵ}_0\ast m/n}$ depending on the heliocentric distance, at the cavity radius (the most constraining cometocentric distance). For large heliocentric distances and at the cavity, both terms become of the same order of magnitude. Farther away from the nucleus, the ordering is verified.

ii ℓ_P  ≫  r_gev_the/u – Assuming an isothermal and spherically symmetric coma, the typical length scale of the electronic pressure gradient is ℓ_P  =  P_e/∂_rP_e  =  r(n_com  +  n_sw)/2n_com. r_ge  =  m_ev_the/(Be) is the electron gyroradius, $v_{\mathrm{the}}=\sqrt{(2e{E}_{\mathrm{e}}/m_{\mathrm{e}})}$ is the electron thermal speed, with E_e  ∼  7.5 eV being the electron thermal energy taken from Eriksson et al. (2017). The thermal speed for both the solar wind electrons and the cometary electrons are of the same order of magnitude, which is greater than the average speed of each electron population. Similarly to the Hall term of the electric field, condition ii is not fulfilled at large heliocentric distances and close to the nucleus.

iii ∂_t ⋅   ≡  0 – Dynamic phenomenons are not modelled, and the system can only reach another state adiabatically, with changing upstream conditions.

iv u_com = u_com and u_sw = u_sw – This is one of the strongest assumptions made in the model, that actually allows one to consider protons as single particles: the trajectories can only be relevant if they are not crossing each other. In the flow presented in Fig. 4, as long as protons do not reach the caustic, neighbouring trajectories are never intersecting: only the density and the bulk velocity change, but the assumed beam distribution is not deformed. We note that this is still true on the caustic itself, where, by definition, all trajectories are aligned. However, immediately after the caustic, two beams of comparable density appear in velocity space. This is a very interesting situation that is not accounted for in the present model. After the caustic, the beam quickly looses density and its particles will experience electric and magnetic fields dictated by the “upstream” beam: they will have a general gyromotion which is not modelled here.

v B  ⊥  u_sw – For an upstream magnetic field along the y-axis, the symmetry of the system guarantees this configuration everywhere in the plane (y = 0). The average Parker spiral angle will break this symmetry, and the two populations will have an additional drift along the y-axis. In a purely parallel case, $\mathit{B}=B\widehat{\mathit{x}}$ and no solar wind deflection nor cometary ion acceleration can happen. At heliocentric distances above 1 au, B becomes closer to 90° from the y-axis on average, and projected in the plane of interest we expect the dynamics to be qualitatively identical to the one depicted here-above.

vi R_ml  ≪  ℓ  ≪  R_d <  R_i – In Sect. 3.5, we have found that R_ml  ∼  10² km, R_d  ∼  10⁶ km and R_i  ∼  10⁸ km, verifying R_ml  ≪  R_d <  R_i. For r  ∼  R_d and beyond, the cometary ion density is already negligible compared to the solar wind density. We remind that only the ratio of the densities matters in the dynamics, and therefore the change of slope in the density profile at these scales barely has any impact on the solar wind dynamics.

vii n_com  ∼  n_new-born – (Discussed with viii )

viii u_com  ≪  u_sw – The observations presented by Nilsson et al. (2017) show that over the mission duration and as seen by the spacecraft, these two assumptions are sound. The flux of pick-up cometary ions is always two to three orders of magnitude above the flux of low-energy cometary ions, a difference that cannot be evened out by the difference of speeds (see also Fig. 5). The solar wind energy was also never close to the cold cometary ions’ energy.

ix $\underline{{\mathit{u}}_{\mathrm{i}}}={\xi }_{\mathrm{sw}}{u}_{\infty }\widehat{\mathit{x}}$ – This is actually inconsistent with the generalised gyromotion, as the ion bulk velocity is seen to change its direction (precisely the interest of the model). But as of now, this simplification seems necessary, and allows us to reduce the proton dynamics to the simple object of Eq. (19). We also mentioned that for most cometocentric distances, u_i ∼ ξ_swu_sw, so as long as the deflection of the protons is not too high, the resulting pile-up is relevantly modelled (within the limits of our description of the coma density).

In summary, different aspects of the flow in Fig. 4 are to be considered lightly, and these are:

• Trajectories of protons some time after they passed the caustic are non-physical (see iv). The local density is however only slightly disturbed.

• The magnetic pile-up will be affected in areas where protons have experienced significant deflection (see ix), typically the region downstream of the caustic, that is, the lower-left quadrant of the top graph in Fig. 4. For these two reasons, the solar wind ion cavity should most likely not be circular.

• In the region closest to the nucleus at large heliocentric distances, electrons and ions are expected to decouple and pressure gradients will be at work (see i and ii), which is not accounted for by the model.

Finally, the entire subject of waves and instabilities has been voluntarily neglected in the pursuit of simplicity, and it is believable that close to the nucleus and close to the Sun, these phenomena begin to play a major role. We refer to the work of Sauer et al. (1996), in which the authors study magneto-acoustic waves propagating transverse to the magnetic field in the frame of the bi-ion fluid theory, in a similar system (the artificial comet experiment conducted by the AMPTE mission).

4. Rosetta data and self-consistent models in the literature

We now have a 2D model to describe the velocity of individual solar wind protons around a comet. The dynamics is governed by a simple force acting on protons that is always perpendicular to their velocity and has an amplitude proportional to 1/r². The resulting flow in Fig. 4 is highly asymmetric and is only scaled with varying heliocentric distances. We now briefly review a selection of studies and results that support the model, both from in situ data and from numerical simulations.

4.1. Rosetta

A force always orthogonal to the solar wind protons and proportional to 1/r² in strength was initially proposed as an empirical model in Behar et al. (2017) to account for high solar wind deflection close to the nucleus together with very low deceleration, which eventually lead to the creation of a solar wind ion cavity. The ion cavity and the diamagnetic cavity observed at 67P and reported by Goetz et al. (2016b, a) are different, the former being much larger than the latter. More detailed and physical differences are discussed by Sauer et al. (1994) and Behar et al. (2017). It was also pointed out that just before the expanding cavity passed the spacecraft position³, the deflection was focusing on a value of around 140°, with proton velocity distributions more stable for a time. This would correspond to the crossing of the caustic, at the vertical (x = y = 0 and z >  0 in the CSE frame, Fig. 3) of the nucleus.

The orbit of the spacecraft over the two years of active mission provided two opportunities to map the solar wind flow over an extended region. The first was a day-side excursion, which ended up being almost entirely within the solar wind ion cavity and therefore of little interest here. During the second excursion, which was conducted at lower activity and in the night-side of the coma, the spacecraft reached distances up to almost 1000 km, and the solar wind was observed during the entire excursion. In Behar et al. (2018a), it is shown that a combination of the spacecraft position and of the upstream electric field orientation results in the spacecraft being fairly close to the plane of the model (y = 0 in the CSE frame) during most of the excursion. Ion data present an excellent match with the modelled flow, especially in the +z_CSE-hemisphere. The value of ν_ml giving the best fit with the data is 0.01 s⁻¹ (with all other parameters taken from other studies), surprisingly close to the value found above, ν_ml = 0.011 s⁻¹. We note that two independent methods based on different ion populations give the same value of ν_ml.

Additionally, no significant deceleration in the night-side of the coma was seen to correlate with the deflection, itself observed from just a few degrees up to 70°. As described by Eq. (19) and by the overall mission analysis of Behar et al. (2017) and Nilsson et al. (2017), the solar wind mostly gives momentum to the coma, without significant loss of kinetic energy. The present work provides a physical explanation for this important observation, and shows under which assumptions the solar wind can indeed be deflected by any angle, with negligible loss of kinetic energy.

4.2. Numerical models

The analytical expressions of the generalised gyromotion were verified with the hybrid simulation results, found in Lidström (2017).

The solar wind deflection pattern and the corresponding caustic (an over-density structure in the solar wind) can be seen in the results of numerical simulations, in several publications. In the context of comet 67P, this curved over-density in the solar wind can clearly be seen in the results of Wiehle et al. (2011; Fig. 3a), Koenders et al. (2016b; Fig. 14a), Koenders et al. (2016a; Fig. 3b), Behar et al. (2016a; Fig. 6), and Deca et al. (2017; Fig. 4c). Many of these results also show the general deflection of the solar wind, in qualitative agreement with the present model.

Such an asymmetric density structure in the solar particle flow can also be spotted in the simulation of the plasma environment at other solar system bodies. A first example is the solar wind dynamics modelled by Delamere (2009) at Pluto; see their Fig. 4. There, for two neutral production rate cases, the flow is highly asymmetric and develops a similar structure along which proton trajectories intersect. An even more familiar result can be found in Kallio & Jarvinen (2012), with the simulation of the solar wind interaction with unmagnetised bodies like Mars or Venus. Figure 3 presents the effect of mass-loading on the solar wind, in a test run where the the body has no physical extent, and newborn ions are created according to a 1/r² law, with a total production rate of 10²⁶ s⁻¹. This is virtually the same system as treated here, and therefore the strong agreement between the flow line of Fig. 3b in Kallio & Jarvinen (2012) and the proton trajectories modelled here-above is natural.

However, all these results present only the bulk velocity of the flow, with the exception of Delamere (2009). This makes it impossible to judge how single particles behave in the structure itself. In Fig. 7, trajectories of single solar wind protons taken from a self-consistent numerical model are given as a first illustrative overview. The hybrid-FLASH model is a hybrid particle-in-cell solver⁴ developed by Holmström (2010, 2013), and used in the context of 67P in Behar et al. (2016a) and Lindkvist et al. (2018). The model of the comet was purposefully kept simple, with a spherically symmetrical outgassing (Haser model), two ion populations (solar wind protons and cometary water ions), no charge exchange between neutral particles and ions, and in this particular simulation, B_∞ is along the y-axis. The cell size is constant and equal to 50 km. The heliocentric distance is 2 au, the production rate is Q = 7.52 × 10²⁶ s⁻¹, the ionisation rate is ν_i = 1.63 × 10⁻⁷ s⁻¹, and the speed of neutral molecules is 700 ms⁻¹. The upstream magnetic field of the solar wind is along the y-axis, B = 2.53 nT, its upstream speed is u_∞ = 430 km s⁻¹ and its upstream density is n_∞ = 1.25 cm⁻³.

The value of ν_ml that gives the best match between the caustics, especially at large distances, is found to be ν_ml  =  0.025 s⁻¹, which is larger than the previous experimental estimates. We find that similarly to in the semi-analytical model, proton trajectories are slowly deflected, and intersect each other in the (−x, −z) quadrant, forming a caustic. Immediately after passing the caustic, the protons are now experiencing the electric field mostly dictated by the denser incoming beam, and therefore starts a more complex gyromotion not accounted for by the analytical model. These protons are accelerated upward and cross the caustic. In both the simulation and the analytical model, immediately below the caustic, the phase space distribution function of the solar wind protons presents two beams (which might be similar to the observations of two proton beams reported by Jones & Coates 1997 at comet Grigg-Skjellerup). One can also note that the presence of the pick-up ions, denser in the +z hemisphere, can be seen in the locally higher deflection of the solar wind protons. We do not discuss the situation for cometocentric distances below 500 km, as with too few cells, one cannot properly resolve the smaller scales in this inner region, where charge exchange is also expected to play a role. A noteworthy observation is that in the simulation, the finite size of the box leads to an underestimation of the deflection, as can be seen at the boundary 2000 km upstream of the nucleus. This is due to an injection of solar wind protons at the upstream boundary with an initial velocity along the x-axis, whereas in the analytical model, protons have already experienced a significant deflection at this cometocentric distance. This issue is pointed out and quantified in Saillenfest et al. (2018; cf. Sect. 3.2). It is also probably one of the reasons why the best-fit value of ν_ml is larger here, since it must compensate for the distortion induced by the finite size of the simulation box.

Downstream of the caustic, the proton density drops in both panels, and accordingly to Eqs. (3) and (4), so does the electric field. In turn, the newborn cometary ions will be less accelerated than immediately upstream from the caustic, and will accumulate: the caustic shields the newborn cometary ions, and in turn a discontinuity in their density is expected to form along the caustic. This issue goes beyond the scope of this article, and is left for further studies.

Finally, the nature of this structure may now be discussed under a new light. Bagdonat & Motschmann (2002) describe the structure (the caustic in our description) as one side of an asymmetric Mach cone formed by the front wave of propagating density and magnetic field disturbances, induced by the obstacle – the newborn cometary ions in this case – in the incident flow. The complete asymmetry of the cone is however not further discussed. It is extremely interesting to note that the present model does not consider the super-magnetosonic character of the solar wind, nor does it propagate any type of disturbance. Based on this 2D approach, this asymmetric structure is not formed by a propagating perturbation, and therefore is not a Mach cone. Furthermore, the cometary newborn ions are only indirectly the obstacle in this picture, as the over-density (the caustic) is formed purely by the geometry of the deflected flow. The solar wind is forming an obstacle to itself, an obstacle with a shape independent of the magnetosonic Mach number. Developing the model to a third dimension, comparing it more thoroughly to self-consistent numerical models, and studying the effect of the plasma pressure on this structure is however necessary to conclude on this topic.

5. Concluding remarks

We have shown how momentum and energy are transferred between two collisionless plasma beams for spatial scales that are large compared to the ion inertial length and in the case of negligible electron pressure gradients. This so-called generalised gyromotion applies to the most general, or arbitrary, 3D configuration of two plasma beams.

There are two possible ways to consider the exchange of energy between the solar wind and the coma. The first way states that at scales that are large compared to the gyroradius, and based on classical fluid concepts, there is necessarily a loss of kinetic energy in the solar wind. The second way considers individual particles of the solar wind, and the present model shows how and under which conditions and assumptions these particles do not lose kinetic energy, at zero order. However, part of the plasma is accelerated in this interaction (pick-up ions); therefore at a higher order, also the solar wind loses kinetic energy. We further note that translated into bulk properties, this model also displays a deceleration of the fluid.

In the plane of symmetry of the classical magnetic field draping at comets, the exchange of momentum through the fields between the solar wind and a comet atmosphere results in a very simple expression of the force applied to the protons: this force is perpendicular to their velocity, with an amplitude proportional to 1/r². The solar particle flow is reduced to a peculiar and highly asymmetric pattern, exhibiting a caustic, which is also seen in numerical models. In situ data from the Rosetta mission show strong support to the semi-analytical model, in terms of deflection and speed, and with the observation of a solar wind ion cavity. We note that these results can be straightforwardly extended to solar wind alpha particles.

The cheapness of the model will allow for an extended and systematic comparison with in situ data, allowing us to distinguish how dominant the motional electric field is during a given activity level. The model may also point to previously unnoticed structures, such as the peculiar distribution function of the solar wind proton close to the caustic. Furthermore, the model should also greatly ease the understanding of complex mass-loaded solar wind kinetic simulations.

Finally, the validity of the model close to the nucleus is expected to crumble closer to the Sun, as a bow shock might form, at least for strong-enough cometary outgassing. There, waves and instabilities will provide additional ways to transfer energy and momentum, and these phenomena will most likely act to transform the described caustic into a bow shock, beyond the domain of validity of this model. Pinning down the conditions for which this transition happens, together with the micro-physics involved, is an obvious direction to explore. However, further out than the potential bow shock, at comets and at unmagnetised bodies as well, the mass-loading is affecting the solar wind flow, and the present model remains relevant.

---

Acknowledgments

This work was supported by the Swedish National Space Board (SNSB) through grants 108/12, 112/13, 96/15 and 201/15. The work at LPC2E/CNRS was supported by ESEP, CNES and by ANR under the financial agreement ANR-15-CE31-0009-01. This work was supported in part by NASA’s Solar System Exploration Research Virtual Institute (SSERVI): Institute for Modelling Plasmas, Atmosphere, and Cosmic Dust (IMPACT). Partial support is also acknowledged by Contract No. JPL-1502225 at the University of Colorado from Rosetta, which is an European Space Agency (ESA) mission with contributions from its member states and NASA. The hybrid solver is part of the openly available FLASH code and can be downloaded from http://flash.uchicago.edu/ developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. The hybrid simulations were conducted using resources provided by the Swedish National Infrastructure for Computing (SNIC) at the High Performance Computing Center North (HPC2N), Umeå University, Sweden. Part of this work was inspired by discussions within International Teams 336: “Plasma Surface Interactions with Airless Bodies in Space and the Laboratory” and 402: “Plasma Environment of Comet 67P after Rosetta” at the International Space Science Institute, Bern, Switzerland.

References

All Tables

All Figures