The Dust-to-Gas Ratio, Size Distribution, and Dust Fall-Back Fraction of Comet 67P/Churyumov-Gerasimenko: Inferences From Linking the Optical and Dynamical Properties of the Inner Comae

In this work, we present results that simultaneously constrain the dust size distribution, dust-to-gas ratio, fraction of dust re-deposition, and total mass production rates for comet 67P/Churyumov-Gerasimenko. We use a 3D Direct Simulation Monte Carlo (DSMC) gas dynamics code to simulate the inner gas coma of the comet for the duration of the Rosetta mission. The gas model is constrained by ROSINA/COPS data. Further, we simulate for different epochs the inner dust coma using a 3D dust dynamics code including gas drag and the nucleus' gravity. Using advanced dust scattering properties these results are used to produce synthetic images that can be compared to the OSIRIS data set. These simulations allow us to constrain the properties of the dust coma and the total gas and dust production rates. We determined a total volatile mass loss of (6.1 ± 1.5) · 10⁹ kg during the 2015 apparition. Further, we found that power-laws with $q=3.7_{-0.078}^{+0.57}$ are consistent with the data. This results in a total of $5.1_{-4.9}^{+6.0}·10^9$ kg of dust being ejected from the nucleus surface, of which $4.4_{-4.2}^{+4.9}·10^9$ kg escape to space and $6.8_{-6.8}^{+11}·10^8$ kg (or an equivalent of $14_{-14}^{+22}$ cm over the smooth regions) is re-deposited on the surface. This leads to a dust-to-gas ratio of $0.73_{-0.70}^{+1.3}$ for the escaping material and $0.84_{-0.81}^{+1.6}$ for the ejected material. We have further found that the smallest dust size must be strictly smaller than ~30μm and nominally even smaller than ~12μm.

1. Introduction

The European Space Agency's (ESA) Rosetta mission escorted comet 67P/Churyumov-Gerasimenko (hereafter 67P) from August 2014 to September 2016 along its orbit through the inner Solar System. It watched as the comet's activity started to develop at large heliocentric distances, come to its culmination at perihelion, and decline as the comet traveled out toward Jupiter's orbit. This long-term continuous monitoring of the comet's activity has provided an unprecedented wealth of data on this comet and its activity.

The observations revealed a complex bi-lobate shape [1, 2] and diverse morphology [3]. As a comet approaches the Sun it is heated and the ices start sublimating and ripping with them dust particles. Thus one of the important questions to be answered was what the bulk of the comet was made of i.e., what the bulk refractory-to-volatile ratio is. In the simplified view where any ejected material is lost to space two measurements are sufficient to determine this ratio. First, the total mass loss during one apparition measured by the Radio Science Investigation (RSI) [4]. Second, the total volatile mass loss which can be indirectly determined by the in-situ measurements of the gas density (e.g., [5–7]) or remote sensing data (e.g., [8–11]). In this simple case, the refractory-to-volatile ratio can be immediately inferred from those two measurements. But the complex surface morphology has revealed large dust deposits [12] that indicate that possibly a large fraction of the ejected dust is re-deposited [13]. If that is indeed the case, then the two above mentioned quantities cannot constrain the total dust mass ejected but rather only the dust mass escaping the nucleus gravity. Further, the process of dust fall-back obscures the emitted dust-to-gas ratio.

One way of constraining the amount of fall-back material would be to attempt to measure the actual change in elevation of the surface as a function of time from local or global digital terrain models (DTM). We cannot assess at this point if that is indeed feasible with the Optical, Spectroscopic and Infrared Remote Imaging System (OSIRIS, [14]) data set. Another way is to couple the scattering properties of the dust with a dynamical model of the dust coma constrained by the brightness of the dust coma. In this work, we have adopted the latter approach and modeled the inner gas and dust comae for the entire Rosetta mission. We use Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA, [15]) data to constrain the gas production rate and OSIRIS data for our dust models. To constrain the dust models we compare the dust coma brightness as measured by OSIRIS to synthetic model images. This process links several dust parameters that are otherwise not easily combined. In particular, we will show how the dust size distribution, the dust-to-gas ratio, the fraction of fall-back and the optical properties are inter-dependant and thus cannot be determined independently.

In section 2, we will describe the method used and lay out the assumptions we have made. Furthermore, we will point out the free parameters of the models, that need constraining through Rosetta data. Some theoretical considerations are presented in section 3. We will discuss the results of our work in section 4 and summarize and conclude our work in section 5.

2. Method

In this work, we have used the modeling approach (and in particular our DRAG3D model for the dust coma) described in detail in Marschall et al. [16]. This approach has been successfully applied for the analysis and interpretation of multiple Rosetta instruments, in particular ROSINA, MIRO (Microwave Instrument for the Rosetta Orbiter), VIRTIS (Visible and Infrared Thermal Imaging Spectrometer), and OSIRIS [16–19]. While in previous work we have applied this approach to specific epochs of the Rosetta mission, we have employed it here to cover the entire mission period to study longer-term processes.

In the following, we will briefly repeat some of the most important parts of the modeling elements and refer to Marschall et al. [16] for a detailed description.

2.1. General Assumptions

The calculation of the 3D gas flow field using the Direct Simulation Monte Carlo (DSMC) method is very computationally expensive and it is therefore currently not feasible to cover the entire escort phase of Rosetta (from August 2014 to September 2016) with a high temporal resolution. It is thus necessary to split the comet's orbit into a number of epochs that are computationally feasible and then interpolate between the results using a linear scaling between epochs. To ensure that the calculated results are representative of the respective epoch we make sure that during each of the epochs neither the total solar energy reaching the surface nor where the energy strikes the surface changes substantially. The amount of energy deposited is driven primarily by the heliocentric distance, R_h, while the location of deposition apart from the rotation of the comet is controlled by the sub-solar latitude, LAT. We thus chose that the inverse square of the heliocentric distance of the comet's location at the start and end time of each epoch shall be within 15% of the location at the center date of each epoch. Furthermore that the difference in sub-solar latitude be less than 5° from the center time of epoch to the start and end of the epoch, respectively. This leads to the 20 epochs listed in Table 1 and illustrated in Figure 1. Simulations were run for the center time of each epoch. This choice also ensures that we cover the exact dates of the in- and outbound equinox (epochs 6 and 18) as well as perihelion (epoch 11) and summer solstice (epoch 12).

TABLE 1

Table 1. Start, center, and end time for each of the epochs as well as the heliocentric distance and sub-solar latitude of the center time of each epoch.

FIGURE 1

Figure 1. Heliocentric distance vs. sub-solar latitude of comet 67P during the escort phase of Rosetta as well as epochs used in this work.

The basis of all simulations is the 3D shape model by Preusker et al. [2]. We use a decimated model with ~440′000 facets due to our computational constraints. To fully define the illumination condition we need to select the sub-solar longitude in addition to the heliocentric distance and sub-solar latitude which are set by the choice of epoch. For each epoch we have run simulations for sub-solar longitudes of 0°, 30°, 60°, 90°, 120°, 150°, 180°, 210°, 240°, 270°, 300°, and 330°. This results in a total of 240 different illumination conditions for the entire mission period.

For each illumination condition, we calculate the incidence angle (angle between the surface normal and the direction of the Sun) of each facet taking into account self-shadowing. This allows calculating the solar energy entering the surface neglecting re-radiation from other facets. By means of a simple energy balance of the incoming solar energy, thermal re-radiation and sublimation we can calculate the sublimation temperature and the sublimation rate of each facet assuming pure water ice.

We do not take into account any emission from shadowed facets, be it due to local night or mutual shadowing by other parts of the nucleus. The calculated pure ice sublimation rate of each facet needs to be scaled to match observed sublimation rates at 67P. Here we assume a pure H₂O ice surface that is areally mixed with inert refractory surface akin to a checkerboard pattern. This surface fraction of the facet covered by ice, which is a priori not known, is a free parameter of the model. We refer to this scaling factor as the effective active fraction (EAF). This factor only has a physical interpretation for a pure ice surface where it would represents the fraction of pure ice of an areally mixed surface needed for a specific sublimation flux. In general though it is not a physical parameter and should not be interpreted as such.

In the next steps, we calculate the gas and dust flow fields in three dimensions. We then perform a column integration along the line-of-sight through the dust coma for a specific viewing geometry of the OSIRIS NAC (narrow-angle) and WAC (wide-angle) cameras [14] and convolve the dust column densities with the optical properties of the dust to arrive at absolute radiance values that can be compared with the OSIRIS images. One major assumption that goes into this approach is that there is no significant back-coupling from the dust to the gas allowing a sequential treatment of the two flows. For low dust-to-gas mass ratios, this is certainly justified [16] but will break down when a lot of dust is released. We will further discuss this limitation later on.

2.2. Gas Kinetic Simulations

The gas flow-field is calculated using the DSMC technique. The code used is called UltraSPARTS¹ and is a commercialized derivative of the PDSC++ code [20] used in previous papers (e.g., [16, 17]). PDSC++ is a C++ based, parallel DSMC code which is capable of simulating 2D, 2D-axisymmetric, and 3D flow fields. The code has been developed over the past 15 years [21–23] and contains several important features including the implementation of 2D and 3D hybrid unstructured grids, a transient adaptive sub-cell method (TAS) for denser flows, and a variable time-step scheme (VTS). In the parallel version, computational tasks are distributed using the Message Passing Interface (MPI) protocol. The improved UltraSPARTS (Ultra-fast Statistical PARTicle Simulation Package) has been applied to 67P [18, 19]. Here we simulate the full 3D gas flow up to a distance of 10 km from the nucleus center.

The sublimation temperature and flux—calculated as described above—for each facet are set as initial conditions of the simulation. This includes implicitly the assumption of the appropriate EAF. We assume here that the EAF of all facets are the same (i.e., homogeneous surface properties) but can change from epoch to epoch. This results in one value for the global EAF per epoch. Though we know from previous works (e.g., [5, 16, 24, 25]) that there are regional inhomogeneities that can be encoded in EAF it is not the focus of this work to constrain these inhomogeneities. Rather we seek a global estimate of the fluxes and dynamical behaviors. Because the EAF is a free parameter it needs to be constrained by data. In our case, we determine the EAF by comparing modeled densities extrapolated to the Rosetta position and actual COmet Pressure Sensor (COPS; [15]) measurements during each epoch. Within each epoch where we match the sub-solar longitude of a measurement, we extract from the respective simulation the gas number density at the position of the spacecraft. If the spacecraft distance is larger than 10 km we extrapolate the value from the 10 km surface to the spacecraft distance assuming free radial outflow. This assumption is well-justified as shown in Marschall et al. [16]. Though this does not capture the detailed structure of the ROSINA/COPS data it does account accurately for the average activity level at each epoch. Table 2 shows the EAF used and the resulting average global H₂O production rate for one comet day.

TABLE 2

Table 2. The effective active fraction (EAF) assumed and the resulting average H₂O production rate for each epoch.

The global gas production rate as a function of time is shown in Figure 11 (purple band in top panel). Because we have used the total ROSINA/COPS data—which also contains the other gas species other than water—to constrain our emission, these values should be understood as a proxy for the entire emission. We have interpolated between the epochs using a local second-order polynomial. The fit for each epoch, i, includes three epochs i−1, i, i+1. The fitting parameters to calculate the mean gas production rate as a function of ephemeris time (ET) is shown in Table 3. The resulting integrated mass loss over the shown period adds up to (6.1 ± 1.5) · 10⁹ kg. This is well in line with values published in other works as e.g., (6.3 ± 2.0) · 10⁹ kg [7] or (5.8 ± 1.8) · 10⁹ kg [6]. The error arises from the uncertainty of the data (up to 15%; [26]) although the relative errors are probably smaller (M. Rubin, pers. comm.) and our model (5–10%; [27]) as well as from the scatter from the comparison of the data and our model. As it was not the goal of this work to constrain as precisely as possible the surface-emission distribution it is nevertheless noteworthy that our estimates come so close to the other published values. This illustrates that it is not necessary to know the surface-emission distribution well to estimate the total global volatile loss. Rather simple assumptions of the surface response is sufficient for such an estimate. Though it is not that surprising, because as pointed out in Marschall et al. [28] the global gas production rate can be fairly well-estimated by even simplified models. Our peak production rate is reached at the summer solstice (epoch 12) and not perihelion (epoch 11) and therefore roughly 22 days post-perihelion. This is in line with dust coma measurements by OSIRIS. Gerig et al. [18] reported peak dust coma brightness 20 days post-perihelion. This also hints at the fact that the obliquity plays an important role in the activity of comets. Though the heliocentric distance still is the main driver of the gas and dust activity (O(1)) it is the obliquity/season that controls the second order. The coincidence of the peak gas activity with the peak dust activity also indicates that the dust activity is mainly driven by H₂O or at least near-surface volatiles without a significant thermal lag.

TABLE 3

Table 3. The mean gas production rate, ${\stackrel{-}{q}}_g$ [kg s⁻¹] as a function of ephemeris time (ET): ${\stackrel{-}{q}}_g(ET)=a·{ET}^2+b·ET+c$.

2.3. Dust Dynamic Simulations

After the gas flow field has been evaluated, we calculate the dust flow field by injecting dust test particles into the flow. We use a typical approach for computing the dust motion in a gas flow-field taking into account gas drag and the comet gravity using our DRAG3D dust coma model detailed in Marschall et al. [16] and the references therein. We assume that the dust mass production rate is proportional to the gas mass production rate and that the dust size distribution does not vary across the surface except in cases where certain dust sizes are no longer lifted because the gas pressure is too low to surpass the local gravity. The dust size distribution is thus only naturally modified by the dynamics and lifting process. It is assumed that the dust particles are at rest on the surface (i.e., the ejection velocity is 0). The dust-to-gas mass ratio as well as the dust size distribution at the surface are free parameters of the model and will be constrained by the data as described below. Due to the presence of gravity, large dust particles may not reach escape speed and eventually return to the surface. The flux of back-fall particles is thus a further output of the model. It is important to note that we assume that the dust particles are desiccated, i.e., contain no significant amounts of volatiles that evaporate while airborne. They may still be wet but do not outgas significantly. This is a consequence of our assumption that there is no significant back-coupling of the dust flow onto the gas flow.

For each epoch (except for two) we have selected one OSIRIS image where the illumination conditions of the image match one of the gas simulations (see Table 1). The images used in this work, as well as some camera and geometric properties, are shown in Table 4. Two main criteria were used to select these images. First, the images needed a large enough field of view such that projected distance (impact parameter, b) in the image plane from the center of the comet to the edge of the image at each side was at least 9 km. Why this is an important constraint will be described in the next paragraph. Second, images need a sufficient signal to noise such that the dust coma brightness could be measured well. These two constraints unfortunately, eliminated all images for epoch 2 and 20. Epoch 2 included the 10 km orbit phase and thus did not provide large enough fields-of-view while the signal-to-noise was bad in epoch 20 due to the very low activity of the comet. Most images we have used were taken by the wide-angle camera (WAC) and filter 18 (central wavelength, λ_c = 612 nm) and at cometocentric distances between 87 and 635 km and phase angles between 37° and 108°.

TABLE 4

Table 4. List of OSIRIS images used in each of the epochs as well as their filter, central wavelength (λ_c), phase angle (α), and cometo-centric distance (D_cc).

The dust field is calculated for each image using 41 different dust sizes from 10 nm to 1 m. The dust sizes are logarithmically spaced with five dust sizes per decade. The particles are assumed to be spherical and have a density of 533 kg m⁻³ matching roughly the bulk density of the nucleus [2]. Even though all dust sizes are simulated, not all of them contribute to the dust brightness in the coma. This is because the particles larger than a certain size might not all be lifted because the gas pressure cannot overcome gravity. Thus the number of dust sizes present in the coma depends on the heliocentric distance (epoch). The upper size limit (largest liftable size) is thus naturally determined and thus an outcome of the simulation. What the smallest dust size should be is unless. The smallest diameter of particle sub-units measured by MIDAS [29] is 100 nm. Whether these could also be the smallest dust particles in the coma or if these measurements have an in-situ collection bias at the spacecraft is not clear. One could imagine that very small particles might not have been collected because of spacecraft and dust charging. Della Corte et al. [30] showed that particles, for which the ratio of the particle charge to its kinetic energy entering the electrostatic field of the space craft $q/E_k>0.24C{J}^{-1}$, will not reach the spacecraft. We will therefore leave this issue open for the moment and examine the impact of the smallest size on the results in section 4.

Once the 3D dust field is simulated we calculate the dust column densities of each size for the specific viewing geometry of the respective OSIRIS image. The final image is composed by weighting the different dust sizes according to a specific dust size distribution and convolving the column densities with the scattering properties described in section 2.4.

For each of the images we compare the integrated radiance of the dust coma along an aperture with impact parameter b = 11 km and compare it to that of the synthetic images. Again it is not our goal in this work to match the structure of the inner dust coma but rather the overall global behavior. Gerig et al. [18] showed from OSIRIS data that the dust flow goes over to free radial outflow at an average impact parameter of b ~ 11 km. This is in line with theoretical considerations of dusty flows [31]. Beyond that point, the dust brightness falls off with 1/b as expected for a freely expanding radial flow. For that reason, we have chosen b = 11 km to be within the free-flow regime. If the field-of-view was not large enough we used the maximum available impact parameter.

2.4. Scattering Model

Previously, we have used a spherical particle model and a Mie scattering code in our modeling pipeline. Here we use a much more sophisticated approach based on the recently introduced radiative transfer with reciprocal transactions framework [32, 33]. The approach allows for scattering analysis of large irregularly shaped particles with wavelength-sized details. Here, the dust particles are considered to be irregular aggregates composed of sub-micrometer-sized organic grains and micrometer-sized silicate grains. Such a particle model has been found to be in good agreement with OSIRIS [34] and VIRTIS [35] phase function measurements. The refractive index for silicate grains is assumed to be m = 1.6048692+i0.0015341 corresponding to magnesium iron pyroxine [36] and for organic grains m = 1.55950+i0.42964 corresponding to amorphous carbon [37]. At 612 nm (WAC filter 18) the resulting scattering phase functions for different particle sizes normalized to the geometric albedo are shown in Figure 2. The figure shows good agreement of the phase function of large particles (>1 cm) with the nucleus phase function as measured by OSIRIS [39–41]. This should indeed be the case because larger particles should behave more and more like “small comets” themselves and thus be representative of the nucleus scattering properties. For small particles, the best agreement of a single dust size with the coma phases function [38] is between 10 and 100 μm. The numerical method of Markkanen et al. [34] is not applicable to particles smaller than 1μm. Thus for the particles smaller than 1μm we use a Mie scattering code to determine the scattering properties [see 16] matching the single scattering albedo of the Mie result with the approach of Markkanen et al. [34] for 1μm particles. This gives us a smooth transition from the large particle region to the Rayleigh scattering region where particle's shape has a negligible effect on its scattering properties. This is a state of the art model and we have thus used its results throughout this work. But because the scattering model does have an effect on the results a re-evaluation of the results can be done if and when a better model arises.

FIGURE 2

Figure 2. Shows the phase functions, S₁₁, as a function of phase angle for different particle sizes (solid lines). The green stars represent the measurements of the coma phase function by Bertini et al. [38] while the blue stars show the nucleus phase function as measured by Fornasier et al. [39], Feller et al. [40], and Masoumzadeh et al. [41].

3. Theoretical Consideration

To put some of our results in the next section (section 4) into context, we present first some general theoretical considerations of what we can expect, in particular with regards to the relationship between the dust-to-gas ratio and the dust size distribution. We thus consider first a simple model where the comet is represented as a sonic (i.e., the gas velocity near the surface is equal to the local sound velocity and defined by the thermodynamic properties of the gas and the surface temperature i.e., R_h) spherical source of ideal perfect gas (with specific heat ratio γ=1.33) accelerating spherical solid grains. The source shall have radius R_N, nucleus mass M_N, total gas production rate Q_g (kg/s). The motion of a spherical grain in a flow from such a source was studied in Zakharov et al. [31] for a wide range of conditions. They defined

$\begin{array}{ll}\text{Iv}=\frac{3Q_g}{32R_{N}r{\rho }_d\pi v_{g0}},& (1)\end{array}$

and

$\begin{array}{ll}\text{Fu}=\frac{G{M}_N}{R_N}\frac{1}{{(v_g^{max})}^2}& (2)\end{array}$

which are dimensionless parameters, where ρ_d is the specific density of the dust particles, v_g0 the gas velocity near the surface, $v_g^{max}$ the theoretical maximal velocity of gas expansion (defined by the thermodynamic properties of the gas and the surface temperature i.e., R_h), and G the gravitational constant. Iv characterizes the ability of a dust particle to adjust to the gas velocity while Fu quantifies the importance of gravity. Zakharov et al. [31] found that for Iv < 0.1 (which is the case of 67P, and dust sizes >1 nm) the dust particles reach 90% of their terminal velocity at about 6 · R_N. The terminal velocity of particles with radius, r, varies as $v_d(r)\propto r^{-0.5}$ for small Fu (i.e., if gravity plays a minor role). The asymptotic dust velocities are given by:

$\begin{array}{ll}{v}_d(r)={\left(\frac{\text{Iv}(r)}{\text{Iv}(r_{*})}\right)}^{1/2},v_d(r_{*})={\left(\frac{r}{r_{*}}\right)}^{-1/2},v_d(r_{*})=r^{-1/2}{C}_{Iv}& (3)\end{array}$

where r_* and v_d(r_*) are some referential size and corresponding terminal velocity, and C_Iv is a constant.

For a dust size distribution given by a power-law, r^−q, the normalized mass distribution, f_md, of particles ejected from the surface is

$\begin{array}{ll}{f}_{md}(r)=\{\begin{array}{cc}\frac{4-q}{r_{max}^{4-q}-r_{min}^{4-q}}{r}^{3-q},& q\ne 4\\ ln\left(-\frac{r_{max}}{r_{min}}\right)r^{3-q},& q=4\end{array}& (4)\end{array}$

where r_min and r_max are the smallest and largest dust sizes ejected from the surface. In the following, we will not considered specially the case of q = 4. The dust production rate, Q_d, of each dust size is

$\begin{array}{ll}{Q}_d(r)=\chi f_{md}(r)Q_{g}dr,& (5)\end{array}$

where χ = Q_d/Q_g is the total dust-to-gas mass loss rate. Therefore, the number density of dust particles with radius, r, at the radial distance, R, from the center of the nucleus is:

$\begin{array}{ll}n(r,R)=\frac{\chi f_{md}(r)Q_{g}dr}{v_d(r)m_d(r)4\pi R^2}.& (6)\end{array}$

The column density at the distance ϱ from the center of the nucleus in the image plane is:

$\begin{array}{ll}{n}_{col}(r,\varrho )={\displaystyle {\int }_{-\infty }^{+\infty }}n(r,R)dz=\frac{\chi f_{md}(r)Q_{g}dr}{v_d(r)m_d(r)4}\frac{1}{\varrho }.& (7)\end{array}$

The total number of dust particles in a column within a circular observing aperture of radius ℜ is:

$\begin{array}{ll}N(r,\Re )={\displaystyle {\int }_0^{\Re }}{n}_{col}(r,\varrho )2\pi \varrho d\varrho =\frac{\chi f_{md}(r)Q_{g}dr}{v_d(r)m_d(r)2}\pi \Re .& (8)\end{array}$

The brightness is proportional to the flux F (W/m²) gathered by an instrument which for an optically thin coma is:

$\begin{array}{ll}F(r,\Re )={F}N(r,\Re )\pi r^2{q}_{sca}(r)\frac{{\phi }_{av}(r)}{4\pi }\frac{1}{{\Delta }^2}& (9)\end{array}$

where ${F}$ is the incident flux, Δ is observational distance, q_sca is scattering efficiency and φ_av is the phase function averaged over phase angle. Substituting Equations (3), (4) and (8) in (9) we get:

$\begin{array}{ll}F(r,\Re )=\frac{3}{32}\frac{{F}{Q}_g\Re }{C_{Iv}{\rho }_d{\Delta }^2}\frac{4-q}{r_{max}^{4-q}-r_{min}^{4-q}}{r}^{\frac{5}{2}-q}\chi q_{sca}(r){\phi }_{av}(r)dr.& (10)\end{array}$

For fixed ${F}$, Q_g, R_N, ρ_d, v_g0, r_*, v_d(r_*), Δ, and ℜ

$\begin{array}{ll}F(r,\Re )=C\frac{4-q}{r_{max}^{4-q}-r_{min}^{4-q}}\chi r^{\frac{5}{2}-q}{q}_{sca}(r){\phi }_{av}(r)dr& (11)\end{array}$

where $C=\frac{3}{32}\frac{{F}{Q}_g\Re }{C_{Iv}{\rho }_d{\Delta }^2}$ is a constant (for the given observational conditions).

For the optical properties, we make some simplifying assumption. First, we approximate the scattering efficiency q_sca(r) to be:

$\begin{array}{ll}{q}_{sca}(r)=\{\begin{array}{cc}0.233·10^{24}·r^{7/2},& 10^{-8}\le r<2·10^{-7}\\ 4.993·10^{-5}·r^{-2/3},& 2·10^{-7}\le r<2·10^{-4}\\ 0.02,& 2·10^{-4}\le r\le 1.0\end{array}& (12)\end{array}$

Figure 3 shows the computed q_sca from section 2.4, the fitted q_fit scattering efficiency and the relative difference. In this fit, we used “round numbers” (i.e., this is a very rough fit). This simplification results in differences of <50%. For the phase function φ_av we assume it to be constant for all sizes. We estimate an error of the order of a factor of 2 from this simplification.

FIGURE 3

Figure 3. Scattering efficiency as a function of dust radius: computed (comp, red) according to section 2.4, and fitted (fit, green). The relative difference between the computed and fitted is shown in blue.

Under the assumption we made the integration of Equation (11) becomes trivial.

For a given gas production rate Q_g the maximum dust size a_max is also constant. Figure 4 shows how the dust-to-gas mass loss rate Q_d/Q_g varies as a function or the power-law exponent of the dust size distribution for the same brightness. With increasing power-law exponent q from minimal value to ≈3 the Q_d/Q_g is slowly decreasing (since in this case practically all dust mass is concentrating in a narrow range of largest sizes), but with increasing q from 3 to 4 the ratio Q_d/Q_g is strongly decreasing. For 4.5 < q < 6 Q_d/Q_g is increasing (within one order of magnitude). This inflection point of Q_d/Q_g occurs at the transition from dust grains distribution with most mass being in the large dust sizes to where most mass is in the small dust sizes.

FIGURE 4

Figure 4. Dust-to-gas mass production rate ratio vs. power-law exponent for constant dust brightness.

The growth of Q_d/Q_g for q > 4.5 is a consequence of the strong decrease of q_sca for small sizes, therefore, in order to maintain the same brightness, it is necessary to eject more dust.

We should remember that this analytical result (Figure 4) assumed several important simplifications:

1. We assumed that the dust expansion is strictly radial;

2. For evaluation of the dust brightness we used simplified optical properties (e.g., isotropic phase function);

3. We assumed that gravity plays only a minor role;

4. We assumed that the dust does not affect the gas flow.

We will discuss in the next section how these assumptions [in particular (1) and (3)] change the result.

4. Results and Discussion

To convolve the results of the dust dynamics model with the scattering properties to arrive at synthetic OSIRIS images we need to assume a dust size distribution. As is commonly done we presume that the number of particles, n, of a certain radius, r, follows a power-law:

$\begin{array}{ll}n(r)\propto r^{-q},& (13)\end{array}$

where q is the differential power-law exponent. Figure 5 illustrates an example of an OSIRIS and synthetic image for epoch 12 (solstice). As described in section 2.3, we extract the integrated brightness along a circle with a constant impact parameter of b = 11 km where possible (illustrated in the figure with the red circles). We should stress here again that it was not the aim of this work to match the emission distribution on the surface and thus all the structures in the coma.

FIGURE 5

Figure 5. Panel (A) shows the OSIRIS image WAC_2015-08-09T09.13.16.574Z of epoch 12 (solstice) with an enhanced contrast to show the dust coma. Panel (B) shows the synthetic image with power-law exponent of q = 4 of our dust model. The crosses in each panel indicate the center of the nucleus and the red circles indicate an impact parameter distance of b = 11 km.

For a given gas production rate, the three major factors controlling the brightness (see Equation (11) for more detail) of the dust coma are:

1. The dust-to-gas mass production rate ratio, Q_d/Q_g, at the surface;

2. The dust size distribution (i.e., the power-law exponent, q) at the surface;

3. The scattering properties of the dust particles.

We should note that although we assume a uniform surface (i.e., globally constant Q_d/Q_g and q) the actual values at each facet vary depending on the local gas flux. If a particular facets' local flux is too low to lift a certain particle size the resulting dust flux and size distribution from that surface facet will differ locally from the nominal values. The three input parameters above are not known a priori and are thus initially free parameters and in need of constraints. We have fixed the scattering properties by using published results that fit OSIRIS data [see sections 2.4 and 34]. This reduces the above parameters from three to two.

There are three further quantities that influence the coma brightness, but as we will show below their influence is small compared to the ones mentioned above. These are:

1. The smallest dust size, r_min;

2. The largest dust size, r_max;

3. The bulk density of the dust, ρ_d.

Of these three we will explore the influence of r_min and ρ_d on our results. We will not be artificially truncate the size range at large sizes by varying r_max. On the contrary, r_max will be naturally regulated due to the balance of forces at the surface. If a given local gas flux is not sufficient for lifting a certain size, that will determine the largest dust size from that surface element.

Apart from the parameters that directly influence the brightness, several indirect factors further constrain the curves Q_d/Q_g and q. We will discuss these constraints at the end of this section.

As has already been shown in Figure 12 of Marschall et al. [16] Q_d/Q_g and q are not independent. Knowing the brightness of the coma these two parameters constrain each other to a limiting set of parameter pairs. Figure 6 shows Q_d/Q_g as a function of q for each of the 18 OSIRIS images of this study. Each line represents an equal brightness curve where the respective q and Q_d/Q_g result in a synthetic image that matches the OSIRIS brightness. Several things are noteworthy. First, all curves show minima between q = 4 and q = 4.5 and thus illustrate the inherent degeneracy between Q_d/Q_g and q. Second, all cases with shallow power-laws (q < 3) require very large Q_d/Q_g of at least 10 but in most cases around 100. Third, steep power-laws (4 < q < 6) in all but one case require much less dust mass, i.e., Q_d/Q_g ≤ 1 to match the brightness of OSIRIS. Fourth, there is a clear trend in the gas production rate. As the gas production rate increases the slope in Q_d/Q_g for shallow power-laws (q < 3) becomes shallow, too. Or conversely for low gas production rates very high Q_d/Q_g are needed to match the OSIRIS brightness when the power-law is shallow. This has to do with the amount of dust that can be lifted and escape the nucleus' gravity.

FIGURE 6

Figure 6. The dust-to-gas mass production rate ratio, Q_d/Q_g as a function of the power-law exponent, q, is shown for each of the 18 OSIRIS images of this study. Each line represents an equal brightness curve where the respective q and Q_d/Q_g result in a synthetic image that matches the OSIRIS brightness. The colors of the curves indicate the global gas production rate of the respective epoch.

Comparing Figure 6 to the analytical solution presented in Figure 4 of section 3 we see that for high gas production rates the model follows the analytical solution rather well. The places where we deviate from the analytical solution illustrate the effect of different physical processes. For the analytical solution we have assumed a minor (but not negligible) role of gravity. The effect of gravity can be seen in the low gas production rate cases with shallow power-laws. There, in contrast to the analytical model which levels off at smaller power-law exponents, the dust coma model results in ever higher Q_d/Q_g. This is caused by the inability to lift large particles from the entire surface and therefore a higher Q_d/Q_g is required to maintain the brightness. Thus, the deviations at low gas production rates and shallow power-laws exhibit the non-minor role of the nucleus' gravity. As in the analytical model for steeper size distributions most mass is in the smallest particles, which are weakly scattering and thus hardly contribute to the brightness. This is compensated by an increase of Q_d/Q_g at these steep power-laws.

Compared to previous work presented in Figure 12 of Marschall et al. [16] the Q_d/Q_g values we find here (in particular for q < 3.5) are much higher while the behavior of the curves for steeper power-laws is within the expected range. The two main reasons we find larger values at shallow slopes are: (1) Marschall et al. [16] assumed the scattering properties of astronomical silicate [42] which is much brighter than we now know; (2) we consider here considerably larger dust sizes as our upper limit. This extension of the size domain increases the Q_d/Q_g by orders of magnitude because of high fall back fractions of dust that is gravitationally bound and weakly scattering.

Two assumptions going into Figure 6 are worth discussing. First, we have assumed that all dust particles have a bulk density equal to the nucleus density (533 kg m⁻³). It is likely that the density of small particles is significantly larger than that and that the density then decreases with size. Because the exact relationship of the density as a function of dust size is currently unknown we have not tested a varying density as a function of size. But, we have varied the bulk density for the entire range of dust sizes between 250 kg m⁻³ and 1, 000 kg m⁻³. The two left panels of Figure 7 show the results for a moderate activity environment (epoch 6—inbound equinox, Q_g = 35 kg s⁻¹) and a high activity environment (epoch 12—solstice, Q_g = 800 kg s⁻¹). For 3 < q < 3.75 the differences between the different dust densities is minimal. For q > 3.75 the differences are larger, in particular in the high activity case. How the bulk dust density impacts the total dust mass loss will be explored later in this section.

FIGURE 7

Figure 7. The four panels show the dust-to-gas ration as a function of power-law index. The left two panels show the results for different bulk dust densities (250 kg m⁻³ [purple], 533 kg m⁻³ [green], 800 kg m⁻³ [blue], 1, 000 kg m⁻³ [orange]) assuming a minimum dust size of 0.01μm and a gas production rate of 800 kg s⁻¹ (top panel, epoch 12) and 35 kg s⁻¹ (bottom panel, epoch 6). The two right panels show the results for varying minimum dust size of the power-law from 0.01μm (red) to 1μm (blue) assuming a bulk dust densities of 533 kg m⁻³ and a gas production rate of 800 kg s⁻¹ (top panel, epoch 12) and 35 kg s⁻¹ (bottom panel, epoch 6).

Second, we have currently assumed that the smallest dust size is 0.01μm. This might not be the preferred choice and a much larger smallest size should be considered. The MIDAS instrument detected 1μm particles (e.g., [29]) and there is indirect evidence of sub-micron particles observed by VIRTIS during outbursts [43]. We have thus explored the range of the smallest sizes between 0.01 and 1μm. The two right panels of Figure 7 explore the effect of the smallest size on the dust-to-gas ratio by varying the smallest size. Compared to the differences seen for different bulk dust density the effect of the smallest size is quite substantial. As we would expect the smallest size does not affect the result for q < 3 as in these cases most of the mass is in the large particles. As q increases from 3 the curves for different smallest sizes start diverging. Two trends can be observed. As the smallest size increases from 0.01μm to ~0.1μm the dust-to-gas ratio starts to flatten out beyond q = 4. This is caused by the fact that the size distribution is no longer dominated by very inefficiently scattering particles. As the smallest size continues to increase to 1μm the overall dust-to-gas ratio increases. This is because the most efficient scatterers (see Figure 3) are being removed from the size distribution and must be compensated by more mass of all other sizes. For very steep power-laws the difference in the dust-to-gas ratio can be up to 1.5 orders of magnitudes. How the smallest size impacts the total dust mass loss will be explored later in this section.

How the fraction of gravitationally bound dust mass, which falls back to the nucleus, varies as a function of global gas production rate is shown in Figure 8 for different power-law exponents. This illustrates that for very low gas production rates and very shallow power-laws (q < 2.5) almost the entire dust mass emitted from the surface will be redeposited. This explains the large increase seen in Q_d/Q_g of Figure 6. Conversely, in the case of steep power-laws (q ≥ 4.5) almost all of the dust escapes the nucleus' gravity field irrespective of the gas production rate. In all cases, the fraction of fall back decreases as the gas production rate increases. Therefore, the fraction of fall back material decreases as the comet approaches the Sun. At large heliocentric distances, large fractions of dust emitted will return to the surface (i.e., >50% for q < 3.5). But the gas and dust production rates are highest at perihelion/solstice, thus the total amount of fall back during one apparition is dominated by the fraction of fall back during that period.

FIGURE 8

Figure 8. The ratio of dust mass falling back, Q_Dfb to total dust mass, Q_D is shown as a function of global gas production rate for different power-law exponents.

The fraction of fall back is also tightly bound to the maximum liftable dust size. Figure 9 shows as a function of global gas production rate the largest dust size that can still be lifted from the surface of the nucleus. The figure also shows the largest dust size that can escape the gravity field of the comet. As the gas production rate increases so do the largest liftable and escaping dust sizes. For Q_g > 300 kg s⁻¹ the largest liftable dust size is larger or equal to 1 m, which is the largest size in our model. Though these sizes, or larger, can be lifted they will not be able to escape the gravity field of the comet and will be redeposited on the surface. The largest size that can escape the comet at peak gas production is roughly 0.6 m. We should note though that this calculation neglects surface cohesion, solar radiation pressure, and heat transport to the subsurface that is needed to eject such large particles. Here we only consider the balance of gas drag and gravity to determine these largest sizes.

FIGURE 9

Figure 9. The largest liftable dust size (red curve) and the largest escaping dust size (blue curve) are shown as a function of the gas production rate.

The discussion of the previous paragraphs illustrates that multiple properties of the dynamical simulation of the dust coma (size distribution, dust-to-gas ratio, and the fraction of fallback) as well as the optical properties of the dust are not independent but mutually constraining. Example, a given fraction of fallback implies a certain size distribution which in turn constrains the possible dust-to-gas ratios for a particular set of scattering properties. Though for a particular single OSIRIS image these parameters can be constrained there still remains a rather large set of parameters that are consistent with the data as presented to this point (including dust coma brightness in the OSIRIS images and local gas densities of ROSINA/COPS).

While we have only considered constraints within each epoch there is one strong constraint covering the entire mission. That is the measurement of the total mass loss during the Rosetta apparition. During the 2 year mission comet 67P had lost (10.5 ± 3.4) · 10⁹ kg [4]. The total mass loss, $M_{tot}=(10.5\pm 3.4)·10^9$kg, is:

$\begin{array}{ll}{M}_{tot}=M_g+M_d^{esc},& (14)\end{array}$

where M_g is the total volatile mass loss, and $M_d^{esc}$ is the total escaping dust mass. For the dust masses, we can further specify that

$\begin{array}{ll}{M}_d=M_d^{esc}+M_d^{fb},& (15)\end{array}$

where M_d is the total dust mass ejected from the nucleus, and $M_d^{fb}$ is the dust mass that falls back to the surface. We have determined the total volatile mass loss to be (6.1 ± 1.5) · 10⁹ kg. Combined with the total mass loss of the nucleus it follows that $M_d^{esc}=(4.4_{-4.4}^{+4.9})·10^9$ kg. Note that within this interval exists the possibility that $M_d^{esc}=0$ kg. Though we know that dust escaped from the nucleus the simple mass balance would not exclude this possibility. We can now integrate the total dust mass loss over the orbit of the comet for different power-law exponents. For the integration, we assume a linear interpolation of the results between epochs. Figure 10 shows the $M_d^{esc}$ as a function of the power-law exponent, q. Cases, where $M_d^{esc}$ exceeds the nominal dust mass loss of 4.4 · 10⁹ kg (horizontal dashed green line) or the maximum dust mass loss of 9.3 · 10⁹ kg (horizontal dashed red line), can be discarded. Furthermore, where the mass loss curve intersects the mass loss indicates the corresponding power-law that fits the data. Figure 10 also illustrates the effect of the smallest size—discussed earlier in more detail for an individual OSIRIS image—on the total mass loss. The effect of the smallest size is rather limited for 0.01μm < r_min < 1μm because the dust mass curves cross between 3.5 < q < 4. This also implies that the effect of the bulk dust density is even smaller than the effect from the choice of the smallest size (see discussion about Figure 7 above). We can also see that for r_min > ~12μm there will no longer be a nominal solution to the constraints. Further, for r_min > ~30μm there is no solution at all because the curve will stay above the maximum escaping mass for all power-law exponents studied here. This means that the minimum dust size must be strictly small than ~30μm and nominally even smaller than ~12μm. Figure 10 illustrates how we can determine the power-law exponent for the nominal and maximum dust mass loss, which in turn determines the dust-to-gas ratio, dust production rates, fraction of dust fallback. As the minimum size grows larger than 1μm the required power-law exponent increases and becomes rather large. There remains the issue of the minimum escaping mass. As discussed above the lower limit according to the total and volatile mass loss is zero. But for our models, the minimum escaping dust mass is never zero. We have thus chosen the smallest possible mass loss of each model as the minimum mass loss. The resulting power-law exponents, dust mass losses, dust-to-gas ratios, and fall back fraction are summarized in Table 5. We have also determined the deposition height, H, that results if the fallback material is spread equally on the smooth deposits (9.43 km²) identified by Thomas et al. [12].

FIGURE 10

Figure 10. The total escaping dust mass, $M_d^{esc}$, is shown as a function of the power-law index for five different minimum dust radii, r_min. The horizontal dashed lines show the nominal ejected dust mass (green) and the maximum ejected dust mass (red).

TABLE 5

Table 5. Power-law exponents, dust mass loss, dust-to-gas ratio, and fall back as a function of the smallest dust size (see Figure 10).

The results in Table 5 show that the integrated quantities are rather insensitive to the choice of the smallest dust size if r_min ≤ 1μm. For minimum sizes larger than 1μm the power-law becomes steeper and thus the amount of dust fall-back goes down. The dust-to-gas ratio is rather stable for all cases and is of the order of 0.8 with an error of the order of 100%. This means that while the nominal case reflects a comet that contains more volatiles than dust the case of a dusty comet lies within the error.

The fallback in all cases is of the order of 10% and results in a deposition height of the order of 10 cm. Because the deposition is likely non-uniform it is therefore easily thinkable that in certain areas dust of the order of meters is deposited while in others only a few centimeters.

This analysis assumes that the dust size distribution does not change along the orbit. There is an indication (e.g., [44]) that this is not the case and that the slope is varying with heliocentric distance. Our model cannot resolve/constrain this. All the quantities here are heavily dominated by the period around perihelion and summer solstice when the emission was the highest. Therefore, the power indexes found here reflect the values for this period.

The power-law we find to be compatible with the data is an independent result based only on the brightness of the dust coma and the total mass loss balance. Because most mass is ejected around perihelion, this power-law mainly reflects this period and deviations of it at larger heliocentric distances [45] would not influence the result. The value we find is in line with other measurements around perihelion e.g., the in-situ measurement of q = 3.7 by GIADA [45], q = 3.8 by COSIMA [44], as well as ground-based estimates for the dust tail of 3.6 < q < 4.3 for sizes smaller than 1 mm and q = 3.6 for sizes larger than 1 mm [46].

A check of our dust dynamics model is the comparison of our model dust speeds with the ones measured by GIADA. For the period between 2.2 AU inbound to 2.0 AU outbound [47] have reported 141 dust particle detection for which a dust speed and mass could be determined. Of these, the smallest particle had a mass of 2.8 · 10⁻⁹ kg, which corresponds to a radius of 108 μm assuming a spherical particle with our nominal dust density. The measured dust speeds varied between 0.3 and 34.8 m s⁻¹ [47]. A further constraint is the fact that Rosetta spent ~65% of the time at phase angles of ~90° and an additional ~25% of the time at phase angles of ~60° which implies that the particles were mainly collected in those locations (see also Figure 5 in [47]). At a phase angle of 90° our model dust particles with radius 100 μm have speeds of (3.5 ± 1.0) m s⁻¹ at the inbound equinox (epoch 6) and (18 ± 5.6) m s⁻¹ at the summer solstice (epoch 12). For phase angles of 60° the model dust particles have speeds of (7.0 ± 1.2) m s⁻¹ at the inbound equinox and (32 ± 5.2) m s⁻¹ at the summer solstice. Our dust speeds are thus well in line with the measured speeds by GIADA given that larger particles will have lower speeds than the ones listed above.

We should highlight that our peak dust production rate (~530 kg s⁻¹) is roughly an order of magnitude lower than those reported by e.g., [46] (~3, 000 kg s⁻¹) or [48] (~8, 300 kg s⁻¹). Furthermore, [46] report a total dust mass loss of 1.4 · 0¹⁰ kg. As neither [46] nor [48] report any error bars on their results, we cannot asses if they are plausible. If taken at face value both results are inconsistent with the RSI measurement of the total mass loss of the comet [4] taking into account the estimates of the volatile mass loss in this work and others [6, 7].

Finally, the determination of the power-law exponent allows us to determine the dust production rate (Figure 11, top panel, green band), dust-to-gas ratio (Figure 11, bottom panel, green line), and fraction of fallback (Figure 11, bottom panel, purple line) as a function of time. The dust production rates are linearly interpolated between the epochs. Unfortunately, our model is rather noisy but the overall trends are robust enough that we feel comfortable making further conclusions. The fraction of dust fallback is highest at large heliocentric distances and then decreases toward perihelion and reaches its minimum at summer solstice where the activity peaked. Though the faction of fallback is smallest at the peak of the activity (solstice), most mass that is falling back will still be from the period of summer solstice because of the high activity. The behavior of the fraction of dust fallback is symmetric for the inbound and outbound part of the comets' orbit. Contrary to that the dust-to-gas ratio is highest (~1.5) at large heliocentric distances inbound and keeps decreasing along the entire orbit and does not significantly increase on the outbound leg but rather flattens out at ~0.1. This might be indicative of the comet shedding its dust mantle, in particular in the northern hemisphere. This trend of decreasing dust-to-gas ratio along the orbit manifests itself also in the asymmetry of the global dust production. To first order, the dust production rate follows the gas production rate during the inbound leg but the dust production rate drops faster than the gas production rate post solstice. This is also observed in ground-based measurements [46]. There is also an intriguing spike in the dust-to-gas ratio after the inbound equinox coinciding with an increase in the total dust production rate. Future in-depth work will be needed to confirm the nature of this feature which does not seem to be present in the observations of the outer coma from ground-based measurements. But if it is truly there it can be understood as the comet shedding its southern dust mantle because the feature coincides with the period when the southern hemisphere receives increasing insolation.

FIGURE 11

Figure 11. The top panel shows the global production rate for the gas (purple) and dust (green) as a function of days to perihelion. The bands indicate the range due to the diurnal variation. The gas production rates have been constrained by ROSINA/COPS measurements while the dust production rates are from combined constraints of OSIRIS and gas fluxes. For the dust a minimum size of r_min = 0.1 μm and power-law exponent of q = 3.75 are assumed. The bottom panel shows the fraction of dust fall-back (purple) and dust-to-gas ratio (green).

5. Summary and Conclusions

In this work, we have simulated the inner gas and dust coma of comet 67P covering the entire Rosetta mission by splitting it into 20 epochs. The gas production rates of each epoch were constrained by in-situ measurements of the gas density by ROSINA/COPS. From that, the total gas mass loss is estimated at (6.1 ± 1.5) · 10⁹ kg. This is in line with values published in other works as e.g., (6.3 ± 2.0) · 10⁹ kg [7] or (5.8 ± 1.8) · 10⁹ kg [6]. It also illustrates that it is not necessary to know the surface-emission distribution well to estimate the total global volatile loss.

By simulating synthetic OSIRIS images of the dust coma we showed how the dynamical and optical properties of the dust can be constrained. In particular, we showed how the dust-to-gas mass production rate ratio, Q_d/Q_g, the power-law exponent, q, the fraction of dust fall back, $Q_d^{fb}$, and the scattering properties are inter-related and constrain each other. Because these parameters are not independent they need to be fit simultaneously. Example, the lowest mass needed to match the brightness of the dust coma as observed by OSIRIS is achieved with power-law distributions with exponents between 4 and 4.5. Using the constraint of the total mass loss of the comet during the 2015 apparition we were able to show that only a narrow parameter set fits all observations. We determined that power-laws with $q=3.7_{-0.078}^{+0.57}$ are consistent with the data. This results in a total of $5.1_{-4.9}^{+6.0}·10^9$ kg of dust being ejected from the nucleus surface, of which $4.4_{-4.2}^{+4.9}·10^9$ kg escape to space and $6.8_{-6.8}^{+11}·10^8$ kg (or an equivalent of $14_{-14}^{+22}$ cm over the smooth regions) is re-deposited on the surface. This leads to a dust-to-gas ratio of $0.73_{-0.70}^{+1.3}$ for the escaping material and $0.84_{-0.81}^{+1.6}$ for the ejected material. Further, the minimum dust size must be strictly smaller than ~30μm and nominally even smaller than ~12μm. We have found that these results are robust with respect to varying the smallest dust size between 0.01 and 1μm and variations in the bulk density of the dust between 250− −1, 000 kg m ⁻³.

It remains an open question as to how dust particles are lifted/ejected from cometary surfaces [see e.g., 50]. Furthermore, a more detailed study of the change in the dust size distribution with heliocentric distance would be of great interest and could refine the work presented here. Finally, comprehensive work on estimating the amount of dust deposition through e.g., local digital terrain modeling [e.g., method by 49] would provide valuable additional constraints.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author. The data of some figures is available on www.spaceMarschall.net.

Author Contributions

RM performed the modeling of the gas and dust coma as well as the comparisons with ROSINA/COPS and OSIRIS and the analysis related to that. JM performed the calculations of the scattering properties of the dust particles. S-BG implemented new features that optimize DRAG3D. OP-R produced the unstructured simulation grid within which all UltraSPARTS and DRAG3D were run. NT wrote the IDL programs for reading and analysis of the OSIRIS images. J-SW provided support for running and optimization of UltraSPARTS. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

RM acknowledges the support from the Swiss National Science Foundation grant 184482.

JM acknowledges the support from ERC Grant No. 757390.

The team from the University of Bern is supported through the Swiss National Science Foundation, and through the NCCR PlanetS.

We thank Frank Preusker and Frank Scholten for providing us the comet shape model SHAP7 [2] used in this work.

We thank Vladimir Zakharov for providing valuable comments on the section of the analytical solution.

We acknowledge the personnel at ESA's European Space Operations Center (ESOC) in Darmstadt, Germany, European Space Astronomy Center (ESAC) in Spain, and at ESA for the making the Rosetta mission possible. Furthermore, we thank the OSIRS and ROSINA instrument and science teams for their hard work. We thank Martin Rubin and Kathrin Altwegg for giving us access and support to/for the ROSINA/COPS data.

Footnote

1. ^http://www.plasmati.tw

References

1. Sierks H, Barbieri C, Lamy PL, Rodrigo R, Koschny D, Rickman H, et al. On the nucleus structure and activity of comet 67P/Churyumov-Gerasimenko. Science. (2015) 347:1044. doi: 10.1126/science.aaa1044

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Preusker F, Scholten F, Matz KD, Roatsch T, Hviid SF, Mottola S, et al. The global meter-level shape model of comet 67P/Churyumov-Gerasimenko. Astron Astrophys. (2017) 601:L1. doi: 10.1051/0004-6361/201731798

CrossRef Full Text

3. Thomas N, Sierks H, Barbieri C, Lamy PL, Rodrigo R, Rickman H, et al. The morphological diversity of comet 67P/Churyumov-Gerasimenko. Science. (2015) 347:440. doi: 10.1126/science.aaa0440

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Pätzold M, Andert TP, Hahn M, Barriot JP, Asmar SW, Häusler B, et al. The Nucleus of comet 67P/Churyumov-Gerasimenko - Part I: the global view - nucleus mass, mass-loss, porosity, and implications. Mthly Notices R Astron Soc. (2019) 483:2337–46. doi: 10.1093/mnras/sty3171

CrossRef Full Text | Google Scholar

5. Fougere N, Altwegg K, Berthelier JJ, Bieler A, Bockelée-Morvan D, Calmonte U, et al. Three-dimensional direct simulation Monte-Carlo modeling of the coma of comet 67P/Churyumov-Gerasimenko observed by the VIRTIS and ROSINA instruments on board Rosetta. Astron Astrophys. (2016) 588:A134. doi: 10.1051/0004-6361/201527889

CrossRef Full Text | Google Scholar

6. Läuter M, Kramer T, Rubin M, Altwegg K. Surface localization of gas sources on comet 67P/Churyumov-Gerasimenko based on DFMS/COPS data. Mthly Notices R Astron Soc. (2018) 483:852–61. doi: 10.1093/mnras/sty3103

CrossRef Full Text | Google Scholar

7. Combi M, Shou Y, Fougere N, Tenishev V, Altwegg K, Rubin M, et al. The surface distributions of the production of the major volatile species, H2O, CO2, CO and O2, from the nucleus of comet 67P/Churyumov-Gerasimenko throughout the Rosetta Mission as measured by the ROSINA double focusing mass spectrometer. Icarus. (2020) 335:113421. doi: 10.1016/j.icarus.2019.113421

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Migliorini A, Piccioni G, Capaccioni F, Filacchione G, Bockelée-Morvan D, Erard S, et al. Water and carbon dioxide distribution in the 67P/Churyumov-Gerasimenko coma from VIRTIS-M infrared observations. Astron Astrophys. (2016) 589:A45. doi: 10.1051/0004-6361/201527661

CrossRef Full Text | Google Scholar

9. Bockelée-Morvan D, Crovisier J, Erard S, Capaccioni F, Leyrat C, Filacchione G, et al. Evolution of CO2, CH4, and OCS abundances relative to H2O in the coma of comet 67P around perihelion from Rosetta/VIRTIS-H observations. Mthly Notices R Astron Soc. (2016) 462(Suppl_1):S170–83. doi: 10.1093/mnras/stw2428

CrossRef Full Text | Google Scholar

10. Marshall DW, Hartogh P, Rezac L, von Allmen P, Biver N, Bockelée-Morvan D, et al. Spatially resolved evolution of the local H₂O production rates of comet 67P/Churyumov-Gerasimenko from the MIRO instrument on Rosetta. Astron Astrophys. (2017) 603:A87. doi: 10.1051/0004-6361/201730502

CrossRef Full Text | Google Scholar

11. Biver N, Bockelée-Morvan D, Hofstadter M, Lellouch E, Choukroun M, Gulkis S, et al. Long-term monitoring of the outgassing and composition of comet 67P/Churyumov-Gerasimenko with the Rosetta/MIRO instrument. Astron. Astrophys. (2019) 630:A19. doi: 10.1051/0004-6361/201834960

CrossRef Full Text | Google Scholar

12. Thomas N, El Maarry MR, Theologou P, Preusker F, Scholten F, Jorda L, et al. Regional unit definition for the nucleus of comet 67P/Churyumov-Gerasimenko on the SHAP7 model. Planet Space Sci. (2018) 164:19–36. doi: 10.1016/j.pss.2018.05.019

CrossRef Full Text | Google Scholar

13. Thomas N, Davidsson B, El-Maarry MR, Fornasier S, Giacomini L, Gracia-Berná AG, et al. Redistribution of particles across the nucleus of comet 67P/Churyumov-Gerasimenko. Astron Astrophys. (2015) 583:A17. doi: 10.1130/abs/2016AM-281342

CrossRef Full Text | Google Scholar

14. Keller HU, Barbieri C, Lamy P, Rickman H, Rodrigo R, Wenzel KP, et al. OSIRIS the scientific camera system onboard Rosetta. Space Sci Rev. (2007) 128:433–506. doi: 10.1007/s11214-006-9128-4

CrossRef Full Text | Google Scholar

15. Balsiger H, Altwegg K, Bochsler P, Eberhardt P, Fischer J, Graf S, et al. Rosina Rosetta orbiter spectrometer for ion and neutral analysis. Space Sci Rev. (2007) 128:745–801. doi: 10.1007/s11214-006-8335-3

CrossRef Full Text | Google Scholar

16. Marschall R, Su CC, Liao Y, Thomas N, Altwegg K, Sierks H, et al. Modelling observations of the inner gas and dust coma of comet 67P/Churyumov-Gerasimenko using ROSINA/COPS and OSIRIS data: first results. Astron Astrophys. (2016) 589:A90. doi: 10.1051/0004-6361/201628085

CrossRef Full Text | Google Scholar

17. Marschall R, Mottola S, Su CC, Liao Y, Rubin M, Wu JS, et al. Cliffs versus plains: Can ROSINA/COPS and OSIRIS data of comet 67P /Churyumov-Gerasimenko in autumn 2014 constrain inhomogeneous outgassing? Astron Astrophys. (2017) 605:A112. doi: 10.1051/0004-6361/201730849

CrossRef Full Text | Google Scholar

18. Gerig SB, Marschall R, Thomas N, Bertini I, Bodewits D, Davidsson B, et al. On deviations from free-radial outflow in the inner coma of comet 67P /Churyumov-Gerasimenko. Icarus. (2018) 311:1–22. doi: 10.1016/j.icarus.2018.03.010

CrossRef Full Text | Google Scholar

19. Marschall R, Rezac L, Kappel D, Su CC, Gerig SB, Rubin M, et al. A comparison of multiple Rosetta data sets and 3D model calculations of 67P/Churyumov-Gerasimenko coma around equinox (May 2015). Icarus. (2019) 328:104–26. doi: 10.1016/j.icarus.2019.02.008

CrossRef Full Text | Google Scholar

20. Su CC. Parallel Direct Simulation Monte Carlo (DSMC) Methods for Modeling Rarefied Gas Dynamics. Taiwan: National Chiao Tung University (2013).

21. Wu JS, Lian YY. Parallel three-dimensional Direct Simulation Monte Carlo method and its applications. Comput Fluids. (2003) 32:1133–60. doi: 10.1016/S0045-7930(02)00083-X

CrossRef Full Text | Google Scholar

22. Wu JS, Tseng KC, Wu FY. Parallel three-dimensional DSMC method using mesh refinement and variable time-step scheme. Comput Phys Commun. (2004) 162:166–87. doi: 10.1016/j.cpc.2004.07.004

CrossRef Full Text | Google Scholar

23. Wu JS, Tseng KC. Parallel DSMC method using dynamic domain decomposition. Int J Num Methods Eng. (2005) 63:37–76. doi: 10.1002/nme.1232

CrossRef Full Text | Google Scholar

24. Bieler A, Altwegg K, Balsiger H, Berthelier JJ, Calmonte U, Combi M, et al. Comparison of 3D kinetic and hydrodynamic models to ROSINA-COPS measurements of the neutral coma of 67P/Churyumov-Gerasimenko. Astron Astrophys. (2015) 583:A7. doi: 10.1051/0004-6361/201526178

CrossRef Full Text | Google Scholar

25. Zakharov V, Crifo JF, Rodionov AV, Rubin M, Altwegg K. The near-nucleus gas coma of comet 67P/Churyumov-Gerasimenko prior to the descent of the surface lander PHILAE. Astron Astrophys. (2018) 618:A71. doi: 10.1051/0004-6361/201832883

CrossRef Full Text | Google Scholar

26. Tzou CY. Calibrations of ROSINA-COPS and Observations at Comet 67P/Churyumov-Gerasimenko. Universitat Bern, Bern (2017).

Google Scholar

27. Finklenburg S, Thomas N, Su CC, Wu JS. The spatial distribution of water in the inner coma of Comet 9P/Tempel 1: comparison between models and observations. Icarus. (2014) 236:9–23. doi: 10.1016/j.icarus.2014.03.032

CrossRef Full Text | Google Scholar

28. Marschall R, Liao Y, Thomas N, Wu J. Limitations in the determination of surface emission distributions on comets through modelling of observational data - a case study based on Rosetta observations. Icarus. (2019) 346:113742. doi: 10.1016/j.icarus.2020.113742

CrossRef Full Text | Google Scholar

29. Mannel T, Bentley MS, Boakes PD, Jeszenszky H, Ehrenfreund P, Engrand C, et al. Dust of comet 67P/Churyumov-Gerasimenko collected by Rosetta/MIDAS: classification and extension to the nanometer scale. Astron Astrophys. (2019) 630:A26. doi: 10.1051/0004-6361/201834851

CrossRef Full Text | Google Scholar

30. Della Corte V, Rotundi A, Zakharov V, Ivanovski S, Palumbo P, Fulle M, et al. GIADA microbalance measurements on board Rosetta: submicrometer- to micrometer-sized dust particle flux in the coma of comet 67P/Churyumov-Gerasimenko. Astron Astrophys. (2019) 630:A25. doi: 10.1051/0004-6361/201834912

CrossRef Full Text | Google Scholar

31. Zakharov VV, Ivanovski SL, Crifo JF, Della Corte V, Rotundi A, Fulle M. Asymptotics for spherical particle motion in a spherically expanding flow. Icarus. (2018) 312:121–127. doi: 10.1016/j.icarus.2018.04.030

CrossRef Full Text | Google Scholar

32. Muinonen K, Markkanen J, Väisänen T, Peltoniemi J, Penttilä A. Multiple scattering of light in discrete random media using incoherent interactions. Optics Lett. (2018) 43:683. doi: 10.1364/OL.43.000683

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Markkanen J, Väisänen T, Penttilä A, Muinonen K. Scattering and absorption in dense discrete random media of irregular particles. Optics Lett. (2018) 43:2925. doi: 10.1364/OL.43.002925

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Markkanen J, Agarwal J, Väisänen T, Penttilä A, Muinonen K. Interpretation of the phase functions measured by the OSIRIS instrument for comet 67P/Churyumov-Gerasimenko. Astrophys J Lett. (2018) 868:L16. doi: 10.3847/2041-8213/aaee10

CrossRef Full Text | Google Scholar

35. Markkanen J, Agarwal J. Scattering, absorption and thermal emission by large cometary dust particles: synoptic numerical solution. Astron Astrophys. (2019) 631:A164. doi: 10.1051/0004-6361/201936235

CrossRef Full Text | Google Scholar

36. Dorschner J, Begemann B, Henning T, Jaeger C, Mutschke H. Steps toward interstellar silicate mineralogy. II. Study of Mg-Fe-silicate glasses of variable composition. Astron Astrophys. (1995) 300:503.

Google Scholar

37. Jäger C, Mutschke H, Henning T. Optical properties of carbonaceous dust analogues. Astron Astrophys. (1998) 332:291–9.

Google Scholar

38. Bertini I, La Forgia F, Tubiana C, Güttler C, Fulle M, Moreno F, et al. The scattering phase function of comet 67P/Churyumov-Gerasimenko coma as seen from the Rosetta/OSIRIS instrument. Mthly Notices R Astron Soc. (2017) 469:S404–15. doi: 10.1093/mnras/stx1850

CrossRef Full Text | Google Scholar

39. Fornasier S, Hasselmann PH, Barucci MA, Feller C, Besse S, Leyrat C, et al. Spectrophotometric properties of the nucleus of comet 67P/Churyumov-Gerasimenko from the OSIRIS instrument onboard the ROSETTA spacecraft. Astron Astrophys. (2015) 583:A30. doi: 10.1051/0004-6361/201525901

CrossRef Full Text | Google Scholar

40. Feller C, Fornasier S, Hasselmann PH, Barucci A, Preusker F, Scholten F, et al. Decimetre-scaled spectrophotometric properties of the nucleus of comet 67P/Churyumov-Gerasimenko from OSIRIS observations. Mthly Notices R Astron Soc. (2016) 462:S287–303. doi: 10.1093/mnras/stw2511

CrossRef Full Text | Google Scholar

41. Masoumzadeh N, Oklay N, Kolokolova L, Sierks H, Fornasier S, Barucci MA, et al. Opposition effect on comet 67P/Churyumov-Gerasimenko using Rosetta-OSIRIS images. Astron Astrophys. (2017) 599:A11. doi: 10.1051/0004-6361/201629734

CrossRef Full Text | Google Scholar

42. Laor A, Draine BT. Spectroscopic constraints on the properties of dust in active galactic nuclei. Astrophys J. (1993) 402:441. doi: 10.1086/172149

CrossRef Full Text | Google Scholar

43. Bockelée-Morvan D, Rinaldi G, Erard S, Leyrat C, Capaccioni F, Drossart P, et al. Comet 67P outbursts and quiescent coma at 1.3 au from the Sun: dust properties from Rosetta/VIRTIS-H observations. Mthly Notices R Astron Soc. (2017) 469(Suppl_2):S443–58. doi: 10.1093/mnras/stx1950

CrossRef Full Text | Google Scholar

44. Merouane S, Stenzel O, Hilchenbach M, Schulz R, Altobelli N, Fischer H, et al. Evolution of the physical properties of dust and cometary dust activity from 67P/Churyumov-Gerasimenko measured in situ by Rosetta/COSIMA. Mthly Notices R Astron Soc. (2017) 469(Suppl_2):S459–74. doi: 10.1093/mnras/stx2018

CrossRef Full Text | Google Scholar

45. Fulle M, Marzari F, Della Corte V, Fornasier S, Sierks H, Rotundi A, et al. Evolution of the dust size distribution of comet 67P/Churyumov-Gerasimenko from 2.2 au to Perihelion. Astrophys J. (2016) 821:19. doi: 10.3847/0004-637X/821/1/19

CrossRef Full Text | Google Scholar

46. Moreno F, Muñoz O, Gutiérrez PJ, Lara LM, Snodgrass C, Lin ZY, et al. The dust environment of comet 67P/Churyumov-Gerasimenko: results from Monte Carlo dust tail modelling applied to a large ground-based observation data set. Mthly Notices R Astron Soc. (2017) 469(Suppl_2):S186–94. doi: 10.1093/mnras/stx1424

CrossRef Full Text | Google Scholar

47. Della Corte V, Rotundi A, Fulle M, Ivanovski S, Green SF, Rietmeijer FJM, et al. (2016) 67P/C-G inner coma dust properties from 2.2 au inbound to 2.0 au outbound to the Sun. Mthly Notices R Astron Soc. 462:S210–9. doi: 10.1093/mnras/stw2529

CrossRef Full Text | Google Scholar

48. Ott T, Drolshagen E, Koschny D, Güttler C, Tubiana C, Frattin E, et al. Dust mass distribution around comet 67P/Churyumov-Gerasimenko determined via parallax measurements using Rosetta's OSIRIS cameras. Mthly Notices R Astron Soc. (2017) 469:S276–84. doi: 10.1093/mnras/stx1419

CrossRef Full Text | Google Scholar

49. Jorda L, Gaskell R, Capanna C, Hviid S, Lamy P, Durech J, et al. The global shape, density and rotation of Comet 67P/Churyumov-Gerasimenko from preperihelion Rosetta/OSIRIS observations. Icarus. (2016) 277:257–78. doi: 10.1016/j.icarus.2016.05.002

CrossRef Full Text | Google Scholar

50. Vincent J.-B., Farnham T, Kuhrt E, Skorov Y, Marschall R, Oklay N, et al. Local manifestations of cometary activity. Space Sci Rev. (2019) 215:30. doi: 10.1007/s11214-019-0596-8

CrossRef Full Text | Google Scholar