© The Authors 2024

1 Introduction

One of the primary scientific objectives of the ESA Rosetta mission was to explore and understand the onset, evolution, and decline of activity exhibited by comet 67P/Churyumov-Gerasimenko (hereafter referred to as 67P) in situ. To achieve this goal, the spacecraft was equipped with a suite of complementary instruments. These instruments employed various techniques and targeted different size ranges of particles, allowing for the analysis and collection of data regarding the refractory material found within the coma of 67P. A synthesis of the main findings for all size ranges can be found in Güttler et al. (2019).

While most of the dust analysis instruments on board Rosetta were designed for particles with typical sizes ≲1 mm, the mission offered a unique opportunity to observe and study larger particles, exceeding this size range, using the Optical, Spec-troscopic, and Infrared Remote Imaging System (OSIRIS), the scientific imaging instrument (Keller et al. 2007). The main challenge of using these data comes from the lack of information about their distances and velocities along the line of sight (LOS) in the images. Several approaches have been taken to address this issue by different authors.

Fulle et al. (2016a) used the parallax method, assuming that particles moved radially from the nucleus. Any deviations from this radial path were attributed to spacecraft motion, enabling them to estimate distances. Several authors (Agarwal et al. 2016; Pfeifer et al. 2022, 2024; Shi et al. 2024) focused on particles appearing to emerge from the comet limb, treating them as being at the same distance as the nucleus. Güttler et al. (2017) studied particles appearing out of focus and correlated their blurring levels with distances. Drolshagen et al. (2017) and Ott et al. (2017) examined particles detected in both the Wide and Narrow Angle Cameras (WAC and NAC), using the differences in measured position between these cameras to calculate distances. In a different approach, Frattin et al. (2021) utilized images where particles can be seen as bright tracks, estimating distances to the camera based on the length of these tracks and making assumptions about typical particle speeds.

In this study we employed a method introduced in Lemos et al. (2023). This method offers an indirect approach based on the comparison between OSIRIS images and synthetic counterparts generated by varying particle properties. The goal is to identify the particle properties that produce the closest matches to the observed images. To achieve this, we analyzed OSIRIS images in which particles are visible as tracks due to their relative motion with respect to the spacecraft. The properties of these tracks, namely their orientation angle, length, and integrated brightness, were extracted from the images.

In parallel, we conducted dynamical simulations to trace the trajectories of particles. Based on the position of the spacecraft and the camera viewing direction, we generated synthetic images for the times at which the OSIRIS data were obtained. Then, the properties of the tracks within the synthetic images were inferred, and their distributions compared to their real counterparts. The simulated particles that generate the tracks that best match the observed data are considered to most closely resemble the real particle population.

The primary objective of this study is to conduct an in-depth analysis ofvarious OSIRIS images employing this method. This approach provides us valuable insights into the physical characteristics of the aggregates responsible for generating the observed tracks. It also sheds light on the source regions of these aggregates on the surface of the comet and offers an estimation of their flux. This information plays a pivotal role in constraining the activity mechanisms driving aggregate ejections.

In Sect. 2, we provide a comprehensive overview of the OSIRIS images utilized in our analysis, along with a description of the track detection algorithm. Section 3 examines the dynamical simulations employed to create synthetic images, which are subsequently compared with the OSIRIS images. Section 4 presents the results derived from this comparison, offering insights into aggregate properties, source regions, and flux. Finally, Sect. 5 summarizes the key findings of this work.

2 Observations

2.1 Image selection

For this work we used images obtained with the OSIRIS/NAC. For more details about the characteristics of the camera we refer the reader to Keller et al. (2007). The images were acquired in five different sets between July 2015 (r_h = 1.32 au inbound) and February 2016 (r_h = 2.39 au outbound) with the number of images per set ranging from 21 to 45 (Table 1). These images were acquired to measure the scattering phase function of the dusty coma, so they share some basic characteristics. The spacecraft was located over the terminator in such a way that the nucleus-Sun vector was roughly perpendicular to the nucleus-spacecraft vector. The observations sampled different phase angles α on a plane that was perpendicular to the nucleus-spacecraft vector, and for each phase angle three images were taken using different filters. In all the sets used in this work, the images were obtained using the broadband Blue F24 (480.7 nm), Orange F22 (649.2 nm), and Red F28 (743.7 nm) filters. Figure 1 shows a sketch of the observation geometry for all sets used in this work.

These images show a more or less homogeneous background, resulting from the light scattered by small dust grains, as well as bright, straight tracks resulting from the projection onto the image plane of the trajectories of individual aggregates with larger typical sizes on the image plane. While the background brightness in such sequences was studied by, for example, Bertini et al. (2017) and Keiser et al. (2024) to constrain the spatial distribution of scattering cross-section and the dust scattering properties, we concentrated on the individual tracks in the foreground, in a dataset overlapping with that used by Bertini et al. (2017).

We analyzed a dataset comprising 189 images. These images were subject to processing, attaining a processing Level 3E according to the OSIRIS scale, equivalent to Level 4 on the Committee on Data Management and Computation (CODMAC) scale. This implies that the images are solar stray light corrected, in-field stray light corrected, radiometrically calibrated and geometric distortion corrected, expressed in radiance units (Tubiana et al. 2015).

Fig. 1 Sketch of the observation geometry for all the image sets used in this work. The dashed blue line indicates the solar direction. The dotted red line indicates the position of the spacecraft. The solid violet lines indicate the different phase angles sampled by the camera when images with the three different filters were obtained. All these sampled directions are located in the green plane, perpendicular to the nucleus-spacecraft direction. Note that the sizes are not to scale.

Fig. 2 Diagram of the main steps involved in the detection algorithm.

2.2 Detection method

We are interested in analyzing the properties of the aggregates with typical sizes ≳1 cm. As mentioned in Sect. 1, these aggregates are seen as bright tracks in the studied images. We used a semiautomatic method that exploits this property to detect these particles. This method, introduced by Lemos et al. (2023) and based on that presented by Frattin et al. (2017), consists of four automatic steps, plus a last manual one:

1. Perform the normalized cross-correlation of the images with track templates in order to generate similarity maps (i.e., a representation of regions in the image that have high probability of containing a track).

2. Create binary images from the similarity maps.

3. Detect tracks on the binary images using the Hough transform method (Hough 1962; Duda & Hart 1972).

4. Refine the results in order to correct imperfections of the algorithm.

5. Manually inspect the results.

A diagram of the automatic steps is presented in Fig. 2.

By applying this procedure to the groups of images described before, we detected 34 616 tracks. The tracks are characterized by three parameters: the orientation angle between the track and the vertical of the image, the total length measured in pixels and the integrated brightness, which was determined using the technique outlined by Güttler et al. (2019). This method involves performing an aperture photometry analysis for each track, but using a stadium-shaped, rather than a circular aperture.

3 Synthetic images

3.1 Dynamical simulations

For an accurate description of the aggregate trajectories, a detailed model of the gas flow in the coma is needed. For this purpose, we used the gas flow results from Marschall et al. (2020), obtained using a direct simulation Monte Carlo (DSMC) method. This code computes the energy input in each facet constituting the nucleus shape model, calculates the surface temperature resulting from the energy balance of incoming light, thermal emission and sublimation, and the gas production rate from each facet, and subsequently computes the steady-state gas distribution in the coma using a DSMC code. In this way, the gas field used for the whole duration of the simulation is the one corresponding to the ejection instant.

The shape model used is a decimated version of the stereo-photogrammetric model SHAP7 (Preusker et al. 2015) formed by ~440 000 facets with typical length-scales of ≃10 m. For more details about this method we refer to Marschall et al. (2016, 2020).

We used a modified version of the DRAG3D code (Marschall et al. 2016, 2020) to solve the dust aggregates equations of motion. This code uses a fourth-order Runge-Kutta solver with a variable time step to solve the equations describing the motion of spherical dust aggregates in the cometary coma.

The original code includes the effects of the gravitational and gas drag forces, which we have extended by the solar radiation pressure and solar tidal forces. The gravitational field is computed by discretizing the volume of the mentioned shape model into volume elements with typical sizes ≃30 m and summing the contributions of each element, assuming an homogeneous density and a total mass M_N = 9.9 × 10¹² kg (Sierks et al. 2015).

The gas drag force F_D is calculated using the equation $$F_{\text{D}}=\frac{1}{2}{C}_{\text{D}}{m}_{\text{g}}{n}_{\text{g}}\sigma \left|v_{\text{g}}-v_{\text{d}}\right|\left(v_{\text{g}}-v_{\text{d}}\right),$$(1)

where m_g, n_g and v_g are the molecular mass, number density and velocity of the gas, respectively, v_d is the dust velocity and C_D is the drag coefficient. Assuming the gas flow is in equilibrium, C_D is defined as $$C_{\text{D}}=\frac{2s^2+1}{\sqrt{\pi }{s}^3}{\text{e}}^{-s^2}+\frac{4s^4+4s^2-1}{2s^4}\mathrm{erf}(s)+\frac{2\sqrt{\pi }}{3s}\sqrt{\frac{T_{\text{d}}}{T_{\text{g}}}},$$(2)

where T_g and T_d are the gas and dust temperatures, respectively (assumed to be equal), and $$s=\frac{\left|v_{\text{g}}-v_{\text{d}}\right|}{\sqrt{\frac{2k_{\text{B}}{T}_{\text{g}}}{m_{\text{g}}}}}.$$(3)

The solar radiation pressure, F_R, and solar tidal forces, F_T, are described as follows: $$F_{\text{R}}=-\frac{C_{\odot }{Q}_{\text{RP}}\pi r_{\text{d}}^2}{r_{\text{h}}^{3}c}{r}_{\text{h}}$$(4) $$F_{\text{T}}=F_{\text{S}}\left(r_{\text{C}}\right)-F_{\text{S}}\left(r_{\text{C}}+r_{\text{D}}\right),$$(5)

where C_⊙ = 1361 W m⁻² is the solar constant, Q_RP is the radiation pressure efficiency (equal to 1), c is the speed of light, F_S(r) is the solar gravity force at r, r_C is the Sun-comet vector, and r_D is the nucleus-dust aggregate vector.

Since the equations of motion are solved in the rotating frame, centrifugal and Coriolis forces are also included in order to account for nucleus rotation. The latter are particularly important for the dust aggregates studied in this work because they have low speed, which implies timescales to traverse the nucleus-spacecraft distance on the order of or much longer than the rotation period of the nucleus. For this work we assumed that the dust aggregates do not sublimate volatiles from their surface. For each simulation, a minimum of one and a maximum of ten aggregates were located in a random position of each facet, with the number of particles per facet scaling linearly with the total water production rate (production rate per unit area times area of the facet). The aggregates were given an initial velocity with a direction aligned with the facet normal and a modulus chosen randomly from a Maxwell-Boltzmann (MB) distribution. The MB distribution is parametrized using its most probable velocity, υ_P. With this, the aggregates were characterized by three parameters: particle radius (r_d), particle density (ρ_d), and most probable initial velocity (υ_P). The values used for these parameters are shown in Table 2. The simulations using different combinations of these parameters were carried out independently, and the solutions were later combined in a procedure explained in Sect. 3.2.

Since the dust flight time needed to reach the camera field of view (FOV) was a priori unknown, we repeated the dynamic simulations using gas solutions for different sub–solar longitudes spanning a whole comet day in steps of 30°. This means that for each simulated epoch, simulations were performed for 12 nucleus rotation states × 60 aggregate parameters combinations =720 cases. The total number of aggregates simulated in each case is variable, depending on the distribution of water vapor production on the surface, but is always on the order of 5 × 10⁵ particles.

3.2 Generation of synthetic images

To generate synthetic images from which a synthetic track population can be extracted, understanding the motion of aggregates in relation to the spacecraft is crucial. Initially, the position of the spacecraft as well as its observation direction were computed using the SPICE system (Acton 1996; Acton et al. 2018). The camera’s FOV is described as a pyramid with the spacecraft at its apex. Consequently, aggregates observable by the camera are those intersecting this shape. However, due to the characteristics of our simulations such as the high number of particles and the discrete rotation states of the nucleus, determining these intersections is not a straightforward task. Thus, a customized method was devised to address this issue.

At large distances, the gas and dust velocities relative to the nucleus permanently increase (Zakharov et al. 2018b). However, it has been shown that the changes of dust velocity above the acceleration region (i.e., the region where the gas efficiently accelerates the dust) on the order of 10 km above the nucleus surface depending on the particle size (Gerig et al. 2018; Zakharov et al. 2018b), are very small. For this reason and to simplify the check for intersections between trajectories and the FOV, we assumed that the trajectories of dust aggregates within the range sampled by the camera are straight lines with constant velocity in a nonrotating frame centered on the nucleus. To address potential influences from nucleus gravity and radiation pressure, an additional examination of the trajectories was conducted. The results consistently affirmed that these trajectories indeed followed straight lines, reinforcing the validity of our approximation.

To ensure a continuous coverage of the entire range of rotation angles, each simulation was deemed valid within the range of [−15, +15]° from its nominal rotation angle. The trajectories were parametrized with respect to an angle θ_I that belonged to that range. Using this parametrization, the minimum distance between the particle trajectory and the LOS was determined as a function of θ_I. Should the angle between the vector connecting the spacecraft to the aforementioned point of minimal distance and the LOS direction be smaller than the angular dimension of the FOV (2.20°), the aggregate in question was said to intersect with the FOV and was then selected for further analysis.

Defining the sides of the FOV as four triangular shapes, we determined the range of angles θ_i ∈ θ_I that provide intersections with each face individually using the algorithm presented by Möller & Trumbore (1997). This intersection provides two points, named entry and exit points, respectively, which define the direction of the track in the synthetic image. The position of the track endpoints on the synthetic image were associated with the exposure time of the real image. To define those endpoints, we selected a random position between entry and exit points. Then, the exposure time of the image was divided into two intervals (not necessarily equal), and the position of this random point was propagated both forward and backward along the trajectory, defining both track endpoints. As already mentioned, the region observed by the FOV is far away from the aggregate acceleration zone, so it is safe to assume that the aggregate velocity remains constant along their trajectories through the FOV. If an endpoint is not located within the FOV, it is replaced by the corresponding entry or exit point. Lastly, these points were projected onto the image plane, defining the synthetic track.

However, an extra condition must be met by the trajectories in order to accept the synthetic tracks: the time at which particles intersect the FOV must fall within the time range when the corresponding image was obtained. Quantitatively, this is the same as saying that the subsolar longitude at the moment of observation that is recorded on the image header λ_obs must be compatible with that of the simulations, obtained by adding the subsolar longitude at the moment of ejection λ_ej to the nucleus rotation during the flight time, t_F, defined as the time elapsed from the aggregates ejection to the moment they intersect the FOV. This condition reads as $${\lambda }_{\text{obs}}={\lambda }_{\text{e}j}+{\omega }_{\text{R}}\times t_{\text{F}},$$(6)

where ω_R is the nucleus angular velocity of rotation.

As mentioned before, λe j does not necessarily match with the subsolar longitude of the gas solutions λ_gas, since the angle θ_i was added to account for discontinuities. Plugging this into Eq. (6) we obtain $${\lambda }_{\text{obs}}={\lambda }_{\text{gas}}+{\theta }_i+{\omega }_{\text{R}}\times t_{\text{F}}$$(7) $${\theta }_i={\lambda }_{\text{obs}}-{\lambda }_{\text{gas}}-{\omega }_{\text{R}}\times t_{\text{F}}.$$(8)

For a trajectory to be accepted, the angle θ_i found from Eq. (8) has to belong to the range defined by the angles θ_I that provide intersections. For the accepted trajectories we generated synthetic images from the projected tracks, and stored their orientation angle, length, and integrated reflectance, B. The integrated reflectance was found using the equation from Agarwal et al. (2016), $$B=\frac{p\text{Φ}(\alpha )I_{\odot }{r}_{\text{d}}^2}{r_{\text{h}}^2{\text{Δ}}^2},$$(9)

where p is the aggregate geometric albedo, Φ(α) its phase function normalized to 1 at zero phase, I_⊙ the solar flux in the corresponding filter, r_d is the dust aggregate radius, r_h the heliocentric distance in au, and Δ is the distance between the aggregate and the camera. If the track was not completely contained inside the synthetic image, we scaled the integrated reflectance B with the fraction of the track belonging to the image. Lemos et al. (2023) found that the phase function of the dust aggregates follow the equation: Φ(α) = exp (−β × α), with a mean β = 8.2 × 10⁻³. We used p = 0.065 (Agarwal et al. 2016) and defined Δ as the mean of the distances between the camera and entry and exit points, respectively.

3.3 Parameter optimization via comparison between observed and synthetic tracks

This comparison was based on the three track properties mentioned in Sect. 2.2: orientation angle, length and integrated brightness. However, since the number of tracks detected in different OSIRIS images was highly variable, we generated an estimation of the probability density function for the distribution of properties from observed tracks using a Gaussian kernel density estimation. Based on this probability density function estimation, we created a new sample of track property triplets that represent the observed ones. Despite the introduction of new statistical uncertainties, this method was preferred since it makes sure that the new set of track properties that will be fitted has a large number of elements for all OSIRIS images studied.

At this point, two groups of tracks exist: the ones representing the observed tracks, and the synthetic ones obtained from the dynamical simulations. For each element in the group that represents the observed tracks, we searched the synthetic track pool for the five closest elements in the orientation angle-length-brightness space. This subset of closest neighbors is selected as the result of the comparison procedure.

4 Results and discussion

With the procedure presented in previous sections, we found the properties of the simulated aggregates that provide the best fits to the tracks observed in OSIRIS images. In this section we present these properties, as well as the regions on the surface from where these aggregates originate and the mass flux escaping from the comet. Additionally, we discuss the implications of these results for ejection and activity models.

4.1 Effects of gas on particle dynamics

While simulations without gas were not conducted, we can compare results across sets obtained at different epochs to gain insights into the influence of gas on particle dynamics. Figure 3 illustrates the proportion of aggregates escaping the integration domain, categorized based on their aggregate parameters, for the image sets STP070 (7 days after perihelion) and STP096 (189 days after perihelion).

In the case of aggregates for which the gas drag forces are highly efficient (small radii and densities), there is a noticeable decrease in the fraction of escaping particles for STP096, which was obtained far from perihelion. In contrast, the fraction of larger, denser aggregates that escape remains relatively constant, indicating a diminished impact of gas on their dynamics.

As expected, the initial velocities also significantly influence the fraction of escaping aggregates. Near perihelion, about 40% of aggregates with zero initial velocity escape, indicating efficient acceleration by gas drag. Conversely, this number decreases to approximately 5% for set STP096, highlighting the diminished influence of gas at that point in the orbit.

This change in efficiency highlights the changing nature of the interactions between gas and aggregates phenomenon occurs because we are observing a region where the impact of gas drag shifts from high to low. This influence depends not only on the aggregate properties but also on the characteristics of the gas flow.

4.2 Aggregate properties

We defined two types of simulated aggregates. First, all those aggregates that generate tracks meeting the requirements explained in Sect. 3.2 (intersecting the FOV and doing so at a time compatible with the observation) are denominated candidates. The subgroup of candidates that generate tracks that are selected as results of the comparison procedure explained in Sect. 3.3 are called fitted aggregates. To assess the relevance of each physical property of the aggregates for replicating the observed tracks, we generated histograms for radius, density and most probable initial velocity of candidates and fitted aggregates, considering all the images in each analyzed set. Figure 4 shows an example of this type of histograms, in this case for the results of all images from STP092. The histograms for the remaining image sets can be found in Appendix A.

The data in Fig. 4 provide insights into the distribution of properties for both candidates and fitted aggregates, with the last column displaying the ratio between the two histograms. The information that can be extracted from them is different. In this particular example, the distribution of candidate radii seems fairly uniform, suggesting that the radius may not be a determining factor for an aggregate to reach the FOV. However, the ratio column reveals a distinct peak around 5 cm, indicating that aggregates of this size are the ones producing tracks similar to the observed ones. In terms of initial velocities, both candidates and fitted particles exhibit similar distributions, emphasizing the importance of this parameter for an aggregate to reach the FOV. Nevertheless, the ratio distribution is approximately uniform, suggesting limited influence on the type of track generated by the aggregates.

Figures 5 and 6 show the mean value for each physical property for all image sets as a function of phase angle, while Table 3 shows their mean considering all images in each set, averaged across all phase angles¹. Two key observations emerge from these figures. Firstly, the sets exhibiting greater parameter dispersion are the ones in which the spacecraft was positioned farther from the nucleus. This diminishes the total count of simulated aggregates that reach the FOV, which, in turn, affects the quality of the results, especially at large phase angles. Secondly, the results are similar to those presented in Lemos et al. (2023), but with a clear difference in particle size. In our case, the mean aggregate sizes lie in the centimeters range, while their results showed aggregates of several decimeters. The source of this difference is a combination of dynamical effects: on the one hand, the model applied here considers additional forces, not taken into account in the mentioned work. On the other hand, the model used for the gas dynamics is different. Lemos et al. (2023) used a fluid model for simulating the gas distribution, while here a DSMC model was used. Zakharov et al. (2018a) showed that the fluid approach yields markedly higher velocity estimates near the nucleus compared to DSMC simulations. Given that these regions exhibit the highest gas densities, the most significant gas drag acceleration occurs there, resulting in significantly greater gas drag expected from the fluid model relative to DSMC. Consequently, according to the predictions of Lemos et al. (2023), dust aggregate velocities are anticipated to be higher than those obtained here for a given size. To reconcile the discrepancies observed in OSIRIS images, Lemos et al. (2023) adjust these velocities by increasing the sizes of the aggregates, which explains the differences between inferred particle sizes among both approaches.

Although smaller than those found in Lemos et al. (2023), the aggregate sizes found in this work are a good match to those found by other authors from the analysis of OSIRIS images at a similar range of distances (Agarwal et al. 2016; Ott et al. 2017; Pfeifer et al. 2022, 2024). These aggregates have sizes larger than the ones expected from activity caused by water vapor sublimation, so more volatile species as CO₂ seem to be responsible for their ejection (Gundlach et al. 2020).

Our results also indicate that the aggregates have a density lower than that measured for dust particles using other instruments as GIADA, ρ_D ≃ 800 kg m⁻³ (Fulle et al. 2016a, 2017), but comparable with the bulk nucleus density ρ_N = 537.8 kg m⁻³ (Preusker et al. 2017). This can be interpreted as the aggregates seen by OSIRIS having the same macro-structure as the nucleus rather than that of pebbles.

Fig. 3 Fraction of aggregates escaping the integration domain for image sets STP070 (top) and STP096 (bottom). Each black dot represents the mean of all particles in a single simulation, whereas the red symbols denote the mean across all aggregates of that particular kind, considering the cumulative data from all simulations.

Fig. 4 Example of the distribution of physical properties for aggregates in STP092 obtained using our method. From top to bottom, rows show the distribution of radius, density, and most probable velocity. The first and second columns show candidates and fitted particles, respectively. The last column shows the ratio between heights of columns for fitted aggregates and candidates.

Fig. 5 Mean aggregate parameters for sets STP063 (top), STP070 (middle), and STP086 (bottom) as a function of phase angle. Each row represents a different image set, while the columns from left to right represent the aggregate radius (rd), density (ρd), and initial velocity (υin). The colors show the results for images taken with filters Blue F24, Orange F22, and Red F28.

4.3 Initial velocities

Our method shows that to reach the FOV, the aggregates require an initial velocity υ_in on the order of 1 m s⁻¹. This comes from the fact that, due to their large size, the aggregates are only weakly affected by gas drag, so their initial energy is the main parameter defining whether they are able to reach altitudes sampled by OSIRIS. Most of the work in the field of dust ejection assumes that particles that are initially resting on the surface are lifted by some mechanism, principally gas drag. However, several studies focusing on laboratory measurements, observational data and computational modeling suggest the potential inclusion of an initial velocity component for dust particles. Yelle et al. (2004) and Huebner et al. (2006) showed that if that dust is ejected from a porous, subsurface layer of the nucleus, the acceleration provided by gas drag before reaching the surface could be interpreted as an initial velocity, and applied this idea to explain the jets observed in comet 19P/Borrely by the Deep Space 1 spacecraft. Continuing along this line of thought, Kramer & Noack (2015, 2016) employed a modeling approach to investigate dust transport and the subsequent formation of structures within the inner coma of 67P assuming a nonzero initial velocity. Examining the coma of comet C/2017 K2, Kwon et al. (2023) conducted an analysis using images obtained by the ground-based Very Large Telescope (VLT). Their findings indicate that the observed structures cannot be adequately accounted for by any of the tested combinations of dust parameters unless an initial velocity is introduced. Laboratory experiments carried out by Bischoff et al. (2019) monitored the trajectories of dust aggregates ejected from a refractory layer situated atop a sublimating water ice block. Using a parabolic function fitting, these authors reported a nonzero initial velocity for the aggregates.

In this work, we investigated an analogous explanation to that proposed by Yelle et al. (2004) and Huebner et al. (2006), which has served as a base for other ejection models, principally dealing with jets (e.g., Belton 2010; A’Hearn et al. 2011; Vincent et al. 2016; Lin et al. 2017; Wesołowski et al. 2020). In this model, initial velocity is interpreted as a result of the ejection process itself. Figure 7 provides a schematic representation of this phenomenon, illustrating the scenario where the aggregate is bound to the nucleus, and the void space beneath it is filled with gas. When the gas pressure surpasses the tensile strength, the bond is broken, resulting in a net force on the aggregate. This force acts temporarily, ceasing when the pressurized gas aligns with the state of the surrounding environment.

In a rough approximation, we assumed that all the internal energy of the gas, U_gas, causing the ejection is transferred and converted into kinetic energy of the aggregate K_agg. It is important to note that in reality, this assumption does not hold true. Thus, this approximation provides an upper limit for the kinetic energy transferred to the aggregate. Assuming that the gas behaves ideally we have $$\begin{array}{cc}{K}_{\text{agg}}& =U_{\text{gas}}\\ \frac{1}{2}\frac{4}{3}\pi r_{\text{d}}^3{\rho }_{\text{d}}{\upsilon }_{\text{in}}^2& =C_{V}nT,\end{array}$$(10)

where C_V , n, and T are the isochoric molar heat capacity, number of moles and temperature of the gas, respectively. Since we assumed the gas is ideal, nT = PV/R, where P and V are gas pressure and volume, respectively, and R is the universal gas constant. With this, Eq. (10) becomes $${\upsilon }_{\text{in}}=\sqrt{\frac{3C_{V}PV}{2\pi r_{\text{d}}^{3}R{\rho }_{\text{d}}}}.$$(11)

At the moment of ejection gas pressure equals the material tensile strength, σ_T. While direct measurements of tensile strength within the specified size range are lacking, various measurements across different sizes and theoretical estimations are available. In this study, the focus on large aggregates emphasizes the importance of estimations related to the macro-structure of the comet, as opposed to smaller components like grains or pebbles. For this, we utilize the approximation ${\sigma }_{\text{T}}=\frac{G}{\sqrt{r_{\text{d}}}}$, where G = 100 Pa m^1/2 (Biele et al. 2022). Lastly, the heat capacity of an ideal gas can be expressed using the equipartition theorem as $C_V=\frac{f}{2}R$, where ƒ represents the number of degrees of freedom of the molecule. Inserting this into Eq. (11) yields $${\upsilon }_{\text{in}}=\sqrt{\frac{3fGV}{4\pi {\rho }_{\text{d}}{r}_{\text{d}}^{7/2}}}.$$(12)

Since thermophysical models suggest that CO₂ is the most plausible candidate for explaining the ejection of chunks with these sizes. If the vibrational modes are not taken into account, ƒ is equal to 5 for carbon dioxide. Figure 8 shows the initial velocity as a function of the volume of the pore containing the gas that causes the ejection for aggregates with density ρ_d =500 kg m⁻³ and radii r_d =5 and 10 cm. For a 5 cm aggregate, the volume needed to provide a 1 m s⁻¹ initial velocity is 1.2 × 10⁻⁴ m⁻³, while for a 10 cm aggregate this volume increases to 1.3 × 10⁻³ m⁻³. The minimal porosity that provides the necessary energy is expressed by the formula ϕ_min = V_void/(V_void + V_agg), yielding values of approximately 20% for both scenarios. The macro-porosity of the nucleus was determined to be 0.4 (Fulle et al. 2016a), indicating that the pores amid these types of aggregates should contain sufficient gas to generate their acceleration.

Fig. 6 Same as Fig. 5 but for image sets STPO92 (top) and STPO96 (bottom).

Fig. 7 Sketch of the breaking process. Sublimated gas occupies the void space below the aggregate. When pressure equals tensile strength, the aggregate is ejected and a net force coming from the gas pressure acts on it, producing the observed initial velocity.

Fig. 8 Initial velocity of aggregates as a function of the volume of the pore containing the gas, as determined by Eq. (12). The blue line corresponds to 5 cm radius aggregates, and the red line 10 cm radius aggregates. The horizontal dashed line indicates an initial velocity of 1 ms−1.

4.4 Source regions

One advantage of our model is its utilization of a complete shape model, which allows a precise identification of the origin locations of the fitted aggregates on the nucleus surface. However, since the flight times are comparable to the nucleus rotation period, minor uncertainties in the determined flight times could significantly impact the determination of the source region, limiting the precision of our method. Therefore, it is essential to note that when referring to source regions, we are indicating a general region rather than pinpointing the exact point from which the aggregates originate.

The trajectory and flight time of the simulated aggregates are influenced not only by the parameters examined in this work, such as gas distribution and velocity, but also by factors not taken into account like particle shape. Fulle et al. (2015) and Frattin et al. (2021) have reported the presence of nonspherical dust particles near 67P, and dynamic analyses have demonstrated that their shape can impact their response to gas drag (Ivanovski et al. 2017b,a; Moreno et al. 2022). As a result, the findings presented in this section are interpreted in a broad, qualitative manner.

By utilizing the initial positions of candidates and fitted particles used in the dynamic simulations, we connected these positions on the surface to specific facets of the SHAP7 shape model (Preusker et al. 2017), which consists of 125 000 facets. It is important to note that this shape model differs from the one employed in the dynamical simulations, as the SHAP7 model is considered a standard model in comparison. This allowed us to introduce the parameters C and F to represent the counts of candidates and fitted aggregates per facet normalized to their total counts, respectively. In this context, C and F can be interpreted as the probability of a candidate or fitted particle originating from a specific facet. It is important to note that both parameters are constrained within the range of 0 < F, C < 1. Combining these parameters, we define the ejection efficiency, T, as $$T=\frac{F-C}{C}.$$(13)

The ejection efficiency coefficient can vary within a range starting from −1. If the value of T is −1, it signifies that aggregates originating from a particular facet can reach the FOV (C > 0), but none of them were identified as valid outcomes of the inversion process (F = 0). In simpler terms, certain aggregates would be able to generate tracks that are visible in the images, but tracks with those properties were not actually observed. This suggests that the region on the nucleus surface linked with this facet is not effectively ejecting particles of the investigated type. Conversely, facets with T > 0 indicate an increase in the count of successfully matched tracks arising from particles originating in that specific surface region.

Figure 9 shows the combined ejection efficiency across all the sets. The colors indicate the value of T, while white zones show regions for which no candidates, and hence no fitted particles, could be found. For facets sampled in more than one image set, the mean value was used. The ejection efficiencies for each individual set can be found in Appendix B. We defined zones characterized by high and low ejection efficiency as clusters, and identified them in the map shown in Fig. 9. Below, we enumerate these clusters, utilizing the regions defined by El-Maarry et al. (2015, 2016) and the subregions defined by Thomas et al. (2018). Additionally, we incorporated the fundamental terrain type, as per the classification by Thomas et al. (2018), to which each region is attributed.

High ejection efficiency

1. Ash (dust) – Apis (consolidated) – Imhotep (smooth) – Khonsu (consolidated): This cluster represents the interface between three very different regions. It comprises subregion c of Imhotep up to its boundary with subregion a, characterized by a intermediate scale roughness and the presence of layers, thought to be exposed by material loss. In the north, Ash-i shows large-scale roughness and layering in cliffs, while Apis is rough in many scales, and similar to subregion b of Khonsu, one of the most complex regions with evidence of significant change, most probably due to activity (El-Maarry et al. 2017).

2. Bes (consolidated) – Imhotep (smooth) – Khonsu (consolidated): This region runs along the interface between Khonsu-c to the northeast, Bes-a to the south and Imhotep-b to the west. This is a rough area, with several boulders where the surface of Bes drops into Imhotep, forming a layered terrain, and containing a steep cliff that defines the border between Bes and Khonsu.

3. Bastet (consolidated) – Hathor (consolidated) –Ma’at (dust): This is a region that was hard to map all at once due to the change in orientation of the head, where it transitions from the northern to the southern hemisphere. It comprises sub-region c of Bastet, a transition zone into the cliff-dominated Hathor, and Ma’at-c, containing many circular depressions that show enhanced activity on a terrain with large-scale roughness.

4. Aker (consolidated) – Bastet (consolidated) – Hapi (smooth) – Sobek (consolidated): This region sits mostly on the neck of the comet, where both head and body adjoin the smooth Hapi terrain. Both Aker-d Bastet-c are fractured cliffs, while Sobek-b shows many boulders.

5. Ash (brittle) – Aten (depression) – Babi (dust) – Khepry (consolidated): This cluster comprises the boulder-covered depression Aten and all the surrounding subregions (dust-covered Ash-a,c,i,j; dusty cliffs Babi-a; rocky, boulder-covered Khepry-a,c). Although formed by diverse terrains, the main characteristic of this cluster is given by the Aten depression, thought to have formed by one or several outburst events (El-Maarry et al. 2015).

6. Aker (consolidated) – Anhur (consolidated) – Bes (consolidated) – Khepry (consolidated): This cluster contains subregions a and b from Aker, similar to each other and separated by a ridge. Subregion a contains tectonic fractures (Thomas et al. 2015b,a), while b is similar to Bes and shows a change in slope and roughness when transitioning into Anhur-a. This is a plateau with extreme roughness, containing isolated ridges and pits. To the east it borders Bes-c, a smooth terrain that shows some boulders from collapsing cliffs in subregion d. Lastly, the northeast area is contained in Khepry-a, a flat, rock-like terrain with ponded deposits.

Fig. 9 Ejection efficiency for all the sets combined. Blue indicates patches with low efficiency, and yellow regions with enhanced efficiency. The white and gray areas represent zones that did not provide any candidate in top and bottom panels, respectively. Red and blue ellipses show the identified high and low efficiency clusters, respectively.

Low ejection efficiency

1. Hapi (smooth) – Seth (brittle): Both Hapi and Seth are fairly homogeneous regions, without any clear subdivision. The only prominent features are the exposure of very limited areas of consolidated material in Hapi, as well as active pits and semicircular depressions in Seth.

2. Anuket (consolidated): As in the last case, this region does not have any subdivision. It is smooth at large and intermediate scales, with a rocky appearance on small scales.

3. Anhur (consolidated) – Bes (consolidated) – Geb (consolidated): Although rough at an intermediate scale, the Anhur region seems quite homogeneous, with some cliffs and ridges, principally where it transitions to the neck. The border between Geb and Bes, formed by subregions Geb-c, a part of Geb-a and Bes-b, complete this cluster. At this point, Geb mainly shows a cliff, with a smooth transition to Bes. On the western end, Geb-a shows fractured terrain.

4. Hatmehit (depression) – Ma’at (dust) – Wosret (consolidated): This is a complex cluster, since it shows a very localized high efficiency spot, surrounded by a low efficiency area. The high efficiency spot is located in Hatmehit-b, at the border with Wosret and containing quasi-circular depressions. Around it there are Hatmehit-a, the flat, dust-covered bottom of the depression, and Hatmehit-c, the transition zone to Ma’at, showing a steep cliff with several fractures and layering. In the north, Ma’at-e is covered by dust, but with some of the consolidated material below exposed, while Wosret-a is located to the south, a flat, smooth surface without any major feature.

5. Ash (brittle) – Imhotep (smooth) – Khepry (consolidated): This cluster is located mainly in Imhotep-a, with some parts in subregions c and d. These are mainly smooth and covered by dust (only subregion c shows a rougher terrain). Although the adjacent areas Ash-i and Khepry-c show more terrain heterogeneity, this cluster only covers them slightly, so the mentioned areas of Imhotep are the most representatives.

Initial examination of the morphological properties within the clusters for both high and low ejection efficiencies, reveals no discernible correlation between specific features and the likelihood of aggregate ejections. The sole fundamental terrain type consistently observed in both cluster classes is consolidated, but this observation can be attributed to its prevalence as the most widespread terrain type, encompassing 25.05 km² of the total 51.74 km² surface area (48.5%).

While high ejection efficiency clusters typically align with regions distant from the axis of rotation, correlating with high centrifugal acceleration, and vice versa, exceptions exist. Notably, low ejection efficiency cluster #4, situated at the head – where centrifugal acceleration peaks – contradicts this trend. This indicates that while centrifugal acceleration may aid in the ejection of aggregates, it is not the sole determinant of this phenomenon.

Nonetheless, zones of high ejection efficiency tend to be concentrated along the boundaries between regions. These boundaries are marked by the diversity of terrains and the presence of cliffs. Conversely, regions with lower ejection efficiency tend to reside in more uniform areas, typically devoid of prominent features. The source regions of the particles analyzed in this study exhibit a distribution akin to that of dust and gas jets (Vincent et al. 2016; Fornasier et al. 2019; Lai et al. 2019). These authors attribute strong, sporadic events as the drivers behind jets and outbursts, which can arise through various mechanisms, including cliff collapses, thermal stress fractures, or the pressurization of volatile substances in deeper subsurface layers.

4.5 Afρ and the mass loss rate

In cometary studies, the Afρ parameter, introduced by A’Hearn et al. (1984), is commonly employed to characterize dust activity. Originally designed for ground-based observations, this parameter is calculated as the product of the albedo of cometary dust A, the filling factor ƒ defined as the ratio of projected area of dust particles to total sampled area, and the projected aperture radius, ρ. Afρ is a valuable parameter for quantifying comet activity because when the ejection and expansion of dust are isotropic and homogeneous, and dust properties remain constant, it is independent of the aperture radius, ρ.

The derivation of the expression for the Afρ parameter outlined in this work follows the approach of Fink & Rubin (2012) and Fink & Rinaldi (2015). While they express Afρ as a function of the total intensity detected by the sensor, we define it here in relation to the count of particles identified in OSIRIS images.

In its most general form, assuming isotropic ejection, and particles having a constant velocity (stationary coma model), the Afρ parameter for a single particle size can be expressed as (Fink & Rubin 2012) $$Af\rho =\frac{A{Q}_{\text{d}}\sigma }{2{\upsilon }_{\text{d}}}\text{,}$$(14)

where A is the dust albedo, Q_d is the production rate of particles of size r_d, σ is the particle cross section and υ_d is the particle velocity. The pivotal component of Eq. (14) is the parameter Q_d, so we formulated an expression for it using insights gained from our observations and simulations.

The first step is to obtain an expression for the aggregate column density as seen by OSIRIS. In the stationary coma model, the volumetric density of particles can be expressed as $$n_{\text{vol}}(r)=\frac{Q_{\text{d}}}{4\pi {\upsilon }_{\text{d}}}\frac{1}{r^2}.$$(15)

From this equation we derived an expression for the column density as seen by OSIRIS as follows. We defined the position of the spacecraft with respect to the nucleus as r_SC, while the camera LOS defines a new direction, ζ. The minimum distance between the LOS and the nucleus is defined as the impact parameter, ρ, and the point of minimum distance marks the origin on the ζ direction, such that the position of the spacecraft is −ζ_SC. Figure 10 shows a sketch of the geometry presented here.

The column density observed by OSIRIS can be expressed as the integral of the volumetric density over all the observed ζ. The expression for n_col reads $$\begin{array}{ll}{n}_{\text{col}}^{sc}\hfill & =\frac{Q_{\text{d}}}{4\pi {\upsilon }_{\text{d}}}\underset{-{\zeta }_{\text{sc}}}{\overset{{\zeta }_{\max }}{{\displaystyle \int }}}\frac{1}{{\zeta }^2+{\rho }^2}\text{d}\zeta \hfill \\ \hfill & =\frac{Q_{\text{d}}}{4\pi {\upsilon }_{\text{d}}}\frac{\arctan \left({\zeta }_{\max }/\rho \right)-\arctan \left(-{\zeta }_{\text{sc}}/\rho \right)}{\rho }.\hfill \end{array}$$(16)

As explained in Sect. 2, all the analyzed images were acquired in such a way that the LOS is roughly perpendicular to r_SC, so the spacecraft is located at the position ζ = 0 and ρ is the spacecraft altitude. With this, the equivalent ejection rate (the extrapolation of the production rate to the whole nucleus surface based on the particles seen by OSIRIS) can be obtained from Eq. (16): $$Q_{\text{d}}=\frac{4\pi {\upsilon }_{\text{d}}{n}_{\text{col}}\rho }{\arctan \left({\zeta }_{\max }/\rho \right)}.$$(17)

By combining the expressions in Eqs. (14) and (17), an expression for Afρ can be obtained. This expression is $$Af\rho =\frac{2\pi n_{\text{col}}\rho A\sigma }{\arctan \left({\zeta }_{\max }/\rho \right)}$$(18)

The primary factor required is the maximum distance, ζ_max, within which OSIRIS can sample particles. However, the sole information from OSIRIS images refers to the minimum detectable brightness. To define the minimum brightness for which the sample of tracks is complete, we first corrected the measured track brightness to zero phase angle using the phase function mentioned in Sect. 3.2, and then plotted the cumulative distribution of brightness. Figure 11 shows an example of this cumulative distribution for all the tracks detected in OSIRIS images acquired using the Red filter in the set STP092. Using a logarithmic scale for the corrected brightness, these plots exhibit a linear segment within the range of intermediate brightness, followed by a plateau for faint tracks. This plateau highlights the incompleteness of the sample. To address this, we fitted a straight line (in logarithmic space) to the linear segment and determined the limit brightness В_lim as the point where the fit separates from the distribution. Figure 11 shows the fit and limit brightness as dashed red and dotted black lines, respectively.

Assuming that the maximum sizes of the aggregates detected in OSIRIS images correspond to the maximum size used for the simulations (r_d = 50 cm), we find the maximum distance, ζ_msx at which the largest aggregates can be located in order to have a brightness higher than the limit, B_lim, from Eq. (9). Nonetheless, it is possible that particles situated closer than ζ_max remain undetected due to their small size. Consequently, an inclusion of a de-biasing factor is required to account for particles present within ζ_max but fainter than the limit brightness.

To determine this de-biasing factor, we simulated a set of aggregates under the conditions of the stationary coma model. We calculated the fraction of these aggregates visible to the camera with a brightness larger than B_lim. This computation yields the aforementioned de-biasing factor, which allowed us to extrapolate the complete size distribution of aggregates ranging from 5 to 50 cm in the vicinity of the spacecraft, which is depicted in Fig. 12, along with the best-fit power law for the differential distribution, characterized by an index of -3.22. This value is in agreement with other measurements in the same size range, obtained from different types of measurements, as particles observed by OSIRIS (Agarwal et al. 2016; Fulle et al. 2016a), boulders on the nucleus surface observed by the lander Philae (Mottola et al. 2015) and trail simulations based on ground-based observations (Agarwal et al. 2010; Fulle et al. 2010).

This de-biasing factor, together with the phase correction factor, can be applied to the size distribution extracted from the combination of OSIRIS images and simulations, in order to provide the corrected number of tracks observed in each image (Fig. 13 and Appendix A). With this procedure, we found the mean number of tracks per image for all the images. We found that images obtained using the Orange F22 filter showed fewer tracks per image, caused by the shorter exposure times used for this type of images compared to the ones acquired with the remaining filters. On the other hand, Blue F24 and Red F28 showed similar numbers of tracks per image, despite the different exposure times used for acquiring images using these filters. This seems to indicate a sort of saturation effect, where the increase in exposure time would not provide a larger number of tracks due to the contribution of light scattered by the diffuse coma. For this reason, the mean track number was calculated using both these filters.

Assuming all the particles are located within the distance ζ_max from the spacecraft, the column density can be found as the ratio of the total number of particles and the area base of the pyramid defining the FOV. With all these elements, we calculated the values of Afρ for each size bin using Eq. (18). These results are shown in Fig. 14, while the integrated values of Afρ for all sizes are presented in Table 4.

Our Afρ-values lie in the millimeter to centimeter range and are considerably smaller than the typical values encountered for 67P derived from ground-based telescope images (Boehnhardt et al. 2016; Snodgrass et al. 2017) or the diffuse coma as seen by Rosetta (Rinaldi et al. 2016; Gerig et al. 2018), which primarily fall in the meter scale. This outcome is expected, given that these studies predominantly focus on smaller dust particles, which are anticipated to be more abundant. However, the Afρ values presented by Fulle et al. (2016b) and Ott et al. (2017) for a similar size range as investigated here are also several orders of magnitude larger than ours. A possible cause of the discrepancy may lie their employing much shorter distances between particles and camera (below 10 km, compared to ζ_max ≃ 60 km here), which significantly augments the particle density.

Fulle et al. (2016b) finds that the aggregate loss rate can be found through the expression $$Q_n=\frac{{(Af\rho )}_b{\upsilon }_{\text{d}}}{2p\sigma },$$(19)

where (Afρ)_b represents the contribution to the cumulative Afρ originating from a particular bin, υ_d denotes the average particle velocity, which can be determined from our dynamical simulations, and p = 0.065 is the geometric albedo at 649 nm (Fornasier et al. 2015). The mass loss rate is straightforwardly defined as Q_m = Q_n× m_b, where m_b represents the characteristic mass corresponding to each bin. Table 4 shows the total particle and mass fluxes (integrated over all the bins) for the five image sets analyzed in this work.

Since the mass loss rate found here accounts for large aggregates, we can compare these values to the mass lost in smaller dust. In their work, Marschall et al. (2020) calculated the mass loss rate for a dust size distribution with a minimum size r_min = 0.1 µm by comparing the measured diffuse coma brightness with that obtained from numerical simulations; this method is more sensitive to small dust. By comparing their values with ours, we determined the ratio of mass lost in large chunks to that lost in fine dust. These results are shown in Fig. 15.

Close to perihelion, the mass lost in large aggregates represents only for around 20% of the total refractory mass lost. However, as the comet recedes in its orbit, this ratio increases up to values ~100. A comparable trend in the ratio of water and carbon dioxide production rates was observed by Läuter et al. (2020). This finding supports the hypothesis of activity in these size ranges being driven by distinct species: while small grains would be ejected by water sublimation, ejection of large chunks would be driven by CO₂ sublimation. In proximity to perihelion, elevated temperatures lead to a substantially greater water production rate. As the comet moves away from the Sun, the decrease in CO₂ production rate is less steep compared to that of water, explaining the observed change in mass loss rates.

Fig. 10 Sketch of the observation geometry for column density determination in the OSIRIS image used in this work.

Fig. 11 Estimation of limit brightness and visible aggregates. Top: cumulative distribution of phase-function-corrected track brightness for images obtained for set STP096 using the Red filter. The dashed line indicates the best fit for the linear part, while the dotted black line shows the minimum brightness for which the sample is complete. Bottom: fraction of aggregates with brightness higher than Blim at a distance ≤ζmax from the spacecraft, assuming a stationary coma model.

Fig. 12 Normalized size distribution obtained from applying the de-biasing factor to the size distribution found in Sect. 4.2. The best-fit power law, with a differential index of 3.22, is indicated with a dashed line.

Fig. 13 Number of tracks per image after phase-angle correction and brightness de-biasing as a function of the phase angle. This example shows the result for STP063. The colors indicate images obtained using different filters.

Fig. 14 Afρ parameter for the range 5-50 cm with bin width equal to 3 cm. Starting from the top left, the panels show the results for image sets STP063, STP070, STP086, STP092, and STP096.

Fig. 15 Ratio of mass lost in large aggregates to that lost in small dust. The mass loss rate for dust dominating the diffuse coma brightness (predominantly small dust) is taken from Marschall et al. (2020).

5 Conclusions

We have examined tracks from 189 images acquired by OSIRIS and organized into five sets, each obtained at a different epoch. Through dynamical simulations, we generated synthetic images for a wide range of dust parameters. By identifying the closest matches between simulations and observations, we deduced properties of the aggregates that generated the tracks and gained insight into the ejection mechanism.

Our main findings are as follows:

• The aggregates have typical mean sizes of 5–10 cm, bulk densities roughly matching that of the nucleus, and initial velocities on the order of 1 m s⁻¹ ;

• Although showing some variation, these values are consistent throughout all the analyzed time periods, ranging from 37 days before perihelion to 189 after it;

• The proposed ejection mechanism, namely CO₂ sublimating from inner layers and overcoming tensile strength, would be able to provide enough energy to the aggregate to reach the required initial velocity;

• Likely source regions on the surface correspond to boundaries between morphological regions, marked by a heterogeneity of terrain. Conversely, zones with low aggregate ejection efficiency are concentrated in smooth, heterogeneous regions. The source regions found here have similar characteristics to those found for jets by Vincent et al. (2016) and Fornasier et al. (2019);

• The Afρ parameter contributed by aggregates between 5 and 50 centimeters falls approximately in the 0.1–1 cm range for the entire studied period, much smaller than the contribution from fine dust;

• The ratio of the mass loss rate in large aggregates versus in fine dust increases with heliocentric distance, supporting the hypothesis that while small grains are ejected by water activity, CO₂ is responsible for the activity of large chunks.

A database that contains information for a grand total of 34616 tracks is available online for public consultation² (Lemos 2023).

Acknowledgements

We thank the anonymous referee for their valuable feedback. We also thank Nick Attree, Manuela Lippi and Johannes Markkanen for the helpful discussions. OSIRIS was built by a consortium of the Max-Planck-Institut für Sonnensystemforschung, Göttingen, Germany; the CISAS University of Padova, Italy; the Laboratoire d’Astrophysique de Marseille, France; the Instituto de Astrofísica de Andalucia, CSIC, Granada, Spain; the Research and Scientific Support Department of the European Space Agency Noordwijk, The Netherlands; the Instituto Nacional de Técnica Aeroespacial, Madrid, Spain; the Univer-sidad Politécnica 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. We thank the Rosetta Science Ground Segment at ESAC, the Rosetta Missions Operations Centre at ESOC and the Rosetta Project at ESTEC for their outstanding work enabling the science return of the Rosetta Mission. PL, JA and MP acknowledge funding by the ERC Starting Grant No. 757390 Comet and Asteroid Re-Shaping through Activity (CAstRA). PL and MP conducted the work in this paper in the framework of the International Max-Planck Research School (IMPRS) for Solar System Science at the University of Göttingen. JA acknowledges funding by the Volkswagen Foundation.

Appendix A Distribution of physical properties

Fig. A.1 Distribution of physical properties of the aggregates for STP063.

Fig. A.2 Distribution of physical properties of the aggregates for STP070.

Fig. A.3 Distribution of physical properties of the aggregates for STP086.

Fig. A.4 Distribution of physical properties of the aggregates for STP096.

Appendix B Source regions

Fig. B.1 Ejection efficiency for image set STP063. Crosses indicate the projected spacecraft position from 12 hours before the start of the acquisition (green) to 12 hours after the end of it (magenta). The solid cyan line indicates the subsolar latitude.

Fig. B.2 Ejection efficiency for image set STP070.

Fig. B.3 Ejection efficiency for image set STP086.

Fig. B.4 Ejection efficiency for image set STP092.

Fig. B.5 Ejection efficiency for image set STP096.

Appendix C Corrected number of tracks per image

Fig. C.1 Corrected number of tracks per image. From the top left we show STP070, STP086, STP092, and STP096.

References

All Tables

All Figures

Fig. 1 Sketch of the observation geometry for all the image sets used in this work. The dashed blue line indicates the solar direction. The dotted red line indicates the position of the spacecraft. The solid violet lines indicate the different phase angles sampled by the camera when images with the three different filters were obtained. All these sampled directions are located in the green plane, perpendicular to the nucleus-spacecraft direction. Note that the sizes are not to scale.

In the text

	Fig. 2 Diagram of the main steps involved in the detection algorithm.
In the text

	Fig. 3 Fraction of aggregates escaping the integration domain for image sets STP070 (top) and STP096 (bottom). Each black dot represents the mean of all particles in a single simulation, whereas the red symbols denote the mean across all aggregates of that particular kind, considering the cumulative data from all simulations.
In the text

	Fig. 4 Example of the distribution of physical properties for aggregates in STP092 obtained using our method. From top to bottom, rows show the distribution of radius, density, and most probable velocity. The first and second columns show candidates and fitted particles, respectively. The last column shows the ratio between heights of columns for fitted aggregates and candidates.
In the text

	Fig. 5 Mean aggregate parameters for sets STP063 (top), STP070 (middle), and STP086 (bottom) as a function of phase angle. Each row represents a different image set, while the columns from left to right represent the aggregate radius (rd), density (ρd), and initial velocity (υin). The colors show the results for images taken with filters Blue F24, Orange F22, and Red F28.
In the text

	Fig. 6 Same as Fig. 5 but for image sets STPO92 (top) and STPO96 (bottom).
In the text

	Fig. 7 Sketch of the breaking process. Sublimated gas occupies the void space below the aggregate. When pressure equals tensile strength, the aggregate is ejected and a net force coming from the gas pressure acts on it, producing the observed initial velocity.
In the text

	Fig. 8 Initial velocity of aggregates as a function of the volume of the pore containing the gas, as determined by Eq. (12). The blue line corresponds to 5 cm radius aggregates, and the red line 10 cm radius aggregates. The horizontal dashed line indicates an initial velocity of 1 ms−1.
In the text

	Fig. 9 Ejection efficiency for all the sets combined. Blue indicates patches with low efficiency, and yellow regions with enhanced efficiency. The white and gray areas represent zones that did not provide any candidate in top and bottom panels, respectively. Red and blue ellipses show the identified high and low efficiency clusters, respectively.
In the text

	Fig. 10 Sketch of the observation geometry for column density determination in the OSIRIS image used in this work.
In the text

Fig. 11 Estimation of limit brightness and visible aggregates. Top: cumulative distribution of phase-function-corrected track brightness for images obtained for set STP096 using the Red filter. The dashed line indicates the best fit for the linear part, while the dotted black line shows the minimum brightness for which the sample is complete. Bottom: fraction of aggregates with brightness higher than Blim at a distance ≤ζmax from the spacecraft, assuming a stationary coma model.

In the text

	Fig. 12 Normalized size distribution obtained from applying the de-biasing factor to the size distribution found in Sect. 4.2. The best-fit power law, with a differential index of 3.22, is indicated with a dashed line.
In the text

	Fig. 13 Number of tracks per image after phase-angle correction and brightness de-biasing as a function of the phase angle. This example shows the result for STP063. The colors indicate images obtained using different filters.
In the text

	Fig. 14 Afρ parameter for the range 5-50 cm with bin width equal to 3 cm. Starting from the top left, the panels show the results for image sets STP063, STP070, STP086, STP092, and STP096.
In the text

	Fig. 15 Ratio of mass lost in large aggregates to that lost in small dust. The mass loss rate for dust dominating the diffuse coma brightness (predominantly small dust) is taken from Marschall et al. (2020).
In the text

	Fig. A.1 Distribution of physical properties of the aggregates for STP063.
In the text

	Fig. A.2 Distribution of physical properties of the aggregates for STP070.
In the text

	Fig. A.3 Distribution of physical properties of the aggregates for STP086.
In the text

	Fig. A.4 Distribution of physical properties of the aggregates for STP096.
In the text

	Fig. B.1 Ejection efficiency for image set STP063. Crosses indicate the projected spacecraft position from 12 hours before the start of the acquisition (green) to 12 hours after the end of it (magenta). The solid cyan line indicates the subsolar latitude.
In the text

	Fig. B.2 Ejection efficiency for image set STP070.
In the text

	Fig. B.3 Ejection efficiency for image set STP086.
In the text

	Fig. B.4 Ejection efficiency for image set STP092.
In the text

	Fig. B.5 Ejection efficiency for image set STP096.
In the text

	Fig. C.1 Corrected number of tracks per image. From the top left we show STP070, STP086, STP092, and STP096.
In the text