11institutetext: Institut für Geophysik und Extraterrestrische Physik,
Technische Universität Braunschweig, Mendelssohnstraße 3, D-38106
Braunschweig, Germany.
11email<EMAIL_ADDRESS>22institutetext: Max Planck
Institute for Solar System Research, Justus-von-Liebig-Weg 3, D-37077
Göttingen, Germany. 33institutetext: CNRS, Laboratoire J.-L. Lagrange,
Observatoire de la Côte d’Azur, Nice, France.

# Ejection and Dynamics of Aggregates in the Coma of Comet 67P/Churyumov-
Gerasimenko

P. Lemos 1122 J. Agarwal 1122 R. Marschall 33 M. Pfeifer 22

(Received ¡date¿; accepted ¡date¿)

###### Abstract

Context. The process of gas-driven ejection of refractory materials from
cometary surfaces continues to pose a challenging question in cometary
science. The activity modeling of comet 67P/Churyumov-Gerasimenko, based on
data from the Rosetta mission, has significantly enhanced our comprehension of
cometary activity. But thermophysical models have difficulties in
simultaneously explaining the production rates of various gas species and
dust. It has been suggested that different gas species might be responsible
for the ejection of refractory material in distinct size ranges.

Aims. This work focuses on investigating abundance and the ejection mechanisms
of large ($\gtrsim$ 1 cm) aggregates from the comet nucleus. We aim to
determine their properties and map the distribution of their source regions
across the comet surface. This can place constraints on activity models for
comets.

Methods. We examined 189 images acquired at five epochs by the OSIRIS/NAC
instrument on board the Rosetta spacecraft. Our goal was to identify bright
tracks produced by individual aggregates as they traversed the camera field of
view. In parallel, we generated synthetic images based on the output of
dynamical simulations involving various types of aggregates. By comparing
these synthetic images with the observations, we determine the characteristics
of the simulated aggregates that most closely resembled the observations.

Results. We identified over 30000 tracks present in the OSIRIS images, derived
constraints on the characteristics of the aggregates and mapped their origins
on the nucleus surface. The aggregates have an average radius of $\simeq 5$
cm, and a bulk density consistent with that of the comet’s nucleus. Due to
their size, gas drag exerts only a minor influence on their dynamical
behavior, so an initial velocity is needed in order to bring them into the
camera field of view. The source regions of these aggregates are predominantly
located near the boundaries of distinct terrains on the surface.

###### Key Words.:

methods: data analysis – methods: numerical – comets: individual:
67P/Churyumov–Gerasimenko

##  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 regarding all size ranges can be found in
Güttler et al. (2019).

While most of dust analysis instruments onboard Rosetta were designed for
particles with typical sizes $\lesssim 1$ mm, the mission offered a unique
opportunity to observe and study larger particles, exceeding this size range,
using the OSIRIS cameras (Keller et al. 2007). The main challenge about using
this data comes from the lack of information about their distances and
velocities along the line of sight 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 employ a method introduced by 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, the method consists on the
analysis of OSIRIS images, where 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, are
extracted from the images.

In parallel, dynamical simulations are conducted in order to trace the
trajectories of particles. Based on the position of the spacecraft and the
camera viewing direction, we generate synthetic images for the times at which
the OSIRIS data were obtained. Then, the properties of the tracks within the
synthetic images are inferred, and their distributions compared to their real
counterparts. The simulated particles generating 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 of
various OSIRIS images employing the previously outlined method. This approach
enables us to extract 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 Section 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, Section 5
summarizes the key findings of this work.

##  2 Observations

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 $\alpha$ 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 broad-band 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.

Table 1: Image sets used for this work. Columns represent the mid and short-term planning cycles, date of acquisition, filters used for the acquisition, heliocentric distances, nucleocentric distances and number of images in the set. Planning cycle (MTP/STP) | Date |  $r_{h}$ (au) |  $r_{S/C}$ (km) | $\\#I$
---|---|---|---|---
018/063 | 2015-07-07 | 1.32 | 153.4 | 45
019/070 | 2015-08-20 | 1.25 | 325.2 | 39
023/086 | 2015-12-14 | 1.89 | 102.6 | 21
025/092 | 2016-01-21 | 2.18 | 79.2 | 39
026/096 | 2016-02-18 | 2.39 | 36.9 | 45
Figure 1: Sketch of the observation geometry for all the image sets used in
this work. The blue dashed line indicates the solar direction. The red dotted
line indicates the position of the spacecraft. The violet solid lines indicate
the different phase angles sampled by the camera, where images with the three
different filters were obtained. All these sampled directions belong to the
green plane, perpendicular to the nucleus-spacecraft direction. Note that the
sizes are not to scale.

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, e.g., Bertini et al.
(2017) and Keiser et al. (forthcoming) to constrain the spatial distribution
of scattering cross-section and the dust scattering properties, we here
concentrate on the individual tracks in the foreground, in a data set
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 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).

###  2.1 Detection method

We are interested in analyzing the properties of the aggregates with typical
sizes $\gtrsim 1$ cm. As mentioned in Section 1, these aggregates are seen as
bright tracks in the studied images. We use a semi-automatic method that
exploits this property in order 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. 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. 2.

Create binary images from the similarity maps.

  3. 3.

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

  4. 4.

Refine the results in order to correct imperfections of the algorithm

  5. 5.

Manually inspect the results.

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

Figure 2: Diagram of the main steps involved in the detection algorithm.

By applying this procedure to the groups of images described before, we
detected 34616 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 use 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 $\sim 440000$ facets with typical
length-scales of $\simeq$ 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) in
order to solve the dust aggregates equations of motion. This code uses a 4th
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 $\simeq 30$ m
and summing the contributions of each element, assuming an homogeneous density
and a total mass $M_{N}=9.9\times 10^{12}$ kg (Sierks et al. 2015).

The gas drag force $\mathbf{F}_{D}$ is calculated using the equation

$\mathbf{F}_{D}=\frac{1}{2}C_{D}\,m_{g}\,n_{g}\,\sigma\,|\mathbf{v}_{g}-\mathbf{v}_{d}|\,(\mathbf{v}_{g}-\mathbf{v}_{d}),$
(1)

where $m_{g}$, $n_{g}$ and $\mathbf{v}_{g}$ are the molecular mass, number
density and velocity of the gas respectively, $\mathbf{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_{D}=\frac{2s^{2}+1}{\sqrt{\pi}s^{3}}\mathrm{e}^{-s^{2}}+\frac{4s^{4}+4s^{2}-1}{2s^{4}}\mathrm{erf}(s)+\frac{2\sqrt{\pi}}{3s}\sqrt{\frac{T_{d}}{T_{g}}},$
(2)

where $T_{g}$ and $T_{d}$ are the gas and dust temperatures respectively
(assumed to be equal) and

$s=\frac{|\mathbf{v}_{g}-\mathbf{v}_{d}|}{\sqrt{\frac{2k_{B}T_{g}}{m_{g}}}}.$
(3)

Solar radiation pressure $\mathbf{F}_{R}$ and solar tidal forces
$\mathbf{F}_{T}$ are decribed as follows

$\displaystyle\mathbf{F}_{R}$ $\displaystyle=-\frac{C_{\odot}Q_{RP}\pi
r_{d}^{2}}{r_{h}^{3}c}\mathbf{r}_{h}$ (4) $\displaystyle\mathbf{F}_{T}$
$\displaystyle=\mathbf{F}_{S}(\mathbf{r}_{C})-\mathbf{F}_{S}(\mathbf{r}_{C}+\mathbf{r}_{D}),$
(5)

where $C_{\odot}=1361$ W m-2 is the solar constant, $Q_{RP}$ is the radiation
pressure efficiency (equal to 1), $c$ is the speed of light,
$\mathbf{F}_{S}(\mathbf{r})$ is the solar gravity force at $\mathbf{r}$,
$\mathbf{r}_{C}$ is the Sun-comet vector and $\mathbf{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 of the order or much longer than the
rotation period of the nucleus. For this work we will assume that the dust
aggregates do not sublimate volatiles from their surface. For each simulation,
a minimum of one and a maximum of ten aggregates are 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 are 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 $v_{P}$. With this, the aggregates are
characterized by three parameters: particle radius $r_{d}$, particle density
$\rho_{d}$ and most probable initial velocity $v_{P}$. The values used for
these parameters are shown in Table 2. The simulations using different
combinations of these parameters are carried out independently, and the
solutions are later combined in a procedure explained in Section 3.2.

Table 2: Values employed for the aggregate parameter in the dynamic simulations. Parameter | Values
---|---
$\rho_{d}$ |  $[100;500;800]$ kg m-3
$r_{d}$ |  $[1;5;10;50]$ cm
$v_{P}$ |  $[0;0.1;0.5;1.0;2.0]$ m s-1

Since the dust flight time needed to reach the camera FOV is 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 $\times$ 60 aggregate parameters combinations $=$ 720 cases.
The total number of aggregates simulated in each case is variable, depending
on the distribution of water vapour production on the surface, but is always
in the order of $5\times 10^{5}$ 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 field of view (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, in 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 in order 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 non-rotating 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 is deemed valid within the range of $[-15,+15]\degree$ from its
nominal rotation angle. The trajectories are parametrized with respect to an
angle $\theta_{I}$ that belongs to that range. Using this parametrization, the
minimum distance between the particle trajectory and the line of sight was
determined as a function of $\theta_{I}$. Should the angle between the vector
connecting the spacecraft to the aforementioned point of minimal distance and
the line of sight direction be smaller than the angular dimension of the FOV
($2.20\degree$), the aggregate in question is 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 $\theta_{i}\in\theta_{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 are associated with
the exposure time of the real image. In order to define those endpoints, we
selected a random position between entry and exit points. Then, the exposure
time of the image was divided in two intervals (not necessarily equal), and
the position of this random point is propagated both forward and backward
along the trajectory, defining both track endpoints. As we mentioned before,
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 are 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 sub-solar longitude at the
moment of observation that is recorded on the image header $\lambda_{obs}$
must be compatible with that of the simulations, obtained by adding the sub-
solar longitude at the moment of ejection $\lambda_{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_{obs}=\lambda_{ej}+\omega_{R}\times t_{F},$ (6)

where $\omega_{R}$ is the nucleus angular velocity of rotation.

As mentioned before, $\lambda{ej}$ does not necessarily match with the sub-
solar longitude of the gas solutions $\lambda_{gas}$, since the angle
$\theta_{i}$ was added in order to account for discontinuities. Plugging this
into Eq. 6 we obtain

$\displaystyle\lambda_{obs}=\lambda_{gas}+\theta_{i}+\omega_{R}\times t_{F}$
(7) $\displaystyle\theta_{i}=\lambda_{obs}-\lambda_{gas}-\omega_{R}\times
t_{F}.$ (8)

For a trajectory to be accepted, the angle $\theta_{i}$ found from Equation 8
has to belong to the range defined by the angles $\theta_{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\Phi(\alpha)I_{\odot}r_{d}^{2}}{r_{h}^{2}\Delta^{2}},$ (9)

where $p$ is the aggregate geometric albedo, $\Phi(\alpha)$ its phase function
normalized to 1 at zero phase, $I_{\odot}$ the solar flux in the corresponding
filter, $r_{d}$ is the dust aggregate radius, $r_{h}$ the heliocentric
distance in au and $\Delta$ is the distance between the aggregate and the
camera. In case 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
$\Phi(\alpha)=\exp{(-\beta\times\alpha)}$, with a mean $\beta=8.2\times
10^{-3}$. We used $p=0.065$ (Agarwal et al. 2016) and defined $\Delta$ as the
mean of the distances between the camera and entry and exit points
respectively.

###  3.3 Parameter optimization by comparison between observed and synthetic
tracks

This comparison was based on the three track properties mentioned in Section
2.1: orientation angle, length and integrated brightness. However, since the
number of tracks detected in different OSIRIS images is highly variable, we
generated an estimation of the probability density function (PDF) for the
distribution of properties from observed tracks using a Gaussian kernel
density estimation. Based on this PDF estimation, we created a new sample of
track property triplets which represents 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 moment, 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 neighbours 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 best fits to the tracks observed in
OSIRIS images. In this Section, we will 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 will 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).

Figure 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 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 Section 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 Section
3.3 are called fitted aggregates. In order 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.

Figure 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, while first and
second columns show candidates and fitted particles respectively. Last column
shows the ratio between heights of columns for fitted aggregates and
candidates.

The data in Figure 4 provides 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.

Table 3: Mean of each dust particle parameter considering all the images in the sets. Image set |  $\rho_{d}$ (kg m-3) |  $r_{d}$ (cm) |  $v_{P}$ (m s-1)
---|---|---|---
STP063 | 458 | 5.2 | 0.65
STP070 | 434 | 7.0 | 0.80
STP086 | 549 | 6.4 | 1.30
STP092 | 468 | 4.5 | 0.68
STP096 | 471 | 8.1 | 1.08

Figures 5-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 angles111Note that while the
velocity shown in Table 3 represents the most probable initial velocity
$v_{P}$, used to define the MB distribution, Figures 5-6 represent the initial
velocity $v_{in}$ obtained from that distribution.. 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.

Figure 5: Mean aggregate parameters for sets STP063 (top), STP070 (middle) and
STP086 (bottom) as a function of phase angle. Each row represent a different
image set, while the columns from left to right represent aggregate radius
$r_{d}$, density $\rho_{d}$ and initial velocity $v_{in}$ respectively. The
colors show the results for images taken with filters Blue F24, Orange F22 and
Red F28.

Although being smaller than the ones found in Lemos et al. (2023), the
aggregate sizes found in this work match well with the ones 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 vapour sublimation, so more volatile species as CO2 seem to be
responsible for their ejection (Gundlach et al. 2020).

Our results also indicate that the aggregates have a density smaller than the
one measured for dust particles using GIADA $\rho_{D}\simeq 800$ kg m-3 (Fulle
et al. 2016a, 2017), but comparable with the bulk nucleus density
$\rho_{N}=537.8$ kg m-3 (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.

Figure 6: Same as Figure 5 for image sets STP092 (top) and STP096 (bottom).

###  4.3 Initial velocities

Our method shows that in order to reach the FOV, the aggregates require an
initial velocity $v_{in}$ on the order of 1 m s-1. 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 in
case that dust is ejected from a porous, sub-surface 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 non-zero initial velocity for the
aggregates.

In this work, we will investigate 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.

Figure 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 a rough approximation, we will assume 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

$\displaystyle K_{agg}$ $\displaystyle=U_{gas}$
$\displaystyle\frac{1}{2}\,\frac{4}{3}\pi\,r_{d}^{3}\,\rho_{d}\,v_{in}^{2}$
$\displaystyle=C_{V}\,n\,T,$ (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, Equation 4.3 becomes

$v_{in}=\sqrt{\frac{3\,C_{V}\,P\,V}{2\,\pi\,r_{d}^{3}\,R\,\rho_{d}}}.$ (11)

At the moment of ejection gas pressure equals the material tensile strength
$\sigma_{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_{T}=\frac{G}{\sqrt{r_{d}}}$, where $G=100$ Pa m1/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 $f$ represents the number
of degrees of freedom of the molecule. Inserting this into Equation 11 yields

$v_{in}=\sqrt{\frac{3\,f\,G\,V}{4\,\pi\,\rho_{d}r_{d}^{7/2}}}.$ (12)

Since thermophysical models suggest that CO2 is the most plausible candidate
for explaining the ejection of chunks with these sizes. If the vibrational
modes are not taken into account, $f$ 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
$\rho_{d}$ =500 kg m-3 and radii $r_{d}=$5 and 10 cm. For a 5 cm aggregate,
the volume needed to provide a 1 m s-1 initial velocity is $1.2\times 10^{-4}$
m-3, while for a 10 cm aggregate this volume increases to $1.3\times 10^{-3}$
m-3. The minimal porosity that provides the necessary energy is expressed by
the formula $\phi_{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.

Figure 8: Initial velocity of aggregates as a function of the volume of the
pore containing the gas, as determined by Equation 12. The blue line
corresponds to 5 cm radius aggregates, while the red line indicates 10 cm
radius aggregates. The horizontal dashed line indicates an initial velocity of
1 m s-1.

###  4.4 Source regions

One advantage of our model is its utilization of a complete shape model,
allowing 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
non-spherical 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 Figure 9. Below, we enumerate these clusters,
utilizing the regions defined by El-Maarry et al. (2015, 2016) and the sub-
regions 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.

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

High ejection efficiency

  1. 1.

Ash (Dust) - Apis (Consolidated) - Imhotep (Smooth) - Khonsu (Consolidated):
this cluster represents the interface between three very different regions. It
comprises the sub-region c of Imhotep up to its boundary with sub-region a,
characterized by a intermediate scale roughness and the presence of layers,
thought to be exposed by material loss. On the north, Ash-i shows large-scale
roughness and layering in cliffs, while Apis is rough in many scales, and
similar to the 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. 2.

Bes (Consolidated) - Imhotep (Smooth) - Khonsu (Consolidated): this region
runs along the interface between Khonsu-c on the northeast, Bes-a on the south
and Imhotep-b on 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. 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 north to south 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. 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. 5.

Ash (Brittle) - Aten (Depression) - Babi (Dust) - Khepry (Consolidated): this
cluster comprises the boulder-covered depression Aten and all the surrounding
sub-regions (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. 6.

Aker (Consolidated) - Anhur (Consolidated) - Bes (Consolidated)-Khepry
(Consolidated): this cluster contains the sub-regions a and b from Aker,
similar to each other and separated by a ridge. Sub-region 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. On the
east it limits with Bes-c, a smooth terrain that shows some boulders from
collapsing of cliffs in sub-region d. Lastly, the northeast area is contained
in Khepry-a, a flat, rock-like terrain with ponded deposits.

Low ejection efficiency

  1. 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 semi-circular depressions in Seth.

  2. 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 at small scale.

  3. 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 sub-regions 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. 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. On the north, Ma’at-e is covered by
dust, but with some of the consolidated material below being exposed, while on
the south Wosret-a is located, a flat, smooth surface without any major
feature.

  5. 5.

Ash (Brittle) - Imhotep (Smooth) - Khepry (Consolidated): this cluster is
located mainly in Imhotep-a, with some parts in sub-regions c and d. These are
mainly smooth and covered by dust (only sub-region 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 km2 of the total 51.74 km2
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 sub-surface
layers.

###  4.5 $Af\rho$ and mass loss rate

In cometary studies, the $Af\rho$ 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 $f$ defined as the ratio
of projected area of dust particles to total sampled area, and the projected
aperture radius $\rho$. $Af\rho$ 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 $\rho$.

The derivation of the expression for the $Af\rho$ parameter outlined in this
work follows the approach of Fink & Rubin (2012) and Fink & Rinaldi (2015).
While they express $Af\rho$ as a function of the total intensity detected by
the sensor, we here define it 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\rho$ parameter for a single
particle size can be expressed as (Fink & Rubin 2012)

$Af\rho=\frac{A\,Q_{d}\,\sigma}{2\,v_{d}},$ (14)

where $A$ is the dust albedo, $Q_{d}$ is the production rate of particles of
size $r_{d}$, $\sigma$ is the particle cross section and $v_{d}$ is the
particle velocity. The pivotal component of Equation 14 is the parameter
$Q_{d}$, so we will now formulate an expression for it using insights gained
from our observations and simulations.

In the stationary coma model, the volumetric density of particles can be
expressed as

$n_{vol}(r)=\frac{Q_{d}}{4\pi v_{d}}\frac{1}{r^{2}}.$ (15)

From this Equation, we derive an expression for the column density as seen by
OSIRIS as follows. We define the position of the spacecraft with respect to
the nucleus as $\mathbf{r}_{SC}$, while the camera line of sight (LOS) defines
a new direction $\mathbf{\zeta}$. The minimum distance between the LOS and the
nucleus is defined as the impact parameter $\rho$, and the point of minimum
distance marks the origin on the $\mathbf{\zeta}$ direction, such that the
position of the spacecraft is $-\zeta_{SC}$. Figure 10 shows a sketch of the
geometry presented here.

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

The column density observed by OSIRIS can be expressed as the integral of the
volumetric density over all the observed $\zeta$

$\displaystyle n_{col}^{sc}$ $\displaystyle=\frac{Q_{d}}{4\pi
v_{d}}\int_{-\zeta_{sc}}^{\zeta_{max}}\frac{1}{\zeta^{2}+\rho^{2}}d\zeta$
$\displaystyle=\frac{Q_{d}}{4\pi
v_{d}}\frac{\arctan{(\zeta_{max}/\rho)}-\arctan{(-\zeta_{sc}/\rho)}}{\rho}.$
(16)

As explained in Section 2, all the analyzed images were acquired in such a way
that the LOS is roughly perpendicular to $\mathbf{r}_{SC}$, so the spacecraft
is located at the position $\zeta=0$ and $\rho$ 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 Equation 4.5

$Q_{d}=\frac{4\,\pi\,v_{d}\,n_{col}\,\rho}{\arctan(\zeta_{max}/\rho)}.$ (17)

Combining Equations 14 and 17, we find the following expression for $Af\rho$

$Af\rho=\frac{2\,\pi\,n_{col}\,\rho\,A\,\sigma}{\arctan(\zeta_{max}/\rho)}.$
(18)

The primary factor required is the maximum distance $\zeta_{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
Section 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 $B_{lim}$ as the point
where the fit separates from the distribution. Figure 11 shows the fit and
limit brightness as red dashed and black dotted lines respectively.

Figure 11: 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 black dotted line
shows the minimum brightness for which the sample is complete. Bottom:
Fraction of aggregates with brightness larger than $B_{lim}$ at a distance
$\leq\zeta_{max}$ from the spacecraft, assuming a stationary coma model.

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, $\zeta_{max}$ at which the largest aggregates can
be located in order to have a brightness higher than the limit, $B_{lim}$,
from Equation 9. Nonetheless, it is possible that particles situated closer
than $\zeta_{max}$ remain undetected due to their small size. Consequently, an
inclusion of a debiasing factor is required to account for particles present
within $\zeta_{max}$ but fainter than the limit brightness.

To determine this debiasing 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 debiasing factor, which allows to
extrapolate the complete size distribution of aggregates ranging from 5 to 50
cm in the vicinity of the spacecraft, which is depicted in Figure 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).

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

This debiasing 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 (Figure 13 and Appendix ). 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 of 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.

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

Assuming all the particles are located within the distance $\zeta_{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\rho$ for each size bin
using Equation 18. These results are shown in Figure 14, while the integrated
values of $Af\rho$ for all sizes are presented in Table 4.

Figure 14: $Af\rho$ parameter for the range 5-50 cm with bin width equal to 3 cm. Starting from top left, the panels show the results for image sets STP063, STP070, STP086, STP092 and STP096 respectively. Table 4: $Af\rho$, particle flux $Q_{n}$ and mass flux $Q_{m}$ integrated for all particle sizes in the range 5-50 cm. Image set | Days from perihelion |  $Af\rho$ (cm) |  $Q_{n}$ (s-1) |  $Q_{m}$ (kg s-1)
---|---|---|---|---
STP063 | -37 | 0.31 | 2.44 | 12.9
STP070 | 7 | 1.73 | 13.27 | 49.5
STP086 | 123 | 0.33 | 1.78 | 8.2
STP092 | 161 | 0.07 | 0.28 | 1.6
STP096 | 189 | 0.15 | 0.84 | 7.2

Our $Af\rho$-values lie in the mm-cm-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, also the $Af\rho$ values presented by Fulle et al.
(2016b) and Ott et al. (2017) for a similar size range as investigated here
are 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 $\zeta_{max}\simeq 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}\,v_{d}}{2\,A_{p}\,\sigma},$ (19)

where $(Af\rho)_{b}$ represents the contribution to the cumulative $Af\rho$
originating from a particular bin, $v_{d}$ denotes the average particle
velocity, that can be determined from our dynamical simulations and
$A_{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}\times 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) calculate the mass loss rate for a dust size
distribution with a minimum size $r_{min}=0.1$ $\mu$m by comparing the
measured diffuse coma brightness with the one obtained from numerical
simulations, method more sensitive to small dust. By comparing their values
with ours, we determine the ratio of mass lost in large chunks to that lost in
fine dust. These results are shown in Figure 15.

Figure 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).

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 $\sim 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 CO2
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 CO2 production rate is less steep compared to that of
water, explaining the observed change in mass loss rates.

##  5 Conclusions

We 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 simulation and observation, we deduced properties
of the aggregates generating the tracks and gained insight into the ejection
mechanism.

Our main findings are listed below.

  * •

The aggregates have typical mean sizes of $5-10$ cm, bulk densities roughly
matching that of the nucleus, and initial velocities of order 1 m s-1.

  * •

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

  * •

The proposed ejection mechanism, namely CO2 sublimating from inner layers and
overcoming tensile strength, would be able to provide enough energy to the
aggregate in order 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 concentrate on 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\rho$ parameter contributed by aggregates between 5 and 50 centimeter
falls approximately in the $0.1-1$ cm range for all the studied period, much
smaller than the contribution by fine dust.

  * •

The ratio of 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, CO2 is responsible for activity of large
chunks.

The database, containing information from 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 Universidad
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.

## References

  * Acton et al. (2018) Acton, C., Bachman, N., Semenov, B., & Wright, E. 2018, Planetary and Space Science, 150, 9, enabling Open and Interoperable Access to Planetary Science and Heliophysics Databases and Tools
  * Acton (1996) Acton, C. H. 1996, Planetary and Space Science, 44, 65, planetary data system
  * Agarwal et al. (2016) Agarwal, J., A’Hearn, M. F., Vincent, J. B., et al. 2016, MNRAS, 462, S78
  * Agarwal et al. (2010) Agarwal, J., Müller, M., Reach, W. T., et al. 2010, Icarus, 207, 992
  * A’Hearn et al. (2011) A’Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2011, Science, 332, 1396
  * A’Hearn et al. (1984) A’Hearn, M. F., Schleicher, D. G., Millis, R. L., Feldman, P. D., & Thompson, D. T. 1984, AJ, 89, 579
  * Belton (2010) Belton, M. J. S. 2010, Icarus, 210, 881
  * Bertini et al. (2017) Bertini, I., La Forgia, F., Tubiana, C., et al. 2017, MNRAS, 469, S404
  * Biele et al. (2022) Biele, J., Vincent, J.-B., & Knollenberg, J. 2022, Universe, 8, 487
  * Bischoff et al. (2019) Bischoff, D., Gundlach, B., Neuhaus, M., & Blum, J. 2019, MNRAS, 483, 1202
  * Boehnhardt et al. (2016) Boehnhardt, H., Riffeser, A., Kluge, M., et al. 2016, MNRAS, 462, S376
  * Drolshagen et al. (2017) Drolshagen, E., Ott, T., Koschny, D., et al. 2017, Planet. Space Sci., 143, 256
  * Duda & Hart (1972) Duda, R. O. & Hart, P. E. 1972, Commun. ACM, 15, 11–15
  * El-Maarry et al. (2017) El-Maarry, M. R., Groussin, O., Thomas, N., et al. 2017, Science, 355, 1392
  * El-Maarry et al. (2015) El-Maarry, M. R., Thomas, N., Giacomini, L., et al. 2015, A&A, 583, A26
  * El-Maarry et al. (2016) El-Maarry, M. R., Thomas, N., Gracia-Berná, A., et al. 2016, A&A, 593, A110
  * Fink & Rinaldi (2015) Fink, U. & Rinaldi, G. 2015, Icarus, 257, 9
  * Fink & Rubin (2012) Fink, U. & Rubin, M. 2012, Icarus, 221, 721
  * Fornasier et al. (2015) Fornasier, S., Hasselmann, P. H., Barucci, M. A., et al. 2015, A&A, 583, A30
  * Fornasier et al. (2019) Fornasier, S., Hoang, V. H., Hasselmann, P. H., et al. 2019, A&A, 630, A7
  * Frattin et al. (2021) Frattin, E., Bertini, I., Ivanovski, S. L., et al. 2021, MNRAS, 504, 4687
  * Frattin et al. (2017) Frattin, E., Cremonese, G., Simioni, E., et al. 2017, MNRAS, 469, S195
  * Fulle et al. (2010) Fulle, M., Colangeli, L., Agarwal, J., et al. 2010, A&A, 522, A63
  * Fulle et al. (2017) Fulle, M., Della Corte, V., Rotundi, A., et al. 2017, MNRAS, 469, S45
  * Fulle et al. (2016a) Fulle, M., Della Corte, V., Rotundi, A., et al. 2016a, MNRAS, 462, S132
  * Fulle et al. (2015) Fulle, M., Ivanovski, S. L., Bertini, I., et al. 2015, A&A, 583, A14
  * Fulle et al. (2016b) Fulle, M., Marzari, F., Della Corte, V., et al. 2016b, ApJ, 821, 19
  * Gerig et al. (2018) Gerig, S. B., Marschall, R., Thomas, N., et al. 2018, Icarus, 311, 1
  * Gundlach et al. (2020) Gundlach, B., Fulle, M., & Blum, J. 2020, MNRAS, 493, 3690
  * Güttler et al. (2017) Güttler, C., Hasselmann, P. H., Li, Y., et al. 2017, MNRAS, 469, S312
  * Güttler et al. (2019) Güttler, C., Mannel, T., Rotundi, A., et al. 2019, A&A, 630, A24
  * Hough (1962) Hough, P. V. 1962, Method and means for recognizing complex patterns, uS Patent 3,069,654
  * Huebner et al. (2006) Huebner, W. F., Benkhoff, J., Capria, M.-T., et al. 2006, Heat and gas diffusion in comet nuclei, Vol. 133 (International Space Science Institute Bern)
  * Ivanovski et al. (2017a) Ivanovski, S. L., Della Corte, V., Rotundi, A., et al. 2017a, MNRAS, 469, S774
  * Ivanovski et al. (2017b) Ivanovski, S. L., Zakharov, V. V., Della Corte, V., et al. 2017b, Icarus, 282, 333
  * Keiser et al. (forthcoming) Keiser, F., Markkanen, J., & Agarwal, J. forthcoming, A&A
  * Keller et al. (2007) Keller, H. U., Barbieri, C., Lamy, P., et al. 2007, Space Sci. Rev., 128, 433
  * Kramer & Noack (2015) Kramer, T. & Noack, M. 2015, ApJ, 813, L33
  * Kramer & Noack (2016) Kramer, T. & Noack, M. 2016, ApJ, 823, L11
  * Kwon et al. (2023) Kwon, Y. G., Opitom, C., & Lippi, M. 2023, A&A, 674, A206
  * Lai et al. (2019) Lai, I. L., Ip, W. H., Lee, J. C., et al. 2019, A&A, 630, A17
  * Läuter et al. (2020) Läuter, M., Kramer, T., Rubin, M., & Altwegg, K. 2020, MNRAS, 498, 3995
  * Lemos (2023) Lemos, P. 2023, Track properties as observed from OSIRIS-Rosetta
  * Lemos et al. (2023) Lemos, P., Agarwal, J., & Schröter, M. 2023, MNRAS, 519, 5775
  * Lin et al. (2017) Lin, Z.-Y., Knollenberg, J., Vincent, J. B., et al. 2017, MNRAS, 469, S731
  * Marschall et al. (2020) Marschall, R., Markkanen, J., Gerig, S.-B., et al. 2020, Frontiers in Physics, 8, 227
  * Marschall et al. (2016) Marschall, R., Su, C. C., Liao, Y., et al. 2016, A&A, 589, A90
  * Möller & Trumbore (1997) Möller, T. & Trumbore, B. 1997, J. Graph. Tools, 2, 21–28
  * Moreno et al. (2022) Moreno, F., Guirado, D., Muñoz, O., et al. 2022, MNRAS, 510, 5142
  * Mottola et al. (2015) Mottola, S., Arnold, G., Grothues, H. G., et al. 2015, Science, 349, 2.232
  * Ott et al. (2017) Ott, T., Drolshagen, E., Koschny, D., et al. 2017, MNRAS, 469, S276
  * Pfeifer et al. (2024) Pfeifer, M., Agarwal, J., Marschall, R., Grieger, B., & Lemos, P. 2024, arXiv e-prints, arXiv:2402.18613
  * Pfeifer et al. (2022) Pfeifer, M., Agarwal, J., & Schröter, M. 2022, A&A, 659, A171
  * Preusker et al. (2017) Preusker, F., Scholten, F., Matz, K. D., et al. 2017, A&A, 607, L1
  * Preusker et al. (2015) Preusker, F., Scholten, F., Matz, K. D., et al. 2015, A&A, 583, A33
  * Rinaldi et al. (2016) Rinaldi, G., Fink, U., Doose, L., et al. 2016, MNRAS, 462, S547
  * Shi et al. (2024) Shi, X., Hu, X., Agarwal, J., et al. 2024, ApJ, 961, L16
  * Sierks et al. (2015) Sierks, H., Barbieri, C., Lamy, P. L., et al. 2015, Science, 347, aaa1044
  * Snodgrass et al. (2017) Snodgrass, C., A’Hearn, M. F., Aceituno, F., et al. 2017, Philosophical Transactions of the Royal Society of London Series A, 375, 20160249
  * Thomas et al. (2015a) Thomas, N., Davidsson, B., El-Maarry, M. R., et al. 2015a, A&A, 583, A17
  * Thomas et al. (2018) Thomas, N., El Maarry, M. R., Theologou, P., et al. 2018, Planet. Space Sci., 164, 19
  * Thomas et al. (2015b) Thomas, N., Sierks, H., Barbieri, C., et al. 2015b, Science, 347, aaa0440
  * Tubiana et al. (2015) Tubiana, C., Güttler, C., Kovacs, G., et al. 2015, A&A, 583, A46
  * Vincent et al. (2016) Vincent, J. B., A’Hearn, M. F., Lin, Z. Y., et al. 2016, MNRAS, 462, S184
  * Wesołowski et al. (2020) Wesołowski, M., Gronkowski, P., & Tralle, I. 2020, Icarus, 338, 113546
  * Yelle et al. (2004) Yelle, R. V., Soderblom, L. A., & Jokipii, J. R. 2004, Icarus, 167, 30
  * Zakharov et al. (2018a) Zakharov, V. V., Crifo, J. F., Rodionov, A. V., Rubin, M., & Altwegg, K. 2018a, A&A, 618, A71
  * Zakharov et al. (2018b) Zakharov, V. V., Ivanovski, S. L., Crifo, J. F., et al. 2018b, Icarus, 312, 121

##  Appendix A Distribution of physical properties

Figure 16: Distribution of physical properties of the aggregates for STP063.
Figure 17: Distribution of physical properties of the aggregates for STP070.
Figure 18: Distribution of physical properties of the aggregates for STP086.
Figure 19: Distribution of physical properties of the aggregates for STP096.

##  Appendix B Source regions

Figure 20: 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). Cyan solid line
indicates the sub-solar latitude. Figure 21: Ejection efficiency for image set
STP070. Figure 22: Ejection efficiency for image set STP086. Figure 23:
Ejection efficiency for image set STP092. Figure 24: Ejection efficiency for
image set STP096.

##  Appendix C Corrected number of tracks per image

Figure 25: Corrected number of tracks per image. From top left are represented
STP070, STP086, STP092 and STP096 respectively.