© R. Marschall et al. 2022

1 Introduction

To date, only seven comets (1P/Halley, 19P/Borrelly, 9P/Tempel 1, 26P/Grigg-Skjellerup [no nucleus imaging], 67P/Churyumov-Gerasimeko, 81P/Wild 2, and 103P/Hartley 2) have been visited by spacecraft. Spacecraft observations provide a detailed, high spatial, and temporal resolution of the surface and surrounding coma, but they are limited to a few target comets (~1–2 per decade).

While these observations are invaluable to understanding the building blocks of our Solar System and how their activity works, exploring them in close proximity bares a potential risk for spacecraft. This is particularly true when a spacecraft passes in close proximity (a few hundred kilometres or closer) of an active comet nucleus during a high velocity flyby.

When comets enter the inner Solar System and heat up sufficiently, their ices begin to sublimate. Additionally, the emitted gas drags dust particles into space and they form the so-called gas and dust comae. Depending on the spatial and temporal scales, the dynamics of the dust particles is primarily governed by the interaction with the gas, gravity from the nucleus and the Sun, and solar radiation pressure. These forces shape the distribution of dust in the vicinity of the comet nucleus as well as far into the tail of a comet millions of kilometres from the comet’s surface.

During a flyby through the coma, dust particles can collide with the spacecraft. Depending on the impact energy, such impacts can pollute the surface of instruments and/or disturb attitude control. For example, four particles between 0.5 and 10 mg impacted the Deep Impact Impactor before it impacted Tempel 1 (A’Hearn et al. 2008). Nine particles (~0.01–~0.1 mg) impacted on the Deep Impact spacecraft during its flyby of Hartley 2 (Hermalyn et al. 2013). Further, several discrete events caused by dust particles in the 1–50 mg range were reported by the Halley Multicolor Camera (HMC) on Giotto (Curdt & Keller 1990). Severe damage to the spacecraft can occur at high impact energies that is to say when the dust mass and/or impact speed are high. An example of the latter is the loss of the HMC camera on Giotto due to a large dust particle impact. It is crucial to quantify the dust environment to (i) guide trajectory design, (ii) inform spacecraft shielding, (iii) assess the performance of star trackers, (iv) evaluate the performance of scientific instruments (e.g. dust sensors and accumulators), and (v) determine the stability of attitude control.

ESA’s Comet Interceptor (CI) mission Snodgrass & Jones (2019) will pass through a potentially hazardous region of the inner coma of a dynamically new comet or long period comet. This is an indispensable condition to get better resolution from remote observations and to collect in situ data on dust and gas environment. It is therefore important to assess the dust impact risk to the spacecraft and their scientific instruments to aid hazard mitigation strategies.

This kind of assessment has a long history dating back to the assessment of the flyby at comet 1P/Halley by the VEGA spacecraft, for example (Sagdeev et al. 1982). Other studies have been performed to quantify the dust environment of comet 9P/Tempel 1 for the Deep Impact mission (Lisse et al. 2005), for the Rosetta mission to comet 67P/Churyumov-Gerasimenko (Agarwal et al. 2007), and most recently for potential missions to active centaurs, for example (Fink et al. 2021).

Though such models describing the dust environment for space missions are not new, the CI mission is in a unique situation. It is the first mission for which the precise target (a specific comet) of the mission could remain unknown until after launch. Naturally this is a particular problem for determining the expected dust coma because of the many parameters with a broad range of possible values. The primary problem is thus not what model to use, but what values for the parameters to assume in a model of the inner dust coma. Furthermore, because the different assumptions about the dust mass loss from cometary nuclei are strongly interdependent Marschall et al. (2020), it is not even obvious a priori which set of parameters represent the best- and worst-case scenarios.

Here, we present the Engineering Dust Coma Model (EDCM) which has a limited number of input parameters, but it retains general physical realism of the dust distribution in the inner coma. The purpose of this model is to make predictions of which dust the three spacecraft of CI will encounter during the active phase of the mission.

It is important to note that due to the uncertainty of the target comet and therefore its parameters, we have to use a new approach to build the model. Instead of defining one set of parameters that shall represent a best-, nominal-, or worst- case scenario, we chose ranges for each parameter of the model based on our knowledge of comets such as 1P/Halley and other comets. All self-consistent combinations within those ranges will be run through our model to give a prediction of all possible coma environments within parameter space. This ensemble of dust environments is subsequently statistically evaluated to determine a probabilistic distribution of possible conditions which the spacecraft might encounter.

We describe the general philosophy behind the model and the EDCM in Sect. 2 and how it was calibrated (Sect. 3). Section 4.1 presents the results for the spacecraft of the CI mission. We also describe how the presented results can be scaled to different flyby trajectories (Sect. 6). We end this work by providing a discussion and conclusions (Sect. 7).

2 The model

2.1 A shift in the modelling approach

For the assessment of the dust coma environment, missions to particular comets (e.g. previous missions to comets 1P/Halley, or 67P/Churyumov-Gerasimenko) had a significant advantage over CI. Their targets were known before launch and ground-based data with information about the targets could inform the expected dust coma environment (Sagdeev et al. 1982; Lisse et al. 2005; Agarwal et al. 2007; Fink et al. 2021).

In contrast the target of CI will be selected shortly before or even after the spacecraft have been launched into space (Snodgrass & Jones 2019). Therefore, it is necessary to assess the dust environment for a range of conditions which CI might encounter.

Two issues are of particular importance in our case. First, because the target is unknown, or at the least it’s activity within the water ice line and in particular at the encounter distance of 1 au, so are the properties about its activity. We don’t know how large the nucleus will be, what the gas and dust production rates are around 1 au, or what the dust size distribution is. Second, more fundamentally, even if some of these properties were known the different dust properties are inter-dependent in a way that makes it unclear a priori which set of parameters represent a best-, nominal-, or worst-case scenario.

We illustrate this with a simple example of just three parameters influencing the dust environment: the gas production rate, the dust-to-gas mass flux ratio, and the dust size distribution. Each of these individually trivially affects the dust densities within the coma. For a given dust-to-gas mass flux ratio and dust size distribution, an increase in the gas production rate will result in an increase of the dust densities in the coma. Likewise, for a given gas production rate and dust size distribution, an increase in the dust-to-gas mass flux ratio will increase the dust densities in the coma. Finally, the abundance of a specific size can be directly determined by the dust size distribution provided a given gas production rate and dust-to-gas mass flux ratio.

While these relationships are intuitive, combinations of them no longer have an obvious result. For example, it is unclear which of the following two combinations has more large particles:

1. a coma with a high gas production rate, low dust-to-gas mass flux ratio, and size distribution dominated by large particles,

2. a coma with a high gas production rate, high dust-to-gas mass flux ratio, and size distribution dominated by small particles.

Of course this depends on the exact values of each of these three parameters, but it illustrates the fundamental dilemma. And this is not the full extent of the problem. In reality, we are not confronted by only three inter-dependent parameters, but, as we discuss below, by at least 11.

It is therefore unclear how to properly define a set of parameters which shall correspond to a certain scenario (best, nominal, worst). For this reason, we have decided to present a new approach to this problem. While previous work could rely on a narrow set of input parameters we turn this around.

For each parameter of the model, we consider a broad range of values which are compatible with known observations of comets. Within this parameter space we select all possible combinations that are self-consistent. For example a parameter combination can be excluded if it were to result in an Afρ higher than the upper limit assumed here. This set of self-consistent parameter combinations will form the basis of the model.

Each self-consistent parameter combination is then run through our dust coma model and results in a unique coma environment. The resulting set of coma environments can then be analysed in a statistical manner to determine, for instance, a nominal case.

Our approach thus turns around the process of predicting the coma environment. Instead of defining a nominal case, we run all plausible parameter combinations and let the model tell us which dust densities are expected to be encountered.

Fig. 1 Sketch of the components of the EDCM and which input parameters feed into which step.

2.2 Overview of the model components

The EDCM is composed of three components:

• the dimensionless dust particle dynamics which calculates the spatial distribution of dust,

• the scaling of dust volume densities which makes the conversion to physical units and performs column integration taking into account dust scattering and a size distribution in order to obtain the brightness Afρ (defined by A’Hearn et al. 1984, it is the product of grain albedo (A), the cross sectional filling factor of the dust within the photometric aperture (f), and the circular aperture’s radius (ρ). For idealised circumstances (constant production, radial outflow) it is proportional to the dust mass loss rate (Fink & Rubin 2012) and make the selection of self-consistent cases,

• the observable extraction that determines the number density encountered along a specific spacecraft trajectory and performs the estimation of the probability of encountering the respective case.

The three components are depicted in magenta, blue, and green, respectively in Fig. 1 and will be described in detail in the following sections. Which parameters are required by which component are illustrated as well.

Table 1 lists the parameters, respective symbols, and ranges used for the target of the CI mission. The ranges used have been chosen such that CI will be able to survive a fly-by similar to the one by Giotto at comet 1P/Halley (see discussion later on).

Some parameters are of note. For now we have fixed the heliocentric distance to 1 au (where the CI fly-by will occur) thus also fixing the nucleus surface temperature at the sub-solar point. The model will be valid for small variations in the heliocentric distance (±20%). This is because of the large range of possible gas production rates we are allowing in the model. To some extent a lower heliocentric distance resulting in a higher production rate is therefore indirectly contained in the model. The nucleus is assumed to be a homogeneous sphere of radius R_N. We assume 2 and 5 km to cover the range of comets such as 67P (Preusker et al. 2017) to that of comet Halley (Keller et al. 1987; Wilhelm 1987). The bulk nucleus density was assumed as low as for objects like the cold classical Arrokoth (Spencer et al. 2020) up to the canonical value of 500 kg m⁻³ for comets (Preusker et al. 2017)

The gas emission distribution, q_g(φ), is assumed to be proportional to cos^k(φ), where φ is the solar zenith angle, and the power k takes the values 1, 2, and 3 (which corresponds to the different mechanisms of gas production from the nucleus). For example Zakharov et al. (2021a) showed that at 1 au sublimation from the surface results in an emission distribution closely following the case where k = 1 (Fig. 2 in Zakharov et al. 2021a). If, on the other hand, sublimation occurs inside the surface layer, then the distribution changes to the case where k = 3. In this sense, k can be interpreted as the parametrisation of where sublimation occurs.

The gas production rate is assumed for an equivalent sphere with a radius of 1 km and then scaled according to the actual nucleus radius used. The lower limit corresponds to a comet with an activity equivalent to comet 67P at 1.4 au inbound (Combi et al. 2020; Marschall et al. 2020) and the upper limit to roughly twice that of comet 1P/Halley during the Giotto fly-by (Krankowsky et al. 1986). Night-side activity is introduced as uniform emission with respect to φ and quantified by the ratio of the gas production rate at the anti-solar point to the production rate at the sub-solar point, $q_{\text{g}}^{180°}/q_{\text{g}}^{0°}$.

We allow for a large range of dust particles from micron to decimetre sized which have been observed at comets (McDonnell et al. 1997; Price et al. 2010; Mannel et al. 2016, 2019; Merouane et al. 2017; Ott et al. 2017; Fulle et al. 2016; Levasseur-Regourd et al. 2018; Della Corte et al. 2019). The slope of the size distribution has also been chosen according to previous measurements (e.g. McDonnell et al. 1997; Price et al. 2010; Merouane et al. 2017). The brightness of the coma, expressed in Afρ was chosen between 2800 cm and 28 000 cm with the upper limit being slightly higher than that of comet 1P/Halley during the Giotto encounter (Fink 1994; Schleicher et al. 1998).

2.3 The dust particle dynamics

We assume that the dusty-gas coma is formed by the gas sublimating from the nucleus (from the surface and/or from the interior), and solid particles (mineral or/and icy) released from the nucleus with zero initial velocity and entrained by the gas flow. It is assumed that the dusty-gas flow is coupled in one way only – the gas drags the dust (i.e. the presence of dust in the coma does not affect the gas motion), and that the dust particles do not collide with each other (the assumption of no feedback on the gas for most comae was articulated in Tenishev et al. 2011). The dust particles are assumed to be isothermal, spherical and internally homogeneous with invariable size and mass. To compute the dust distributions we adopt a frame attached to the nucleus but not rotating with it. Under these assumptions, dust particle motion is governed by the following equations (Crifo et al. 2005): $$m_{\text{d}}\frac{\text{d}{\upsilon }_{\text{d}}}{\text{d}t}=F_{\text{A}}+F_{\text{G}}+F_{\text{S}}+F_{\text{I}}$$(1) $$\frac{dr}{dt}={\upsilon }_{\text{d}},$$(2)

where t is the time, r and υ_d are the particle radius and velocity vectors, m_d is the particle mass, F_A is the aerodynamic force, F_G is the nucleus gravity force, F_S is the solar force, and F_I is the inertia force.

For isothermal, spherical, and internally homogeneous non-sublimating grains, the aerodynamic force degenerates to a drag force of $$F_{\text{A}}=\frac{{\left({\upsilon }_{\text{g}}-{\upsilon }_{\text{d}}\right)}^2}{2}\frac{\left({\upsilon }_{\text{g}}-{\upsilon }_{\text{d}}\right)}{\left|{\upsilon }_{\text{g}}-{\upsilon }_{\text{d}}\right|}{\rho }_{\text{g}}{\sigma }_{\text{d}}{C}_{\text{D}}\left(s,\frac{T_{\text{d}}}{T_{\text{g}}}\right),$$(3)

where σ_d and T_d are the particle cross section and temperature; ρ_g, υ_g, and T_g are the gas density, velocity, and temperature; and C_D is the drag coefficient. For a spherical particle σ_d = πa² and $m_{\text{d}}=\frac{4}{3}\pi a^3{\varrho }_{\text{d}}$, where a is the particle radius and ϱ_d is the particle specific density. If the particle size a is much smaller than the mean free path of the gas molecules (typical case in the coma), the drag coefficient C_D can be estimated by the following ‘free-molecular’ expression for a spherical particle (see Bird 1994): $$\begin{array}{lll}{C}_{\text{D}}\left(s,\frac{T_{\text{d}}}{T_{\text{g}}}\right)\hfill & =\hfill & \frac{2s^2+1}{s^3\sqrt{\pi }}\text{exp}\left(-s^2\right)\hfill \\ \hfill & \hfill & +\frac{4s^4+4s^2-1}{2s^4}\text{erf}\left(s\right),\hfill \\ \hfill & \hfill & +\frac{2\sqrt{\pi }}{3s}\sqrt{\frac{T_{\text{d}}}{T_{\text{g}}}}\hfill \end{array}$$(4)

where $s=\left|{\upsilon }_{\text{g}}-{\upsilon }_{\text{d}}\right|/\sqrt{2T_g{k}_{\text{B}}/m_{\text{g}}}$ is a speed ratio, m_g is the mass of the molecule, and k_B is the Boltzmann constant.

The nucleus gravitational force F_G is $$F_{\text{G}}=m_{\text{d}g},$$(5)

where g is the gravitational acceleration in the field of the nucleus. For the spherical nucleus g = −GM_Nr/|r³| where G is the gravitational constant, M_N is the mass of the nucleus.

The solar force F_S which includes solar gravity and solar radiation pressure is given, for spherical grains, by $$F_{\text{S}}=\left(\frac{G{M}_{\odot }{m}_{\text{d}}}{r_{\text{h,d}}^2{r}_{\oplus }^2}-{\sigma }_{\text{d}}\frac{Q_{\text{ef}}{c}_{\odot }\xi \left(r\right)}{c{r}_{\text{h,d}}^2}\right)i_{\odot ,d},$$(6)

where M_⊙ is the mass of the Sun, Q_ef is the radiation pressure efficiency, c is the speed of light, c_⊙ is the solar energy flux at 1AU, r_h is the heliocentric distance of the nucleus (in AU), r_h,d is the dust particle heliocentric distance (in AU), r_⊕ = 1.496 × 10¹¹ m is the Earth orbit semimajor axis, i_⊙,d is the unit vector directed from the dust particle to the Sun, ξ(r) = 0 if r is inside the shadow, and ξ(r) = 1 if r is in sunlight.

In the adopted nucleus-attached (but non rotating) frame, the inertia force F_I is given (Landau & Lifshitz 1976) by the following: $$F_{\text{I}}=m_{\text{d}}\frac{G{M}_{\odot }}{r_{\text{h}}^2{r}_{\oplus }^2}{i}_{\odot ,\text{n}},$$(7)

where i_⊙,n is the unit vector from the nucleus centre of mass to the Sun.

Following Zakharov et al. (2018b) and Zakharov et al. (2021b) we introduce the following dimensionless variables: θ = T_d/T_s, ${\tilde{\upsilon }}_{\text{g}}={\upsilon }_{\text{g}}/{\upsilon }_{\text{g}}^{\text{max}},{\tilde{\upsilon }}_{\text{d}}={\upsilon }_{\text{d}}/{\upsilon }_{\text{d}}^{\text{max}},\tilde{r}=r/R_N,\tilde{t}=t/\text{Δ}t$, τ = T_d/T_g, ${\tilde{\rho }}_{\text{g}}={\rho }_{\text{g}}/{\rho }_{\text{s}}$. Here ${\upsilon }_{\text{g}}^{\text{max}}=\sqrt{\gamma \frac{\gamma +1}{\gamma -1}\frac{k_{\text{B}}}{m_{\text{g}}}{T}_s}$ is the theoretical maximum velocity of gas expansion, ρ_s and T_s are the gas density and temperature on the sonic surface (i.e. where the gas velocity is equal the local sound velocity $\sqrt{\gamma k_{\text{B}}{T}_{\text{g}}/m_{\text{g}}}$) at the sub-solar point, γ is the specific heat ratio, R_N is the characteristic linear scale (e.g. the radius of the nucleus), and $\text{Δ}t=R_{\text{N}}/{\upsilon }_{\text{g}}^{\text{max}}$. In addition, we introduce $\tilde{g}$ via normalising g on $G{M}_{\text{N}}/R_{\text{N}}^2$ (therefore, in the case of a spherical nucleus $\tilde{g}=-\tilde{r}/\left|{\tilde{r}}^3\right|$).

Using these dimensionless variables, Eqs. (1) and (2) can be rewritten as follows: $$\begin{array}{lll}\frac{\text{d}{\tilde{\upsilon }}_{\text{d}}}{\text{d}\tilde{t}}\hfill & =\hfill & {\tilde{\rho }}_{\text{g}}{\left({\tilde{\upsilon }}_{\text{g}}-{\tilde{\upsilon }}_{\text{d}}\right)}^2\frac{{\tilde{\upsilon }}_{\text{g}}-{\tilde{\upsilon }}_{\text{d}}}{\left|{\tilde{\upsilon }}_{\text{g}}-{\tilde{\upsilon }}_{\text{d}}\right|}{C}_{\text{D}}\left(s,\tau \right)\text{Iv}+\text{Fu\hspace{0.17em}}\tilde{g}\hfill \\ \hfill & \hfill & +\frac{R_{\text{N}}G{M}_{\odot }}{{({\tilde{\upsilon }}_{\text{g}}^{\max })}^2{r}_{\text{h}}^2{r}_{\oplus }^2}\times \left[\left(1-\frac{3c_{\odot }{r}_{\oplus }^2}{4G{M}_{\odot }c}\frac{Q_{\text{ef}}}{{\varrho }_{\text{d}}a}\xi \left(\tilde{r}\right)\right)\frac{r_{\text{h}}^2}{r_{\text{h},\text{d}}^2}{i}_{\odot ,d}-i_{\odot ,n}\right]\hfill \end{array}$$(8) $$\frac{\text{d}\tilde{r}}{\text{d}\tilde{t}}={\tilde{\upsilon }}_{\text{d}}.$$(9)

If the linear size of computational domain is much smaller than r_h, then r_h ≃ r_h,d and i_⊙,d ≃ i_⊙,n. In this case, Eq. (8) may be rewritten as: $$\frac{\text{d}{\tilde{\upsilon }}_{\text{d}}}{\text{d}\tilde{t}}={\tilde{\rho }}_{\text{g}}{\left({\tilde{\upsilon }}_{\text{g}}-{\tilde{\upsilon }}_{\text{d}}\right)}^2\frac{{\tilde{\upsilon }}_{\text{g}}-{\tilde{\upsilon }}_{\text{d}}}{\left|{\tilde{\upsilon }}_{\text{g}}-{\tilde{\upsilon }}_{\text{d}}\right|}{C}_{\text{D}}\left(s,\tau \right)\text{Iv}+\text{Fu\hspace{0.17em}}\tilde{g}-\text{Ro}\xi (\tilde{r})i_{\odot ,n}.$$(10)

As a consequence of this assumption, this paper neglects the solar tidal forces that affect significantly the dust motion out of the acceleration zone over long times (weeks and longer). For particles equal and larger than mm-sized these effects are important also just outside the dust acceleration zone. The rigorous computation in the heliocentric reference frame of the dust motion out of cometary Hill’s sphere can be found in Fulle et al. (1995).

The Eqs. (8) and (10) contain three dimensionless parameters: $$\text{Iv}=\frac{1}{2}\frac{{\rho }_{\text{s}}{\sigma }_{\text{d}}{R}_{\text{N}}}{m_{\text{d}}}=\frac{3}{8}\frac{{\rho }_{\text{s}}{R}_{\text{N}}}{a{\varrho }_{\text{d}}}=\frac{3Q_{\text{g}}^{\text{sph}}{m}_{\text{g}}}{32R_{\text{N}}a{\varrho }_{\text{d}}\pi \sqrt{T_{\text{s}}\gamma k_{\text{B}}/m_g}}$$(11)

and $$\text{Fu}=\frac{G{M}_{\text{N}}}{R_{\text{N}}}\frac{1}{{({\upsilon }_{\text{g}}^{\max })}^2}$$(12)

and $$\text{Ro}=\frac{1}{m_{\text{d}}{({\upsilon }_{\text{g}}^{\max })}^2}{R}_{\text{N}}\frac{{\sigma }_{\text{d}}{Q}_{\text{ef}}{c}_{\odot }}{c{r}_{\text{h}}^2}.$$(13)

Here $Q_{\text{g}}^{\text{sph}}$ is the total gas production rate (s⁻¹) of a spherical nucleus of radius R_N with the local gas production uniform over the surface and equal the gas production in sub-solar point.

We note that the theoretical maximum velocity of gas expansion could be expressed also in terms of a stagnation temperature T₀ and a heat capacity C_p as ${\upsilon }_{\text{g}}^{\text{max}}=\sqrt{2C_{\text{p}}{T}_0}$. In this way the gas velocity on the sonic surface is ${\upsilon }_{\text{s}}=\sqrt{2C_{\text{p}}{T}_0}\sqrt{\frac{\gamma -1}{\gamma +1}}$ and the local gas production is q = ρ_sυ_s. Therefore we can rewrite Iv as: $$\text{Iv}=\sqrt{\frac{\gamma +1}{2\left(\gamma -1\right)}}\frac{3q{R}_{\text{N}}}{8\sqrt{C_{\text{p}}{T}_0}a{\varrho }_{\text{d}}}.$$(14)

The second term in this equation is the reciprocal of a dimensionless similarity parameter characterising the ability of a dust particle to adjust to the local gas velocity introduced in Probstein (1969).

These parameters have the following meaning:

1. Iv represents the ratio of the gas mass present in a flow tube, with the cross section of the particle and a characteristic length R_N, to the particle mass. This parameter characterises the efficiency of entrainment of the particle within the gas flow (i.e. the ability of a dust particle to adjust to the gas velocity);

2. Fu represents the ratio of the comet surface gravitational potential to the flow thermodynamic potential (enthalpy, C_pT₀). This parameter characterises the efficiency of gravitational attraction;

3. Ro represents the ratio of the specific work done by the solar pressure force on the characteristic length R_N to the flow thermodynamic potential. This parameter characterises the contribution of solar radiation pressure.

In order to define Iv, Fu, Ro it is necessary to know: m_g, γ, ρ_s, T_s, R_N, M_N, a and ϱ_d (or σ_d and m_d), Q_ef and r_h.

We have computed the dust distribution within the cometary coma up to 1000 nucleus radii (R_N). It uses a minimal number of parameters for the description of a cometary dust coma, while keeping it physically realistic. This model physically consistently takes into account the expanding nature and asymmetry of the gas coma (caused by gas production modulated by solar radiation) and considers the dust dynamics driven by the gas drag force, force from the nucleus gravity, and solar radiation pressure.

A series of general assumptions were made to simplify the model. For each of the simplifications we also outline the resulting limitations that need to be appreciated.

• The nucleus shape is assumed to be spherical. The dust dynamical model can therefore not reproduce possible heterogeneity within the coma that arise due to topography (and often occurs only at particular orientation of the nucleus with respect to the Sun). Nevertheless, when the parameters of the coma of a complex shape nucleus are averaged over a full comet rotation the resulting average coma distribution resembles closely the one from a spherical nucleus (e.g. Zakharov et al. 2021a).

• The gas is assumed to be an ideal perfect gas. The gas flow in the coma is described by the gas-dynamic Euler equations which express the conservation of the mass, momentum and energy in the flow of an ideal perfect gas.

In other words these equations describe the motion of an equilibrium fluid flow without viscous dissipation and heat conductivity. The real atmosphere of a comet contains rarefied non-equilibrium regions as well.

The transfer of thermal energy into kinetic energy in a rarefied flow is less efficient than in a fluid flow, therefore the rarefied flow accelerates slower. Due to the presence of viscous dissipation flow structures like shocks are diffused in a rarefied flow.

Nevertheless, as was shown in many publications, the description of the flow from a spherical nucleus based on the fluid equations preserves general physical realism of the flow. As was noticed in Crifo et al. (2004), the existing results from comparisons of kinetic and fluid approaches (e.g. Crifo et al. 2002, 2003) show that the Navier-Stokes equations and even the Euler equations provide acceptable solutions over practically the whole day-side coma of observable comets. Two restrictions are to be made, however: (1) the immediate vicinity of the nucleus surface (several R_N), and (2) the outer reaches of the coma, where dissociation products are dominant, it is presently not known to which accuracy these equations represent the real situation.

An important merit of the Euler equations is that they don’t depend on the flow rarefaction and therefore the solutions can be scaled for different production rates, the main reason why we use them in the present model.

• It is assumed that the dust does not influence the gas flow (i.e. no back-coupling of the dust to the gas flow). This allows the separate, and sequential determination of the gas and dust flows. This is given for comets with low dust-to-gas ratios [~<10]. Nevertheless, for a comet with a high dust-to-gas ratio the model will be able to produce reliable predictions.

• It is assumed that the gas coma is constituted of one single species, H₂O. This assumption will be reasonable well satisfied for comets where H₂O is the dominant gas species (as e.g. for comet 67P, most comets at 1 au).

• It is assumed that there is no extended gas/dust source/sink in the coma. The model should not be applied to comets with a significant extended gas/dust source/sink, for example, emission of gas from dust particles within the coma, sublimation and destruction of dust particles etc.

• It is assumed that the gas and dust emission is smooth across the surface or that any inhomogeneities are blurred within the drag acceleration region r ≤ 10R_N. The gas and dust emission is modulated, for instance, by the solar zenith angle. This means that the emission is not dominated by a few very localised jets.

• The model is reliable for cometo-centric distances of 10R_N < r < 1000 R_N. The inner boundary is defined by the possible existence of ‘fine structures’ of the flow due to particularities of the surface topography. The outer boundary is defined by the size of the computational region. It is possible to extrapolate the data beyond this upper limit via introducing additional assumptions (e.g. sphericity of expansion etc.).

• The model is run for a heliocentric distance of 1 AU, the approximate distance of CI’s comet during the encounter. Result shall thus not be attempted to be scaled to vastly different heliocentric distances. The results remain valid for variations of the heliocentric distance that do not change the solar flux to the nucleus surface by more than a factor of 2 (i.e. the model is valid for the heliocentric distance range of CI).

• It is assumed that dust particles are spherical, homogeneous with invariable size and mass.

This simplification allows us to solve the equation of motion generally in (Iv, Ro, Fu) space and then re-scale these results to the respective physical values in a later step. For all details on this model we refer the reader to Zakharov et al. (2018b) and Zakharov et al. (2021b).

For the underlying gas dynamics model we used the results by Zakharov et al. (2021a) who have calculated the gas field by solving the Euler equations. The gas results are used to calculate the dynamics of spherical dust particles taking into account gas drag, nucleus gravity, and solar radiation pressure (Eq. (1); Zakharov et al. 2018b, 2021b). An important implicit assumptions to reiterate is that we assume that the dust does not have a back reaction effect on the gas flow.

The dust dynamics makes use of dimensionless variables (derived by combining the respective physical parameters from Table 1) to parametrise the dust dynamics as described in Zakharov et al. (2018a, 2021b) and reduce the number of numerical solutions. The result of this step is the dust number density and velocity in three-dimensional space.

At this point there is also the implicit assumption that the dust-to-gas mass production rate ratio, χ, is unity. These solutions therefore do not yet have an absolute scale and thus requiring the second component of the EDCM.

2.4 Scaling the model to physical parameters

In a first step we convert the dimensionless results given in (Iv, Ro, Fu)-space into physical units. This is done for each of the parameter combination in our parameter set.

Once this is done only one parameter has yet to be fixed. Until now the dust-to-gas ratio, χ, is still assumed to be one. χ will be chosen such that a certain brightness of the coma, Afρ, is reached. The larger χ the larger Afρ. But we need to identify what value of χ results in which Afρ. The calculation of Afρ for each parameter combination follows the approach described in Marschall et al. (2016).

First, the dust column density of an aperture of 20R_N is calculated. For points outside the simulation domain (>1000R_N) a 1/r² extrapolation is applied. The column densities are then weighted with a power law (n ~ a−^β) and converted into reflectance using the scattering model of Markkanen et al. (2018) as shown in Marschall et al. (2020).

Second, the resulting reflectance can then be used to calculate the Afρ as explained in Gerig et al. (2018). The absolute scaling χ can then be determined by linearly scaling the model Afρ to the desired Afρ, i.e. if the model Afρ = 100 cm then an actual coma with Afρ = 200 cm is achieved with χ = 2.

Two additional constraints are added at this point:

1. χ < 10, and

2. −29.7 < $\text{log}\left(\text{Af}\rho /{\text{Q}}_{{\text{H}}_2\text{O}}\right)$ < −27.65.

The second constraint comes from ground-based observations of comets (A’Hearn et al. 1995; Schleicher & Bair 2016). These measurements show that comets don’t exhibit an arbitrary range of $\text{Af}\rho /{\text{Q}}_{{\text{H}}_2\text{O}}$.

Figure 2 shows the ranges for typical (yellow) and depleted (orange) comets according to Table VI in A’Hearn et al. (1995). The range covered in this work corresponds to a range spanning the ones from A’Hearn et al. (1995); Schleicher & Bair (2016), and Schleicher & Bair (private communication; Oort cloud comets).

We have chosen to calibrate the model on the Giotto data (see Sect. 3) for the comparison but employed a scattering model calibrated to comet 67P (Markkanen et al. 2018). This introduces some bias because the phase functions of 1P and 67P differ significantly, especially at large phase angles (Fig. 3). We note that these two phase curves have been measured differently (1P from ground, 67P from inside the coma). The Halley–Marcus phase curve (or Schleicher-Marcus phase curve) is a composite. The part originating from Halley only extended up to 55° (Schleicher et al. 1998), and the remainder is based on a Henyey–Greenstein function that is consistent with some near-Sun comets (Marcus 2007)¹.

Because we employ the scattering model by Markkanen et al. (2018) for the calculation of the model Afρ, a global scaling factor of 23 is applied to bring in line the Afρ measurements of 1P with the Giotto dust densities. Therefore, by definition our model fits the Halley data when the appropriate input parameters for Halley are chosen. Thus all results can be understood as relative to Halley. In this sense if the model is applied to, for example, a comet that is otherwise identical to Halley but has a dust-to-gas ratio that is double that of Halley, the model will predict densities that are double those measured with Giotto.

Fig. 2 $\text{Log}\left(\text{Af}\rho /{\text{Q}}_{{\text{H}}_2\text{O}}\right)$ from A’Hearn et al. (1995) for typical (yellow) and depleted (orange) comets as well as the range of the EdCM.

Fig. 3 Phase curves of 1P (Schleicher) and 67P (Bertini et al. 2017) are shown.

Fig. 4 Giotto flyby trajectory relative to the cometary nucleus and Sun is shown, including offset angles and positions forthelast4usefulHMC images (from Keller et al. 1987).

2.5 Extracting observables

Now that our set of dust comae have been scaled to the respective physical units we are left with the final component of the EDCM. We can now extract the number density encountered along a specific spacecraft trajectories for coma in our set of input parameters. Again, for points outside the simulation domain (>1000R_N) a 1/r² extrapolation is applied.

The trajectories are defined by three parameters (Table 1). The closest approach (CA) distance, the trajectory angle in the meridional plane, α_T, and the trajectory angle in the equatorial plane, β_T. The nominal CA distances for the three CI spacecraft A, B1, and B2 are 1000, 500, and 200 km respectively. The two angles are depicted as α and β in Fig. 4.

For our nominal cases we have assumed that the fly-bys of all three spacecraft occur in the same plane, cross the Sun-comet line at 90°, and thus go through 0° phase angle. Though these will not be the actual trajectories of the three spacecraft our assumed trajectories are the more conservative assumption. Trajectories out of the assumed plane result in smaller dust densities (see Sect. 6). As soon as the trajectories are known the model can be adjusted to the actual trajectory. We have extracted the dust number density along the trajectory as well as the integrated number of particles along the trajectory (i.e. the column density) for each dust size separately.

In a final step, we analyse the results from the set of dust comae. In particular we have calculated the median column density for each dust size as well as certain percentiles (specifically 5%, 10%, 25%, 75%, 90%, 95%). This gives a measure for the scatter of cases and thus a measure for the uncertainty of possible outcomes.

3 Model calibration

To calibrate the EDCM we have used the case of comet 1P/Halley which was visited, amongst others, by the Giotto mission.

The three-dimensional shape of the comet Halley’s nucleus was derived from the combination of Vega and Giotto images. A self-consistent model has been produced using an ellipsoid with 8.0, 8.2, and 16 km for the three axes Wilhelm (1987). The surface is about 400 km² and volume 580 km³. The spin vector points about 60° south of the comet to sun direction and it is about perpendicular to the long axis of the nucleus. A spin period of 54 ± 1 h matches the data from the three spacecraft.

In situ measurements Krankowsky et al. (1986) showed a total gas production rate of 6.9 × 10²⁹ molecules s⁻¹ with the predominance of water vapour ~80% (and 20% more volatile compounds) by volume and the photodestruction scale length for H₂O of 3.9 × 10⁴ km. Keller et al. (1987) noted that the Halley Multicolour Camera (HMC) images have revealed that the activity of comet Halley’s nucleus was concentrated towards the sun and that the variability of the comet during Multi Detector Mode was seen to be small (the minimum-maximum variation is considerably less than a factor of 2).

The results of analysis in Rickman (1989) provide the mass of the nucleus 1.3 ± 0.3–3.1 ± 0.4 × 10¹⁴ kg, the mean density of the nucleus 280 ± 100–650 ± 190 kg m⁻³, and the dust-to-gas ratio 0.3–1.

The EDCM uses axially-symmetric numerical solutions. Therefore, some trajectories which are different in three-dimensional space might be identical in the axially-symmetric frame. For the model with a spherical nucleus we use the radius of 5 km, a nucleus mass of 2.7 × 10¹⁴ kg, the total gas production rate (5.1–5.5) × 10²⁹ s⁻¹, and Afρ of 28 000 cm (Fink 1994; Schleicher et al. 1998). The other parameters have been left as free parameters which results in 84 parameter variations.

The results of the EDCM for 1P/Halley during the Giotto fly-by are shown in Fig. 5. The data for dust masses below 10⁻¹² kg is an lower limit measurement and should thus not be considered in the comparison. Apart from that the data lies within the confidence interval of the model for most of the mass range. In the mass range from 10⁻¹² to 10⁻⁸ kg the model matches well the slope of the Giotto data. Because the size distribution for Halley has a break at ~3 × 10⁻⁸ kg our model cannot match that behaviour as it assumes an unbroken power law. But such a power law is implicitly contained in the model because a variety of slopes are simulated and thus are reflected in the uncertainty.

Fig. 5 Total number of particles during encounter as a function of mass for the EDCM Halley case and the data by McDonnell et al. (1997).

4 Results

4.1 Predictions for the Comet Interceptor spacecraft

The different values for the input parameters listed in Table 1 result in ~23 000 self-consistent combinations (i.e. possible cases of the dust coma). To arrive at this set of cases, each possible combination of parameter space was checked, for example, for each nucleus radius, two different nucleus bulk densities; for each of those combinations, three different gas emission distributions; for each of those combinations, three different night side activity levels; for each of those combinations, twelve different dust size distributions; among others. Multiplying all degrees of freedom in Table 1 results in a total of over 9 million combinations. But most of these combinations are not self-consistent and were thus discarded, if they violated a specific constraint, such as χ > 10 or, for instance, a given dust production rate resulted in larger Afρ than we consider for this work.

As described in Sect. 2, we have extracted for each dust size bin the number of particles encountered along the spacecraft trajectory.

At each point along the trajectory we calculated the median number of particles predicted by the set of all parameter combinations as well as the 5th, 10th, 25th, 75th, 90th, and 95th percentile. Each combination within parameter space of the model is assumed to be equally likely. The results for four size bins and spacecraft A is shown in Fig. 6. The shaded areas illustrate the variation in predicted number of particles based on the variation of the input parameters. These ranges thus reflect to a large degree the uncertainty of our knowledge of the future target of CI. As the dust size increases the expected number of particles decreases but the uncertainty increases. Further, the spike in particles around CA highlights that most particles are encountered very close to CA. For example, from cometo-centric distances of 10 000 km to CA at 1000 km the dust densities increase by roughly 2.5 orders of magnitude.

An important thing to note for the number density is that there does not exist a global median case along the trajectory. This means what is the median case at far distances may not be the median case at CA. This also means that the median column density is not equal to the integration of the median number density. Therefore, we have calculated the percentiles at each distance from the nucleus independently.

We have further integrated the total number of particles per square meter expected along the entire fly-by. The results for spacecraft A are shown in Fig. 7. The figure shows the number of particles per 1.5 decade in dust mass. This means that the value, for instance, given at 1 μm or 3 mm includes the number of particles in the interval from 1–3 μm or 3–10 mm respectively. As another example we expect in the median case ~1 particle per square meter along the entire fly-by trajectory of spacecraft A in the size range of 300–1000 μm.

Additionally, Fig. 7 indicates in orange the median case for an encounter of spacecraft A at comet Halley (as defined Sect. 3). The main reason that the median case for spacecraft A lies below a Halley encounter is that the EDCM model assumes a larger range of possible coma environments which turns out to slightly skew to scenarios that produce dust comae with lower dust densities than Halley exhibited. In a simplified picture, the median target comet in the EDCM is slightly less active than Halley.

Finally, it is worth mentioning that in contrast to the number density there do exist global percentile cases (e.g. a global median case). All results for all three CI spacecraft and additional supplementary results are available online².

Fig. 6 Number of dust particles according to the EDCM along the spacecraft trajectory of spacecraft A as function of cometo-centric distance. The shaded areas show different percentile ranges within which cases fall. Closest approach for S/C A is assumed at 1000 km.

Fig. 7 Total number of dust particles encountered according to the EDCM along the spacecraft trajectory of spacecraft A as function dust mass. The shaded areas show different percentile ranges within which cases fall. Additionally the orange curve shows the predicted median densities for a Halley type comet (according to the definition in Sect. 3).

4.2 Predictive strength of parameters

Based on the ensemble of self-consistent cometary comae described in Sect. 4.1 we can explore which parameters are predictive of the dust densities within the coma. Generally speaking, the parameters we have found to be most predictive of the fly-by column density are the combination of: Afρ, the power law exponent of the dust size distribution, β, and the dust production rate. Importantly the predictive power of these three parameters is not equal for all dust sizes.

Figure 8 shows the dust column density for dust particles with radii between 1 and 3 μm at the Comet Interceptor spacecraft A for different combinations of Afρ, the power law exponent, β, and the dust production rate. Afρ by itself has the highest predictive power in this case. The dust column density for a given Afρ is constrained to ~half an order of magnitude in all cases except for the lowest Afρ (left panels in Fig. 8). For any given dust production rate (centre panels in Fig. 8) or power law exponent (right panels in Fig. 8) by themselves the possible column densities span at least one order of magnitude in most cases. Afρ in combination with either the power law index or the dust production rate will tightly constrain the column density.

The behaviour for large particles (e.g. 10–30 cm) is very different from the one described above for small particles. In this case Afρ is the most unpredictive of the three parameters (left panels in Fig. 9). For a given Afρ the possible dust column densities span ~5 orders of magnitude. For large particles the power law index is generally the most predictive parameter. For a given β alone (right panels in Fig. 9) the expected dust column densities can be predicted to ~1.5 orders of magnitude accuracy. The dust production rate by itself (centre panels in Fig. 9) has variable predictive power. The predicted range of dust column densities is largest at small dust production rates and then decreases with increasing dust production rates. The best predictions of the dust column density of large particles can be achieved with a combination of the power law exponent and either Afρ, or the dust production rate (left panels in Fig. 9). By adding one of these constraints to the power law exponent the predicted dust column densities can be determined down to an accuracy of ~half an order of magnitude.

Often Afρ is used to estimate the dust production rate of a comet (e.g. Afρ[cm] ~Q (kg s⁻¹), A’Hearn et al. 1995). With our ensemble of dust comae we can test how large the associated uncertainties are. Figure 10 shows the production rate as a function of Afρ and the power law index. First, it is quite clear that Afρ is a rather poor predictor of the dust production rate. For a given Afρ the dust production rate spans ~2.5 to ~4 orders of magnitude. This illustrates why Afρ by itself should not be used to estimate a comet’s dust production rate. The estimates of the dust production rate can however be greatly improved when the power law index of the size distribution is known. In this case the uncertainty in the dust production rate can be reduced to ~1 order of magnitude.

Finally, for completeness we show in Fig. 11 that the other parameters in the model are not predictive of the dust column density. There are some minor trends, for instance, in the fraction of night side activity (bottom centre panel of Fig. 11). In this case this minor correlation is easily understood. The column density considered here is located on the day side (i.e. where the CI fly-by of spacecraft A will occur). Therefore, naturally as more activity occurs on the night side the relative amount of dust on the day side will decrease.

Fig. 8 Dust column density for dust particles with radii between 1 and 3 μm at the Comet Interceptor spacecraft A is shown as a function of different combinations of Afρ, the power law exponent, β, and the dust production rate.

Fig. 9 Dust column density for dust particles with radii between 10 and 30 cm at the Comet Interceptor spacecraft A is shown as a function of different combinations of Afρ, the power law exponent, β, and the dust production rate.

4.3 Dust speeds and impact angles

We would be remiss not to comment on the dust speeds and what those would mean for the impact angle on the spacecraft. With a simple estimates we want to illustrate that the spacecraft motion with respect to the comet will dominate the flux direction of the dust onto the spacecraft in almost all cases. The terminal dust speeds cannot exceed the gas speed with an upper limit of ~ 1000 m s⁻¹. In the case of isotropic gas expansion from a spherical nucleus the terminal dust velocity, υ_d,term, can be calculated from the approximation of the numerical solution: $${\upsilon }_{\text{d},\text{term}}=\frac{58.903\sqrt{T_{\text{N}}}}{1+0.6\sqrt{1.5/\text{Iv}}},$$(15)

where $$\text{Iv}=4.038\times {10}^{-29}\frac{Q_{\text{g}}}{R_{\text{N}}{\varrho }_{\text{d}}a\sqrt{T_{\text{N}}}},$$(16)

and where T_N is the nucleus temperature (to be set to 300 K for our purposes here), Q_g is the gas production rate (molecules s⁻¹), R_N is the radius of the nucleus (m), and ϱ_d is the bulk dust density. The left panel of Fig. 12 shows the terminal dust speeds for the maximum and minimum gas production rates considered in the EDCM for CI. Particles larger than ~600 μm are always slower than 100 m s⁻¹.

From the terminal dust speed we can derive the maximum impact angle of a dust particle with respect to the ram-direction of the spacecraft given a certain fly-by velocity. Currently fly-by velocities between 10 km s⁻¹ and 70 km s⁻¹ are considered for CI. The most likely fly-by speed is expected to be 50 km s⁻¹.

The right panel of Fig. 12 shows this angle which is always smaller than ~0.5° for particles with masses larger than 1 mg. The maximum impact angle occurs for the smallest particles (1 μm) in the highest gas production rate case and is ~4.3°. This illustrates that the impact angles is dominated by the spacecraft fly-by velocity. Such impacts can nevertheless cause some issues for attitude control, as is shown in Haslebacher et al. (2022).

Fig. 10 Dust production rate as a function of Afρ and the power law exponent, β, of the dust size distribution. The values of Afρ are slightly offset depending on the power law exponent to show the overlap resulting from β. By itself Afρ is a poor predictor of the dust production rate of a comet.

5 Guiding mission design and safety

In the introduction, we highlighted five areas that motivate this modelling approach. Here we revisit these points to show how the EDCM is/will be used to guide mission design and safety.

5.1 Trajectory design and spacecraft shielding

The choices of spacecraft trajectory and spacecraft shielding are interwoven. A more distant closest approach distance or an encounter through the anti-solar point rather than the sub-solar point would result in fewer particle impacts. How much the dust densities can be reduced through changes in the fly-by trajectory is detailed in Sect. 6. Therefore, there is a trade-off between the amount of shielding needed and how close the spacecraft can fly to the comet. The CI spacecraft are designed to survive a Halley/Giotto type comet encounter. This corresponds to the median case of the EDCM. The model is being used to assess the risk of particle hits, in particular the number of particles in the different size bins. The shielding of the CI spacecraft have been sized following this criterion.

5.2 Performance of scientific instruments and star trackers

One of the issues during operations of Rosetta at comet 67P was triggered by star trackers falsely identifying dust particles as stars (Accomazzo et al. 2017). Modern star trackers in combination with the EDCM dust densities will allow for in-depth bench-marking and thus safe operations.

Further, the dust accumulation rate, in the ram-direction, can be calculated using the column densities of the model. This is particularly important for instruments such as the Dust Impact Sensor and Counter (DISC) which are placed on spacecraft A and probe B2. DISC will measure the momentum and mass of the particles when they impact its sensitive surface.

5.3 Attitude control

The impact of dust particles onto the spacecraft’s surface will impact attitude control and thus the pointing of the instruments. As shown by Haslebacher et al. (2022) this can be significant and is a potential concern, for example, for keeping the nucleus within the field of view of the COmet CAmera (CoCa) on spacecraft A. They found that a total change in angular velocity that needs to be corrected can reach 1–10° s⁻¹ resulting in a median shift of the target object on the CoCa detector of more than 10 pixels. Haslebacher et al. (2022) used an idealised coma model and thus the updated values of the EDCM will allow for a refined assessment of this issue and thus provide better information on the need and scope of mitigating strategies.

6 Scaling the results

The nominal trajectories (Table 1) of the model only vary in CA distance but not in the two angles α_T and β_T (α and β in Fig. 4) which are both assumed to be zero. Because this will not be the case for all three spacecrafts of CI we have tested the model to the sensitivity of changes in α_T, β_T, as well as the CA distance, r_CA. Furthermore, the results might be needed for a different dust size or size interval than we have used in this work.

6.1 Variations in CA distance

First, we explored variations in the CA distance. For a force free, radially expanding coma the column density scales with 1/CA. Or put another way, the product of the column density and the CA distance is constant. Figure 13 shows this for the three S/C and reveals that this scaling works well for S/C A & B1 but not perfectly for S/C B2 which results in densities roughly a factor of 2 higher than expected from that scaling of the other two S/C. This is likely due to the fact that S/C B2 gets very close to the nucleus (down to 20 R_N in some of the cases considered).

In a sense S/C B1 is the nominal case as it is not too close to the nucleus and not too far from the nucleus (where we need to extrapolate the dust results from our modelling domain and thus ‘force’ a 1/CA behaviour). The medians of the different S/C are all still within the centre 50 percentiles of the cases. Thus the uncertainty from the unknown parameters of the model is much larger than the uncertainty of this simple scaling.

Also the assumed differences compared here are very large, therefore for small changes in the CA distance the results of the model can be scaled with $$N_{\text{co}1}\left(r_{\text{CA}}\right)=N_{\text{co}1}\left(r_{\text{CA};0}\right)\frac{r_{\text{CA};0}}{r_{\text{CA}}},$$(17)

where r_CA;0 is the reference CA distance. The number densities, Nnum, scale respectively with the square of the closest approach distance, i.e. $$N_{\text{num}}\left(r_{\text{CA}}\right)=N_{\text{num}}\left(r_{\text{CA};0}\right){\left(\frac{r_{\text{CA};0}}{r_{\text{CA}}}\right)}^2.$$(18)

Fig. 11 Dust column density for dust particles with radii between 100 and 300 μm at the Comet Interceptor spacecraft A is shown as a function of different combinations of the gas production rate, dust bulk density, surface activity distribution, night time activity, and radiation pressure efficiency. These quantities are not predictive of the dust column density.

Fig. 12 Dust speed and impact angles for different dust sizes and activity levels. Left panel: terminal dust speed, υd,term, according to Eq. (15) for two gas cases. The lowest gas production rate case considered in the EDCM in purple (min Qg case, where RN = 2 km and Qg = 155 kg s−1) and the highest gas production rate case in orange (max Qg case, where RN = 5 km and Qg = 16 434 kg s−1). Right panel: maximum angle of impact of dust particles with respect to the ram-direction of the spacecraft. The solid lines are for a fly-by speed of 10 km s−1 and the dashed lines for 70 km s−1. The colours are the min. and max. gas production cases as in the left panel. The corresponding radii of particles with masses 1 mg, 10 mg, and 100 mg are indicated by the grey lines.

Fig. 13 Column density multiplied with the CA distance as a function of dust mass.

6.2 Variations in α_T and β_T

Next we explore changes in α_T (α in Fig. 4). Changes to only α_T (β_T = 0°) merely change the phase angle at CA but retain a symmetric trajectory around CA. We have run the EDCM for α_T = 0°, 15°, 30°, 45°, 60° for S/C B2. Figure 14a shows the normalised column density to the values of each dust sizes column density at α_T = 0° as well as the size averaged relative column density. The column densities scale roughly as a function of the cosine but barely reach a factor of 2 decrease at α_T = 60°, which again is much less than the uncertainties inherent to the EDCM from the large uncertainties of the input parameters. We can therefore use the rough empirical relation $$N_{\text{co}1}\left({\alpha }_{\text{T}}\right)=N_{\text{co}1}\left(0\right){\cos }^{0.8}\left({\alpha }_{\text{T}}\right)$$(19)

to scale the nominal results to different values of α_T.

Finally, we varied β_T (β in Fig. 4 and not to be confused with β the power law) which not only changes the phase angle at CA to β_T (if α_T = 0°) but also skews the trajectory such that the asymptotic phase angles are 90 ± β_T. If α_T = 0°, as assumed here, the S/C will usually go through zero phase shortly after/before CA. Figure 14b shows the normalised column density to the values of each dust sizes column density at β_T = 0° as well as the size averaged relative column density. The dependency of the column density on β_T is weaker than that of α_T. The maximum decrease, at β_T = 60°, is ~15%. We can use the rough empirical relation $$N_{\text{co}1}\left({\beta }_{\text{T}}\right)=N_{\text{co}1}\left(0\right){\cos }^{0.24}\left({\beta }_{\text{T}}\right)$$(20)

to scale the nominal results to different values of β_T.

But the change in β_T will result in an asymmetric dust distribution along the S/C trajectory. This is illustrated for S/C B2 in Fig. 15 showing the examples where β_T is 0° (top row) and 60° (bottom row). In the case where β_T = 60° the inbound densities are lower than for the symmetric trajectory because the inbound asymptotic phase angle is 150° instead of 90°. Conversely, the outbound densities are higher than for the symmetric trajectory because the outbound asymptotic phase angle is 30° instead of 90°. Further the maximum number density will occur shortly after CA instead of at CA.

Finally, we should state that all changes in the CA distance, α_T, and β_T result in small changes in the column densities compared to the inherent uncertainties resulting from our poor knowledge of the input parameters.

6.3 Scaling of results to other dust sizes

In the case that the model does not provide the results for the desired dust radius or size of the dust bin (by default half a decade) simple scalings can be applied. We should point out to use these scaling with caution as they might break down if pushed to their extremes (e.g. to very short size intervals). All scaling we provide here apply only to the respective median cases of the column density of the respective spacecraft.

First, we define some general functions to calculate the respective scaling later on. The column density as a function of the dust radius, a, follows a double asymptotic power law of the form $$N_{\text{co}1,\text{fit}}\left(a\right)=C_1{a}^{-{\beta }_1}+C_2{a}^{-{\beta }_2},$$(21)

where β_i are the respective asymptotic power law exponents, and $$C_i\left(a_{0,i},N_{0,i},{\beta }_i\right)=N_{0,i}{a}_{0,i}^{{\beta }_i}$$(22)

is the normalisation constant given where My is the column density at reference size a_0,i. The values of these parameters for the three spacecraft are given in Table 2.

Next we define the length of our dust size bins as $$\text{Δ}a={10}^d,$$(23)

where d is a decade fraction, for example, 0.5 for half a decade or 0.25 for a quarter decade. From this it follow that a specific dust bin starting at a will extend to a · Δa.

To interpolate the results from one dust bin radius to another while keeping the length of the dust bin (Δa) unchanged we can proceed in the following manner. When the target dust size is a_t we lookup the dust bin size, a_r, which is just smaller than our target size. This reference dust bin, a_r, has a certain column density, N_col(a_r), according to the EDCM. Using Eq. (21) the column density for our target size, a_t (covering the dust range from a_t to a_t · Δa), will be $$N_{\text{co}1}\left(a_{\text{t}}\right)=N_{\text{co}1}\left(a_{\text{r}}\right)\frac{N_{\text{co}1,\text{fit}}\left(a_{\text{t}}\right)}{N_{\text{co}1,\text{fit}}\left(a_{\text{r}}\right)}.$$(24)

Fig. 14 Dependence of the column density with αT and βT. αT is the angle between the orbital plane of the comet and the direction from the comet to the spacecraft at closet approach (α in Fig. 4). βT is the angle between the direction of the Sun and the direction from the nucleus to the closest approach projected into the orbital plane (β in Fig. 4 and not to be confused with β the power law). Left panel: normalised column density as a function of αT (°). Right panel: normalised column density as a function of βT (°)

Fig. 15 Number densities as a function of cometo centric distance for S/C B2 for different dust sizes. Top left: case where αT = βT = 0°. Top right: same as top left but zoomed view. Bottom left: case where αT = 0°;βT = 60°. Bottom right: same as bottom left but zoomed view.

6.4 Interpolation of dust bin length

Next we discuss how the length of the dust bins can be changed. The results of the EDCM are provided for bins with a length of half a decade in dust radius. If one wants to know the results for a bin size of a quarter decade instead of a half a decade the following procedure can be applied. We should caution the reader at this point that we have tested this scaling only down to a quarter decade. While we believe that it also holds for smaller intervals the user should be aware of this limitation and that smaller bin lengths could be associated with additional uncertainty.

The provided solution have dust bin lengths which are half a decade in dust radius. Or according to Eq. (23) this means that our reference decade fraction is d_r = 0.5. The target decade fraction shall be d_t. First we define the helper function $$H\left(a_0,N_0,\beta ,d,a\right)=\frac{C\left(a_0,N_0,\beta \right)}{\left(1-\beta \right)}{a}^{1-\beta }\left({10}^{d\left(1-\beta \right)}-1\right),$$(25)

which allows us to calculate the integral, I, of Eq. (21) for one dust bin: $$\begin{array}{lll}I\left(a_{0,1},N_{0,1},{\beta }_1,a_{0,2},N_{0,2},{\beta }_2,d,a\right)\hfill & =\hfill & H\left(a_{0,1},N_{0,1},{\beta }_1,d,a\right)\hfill \\ \hfill & \hfill & +H\left(a_{0,2},N_{0,2},{\beta }_2,d,a\right).\hfill \end{array}$$(26)

From this, we can now calculate the scaling constant, S, for a bin that was initially covered the decade fraction d_r to one that covers the decade fraction d_t. This scaling constant can be written as $$\begin{array}{l}S\left(a_{0,1},N_{0,1},{\beta }_1,a_{0,2},N_{0,2},{\beta }_2,d_{\text{r}},d_{\text{t}},a\right)\\ =\frac{I\left(a_{0,1},N_{0,1},{\beta }_1,a_{0,2},N_{0,2},{\beta }_2,d_{\text{t}},a\right)}{I\left(a_{0,1},N_{0,1},{\beta }_1,a_{0,2},N_{0,2},{\beta }_2,d_{\text{r}},a\right)}\end{array}$$(27)

The column density of the new bin length, N_col(a, d_t), can now be calculated from the column density given by the EDCM, N_col(a, d_r): $$\begin{array}{lll}{N}_{\text{co}1}\left(a,d_{\text{t}}\right)\hfill & =\hfill & N_{\text{co}1}\left(a,d_{\text{r}}\right)\hfill \\ \hfill & \hfill & \times S\left(a_{0,1},N_{0,1},{\beta }_1,a_{0,2},N_{0,2},{\beta }_2,d_{\text{r}},d_{\text{t}},a\right).\hfill \end{array}$$(28)

Though we have kept all equations general, we remind the reader that currently d_r = 0.5.

7 Discussion and conclusions

In this work we have considered a new approach to determine the dust environment of a currently unknown comet as the one to be intercepted by ESA’s CI mission. We have argued that it is not clear – a priori – how a best-, nominal-, or worst-case scenario can be defined because of the non-linear interdependence of the parameters in any dust coma model.

We have therefore taken the opposite approach of defining ranges for all input parameters to the dust coma model and identifying all possible parameter combinations. The collection of all these combinations provides a set of coma environments which can then be analysed statistically to determine a median case, for instance.

The results we have presented show that the range of possible dust densities, which can be encountered by a spacecraft, span many orders of magnitude. This large uncertainty is dominated by the unknowns of the parameters of the nucleus, the dust particles, and the coma. As long as the uncertainties in these parameters remain large, so does the uncertainty of the predicted dust environment. However, as we learn more about a target, these uncertainties can be successively diminished.

We should also point out that even though it is possible to define a clear median case, this obviously does not imply that it is the most probable case. In this work we have made no attempt to quantify which input parameters are more likely than others. The EDCM assumes that each parameter is equally likely. Thus the model is agnostic as to the likelihood of any parameter taking a specific value and assigns it an equal probability. We should therefore not expect the median, but rather any value within a certain confidence interval around it.

Our model shows that the brightness Afρ (as defined in A’Hearn et al. 1984, where A is the bond albedo, f is the filling factor of grains within the field of view, and ρ is the radius of the field of view) alone is a poor indicator of the dust production rate. In combination with a known dust size distribution, the expected uncertainty can be reduced to ~1 order of magnitude.

The parameters most predictive of the column density along the spacecraft flyby are a combination of Afρ, the power law exponent, and the dust production rate. For small particles, Afρ in combination with either the power law exponent or the dust production rate is the most predictive. In contrast, for large particles, the power law exponent in combination with either Afρ or the dust production rate is the most predictive.

Unfortunately, only Afρ (or an equivalently measured brightness in flux or magnitudes) is an observable. The size distributions and dust production rates are always model dependent. Any assessment of the coma environment, therefore, needs to be rephrased to what range of model parameters are possible to match a given Afρ.

A somewhat independent constraint on the size distribution can be achieved using Finson-Probstein modelling of the dust tails (e.g. Agarwal et al. 2010). However, even this size distribution is a model-dependent interpretation rather than a direct observational measurement. Thus, because only Afρ is available as a direct measurement of the three most predictive parameters of the dust abundance in the coma, the only remaining approach is the one we have outlined in this work. For a certain Afρ, one can assess what comet or coma properties are consistent with that measurement.

For typical encounter velocities, we have shown that the angle of impact of dust particles onto the spacecraft are dominated by the spacecraft motion. This means that shielding of the spacecraft can be concentrated to the ram side of the spacecraft. Though there might still be implications for attitude control which need to be checked separately.

Finally, we provide all results to the public in easy-to-read formats³. Should the standard results we present here not satisfy a reader’s needs, we have provided some scaling relationships in Sect. 6 to adapt our results.

Acknowledgements

We acknowledge the support from the International Space Science Institute (ISSI) through the team “Closing The Gap Between Ground Based And In-Situ Observations Of Cometary Dust Activity: Investigating Comet 67P To Gain A Deeper Understanding Of Other Comets”. R.M. acknowledges the support from the Swiss National Science Foundation (SNSF) under the grant P2BEP2_184482, and funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement no. 101019380). VDC acknowledges the support by the Italian Space Agency (ASI) within the ASI-INAF agreements I/032/05/0, I/024/12/0 and 2020-4-HH.0.

References

All Tables

All Figures

	Fig. 1 Sketch of the components of the EDCM and which input parameters feed into which step.
In the text

	Fig. 2 $\text{Log}\left(\text{Af}\rho /{\text{Q}}_{{\text{H}}_2\text{O}}\right)$ from A’Hearn et al. (1995) for typical (yellow) and depleted (orange) comets as well as the range of the EdCM.
In the text

	Fig. 3 Phase curves of 1P (Schleicher) and 67P (Bertini et al. 2017) are shown.
In the text

	Fig. 4 Giotto flyby trajectory relative to the cometary nucleus and Sun is shown, including offset angles and positions forthelast4usefulHMC images (from Keller et al. 1987).
In the text

	Fig. 5 Total number of particles during encounter as a function of mass for the EDCM Halley case and the data by McDonnell et al. (1997).
In the text

	Fig. 6 Number of dust particles according to the EDCM along the spacecraft trajectory of spacecraft A as function of cometo-centric distance. The shaded areas show different percentile ranges within which cases fall. Closest approach for S/C A is assumed at 1000 km.
In the text

	Fig. 7 Total number of dust particles encountered according to the EDCM along the spacecraft trajectory of spacecraft A as function dust mass. The shaded areas show different percentile ranges within which cases fall. Additionally the orange curve shows the predicted median densities for a Halley type comet (according to the definition in Sect. 3).
In the text

	Fig. 8 Dust column density for dust particles with radii between 1 and 3 μm at the Comet Interceptor spacecraft A is shown as a function of different combinations of Afρ, the power law exponent, β, and the dust production rate.
In the text

	Fig. 9 Dust column density for dust particles with radii between 10 and 30 cm at the Comet Interceptor spacecraft A is shown as a function of different combinations of Afρ, the power law exponent, β, and the dust production rate.
In the text

	Fig. 10 Dust production rate as a function of Afρ and the power law exponent, β, of the dust size distribution. The values of Afρ are slightly offset depending on the power law exponent to show the overlap resulting from β. By itself Afρ is a poor predictor of the dust production rate of a comet.
In the text

	Fig. 11 Dust column density for dust particles with radii between 100 and 300 μm at the Comet Interceptor spacecraft A is shown as a function of different combinations of the gas production rate, dust bulk density, surface activity distribution, night time activity, and radiation pressure efficiency. These quantities are not predictive of the dust column density.
In the text

Fig. 12 Dust speed and impact angles for different dust sizes and activity levels. Left panel: terminal dust speed, υd,term, according to Eq. (15) for two gas cases. The lowest gas production rate case considered in the EDCM in purple (min Qg case, where RN = 2 km and Qg = 155 kg s−1) and the highest gas production rate case in orange (max Qg case, where RN = 5 km and Qg = 16 434 kg s−1). Right panel: maximum angle of impact of dust particles with respect to the ram-direction of the spacecraft. The solid lines are for a fly-by speed of 10 km s−1 and the dashed lines for 70 km s−1. The colours are the min. and max. gas production cases as in the left panel. The corresponding radii of particles with masses 1 mg, 10 mg, and 100 mg are indicated by the grey lines.

In the text

	Fig. 13 Column density multiplied with the CA distance as a function of dust mass.
In the text

Fig. 14 Dependence of the column density with αT and βT. αT is the angle between the orbital plane of the comet and the direction from the comet to the spacecraft at closet approach (α in Fig. 4). βT is the angle between the direction of the Sun and the direction from the nucleus to the closest approach projected into the orbital plane (β in Fig. 4 and not to be confused with β the power law). Left panel: normalised column density as a function of αT (°). Right panel: normalised column density as a function of βT (°)

In the text

	Fig. 15 Number densities as a function of cometo centric distance for S/C B2 for different dust sizes. Top left: case where αT = βT = 0°. Top right: same as top left but zoomed view. Bottom left: case where αT = 0°;βT = 60°. Bottom right: same as bottom left but zoomed view.
In the text