© ESO 2018

1. Introduction

Minor bodies, especially comets, are invaluable sources of information that allow us to better understand how the solar system formed. They are considered to be time capsules and planetary building blocks, and are the oldest and least-evolved bodies left over from the primitive proto-solar system. They therefore represent the earliest record of material from this epoch. In addition, they are the most organic-rich bodies in the proto-solar system, and the fully formed molecules contained within their nuclei could have played a key role in the origin of life on Earth. Moreover, it is thought that minor bodies also played an important role in the hydration process of early Earth; see e.g. Hartogh et al. (2011) and Jewitt et al. (2007). Our understanding and knowledge of comets have been revolutionised by the outstanding Rosetta mission (Taylor et al. 2015) to the Jupiter family comet (JFC) 67P/Churyumov–Gerasimenko (hereafter 67P). Snodgrass et al. (2017) concluded from the observing campaign of 67P that Rosetta was seeing a typical JFC object. This allows the conclusions from Rosetta measurements to be taken as generally true for JFCs. During April 2017, the close encounter between 41P/Tuttle–Giacobini–Kresak (hereafter 41P) and Earth offered a unique opportunity to apply the lessons learned from Rosetta to other JFCs. To date, comet 41P has been poorly studied, where the only published literature was the work of Kresak (1974) who reported on two large outbursts suffered by the comet in 1973. However, 41P came to international attention because of its surprisingly fast rotational variation reported by Bodewits et al. (2018) as gleaned from observations obtained from March to May, 2017. As well as being discovered in 1858 by Horace Parnell Tuttle at Cambridge, 41P was also independently identified by Michael Giacobini in 1907 in Nice. However, it was not until 1951 when Lubor Kresák performed a third, independent identification at Skalnaté Pleso, that enough information was obtained for a proper orbit characterisation. The comet was classified as a member of the JFCs, which are assumed to come from the trans-Neptunian region, where the comets are dynamically controlled by Neptune, and in some cases are injected inwards passing to the status of Centaurs, the direct progenitors of JFCs. When Centaurs finally fall under the gravitational control of Jupiter, they become JFCs (Duncan et al. 1988; Levison & Duncan 1997). On the other hand, it has been suggested that JFCs could have other source regions closer to the Sun, such as the main asteroid belt (Fernández & Sosa 2015), the Hilda family in the 3:2 Mean Motion Resonance with Jupiter Di Sisto et al. 2005; Toth 2006, or, with an extremely low contribution, Jupiter’s Trojans (Volk & Malhotra 2008).

In this paper we present the results of a dust analysis using observations obtained with a large monitoring campaign carried out with the TRAPPIST telescopes (Jehin et al. 2011). In a second paper, an analysis of the gas content of the coma will be presented by Moulane et al. (in prep.). The observations were obtained between January and July, 2017, which encompasses both pre- and post-perihelion epochs. In the first part of this paper, we determine the evolution of the dust environment of 41P throughout the time the comet was active using Monte Carlo dust tail simulations (Moreno et al. 2012). This provides us information such as the dust production rate, the size distribution and ejection velocities of the dust particles, and the emission pattern. We compare our results with those obtained by other authors for similar comets, especially with those obtained from the Rosetta mission; for example, Rotundi et al. (2015), Fulle et al. (2016b), and Moreno et al. (2017a). In the second part of the study, we analyse the comet’s orbital stability using numerical simulations. We also characterize its dynamical nature, which allows us to constrain its origin and future evolution. By merging our dust analysis results with its dynamical characteristics, we have developed a better understanding of 41P’s current and future behaviour.

2. Observations and data reduction

We performed long-term, high-cadence monitoring of comet 41P using the TRAPPIST network, that is, TRAPPIST-South at ESO-La Silla Observatory in Chile and TRAPPIST-North at Oukaimeden Observatory in Morocco at several epochs pre- and post-perihelion. Both are 60 cm Ritchey-Chretien telescopes, which have thermoelectrically cooled detectors: a 2K × 2K FLI Proline CCD camera with a field of view of 22′ × 22′ at TRAPPIST-South (Jehin et al. 2011) and an Andor IKONL BEX2DD CCD camera at TRAPPIST-North. For our dust modelling purposes, we used the broad-band R Johnson-Cousins filter to minimise gas contamination due to intense emission bands located in the UV and the blue portion of the comet’s spectrum.

Our observational campaign started on January 20, 2017 (r_h = 1.49 au, inbound), when the comet was bright enough to be observed with TRAPPIST-South. We monitored the comet, mostly with TRAPPIST-North, until it was too faint to be detected on July 27, 2017 (r_h = 1.69 au, outbound). During that period, we obtained more than 40 nights of observations. The final data set is composed of 30 photometric nights, as we discarded observations obtained on cloudy nights. Our observational log is shown in Table 1.

In order to improve the signal-to-noise ratio, the comet was imaged several times each night using integration times in the range 60–120 s. The individual images were flat-fielded and bias subtracted using standard techniques, then a median stack was obtained from the available images. The flux calibration was done using the USNO-B1.0 star catalogue (Monet et al. 2003), so that each image we acquired was calibrated in mag arcsec⁻², and then converted to Solar Disk Units (SDUs). We oriented each image so that north is up and east is to the left. When necessary, the images were rebinned to have physical dimensions small enough as to be analysed with the Monte Carlo dust tail code. The complete set of observations is displayed in Fig. 1.

Fig. 1. Complete set of observations using both TRAPPIST-South and TRAPPIST-North telescopes at La Silla Observatory and Oukaimeden Observatory, respectively. Isophote levels in Solar Disk Units (SDUs) are displayed for each row on the right-hand side. In all cases, north is up and east to the left. The observational conditions of each image is given in Table 1. Images 1, 7, 10, 11, 18, 21, 24, 25, 28, and 30 are marked with †, which are used for comparison with the model in Fig. 8.

3. Dust analysis

3.1. The Monte Carlo dust tail model

The dust analysis was performed using the Monte Carlo dust tail code described in Moreno et al. (2012). The code computes the brightness of a cometary tail by generating synthetic images that can be directly compared with observations. The code has been extensively used to characterise the dust environments of comets as a function of the heliocentric distance; see, for example, Moreno et al. (2014b), Pozuelos et al. (2014b, 2015), and activated asteroids; see, for example, Moreno et al. (2014a, 2016b, 2017b).

The model computes the trajectory and the scattering of a large number of particles ejected from the object’s surface. The particles, after being ejected from the nucleus, are subjected to the gravitational force of the Sun and radiation pressure. The gravity of the comet itself is neglected, which is a valid approximation for small size objects. Any gas molecules, arising from sublimated ice, drag and accelerate the dust particles to their terminal velocities. These velocities and their physical properties are responsible for the final Keplerian motion around the Sun. The ratio of the radiation pressure to the gravity exerted on each particle is given by the β parameter, which is expressed as (Finson & Probstein 1968):$$\beta =\frac{C_{\text{pr}}{Q}_{\text{pr}}}{2\rho r},$$(1)

in this equation, C_pr is given by:$$C_{\text{pr}}=\frac{{3E}_{\odot }}{{8\pi cGM}_{\odot }},$$(2)

where E_⊙ is the mean solar radiation, c is the light speed, G is the gravitational constant, and M_⊙ is the solar mass. Q_pr is the radiation pressure coefficient, which is ~1 for particles of radius r ≥ 1 μm (Moreno et al. 2012, their Fig. 5), ρ is the particle density, and r is the radius.

Since the model has many parameters, a number of assumptions were made to make the problem more tractable. With this aim, we followed the recent application of the Monte Carlo dust model for 67P by Moreno et al. (2017a), where the authors developed a detailed dust analysis based on their large groundbased observational data set and a useful compilation of the results of the Rosetta mission regarding the properties of its dust particles. Here we briefly summarise these assumptions, and we refer the reader to Moreno et al. (2017a) and the references therein for further information.

From the Grain Impact Analyser and Dust Accumulator (GIADA; Colangeli et al. 2007) and the Optical, Spectroscopic, and Infrared Remote Imaging System (OSIRIS; Keller et al. 2007) measurements, the density of the particles and the geometric albedo were set to values of ρ = 800 kg m⁻³ and p_v = 0.065, respectively (Fulle et al. 2016b; Fornasier et al. 2015). From the Micro-Imaging Dust Analysis System (MIDAS; Riedler et al. 2007), we assumed that the minimum size of the particles is 10 μm, which is time-independent. On the other hand, the maximum size of particles is considered time-dependent ranging from 1 cm at large heliocentric distances (Rotundi et al. 2015) to decimeter-sized aggregates during the perihelion passage (Fulle et al. 2016a). The size distribution was assumed to be a power-law function, n(r) ∝ r^δ(t), where the timedependent parameter δ(t) is set to vary between −4.2 and −2; see, for example, Rotundi et al. (2015), Fulle et al. (2016a), and Ott et al. (2017). The terminal velocities of the particles are parametrised by:$$\upsilon (t,\beta )={\upsilon }_0\left(t\right)\times {\beta }^{\gamma },$$(3)

where γ = 1/2. This assumption is commonly accepted for hydrodynamical drag from sublimating ice; see, for example, Moreno et al. (2011) and Licandro et al. (2013), which also agrees with the value range of 0.42–1.5 reported from GIADA measurements (Della Corte et al. 2015, 2016). On the other hand, ν₀(t) is a time-dependent term which is determined during the modelling process.

In addition to ν₀(t), the other time-dependent parameters in the model are the dust mass loss rate, Q_dust(t), the maximum size of the particles, r_max(t), and the power-law index of the particle size distribution δ(t). Due to the large number of free parameters used in the model, the obtained solution is not always unique, and it may be possible to find an alternative set of parameter values that also fit the observational data. However, this indetermination is reduced considerably when the available observations cover a significant orbital arc. The modelling method consists of a trial-and-error procedure starting from the simplest scenario, i.e., symmetric behaviour of the time-dependent parameters with respect to perihelion, assuming an isotropic ejection pattern. From this starting point, we subsequently varied the parameter values. If after many trials a good match with the observations is not found, we switch to anisotropic ejection pattern, where the emission of the particles is characterised by active areas on the comet’s surface and the rotational state is defined by by two angles (Sekanina 1981): the obliquity of the orbit plane to the equator, I, and the argument of the subsolar meridian at perihelion, ϕ. The obliquity determines the sense of the rotation, which is prograde when 0° ≤ I < 90° and retrograde when 90° < I ≤ 180°. On the other hand, when 0° < ϕ < 180°, the northern pole experiences sunlight at perihelion; the southern pole when 180° < ϕ < 360°. However, there is an ambiguity in the formalism, and equivalent solutions are given by 180° −I and ϕ ± 180°, thus the sense of the rotation remains indeterminate.

To determine which is the best model, we followed the method of Moreno et al. (2016a), who computed the goodness-of-fit through the quantity χ = Σσ_i, where the summation extends to all images under consideration, and:$${\sigma }_i=\sqrt{\mathrm{\Sigma }{[\text{log}\left(I_{\text{obs}}\left(i\right)\right)-\text{log}\left(I_{\text{mod}}\left(i\right)\right)]}^2/N\left(i\right)},$$(4)

where N(i) is the number of pixels in the image i, and I_obs(i) and I_mod(i) are the observed and modelled brightnesses, respectively. For every trial, we calculated σ for each image looking for the minimum χ. The choice of work in the logarithm of the intensities instead of the original intensities, allow us to give more appropriate weights to the outermost and innermost isophotes.

Due to a lack of information about the nucleus itself, we had to assume its bulk density, ρ_N, and size, R_N. Pätzold et al. (2016) estimated that the average bulk density of the nucleus of 67P was ρ = 533 ± 6 kg m⁻³, based on Rosetta measurements. From the Deep Impact mission, A’Hearn et al. (2005) estimated the nucleus density of comet Tempel 1 to be ρ_N ~ 600 kg m⁻³. Since comet 41P belongs to the same family of comets, some properties are expected to be roughly similar; therefore, we assumed an intermediate value of ρ_N ~ 550 kg m⁻³. On the other hand, Tancredi et al. (2000, 2006) estimated the size of 41P’s nucleus to be R_N ~ 0.7 km, and Howell et al. (2017) reported a minimum size of R_N ~ 0.9 km from radar observations. We adopted the latter value in this study.

3.2. Results and discussion

We found that we could not explain the observational data set using a full isotropic ejection model, and it was necessary to switch to a more complex scenario. In fact, we find that the emission pattern experienced two transformations: from full isotropic to anisotropic dominated by two strongly active areas, and then from anisotropic to full isotropic again. We named this the hybrid model. A comparison between full isotropic emission and the hybrid model is shown in Fig. 2.

Fig. 2. Comparison between the full isotropic model (panel a), and the hybrid model (panel b). In both cases, the red contours correspond to the models and the black ones to the observations. The observation date is March 14, 2017. For modelling purposes, the image was rebinned ×2 with respect to the value given in Table 1. Therefore, the resolution is 612.6 km pixel−1. In all cases the isophote levels are: 4 × 10−14, 7 × 10−14, 1.25 × 10−13 and 4 × 10−13 SDU. The validity of the models, determined via the χ parameter, confirms that the isotropic model offers a poorer fit (χisotropic = 3.9) than the hybrid model (χhybrid = 2.4). The plot is orientated so that north is up and east is to the left.

Chronologically, the comet started with a full isotropic ejection pattern, when the level of activity was very low and it had dust production rates of Q_dust = (3−20) kg s⁻¹. When the comet was ~(1.218–1.116) au inbound (February 24–March 14, 2017), the emission started to switch to anisotropic. Two active areas took over and dominated the emission of the dust particles, ejecting ~90% of them. The latitudes of these active areas were found to be located both in the northern hemisphere (from (45 ± 10)° to $\left({90}_{-10}^{+0}\right)$°) and in the southern hemisphere (from (−35 ± 10)° to $(-{90}_{-10}^{+0})$°). The combination of these areas represents about 23.7% and 51.2% of the total cometary nucleus surface, respectively. In principle, from our models, we are also able to constrain the longitude of the active areas. However, a well defined area in both longitude and latitude on the surface of a rotating nucleus, after many rotations, will produce the same effect as if the nucleus was ejecting particles across the entire longitude range (from 0° to 360°). Therefore, we were only able to properly constrain the latitudes of the active areas, and we are missing information about their longitudes.

This ejection pattern lasted until ~(1.20–1.45) au outbound (June 7–28, 2017), when the activity decreased and switched again to full isotropic. During the anisotropic period, the level of activity reached a maximum a few days before perihelion, and it was possible to characterise the rotational state of the comet via rotational angles ϕ and I. The subsolar meridian at perihelion was found to be ϕ = (5 ± 3)°, and the obliquity I = (25 ± 15)°. A schematic vision of the hybrid model is displayed in Fig. 3.

Fig. 3. Schematic vision of the evolution of the dust emission pattern of 41P along its orbit. The four inner planets in the solar system are also included. The configuration corresponds to the time of our last observation, that is, July 27, 2017. The two pairs of red crosses point out the changes in the emission pattern: first is from full isotropic to anisotropic during the inbound journey (February 24–March 14, 2017) and the second is from anisotropic to full isotropic during the outbound journey (June 7–28, 2017). The estimated start date of the activity (October 14, 2016) is marked with green dots, and perihelion (April 12, 2017) is marked with an orange triangle.

From aperture photometry performed on CN narrowband imaging, Bodewits et al. (2018) found a very rapid change in the apparent periodicity of the jets from 20 to 50 h between March and May 2017. In Moulane et al. (in prep.) we confirm this behaviour through TRAPPIST measurements by the comparison of coma features exhibited by the CN gas species between March and April 2017. Rotational period changes are not uncommon in comets, however, it is expected that such changes take place over a time period of a number of years, being observable from one orbit to another (Samarasinha et al. 2004). Bodewits et al. (2018) found two jets associated with 41P, which were proposed as the main cause of the spin down.

In this context, our dust model with two strongly active areas seems to match well with their results. In order to explain the spin down, these powerful dust jets, which imply strong outgassing, should be located at different longitudes in the north and in the south, in such a way that when combined with the small size of the nucleus, they exert torques that act as brakes due to reactional forces. However, as explained before, we are unable to properly constrain the longitude of the active areas and hence unambiguously confirm this assertion.

On the other hand, in order to provoke a spin down in a cometary nucleus, the most efficient place would be at the equator. Unfortunately, all of our attempts to place an active area at the equator were unsuccessful. In addition, the local topography and/or the real physical shape of the cometary nucleus could also greatly affect the final net torque. However, the lack of information about the nucleus itself prevents us from confirming this hypothesis, and leads us to consider that other unknown factors could also be affecting the final torque. Moreover, since the rotational periods found by Bodewits et al. (2018) are much shorter than the age of the dust tail, it is not possible to confirm these values with our dust model.

The evolution of the dust parameters that best describe the dust environment along its orbit are displayed in Fig. 4, the dust production rate; Fig. 6, power index of the size distribution and maximum size of particles; and Fig. 7, ejection velocities of the particles for several sizes.

Fig. 4. Dust production rate (solid black line) and its error (red border) given by the hybrid model as a function of day relative to perihelion (lower x-axis) and heliocentric distance (upper x-axis). The water production rates computed from OH observations with TRAPPIST telescopes are shown as blue dots (Moulane et al., in prep.).

Overall, we found that the activity started about 180 ± 10 days before perihelion (October 14, 2016), which corresponds to 2.31 ± 0.08 au inbound, that is, between October 4 and October 24, 2016. The activity increased gently until −1.64 au, whereby it began to increase to a faster rate until the peak of the activity at ~17 days pre-perihelion with a dust production rate of Q_dust = 110 kg s⁻¹. The total dust ejected from the estimated activation date to our last observation is ~7.5 × 10⁸ kg. Since our observations cover a range that spans the beginning of the activity until when the comet nearly switched off, we can consider that this is a good estimate of the total dust ejected during the complete orbit. This value is significantly lower than previously reported for other comets of the same family also considering a complete orbit; see, for example, 67P, 1.4 × 10¹⁰ kg given by Moreno et al. (2017a), 81P/Wild 2 (hereafter 81P), 1.1 × 10¹⁰ kg given by Pozuelos et al. (2014a) or 22P/Kopff, with a total of 8 × 10⁹ kg estimated by Moreno et al. (2012). Comet 41P is therefore a dust-poor comet (Fig. 4). However, one has to keep in mind that 41P is remarkably smaller in size than these comets, where its total mass is estimated to be M ~ 1.6 × 10¹² kg from the current data available and the hypotheses made in this study regarding its nucleus density and spherical shape. This means a small but non-negligible amount of erosion occurs per perihelion passage. The dust-to-gas production rate ratio provides important clues about the formation mechanism and evolution of comets. Using the TRAPPIST telescopes, we also carried out observations of the main gases: OH, CN, NH, C₂, and C₃. Those observations and the gases’ production rates will be presented in a separate paper (Moulane et al., in prep.). In particular, from the production rates of OH, we computed the water production rates using the relation:$$Q_{{\mathrm{H}}_2\mathrm{O}}=1.361\times r^{-0.5}\times Q_{\text{OH}},$$(5)

given by Cochran & Schleicher (1993). In Fig. 4, both the dust production and water production rates are displayed, while in Fig. 5 the evolution of the dust-to-water mass ratio is shown, which varies from 0.5 at ~(1.3–1.10) au inbound to 1.5 when the peak of activity reached a maximum a few days before perihelion at 1.07 au. Post-perihelion, the dust-to-water mass ratio stayed at values close to 1 until the last observation available, when it dropped to 0.25. This low dust-to-water mass ratio implies, and confirms, the aforementioned assertion that 41P is a dust-poor comet, and its nucleus is richer in volatiles. These values are notably smaller than those obtained for 67P, where the dust-to-water mass ratios ranged from 6 to 100 at perihelion and immediately afterwards from in situ measurements obtained with the GIADA and OSIRIS instruments inferred by individual particle detections (Fulle et al. 2016a), and 4 ± 2 in Rotundi et al. (2015). Comets with low dust-to-gas mass ratios are more prone to outbursts, which is attributed to a mechanism based on the phase transitions of volatile species between solid, liquid, and gaseous states. Moreover, high mixing ratios of volatile ices may indicate low physical strengths, where dust fails to act as an efficient binder, thus provoking an unstable nucleus prone to fragmentation events (Miles 2016). An example of this is comet 73P/Schwassmann-Wachmann, which split into many fragments in 2006 (Ishiguro et al. 2009); its dust-to-gas ratio was 0.12 (Lisse et al. 2002). 41P has suffered several outbursts in its recent past: two violent outbursts in 1973, with a brightness amplitude of nine magnitudes (Kresak 1974) and two other outbursts in 2001, with a brightness amplitude of four magnitudes (Bodewits et al. 2018). Due to the dust erosion per orbit found in our model for such a small nucleus, the low dust-to-gas ratio and its propensity to outbursts, 41P is a good candidate for a split event in the future.

Fig. 5. Dust-to-water mass ratio as a function of day relative to perihelion (lower x-axis) and heliocentric distance (upper x-axis).

The power-law index of the size distribution at the beginning of the activity was δ = −3.50. It reached a minimum around perihelion as δ = −3.75 and finished as δ = −3.60 (Fig. 6, panel a). This behaviour of having a minimum value at perihelion was reported previously in other comets; see, for example, 67P in Moreno et al. (2017a), and 81P and 103P/Hartley 2 (103P) in Pozuelos et al. (2014a), which indicates that the population of small particles ejected increase and dominate the size distribution and the observed brightness. The minimum size of particles was a constant in the model presented here (set to 10 μm), and it is these particles that were the fastest, reaching ν_ejec ~200 m s⁻¹ at perihelion. On the other hand, the slowest velocity in the model corresponds to the largest particles, which increase from 1 cm at the start of the activity to 40 cm during the peak of the activity, and decrease to 2 cm at the end. We called this velocity V_min (Fig. 6, panel b and Fig. 7). The maximum size of the particles found during perihelion agrees very well with the results from in situ missions such as Rosetta, where the authors found particles up to 40 cm at perihelion distances (Rotundi et al. 2015). Values of 20 cm were reported by Harmon et al. (2011) during the fly-by of EPOXI mission (A’Hearn et al. 2011; Meech et al. 2011) around 103P. The minimum velocity in the model, V_min, has to fulfil the condition V_min ≳ V_esc, where V_esc is the escape velocity of the nucleus given by$$V_{\text{esc}}=R_{\mathrm{N}}\sqrt{\frac{2}{15}\pi \rho G},$$(6)

Fig. 6. Left panel: evolution of the power-law index of the particle-size distribution, δ(t). Right panel: evolution of the maximum size of the ejected particles, rmax(t). Both are presented as functions of day relative to perihelion (lower x-axis) and heliocentric distance (upper x-axis).

Fig. 7. Ejection-velocity field of the hybrid model, as a function of day relative to perihelion (lower x-axis) and heliocentric distance (upper x-axis). The velocity is parametrised in the model by Eq. (3). The fastest particles correspond to a size of 10 μm. In addition, the velocities of the 100 μm, 0.1 cm, and 1 cm sized particles are shown. The slowest velocity in the model is red-labelled as Vmin, which corresponds to the largest particles (from 1 cm to 40 cm; see Fig. 6, right panel).

for spherical shape, at distance ~20 R_N where the gas drag vanishes. We find V_esc ~ 0.11 m s⁻¹, the condition being overcome all the time. The velocities reported by Della Corte et al. (2016) from GIADA measurements, adopted by Moreno et al. (2017a, their Fig. 2), showed particles of 100 μm at ~1.5 au with velocities of ~22 m s⁻¹. In this study we derive that, at the same heliocentric distance and for the same size of particles as 67P, 41P ejected particles at about ~10 m s⁻¹ (see Fig. 7). This is consistent with the fact that 41P is smaller in size, and it has both weaker dust production and gas production rates than 67P. This can be concluded by comparing the values obtained for 67P by Opitom et al. (2017) and the values obtained for 41P by Moulane et al. (in prep.). Both studies used TRAPPIST telescopes.

A comparison of ten selected images from the observational data set (panels marked with a † in Fig. 1) with the corresponding synthetic images generated by the model are shown in Fig. 8. In addition to the images, our model was forced to match the observationally derived A(θ)fρ parameter (A’Hearn et al. 1984). This quantity is related to the dust coma brightness, where A(θ) is the geometric albedo as a function of phase angle (θ), f is the filling factor in the aperture of field of view, and ρ is the projected distance from the nucleus. We computed A(θ)fρ for our complete set of observations at ρ = 10⁴ km, which allowed us to track its time evolution.

Fig. 8. Comparison of a subset of ten observed and modelled images. The selected observed images are marked with a † in Fig. 1; here we keep the same numeration. For modelling purposes, those images labelled with an “*” were rebinned ×2 with respect to the values given in Table 1. The isophote levels in SDU in each case are: (1†) 0.5 × 10−14, 1.0 × 10−14, 2.0 × 10−14 and 6.0 × 10−14. (7†,*) 1.7 × 10−14, 3.0 × 10−14, 7.0 × 10−14 and 2.0 × 10−13. (10†,*) 4.0 × 10−14, 7.0 × 10−14, 1.2 × 10−13 and 4.0 × 10−13. (11†,*) 4.0 × 10−14, 7.0 × 10−14, 1.2 × 10−14 and 4.0 × 10−13. (18†,*) 4.0 × 10−14, 7.0 × 10−14, 1.2 × 10−14 and 4.0 × 10−13. (21†,*) 4.0 × 10−14, 7.0 × 10−14, 1.2 × 10−14 and 4.0 × 10−13. (24†,*) 4.0 × 10−14, 7.0 × 10−14, 1.2 × 10−14 and 4.0 × 10−13. (25†) 4.0 × 10−14, 7.0 × 10−14, 1.2 × 10−14 and 4.0 × 10−13. (30†) 1.7 × 10−14, 3.5 × 10−14 and 9.0 × 10−14. In all cases, the black contours correspond to observations and the red ones to the model. The black vertical lines correspond to 2 × 104 km in the sky plane, and the y and x axes are given in pixels.

Since the phase angle varied significantly from 0° to 70° (see Table 2), it was necessary to adopt a correction to distinguish between enhancements due to the actual behaviour of the comet and those arising from phase angle effects, such as backscattering at low phase angles. This effect has been observed for a number of comets. In particular, for 67P it was observed through ground-based observations made by Moreno et al. (2017a), where the authors reported two enhancements at low phase angles. The first of them occurred 400 days pre-perihelion, when it was almost null because of the low level of activity of the comet at that moment. The second peak occurred 200 days post-perihelion, when the level of activity was higher and the enhancement was better resolved. Moreover, this effect was confirmed by Bertini et al. (2017), where the authors studied the scattering phase function at a wide range of phase angles both pre- and post-perihelion in situ using the OSIRIS instrument on board Rosetta. However, in our A(θ)fρ measurements the backscattering effect does not seem appreciable. There may be several reasons for this behaviour: it may be due to the very low dust production rate at the time of the minimum, which occurred at the beginning of the activity, and therefore the total amount of dust was very low and the precision of the data are not good enough. It may also be due to the low number of observations available around the minimum. These reasons suggest that the expected backscattering enhancement was not detectable in our dataset. Moreover, it may have been hidden because of the rapid increase of cometary activity a short while afterwards. Notwithstanding, we adopted a correction following the combined phase function computed by D. Schleicher¹ from observations of comets at all phase angles.

The values of A(θ)fρ directly computed from the observations, and the conversion to A(θ = 0°)fρ, are compared with the synthetic values obtained from the hybrid model and displayed in Fig. 9. In general terms, the hybrid model proposed here fits very well both images and the A(θ = 0°)fρ parameter.

Fig. 9. Upper panel: A(θ)fρ values computed directly from the observational data set (left y-axis) and phase angle, θ (right y-axis). Lower panel: A(θ = 0°)fρ observational values and the synthetic A(θ = 0°)fρ values obtained from the hybrid model. Both panels show the data as a function of day relative to perihelion (lower x-axis) and heliocentric distance (upper x-axis).

4. Dynamical analysis

The main reason for studying the dynamical evolution of comet 41P is to understand how long this comet has been suffering its current rate of erosion, and for how long it will continue. Comet 41P is a near-Earth Jupiter family comet (NEJFC), that is, a JFC with a perihelion distance of q < 1.3 au. The recent work by Fernández & Sosa (2015) revealed a subgroup among NEJFCs that reside in highly stable orbits, with a likely origin in the main asteroid belt. This new class of objects could be the counterparts to the Main Belt Comets (Jewitt et al. 2015), that is, they may be asteroids disguised as comets. In order to clarify the dynamical nature of 41P, we performed numerical integrations like many other authors before; see, for example, Di Sisto et al. (2009), Ye et al. (2016), and Fernández et al. (2017).

4.1. Numerical integrations

In order to study the dynamical evolution of 41P, we used numerical integrations in the heliocentric frame. A quick first inspection consisted of integrations over 2 × 10⁴ yr: from current time to 10⁴ yr backward and 10⁴ yr forward. The current time was set as October 25, 2017. The initial conditions of the orbital parameters were extracted from the NASA/JPL Small-Body Database, and can be consulted in Table 2. The numerical integrations were performed twice using different numerical packages, MERCURY6 (Chambers 1999) and REBOUND (Rein & Liu 2012; Rein & Spiegel 2015) with equivalent conditions. The results obtained were the same; therefore here we only describe the ones obtained with MERCURY6, which allow for a more meaningful comparison with analyses performed by other authors. The Sun and the eight planets were included in the simulations. We used the hybrid algorithm, which combines a second-order mixed-variable symplectic algorithm with a Bulirsch-Stoer integrator to manage the close encounters. The initial time-step was set to 8 days, and the computed orbital evolution was stored every year. We considered as ejected the particles with heliocentric distance r_h > 100 au. Any close encounters with planets at distances smaller than 3× Hill radii were also registered. From the dust analysis, we derived the action of two strongly active areas which seem to be related to the fast rotational period variations. We decided to perform two trials: the first one considered only a pure gravitational model, and the second included non-gravitational forces (see Table 2). In the first experiment, the results were exactly the same, therefore only the results for the pure gravitational model are shown in Fig. 10. From these results, we find that its current state as an NEJFC was obtained in the recent past after a close encounter with Jupiter in the period of time studied (~−1572 yr), at a distance of d_{41P−Jupiter} ~ 0.12 au, which is well below the Hill radius. Due to this event, all the orbital parameters changed dramatically. In addition, it was observed that currently the comet is in a transitional state, and it may reach a more stable orbit with some of its orbital parameters coupled. The perihelion distance would be even smaller in that future scenario, with a value of q ~ 0.8 au. Therefore, we conclude that being a NEJFC seems to be a relatively new status for 41P, and it will remain as part of this population of comets for at least the next 10⁴ yr.

Fig. 10. Orbital evolution of 41P during 2 × 104 yr: from current time to 104 yr backward and 104 yr forward in time. From top to bottom: closest approaches with Jupiter, semimajor axis, eccentricity, inclination, perihelion and aphelion distance, and Tisserand parameter. The initial orbital elements are displayed in Table 2, which were taken from the JPL Small-Body Data Browser.

However, the chaotic nature of the orbits of the minor bodies in the solar system is well known (Levison & Duncan 1994). Therefore, in order to make more robust conclusions, it is necessary to perform statistical studies. With this aim, we carried out the set of tests proposed by Fernández et al. (2014) and Fernández & Sosa (2015) to obtain the instability degree of comet 41P.

4.2. Characterisation of the degree of instability of 41P’s orbit

While typical JFCs have unstable orbits, probably coming from trans-Neptunian regions, a small group of them reside in asteroidal orbits, and may even originate from the asteroid belt, like near-Earth asteroids. Fernández & Sosa (2015) hypothesised that some NEJFCs could be interlopers from the asteroid belt. These authors analysed a sample of 58 NEJFCs using numerical simulations of their orbits in order to characterise their stability, and they performed statistical studies using clones of those NEJFCs. They found different degrees of stability, which allowed them to distinguish three different categories: Highly asteroidal, whose orbits are highly stable with a possible source region being the outer main asteroid belt; Moderately asteroidal, despite generally having highly stable orbits, a small fraction of the clones were unstable, and therefore gave less confidence to their origin; and Maybe asteroidal, predominantly displaying stable orbits, but an important fraction of the clones had highly unstable orbits due to the occurrence of encounters with Jupiter at distances of less than 0.2 au, more in consonance with the evolution of JFCs. Therefore, a possible asteroidal origin is more uncertain.

Fernández & Sosa (2015) included 41P in their sample. However, they did not report any peculiar result regarding it; therefore, we infer that they did not find anything suspicious to suggest that this comet originated in the asteroid belt. However, in the last perihelion passage in April 2017, the extremely favourable conditions for its observation attracted international attention, increasing the number of observations available, and improving the quality of its determined orbital parameters. Its current condition code is 2². Therefore, we decided to reanalyse its orbital stability following the same steps given by Fernández & Sosa (2015), where a likely dynamical path is defined as the average of the set of results obtained for a given object and its clones, characterised by the f_q index, f_a index, the capture time (here after t_cap), and the closest approach to Jupiter, d_min. Here, we briefly describe these parameters, and we refer the reader to Fernández & Sosa (2015) and Fernández et al. (2014), and references therein, for further information.

The f_q index is computed as:$$f_q=\frac{{\mathrm{\Sigma }}_{j=1}^{N+1}{\mathrm{\Delta }t}_j}{(N+1)\times {10}^4},$$(7)

where Δt_j is the fraction of time in the last 10⁴ yr in which a given JFC or its clones ( j = 1, … , N) describe an orbit with q > 2.5 au or reaches heliocentric distances r_h > 100 au.

The f_a index is computed as:$$f_q=\frac{{\mathrm{\Sigma }}_{j=1}^{N+1}{\mathrm{\Delta }t}_j^{\prime }}{(N+1)\times {10}^4},$$(8)

where Δt′_j is the fraction of time in the last 10⁴ yr in which a given JFC or its clones ( j = 1, … , N) move along an orbit with semi-major axis a > 7.37 au, i.e., an orbital period >20 yr. When reaching these criteria, an object is no longer considered to be a JFC.

In general terms, comets in unstable orbits have f_q and f_a values well above zero, which means that they spend an important fraction of the 10⁴ yr with q > 2.5 and/or a > 7.37 au. On the other hand, when f_q ~ f_a ~ 0, the comets move in stable orbits.

The t_cap parameter is defined as the time in the past at which the average perihelion distance of a given NEJFC and its clones at a certain time, that is,$$\overline{q}\left(t\right)=\frac{{\mathrm{\Sigma }}_{j=1}^{N+1}\mathrm{\Delta }{q}_j\left(t\right)}{(N+1)},$$(9)

increased by 1 au with respect to the observed value at the discovery time, t_disc, namely:$$q\left(t_1\right)=q\left(t_{\text{disc}}\right)+1\Rightarrow t_{\text{cap}}=t_{\text{disc}}-t_1\mathrm{.}$$(10)

An increase of 1 au in q means that the comet is q ~ 2 au farther from the Earth’s vicinity. Thus, the concept of t_cap is related to the time span during which the comet has been in the Earth’s region, i.e. making it a NEJFC. For instance, typical JFCs have t_cap ranging from a few years to a few times 10³ years. On the other hand, the NEJFCs in the category defined by Fernández & Sosa (2015) as Highly asteroidal show values for t_cap that largely exceed 10⁴ yr.

4.3. Results and discussion

In order to achieve a balance between computational cost and good coverage of the physical space around 41P, we generated 200 clones of the original 41P comet. The orbital parameters a, e, and i were chosen randomly from a Gaussian distribution, whose mean values and standard deviations were their osculating values and their 3× osculating uncertainties, σ, respectively (see Table 2). With this choice, we ensured that most of our clones (~70%) were in the 3 × σ area around the osculating orbital parameters. The total integration time was 10⁵ yr: from current time to 5 × 10⁴ yr backward and 5 × 10⁴ yr forward. We applied the concepts described above to our simulations with 200 clones for both models (pure gravitational and one that includes non-gravitational forces). We found very similar results for both of them. In the first instance, for a purely gravitational model, we obtained: f_q = 0.025, f_a = 0.007, and for that with non-gravitational forces, f_q = 0.027, f_a = 0.004. The t_cap parameters were computed for both models from the average perihelion in the last 10⁴ yr given by Eq. (9), and they are displayed in Fig. 11. We find for the pure gravitational model t_cap = 1.09 × 10⁴ yr, and t_cap = 1.10 × 10⁴ yr with non-gravitational forces. Finally, the closest approach to Jupiter is respectively found to be d̄_min = 0.19 au (purely gravitational) and d̄_min = 0.20 au (including non-gravitational forces).

Fig. 11. Average perihelion distance as a function of time for both the pure gravitational model (magenta solid line) and the model that includes non-gravitational forces (yellow solid line). The red circles indicate the time when q̄ = q disc + 1, which is necessary to compute the capture time (see text).

These values seem to be different from the ones obtained by Fernández & Sosa (2015), where the authors did not report a stable orbit; we therefore assume that they obtained f_q > 0.2 and f_a > 0.1. This discrepancy could be explained by the different quality of the orbits (now being more accurate) or the integrator algorithm used in both studies. We used a hybrid code that combines a symplectic algorithm with the Bulirsch-Stoer code, while Fernández & Sosa (2015) used only Bulirsch-Stoer code. Another reason for the discrepancy could be how the clones were built; they generated 50 clones using a Gaussian distribution with a standard deviation given by the osculating uncertainties, while we generated 200 clones, using a Gaussian distribution with a standard deviation given by the 3× osculating uncertainties.

Therefore, according to the new results, comet 41P belongs to the Moderately asteroidal category (see Fernández & Sosa 2015, their Table 3). The authors included in this category the following comets: 197P/LINEAR, 207P/NEAT, 209P/LINEAR, 210P/Christensen, 217P/LINEAR, and 317P/WISE. Therefore, comet 41P is the seventh. All of them generally show little activity compared to normal JFCs. This is in contrast with the comets in the Highly asteroidal category, which were reported to have extremely low levels of activity and even stellar-like appearances in some cases. According to Fernández & Sosa (2015), Moderately asteroidal category still could imply a main asteroid belt origin. However, comet 41P has shown a typical cometary activity based on ice sublimation; therefore, an asteroidal belt origin sounds unlikely. On the other hand, our set of experiments confirms that its orbit is more stable than usual for typical JFCs. Most of them are indeed moving in highly unstable orbits, with capture times in their current near-Earth orbits being lower than a few times 10³ yr.

After completing the set of tests proposed by Fernández et al. (2014) and Fernández & Sosa (2015), one wonders about the rough period of time that comet 41P will stay in its current orbit. Using the same sample of 200 clones, we performed a couple of extra experiments to address this question. First, we divided the possible locations of the comet and its clones into four regions depending on their dynamical properties, and we computed the time spent on them, namely:$$a<\frac{a_{\mathrm{s}}}{(1+e)}\Rightarrow \text{JFC} \text{type},$$(11) $$\begin{array}{r}\frac{a_{\mathrm{s}}}{(1+e)}<a_{\mathrm{N}}\\ e<0.8\\ \end{array}\}\Rightarrow \text{Centaur} \text{type},$$(12) $$\begin{array}{r}\frac{a_{\mathrm{s}}}{(1+e)}<a_{\mathrm{N}}\\ e>0.8\\ \end{array}\}\Rightarrow \text{Halley} \text{type},$$(13) $$a>a_N\Rightarrow \text{Trans}-\text{Neptunian} \text{type}\mathrm{.}$$(14)

In all these equations, a_S and a_N are the semi-major axes of Saturn and Neptune, respectively. The results of this experiment are displayed in Fig. 12, where the top panel shows the pure gravitational model and the bottom panel shows the model that includes non-gravitational forces. In this figure, the red area corresponds to JFC types, the clear blue corresponds to Centaur types, the yellow to Halley types and the dark blue to trans-Neptunian types. In particular, we are interested in the amount of time spent in the JFC region, where comets reach a temperature high enough to be periodically active. To extract this information from the figure, we computed the time for different confidence levels (CL), based on the fraction of surviving clones in the JFC region, namely: 100% CL (A); 90% CL (B); 80% CL (C), and 70% CL (D).

Fig. 12. Statistical orbital evolution of comet 41P and its 200 clones over 105 yr: 5 × 104 backward in time and 5 × 104 forward. The top panel corresponds to a pure gravitational model, while in the bottom panel non-gravitational forces are included. In both panels different colours refer to different regions in the solar system, namely: red for those in the Jupiter family region, clear blue are Centaurs, yellow are Halley types, and dark blue are trans-Neptunians. Labels A, B, C, and D mean the time spent in the Jupiter family region at different confidence levels, i.e., % of clones in that region. Therefore, A is 100%, B is 90%, C is 80%, and D is 70%.

Therefore, we find that the time spent by 41P in the JFC region is: 3700 ± 100 yr (A); 18 800 ± 100 yr (B); 34 000 ± 100 yr (C); and 67 300 ± 100 yr (D) for the pure gravitational model. For the model that includes non-gravitational forces, we find: 5500 ± 100 yr (A); 15 000 ± 100 yr (B); 33 100 ± 100 yr (C); and 56 300 ± 100 yr (D). In both models we observe that at the end of the simulations, more than half of the clones still remain in the JFC region, with less than 8% of the initial particles being ejected. In the upper panel, we observe that 100% of the particles became JFCs at −1600 ± 100 yr, which matches the result found in Fig. 10 regarding the incursion to its current region after a strong close encounter with Jupiter. In contrast, when we include non-gravitational forces, this incursion occurred earlier. With a 100% CL we obtain that the comet will be in the JFC region for a time ranging between 3700 and 5500 yr. This confirms the stable orbit found in the previous experiments, but favours the trans-Neptunian region as a more plausible origin.

However, 41P is a special case of a JFC in that it belongs to the NEJFC subgroup, that is, a JFC with a perihelion value of q < 1.3 au. Our definition of a JFC given in Eq. (11) does not differentiate between JFCs and NEJFCs. Therefore, our second experiment consisted of computing the time the comet was an NEJFC, and focused on the period in which the comet remains in the JFC region, which was from −8400 to 12 700 yr, for a 90% CL.

Since from both models we obtained very similar results, we decided to compute the average of them. The result is shown in Fig. 13. In the figure, the perihelion evolution of the complete set of clones is displayed in grey, which includes both models, and the average perihelion evolution, q̄, shown in red. We divided the physical space between the JFC and NEJFC regions. We observe that the dynamical evolution will carry 41P to a minimum perihelion distance of q = 0.82 au 950 yr from now. The evolution of the clones is extremely compact during a period of 750 yr around the current time (from −300 to 450 yr), the equivalent to ~140 orbits. During this period, there is no divergence in the orbits. From the current epoch until 3600 yr, q̄ < 1.3 au, therefore we considered this time as the period in which the comet 41P will be in Earth’s neighbourhood. All these results suggest that the orbit of 41P is more stable than typical JFCs. The time in the NEJFC region will last roughly 3600 yr, which is equivalent to ~665 orbits.

Fig. 13. Perihelion evolution during the time in which comet 41P is in the JFC region with a 90% of CL. Grey orbits correspond with the whole set of clones, which take into account both the pure gravitational model and the model that includes non-gravitational forces. The red line is the average perihelion q̄. The dashed horizontal line differentiates between the JFC and NEJFC regions.

5. Summary and conclusions

In the first part of this work, we describe and analyse an extensive observational data set of images of comet 41P obtained with the TRAPPIST telescopes using dust tail models. Those observations cover the main portion of the orbital arc in which the comet displayed activity, from 1.5 au pre-perihelion to 1.7 au post-perihelion. In our dust models we followed the assumptions and guidelines given by Moreno et al. (2017a) for the study of comet 67P, which was the Rosetta target.

Our main conclusion is that it is not possible to explain the complete set of observations using a full isotropic ejection model. In fact, we find that a complex ejection pattern which switched from full isotropic to anisotropic (February 24–March 14), and then back from anisotropic to full isotropic again on June 7–28 provides the best description of the observations. During the anisotropic period, we find that ~90% of the ejected particles came from two strongly active areas, one located in the northern hemisphere and the other in the southern. This model is in agreement with the recent discovery of the fast rotational period variation reported by Bodewits et al. (2018) from March to May, 2017, in the sense that the two powerful active areas could have acted as brakes, increasing the nucleus rotation period. However, the location found in our model for these active areas prevents us from giving final confirmation, and leads us to the consideration that other factors may be acting and possibly affecting the fast spin-down observed. Further investigations are encouraged.

In general terms, from the dust model we obtain that the total dust mass ejected is ~7.5 × 10⁸ kg. This quantity is roughly the total dust ejected by the comet during the whole orbit. This amount of dust is low compared to other comets of the same family; however, 41P is a small comet, and this quantity represents a non-negligible fraction of its total mass. This implies that 41P suffered a substantial amount of erosion during its last incursion to perihelion. From observations of gases also performed with TRAPPIST telescopes (Moulane et al. in prep.), we find that the dust-to-water mass ratio is low ranging from 0.25 to 1.5.

The complete set of dust parameters, which best describe the evolution of its dust environment, is also reported, which includes the maximum particle size, the power-law index of the size distribution, and the ejection velocity field of the particles. All these results confirm the dust-poor nature of 41P.

In the second part of this work, we explored the dynamical nature of 41P with numerical simulations. We followed the set of tests proposed by Fernández et al. (2014) and Fernández & Sosa (2015) to constrain the degree of instability of its orbit. In those experiments, we always considered two models: a pure gravitational model and a model that included non-gravitational forces. No significant differences between them were found. We obtain that 41P is more stable than typical JFCs, its orbit being in the category Moderately asteroidal, which could imply a main asteroid belt origin. However, the complete set of dynamical experiments performed, and the activity being driven by ice-sublimation, favour the trans-Neptunian origin hypothesis. The status of NEJFC, that is, a JFCs with a perihelion distance of q < 1.3 au, seems to be relatively new for this comet. The expected period of time during which the comet will remain in this region is roughly ~3600 yr. A minimum perihelion distance will be reached in 950 yr, with a value of q ~ 0.8 au.

With the information currently available, we estimate the total mass of the comet to be 1.6 × 10¹² kg. If the dust production rate per orbit of ~7.5 × 10⁸ kg remains constant or similar, every incursion to perihelion during which the comet is in Earth’s vicinity will mean a prolonged period of nucleus erosion, and it may lose up to 30% of its mass within 3600 yr; even more if we also consider the gas production rates. This erosion could be even larger since the perihelion distance will decrease by ~0.2 au over the next 950 yr. This fact, combined with its low dust-to-gas mass ratio, and its propensity to undergo outbursts (Kresak 1974; Bodewits et al. 2018) could provoke the disruption of the nucleus in a relatively short period of time.

Acknowledgments

This work is supported by a Marie Curie CO-FUND fellowship, co-founded by the University of Liège and the European Union. TRAPPIST-South is funded by the Belgian Fund for Scientific Research (Fond National de la Recherche Scientifique, FNRS) under the grant FRFC 2.5.594.09.F, with the participation of the Swiss FNS. TRAPPIST-North is a project funded by the University of Liège, and performed in collaboration with Cadi Ayyad University of Marrakesh. E.J. is a Belgian FNRS Senior Research Associate, and M.G. is an FNRS Research Associate. We thank the anonymous referee for having significantly improved the paper and Fernando Moreno and Damien Hutsemékers for their valuable discussions. Some of the simulations in this paper made use of the REBOUND code which can be downloaded freely at http://github.com/hannorein/rebound. We thank Hanno Rein and Daniel Tamayo for their help using REBOUND.

References

All Tables

All Figures

Fig. 1. Complete set of observations using both TRAPPIST-South and TRAPPIST-North telescopes at La Silla Observatory and Oukaimeden Observatory, respectively. Isophote levels in Solar Disk Units (SDUs) are displayed for each row on the right-hand side. In all cases, north is up and east to the left. The observational conditions of each image is given in Table 1. Images 1, 7, 10, 11, 18, 21, 24, 25, 28, and 30 are marked with †, which are used for comparison with the model in Fig. 8.

In the text

Fig. 2. Comparison between the full isotropic model (panel a), and the hybrid model (panel b). In both cases, the red contours correspond to the models and the black ones to the observations. The observation date is March 14, 2017. For modelling purposes, the image was rebinned ×2 with respect to the value given in Table 1. Therefore, the resolution is 612.6 km pixel−1. In all cases the isophote levels are: 4 × 10−14, 7 × 10−14, 1.25 × 10−13 and 4 × 10−13 SDU. The validity of the models, determined via the χ parameter, confirms that the isotropic model offers a poorer fit (χisotropic = 3.9) than the hybrid model (χhybrid = 2.4). The plot is orientated so that north is up and east is to the left.

In the text

Fig. 3. Schematic vision of the evolution of the dust emission pattern of 41P along its orbit. The four inner planets in the solar system are also included. The configuration corresponds to the time of our last observation, that is, July 27, 2017. The two pairs of red crosses point out the changes in the emission pattern: first is from full isotropic to anisotropic during the inbound journey (February 24–March 14, 2017) and the second is from anisotropic to full isotropic during the outbound journey (June 7–28, 2017). The estimated start date of the activity (October 14, 2016) is marked with green dots, and perihelion (April 12, 2017) is marked with an orange triangle.

In the text

	Fig. 4. Dust production rate (solid black line) and its error (red border) given by the hybrid model as a function of day relative to perihelion (lower x-axis) and heliocentric distance (upper x-axis). The water production rates computed from OH observations with TRAPPIST telescopes are shown as blue dots (Moulane et al., in prep.).
In the text

	Fig. 5. Dust-to-water mass ratio as a function of day relative to perihelion (lower x-axis) and heliocentric distance (upper x-axis).
In the text

	Fig. 6. Left panel: evolution of the power-law index of the particle-size distribution, δ(t). Right panel: evolution of the maximum size of the ejected particles, rmax(t). Both are presented as functions of day relative to perihelion (lower x-axis) and heliocentric distance (upper x-axis).
In the text

Fig. 7. Ejection-velocity field of the hybrid model, as a function of day relative to perihelion (lower x-axis) and heliocentric distance (upper x-axis). The velocity is parametrised in the model by Eq. (3). The fastest particles correspond to a size of 10 μm. In addition, the velocities of the 100 μm, 0.1 cm, and 1 cm sized particles are shown. The slowest velocity in the model is red-labelled as Vmin, which corresponds to the largest particles (from 1 cm to 40 cm; see Fig. 6, right panel).

In the text

Fig. 8. Comparison of a subset of ten observed and modelled images. The selected observed images are marked with a † in Fig. 1; here we keep the same numeration. For modelling purposes, those images labelled with an “*” were rebinned ×2 with respect to the values given in Table 1. The isophote levels in SDU in each case are: (1†) 0.5 × 10−14, 1.0 × 10−14, 2.0 × 10−14 and 6.0 × 10−14. (7†,*) 1.7 × 10−14, 3.0 × 10−14, 7.0 × 10−14 and 2.0 × 10−13. (10†,*) 4.0 × 10−14, 7.0 × 10−14, 1.2 × 10−13 and 4.0 × 10−13. (11†,*) 4.0 × 10−14, 7.0 × 10−14, 1.2 × 10−14 and 4.0 × 10−13. (18†,*) 4.0 × 10−14, 7.0 × 10−14, 1.2 × 10−14 and 4.0 × 10−13. (21†,*) 4.0 × 10−14, 7.0 × 10−14, 1.2 × 10−14 and 4.0 × 10−13. (24†,*) 4.0 × 10−14, 7.0 × 10−14, 1.2 × 10−14 and 4.0 × 10−13. (25†) 4.0 × 10−14, 7.0 × 10−14, 1.2 × 10−14 and 4.0 × 10−13. (30†) 1.7 × 10−14, 3.5 × 10−14 and 9.0 × 10−14. In all cases, the black contours correspond to observations and the red ones to the model. The black vertical lines correspond to 2 × 104 km in the sky plane, and the y and x axes are given in pixels.

In the text

	Fig. 9. Upper panel: A(θ)fρ values computed directly from the observational data set (left y-axis) and phase angle, θ (right y-axis). Lower panel: A(θ = 0°)fρ observational values and the synthetic A(θ = 0°)fρ values obtained from the hybrid model. Both panels show the data as a function of day relative to perihelion (lower x-axis) and heliocentric distance (upper x-axis).
In the text

	Fig. 10. Orbital evolution of 41P during 2 × 104 yr: from current time to 104 yr backward and 104 yr forward in time. From top to bottom: closest approaches with Jupiter, semimajor axis, eccentricity, inclination, perihelion and aphelion distance, and Tisserand parameter. The initial orbital elements are displayed in Table 2, which were taken from the JPL Small-Body Data Browser.
In the text

	Fig. 11. Average perihelion distance as a function of time for both the pure gravitational model (magenta solid line) and the model that includes non-gravitational forces (yellow solid line). The red circles indicate the time when q̄ = q disc + 1, which is necessary to compute the capture time (see text).
In the text

Fig. 12. Statistical orbital evolution of comet 41P and its 200 clones over 105 yr: 5 × 104 backward in time and 5 × 104 forward. The top panel corresponds to a pure gravitational model, while in the bottom panel non-gravitational forces are included. In both panels different colours refer to different regions in the solar system, namely: red for those in the Jupiter family region, clear blue are Centaurs, yellow are Halley types, and dark blue are trans-Neptunians. Labels A, B, C, and D mean the time spent in the Jupiter family region at different confidence levels, i.e., % of clones in that region. Therefore, A is 100%, B is 90%, C is 80%, and D is 70%.

In the text

	Fig. 13. Perihelion evolution during the time in which comet 41P is in the JFC region with a 90% of CL. Grey orbits correspond with the whole set of clones, which take into account both the pure gravitational model and the model that includes non-gravitational forces. The red line is the average perihelion q̄. The dashed horizontal line differentiates between the JFC and NEJFC regions.
In the text