Research

^{Image: NAC_2016-01-27T16.20.08.963Z_ID10_1397549000_F22; Credit: ESA/Rosetta/MPS for OSIRIS Team MPS/UPD/LAM/IAA/SSO/INTA/UPM/DASP/IDA}

Cometary gas and dust outflow

In recent years I have focussed my research on building a self-consistent simulation pipeline to model the gas and dust emission from cometary nuclei taking into account the relevant physical processes within the inner comae (~ tens of kilometres above the cometary surface). Although the described methods can be applied to any comet we have focussed on understanding multi-instrument results from the European Space Agency’s (ESA) Rosetta mission to comet 67P/Churyumov–Gerasimenko (hereafter 67P). The results using this modelling approach have been published in Marschall et al. (2016), Marschall et al. (2017), Gerig et al. (2018), Marschall et al. (2020a), and Marschall et al. (2020b). On this page, you can find a brief overview of the different steps involved in our modelling pipeline.

The physics of the outflow from cometary surfaces has been the subject of intense study for more than 40 years. Ice sublimating into vacuum forms a non-equilibrium boundary layer, the “Knudsen layer” (Kn-layer), with a scale height of ~20 mean free paths (MFP). If the production rate is low, the Kn-layer becomes infinitely thick and the velocity distribution function (VDF) remains strongly non-Maxwellian. Nowadays, the state of the art method to study gas flows inside non-equilibrium regions is Direct Simulation Monte Carlo (DSMC). After obtaining a steady-state gas solution we introduce dust particles into the flow to determine the physical parameters of the dust coma (i.e. dust number densities).

1. THE SHAPE MODEL

The starting point of all our simulations is a shape model of the comet. We can use arbitrary shapes as input as long as the surface is closed. In the case of 67P, we use the SHAP7 model by Preusker et al. (2017) which is shown in the animation. These shape models are built of triangles, which we refer to as surface facets. The full resolution shape model of SHAP7 has more than 40 million facets which is far greater than we can use in our simulation due to computational costs. We thus use a reduced model with around 400’000 facets instead.

2. THE ILLUMINATION MODEL

Solar insolation is the main driver of cometary activity by providing the energy to sublimate ices at or close to the surface. Before we can calculate the gas production rate from the cometary surface we first need to determine the heat input to it. The two things having the largest impact on the energy available at the surface are the distance to the Sun and the angle between the normal of the local surface and the direction to the Sun, the so-called incidence angle, $i$ . For a specific point in time, we can calculate the position of the Sun relative to the comet. Or illumination model then calculates the incidence angles for every facet taking into account shadows cast by other parts of the cometary shape onto itself. The animation on the left shows one of the results that can be obtained. It can be clearly seen how the illumination conditions including the shadows change as the nucleus rotates. The Sun in this animation is located toward the top centre of the frame.

3. THE THERMAL MODEL

Now that we have calculated the incidence angle of the Sun light to the surface of the comet we can now use a thermal model to calculate the surface temperature and gas production rate from each of the surface facets. We have decided to use a simple thermal model omitting thermal conductivity (i.e. the thermal inertia was set to zero. We note that Gulkis et al (2015), Schloerb et al. (2015), and Choukroun et al. (2015) measured low values consistent with values seen at 9P/Tempel1) but including sublimation of water ice. We have thus imposed a simple thermal balance for each surface facet: $0 = \frac{S \left(1-A_H \right) \cos\left( i \right)}{{R_h}^2} - \epsilon \sigma T^4 - L Q_{gas} \quad,$ where $A_H$ is the directional–hemispheric albedo (set to 0.04, Fornasier et al. (2015)), $S$ the solar constant at 1 AU (taken to be 1384 W^-2), $i$ the angle of incidence, $R_h$ the heliocentric distance of the comet in AU, $\epsilon$ the IR emissivity (set to 0.9), $\sigma$ Stefan-Boltzmann’s constant, $L$ the latent heat of sublimation of water ice (2.84 MJ kg^-1; Huebner 2006), and $Q_{gas}$ the sublimation rate (i.e. the gas mass loss rate per unit time and unit area). The animation on the right shows how the temperature at the surface and the local production rate of each facet changes as the comet rotates.

4. THE GAS DYNAMICS MODEL

Now that we know which facet emits how many gas particles at which temperature we can now apply a gas dynamics calculation to determine where the molecules flow to in 3D space. Because the production rate at comet 67P spans large ranges from low to high the gas flow in principle transitions from either fluid or collisional to a free adiabatic outflow. Thus our preferred method is Direct Simulation Monte Carlo (DSMC) because it intrinsically covers all of these regimes. For high production rates, a fluid approach might be computationally more economical though. The code we are using, named PDSC++ (Su (2013)), has previously been used to model the water vapour distribution in the vicinity of comet 9P/Tempel 1 (Finklenburg et al. (2014)) and is based on the PDSC code developed by Wu and co-workers (Wu et al (2003, 2004, 2005)). PDSC++ allows the simulation of 2D, 2D-axisymmetric, and 3D flows on hybrid unstructured grids. The code has been parallelised, allowing a much larger number of cells, and has been implemented on several clusters in Bern (Switzerland) and Taiwan. The code is especially useful because it is able to treat the large density gradients by the implementation of a variable time step and a transient adaptive sub-cell technique to increase computational speed and accuracy in the regions of high density (Finklenburg et al. (2014), Finklenburg (2014)).

Within the DSMC code, the real gas molecules are represented by a smaller number of simulation particles that are used to generate a statistical representation of the flow. The microscopic behaviour of the simulation molecules is separated into a translational step and a collision step. In the collision step, the collision pairs are chosen randomly from all the molecules that are within one cell at the time of the collision step. Macroscopic gas properties such as number density, velocity, and temperature are then calculated by averaging the appropriate quantity over all simulation molecules in the respective sampling cell.

Even though the code is fully capable to treat multi-species flows, in this work we consider water only because it is the dominant species in the coma of 67P (Haessig et al. (2015)). In addition, we are interested in the dynamics within the first kilometres of the gas coma above the cometary surface and thus do not need to include any chemistry into our gas model. Processes such as photo-dissociation and ionisation occur on much longer time scales.

The results of this model provide gas number densities, temperatures, and speeds in 3D space around to comet out to the edge of our simulation domain (usually 10-20 km from the nucleus centre). An example result can be seen in the figure illustrating two slices through the 3D gas field showing the gas speed [m s-1]. The nucleus shows the illumination of the surface at the time of this simulation. The results from these simulations can be compared to several instruments on-board Rosetta incl. the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA), the Microwave Instrument for the Rosetta Orbiter (MIRO), and the Visible and Infrared Thermal Imaging Spectrometer (VIRTIS). By comparing simulation results and measurements we are able to constrain the gas emission regions at the surface and the physical processed occurring during the expansion of the gas from the surface into the vacuum. Furthermore, these kinds of simulations also provide a powerful tool to determine the limits of data interpretation possible.

5. THE DUST DYNAMICS MODEL

Given a gas flow field, we can now study the motion of dust particles in this flow. The motion of the dust is governed by these two dominant forces: the force produced by the gas drag on the dust particles to accelerate them away from the surface, and the gravity of the nucleus as an opposing force. The dust follows, though in a complex way, the gas flow from the nucleus. If the gas flow is perpendicular to the surface so is the motion of the dust particles within the first meters above the surface. But lateral flows in areas with large gas production rate gradients (such as e.g. close to the terminator) are possible and thus also result in lateral dust flows close to the surface. We treat the dust grains as test particles and consider them to be spherical, although there are strong indications that especially the larger dust grains are porous and fluffy aggregates (Kolokolova et al. (2010), Schulz et al. (2015), Rotundi et al. (2015), Langevin et al. (2016), Bentley et al. (2016), Mannel et al. (2016)). We justify our choice of treating the dust as spherical particles with the fact that the main parameters influencing the dust trajectory are the mass and cross-section. Therefore our dust particles must be understood as effective spheres with the respective cross-section and mass of the dust particles they represent. At this point, we still know little about the physical properties of the dust particle (such as their shapes, density) and, perhaps more importantly, about how they orient themselves in the gas flow (Ivanovski et al. (2017).
The equation of motion for each dust grain which has mass $m_d$ , radius $r_d$ , and geometric cross-section $\sigma_d = \pi r_d^2$ at position $\vec{x}$ and with velocity $\vec{v}_d = \frac{d\vec{x}}{dt}$ is
$m_d \frac{d^2 \vec x}{dt^2} = \vec F_G + \vec F_D$
$= \vec F_G + \frac{1}{2} C_D m_g n_g \sigma_d \left| \vec v_g - \vec v_d \right| \left( \vec v_g - \vec v_d \right) \quad ,$
where $\vec F_G$ is the gravitational force, $m_g$ the mass of the gas molecule considered (in our case molecular water), and $n_g$ and $\vec v_g$ are the number density and velocity of the gas. If we assume an equilibrium gas flow where the mean free path of the gas is much larger than the dust size (free molecular aerodynamics), the drag coefficient $C_D$ is defined as $C_D = \frac{2 \zeta^2 + 1}{\sqrt{\pi} \zeta^3} e^{-\zeta^2} + \frac{4 \zeta^4 + 4 \zeta^2 - 1}{2 \zeta^4} \text{erf}(\zeta) + \frac{2 \left( 1 - \varepsilon \right) \sqrt{\pi}}{3 \zeta} \sqrt{\frac{T_d}{T_g}}$ , with the gas temperature $T_g$ , the dust temperature $T_d$ (chosen to be equal to $T_g$ ), and $\varepsilon$ is the fraction of specular reflection (set to 0), and $\zeta = \frac{ \left| \vec v_g - \vec v_d \right| }{\sqrt{\frac{2kT_g}{m_g}}}$ .

By injecting millions of dust particles of different sizes into the gas flow and solving their equation of motion to determine their trajectories we can determine the local properties of the dust flow within the inner coma (i.e. dust number densities). How dust particles travel through space is illustrated in the animation on the right. These results can be compared to several instruments on-board Rosetta incl. the Optical, Spectroscopic, and Infrared Remote Imaging System (OSIRIS), the Grain Impact Analyser and Dust Accumulator (GIADA), the Cometary Secondary Ion Mass Analyser (COSIMA), and the Micro-Imaging Dust Analysis System (MIDAS).