© ESO, 2016

1. Introduction

The European Space Agency’s Rosetta spacecraft entered an irregular orbit about the Jupiter family comet 67P/Churyumov-Gerasimenko (hereafter 67P) on 6 August 2014. Rosetta carries 11 instruments for investigating the comet and its environment (Glassmeier et al. 2007; Schulz 2010). In particular, the payload has been designed to comprehensively investigate the gas and dust emission from the nucleus as the comet approaches the Sun. In this work, we focus on making a first comparison between observations of the nucleus and the innermost dust coma using the Optical, Spectroscopic, and Infrared Remote Imaging System (OSIRIS; Keller et al. 2007) and observations of the spacecraft in situ gas density using the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA; Balsiger et al. 2007). Our goal in this work is to use the observations of the nucleus, the innermost dust coma, and the local gas density to constrain the outgassing distribution on the nucleus surface and to provide a preliminary model of the innermost gas and dust coma of the comet when it was at 3.4 AU pre-perihelion from the Sun. In particular, we are looking for models in support of the supposed inhomogeneous surface ice distribution inferred from date from the Microwave Instrument for the Rosetta Orbiter – MIRO (Biver et al. 2015; Lee et al. 2015; Gulkis et al. 2015), the Ultraviolet Imaging Spectrometer – ALICE (Feldman et al. 2015), the Visible and Infrared Thermal Imaging Spectrometer – VIRTIS (Bockelée-Morvan et al. 2015), and ROSINA (Bieler et al. 2015) for the gas and reproducing the dust jet features observed by OSIRIS (Lara et al. 2015; Lin et al. 2015).

Fig. 1 Diagram of our modelling and verification approach.

Our modelling approach is illustrated in Fig. 1 and is similar to the work of Tenishev et al. (2011). Each tile represents the respective self-contained model. We begin by using the full 3D shape model derived from the OSIRIS observations of the nucleus. We compute the surface insolation taking the orientation of the nucleus with respect to the Sun into account. This requires using a ray-tracing algorithm because of the complex shape of the body (Sierks et al. 2015; Preusker et al. 2015) because shadowing plays an important role in determining the surface energy budget. This energy budget is then calculated by taking sublimation from the surface into account. We adopt here a simple approach that ignores thermal conductivity (Huebner et al. 2006). The thermal inertia of cometary surfaces has been shown in the past to be extremely low (e.g. Groussin et al. 2013) and preliminary indications suggest that 67P has similar surface thermal properties (Gulkis et al. 2015). As a result, this approximation is probably adequate as a first approach. However, we recognize that more sophisticated models may be necessary because conductivity into the interior at these heliocentric distances can produce strong relative reductions in gas emission rates.

The thermal model results in a gas production rate from each surface facet, which is reduced to observed global production rate values by using the concept of an active fraction for each facet. This is used, together with a gas temperature, to compute a gas flow field using a 3D parallel Direct Simulation Monte Carlo (DSMC) code. The modelling of the gas dynamics in the immediate vicinity of a cometary nucleus has made great strides forwards in recent years because the available computational power increased. It has been shown that the DSMC method converges to the Boltzmann equation in the limit of many simulation particles and collisions (Nanbu 1980; Wagner 1992). Combi & Smyth (1988) first studied the DSMC approach for comets, and later, 2D axially symmetric DSMC coma models were published by several authors (Zakharov et al. 2009; Tenishev et al. 2008, 2011; Skorov et al. 2004, 2006; Crifo et al. 2002, 2003, 2005).

Here we have followed approaches described by, for example, Tenishev et al. (2008) to produce the gas flow field from the initial boundary conditions. Once this flow field has been determined, it can be compared to measurements at the spacecraft by the ROSINA experiment through its COmet Pressure Sensor (COPS) instrument. This provides our first constraint on the model in this work, and we show that our results broadly agree with similar modelling efforts by Bieler et al. (2015). Following this, we computed trajectories of dust particles, distributed in size, emitted from the surface using a Monte Carlo approach combined with standard gas drag equations to derive accelerations and local dust densities in the innermost coma. Finally, we determined column densities of the dust when viewed from specific geometries and input this into a Mie scattering code (Bohren & Huffman 1983; Fink & Rubin 2012) to compute brightnesses of the inner coma.

These results can then be compared to OSIRIS observations of the dust coma to give a further constraint on the model.In our discussion of the results, we refer to selected regions of the nucleus that have been defined by Thomas et al. (2015b) and El-Maarry et al. (2015). In the following section, we describe the instruments and observations that are used for this study. In Sect. 3, we present the models we use (nucleus shape, thermal, gas dynamics, dust motion and dust scattering) in more detail. In Sect. 4, we describe the simulation results and how they indicate the need for inhomogeneous outgassing (i.e. not directly in relation to the insolation) and present a solution for such inhomogeneous outgassing. We conclude with a discussion of the simplifying assumptions made and the need for further constraints.

2. ROSETTA instruments considered for this study

2.1. ROSINA/COPS

The COmet Pressure Sensor (COPS) consists of two extractor-type ionization gauges combined with high sensitivity electrometers: the Nude Gauge (NG) measures the total neutral density, and the Ram Gauge (RG) derives the dynamic pressure of the outflowing gas. Combining the two gives the outflow velocity and the temperature of the expanding neutral gas coma (Balsiger et al. 2007). In nominal operation mode, the nude gauge obtains ten-second averages every minute in the operational pressure range of 10^-11–10^-5 mbar. It has been shown that the Rosetta spacecraft itself produces a background of neutral gas through desorption of volatiles, diffusion of trapped molecules, and spacecraft material decomposition (Schläppi et al. 2010), which is taken into account when deriving the neutral gas density in the coma. Furthermore, the background can also be affected on short time scales, when e.g. the spacecraft attitude changes and previously cold surfaces come into sunlight.

The data set we use was obtained between 21 August 2014 and 22 September 2014 when the comet was between 3.51 and 3.32 AU from the Sun. We have split this data set into a series of zones for more detailed comparisons. We have eliminated data acquired during major slews and thruster firings, both of which influence the local gas density measurements.

2.2. OSIRIS

OSIRIS consists of two high-resolution cameras, one a narrow angle camera (NAC) with an angular resolution of 18.6 μrad px^-1 and a field of view of 2.20 × 2.22°, and the other a wide angle camera (WAC) with 100 μrad px^-1 angular resolution and a field of view of 11.35 × 12.11° (Keller et al. 2007). Both cameras operate in the optical wavelength regime (250–1000 nm) with 12 filters for the NAC and 14 filters for the WAC. The images used for comparisons in this work were acquired with the WAC using a dust continuum filter with a central wavelength of 612.6 nm and a bandwidth of 9.8 nm.

3. Gas and dust models

3.1. Thermal model

3.1.1. Shape model and incident angle calculation

For this work we adopted the shape model SHAP4S of Preusker et al. (2015), which is based on a stereophotogrammetric technique. The original model has a horizontal spacing of 2 m and a vertical accuracy on the decimetre scale. With over 16 million facets, this shape model is currently beyond our computing capacity. We therefore used a decimated version of this model with 50 000 facets, providing a dimension of roughly 32 m per facet.

For each surface facet and for a given position of the Sun, we compute the angle of incidence of the light, i.e. the angle between the direction to the Sun and the respective surface normal. By means of ray tracing, we determine whether a facet is occluded by another part of the comet and is thus in shadow. At this stage we do not account for self-heating.

3.1.2. Energy balance, surface temperature, and gas production rate

A simple thermal model was constructed omitting thermal conductivity (i.e. the thermal inertia was set to zero. We note that Gulkis et al. (2015), Schloerb et al. (2015), and Choukroun et al. (2015) measured low values consistent with values seen at 9P/Tempel 1) but including sublimation of water ice. The sublimation coefficient was set to 1 for simplicity. The thermal balance for each surface facet was produced by $\mathrm{0}\mathrm{=}\frac{\mathit{S}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{A}}_{\mathrm{H}}\mathrm{)}\cos \mathrm{(}{\mathit{i}}^{\mathrm{)}}}{{{\mathit{R}}_{\mathrm{h}}}^{\mathrm{2}}}\mathrm{-}\mathit{ϵ\sigma }{\mathit{T}}^{\mathrm{4}}\mathrm{-}\mathit{L}{\mathit{Q}}_{\mathrm{gas}}\mathit{,}$(1)where A_H is the directional–hemispheric albedo (set to 0.04), S the solar constant at 1 astronomical unit (AU) (taken to be 1384 W m^-2), i the angle of incidence, R_h the heliocentric distance of the comet (set to 3.4 AU), ϵ the IR emissivity (set to 0.9), σ Stefan-Boltzmann’s constant, L the latent heat of sublimation of water ice (2.84 MJ kg^-1; Huebner et al. 2006), and Q_gas the sublimation rate (i.e. the gas mass loss rate per unit time and unit area). We note that Keller et al. (2015) presented a thermal model including self heating which, together with the observed thermal inertia, needs to be assessed in a later step.

The sublimation rate was then computed from the surface temperature, T, using the equation ${\mathit{Q}}_{\mathrm{gas}}\mathrm{=}{\mathit{p}}_{\mathrm{evp}}\sqrt{\frac{{\mathit{M}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}{\mathrm{2}\mathit{\pi kT}}}\mathit{,}$(2)where k is the Boltzmann constant, M_H2O the molecular mass of water, and the equilibrium vapour pressure of water vapour (p_evp) was computed from values given by Huebner et al. (2006). This scheme provided a sublimation flux and a gas temperature for each facet and limits the number of free parameters. Essentially only ϵ is a free parameter because varying A_H within reasonable bounds has only a small effect. For unilluminated surfaces, the gas flux was set to zero and the nominal surface temperature to 100 K. Using this scheme would normally produce gas production rates far in excess of what is currently observed. We therefore scale the results of the production rates to produce values that are closer to those observed at 67P. One can visualize this as being equivalent to only a fraction of each surface facet being active with the rest being inert (akin to a chessboard pattern; Keller et al. 2015). We refer to this as the effective active fraction with required values at this time of typically 1% (Gulkis et al. 2015) for the purely insolation driven models.

Assuming a half Maxwellian velocity distribution at the surface, we can convert the production rate per unit area, Q_gas, to an initial number density at the surface. By using ${\mathit{v}}_{\mathrm{g}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}\sqrt{\frac{\mathrm{8}\mathit{kT}}{\mathit{\pi }{\mathit{m}}_{\mathrm{g}}}}\mathit{,}$(3)for the speed of the gas, v_g, with m_g being the molecular mass of the gas, we can convert the mass production rate to a number density n_g (Bird 1994) due to flux conservation ${\mathit{Q}}_{\mathrm{gas}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}{\mathit{n}}_{\mathrm{g}}{\mathit{m}}_{\mathrm{g}}{\mathit{v}}_{\mathrm{g}}\mathit{.}$(4)This scheme implies that surface temperatures where sublimation is occurring never exceed approximately 200 K (i.e. the free sublimation temperature of water), which is formally inconsistent with measurements of the surface temperatures of comets (Emerich et al. 1987). However, the effective active fractions have been shown to be small (as confirmed here), and this approach does not exclude a hotter surface being adjacent to our active surface.

For the simulations presented here, we assume the gas equilibrates with the subliming surface temperature at the source. This is an intermediate assumption between two extremes. On the one hand, the energy in the gas may be lower than this assumption if the rotational degrees of freedom are in disequilibrium after emission from the icy surface. On the other hand, if the gas passes through a hot, inert surface layer, then it will be heated above the sublimation temperature before leaving the surface.

3.2. DSMC gas model

3.2.1. DSMC method

Our preferred method for gas dynamics simulations of the coma is direct simulation Monte Carlo DSMC (Bird 1994). Our code for modelling the coma, named PDSC⁺⁺ (Su 2013), has previously been used to model the water vapour distribution in the vicinity of comet 9P/Tempel 1 (Finklenburg et al. 2014) and is based on the PDSC code developed by Wu and co-workers (Wu & Lian 2003; Wu et al. 2004; Wu & Tseng 2005). PDSC⁺⁺ allows a simulation of 2D, 2D-axisymmetric, and 3D flows on hybrid unstructured grids. The code has been parallelized, allowing a much larger number of cells and has been implemented on several clusters in Bern (Switzerland) and Taiwan. The code is especially useful in that it is able to treat the large density gradients by implementation of a variable time step and a transient adaptive sub-cell technique to increase computational speed and accuracy in the regions of high density (Finklenburg et al. 2014). The DSMC code represents the real gas molecules by a smaller number of simulation particles that are used to generate a statistical representation of the flow.

The microscopic behaviour of the simulation molecules is separated into a translational step and a collision step. In the collision step, the collision pairs are chosen randomly from all the molecules that are within one cell at the time of the collision step. Macroscopic gas properties such as number density, velocity, and temperature are then calculated by averaging the appropriate quantity over all simulation molecules in the respective sampling cell.

For the simulation to yield reliable results it is important that in each cell, the mean collisional separation (mcs, i.e. the mean distance between the collision of two simulation particles) divided by the physical mean free path (mfp) be less than one. Additionally, the total number of simulation particles per cell should be at least 10 to allow meaningful sampling of the gas properties. These constraints have usually been met with the simulations shown here.

3.2.2. The unstructured grid and simulation domain

A substantial amount of time goes into building the grid. On the basis of the shape model described in Sect. 3.1.1, we constructed an unstructured 3D grid with the software Gridgen^TM by Pointwise¹. Since we are primarily interested in the gas flow immediately above the nucleus, the domain has been extended out to only 10 km from the centre of the nucleus. The unstructured grid uses tetrahedron cells, and the cell size must be adjusted to satisfy the condition that mcs/mfp< 1, which requires some experience. The grid used for the simulations in this study consists of just below one million cells that are then split automatically during runtime into much smaller sub-cells in regions with short mfp close to the active part of the surface. It must be noted that the size of the cells increases with the distance from the surface to compensate for the reduction in gas density caused by the expansion of the gas. The surface cell size is close to that of the decimated version of the shape model (see above). The cell size contributes to the error bar in determining the gas density at the spacecraft shown later in our results.

3.2.3. Input parameters and assumptions

For the inlet boundary condition of the gas model, we set the surface temperature and number densities calculated according to Sect. 3.1.2. The outer boundary is taken to be a vacuum boundary. In the code we can set each inlet boundary facet to be either absorbing or reflecting. For the reflecting surface option, we have the additional choice of setting specular or diffuse reflection where, in the latter case, the gas hitting the comet surface is thermalized to the respective surface facet temperature and is then re-emitted. We have performed simulations with both specular and diffuse reflection and found that the final result is not significantly affected by this choice as far as local gas densities and gas column densities are concerned. As these are the initial results that we seek to compare with ROSINA/COPS data, we only present the diffuse reflection results because this is, to our mind, the more physically appropriate representation of a reflecting surface, i.e. the assumption that any surface can be considered rough on the level of the size of individual atoms or molecules. Obviously a temperature-dependent behaviour of the surface reflectance seems the most physically plausible. Cold surfaces could thermally trap gas and thus act as absorbing surfaces (Rubin et al. 2014), and hot surfaces would diffusely reflect any gas back flux. This can be included in our code but has not yet been studied, especially as we currently lack the data that would be able to differentiate between all of these cases. We have therefore taken the approach of studying the two extreme cases of total absorption and total diffuse reflection at the inlet boundary.

To simplify the model we have made some additional assumptions. First, we neglect the rotation of the nucleus and thus consider only steady-state solutions for both the gas and the dust. In the case of the gas flow, this assumption is easily justifiable because the high velocities of the gas (~ 500 m s^-1) result in the gas crossing our simulation domain within seconds to minutes, which is much less than the rotation period of the comet (12.4 h, Mottola et al. 2014). For the dust this is equally true only for the small and thus fast particles. For large particles (>300 μm), a time-dependent solution is necessary and will have to be dealt with at a later stage. On the other hand, in the case of a power-law distribution for the dust number densities with power-law indices, b, larger than 3.5, we see in Sect. 4.3.1 that the contribution to the observed brightness from large particles that are affected by this is very minor.

Second, we assume that although the dust is accelerated by the gas, the back reaction of the dust onto the gas is negligible, allowing us to model the gas and dust dynamics separately and in series. Tenishev et al. (2011) note that heating of the gas by collisions with the dust should be negligible and that kinetic back coupling of the dust in the coma does not affect the gas field in the case of 67P. Our models show the total kinetic energy of the dust to be two orders of magnitude lower than that of the gas. We thus agree with this statement.

Third, the initial velocity distribution function (VDF) of the gas at the surface is set to be a half Maxwellian. Although it is conceivable that the initial flow is more collimated, the effects of this as seen by the in situ measurements at the spacecraft position are unclear. Our preliminary investigations regarding this show a variance of a few percent at a distance of 10 km from the nucleus centre when comparing a half Maxwellian VDF with a cosine law with a high n if $\mathit{f}{}^{\mathrm{\right(}}\mathit{\theta }^{\mathrm{\left)}}\mathrm{\propto }{\cos }^{\mathit{n}}{}^{\mathrm{\right(}}\mathit{\theta }^{\mathrm{\left)}}$ (Liao et al. 2015). Furthermore, we consider water as the only gas species for this investigation. We are aware of the high relative abundance of carbon dioxide (particularly above southern latitudes, i.e. the winter hemisphere at the time of the observations discussed here) in the ROSINA dataset (Hässig et al. 2015), but we draw attention to the fact that H₂O was always the dominant species at the spacecraft for positive spacecraft latitudes in the time span considered in this study.

3.3. Dust model

3.3.1. Dust tracking model

The motion of the dust is driven by two forces: the force produced by the gas drag on the dust particles to accelerate them away from the surface, and the gravity of the nucleus as an opposing force. Solar radiation pressure is assumed to be negligible this close to the nucleus (see also Tenishev et al. 2011). We treat the dust grains as test particles and consider them to be spherical, although there are strong indications that especially the larger dust grains are porous and fluffy aggregates (Fig. 1 in Schulz et al. 2015; Kolokolova & Kimura 2010). We can justify our choice of treating the dust as spherical particles with the fact that the main parameters influencing the dust trajectory are the mass and cross-section. Since we still know little about the particle shapes and, perhaps more importantly, about how they orient themselves in the flow, we interpret our parameters as representative of effective spheres. We have seen that the mfp of the dust particles in our models for reasonable setups is at least 2000 m, so we can safely neglect dust-dust collisions.

As a result, the equation of motion for each dust grain with mass m, radius r, and geometric cross-section σ = πr² at position x and with velocity ${{v}}_{\mathrm{d}}\mathrm{=}\frac{\mathrm{d}{x}}{\mathrm{d}\mathit{t}}$ is $\begin{array}{ccc}\mathit{m}\frac{{\mathrm{d}}^{\mathrm{2}}{x}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}& \mathrm{=}& {{F}}_{\mathrm{G}}\mathrm{+}{{F}}_{\mathrm{D}}\\ & \mathrm{=}& {{F}}_{\mathrm{G}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{C}}_{\mathrm{D}}{\mathit{m}}_{\mathrm{g}}{\mathit{n}}_{\mathrm{g}}\mathit{\sigma }\left|{{v}}_{\mathrm{g}}\mathrm{-}{{v}}_{\mathrm{d}}\right|\mathrm{(}{{v}}_{\mathrm{g}}\mathrm{-}{{v}}_{\mathrm{d}}\mathrm{)}\mathit{,}\end{array}$(5)where F_G is the gravitational force, m_g the mass of molecular water, and n_g and v_g are the number density and velocity of the gas. If we assume an equilibrium gas flow, the drag coefficient C_D is defined as ${\mathit{C}}_{\mathrm{D}}\mathrm{=}\frac{\mathrm{2}{\mathit{s}}^{\mathrm{2}}\mathrm{+}\mathrm{1}}{\sqrt{\mathit{\pi }}{\mathit{s}}^{\mathrm{3}}}{\mathrm{e}}^{\mathrm{-}{\mathit{s}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{4}{\mathit{s}}^{\mathrm{4}}\mathrm{+}\mathrm{4}{\mathit{s}}^{\mathrm{2}}\mathrm{-}\mathrm{1}}{\mathrm{2}{\mathit{s}}^{\mathrm{4}}}\mathrm{erf}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{+}\frac{\mathrm{2}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\epsilon }}^{\mathrm{)}}\sqrt{\mathit{\pi }}}{\mathrm{3}\mathit{s}}\sqrt{\frac{{\mathit{T}}_{\mathrm{d}}}{{\mathit{T}}_{\mathrm{g}}}}\mathit{,}$(6)with the gas temperature T_g, the dust temperature T_d (chosen to be equal to T_g), and ε is the fraction of specular reflection (set to 0), and $\mathit{s}\mathrm{=}\frac{\left|{{v}}_{\mathrm{g}}\mathrm{-}{{v}}_{\mathrm{d}}\right|}{\sqrt{\frac{\mathrm{2}\mathit{k}{\mathit{T}}_{\mathrm{g}}}{{\mathit{m}}_{\mathrm{g}}}}}\mathrm{·}$(7)The gravitational acceleration has been calculated according to the numerical approach described in Sect. 3.3.2.

Two further assumptions have been made. Firstly, as the drag force depends on the relative velocity, a non-Maxwellian velocity distribution may affect the result, an effect we have ignored. Skorov & Rickman (1999) have derived more complex expressions for the drag force taking the non-LTE (local thermodynamic equilibrium) in the Knudsen layer into account. They have found that the difference in grain speeds can be up to a factor of two. Secondly, concerning the liftability of dust grains, our model neglects surface cohesion. This may be justifed for certain ejection processes but can also be accounted for by modifying the ejected size distribution.

The equation of motion is solved by means of a fourth-order Runge-Kutta Method with an adaptive time step, ensuring that the spatial displacement is either less than one tenth of the characteristic size, l_c, of the cell (${\mathit{l}}_{\mathrm{c}}\mathrm{=}{}^{\mathrm{3}}\sqrt{{\mathit{V}}_{\mathrm{c}}}$, where V_c is the cell volume), or the change in velocity is less than 10% of the current velocity, whichever is lower.

We simulated 40 dust sizes in the range between 8 nm and 0.3 mm. The reason for this range will become apparent in Sect. 3.3.4. Approximately four million test particles were simulated in each run. The particles were placed randomly on the surface with initial velocity ${{v}}_{\mathrm{d}}^{\mathit{t}\mathrm{=}\mathrm{0}}\mathrm{=}\mathrm{0}$ and were subsequently tracked though the unstructured grid of the simulation domain (identical to the grid in the gas model) until they either reached the outlet surface at 10 km distance from the centre of the nucleus or “redeposited” on the inlet surface, i.e. the cometary surface.

The density of the dust particles was assumed to be 440 kg m^-3 for all dust sizes. Fulle et al. (2015) have estimated the dust density of very fluffy particles of sizes ranging between 0.2 to 2.5 mm detected with the Grain Impact Analyser and Dust Accumulator (GIADA) to <1 kg m^-3. We have considered smaller more compact particles in this study and thus chosen to assume a dust density that is close to the nucleus bulk density.

The simulation particle i is weighted with w_i to represent the physical number of particles leaving the surface facet per second, when presuming a dust-to-gas mass production rate ratio, Q_d/Q_g. To calculate, for example, the number density in a specific cell k, one simply has to measure the time ${\mathit{t}}_{\mathit{k}}^{\mathit{i}}$ that each simulation particle i spends in that cell k. The total number density n_k of cell k can thus be calculated by ${\mathit{n}}_{\mathit{k}}\mathrm{=}\frac{{\mathit{Q}}_{\mathrm{d}}}{{\mathit{Q}}_{\mathrm{g}}}\sum _{\mathit{i}}{\mathit{t}}_{\mathit{k}}^{\mathit{i}}{\mathit{w}}_{\mathit{i}}\mathit{.}$(8)Similarly the mean velocity of the dust particles in cell k can be written as ${{v̅}}_{\mathit{k}}\mathrm{=}\frac{{\sum }_{\mathit{i}}{\mathit{t}}_{\mathit{k}}^{\mathit{i}}{\mathit{w}}_{\mathit{i}}{{v}}_{\mathit{k}}^{\mathit{i}}}{{\sum }_{\mathit{i}}{\mathit{t}}_{\mathit{k}}^{\mathit{i}}{\mathit{w}}_{\mathit{i}}}$(9)and equally for any other dust parameter that one would want to be tracked.

3.3.2. Gravity field calculation

To determine the gravity field of a complex shape like comet 67P, we used a simple model. The gravitational acceleration of an arbitrary object exhibited at any point in space can be written as ${{a}}_{\mathrm{G}}\mathrm{\left(}{y}\mathrm{\right)}\mathrm{=}\mathit{G}{\mathrm{\int }}_{\mathit{V}}\mathit{\rho }\frac{{r}}{\mathrm{|}{r}{\mathrm{|}}^{\mathrm{3}}}\mathrm{d}\mathit{V}$(10)where r is the vector of point y to the volume element dV, G the gravitational constant, and ρ the local mass density of the body. Since this integral only has an analytical solution for the homogeneous sphere or point source, the value of a_G must be determined numerically in our case. We did this by discretizing the volume of the shape model “SHAP4S” of Preusker et al. (2015) with a resolution of 30 m resulting in 801′757 volume elements ΔV and assuming a constant density of 462 kg m^-3, which results in a total mass of 9.9 × 10¹² kg (Sierks et al. 2015). The integral thus reduces to a sum over all these elements: ${{a}}_{\mathrm{G}}\mathrm{\left(}{y}\mathrm{\right)}\mathrm{=}\mathit{G\rho }\mathrm{\Delta }\mathit{V}\sum _{\mathit{n}}\frac{{r}}{\mathrm{|}{r}{\mathrm{|}}^{\mathrm{3}}}\mathrm{·}$(11)This has been done for over 21 million points of a cubical grid up to a distance of 10 km from the centre of the nucleus. The resolution of the cubical grid is 30 m up to a distance of 3 km, 50 m between 3 and 5 km, and 100 m beyond 5 km distance to the nucleus centre.

3.3.3. Column densities from line-of-sight integration

The dust column density can be calculated from the number densities along a specific line of sight as an intermediate step towards an artificial image. At this point it is necessary to define a specific viewing geometry. The position of the spacecraft from the centre of the nucleus x_s/c, the Euler rotation angles (α,β,γ) of the virtual camera with respect to the spacecraft coordinates, and the pointing direction p of the spacecraft are all derived from calibration and the reconstructed spacecraft SPICE kernels. For calculating each pixel, we are using a planar projection onto the image plane compared to the actual angular field of view of the OSIRIS cameras. The artificial column density maps have a field of view of 20 km × 20 km in the image plane and are comprised of 1024 × 1024 pixels, resulting in a resolution of 19.5 m px^-1. The depth of the line-of-sight integration extends 200 km beyond the nucleus centre. Because of the quadratic footprint of our integration and the fact that we are integrating beyond our simulation domain of the gas-and-dust dynamics model, we need to extrapolate the dust number densities to up to a distance of 200 km from the comet. We do this by interpolating the 14′000 values at the outlet surface at 10 km to a 2D longitudinal and azimuthal coordinate grid with a angular resolution of 0.2°, yielding more than 1.6 million grid points. For the extrapolation to any point in 3D space, we assume that the number density falls off with the distance squared. The density at a point y = y(x,y,z) = y(r,θ,φ) is thus calculated by ${\mathit{n}}_{\mathrm{d}}\mathrm{\left(}{y}\mathrm{\right)}\mathrm{=}{\mathit{n}}_{\mathrm{d}}\mathrm{\left(}\mathit{r,\theta ,\phi }\mathrm{\right)}\mathrm{=}{\mathit{n}}_{\mathrm{d}}\mathrm{(}\mathrm{10}\hspace{0.17em}\mathrm{km}\mathit{,\theta ,\phi }\mathrm{)}{\left(\frac{\mathit{r}}{\mathrm{10}\hspace{0.17em}\mathrm{km}}\right)}^{\mathrm{2}}\mathrm{·}$(12)This method only works well for dust that has reached terminal velocity and expands primarily radially, so that it does not perform well in situations with dust jets that have significant non-radial components. However, this appears to be a good approximation for most of the cases we have considered so far. In any case, for small impact parameters, the column density is dominated by cells close to the nucleus.

This line-of-sight integration is executed for each size bin individually under the assumption of low optical depth.

3.3.4. Scattering model

For the scattering of the dust, we use Mie theory for spherical particles and the algorithm of Bohren & Huffman (1983). Mie theory and its application to comet observations have been extensively investigated by Fink & Rubin (2012), who provide plots for different values of the refractive index of the scattering efficiency versus the size parameter, x, defined as $\mathit{x}\mathrm{=}\frac{\mathrm{2}\mathit{\pi r}}{\mathit{\lambda }}\mathit{,}$(13)where λ is the wavelength and r the particle radius. The scattering efficiency typically shows a strong peak at values of x = 1−3. We compute the scattering functions for 40 discrete sizes, quasi-logarithmically distributed (10 linear bins per decade) between 8 nm to 0.3 mm. Larger particle sizes are computationally too time consuming for Mie theory, although their contribution to the total dust brightness from the inner coma is likely to be small (see below). Discretization can lead to oscillation in the scattering function in Mie calculations (see also Fink & Rubin 2012) and hence, we avoid this by computing a mean scattering function over a size range that is ±10% of the nominal discrete size. As refractive index we have chosen that of astronomical silicate (Laor & Draine 1993), which corresponds to n = 1.81 + 0.1012i for λ = 612 nm of the WAC dust continuum filter.

Under the assumption of zero optical depth, the observed radiance can be summed and compared to the expected radiance of a column of dust with a specified size distribution of unit mass. When referring to the size distribution of the dust, we always refer to the distribution of dust sizes of the injected flux produced at the surface prior to any acceleration. Two approaches have been adopted. Firstly, a power law distribution for the size is set at the surface using $\mathit{n}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{~}{\mathit{r}}^{\mathrm{-}\mathit{b}}\mathit{,}$(14)and we computed the resulting brightness “images” of the inner coma at the wavelength of the OSIRIS WAC 18 filter for direct comparison with the absolute calibrated OSIRIS data (Fornasier et al. 2015). Secondly, a single dust size distribution, n(r) ~ δ(r−r_i), with one large dust radius r_i was studied.

The brightness for each size bin is calculated according to $\mathit{I}\mathrm{=}\frac{\mathit{F}\mathrm{(}{\mathit{\lambda }}^{\mathrm{)}}}{{\mathit{R}}_{\mathrm{h}}^{\mathrm{2}}}{\mathit{n}}_{\mathrm{col}}{\mathit{\sigma }}_{\mathrm{geo}}{\mathit{q}}_{\mathrm{eff}}\frac{\mathit{p}\left(\mathit{\phi }\right)}{\mathrm{4}\mathit{\pi }}\mathit{,}$(15)where F is the solar flux at 1 AU, n_col the dust column density, σ_geo the geometric cross-section of the dust grain, q_eff the scattering efficiency, and p the phase function for the phase angle, φ. Individual images for each particle size can also be extracted to investigate the relative influence of each particle size on the result. It should be noted that we make no statement about how this initial power law is produced by the surface ejection mechanism. We recognize that many effects (e.g. charging or van der Waals forces) may lead to deviations from a simple power law distribution and, in particular, that there may be evidence that a population of larger particles is being ejected from the nucleus at the time of the observations used herein (Rotundi et al. 2015). Also potential fragmentation of grains will certainly change the size distribution.

4. Simulation results and comparison with data (COPS and OSIRIS)

4.1. Investigated geometries

The illumination geometries considered in this study have been selected for when the comet was at a heliocentric distance of 3.4 AU. At that time pre-perihelion, the Sun was at a northern latitude of 48°. Steady state gas solutions were run for longitudinal positions of the Sun at 0°, 50°, 90°, 140°, 180°, 230°, 270°, and 320°. The 0° longitude is aligned with the −x direction of the “SHAP4S” model and goes through the Imhotep region on the “body” lobe. We presumed that the models are valid within the time span between 16 August to 25 September 2014 because during this period, the heliocentric distance changed by only ± 4%, and the Sun direction was within a 4° cone of the used value.

Fig. 2 Cut through the 20 km simulation domain though the long axis of the comet of the insolation driven outgassing case with an absorbing surface where the Sun is at the solar longitude of 140°. The left panel shows the total gas speed and the right panel the radial component of the gas velocity relative to the total speed.

As in the case of the initial conditions on the cometary surface, we studied two specific issues. The first is whether the surface is reflecting or absorbing to any back flux striking it. Secondly we have investigated two different initial production rates. The first is solely insolation driven outgassing. As described in Sect. 3.1.2, we scale the inlet number density by a fixed factor in order to meet the observational constraint and to match the variation in the ROSINA measurements as the spacecraft moves and the comet rotates. For the absorbing insolation driven model,for example, the global scaling constant is of the order of 0.01 (i.e. 1% of the sublimation from a pure water ice surface at the given albedo) was applied to all facets to provide a number density at the spacecraft that agrees with the ROSINA/COPS data. One can envision this as the active fraction within a facet, where the rest of the facet is inert, so covered by an insulating layer of regolith. The second inlet condition we refer to as inhomogeneous outgassing. In contrast to the first case, the scaling constant is now regionally defined and the regional production rate scaling constants vary. In some sense the inhomogeneous outgassing is simply a regionally varying enhancement superimposed on the insolation driven outgassing. We present the results of our simulations of the absorbing surface in Sect. 4.2.1 and those of the diffusely reflecting surface in Sect. 4.2.2.

For the dust simulations, we assume a global dust-to-gas mass production rate ratio Q_d/Q_g applied to all surface facets. The viewing geometries are defined by the respective OSIRIS images for which we extracted the necessary geometrical parameters from the SPICE kernels (Acton 1996) provided by the European Space Operations Centre (ESOC).

4.2. Results of the gas simulations and comparison with COPS

All our PDSC⁺⁺ simulations were run successfully and resulted in the mean collisional separation (mcs) to mean free path (mfp) ratio being less than 1 at most points in the simulation domain and therefore acceptable according to recognized criteria (Bird 1994). Only for some cells in inhomogeneous outgassing cases did mfp/mfp exceed 1 but was well below 10, which we still consider to be reasonable.

We can use the 3D gas number density distribution within the first 10 km of the nucleus centre to extract the flow properties in our domain to compare with the in situ measurements performed on-board Rosetta by the ROSINA experiment, in particular the COPS nude gauge (NG) density sensor, when the spacecraft is within 10 km of the nucleus. However, this was not the case within the time frame we investigated in this study, so there is a need to extrapolate. When the spacecraft was at cometocentric distances greater than 10 km, we chose the corresponding value at 10 km and scale this, assuming the density falls off with the square of cometocentric distance. This is equivalent to saying that the flow is radial and that the gas speed does not increase beyond 10 km. We have found that the radial component of the gas velocity on the day side is larger than 99% of the total gas velocity at the 10 km boundary. Even on the night side, the radial component accounts for 95% of gas velocity. This is illustrated in the righthand panel of Fig. 2. We can thus consider the flow to be radial. The second condition was checked in several test cases with a 20 km simulation domain and shows only a minor increase in the gas speed beyond 10 km from the nucleus centre as seen in the lefthand panel of Fig. 2. This is in line with the findings of Tenishev et al. (2008; see Fig. 7 in Tenishev et al. 2008), and the subsequent comparison with the COPS data also showed that a constant radial outflow beyond 10 km is a fair assumption to make for a comet such as 67P.

4.2.1. Absorbing surface models

We first discuss our simulation results of the absorbing surface models. Any gas back flux to the cometary surface will be absorbed fully. The solid lines in Fig. 3 show the integrated production rates for a full revolution of the nucleus with the insolation driven (red) and inhomogeneous (blue) outgassing condition. For the insolation driven model, a global scaling constant (effective active fraction) of 0.014 (i.e. 1.4% of the sublimation from a pure water ice surface at the given albedo) has been applied to all surface facets. In contrast to the insolation driven case, the scaling constant for the inhomogeneous model is now regionally defined as shown in Fig. 4. The fluxes have been adjusted such that we have higher activity (0.1) in the Hapi and Hathor regions. This was identified as a location of a potentially strong source region by Lin et al. (2015) for the dust and Bieler et al. (2015) for the gas. Low emission (0.005) has been set in the Hatmehit, Maftet, Nut, and parts of Ma’at and Bastet regions on the “head” of the nucleus. Similarly very low values have been assigned to the Imhotep region, and parts of Khepry, Aten, and Ash on the “body”. Medium activity (0.015) has been assigned to the rest of the comet. Within these regions the outgassing is still driven by insolation. The injected gas production rates Q_g for these two models range from 1.32 to 2.19 kg s^-1 (3.94 × 10²⁶ to 6.56 × 10²⁶ molecules s^-1) for the insolation-driven and 2.39 to 3.59 kg s^-1 (7.15 × 10²⁶ to 1.07 × 10²⁵ molecules s^-1) for the inhomogeneous case over the entire revolution of the nucleus. The modulation in the production rate seen in Fig. 3 is a direct consequence of the shape of the nucleus and thus the variation in the cross section of the comet exposed to the Sun. For the inhomogeneous models, the differences in active areas modify the modulation further. Trials with other distributions have been tested, but we show here only the best result obtained so far. The criterion for setting the effective active fraction, hence the gas production rate, was to match the COPS/ROSINA data as well as possible while modulating the outgassing within the morphological units of Thomas et al. (2015a,b).

Figure 5 shows the comparison of the in situ measured number density and the two sub-models, the insolation-driven and inhomogeneous outgassing cases for an absorbing surface. The drawn line in our simulated data is solely for the purpose of readability and is not an interpolation of the model results. This must be kept in mind especially when we have no simulated data point at the centre of a peak, and thus the line will not show the actual maximum for when our simulation had been run at the exact geometry of the peak. The error bars drawn are ± 10% of the value and correspond to our estimate for the accuracy of our DSMC code and the applied extrapolation. We also truncated the time axis to show the main zones of interest that are, to the greatest extent, representative of the entire time interval.

Fig. 3 Global gas production rates Qg in kg s-1 for the insolation-driven (red) and inhomogeneous (blue) outgassing models with an absorbing (solid) and reflecting (dashed) surface and the eight solar geometries run that we modelled (vertical lines).

Fig. 4 Regional effective active fraction on the comet as seen from a north polar view used for the setting up the absorbing inhomogeneous outgassing model.

Fig. 5 Top panel: comparison of the COPS NG data with our insolation driven and inhomogeneous outgassing models with an absorbing nucleus surface over the period from the 21.08. to 23.09.2014. Middle up panel: cometocentric distance and the phase angle of the observations shown on the same scale showing how the spacecraft approached the comet towards Zone D but at relatively high phase. Middle lower panel: sub spacecraft latitude (left axis) and longitude (right axis) showing how in Zone D the spacecraft was moving towards the northern pole which is located in the Hapi region. Bottom panel: sub spacecraft local time (SCLT) and the local time at the position of the zero longitude meridian (CLT) which runs through Imhotep.

At this point, we highlight the main structure of the COPS NG data for this time period. First, the density is clearly increasing as the distance to the comet decreases. Second, the geometry plots show that the orbital velocity of Rosetta is low and that 67P rotates underneath the spacecraft. The COPS data have what can be interpreted as a diurnal structure that has two periods per nucleus rotation. One nucleus revolution of 12.4 h is indicated by the vertical lines in Fig. 5. The two maxima correspond to the instances when the spacecraft, if nadir pointing, would be observing the nucleus neck region. The more pronounced peak usually corresponds to the spacecraft observing the “neck” from the Anuket side of Hapi. The lower maximum corresponds to the spacecraft being on the Aker side of Hapi. The minimum after the higher peak corresponds to when the spacecraft is on the Imhotep side of the comet, and the second minimum corresponds to a spacecraft location above the Hatmehit side. Third, the magnitude of the number density increases with increasing spacecraft northern latitude as identified by Bieler et al. (2015).

When we compare the data with the models, we can see that for Zone A, both models fit fairly well. The insolation-driven model fits the data almost perfectly, while the inhomogeneous model is slightly inferior with the density systematically overestimated in this case. In Zone B we can see that the insolation-driven model underestimates the magnitude of the rotational oscillation with the minima very weak. This is an indication of low activity in the larger surroundings of Hatmehit. This is confirmed by our inhomogeneous model that is in fact reproducing this feature to higher accuracy.

Zone C shows three interesting trends in the insolation-driven model. First, this model is systematically underestimating the number density by around 30%. Second, while we should see two peaks per cometary day, we only see one. Third, the peaks seem to be shifted by approximately minus one hour. The inhomogeneous model shows improvement in the magnitude of the measurements but primarily in the shape of the curve that now exhibits the expected frequency over the nucleus rotation. There is no improvement in the shift already observed in the insolation-driven case. The shift in Zone C is most likely linked to the viewing geometry. The spacecraft is just in the process of crossing over the north pole, so this shift may be an indication that the structure of the more active polar region that we have placed in the Hapi region is not yet defined accurately enough. Thermal inertia or self-heating effects may also be relevant. Further investigation is needed to resolve this.

The insolation-driven model exhibits two trends in Zone D. As in Zone C the model underestimates the actual data, but especially in this case the minima are low. In the COPS data we see two peaks during one revolution of the nucleus with one being more pronounced than the other. The more pronounced peak corresponds to geometries where Hathor is illuminated and additionally in the field of view. In our insolation-driven model, what should be the higher peak is the lower one and vice versa. This inversion of peaks is an indicator of stronger activity from the Hathor region. The inhomogeneous model shows that this is indeed a viable explanation because the relative size of the maxima is correctly modelled in the inhomogeneous case.

In general, there is significant improvement in the fit to the data with the inhomogeneous model and especially the form of the curve matches the actual data considerably better. Neither the insolation-driven nor the inhomogeneous models fit the minima in Zone D. This presumably indicates night-side outgassing, which our thermal model cannot fit under the assumptions adopted. Bieler et al. (2015) use an ad hoc assumption of 7–10% of activity relative to the maximum flux coming from the night side and shadowed areas.

The downside of these models with an absorbing surface is that strictly speaking, we lose self consistency with Eq. (1)since we are effectively reducing the total flux from the inlet surface. Increased back flux is observed in geometries where either the Hathor or Seth regions are not illuminated, which shows that a part of the absorbed flux is due to the concave shape of the nucleus around the Hapi region and not to immediate reabsorption after sublimation on the original facet. When both Seth and Hathor are illuminated simultaneously, the outgassing prevents additional back flux from other regions and thus we see lower back flux rates in these illumination geometries. Because of the fairly high back flux in the inhomogeneous case, it becomes clear that the production rate must be higher than in the insolation-driven case to produce equivalent number densities at the spacecraft position, especially for the absorbing case. This motivated us to investigate whether reflecting surfaces can achieve similar results and at the same time fulfil self-consistency.

4.2.2. Reflecting surface models

As we have just seen in Sect. 4.2.1, high gas fluxes back to the cometary surface warrant the examination of reflecting surfaces. With a reflecting surface, none of the gas is absorbed when hitting the cometary surface so that compared to the absorbing surface models with identical setups, all subliming gas reaches the outlet surface. To reproduce the in situ measurements of the gas number density by COPS less gas thus needs to be produced, and consequently the effective active fraction is lower compared to the absorbing surface models. In the insolation-driven case, we used an effective active fraction of 0.012 compared to 0.014 with the absorbing surface. Unilluminated surfaces were set to a nominal temperature of 100 K. We can see in Fig. 6 the effect of the reflecting surface in a cut through the coma comparing the insolation-driven absorbing case with the diffusely reflecting case. The Sun is at 230° longitude so that Hathor is in shadow and Seth is illuminated. One can see that the coma is distorted in the absorbing case. The reflecting surface means we see higher densities at the surface and the coma being pushed back towards the “body” lobe.

Fig. 6 Cuts though the long axis of the comet with a solar longitude of 230° comparing the gas number density [m-3] of the insolation driven model with an absorbing surface (left column) and a diffuse reflecting surface (right column).

The comparison of the COPS data with the insolation-driven outgassing model with a diffusely reflecting surface shows that in some regions such as Zones A and C, the results are similar to the insolation-driven outgassing model with an absorbing surface. Figure 7 shows that the insolation-driven reflecting surface models exhibit similar problems to those with the absorbing surface (e.g. the missing peaks in Zone C). As in Sect. 4.2.1, we have adopted the approach of improving the shortcomings of the insolation-driven model by introducing regional heterogeneity. The fluxes were adjusted such that we have higher activity in the Hapi (0.075) and Hathor (0.04) regions. Low emission (0.005) were set in the Hatmehit, Maftet, Nut, and parts of Ma’at and Bastet regions on the “head” of the nucleus. Similarly very low values were assigned to the Imhotep region and parts of Khepry, Aten, and Ash on the “body”. Medium activity (0.0095) was assigned to the rest of the comet. Within these regions the outgassing is still insolation-driven though the outgassing power varies regionally. Compared to the absorbing surface, the Hathor region’s activity is weaker than the one in Hapi. The gas production rates Q_g for these two models range from 1.28 to 2.20 kg s^-1 (3.82 × 10²⁶ to 6.59 × 10²⁶ molecules s^-1) for the insolation-driven and 1.54 to 2.58 kg s^-1 (4.61 × 10²⁶ to 7.71 × 10²⁶ molecules s^-1) for the inhomogeneous case over the entire revolution of the nucleus. The global production rates for the entire nucleus rotation can be seen in Fig. 3.

As with the absorbing surface, this inhomogeneous model improves the fit to the ROSINA/COPS data in all regions.

4.2.3. Summary of the gas models

We ran insolation-driven and regionally inhomogeneous outgassing models with absorbing and diffusely reflecting nucleus surface boundary conditions and compared the simulations with the in situ gas number density measurements of the COPS NG. For both surface behaviours, we have seen that insolation driven outgassing does not reproduce the data to a satisfactory degree. The shape of the ROSINA/COPS curves in particular indicate the inadequacy of the model. Regionally inhomogeneous outgassing models have, both in the absorbing and the reflecting cases, improved the fit to the data and overcome most shortcomings of the insolation-driven outgassing.

To place this on a more objective footing, we use the approach of Bieler et al. (2015) and use the Pearson product-moment correlation coefficient (PPMCC) which shows that the inhomogeneous models (PPMCC = 0.87) are statistically better fits to the data than the insolation-driven models (PPMCC = 0.816/0.834 for the absorbing/reflecting surface) and compare well to the values in Bieler et al. (2015). We must stress, though, that compared to Bieler et al. (2015), we do not need any post simulation correction of the simulation results to achieve these fits. Additionally,we do not consider any artificially introduced night-side activity, which leads to underestimating the number density over lower latitudes where CO₂ emission becomes relevant. The physical process of the CO₂ production, especially in the southern hemisphere, has not been tackled and understood yet so certainly warrants further study.

That we can achieve a good correlation with both the absorbing and reflecting surfaces shows that ROSINA/COPS alone cannot distinguish between these two types of surface properties. On the other hand, line-of-sight instruments, such as MIRO or VIRTIS, should be able to see such a difference in surface reflectance, as can be suspected from Fig. 6. We look at this in more detail in Sect. 4.2.4.

One further point in need of attention is that we have only looked at heterogeneities on a regional scale. The results shown can be seen as a hint of more local inhomogeneities than assumed in our models. The measured gas densities in viewing geometries that are at high northern latitudes, especially directly over the pole, seem to indicate this. Thus locally inhomogeneous models may need to be studied.

Fig. 7 Top panel: comparison of the COPS NG data with our insolation-driven and inhomogeneous outgassing models with an absorbing and a diffusely reflecting nucleus surface over the period from the 29.08. to 23.09.2014. Upper middle panel: cometocentric distance and the phase angle of the observations shown on the same scale showing how the spacecraft approached the comet towards Zone D but at relatively high phase. Middle lower panel: sub-spacecraft latitude (left axis) and longitude (right axis) showing how in Zone D the spacecraft was moving towards the northern pole, which is located in the Hapi region. Bottom panel: sub-spacecraft local time (SCLT) and the local time at the position of the zero longitude meridian (CLT) that runs through Imhotep.

4.2.4. Gas column densities

Having the 3D gas number density distribution within the first 10 km allows us to produce a water column density within the field of view using a line-of-sight integration. We did this in a similar way to the extrapolation described in Sect. 3.3.3. The results shown in Fig. 8 for 7 September 2014 at 12:30:00 UTC (identical viewing geometry as in Fig. 6 in Biver et al. 2015). We clearly see the influence in the case of inhomogeneous outgassing boosting the column density by at least a factor of 2. The influence of the reflecting surface can be seen especially in the inhomogeneous case where the coma is pushed away from the “head” lobe towards the “body” lobe (top of the image), and thus the main direction of the gas flow is clearly deflected by approximately 15°. The inhomogeneous outgassing models with increased activity from the “neck” structurally agree with the conclusion of Biver et al. (2015). However, our column densities are roughly one order of magnitude higher even though our production rates are in general agreement with their derived production rate of 1.46 ± 0.75 kg s^-1. (Our values are within that interval for the insolation-driven outgassing and slightly higher for the inhomogeneous outgassing models as seen in Fig. 3.) The observations of MIRO were performed during a time period of 4 h (or 116° of comet rotation) compared to our steady state solution in Fig. 8. Therefore a more detailed analysis will have to be performed when the datasets can be combined to resolve this discrepancy. VIRTIS observations could provide an additional constraint. A first comparison has shown that our column densities resemble the ones reported by VIRTIS in Bockelée-Morvan et al. (2015).

Fig. 8 Comparison view of the water column density in log 10(cm-2) for the insolation-driven (left) and inhomogeneous (right) outgassing cases with an absorbing (top) and a reflecting (bottom) surface using the viewing geometry on 7 September 2014 at 12:30:00 UTC. The Sun is at a longitude of 180°.

4.3. Results of the dust simulations and comparison with OSIRIS

Although we have made considerable effort to fit the COPS data, our main objective was to take this a step further and fit the OSIRIS dust brightness measurements. The dust results shown here have been created in an attempt to reproduce the OSIRIS image WAC_2014-09-05T09.19.13.810Z_ID30_1397549700_F18 taken with the wide angle camera (WAC) on 5 September 2014 at 9:20:23 UTC with filter 18 from a distance of 42.5 km to the nucleus centre, a phase angle of 59°, and an exposure time of 0.469 s. In this image, seen in Fig. 9, the Sun is at 140° longitude located in the top left of the image. The image on the right has been stretched to show the dust emission. The data have been absolute-calibrated (Fornasier et al. 2015) and are in physical units of spectral radiance. This processed form of the image clearly shows a fan-like structure of the dust brightness. Other than this feature, we do not see any other major dust activity.

The coma brightness depends on the dust size distribution. We have thus included two types of distribution in our analysis: (1) a power law distribution over a large interval of dust sizes from tens of nanometres to tenths of millimetres favouring a coma dominated by small particles and (2) a single-size coma dominated by large particles of tens or hundreds of micrometres.

Fig. 9 Wide angle camera OSIRIS image (WAC_2014-09-05T09.19.13.810Z_ID30_1397549700_F18) taken with filter no. 18 on the 5. September 2014 at 9:20:23 UTC from a distance of 42.5 km distance to the nucleus centre, a phase angle of 59°, and an exposure time of 0.469 s. This image shows the spectral radiance in [W m-2 nm-1 sr-1] with the left one being on a linear scale with a maximal radiance of 1 × 10-3 W m-2 nm-1 sr-1 and the right one being stretched to the range 3 × 10-7 to 3 × 10-4 W m-2 nm-1 sr-1.

4.3.1. Coma with a power law size distribution

With a power law for the dust size distribution where the number of particles of a given size n(r) is proportional to r^{− b} (Eq. (14)) at the inlet boundary, we ran cases for different values of the power law exponent, b, and for the insolation-driven and inhomogeneous outgassing solutions for the two surface types presented in Sect. 4.2. A power law of b = 3 amounts to an equal mass distribution for any dust size and b = 4 an equal mass distribution over each dust size decade. Figure 10 shows the spectral brightness with insolation-driven outgassing and the reflecting surface of the individual dust size bin run with b = 4.5. The final image is simply the sum over these 40 partial brightnesses.

Fig. 10 Spectral radiance in [W m-2 nm-1 sr-1] for all the 40 dust size bins between 8 × 10-8 and 3 × 10-4 m for the power law exponent b = 4.5 of the insolation driven outgassing case with a reflecting surface. The individual size ranges are indicated in each plate.

The main contribution to the total coma brightness for such a coma with b = 4.5 comes from the size bins in tenths of μm range, and the small and large dust sizes only provide a small contribution. This is a combined effect of the power law and the scattering efficiency. Even though large particles are moving slowly, and thus in principle producing higher densities, this is more than compensated for by the combined effect of the power law and low scattering efficiency of large particles both of which drive the radiance down. For very small particles, the situation is slightly different. Owing to the power law, they are much more numerous, which might suggest a high brightness. But because they are much smaller than the observation wavelength, their scattering efficiency drops to almost zero, which results in a small contribution to the total radiance from these sizes. This is also illustrated in Fig. 11, which shows the phase function, p, multiplied by the scattering efficiency, q, of spherical particles for different dust radii as a function of scattering angle in a polar diagram. Small particles are very much fainter than all other sizes because their scattering efficiency is low, and they scatter uniformly in all directions.

Because the power law exponent increases from low values to 4, we have an increasing number of particles in the efficient scattering regime for the same dust mass, so the total spectral radiance increases. When the power law exponent goes beyond 4, most of the particles are found in the very small size bins, which have an almost negligible scattering efficiency, resulting in a drop in the total brightness. We can put this differently by saying that for a constant (measured) brightness of the coma, the dust-to-gas mass production rate ratio Q_d/Q_g must decrease because the power law exponent increases in the range of low values up to 4 (i.e. there is degeneracy between b and Q_d/Q_g). This is illustrated in Fig. 12 where we can see Q_d/Q_g as a function of the power law exponent, and the figure shows a coma of constant brightness (in this case the one of the OSIRIS image WAC_2014-09-05T09.19.13.810Z_ID30_1397549700_F18 we have chosen). This figure furthermore shows the degeneracy we have just mentioned. For a given Q_d/Q_g, up to two power law exponents will produce the desired brightness. Recent results (Rotundi et al. 2015) suggest that b ~ 3 for large particles (>mm) and closer to 2 for sub-millimetre particles and that Q_d/Q_g = 4 ± 2, which would, however, constrain this result and lift the degeneracy (light red area in Fig. 12). Our models are able to reproduce results within these constraints, though we are at the lower end of the constraining interval. It must be stressed, however, that for b< 3.5 the brightness is dominated by the large size bin, hence the upper limit of our considered size interval, such that the size of the interval can become important by influencing the trend of the curve in Fig. 12. We have therefore also studied single-size comae in Sect. 4.3.2 to circumvent this problem.

We should highlight at this point that the power law exponent is set at the surface. As the different dust sizes accelerate at different rates, i.e. the small particles undergo higher acceleration compared to the smaller ones, the small particles are diluted compared to the large particles, thereby changing the power law. The effective or mean power-law exponent of the total coma will thus always be lower than the one at the surface (McDonnell et al. 1987) because the dust speeds as a function of the dust size, r, typically scale with r^-0.5.

Fig. 11 Polar diagram of the phase function, p, multiplied by the scattering efficiency, q, of spherical particles for different dust size radii as a function of scattering angle. Large particles exhibit a strong forward scattering peak, while smaller particles scatter more uniformly.

Fig. 12 Dust-to-gas mass production rate ratio as a function of the power law exponent, b. The solid blue line represents the mean value of the four models run (insolation-driven/inhomogeneous outgassing and reflecting/absorbing surface), and the blue band indicates the maximum and minimum values obtained. The red area shows the constrained area by Rotundi et al. (2015).

The outgassing distribution also affects the dust distribution. Figure 13 shows a comparison of the total brightness of the dust coma for the insolation-driven and inhomogeneous outgassing solutions. By comparing with the actual image, we can see immediately that there is considerably more dust visible in our models on the entire illuminated part of the comet. The dust filaments seen are a product of the nucleus shape, including the varying illumination conditions across the surface and not of any manually introduced jet sources.

None of these models reproduces the data well. This indicates that we either (1) are dealing with very locally varying outgassing rather than the broad regional inhomogeneities we were considering in Sect. 4.2 (this could also include small scale shadowing effects below the currently used resolution); (2) have regionally or locally varying Q_d/Q_g; or (3) have regionally or locally varying power law exponents, b. Comparing the insolation-driven and inhomogeneous solutions, we note that the structure of the coma varies substantially.Since the dust is tracing the motion of the gas to some extent, we immediately see the differences caused by the inhomogeneous outgassing and thus higher activities in the Hapi and Hathor regions. The coma filaments have clearly shifted to the Seth side of the comet from the insolation-driven to the inhomogeneous outgassing case. The inhomogeneous models come closer to the data because they show a more pronounced main filament originating in the Hapi region and more faint emission everywhere else. Additionally, the influence of the reflecting surface can be seen especially with the main filament turning towards the “body” lobe.

Of the four models shown here, the inhomogeneous model with a reflecting surface resembles the actual data the closest, particularly since the main filament coming from the Hapi region has the correct direction. We must stress that dust filaments coming from the Hatmehit region on the “head” lobe and the Seth, Ash, Aten, and Babi regions on the “body” lobe can still be seen even though they are not observed in our specific OSIRIS image. We can also see bent jet filaments near the surface in the “neck” part of the comet, especially in the insolation-driven outgassing models. Because the dust speeds are too high (see Table 1), this bending is not caused by gravity but rather by the gas drag.

Fig. 13 Spectral brightness of a power law dust coma for insolation-driven (left column) and inhomogeneous (right column) outgassing for the absorbing (top row), and reflecting (bottom row) surface with a power law exponent of b = 4.5 and dust to gas mass production rate ratios Qd/Qg = 0.075.

4.3.2. Coma dominated by large particles

Finally we examined a coma dominated by large particles in contrast to a power law dust size distribution. This was motivated by the inferred dominance of large particles from the GIADA and COSIMA experiments (Rotundi et al. (2015); Langevin presentation at SWT 2015). We looked especially at a coma consisting solely of particles with radii of 3.18 × 10^-4 m. Figure 14 shows the results for a reflecting surface with insolation-driven and inhomogeneous outgassing.

Fig. 14 Spectral brightness of the dust for insolation-driven (left) and inhomogeneous (right) outgassing for a coma consisting only of 318 μm sized particles.

What immediately stands out is that the insolation-driven model in no way resembles the real coma. The inhomogeneous outgassing model is a much better approximation and structurally resembles the actual coma quite well. In particular, the inhomogeneous coma dominated by 0.318 mm particles no longer shows the filament from Hatmehit and only retains a few excess jets coming from the Seth, Ash, Aten, and Babi regions. This case structurally resembles the observed coma most closely. The dust-to-gas mass production rate ratio in this case is Q_d/Q_g = 2.2, so it is within the expected range. A coma dominated by large particles thus seems plausible but is unlikely to be a unique solution.

The travel time needed to reach the edge of the domain for the large particles is not negligible (90 min worst case for 318 μm radius particles in the insolation driven flow field), so that ultimately a time-dependent solution will be required to account for the rotation of the comet. Figure 15 shows the amount of dust mass reaching the outlet surface for different dust sizes in the range of 1 μm to 1 cm. Approximately 27% of the particles do not reach the outlet surface regardless of the size.

These are the particles that mainly get reabsorbed in the concave shape of the comet’s “neck” region. Additionally, some of these particles are on ballistic trajectories because they originate near the terminator. Since the dust size increases beyond 0.1 mm, the number of particles reaching the outlet surface decreases further as areas with lower activity can no longer lift the dust particles. For the insolation-driven outgassing model with an absorbing surface, the maximum liftable mass of our particles with a density of 440 kg m^-3 lies at 1.5 mm. For the inhomogeneous case, the largest liftable size is 5 mm. It is higher because we have stronger gas emission from the “neck” region compared to the insolation-driven case. It needs to be stressed that localized activity within the facet might lead to larger particles being lifted.

Fig. 15 Mass fraction of dust reaching the outlet boundary as a function of the dust size in the size range of 1 μm to 1 cm for the insolation driven (red triangles) and inhomogeneous (blue circles) outgassing model with an absorbing surface. The dust density is 440 kg m-3 for all sizes.

For a point source the brightness drops off as $\frac{\mathrm{1}}{\mathit{r}}$, where r is the impact parameter of the line of sight with respect to the nucleus in the image plane. In all of our simulations we have observed deviations from this behaviour when studying profiles along jets and total radial brightness behaviour, especially close to the nucleus. As a result, the observation of non $\frac{\mathrm{1}}{\mathit{r}}$ behaviour (Lin et al. 2015) can be explained by a combination of non-point source geometry and acceleration (e.g. Thomas & Keller 1990) although a detailed comparison should still be performed. The signal-to-noise in the August-September 2014 observation is, however, low in the data and not optimal. Extension of this work to times of higher activity (near perihelion) would offer better prospects for evaluation and comparison. The deviation from a $\frac{\mathrm{1}}{\mathit{r}}$ drop of the brightness in our simulations is illustrated in Fig. 16. The lefthand panel shows a polar transform of the dust brightness of the smallest size bin of the insolation-driven outgassing model with an absorbing surface. The righthand panel shows the integrated brightness multiplied by the radial distance to the nucleus centre, I(r)·r, as a function of the radial distance, r, given by $\mathit{I}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{r}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathit{I}\mathrm{\left(}\mathit{r,\phi }\mathrm{\right)}\mathit{r}\mathrm{d}\mathit{\phi ,}$(16)with I(r,φ) being the spectral radiance at the respective point in the image plane given in polar coordinates, (r,φ).

The two curves shown correspond to the smallest (8 nm, green line) and largest bins (0.32 mm, dashed red line) we ran. The values have been scaled relative to each other at a radial distance of 4000 m from the comet centre. A freely outflowing coma from a point source would produce a drop in number density of $\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}$ and thus a drop in brightness of $\frac{\mathrm{1}}{\mathit{r}}$. Were this to be multiplied by the radial distance, r, we would thus expect a constant line in the righthand panel of Fig. 16. But we observe a drop in brightness that is steeper than $\frac{\mathrm{1}}{\mathit{r}}$. This deviation indicates that our flow is not radial and the dust is still accelerating. For radial distances larger than four kilometres the brightness drops comparably for the large and small particles and flattens out towards the edge of our simulation domain. This indicates that the coma indeed goes over to a radial flow with constant dust speed when our particles start reaching their terminal velocity. More interesting is the regime between 2.5 and 4.0 km. Here we would expect the small particles to exhibit a higher relative brightness compared to the large particles owing to their larger relative acceleration. But the exact opposite is the case. Our analysis has shown that this is caused by the fact that the small particles go over to a radial outflow faster than the large particles. The large particles have a more significant non-negligible tangential motion close to the surface and thus provide higher relative brightnesses. This non-radial motion is primarily driven by the non-spherical shape of the nucleus itself.

Fig. 16 Left panel: polar transform of the spectral brightness of the smallest dust size bin of the insolation-driven outgassing model with an absorbing surface. Right panel: integrated brightness at a given radial distance multiplied by the radial distance, I(r)r, as a function of the radial distance, r, calculated according to Eq. (16). The two curves shown correspond to the smallest (8 nm, green line) and largest size bins (0.32 mm, dashed red line) that we ran. The values have been scaled relative to each other at a radial distance of 4000 m from the comet’s centre

5. Discussion

Inevitably, given the number of free parameters implied by the complete scheme shown in Fig. 1, combined with the resolution of the simulation, many assumptions are needed to obtain a tractable and useful result at this stage of the analysis. We have deliberately used the simplest thermal model for setting the initial boundary condition in this first work. This is in part because the model used has essentially no free parameters. (Only the hemispherical albedo and the emissivity can be chosen, the albedo having very little influence within the expected range.) The low thermal inertia shown by other experiments would indicate that the approximation is suitable. However, the assumption of surface sublimation provides considerably more uncertainty. Conductivity provides heat for the sub-surface, which would lead to sublimation and gas transport through a porous surface layer. We have assumed initial local thermal equilibrium (LTE) at the source. This might be envisaged as gas equilibrating while flowing through the porous surface layer. However, if there is direct ejection, then there will be less energy for the expansion because the rotational degrees of freedom will not have had time to equilibrate. On the other hand, the gas could be warmer than the free sublimation temperature of water ice if there is sub-surface sublimation followed by interaction with the hot inert surface layer. These processes affect the total energy of the gas and thus influences the final gas speed and the lateral expansion. Their effect on the solutions presented need to be assessed.

The outgassing distribution on the surface has been assumed, even in the inhomogeneous case, to be insolation-driven on a regional scale. The morphological evidence from OSIRIS (Thomas et al. 2015a,b) suggests that deposition of larger particles back onto the surface has produced a coating of “airfall” in several regions. There is no strong evidence for dust and/or gas emission from these coated areas. As a result, locally, the outgassing may also be inhomogeneous. This will affect the properties of the gas flow in several ways. For example, local activity will result in lateral expansion close to the source, reducing the energy available for radial expansion and lowering the final gas speed. The interaction between the molecules will have raised their temperature relative to a single source. It remains to be assessed whether the resolution of the grid and the description of active areas can be adapted to provide a realistic simulation of the varying surface outgassing properties.

At this stage we assumed a single species. We are aware that CO₂ has a high mixing ratio above the southern hemisphere (Hässig et al. 2015). We have studied absorbing and reflecting surfaces for the gas dynamics and seen an influence of this choice that can be distinguished not by the ROSINA/COPS NG but in the structure of the dust coma. We also assumed two extremes for the dust with a specific refractive index (that of astronomical silicate, Laor & Draine 1993), a power law distribution and a single fixed size. Although these are common assumptions, they may be violated for numerous reasons in order to match the integrated brightness seen by OSIRIS. Perhaps more critical is that we have assumed as an input condition that the local dust-to-gas production function at the surface is a global constant only modified naturally by means of dust not being liftable depending on the illumination. Given the remarkably diverse surface morphology, this appears to be unjustifiable, but we have no clear evidence to the contrary.

Finally, to keep the computation time down, we have used a small domain and extrapolated where necessary. Ultimately, once the computed result in the 10 km domain appears close to final, a longer run with a larger domain has to be performed to verify our extrapolations. One example of that has, however, been presented here.

6. Conclusions

We have studied a wide range of initial conditions to examine the neutral gas coma for a heliocentric distance of 3.4 AU pre-perihelion, ranging from different surface properties (absorbing to reflecting) to outgassing distributions (insolation-driven to regionally inhomogeneous). This analysis has shown that in the considered time span, the most straightforward assumption of purely insolation-driven outgassing has no solutions fitting the in situ COPS/ROSINA measurements regardless of which surface property is chosen. Furthermore, based on these insolation-driven outgassing models,dust comae do not reproduce the observed dust coma structure either. In this sense the data sets of COPS/ROSINA and OSIRIS are consistent with insolation-driven outgassing not being a viable solution.

On the other hand, we have found fits using an absorbing and reflecting surface with increased activity from the Hapi and Hathor regions that come remarkably close to the observations of COPS/ROSINA. This kind of regionally inhomogeneous outgassing model is also consistent with measurements presented by the MIRO (Biver et al. 2015) and VIRTIS (Bockelée-Morvan et al. 2015) teams. Furthermore, our conclusion of increased activity in the “neck” region is consistent with the findings of Bieler et al. (2015), who have found an improvement in their fits to the COPS NG date with an a posteriori skew of the activity in higher northern latitudes. The shortcomings of our fits point to stronger locally varying outgassing especially in the northern polar region. This and the fact that multiple inputs can lead to equally good fits to the data show the importance of including data from other instruments, such as MIRO and VIRTIS, and other data from ROSINA to further constrain some of the free parameters.

Regarding the dust coma we have learnt that for a dust-to-gas mass production rate ratio, Q_d/Q_g, in the range of 1–10, we require a power law exponent, b, in the range of 2–3 to match the observed dust coma brightness. This is consistent with the findings of Rotundi et al. (2015). However these power law models exhibit too many dust filaments and thus do not match the observed coma structure particularly well. This implies either a locally varying Q_d/Q_g or a locally varying dust emission distribution. On the other hand, we have seen that single-sized comae with large particles (r_d ~ 100 μm) and dust densities of 440 kg m^-3 fit the OSIRIS data rather well. As a by-product, we showed how the assumed surface boundary condition (i.e. absorbing vs. reflecting) can change the direction of the dust jet originating in the Hapi region. Additional data from dust instruments such as GIADA and COSIMA should be used to check our models and help constrain some of the dust parameters.

Acknowledgments

OSIRIS was built by a consortium of the Max-Planck-Institut für Sonnensystemforschung in Göttingen, Germany; CISAS-University of Padova, Italy; the Laboratoire d’Astrophysique de Marseille, France; the Instituto de Astrofisica de Andalucia, CSIC, Granada, Spain, the Research and Scientific Support Department of the European Space Agency, Noordwijk, The Netherlands; the Instituto Nacional de Tecnica Aeroespacial, Madrid, Spain; the Universidad Politechnica de Madrid, Spain the Department of Physics and Astronomy of Uppsala University, Sweden; and the Institut für Datentechnik und Kommunikationsnetze der Technischen Universität Braunschweig, Germany. The support of the national funding agencies of Germany (DLR), France (CNES), Italy (ASI), Spain (MEC), Sweden (SNSB), and the ESA Technical Directorate is gratefully acknowledged. Work on ROSINA at the University of Bern was funded by the State of Bern, the Swiss National Science Foundation, and the ESA PRODEX Programme. The team from the University of Bern is supported through the Swiss National Science Foundation and through NCCR PlanetS. We thank Marco Fulle for pointing out a verification method that led to identiying an error.

References

All Tables

All Figures

	Fig. 1 Diagram of our modelling and verification approach.
In the text

	Fig. 2 Cut through the 20 km simulation domain though the long axis of the comet of the insolation driven outgassing case with an absorbing surface where the Sun is at the solar longitude of 140°. The left panel shows the total gas speed and the right panel the radial component of the gas velocity relative to the total speed.
In the text

	Fig. 3 Global gas production rates Qg in kg s-1 for the insolation-driven (red) and inhomogeneous (blue) outgassing models with an absorbing (solid) and reflecting (dashed) surface and the eight solar geometries run that we modelled (vertical lines).
In the text

	Fig. 4 Regional effective active fraction on the comet as seen from a north polar view used for the setting up the absorbing inhomogeneous outgassing model.
In the text

Fig. 5 Top panel: comparison of the COPS NG data with our insolation driven and inhomogeneous outgassing models with an absorbing nucleus surface over the period from the 21.08. to 23.09.2014. Middle up panel: cometocentric distance and the phase angle of the observations shown on the same scale showing how the spacecraft approached the comet towards Zone D but at relatively high phase. Middle lower panel: sub spacecraft latitude (left axis) and longitude (right axis) showing how in Zone D the spacecraft was moving towards the northern pole which is located in the Hapi region. Bottom panel: sub spacecraft local time (SCLT) and the local time at the position of the zero longitude meridian (CLT) which runs through Imhotep.

In the text

	Fig. 6 Cuts though the long axis of the comet with a solar longitude of 230° comparing the gas number density [m-3] of the insolation driven model with an absorbing surface (left column) and a diffuse reflecting surface (right column).
In the text

Fig. 7 Top panel: comparison of the COPS NG data with our insolation-driven and inhomogeneous outgassing models with an absorbing and a diffusely reflecting nucleus surface over the period from the 29.08. to 23.09.2014. Upper middle panel: cometocentric distance and the phase angle of the observations shown on the same scale showing how the spacecraft approached the comet towards Zone D but at relatively high phase. Middle lower panel: sub-spacecraft latitude (left axis) and longitude (right axis) showing how in Zone D the spacecraft was moving towards the northern pole, which is located in the Hapi region. Bottom panel: sub-spacecraft local time (SCLT) and the local time at the position of the zero longitude meridian (CLT) that runs through Imhotep.

In the text

	Fig. 8 Comparison view of the water column density in log 10(cm-2) for the insolation-driven (left) and inhomogeneous (right) outgassing cases with an absorbing (top) and a reflecting (bottom) surface using the viewing geometry on 7 September 2014 at 12:30:00 UTC. The Sun is at a longitude of 180°.
In the text

Fig. 9 Wide angle camera OSIRIS image (WAC_2014-09-05T09.19.13.810Z_ID30_1397549700_F18) taken with filter no. 18 on the 5. September 2014 at 9:20:23 UTC from a distance of 42.5 km distance to the nucleus centre, a phase angle of 59°, and an exposure time of 0.469 s. This image shows the spectral radiance in [W m-2 nm-1 sr-1] with the left one being on a linear scale with a maximal radiance of 1 × 10-3 W m-2 nm-1 sr-1 and the right one being stretched to the range 3 × 10-7 to 3 × 10-4 W m-2 nm-1 sr-1.

In the text

	Fig. 10 Spectral radiance in [W m-2 nm-1 sr-1] for all the 40 dust size bins between 8 × 10-8 and 3 × 10-4 m for the power law exponent b = 4.5 of the insolation driven outgassing case with a reflecting surface. The individual size ranges are indicated in each plate.
In the text

	Fig. 11 Polar diagram of the phase function, p, multiplied by the scattering efficiency, q, of spherical particles for different dust size radii as a function of scattering angle. Large particles exhibit a strong forward scattering peak, while smaller particles scatter more uniformly.
In the text

	Fig. 12 Dust-to-gas mass production rate ratio as a function of the power law exponent, b. The solid blue line represents the mean value of the four models run (insolation-driven/inhomogeneous outgassing and reflecting/absorbing surface), and the blue band indicates the maximum and minimum values obtained. The red area shows the constrained area by Rotundi et al. (2015).
In the text

	Fig. 13 Spectral brightness of a power law dust coma for insolation-driven (left column) and inhomogeneous (right column) outgassing for the absorbing (top row), and reflecting (bottom row) surface with a power law exponent of b = 4.5 and dust to gas mass production rate ratios Qd/Qg = 0.075.
In the text

	Fig. 14 Spectral brightness of the dust for insolation-driven (left) and inhomogeneous (right) outgassing for a coma consisting only of 318 μm sized particles.
In the text

	Fig. 15 Mass fraction of dust reaching the outlet boundary as a function of the dust size in the size range of 1 μm to 1 cm for the insolation driven (red triangles) and inhomogeneous (blue circles) outgassing model with an absorbing surface. The dust density is 440 kg m-3 for all sizes.
In the text

Fig. 16 Left panel: polar transform of the spectral brightness of the smallest dust size bin of the insolation-driven outgassing model with an absorbing surface. Right panel: integrated brightness at a given radial distance multiplied by the radial distance, I(r)r, as a function of the radial distance, r, calculated according to Eq. (16). The two curves shown correspond to the smallest (8 nm, green line) and largest size bins (0.32 mm, dashed red line) that we ran. The values have been scaled relative to each other at a radial distance of 4000 m from the comet’s centre

In the text