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 \pm 1.5) \cdot 10^9$~kg during the 2015 apparition. Further, we found that power-laws with $q=3.7^{+0.57}_{-0.078}$ are consistent with the data. This results in a total of $5.1^{+6.0}_{-4.9}\cdot10^9$~kg of dust being ejected from the nucleus surface, of which $4.4^{+4.9}_{-4.2}\cdot10^9$~kg escape to space and $6.8^{+11}_{-6.8}\cdot10^8$~kg (or an equivalent of $14^{+22}_{-14}$~cm over the smooth regions) is re-deposited on the surface. This leads to a dust-to-gas ratio of $0.73^{+1.3}_{-0.70}$ for the escaping material and $0.84^{+1.6}_{-0.81}$ for the ejected material. We have further found that the smallest dust size must be strictly smaller than $\sim30\mu$m and nominally even smaller than $\sim12\mu$m.


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 travelled out towards 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 (Sierks et al., 2015;Preusker et al., 2017) and diverse morphology (Thomas et al., 2015b). 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) (Pätzold et al., 2019). Second, the total volatile mass loss which can be indirectly determined by the in-situ measurements of the gas density (e.g. Fougere et al., 2016;Läuter et al., 2018;Combi et al., 2020) or remote sensing data (e.g. Migliorini, A. et al., 2016;Bockelée-Morvan et al., 2016;Marshall et al., 2017;Biver et al., 2019). 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 ) that indicate that possibly a large fraction of the ejected dust is re-deposited (Thomas et al., 2015a). 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, Keller et al., 2007) 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 modelled the inner gas and dust comae for the entire Rosetta mission. We use Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA, Balsiger et al., 2007) 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 Sec. 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 Sec. 3. We will discuss the results of our work in Sec. 4 and summarise and conclude our work in Sec. 5.

METHOD
In this work, we have used the modelling approach (and in particular our DRAG3D model for the dust coma) described in detail in Marschall et al. (2016). 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 (Marschall et al., 2016(Marschall et al., , 2017Gerig et al., 2018;Marschall et al., 2019b). 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 modelling elements and refer to Marschall et al. (2016) for a detailed description.

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 centre date of each epoch. Furthermore that the difference in sub-solar latitude be less than 5 • from the centre time of epoch to the start and end of the epoch, respectively. This leads to the twenty epochs listed in Tab. 1 and illustrated in Fig. 1. Simulations were run for the centre 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). The basis of all simulations is the 3D shape model by Preusker et al. (2017). 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 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 2 O ice surface that is areally mixed with inert refractory surface akin to a chequerboard 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 (Keller et al., 2007) 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 (Marschall et al., 2016) but will break down when a lot of dust is released. We will further discuss this limitation later on.

Gas kinetic simulations
The gas flow-field is calculated using the DSMC technique. The code used is called UltraSPARTS 1 and is a commercialized derivative of the PDSC++ code (Su, 2013) used in previous papers (e.g Marschall et al., 2016Marschall et al., , 2017. 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 (Wu and Lian, 2003;Wu et al., 2004;Wu and Tseng, 2005) 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 Gerig et al. (2018); Marschall et al. Here we simulate the full 3D gas flow up to a distance of 10 km from the nucleus centre.
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 Bieler et al., 2015;Marschall et al., 2016;Fougere et al., 2016;Zakharov et al., 2018a) 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 behaviours. Because the EAF is a free parameter it needs to be constrained by data. In our case, we determine the EAF by comparing modelled densities extrapolated to the Rosetta position and actual COmet Pressure Sensor (COPS; Balsiger et al. (2007)) 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. (2016). 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. Tab. 2 shows the EAF used and the resulting average global H 2 O production rate for one comet day.
The global gas production rate as a function of time is shown in Fig. 2 (purple band in top panel). Because we have used the total ROSINA/COPS data -which also contains the other gas species other than waterto 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 Tab. 5. The resulting integrated mass loss over the shown period adds up to (6.1 ± 1.5) · 10 9 kg. This is well in line with values published in other works as e.g. (6.3 ± 2.0) · 10 9 kg (Combi et al., 2020) or (5.8 ± 1.8) · 10 9 kg (Läuter et al., 2018). The error arises from the uncertainty of the data (up to 15%; Tzou (2017) although the relative errors are probably smaller (M. Rubin, pers. comm.)) and our model (5-10%; Finklenburg et al. (2014)) 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. (2019a) 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. (2018) 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 2 O or at least near-surface volatiles without a significant thermal lag. Figure 2. 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). Table 3. 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

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. (2016) 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 Tab. 1). The images used in this work, as well as some camera and geometric properties, are shown in Tab. 3. 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 centre 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 • .
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 −3 matching roughly the bulk density of the nucleus (Preusker et al., 2017). 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 (Mannel et al., 2019) 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. (2019) 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.24CJ −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 Sec. 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 Sec. 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 behaviour. Gerig et al. (2018) 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 (Zakharov et al., 2018b). 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.

Scattering model
Previously, we have used a spherical particle model and a Mie scattering code in our modelling pipeline. Here we use a much more sophisticated approach based on the recently introduced radiative transfer with reciprocal transactions framework (Muinonen et al., 2018;Markkanen et al., 2018b). 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-micrometre-sized organic grains and micrometre-sized silicate grains. Such a particle model has been found to be in good agreement with OSIRIS (Markkanen et al., 2018a) and VIRTIS (Markkanen and Agarwal, 2019) phase function measurements. The refractive index for silicate grains is assumed to be m = 1.6048692 + i0.0015341 corresponding to magnesium iron pyroxine (Dorschner et al., 1995) and for organic grains m = 1.55950 + i0.42964 corresponding to amorphous carbon (Jäger et al., 1998). At 612 nm (WAC filter 18) the resulting scattering phase functions for different particle sizes normalized to the geometric albedo are shown in Fig. 3. The figure shows good agreement of the phase function of large particles (> 1 cm) with the nucleus phase function as measured by OSIRIS (Fornasier et al., 2015;Feller et al., 2016;Masoumzadeh et al., 2017). 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 (Bertini et al., 2017) is between 10 and 100 µm. The numerical method of Markkanen et al. (2018a) 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 Marschall et al., 2016) matching the single scattering albedo of the Mie result with the approach of Markkanen et al. (2018a) 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.

THEORETICAL CONSIDERATION
To put some of our results in the next section (Sec. 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. (2018b) for a wide range of conditions. They defined and which are dimensionless parameters, where ρ d is the specific density of the dust particles, v g0 the gas velocity near the surface, v max g 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 characterises the ability of a dust particle to adjust to the gas velocity while F u quantifies the importance of gravity. Zakharov et al. (2018b)) 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) ∝ r −0.5 for small F u (i.e. if gravity plays a minor role). The asymptotic dust velocities are given by: 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 normalised mass distribution, f md , of particles ejected from the surface is 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 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 centre of the nucleus is: The column density at the distance from the centre of the nucleus in the image plane is: The total number of dust particles in a column within a circular observing aperture of radius is: The brightness is proportional to the flux F (W/m 2 ) gathered by an instrument which for an optically thin coma is: 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 Eqs. (3), (4) and (8) in (9) we get: For fixed F, Q g , R N , ρ d , v g0 , r * , v d (r * ), ∆ and F (r, ) = C 4 − q r 4−q max − r 4−q min χr 5 2 −q q sca (r)ϕ av (r)dr (11) where C = 3 32 FQg C Iv ρ d ∆ 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: q sca (r) =    0.233 · 10 24 · r 7/2 , 10 −8 ≤ r < 2 · 10 −7 4.993 · 10 −5 · r −2/3 , 2 · 10 −7 ≤ r < 2 · 10 −4 0.02, 2 · 10 −4 ≤ r ≤ 1.0 (12) Figure 4 shows the computed q sca from Sec. 2.4, the fitted q f it 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.
Under the assumption we made the integration of Eq. (11) becomes trivial. Figure 5. Dust-to-gas mass production rate ratio vs. power-law exponent for constant dust brightness.

Frontiers
For a given gas production rate Q g the maximum dust size a max is also constant. Fig. 5 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.
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 (Fig. 5) 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.

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: where q is the differential power-law exponent. Fig. 6 illustrates an example of an OSIRIS and synthetic image for epoch 12 (solstice). As described in Sec. 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. For a given gas production rate, the three major factors controlling the brightness (see Eq. (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 Sec. 2.4 and Markkanen et al., 2018a). 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. Figure 7. 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 colours of the curves indicate the global gas production rate of the respective epoch.
As has already been shown in Fig. 12 of Marschall et al. (2016) 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. Fig. 7 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.
Comparing Fig. 7 to the analytical solution presented in Fig. 5 of Sec. 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 Fig. 12 of Marschall et al. (2016) the Q d /Q g values we find here (in particular for q < 3.5) are much higher while the behaviour 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. (2016) assumed the scattering properties of astronomical silicate (Laor and Draine, 1993) 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 Fig. 7 are worth discussing. First, we have assumed that all dust particles have a bulk density equal to the nucleus density (533 kg m −3 ). 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 −3 and 1000 kg m −3 . The two left panels of Fig. 8 show the results for a moderate activity environment (epoch 6 -inbound equinox, Q g = 35 kg s −1 ) and a high activity environment (epoch 12 -solstice, Q g = 800 kg s −1 ). 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.
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 Mannel et al., 2019) and there is indirect evidence of sub-micron particles observed by VIRTIS during outbursts . We have thus explored the range of the smallest sizes between 0.01µm and 1µm. The two right panels of Fig. 8 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 Fig. 4) 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 Fig. 9 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 Fig. 7. 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 9. The ratio of dust mass falling back, Q D f b 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. Fig. 10 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 −1 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.
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. E.g. 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 Figure 10. 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 2 year mission comet 67P had lost (10.5 ± 3.4) · 10 9 kg (Pätzold et al., 2019). The total mass loss, M tot = (10.5 ± 3.4) · 10 9 kg, is: where M g is the total volatile mass loss, and M esc d is the total escaping dust mass. For the dust masses, we can further specify that where M d is the total dust mass ejected from the nucleus, and M f b d 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 9 kg. Combined with the total mass loss of the nucleus it follows that M esc d = (4.4 +4.9 −4.4 ) · 10 9 kg. Note that within this interval exists the possibility that M esc d = 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 11 shows the M esc d as a function of the power-law exponent, q. Cases, where M esc d exceeds the nominal dust mass loss of 4.4 · 10 9 kg (horizontal dashed green line) or the maximum dust mass loss of 9.3 · 10 9 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 11 also illustrates the effect of the smallest size -discussed earlier in more detail for an individual OSIRIS imageon 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 Fig. 8 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. Fig. 11 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 summarised in Tab. 4. We have also determined the deposition height, H, that results if the fallback material is spread equally on the smooth deposits (9.43 km 2 ) identified by Thomas et al. (2018). Figure 11. The total escaping dust mass, M esc d , 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).
The results in Tab. 4 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. Table 4. Power-law exponents, dust mass loss, dust-to-gas ratio, and fall back as a function of the smallest dust size (see Fig. 11) . Dashed entries mean that no solution is possible for this size.
This analysis assumes that the dust size distribution does not change along the orbit. There is an indication (e.g. Merouane et al., 2017) 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  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 , q = 3.8 by COSIMA (Merouane et al., 2017), 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 .
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 Della Corte et al. (2016) 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 −9 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 −1 . 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 Fig. 5 in Della Corte et al., 2016). At a phase angle of 90 • our model dust particles with radius 100 µm have speeds of (3.5 ± 1.0) m s −1 at the inbound equinox (epoch 6) and (18 ± 5.6) m s −1 at the summer solstice (epoch 12). For phase angles of 60 • the model dust particles have speeds of (7.0 ± 1.2) m s −1 at the inbound equinox and (32 ± 5.2) m s −1 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 −1 ) is roughly an order of magnitude lower than those reported by e.g. Moreno et al. (2017) (∼ 3000 kg s −1 ) or Ott et al. (2017) (∼ 8300 kg s −1 ). Furthermore, Moreno et al. (2017) report a total dust mass loss of 1.4 · 0 10 kg. As neither Moreno et al. (2017) nor Ott et al. (2017) 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 (Pätzold et al., 2019) taking into account the estimates of the volatile mass loss in this work and others (Combi et al., 2020;Läuter et al., 2018).
Finally, the determination of the power-law exponent allows us to determine the dust production rate (Fig. 2, top panel, green band), dust-to-gas ratio (Fig. 2, bottom panel, green line), and fraction of fallback (Fig. 2, 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 towards 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 behaviour 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 . 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.

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 9 kg. This is in line with values published in other works as e.g. (6.3 ± 2.0) · 10 9 kg (Combi et al., 2020) or (5.8 ± 1.8) · 10 9 kg (Läuter et al., 2018). 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 f b d , and the scattering properties are inter-related and constrain each other. Because these parameters are not independent they need to be fit simultaneously. E.g. 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.57 −0.078 are consistent with the data. This results in a total of 5.1 +6.0 −4.9 · 10 9 kg of dust being ejected from the nucleus surface, of which 4.4 +4.9 −4.2 · 10 9 kg escape to space and 6.8 +11 −6.8 · 10 8 kg (or an equivalent of 14 +22 −14 cm over the smooth regions) is re-deposited on the surface. This leads to a dust-to-gas ratio of 0.73 +1.3 −0.70 for the escaping material and 0.84 +1.6 −0.81 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 − 1µm and variations in the bulk density of the dust between 250 − 1000 kg m −3 .
It remains an open question as to how dust particles are lifted/ejected from cometary surfaces (see e.g. Vincent et al., 2019). 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 modelling (e.g. method by Jorda et al., 2016) would provide valuable additional constraints.

AUTHOR CONTRIBUTIONS
RM performed the modelling 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. SBG implemented new features that optimise DRAG3D. OPR 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. JW provided support for running and optimisation of UltraSPARTS.

FUNDING
Raphael Marschall acknowledges the support from the Swiss National Science Foundation grant 184482. Johannes Markkanen 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. Table 5. The mean gas production rate,q g [kg s −1 ] as a function of ephemeris time (ET): q g (ET ) = a · ET 2 + b · ET + c epoch a b c