The dust and gas environment of comet 8P/Tuttle


 Comet 8P/Tuttle has been selected as a possible backup target for the Comet Interceptor mission (ESA). This comet was observed intensively during its previous perihelion passage, in 2008 January. From those observations, important information was obtained about the physical properties of the nucleus and coma. This study focuses on the coma of 8P/Tuttle using visible spectra and images to derive gas and dust production rates. The production rates obtained suggest that this comet can be considered as ‘typical’ concerning the C2/CN and C3/CN ratios, although, depending on the criteria adopted, it could be defined as C3 depleted. NH2 production rates suggest an enrichment of this molecule. Visible and infrared images have been analysed using a Monte Carlo dust tail model. At comparatively large heliocentric distances, the coma is characterized by a dust-to-water ratio around or less than 1. Nevertheless, when the comet approaches perihelion, and the subsolar latitude crosses the equator, the coma dust-to-water ratio increases significantly, reaching values larger than six. Such a high dust-to-gas ratio around perihelion suggests that the nucleus of 8P/Tuttle is also ‘typical’ regarding the refractory content, considering the comparatively high values of that magnitude estimated for different comets.

2008, and also from distant observations in the visible (at 7.42 au inbound), Snodgrass, Lowry & Fitzsimmons (2008) derived an equally large size for the nucleus, obtaining a radius of 7.6 km. Weissman, Choi & Lowry (2008) also estimated a large size for the nucleus of 8P, obtaining a radius of approximately 6 km, as previously, from observations at large heliocentric distances, according to the information given by Groussin et al. (2019) The close approach to Earth of comet 8P in 2008 allowed Harmon et al. (2008a,b,c) to characterize the nucleus from radar observations. These observations soon revealed that the nucleus of 8P was formed by two lobes with approximate diameters of 5.6 and 4.4 km (Harmon et al. 2008b). A detailed analysis of those radar observations was later given in Harmon et al. (2010). By fitting radar images with a shape model formed by two prolate spheroids, Harmon et al. (2010) estimated that the two lobes had dimensions of 5.75 x 4.11 km and 4.25 x 3.27 km, leading to an effective radius of 3.1 km if calculated from its projected area. Taking advantage of its close approach to Earth, 8P was also observed with the HST by Lamy et al. (2008). These authors concluded that the visible light curve was also best fitted by considering a nucleus composed of two lobes. Assuming an albedo of 0.04, Lamy et al. (2008) estimated that the nucleus consisted of two spheres with radii 1.2 and 2.8 km, leading to a mean radius of 3.0 km from the equal projected area. This estimate is in agreement with the size derived from radar observations. During the 2007-2008 apparition, comet 8P was also observed in the infrared with Spitzer by Groussin et al. (2019). These authors studied the 24μm infrared emission of 8P considering the slightly different results on the shape (prolate shapes versus spheres for the two lobes) and constraints on the spin axis orientation from Harmon et al. (2010) and Lamy et al. (2008). Groussin et al. (2019, table 3) summarize the differences in the characterization of the nucleus of 8P from the observations by Lamy et al. (2008) and Harmon et al. (2010), including the derived spin axis orientations compatible with each observation. From their detailed and careful analysis, Groussin et al. (2019) concluded that the thermal emission was best reproduced by the HST main results, refining the size of the nucleus, i.e. a bilobate nucleus formed by two spheres with radii 1.1 ± 0.1 and 2.7 ± 0.1 km and a spin axis orientation of (RA, Dec.) = (285 ± 12, +20 ± 5) • . By combining their size estimate with the HST magnitude, Groussin et al. (2019) determined that the R geometric albedo of 8P was 0.042 ± 0.008.
Regarding the spin period of the nucleus, Lamy et al. (2008) and Harmon et al. (2010) obtained a very similar value around 11.4 h. Half that spin period has been favoured by Waniak et al. (2009) (as well as by Woodney, Schleicher & Bair 2008) from their analysis of the repeatability and kinematics of CN shells. Nevertheless, the authors recognized that their analysis was based on a simplified model and that the 11.4-h period could be confirmed from CN structures with more sophisticated modelling.
Volatile production rates were also subject of intensive study (e.g. Biver et al. 2008;Bockelée-Morvan et al. 2008;Böhnhardt et al. 2008;Bonev et al. 2008;Borisov et al. 2008;Crovisier et al. 2008;Lovell & Howell 2008;Barber et al. 2009;Jehin et al. 2009;Lippi et al. 2009;Kobayashi et al. 2010 etc.) during the last closest approach to Earth and near perihelion by using different techniques and telescopes. These studies made it possible to estimate the water production rate around perihelion, being within the range of 1 to 6 × 10 28 mol sec −1 . Likewise, the production rates of different molecules (CO, HCN, C 2 H 2 , CH 4 , C 2 H 6 , CH 3 OH, CN, C 3 , and C 2 ) were also obtained. Observations also suggested that 8P could have a low dust-to-gas ratio (e.g. Schleicher 2007a). Some discrepancies exist concerning the singularity of this comet. On the one hand, Bonev et al. (2008) and Böhnhardt et al. (2008) concluded that the mixing ratios of parent molecules in 8P were unusual (with a strong enrichment of methanol and being depleted in other molecules), suggesting that the two lobes could have different volatile composition. On the other, Kobayashi et al. (2010) obtained that the relative abundances of 8P, even if depleted in some molecules with regard to Halley-type comets, were not atypical when compared to long-period comets, finding no evidence of chemical heterogeneity in their observations. Regarding daughter molecules, a typical composition has been reported (e.g. Schleicher 2007a).
Knowledge of 8P is very valuable for the planning of the future Comet Interceptor (hereinafter CI) mission. CI is an ESA F-class mission expected to be launched in 2028 (Snodgrass & Jones 2019) and aimed at exploring a dynamically new comet. The mission is being designed and will be launched without a specific main target. The spacecraft will travel up to the L2 Lagrangian point where it will wait in hibernation until a suitable target is found. In the event of an appropriate target not being discovered before CI had to leave L2, several comets have been selected as scientific valuable backup targets for the mission (Schwamb et al. 2020). Comet 8P/Tuttle is presently considered as a potential backup target for CI. The CI team encouraged the community to observe those comets and analyse archival data in order to assist its scientific prioritisation (Schwamb et al. 2020), and the safe and scientifically fruitful flyby of the spacecraft if finally one of those comets is selected as target. Among the selected backup comets, 8P is the only one probably coming from the Oort cloud, which makes it particularly interesting.
Comet 8P was again at its perihelion (the last one before CI launching) in late 2021 August. Nevertheless, its observability was poor given the low solar elongation. This unfortunate circumstance makes all the data and information already available on this comet of great importance. All the above studies provide the community with invaluable knowledge on comet 8P, particularly regarding the characteristics of the nucleus and the volatile composition of the coma. Information regarding the refractory component is still scarce. In this work, we focus on the gas and dust characteristics. In the next section, a summary of our observations is presented. Section 3 is devoted to our spectroscopic observations, aimed at deriving gas productions rates in the visible. In Section 4, images in the visible, together with archived images taken in the infrared with Spitzer under the programme ID 40270 (PI O. Groussin), are analysed by using a Monte Carlo model to derive dust characteristics and production rates. In Section 5, the main results on both gas and dust are discussed.

O B S E RVAT I O N S A N D DATA R E D U C T I O N
8P was observed on 2007 December 15 from Calar Alto Observatory (Spain) using the 2.2-m telescope equipped with the CAFOS instrument. On that night, several images through the Cousins Rfilter as well as medium-resolution spectra (see Fig. 1) with grism B200 were acquired. The grism B200 provides us with an observable spectral range between 328 and 843 nm with a wavelength scale of 0.447 nm per pixel and a spectral resolving power of 225 (at 500 nm). The slit of the spectrograph is oriented along the northsouth and east-west directions, as projected on the plane of the sky, giving dust and gas radial profiles in those directions to study the deviations from isotropic gas and dust emission. For absolute calibration, observations of appropriate spectrophotometric standard stars were acquired. For the comet observations, the slit width is 2.5 arcsec whereas the usable selected length is 10.6 arcmin, which provides us with radial profiles along selected directions up to projected cometocentric distances of ∼10 5 km along both directions. The spectrophotometric standard star was observed with a width of 5 arcsec and the same slit length.
During and shortly after its closest approach to Earth, a second run to observe 8P was developed from El Roque de los Muchachos Observatory (Spain) with the ALFOSC instrument attached to the NOT telescope. Unfortunately, either technical problems or nonphotometric conditions prevented us to scientifically use most of the observing time. 8P was successfully observed in this second run on the night of 2008 January 4. Broad-band images were taken through the Bessell R-filter and spectroscopic measurements were obtained with grism #3 with initial and final wavelengths of ∼300 and 670 nm, respectively, a spectral resolution of 0.23 nm per pixel and a spectral resolving power of 190 (at 500.0 nm). The slit, of 2.5 arcsec of width and 6.3 arcmin length, was positioned along the N-S direction. Tracking on the comet was used and corresponding bias, flats, lamp spectra as well as calibration stars were also taken to perform standard calibration.
The spectra were reduced using the ESO-MIDAS standard reduction context long for long-slit spectra. The spectra were bias subtracted, flatfielded, wavelength calibrated, extinction corrected (using the standard extinction curve for both observatory sites), and flux calibrated. Since the spectrum of 8P/Tuttle covers the entire length of the slit, it is not possible to extract the sky contamination directly from the spectrum itself. Instead, this subtraction has been estimated from spectral regions free of gas emissions nearby the CN, C 2 , C 3 , and NH 2 bands. More details on the images and spectra reduction and calibration can be found in Lara et al. (2001Lara et al. ( , 2004a and Bertini et al. (2009). Table 1 summarizes observations as well as comet geometrical circumstances. Our observations were aimed at retrieving production rates of CN, C 2 , C 3 , and NH 2 and at determining dust properties and production rate using Monte Carlo models as described in e.g. Moreno et al. (2016). As our images were close in time, and the performance of Monte Carlo models may benefit from distant observations, we decided to include publicly available data in our study. On the one hand, our Afρ determined from the spectroscopic measurements in clear dust continuum regions were combined with amateur determinations of the Afρ parameter as a function of the heliocentric distance to constrain the dust model, as described in Section 4. On the other hand, considering the recent publication by Groussin et al. (2019) and the suitability of including infrared observations in our Monte Carlo study, images of 8P available at the Spitzer Heritage Archive 2 corresponding to the programme with ID 40270, 3 PI: O. Groussin, were used. In particular, we utilized MIPS level 2 24-μm images 4 taken on 2008 June 22, when the comet was at 2.25 and 1.58 au from the Sun and Spitzer, respectively. Shadow observations, showing consistency with the background sky, were subtracted from the mosaicked images to extract comet flux.

T H E G A S C O M A
The spectra of the comet are used to investigate the CN, C 2 , C 3 , and NH 2 , in the directions probed by the slit. The spectral regions are specified in Table 2. For the subtraction of the underlying continuum in the gas emission bands, we have measured the continuum bordering each emission band (see Table 2) and approximated the continuum contribution to the band by interpolating the left-and right-hand continua (Lara et al. 2001). The conversion of the emission band fluxes into column densities makes use of constant g-factors for C 2 , C 3 , and NH 2 (Cochran et al. 1992, A'Hearn et al. 1995, Schulz et al. 1998, whereas the g-factor of the CN molecule is calculated for the heliocentric distance and velocity of 8P on every date from the set of values given by Schleicher (2010). Considering all the procedures involved, our estimate of the error in column densities determination is less than 5 per cent, at least up to moderate cometocentric distances (∼5000 km). The outflow velocity of the gas adopted in this study is v = 0.85 r h −0.5 (e.g. Cochran, Barker & Gray 2012), where r h is the heliocentric distance.
The gas profiles are asymmetric in the directions probed by the slit, indicating a non-isotropic gas emission. To obtain the production rates, Q, we used the Haser model (Haser 1957). We performed fits of the column density profiles, N, as a function of the projected cometocentric distance, ρ. Specifically, we fitted log N − log ρ profiles in the E-W and N-S directions on 2007 December 15, and in the N-S direction on 2008 January 4. By considering the standard scale lengths given in Cochran et al. (1992), A'Hearn et al. (1995), and in Schulz et al. (1998), we produced theoretical column density profiles for each species by varying the production rate Q until the best match between observations and theoretical predictions is achieved. However, the parent and daughter scale lengths (l p , l d ) do not satisfactorily fit the observed N(ρ) regardless of the Q value (some examples are shown in Fig. 2). In general, the observed CN, C 3 , and C 2 column density profiles versus projected cometocentric distance are noticeably flatter than the results of the Haser modelling with classical parent and daughter scale lengths. This fact points to an l p , or l d or both larger than the standard values. Only slightly better results are obtained if the standard scale lengths are corrected for heliocentric distance by using a power-law index of 1.5 (subtracting our heliocentric dependence of the outflow velocity to the theoretically expected increase of the photodissociation rates).
Therefore, to retrieve more confident production rates, both parent (l p ) and daughter (l d ) scale lengths have to be fitted as well. Given the dispersion of the column densities profiles (see examples in Fig. 2), to estimate the best-fitting parameters, their significance, and errors bars, the following procedure, inspired by the work performed by Fray et al. (2005) and Langland-Shula & Smith (2011), among others, has been implemented. Column densities have been first binned in a cometocentric logarithmic scale. The error bar for each bin is defined by comparing the 5 per cent of the median with the standard deviation of the data included in the bin, and choosing the larger of the two quantities. Then, a grid of plausible parent (l p ) and daughter (l d ) scale lengths is defined, and the cometocentric Haser profile for each (l p , l d ) pair is calculated. The cometary gas production rate, Q, is then obtained by scaling all those profiles to minimize the chi-square between the theoretical and observational column densities. This allows to obtain a map of the chi-square (see examples in Fig. 3) as a function of both scale lengths and to build the corresponding chisquare distribution, which is used to estimate confidence regions. Values and errors of the best-fitting parameters are defined as the medians and standard deviations, respectively, of the parameters within the region of the 99 per cent of confidence, i.e. the sets of parameters that have a probability of less than 1 per cent of finding another set with a smaller chi-square.
The region of the 99 per cent of confidence is well constrained for most of the observational column densities (see as an example, the left-hand panel of Fig. 3). Nevertheless, some column densities do not allow to define the daughter scale length appropriately. The chi-square slightly drifts towards smaller values when the daughter scale length increases around a nearly constant parent scale length. This circumstance indicates that either the quality of the data is not sufficient, or the Haser model does not describe the underlying physics adequately, or both. An example of the chi-square confidence intervals when the daughter scale length cannot be constrained is shown in the right-hand panel of Fig. 3. In these circumstances, we have verified that a change of a factor of ∼2 in the daughter Note. Infrared images are obtained from the Spitzer Heritage Archive from the programme developed by O. Groussin as described in the text. Details on the characteristics of the grisms are given in the text.  . It can be seen that those profiles for CN, C 2 , and C 3 are considerably steeper than observations. The values of scale lengths and productions rates are shown in Table 3. scale length implies a variation of the chi-square of the order of ∼1 per cent or less. Fortunately, that circumstance, i.e. a poorly constrained daughter scale length, has a limited impact on the retrieved production rate. In regions with high confidence levels, the production rate best-fitting the data is controlled mainly by the parent scale length. Thus, in those cases, a confidence region of 95 per cent is used to determine the best-fitting parameters and their corresponding errors. The level of 95 per cent of confidence is admittedly arbitrary but, given the aforementioned slow drift of the chi-square, smaller confidence values may artificially increase the error bars with poorer representations of the data. As previously described, within the 95per cent confidence region, most of the dispersion of the Q values is due to the variation of the parent scale length within the region, regardless of the daughter scale length. Best-fitting parameters and errors defined from median values and standard deviations in the confidence region may introduce a bias, which depends on the upper limit of the daughter scale length. In order to reduce the bias introduced by the drift of the chi-square when the daughter scale length increases (which also favours smaller parent scale lengths and production rates), just in these circumstances, Q and l p are defined as the mean of the maximum and minimum values found within the 95-per cent confidence region, and the corresponding errors are calculated as half the distance between the maximum and the minimum. In these cases, l d is left unconstrained. The sets of l p , l d , and Q, and corresponding errors, that allow a good fit to the observations are listed in Table 3. Only a lower limit, obtained from the 95-per cent confidence region, is given for l d when it is unconstrained. These cases are generally characterized by large error bars both in Q and l p , although they are usually not larger than 20 per cent of the nominal value. The worst case is C 2 on 2008 January 4, in which the errors are estimated somewhat below 50 per cent. To facilitate comparison, Table 3 also shows production rates derived for standard (l p , l d ), both uncorrected and corrected for heliocentric distance. Fig. 2 shows the gas column density radial profiles observed on 2007 December 15 in the N direction together with the Haser fit using the standard (l p , l d ) and those retrieved in this work. Table 4 lists the log [Q(C 2 )/Q(CN)], log [Q(C 3 )/Q(CN)], log [Q(NH 2 )/Q(CN)], and also log [Afρ/Q(CN)] obtained from our study. As the gas-production rates are determined for several slit directions, we have averaged them before calculating the ratios.
The best-fitting procedure results in averaged Q values for CN, C 2 , and NH 2 that are significantly larger (more than a factor 2) than those derived using standard scale lengths on 2007 December 15 (r h = 1.21 au). This difference becomes smaller on 2008 January 4 (r h = 1.08 au), when the best-fitting procedure yields production rates for CN, C 2 , and NH 2 that are 74, 44, and 87 per cent, respectively, larger than those obtained with standard scale lengths. Differences for C 3 are significantly smaller, being ∼25 per cent on both dates. When the standard scale lengths are corrected for heliocentric distance, differences diminish, but they are still significant.
Asymmetries in production rates along the directions probed by the slit are seen for most molecules. On 2007 December 15, NH 2 shows the highest asymmetry, ranging from 6 to 50 per cent of the average value. CN and C 2 display more moderate asymmetries, with values between 16 and 37 per cent for the former, and between 0 and 16 per cent for the latter. C 3 shows a low asymmetry, with values ranging from 1 to 10 per cent. On 2008 January 4, asymmetries tend to reduce. NH 2 displays a 37 per cent of asymmetry, CN a 22 per cent, C 2 a 12 per cent, and C 3 a 4 per cent. Asymmetries for CN and C 2 are similar to those obtained by Jehin et al. (2009) at perihelion, about four weeks after our last observation. The main difference is found in C 3 . In our case, this molecule shows the lowest asymmetry, while Jehin et al. (2009) obtained a ∼30 per cent. Given the different slit orientations (Tail-Sun in Jehin et al. 2009), a direct comparison of the asymmetries seems difficult.
In order to compare our best-fitting scale lengths with those available in the literature, they must be converted to a common heliocentric distance. Scale lengths are generally assumed to follow a power law of the heliocentric distance, r h , i.e. r h α . Presently, different α indexes have been retrieved from various observational data sets, frequently showing discrepancies when the same molecule is considered (see e.g. Rauer et al. 2003;Langland-Shula & Smith 2011). This fact shows the complexity of the physicochemical evolution of the coma. Theoretically, if the impact of coma dynamics is minimized, an α = 2 is expected due to the dependence of the solar radiation on heliocentric distance. As in our study, an expansion velocity which depends on the heliocentric distance as r h −0.5 is considered, scale lengths may be approximately transformed to 1 au using α = 1.5. Averaging the values retrieved for the different slit orientations (when l d is constrained) and assuming α = 1.5, all results Table 3. Gas-production rates, Q, of different molecules, and the best-fitting parent and daughter scale lengths (l p , l d ) in the coma of 8P/Tuttle.  Note. For comparison, although the fit to observations is not optimal, the production rates obtained using standard (l p , l d ), both uncorrected and corrected for heliocentric distance as described in the text, are also listed. Those latter data are identified with 'st' at the corresponding line, and the solutions for the uncorrected and corrected scale lengths are separated by '//'. 2) × 10 4 km at 1 au. This length is smaller than the lower limit of 19.9 × 10 4 km for l d at 1 au of the CN scale lengths compiled by Langland-Shula & Smith (2011). Actually, our nominal value at r h = 1.21 au is already slightly under that lower limit for r h = 1 au. Given our large error bar, the discrepancy could be marginally explained by the chosen heliocentric dependence because a value of α = 1.5 may be excessive based on the results by Langland-Shula & Smith (2011). These authors found that the heliocentric dependence of the l d of CN may be characterized by an uncertain but shallow slope of 0.6 ± 0.5. Thus, in principle, our values for the l d of CN could be marginally compatible with those found in the literature if a very low slope is considered to transform them to 1 au. Nevertheless, it is also possible that CN in comet 8P, when described by Haser modelling, actually pushes down the lower limit for l d generally found in the literature. Jehin et al. (2009), on 2008 January 28, when the comet was at perihelion, obtained a value of l d = 16 × 10 4 km, already under the lower limit of the compilation by Langland-Shula & Smith (2011), as in our case. If that value is transformed to 1 au with our α = 1.5, we would obtain an l d = 15.3 × 10 4 km, compatible with our (14.6 ± 3.2) × 10 4 km, within the error bar.

T H E D U S T M O D E L
The available observations to fit with the model are the CAFOS and ALFOSC images, on 2007 December 15, and 2008 January 4, respectively, the Spitzer/MIPS 24-μm image on 2008 June 23, and the Afρ observations by CAFOS, ALFOSC, and CARA 5 and Cometas Obs 6 amateur teams, all converted to R-band photometry using an aperture of ρ = 10 4 km.
To retrieve the physical parameters of the dust from the observations, we used our Monte Carlo dust tail code. This code has been already used to characterize the dust environment in several dustemitting objects, including comets and active asteroids (for recent applications of the model, see e.g. Moreno et al. 2019;de León et al. 2020). This model simulates a dust tail from a comet, or an active asteroid, by adding up the individual contribution to the brightness of each particle ejected from the parent nucleus. The particles, after leaving the object's surface and at a distance where they become decoupled from the gas drag (some 20 nuclear radii, R n ), move mainly under the influence of the Sun's gravity and solar radiation pressure. The nucleus gravity force is neglected, an approximation valid for small-sized nuclei, such as 8P. In such conditions, the trajectory of the particles becomes Keplerian, and their orbital elements are functions of their physical properties and velocities when the dust becomes decoupled from the gas (e.g. Fulle 1989). To build up sufficiently populated synthetic images with the Monte Carlo procedure, we usually launch from 2 × 10 6 to 10 7 particles for each observation date. Since the computing time of each Monte Carlo run is proportional to the dimensions of the image to be simulated, the CAFOS and ALFOSC images have been rebinned by factors of 4 and 8, resulting in pixel sizes of 616.7 and 321.5 km px −1 , respectively.
The ratio of radiation pressure force to the gravitational force exerted on the particles is given by the parameter β, defined as In the previous expression, ρ is the particle density, and r is the radius of the particle. Q pr is the scattering efficiency for radiation pressure, here taken as 1, as it converges to that value for absorbing spherical particles of radius r 1 μm from Mie theory (see e.g. Moreno et al. 2012, their fig. 5). Finally, C pr is the radiation pressure coefficient, which is computed as where E is the total solar radiation, c is the speed of light, G is the gravitational constant, and M is the mass of the Sun. 5 Cometary Archive for Afrho, http://cara.uai.it/. 6 http://www.astrosurf.com/cometas-obs/ The brightness, m, of a particle of radius r, expressed in mag arcsec −2 , is computed from the expression πr 2 = 2.24 × 10 22 πr h 2 2 10 0.4(m −m) where r h is the comet heliocentric distance in au, is the Earth-tocomet distance in au, α is the phase angle, φ(α) is the linear phase coefficient, and m is the magnitude of the Sun in the corresponding red bandpass filter. For the simulations of the Spitzer/MIPS 24-μm image, we computed the particle flux using the equation where λ is the wavelength, is the grain emissivity, B λ is the Planck function, T(r h ) is the particle equilibrium temperature, and S is the Spitzer-to-comet distance. The Spitzer heliocentric coordinates at the time of the 2008 June observation are obtained from the JPL Ephemeris (Giorgini et al. 1996). The equilibrium temperature can be computed by the balance between the absorbed solar and emitted thermal radiation as In this equation, A B is the Bond albedo. The dust grain temperature at 1.6 au has been estimated at 258 ± 10 K by Groussin et al. (2019). Hence, the quantity ( 1−A B ) 1/4 = 1.17, so that the grain equilibrium temperature is computed as T (r h ) = 326.2/ √ r h .
Given the large number of free parameters in the model, some assumptions must be made to make the problem tractable. Thus, the particle density is assumed at ρ p = 800 kg m −3 , and a geometric albedo at red wavelengths of 0.065 is considered. These values correspond to those determined for comet 67P particles from GIADA (the Grain Impact Analyser and Dust Accumulator; Colangeli et al. 2007) and OSIRIS (the Optical, Spectroscopic and Infrared Remote Imaging System; Keller et al. 2007) measurements, on Rosetta spacecraft (Fornasier et al. 2015;Fulle et al. 2016). Further, the particle brightness is corrected for the phase function using a linear phase coefficient of φ(α) = 0.03 mag deg −1 . This value has been determined for several comets by Meech & Jewitt (1987).
Following a common practice, the particle size distribution is assumed to be described by a power-law function, i.e. n(r) ∝ r −κ . The limiting particle sizes, r min and r max , are assumed at 5 μm and 5 cm, respectively. The value of r min is found after some experimentation with the code. We found that the influence on the tail brightness of particles smaller than that size becomes negligible for the available observations. This result is in line with the findings of Groussin et al. (2019), who determined an excess of 10 μm particles compared with submicron ones from the temperature of dust grains in 8P. On the other hand, the maximum radius was set as the order-ofmagnitude size of the particles associated with the backscattered signal as detected by radar observations by Harmon et al. (2010), which were of the order of a few centimetres in size. The large upper limit for the particle radius is also justified from observations of the Ursid meteor shower, linked to comet 8P (e.g. Jenniskens et al. 2002). Moreno-Ibanez (2018) analysed several Ursid meteors finding a radius size between 1.7 and 2.4 cm, assuming that they are made of refractory material with a density of 3500 kg m −3 . The authors justified that high density by arguing that the particles were emitted from the nucleus a long time ago (several centuries), enough to lose all the volatile material initially contained in them. Assuming that the particles had the same amount of volatile and refractory material when they were released, and that the density was that considered in the present study (800 kg m −3 ), the initial radius of those particles would be between 3 and 5 cm.
The particle-size distribution was assumed constant with the comet heliocentric distance. The value of the power index κ is to be found in the parameter fitting procedure as described below.
The dust velocity after gas decoupling is assumed to follow the typical expression v(β, t) = β γ u(t), with β standing for the ratio between radiation and gravity forces (equation 1), γ is set to γ = 0.5 from gas drag hydrodynamical modelling, and u(t) is a function describing the time evolution of the velocity, which needs to be found from the fitting procedure. The velocity v(β, t) is constrained to be larger than the nucleus escape velocity. Groussin et al. (2019) (see also Lamy et al. 2008 andHarmon et al. 2010) derived a bilobated nucleus shape consisting of two contact spheres of radii 2.7 and 1.1 km. Assuming a reasonable nucleus density of 500 kg m −3 , this results in a total mass of M = 4.4 × 10 13 kg. The size of an equivalent sphere would have a radius of 2760 m, and this results in an escape velocity of v esc = 1.46 m s −1 . In the simulations concerning the visible data, we added the nucleus brightness to the synthetic dust tails generated. For this purpose, the nucleus is assumed to be a sphere of 2.76 km of radius, having a geometric albedo of 0.065 and a linear phase coefficient of 0.047 mag deg −1 , as was found for comet 67P (Fornasier et al. 2015). The effect of the nucleus brightness on the computed synthetic tail images becomes important only for large heliocentric distances. For the simulation of the Spitzer image, instead of adding the nucleus contribution, we subtracted the nucleus point spread function profile from the Spitzer image, as was estimated by Groussin et al. (2019).
In addition to the time-dependent velocity profile u(t), the dust mass-loss rate profile, dM/dt(t), must also be determined from the fitting procedure. We proceeded as follows. We approximate both time-dependent functions by Gaussian profiles. However, asymmetric branches with respect to the maxima have been found in both functions for a wide variety of comets (e.g. Cremonese & Fulle 1994;Fulle et al. 2000;Moreno et al. 2017). In order to account for these asymmetries, we introduce skew parameters. The fitting functions have the general form where A is the maximum value of f(t), t 0 is the time of the maximum, σ accounts for the width of the Gaussian, and s is the skew parameter. The nominal Gaussian function has s = 1. The Gaussian profile can be skewed to the left or to the right, so we introduce a parameter s = s l when t < t 0 , and s = s r for t > t 0 . In this way, we can make the profile narrower or broader than the nominal Gaussian, at either side of the maximum, by playing with s l and s r (narrower profiles for s < 1, and broader profiles for s > 1). Since we need the skewed Gaussians for both the dust loss rate and the velocity, this implies eight fitting parameters. The size distribution power index κ completes the set of nine parameters to be found by the multidimensional fitting method. The method used is the downhill simplex described by Nelder & Mead (1965) and is aimed at finding the best-fitting solution to the available observational data, namely, the optical images, the SpitzerIR-image, and the Afρ data from the amateur observers. In this regard, it is necessary to point out that Afρ data from the Cometas obs team between −60 and −20 d from perihelion show a strong discrepancy with the rest of Afρ data (both from the CARA team and from our own analysis, see Fig. 4). Presently, there is no satisfactory explanation for that discrepancy.  For that reason, regarding Cometas obs data, only those between −200 and −60, and between −20 and +150 d from perihelion are considered during the fitting procedure. Data from Cometas obs between −60 and −20 d from perihelion are disregarded. The bestfitting parameters are found by minimizing the squared sum of the differences between the modelled and measured tail brightness and Afρ values. The time evolution of the Afρ parameter is depicted in Fig. 4. The best-fitting time-dependent dust mass-loss rates and velocities (for particles having r = 0.1 mm) are shown in Fig. 5. The power-law index of the size distribution was found as κ = 3.2 ± 0.2. The resulting fits to the R-band images and the Spitzer IR-image are shown in Figs 6 and 7, respectively.   Fig. 8 shows CN, C 2 , C 3 , and NH 2 productions rates of 8P estimated by various authors at different pre-perihelion heliocentric distances.

Comments on the gas results
Our results for CN, C 2 , and C 3 , obtained at heliocentric distances of 1.21 and 1.08 au pre-perihelion, are compatible with those of Jehin et al. (2009), obtained at perihelion (1.027 au) with a similar procedure, i.e. applying a Haser model and allowing the scale lengths to be different from the customary values to obtain the best-fitting. CN, C 2 , and C 3 production rates derived in this work, together with the results by Jehin et al. (2009), suggest a slight increase of activity Figure 8. Determinations of CN (black), C 2 (blue), C 3 (orange), and NH 2 (green) production rates of 8P by different authors at different heliocentric distances. Lines joining symbols are intended to guide the eye to connect production rates corresponding to the same molecule.
as the comet approaches perihelion. Our NH 2 estimates also show an increase towards perihelion but, unfortunately, Jehin et al. (2009) did not obtain the NH 2 production rate. By contrast, Borisov et al. (2008), using low-resolution spectra and the scale lengths from A' Hearn et al. (1995), obtained much higher values. Although production rates estimated by different authors or using different techniques frequently show some discrepancies (see e.g. Langland-Shula & Smith 2011), the production rates obtained by Borisov et al. (2008) seem to be very high compared with production rates estimated at nearby heliocentric distances. Transients effects, such as outbursts, cannot be invoked to explain this discrepancy as no event of that type was reported during the continuous observations of the comet in its closest approach to Earth. 7 A similar situation occurs with the production rates obtained at larger pre-perihelion heliocentric distances by Langland-Shula & Smith (2011), being higher than our estimates and those by Jehin et al. (2009). The results of Schleicher (2007a,b) on CN, C 2 , and C 3 at 1.6 and 1.3 au, obtained from narrow-band imaging, are also displayed in Fig. 8. Their results on CN and C 2 , together with those of Jehin et al. (2009), and those reported in this work, show an increase of CN and C 2 towards perihelion. C 3 shows a different behaviour. Schleicher (2007b) obtains a C 3 production rate at 1.3 au similar to that obtained by Jehin et al. (2009) at perihelion, suggesting a constant production rate as the comet approaches perihelion. Our C 3 production rate is a factor between 2 and 3 below that level. In principle, that difference cannot be attributed to the fitting procedure and the determination of the scale lengths. That molecule showed little sensitivity to variations in scale lengths, as shown by the small differences between our estimates and the values obtained when the standard scale lengths were used. It is implausible to invoke a change in the composition of the region originating activity as the subsolar point moves towards the North Pole (see Fig. 9) affecting just C 3 . Being our data compatible with those of Jehin et al. (2009), currently, we do not have an explanation for the difference with the comparatively high C 3 production rates obtained farther from perihelion.
To compare 8P with other comets, it is convenient to consider the ratios 8 of the different molecules with respect to CN (see Table 4). Regarding the C 2 /CN ratio, the previously published results show some discrepancies. On the one hand, Schleicher (2007b) obtained a ratio of 0.17 from observations at 1.3 au, a value that can be considered in agreement with the ratio of 0.11 obtained by Jehin et al. (2009) at perihelion (averaging their productions in the Tail/Sun directions). Both values are in agreement with the estimated one by A'Hearn et al. (1995) at the 1994 perihelion passage, defining 8P as a 'typical' comet regarding C 2 /CN ratio. On the other hand, Langland-Shula & Smith (2011) obtained a comparatively high ratio of 0.38 (at 1.24 au) while Borisov et al. (2008) reported a C 2 /CN ratio of −0.96 (at 1.06 au), a very low value, more typical of comets at large heliocentric distances (A'Hearn et al. 1995). 8 All daughter production ratios in this section are given as the log 10 [ratio]. At 1.21 au (2007 December 15), our nominal C 2 /CN ratio (see Table 4) is the same as that of Jehin et al. (2009), and therefore compatible with those obtained by Schleicher (2007b) and A'Hearn et al. (1995). At 1.08 au (2008 January 4), the uncertainly in the C 2 /CN ratio is much larger, inherited mainly from the uncertainty in C 2 production rates. Thus, our results for that heliocentric distance cannot be interpreted with confidence.
Considering just the value obtained on 2007 December 15, i.e. 0.11 ± 0.10, it is compatible with the average value for Halley-type comets reported by Cochran et al. (2012), 0.16 ± 0.13, which is also coincident with the average value from their complete set of comets, 0.15 ± 0.20. Therefore, our results support that 8P can be considered as 'typical' regarding that ratio.
As in the case of C 2 , C 3 /CN ratio also shows some discrepancies when comparing previously published values. For this ratio, Schleicher (2007b) obtained an unusually high value of −0.40. That value is much higher than those considered by A' Hearn et al. (1995) as 'typical', around −1.0, and also higher than the average value of That change would suggest that the C 3 production remained approximately constant while that of CN increased as 8P approached perihelion. Borisov et al. (2008), from their observations at 1.06 au, also obtained a small value for the ratio, −1.04. Our results give a ratio between −1.46 ± 0.09 and −1.31 ± 0.11 for the range 1.21-1.08 au of heliocentric distances, close to the value of Jehin et al. (2009). According to the definition given by Cochran et al. (2012), our estimate for 8P would meet the condition on the C 3 /CN ratio to be considered as depleted (log [Q(C 3 /Q(CN)] < −0.86). Actually, our ratio is less than the mean value estimated from the sample of depleted comets in Cochran et al. (2012), i.e. −1.06 ± 0.39. The C 2 /CN ratio obtained implies that 8P cannot be defined as 'depleted' in the strict sense (C 2 and C 3 ), but our results, supported by those of Jehin et al. (2009), would suggest that 8P is C 3 depleted according to the criterion by Cochran et al. (2012). The situation is different if the data set by A' Hearn et al. (1995) is considered. A'Hearn et al. (1995, their table 6) provide the range of ratios characterizing comets referred to OH (rather than to CN), and the CN/OH ratio slightly depends on the taxonomic class. Schleicher (2007b) found that the CN/OH ratio at 1.3 au was −2.58. Using that value to estimate our C 3 /OH ratio, we obtain −4.04 (2007 December 15) and−3.89 (2008 January 4). Both values are within the range of the C 3 /OH ratios for 'typical' comets, [−4.26, −3.09], according to A'Hearn et al. (1995), if the C 2 /CN ratio is higher than −0.18, as in our case. Thus, classification from A' Hearn et al. (1995) allows us to define 8P as a 'typical' comet.
Regarding NH 2 /CN ratio, Langland-Shula & Smith (2011) obtained a value of 0.95 when the comet was at 1.24 au. Our results, 0.81 ± 0.10 and 0.93 ± 0.11, may be considered in comparatively good agreement with those of Langland-Shula & Smith (2011). Those high values are much larger than the mean value of −0.26 ± 0.40 that Cochran et al. (2012) estimated for 8P from data corresponding to the 1980 perihelion passage. The high values for the NH 2 /CN ratio found by Langland-Shula & Smith (2011) and in this study are, in fact, higher than any of the ratios found by Cochran et al. (2012) from their complete data set. From our results and those of Langland-Shula & Smith (2011), although statistics on NH 2 /CN ratio do not allow a meaningful conclusion, 8P would be defined as NH 2 enriched. Interestingly, Cochran et al. (2012) found that Halley-type comets, group to which 8P belongs, showed a slight tendency to enrichment in NH 2 (with a positive average ratio of 0.10 ± 0.34) compared with JFCs (−0.03 ± 0.28) and long-period comets (−0.17 ± 0.22). Therefore, even if 8P is rich in NH 2 , this circumstance seems to be common among Halley-type Comets.
Concerning parents molecules, as in Jehin et al. (2009) andBonev et al. (2008), the scale length measured for CN indicates that HCN cannot be the main parent species. The photodissociation rate of HCN at the comet heliocentric distances ranges from 8.6 × 10 −6 to 1.1 × 10 −5 s −1 , whereas those for CN span from 2.2 × 10 −6 to 2.7 × 10 −6 s −1 (Huebner, Keady & Lyon 1992). Using 0.8 km s −1 for the HCN expansion velocity (Fray et al. 2005) and 1 km s −1 for the CN ejection velocity, the equivalent parent and daughter scale lengths are in the range of 49 400−57 700 and 368 000−462 000 km, respectively. From the values found in this work (l p ranging from 27 000 to 63 000 km, and l d poorly constrained but with values as low as 140 000 km), we may conclude that HCN is not the main, or at least the unique, parent of the cyano radical. Other species whose photolysis gives rise to CN are C 2 N 2 , HC 3 N, and CH 3 CN. The cases of cyanogen and cyanoacetylene as parents of CN would have equivalent parent scale lengths within 20 000 and 32 000 km, shorter than the average value measured in this work, although compatible with our lower limit. Additionally, from our study of the (l p , l d , Q) parameter space, we can conclude that both C 2 and C 3 scale lengths (parent and daughter) are larger than customary values. Constraining the CN, C 2 , and C 3 parent species requires the development of complex chemical models beyond the scope of this paper.

Comments on the dust results.
From observations at 1.63 au pre-perihelion (at approximately 87 d before perihelion), Schleicher (2007a) found that, although the daughter molecules abundances of 8P were 'typical', the comet showed a very low dust-to-gas ratio. That low ratio was in agreement with estimates at the 1994 perihelion passage obtained by A'Hearn et al. (1995). Schleicher (2007a) concluded that 8P had a very low dust-to-gas ratio based on the quotient Afρ/Q(H 2 O), which placed the comet among those with the lowest ratio in the data set of A' Hearn et al. (1995). Combining our dust mass-loss rate with the water production estimated by Schleicher (2007a), a dust-to-water ratio of 0.98 is obtained when the comet was at 1.6 au pre-perihelion, a ratio close to the canonical value of 1. That value seems to indicate that 8P certainly had a low dust-to-water ratio, particularly if compared with the results from Sykes & Walker (1992), who estimated an average dust-to-gas ratio of 3 from the trails of an ensemble of comets, with that of C/1995 O1 (Hale-Bopp) (e.g. Weaver et al. 1999;González, Gutiérrez & Lara 2014) or with that of comet 1P/Halley (e.g. McDonnell, Lamy & Pankiewicz 1991;Fulle et al. 2000), among others. Today, it is known that the dust-to-gas ratio varies with heliocentric distance. In fact, as described in the thorough review by Choukroun et al. (2020) as well as in Fulle et al. (2019), estimating the refractory-to-ice, dust-to-gas or dust-to-water ratios is a complex task even in the case of comet 67P/Churyumov-Gerasimenko, for which a wealth of data is available thanks to both Rosetta mission and ground-based observations. Choukroun et al. (2020) estimated the mission integrated dust-to-gas ratio for 67P as 2.3. Nevertheless, this dust-to-gas ratio, which can be considered the best constrained ratio, is affected by a large uncertainty, ranging from 0.64 (Choukroun et al. 2020) up to values as high as 6 (Rotundi et al. 2015), showing a comparatively large variability if estimated from single snapshots at particular heliocentric distances. Fulle et al. (2019) convincingly argued that discrepancies between different estimates of the coma dust-to-gas ratio may converge into a comparatively high value of the nucleus refractory-to-ice ratio if phenomena, such as mass transfer between hemispheres, are considered.
In any case, 8P seems to have a low dust-to-water ratio when estimated from production rates at distances larger than at least 1.6 au pre-perihelion. Nevertheless, from our Monte Carlo modelling of the images of 8P, the dust-to-water ratio increases as the comet approaches perihelion, reaching values between 6 and 11.5, depending on the estimate of the water production rate around perihelion (5.97 x 10 28 molecule s −1 by Böhnhardt et al. 2008 or 3.0 x 10 28 molecule s −1 by Kobayashi et al. 2010, obtained just after perihelion at, 1.03 au).
It is worth mentioning that caution is necessary when interpreting the dust content of a comet from Afρ estimates. Certainly, comet 8P, from the Afρ/Q(H 2 O) or Afρ/Q(CN) ratios (see Table 4), seems to be among those with the lowest ratios in the data set by A'Hearn et al. (1995, their fig. 3). If the Afρ at perihelion is considered (approximately 240 cm from Fig. 4) together with the CN production obtained by Jehin et al. (2009), a value of −23.70 is obtained. That value is similar to those in Table 4, indicating little variation of the ratio from 1.21 au up to perihelion. Although the low dust content of the coma derived from Afρ is consistent with our modelling at large heliocentric distances, estimates based on Afρ close to perihelion show conflicting results when compared to those obtained using Monte Carlo modelling. In this regard, it is necessary to realize that Afρ is proportional to the dust loss rate times the cross-section of the dust grains divided by the mean dust velocity of the observed grains (Fulle 2000).
The increase of the dust-to-water ratio is mainly due to the relative increase of the dust mass ejection. Indeed, while water production increases by a factor between 7.5 and 15 from 1.6 au up to perihelion, dust ejected mass increases by a factor approximately an order of magnitude larger. After perihelion, dust mass ejection decreases much faster than water production. At 2.24 au post-perihelion, considering the water production estimated by Groussin et al. (2019), the dust-to-water ratio is only 0.20. Thus, the dust mass-loss rate obtained from the Monte Carlo dust tail modelling, together with the water production rates available from different observations, indicates that comet 8P shows a coma with a comparatively low dust-to-water ratio at large heliocentric distances. By contrast, close to perihelion, the ratio reaches values that can be considered high, between 6 and 12. These high values are in line with those considered by Fulle et al. (2000) to analyse 1P/Halley data from Giotto (ESA) and with the ratios found for comet C/1995 O1 (Hale-Bopp) (see e.g. González et al. 2014, their fig. 1). Unfortunately, available water production rates do not allow us to integrate the productions to give an average value for the coma dust-to-water ratio. In any case, given the behaviour of the comet near perihelion, it cannot be concluded that 8P has a low dust-to-water ratio. The ratio reaches values that can be described as 'typical'. It is necessary to notice that the ratios found in this study may be considered as lower limits if the nucleus emits particles larger than 5 cm. Very large particles would contribute efficiently to the total mass-loss, but their effect on the coma brightness would be comparatively small, making it difficult to assess their contribution with the dust model. Fig. 5 shows that peaks of particles velocity and dust production do not occur at the same time. While particle velocity peaks at perihelion, dust loss rate peaks approximately 20−25 d after. Interestingly, the same behaviour was found for comet 67P (Moreno et al. 2017). Dust production curve and asymmetry resemble those of the subsolar latitude for the spin axis determined by Groussin et al. (2019), that is (RA, Dec.) = (285 ± 12, +20 ± 5) • . For that pole solution, subsolar latitude is −40 • at −200 d from perihelion, reaches the northward equinox at −73 d, and peaks at +80 • two weeks after perihelion. Fig. 9 displays the dust mass production rate as a function of the subsolar latitude, clearly showing a correlation between both quantities. In principle, that figure may suggest different activity characteristics for both hemispheres of 8P, as was the case for 67P. As it is known, the Northern hemisphere of 67P was fed with icy particles that fell from those released from the Southern hemisphere when it was illuminated at perihelion (e.g. Thomas et al. 2015;Pajola et al. 2019;Davidsson et al. 2021). As Fulle et al. (2017) pointed out, mass transfer may apply to all comets because the asymmetric gas density cannot prevent 'airfall' 9 . This would be particularly true for comets with strong seasonal effects, as 8P with the spin axis inclination determined by Groussin et al. (2019). Indeed,  already pointed out that 8P could be affected by significant seasonal effects. The main difference between 67P and 8P is that, given the different arguments of their spin axes, in the latter comet the Northern hemisphere would be the producer of material falling back on to the Southern hemisphere.
Leaving aside the mass transfer as a process contributing to the seasonal dependence of the activity, other factors may also have contributed to the increased dust mass ejection when the subsolar latitude moves across the surface of the nucleus. For example, the activation of new source regions (as Schleicher & Bair 2008 suggested) or increased sublimation in a permanently illuminated polar region.
Another quantity that may show seasonal dependence is the slope of the dust size distribution (κ index). Comet 67P showed a transition from a slope of 3 when the Northern hemisphere was illuminated to a steeper one, between 3.6 and 4.3, when the illumination was over the Southern hemisphere (e.g. Della Corte et al. 2016;Moreno et al. 2017). Regarding 8P, given the potentially high spin axis inclination, its coma could also hide a strong seasonal variation of the dust size distribution. Unfortunately, the scarcity of available observations, particularly imaging data, prevents us from performing a more complete analysis, including time-dependent size distribution functions. The slope of the size distribution found for 8P, i.e. 3.2 ± 0.2, can be considered as 'typical' among comets, with 3.3 ± 0.2 being the average value obtained from inverse tail models (Fulle et al. 1995). Given the comparatively small slope found with our procedure, most of the dust mass in the perihelion coma of 8P is released in the form of large grains.

S U M M A RY A N D C O N C L U S I O N S
Comet 8P/Tuttle, selected as a possible backup target for the CI mission, was observed during its previous perihelion passage to obtain visible images and low-resolution spectra. Observations were performed when the comet was between 1.21 and 1.08 au from the Sun pre-perihelion, during its close approach to Earth in 2007-2008. These data help fill the gap between similar previously reported observations (e.g. Schleicher 2007a; Borisov et al. 2008;Jehin et al. 2009).
Low-resolution spectra were used to estimate CN, C 2 , C 3 , and NH 2 production rates. From our estimates and the corresponding ratios, we conclude that 8P can be defined as 'typical' regarding C 2 and C 3 according to A'Hearn et al. (1995) taxonomic classification. Under Cochran et al. (2012) criteria, 8P is found to be C 3 depleted, although it cannot be defined as a 'depleted' comet in the strict sense because it does not meet the C 2 condition. Regarding NH 2 , 8P could be considered, as found for other Halley-type Comets, NH 2 enriched according to the statistics of Cochran et al. (2012).
Our visible images, together with infrared images available at the Spitzer Heritage Archive corresponding to the programme lead by O. Groussin (Groussin et al. 2019) and the Afρ estimates by both CARA and Cometas Obs amateur teams, have been analysed by using a Monte Carlo dust tail model. From this analysis, it has been obtained that 8P/Tuttle has a dust size distribution characterized by a slope of 3.2 ± 0.2, a value that can be considered 'typical' among those retrieved using the same technique. Combining the dust loss rate obtained in this study with available water production rates, it is obtained that the dust-to-gas ratio increases from comparatively low values when the comet is at 1.6 au pre-perihelion up to values above 6 when the comet is at perihelion. After perihelion, the dust-towater ratio decreases again, reaching values as low as 0.2 at 2.24 au post-perihelion. The significant increase observed in the coma dustto-water ratio may be a consequence of the high inclination of the spin axis, favouring the existence of a region quasi permanently illuminated when the comet approaches the Sun. Seasonal mass transfer between hemispheres may also contribute to the evolution of the dust-to-water ratio over time. Our results at large heliocentric preperihelion distances confirm those obtained by Schleicher (2007a), defining the coma of 8P as one with a low dust-to-water ratio. Nevertheless, as the comet approaches the Sun, the coma dust-towater ratio increases up to values that can be considered 'typical'. A high refractory-to-ice ratio would characterize the nucleus of 8P if, as Fulle et al. (2019) pointed out, that ratio is an order of magnitude larger than that estimated from the coma.
One of the 8P most noticeable characteristics may be the high inclination of its spin axis if the determination by Groussin et al. (2019) is confirmed. That high inclination would bring the nucleus under strong seasonal effects worth studying to continue learning on the development of cometary activity. Should this comet be finally selected as backup target for the CI mission, we think that an encounter at the pre-perihelion equinox would be very valuable, and at the same time still safe, to contribute to our understanding of cometary activity. and PJG acknowledge funding from project PGC2018-099425-B-I00 (MCI/AEI/FEDER, UE). FM acknowledges support from projects P18-RT-1854 (Junta de Andalucía) and RTI-2018-095330-B-100 (MCIU). This work has made use of NASA's Astrophysics Data System Bibliographic Services and of the JPL's Horizons system.

DATA AVA I L A B I L I T Y
The data underlying this article will be shared on reasonable request to the corresponding author.