© ESO 2021

1 Introduction

Comet P/2019 LD₂ (ATLAS), hereafter LD₂, was discovered in early June 2019 by the Asteroid Terrestrial-impact Last Alert System (ATLAS; Tonry et al. 2018) as a faint asteroidal object. It was initially classified as a Jupiter trojan. In-depth inspection of images obtained in 2019 revealed a faint tail, suggesting that it presented comet-like activity¹. This cometary nature was explicitly acknowledged with the publication of MPEC 2020-K134² and CBET 4780. Prediscovery observations made by the Dark Energy Camera (DECam) on August 10, 2018, were reported on August 2, 2020, with the publication of CBET 4821 (Schambeau et al. 2020). They also reported that the object was not detected in DECam images acquired in 2017, suggesting that activity started sometime between 2017 and 2018 and that a compact coma was present at the time when the August 2018 images were taken. These observations are consistent with an upper limit of the nucleus radius of around 3 km.

As no signs of comet-like activity have been detected on a Jupiter trojan, even though it is widely accepted that they are captured objects from the outer solar system (e.g., Morbidelli et al. 2005; Nesvorný et al. 2018), the discovery of persistent activity on a putative Jupiter trojan asteroid is a very important result as it may suggest that some of them contain volatile material (likely water ice) on their surfaces. They are originally expected to consist of a combination of rock, dust, ice, and frozen gases. However, so far no activity has been detected by sublimation of water ice or other volatiles in any of them, which would confirm this hypothesis.

For this reason, we scheduled observations of LD₂ with the world’s largest optical telescope, the 10.4 m Gran Telescopio Canarias (GTC), to study the possible comet-like activity of LD₂ as soon as it became visible in May 2020. We also started a numerical exploration of its dynamical properties to determine if it could be a true member of the primordial Trojan population.

In the meantime, Kareta et al. (2020) used the orbital elements determined adding new astrometric data obtained in 2020 to show that the comet had experienced a close encounter with Jupiter on February 17, 2017, at 0.092 au, well inside the Hill radius of the planet (0.338 au). These authors concluded that LD₂ is a recently captured centaur, not a Jupiter trojan. As the value of its Tisserand parameter relative to Jupiter is T_J = 2.94, LD₂ can be classified as a Jupiter-family comet (JFC) according to Levison & Duncan (1997). Even if LD₂ is not a Jupiter trojan, it is certainly a very interesting object that could help us understand better the transition from centaur to JFC, the activation mechanisms of these bodies at large heliocentric distances, and how activity affects the surface of centaurs. LD₂ orbits the Sun just beyond Jupiter, and its activity is similar to that of 29P/Schwassmann-Wachmann-1, an object considered a prototypical “gateway” between the centaurs and JFCs by Sarid et al. (2019).

In this paper, we present the results of the observations (visible images and spectra) obtained with the 10.4 m GTC, and the dynamical properties of LD₂ derived from the analysis of an extensive sample of N-body simulations. In Sect. 2, we describe the observations and data reduction, derive the absolute magnitude and colors obtained from the images, and present the gas production rate upper limits derived from the spectra. In Sect. 3, we present the analysis of the activity based on the observed dust tail using a Monte Carlo dust scattering model. In Sect. 5, we present the results of the N-body simulations, and describe the past, present, and future dynamical evolution of LD₂. Our conclusions are laid out in Sect. 6.

2 Observations

We obtained CCD images of LD₂ on May 16, 2020, and low-resolution visible spectra on May 17, 2020, using the Optical System for Imaging and Low Resolution Integrated Spectroscopy (OSIRIS) camera-spectrograph (Cepa 2010) at the 10.4 m GTC. The observational circumstances are shown in Table 1. The OSIRIS detector is a mosaic of two Marconi 2048 × 4096 pixel CCDs. The total unvignetted field of view is 7.8′ × 7.8′, and the plate scale is 0.127′′ pix⁻¹. Standard operation mode consists of a 2 × 2 binning, with a readout speed of 200 kHz (with a gain of 0.95 e⁻/ADU and a readout noise of 4.5 e⁻). On May 16, we obtained individual images using the Sloan g’,r’,i’,z’ filters with individual exposure times of 60 s. We did one r’, g’, r’, i’, r’, z’, r’ sequence of images with the telescope tracking at the comet’s proper motion. Images were bias and flat-field corrected (using sky flats). The comet presents a conspicuous coma and tail as seen in Fig. 1. The full width at half maximum (FWHM) of the point spread function (PSF) of the comet, measured in one of the r’ 60-s images, is wider than that of the stars, 2.6′′ versus 1.5′′ (see Fig. 2) and the tail is > 1′ long (>1.9 × 10⁵ km at the comet distance).

Aperture photometry was computed using standard tasks in the Image Reduction and Analysis Facility (IRAF³), following a procedure similar to that described in Licandro et al. (2019). Using an aperture diameter equivalent to the comet’s FWHM (2.6′′), we obtained a magnituder′ = 19.34 ± 0.02 and the colors (g′− r′) = 0.78 ± 0.03, (r′ − i′) = 0.31 ± 0.03, and (i′ − z′) = 0.26 ± 0.03. The background sky was measured (and subsequently subtracted) as the median value in a region close to the comet free of coma, tail, and background stars. Flux calibration was done using GTC zero-points computed for each observing night and provided by the telescope support astronomer. The colors of the comet, which correspond to the colors of the dust coma in the aperture used, are much larger than the solar values, (g′ − r′) = 0.44 ± 0.02, (r′ − i′) = 0.11 ± 0.02, and (i′ − z′) = 0.02 ± 0.03⁴, showing that the dust coma is redder than the Sun.

From the apparent magnitude we derived the absolute magnitudes in the g’ filter using Eq. (1) from Jewitt & Luu (2019), obtaining H_g = 13.10 ± 0.03 mag. Assuming a visible geometric albedo between 0.1 and 0.04, this provides an upper limit for the nucleus radius R_N between 5.0 and 8.0 km. Considering the conspicuous activity observed, the diameter of the comet nucleus should be much smaller than these values, in agreement with the Schambeau et al. (2020) results.

We also obtained two visible spectra of P/2019 LD₂ on May 17, 2020, with the aim of looking for signatures of the typical gas species observed in comets. Each individual spectrum consisted of an exposure of 600 seconds using the R300B grism and the 1.49′′ slit width, covering a wavelength range from 3600 to 7500 Å, and with a dispersion of 4.96 Å pix⁻¹ for a 0.6′′ slit. As the aim was to look for gas species, the slit was oriented in the direction of the comet tail, not in parallactic angle. This allows a better study of the gas emission along the tail, but introduces a significant error in the slope of the spectrum due to the differentialatmospheric refraction. For this reason, we did not use this spectrum to compute the spectral slope (color) of the comet.

Spectral images were bias and flat-field corrected, using lamp flats. The 2D spectra were background subtracted, wavelength calibrated (using Xe + Ne + HgAr lamps), and flux calibrated using the spectrophotometric standard star Ross 640. The spectra were then extracted and collapsed to one dimension, using an aperture of ± 6 pixels centered at the maximum of the intensity profile of the comet. The value for the extraction aperture corresponds to the distance from the center to where the intensity decreases to 10% of the maximum. Finally, both spectra were averaged to obtain the final spectrum (see Fig. 3).

The final 2D flux-calibrated spectrum is used to analyze the gas emission of comet P/2019 LD₂. In order to visualize any emission band associated with cometary species, we used a visible spectrum of the Sun from the CALSPEC compilation (Bohlin et al. 2014) to remove the solar continuum. We scaled the Sun’s spectrum to account for the redness of the spectrum of the target, and then we subtracted it from the comet spectrum (see Fig. 3). We did not detect any evidence of CN, C₂ or C₃ emission within the 3σ level. In particular, there were no signs of the CN (0–0) emission at 3880 Å, which is usually the strongest emission observed in comets. Nevertheless, we were able to provide an upper limit to the gas production rate of CN using the same procedure as described in de León et al. (2020), obtaining Q(CN) < (1.4 ± 0.7) × 10²⁴ mol s⁻¹. We did it by using the two regions that border the CN emission band and fitting a linear continuum that was then subtracted from the spectrum. The band flux was measured and converted into column density using the g-factor in Schleicher (2010) scaled to both the heliocentric distance and velocity of the comet. Then we computed the gas production rate assuming the Haser modeling with the outflow velocity v_p scaled with r_h ( km s⁻¹), customary values for the daughter velocity v_d = 1 km s⁻¹, and scale lengths given in A’Hearn et al. (1995). Theoretical column density profiles for CN were produced for the corresponding set of parameters in the Haser modeling, varying the production rate until the best match between observationsand theoretical predictions was achieved. The non-detection of gas emission in the visible spectrum of LD₂ is not that surprising as it is very hard to detect them at large heliocentric distances, even for rather large active comets. As an example, in their large compilation of cometary gas production rates, A’Hearn et al. (1995) reported a Q(CN) = 1.7 × 10²⁴ mol s⁻¹ for comet 74P/Smirnova-Chernykh at 3.56 au. Therefore, it is expected that at the heliocentric distances of LD₂, the CN production rate, if any, should be below the value 74P at r_h = 3.56 au.

Fig. 1 Imageof P/2019 LD2 (ATLAS) obtained on May 16, 2020. The image shown is the sum of four images, 60 s exposure time each, obtained using the r′ -band filter. The field is 150′′ × 60′′; north is up, east to the left. The object presents a conspicuous comet-like coma and tail almost aligned with the extended Sun-to-target radius vector (θ).

Fig. 2 Brightness normalized radial profile of P/2019 LD2 (ATLAS), plotted as stars, compared to the profile of a field star, plotted as open circles. The profiles were obtained from one single 60 s exposure time r-band image taken on May 16, 2020. The comet profile is wider than that of the field stars because of the presence of a conspicuous coma.

Fig. 3 Average of the two individual and flux calibrated spectra of comet P/2019 LD2 (ATLAS) obtained with the 10.4 m GTC (upper panel). The spectrum of the Sun was used to obtain the continuum shown in red that was then used to remove the dust reflected spectrum and obtain the emission spectrum (bottom panel). No CN emission band was detected at the 3σ level.

3 Dust tail modeling

To gain insight into the dust physical properties of this object, we used our Monte Carlo dust tail fitting code. In addition to the GTC image obtained on May 16, 2020, and in order to place stronger constraints on the dust parameters, we have included in the analysis some of the photometric observations available, namely the object magnitudes published in the Minor Planet Center (MPEC 2020-K134). Specifically, we considered the reported ATLAS magnitudes at the time of discovery, on 2019 June 10.4, and those corresponding to Pan-STARRS1 precovery observations dated as early as May 21, 2018 and June 10, 2018. The ATLAS magnitudes correspond to the orange filter (o), while the Pan-STARRS magnitudes correspond to the wide-band purple filter (w). The measured magnitudes are given in Table 2. Our model estimates refer to r-Sloan magnitudes, but taking into account the scatter in the measured values and that both measurements refer to wide red bandpasses, we did not apply any photometric correction.

The dust tail fitting code has already been described in various papers (see, e.g., Moreno et al. 2016, 2017b, and references therein), and it has proven to be useful for the characterization of the dust environment of normal comets and main belt comets. Briefly, we assume that the ejected particles (assumed spherical) describe a trajectory that is function of the β parameter. This parameter is the ratio of the solar pressure to the solar gravitational force, and is given by (Finson & Probstein 1968), where C_pr = 1.19 × 10⁻⁴ g cm⁻², ρ_d is the particle density, d is its diameter, and Q_pr is the scattering efficiency for radiation pressure, which is Q_pr ~ 1 for large absorbing particles (e.g., Burns et al. 1979). The tail brightness is calculated as the contribution to the brightness of individual particles, whose trajectory is computed from the ejection time until the time of the observation. Those trajectories depend on β and the terminal velocities. We assume isotropic ejection from a spherical nucleus having a certain radius R_N. The particles are assumed to be distributed in size following a power-law distribution with power index κ. We assume a broad size distribution with limiting radii given by r_min = 10⁻⁴ cm and r_max = 1 cm, and κ = –2.9. This power-law exponent is a bit larger than the typical time-averaged values found in many comets (–4.1 to –3.0; see Fulle 2004), indicating a higher relative abundance of larger particles. This value was found as the one that best captures the brightness distribution along the tail (see Fig. 5, right panel). The particle density is set nominally to ρ_d = 1000 kg m⁻³, although we have also considered a higher density of ρ_d = 2500 kg m⁻³ (see Poppe 2019). The geometric albedo is also uncertain. Most centaurs have albedos in the 0.05–0.12 range (Peixinho et al. 2020), we thus set p_v = 0.07 as a typical value. A linear phase function coefficient of 0.03 mag deg⁻¹ is further assumed. For the nucleus, we also assumed p_V = 0.07, and the same density (ρ_N) as for the particles (i.e., ρ_N = ρ_d), but a steeper linear phase coefficient of 0.047 mag deg⁻¹, based on the precise estimates for the nucleus of comet 67P from Rosetta/OSIRIS measurements (Fornasier et al. 2015).

As in our previous works (see, e.g., Moreno et al. 2019), the terminal velocities are parameterized as v = v₀β^γ, where v₀ is a time-independent speed and the constant γ controls the size dependence of the speed. For activity driven by ice sublimation, hydrodynamical models predict γ ~ 0.5. However, in situ measurements of individual dust particles by Rosetta/OSIRIS and GIADA in the vicinity of comet 67P by Rotundi et al. (2015) show no dependence of particle speed on size at large heliocentric distances (i.e., γ ~ 0). Given the lack of further information, we left γ and v₀ as free parameters of the model. The remaining model parameters refer to the dust loss rate distribution. This function is parameterized by a Gaussian function with peak dust loss rate , full width at half maximum (FWHM), and time of the peak loss rate relative to the GTC observation t₀. In order to perform an adequate comparison with the observed tail, the modeled tails are convolved with a Gaussian function with a full width at half maximum equal to the prevailing seeing conditions during the night (~ 1.5′′).

One of the parameters that influence the innermost isophote levels in the computed tails is the assumed nuclear radius. This is best guessed from the Pan-STARRS 1 precovery observations owing to its faintness, which suggest little or no activity. Assuming complete inactivity, for a nuclear radius of R_n = 3 km we found r′ = 21.3 for the assumed albedo parameters, which is in line with the measured magnitudes (see Table 2). The remaining best-fit parameters are found by minimizing the squared sum of the differences between the modeled and measured tail brightness for the GTC image, and the squared sum of the differences between synthetic and measured magnitudes at the dates given in Table 2. To do this, at each iteration of the model we computed synthetic tail images at the corresponding dates, and calculated the r-Sloan synthetic magnitudes from them. For the GTC image on May 16, 2020, the region of the CCD affected by the strong brightness by the field star located near RA = 35 000 km, Dec = 5000 km (see Fig. 4, left panel) is avoided in the fitting procedure. The best fit was found by the downhill simplex method (Nelder & Mead 1965), resulting in v₀ = 0.8 m s⁻¹, γ = 0.04, = 60 kg s⁻¹, t₀ = 275 days, and FWHM = 354 days. This implies a total dust mass loss of 1.9 × 10⁹ kg since the start of the dust emission till the date of observation. With respect to the GTC image, a comparison of the isophote fields of the observation and the model is depicted in Fig. 5, left panel. The good agreement between the brightness along the tail for the observation and modeled images is shown in the right panel of Fig. 5. On the other hand, the modeled values of the magnitudes agree well with the discovery and the precovery observations (see Table 2). The time evolution of the dust loss rate is depicted in Fig. 4 where some relevant dates are indicated. As can be seen, themodel predicts that the object was already active, although with a very limited dust production of ~ 0.6 kg s⁻¹, at the time of the precovery Pan-STARRS 1 observations on May 21, 2018, and continues to be active as of the current epoch (May 16, 2020) (at 11 kg s⁻¹), spanning some 2 yr of continuous activity. The maximum level of activity of 60 kg s⁻¹ is lower than that found for other centaurs at comparable heliocentric distances. Thus, Mazzotta Epifani et al. (2006) found 100 kg s⁻¹ for P/2004 A1 (LONEOS), while Fulle (1992) and Moreno (2009) found 300–900 kg s⁻¹ for 29P/Schwassmann-Wachmann. For 174P/Echeclus, values of 20–40 kg s⁻¹ have been reported (Bauer et al. 2008), but at much larger heliocentric distance (13 au). Our reported maximum production rate is, however, much higher than typical JFCs at the same heliocentric distance. For instance, for 67P/Churyumov-Gerasimenko at 4.5 au, the dust production rate is estimated at less than 1 kg s⁻¹ (Moreno et al. 2017a). Regarding particle sizes, and as described above, we found a power-law exponent of κ = −2.9, which implies the presence of a higher amount of large particles compared with those usually found in most comets. This is in line with the large particle sizes estimated for centaurs P/2004 A1 (LONEOS) (Mazzotta Epifani et al. 2006) and 174P/Echeclus (Bauer et al. 2008).

The derived particle speeds are only weakly dependent on size (γ = 0.04), as was found in situ for comet 67P at 3.7 au (Rotundi et al. 2015). The speed (~ 0.8 m s⁻¹) corresponds to the escape velocity of a R_N = 1.1 km object with the assumed nominal density of ρ_N = 1000 kg m⁻³. This is smaller than the radius estimated above, (R_N = 3 km), but it is in line with it when taking into account the uncertainties on the model parameters, particularly density and albedo.

As stated above, we also ran the model considering a higher density for both the nucleus and the particles of 2500 kg m⁻³. In that case, a fit of similar quality to that depicted in Fig. 5 was found, but with a steeper size distribution function (κ = −3.2) and slightly smaller terminal velocities (v₀ = 0.75 m s⁻¹), the remaining model parameters being the same in both cases.

Fig. 4 Time evolution of the modeled dust loss rate from 2019 LD2 relative to the GTC observation time (2020 May 16.19). The perihelion date and the peak emission time are indicated, as well as the time of the ATLAS discovery (June 10, 2019) and the precovery observations from Pan-STARRS 1 (May 21, 2018).

Fig. 5 Fitting of the dust model with the observations. Left panel: comparison of the observed (black contours) and modeled (red contours) tail brightness isophotes. The innermost contour corresponds to 8 × 10−14 solar disk units, and the brightness decreases in factors of two outward. One solar disk intensity unit corresponds to 2.4 × 106 erg s−1 cm−2 sr−1 Å−1. The images are rotated to the conventional orientation (north up, east to the left). The x- and y-axes are labeled in kilometers projected on the sky at the object distance. The blob observed to the right of the comet optocenter is a bright star. Right panel: comparison of observed (black line) and modeled (red line) tail brightness along the direction described by the blue dotted line in the left panel.

4 Implications of the observed activity

At the time of the discovery of LD₂ on June 10.4, 2019, the model in Sect. 3 shows a dust emission of ~ 50 kg s⁻¹, which increases to its maximum value (d M∕dt = 60 kg s⁻¹) on August 15, 2019 (i.e., roughly eight months before perihelion), and then decreases again to a dust loss rate of 11 kg s⁻¹ on May 16, 2020. The fact that this object became active when approaching perihelion in the present orbit after the latest Jupiter encounter, as well as the long-lasting character of the emission pattern, strongly suggest a thermally driven process as the mechanism responsible for the activity, but the specific mechanism at play in this case is unknown.

During the activity period the r_h ranging from 4.5 to 5 au, and CO and other volatiles might be playing a role (see Womack et al. 2017, for a review on the activity of distant objects). However, if CO ice sublimation were the dominant driver, and owing to the high volatility of CO ice, we should have observed a more prominent coma around the object at the time of Pan-STARRS precovery images on May 21, 2018 (at r_h ~ 5 au), which is notthe case. This is in line with the fact that most centaurs do not show activity at large heliocentric distances, so their activity cannot be associated solely with the strong volatility of CO ice (Jewitt 2009). On the other hand, crystallization of amorphous water ice has been proposed as a plausible mechanism to trigger the outburst activity in centaurs (Jewitt 2009) and distant comets (Prialnik & Bar-Nun 1992; Capria et al. 2002), the most paradigmatic case being comet 29P/Schwassmann-Wachmann (see Gronkowski 2014, for a discussion). However, this mechanism is obviously linked to the presence of amorphous ice in the nuclei of such objects. Observations of comet 67P with the Rosetta/ROSINA instrument show that the outgassing pattern is not consistent with the presence of amorphous ice (Luspay-Kuti et al. 2016). Instead, these observations show that the nucleus contains crystalline water ice and clathrates. The activity observed in LD₂, with a peak well before perihelion, could then be associated with the sublimation of such ices, being exposed at a certain portion of the orbit according to the object’s seasons. Interestingly, as in the case of many other centaurs, the activity occurs several months before perihelion (Jewitt 2009). However, we should also consider that the current level of activity may have been triggered by a sub-catastrophic collision with a smaller body that exposed fresh volatiles from layers below the mantle. Objects moving along co-orbital or nearly co-orbital paths face an increased risk of collisions (see, e.g., de la Fuente Marcos & de la Fuente Marcos 2018). During the last 3000 yr the only other close encounters with Jupiter inside the Hill radius might have taken place nearly 170 yr ago (at perhaps ~ 0.06 au), 370 yr ago (at ~ 0.17 au), and 2642 yr ago (at ~0.27 au), for the nominal orbital evolution in Table 3 (see also Fig. 6 and the discussion in the next section regarding the predictability horizon of these calculations).

Fig. 6 Short-term evolution of comet P/2019 LD2 (ATLAS). The top three panels show the evolution of some representative parameters computed from the orbit determination in Table 3 (in black) and relevant control orbits with Cartesian vectors separated +3σ (in green), −3σ (in cyan), +9σ (in brown), and −9σ (in orange) from the nominal values in Table 4. The top panel displays the evolution of the distance from Jupiter to the comet; the red dashed line indicates the Hill radius of Jupiter at 0.338 au, and the red thick line indicates a distance of 10 Jovian radii. The second panel shows the variation over time of the value of the Tisserand parameter with respect to Jupiter TJ for the same set of reference orbits. The third panel shows the evolution of the resonant angle λr for the same sample of control orbits. The bottom panel displays the trajectory followed by P/2019 LD2 in a frame of reference centered at the Sun and rotating with Jupiter, projected on to the ecliptic plane. The diagram also includes the orbit of Jupiter, its position at (5.2, 0) au, and the Sun at (0, 0) au. This panel shows the evolution of the nominal orbital solution as in Table 3. The zero time instant in the top three panels corresponds to epoch JD 2459000.5 TDB. The blue ellipse in the bottom panel is described by Jupiter that has a small but non-zero value of the orbital eccentricity. The time span plotted in the bottom panel is the same as shown in the previous panels.

5 Past, present, and future dynamical evolution

The assessment of the dynamical evolution of LD₂ requires the analysis of an extensive sample of N-body simulations. In this work, we use the approach discussed in de la Fuente Marcos & de la Fuente Marcos (2019) and Licandro et al. (2019). The calculations were performed using the Hermite integration scheme described by Makino (1991) and implemented by Aarseth (2003). The standard version of this direct N-body code is publicly available from the web site of the Institute of Astronomy of the University of Cambridge⁵. Relative errors in the total energy for the longest integrations presented here are as low as 2 × 10⁻¹² or lower; for the shorter integrations in Fig. 6 the relative errors in the total energy are always below 5 × 10⁻¹⁷. These values are as good as those in Fig. 4 of Rein & Tamayo (2015) or better. The relative error in the total angular momentum is several orders of magnitude smaller. As pointed out by de la Fuente Marcos & de la Fuente Marcos (2012), the results from this code compare well with those from Laskar et al. (2011) among others. In order to generate the initial conditions (control orbits or clones) used in our calculations, we used the orbit determination in Table 3, which is the most recent one (as of January 27, 2021) and was released by the Jet Propulsion Laboratory’s Solar System Dynamics Group Small-Body Database (JPL’s SSDG SBDB, Giorgini 2015)⁶.

5.1 Current dynamical status

Comet LD₂ was initiallyclassified as a Jupiter trojan even though the ephemerides showed that it had experienced a close encounter with Jupiter on February 17, 2017, at 0.092 au, well inside the Hill radius of the planet (0.338 au). Kareta et al. (2020) used an early orbit determination, less precise than the one considered here (123 observations spanning 704 days versus 555 observations spanning 960 days for the orbit determination; see Table 3) to conclude that LD₂ is an active centaur instead of a Jupiter trojan. Hsieh et al. (2021) arrived at similar conclusions: it was a centaur prior to July 2018, then aJovian co-orbital, before returning to centaur after February 2028. They predicted that it will become a JFC after February 2063. The orbit determination used by Hsieh et al. (2021) included 168 observations for a data-arc span of 741 days; therefore, it is also less precise than the orbit investigated here. Steckloff et al. (2020) used the same orbit determination considered by Kareta et al. (2020) to conclude that the object will become a member of the JFC dynamical group after 2063.

Figure 6 shows the evolution of representative control orbits with Cartesian vectors separated ±3σ and ±9σ from the nominal values in Table 4. Figure 6, third panel, confirms that LD₂ is not a Jupiter trojan; the evolution of the resonant angle λ_r (relative mean longitude with respect to Jupiter) does not exhibit an oscillation about ± 60deg that is the condition to pursue a tadpole orbit in a frame of reference rotating with Jupiter. Figure 6, third and fourth panels, show that its current dynamical status is similar to that of a quasi-satellite (for about 11 yr) in which the minor body seems to orbit around the planet, although it is not gravitationally bound to it. Figure 6, first, second and third panels, show that the short-term orbital evolution into the past of LD₂ is robust asthe object remains within a safe distance from Jupiter (during the last ~170 yr), well beyond the Hill radius of the planet. In sharp contrast, the future orbital evolution beyond 2063 is highly uncertain due to a very close encounter with Jupiter (see below). In summary, we conclude that LD₂ is now an ephemeral co-orbital of Jupiter, following what looks like a short arc of a quasi-satellite cycle that started in 2017 and will end in 2028. A number of less unstable Jovian quasi-satellites have already been documented (see, e.g., Wajer & Królikowska 2012; de la Fuente Marcos & de la Fuente Marcos 2016); therefore, the current dynamical status of LD₂ is not at all surprising.

Figure 6 indicates that the predictability horizon as defined originally by Lighthill (1986) for P/2019 LD₂ spans about 43 yr into the future and about 170 yr into the past. Due to its close encounters with Jupiter, it is not possible to predict the evolution of this object for an arbitrarily long amount of time, its trajectory is essentially unstable. This also affects the value of the Tisserand parameter (T_J; see Fig. 6, second panel). Although the present-day value of T_J is < 3 and consistent with JFC membership according to Levison & Duncan (1997), after 2063, the probability of becoming a centaur is slightly higher than that of remaining as a JFC, but the nominal evolution (in black) is consistent with JFC status.

5.2 Future orbital evolution

As pointed out above, the orbital evolution of P/2019 LD₂ a few hundred years into the past and up to a few decades into the future can be precisely predicted. However (as shown in Fig. 6), on January 23, 2063, this object will experience a very close encounter with Jupiter at about 0.016 au. Considering the current uncertainty of its orbit determination, the computed range of minimum approach distance excludes that the comet may get as close to Jupiter as ten Jovian radii. This very close and slow flyby strongly affects our ability to make reliable dynamical predictions beyond early 2063. In fact, P/2019 LD₂ may not survive the close approach in one piece due to the strong tidal forces that will presumably occur during the event (tidal breakup may require an approach to 0.001 au in the case of Jupiter); on the other hand, such a close flyby may even lead to a collision with one of the Jovian moons (see, e.g., Melosh & Schenk 1993). Assuming that P/2019 LD₂ survives its flyby in 2063, integrations indicate that its probability of escaping the solar system during the next 0.5 Myr is 0.53 ± 0.03. In general, comets following trajectories similar to that of P/2019 LD₂ are expected to collide with either the Sun or one of the planets, or to abandon the solar system, either reaching the Oort Cloud or venturing into interstellar space within a timescale of a few million years (see, e.g., Levison & Duncan 1997; Di Sisto & Rossignoli 2020). Our calculations appear to confirm a similar outcome for P/2019 LD₂.

For these longer calculations into the future and those into the past discussed in the next section, we used the Monte Carlo using the Covariance Matrix (MCCM) methodology described by de la Fuente Marcos & de la Fuente Marcos (2015) in which a Monte Carlo process generates control or clone orbits (500) based on the nominal orbit, but adds random noise on each orbital element by making use of the covariance matrix that was retrieved from JPL’s SSDG, Horizons On-Line Ephemeris System.

5.3 Past orbital evolution: possible origin

Comet P/2019 LD₂ has a current value of the Tisserand parameter of 2.94; therefore, and following Levison & Duncan (1997), it is a JFC. These comets move in very unstable orbits as they experience slow close encounters with Jupiter; although originally thought to have an origin in the trans-Neptunian belt (see, e.g., Fernandez 1980; Levison & Duncan 1997), it is now widely assumed that this population has its source in the scattered belt (see, e.g., Di Sisto et al. 2009; Brasser & Wang 2015). We performed integrations backward in time using MCCM to generate control orbits to find that the probability of this comet having been captured from interstellar space during the last 0.5 Myr is 0.49 ± 0.02 (average and standard deviation). This probability increased to 0.67 ± 0.06 for integrations backward in time for 1 Myr, 0.83 ± 0.06 for 3 Myr, and to 0.91 ± 0.09 for 5 Myr integrations. The most simple interpretation of these results is that P/2019 LD₂ almost certainly arrived from interstellar space during the last few million years. It is therefore a dynamically young object, not an object that has remained in its present trajectory since the formation of the solar system. For this reason, the situation described here is very different from that discussed by Morbidelli et al. (2020), who showed that any interstellar planetesimals trapped during the formation of the solar system are highly unlikely to remain with us.

Figure 7, left panel, shows the distribution of inbound velocities for virtual objects (control orbits, see above) in hyperbolic paths with respect to the barycenter of the solar system 0.5 Myr (black), 1 Myr (violet), 3 Myr (blue), and 5 Myr (green) into the past. The median, and the 16th and 84th percentiles of the velocity distributions are km s⁻¹ (0.5 Myr), km s⁻¹ (1 Myr), km s⁻¹ (3 Myr), and km s⁻¹ (5 Myr). InFig. 2 of Hands & Dehnen (2020) the distribution of hyperbolic excess velocities for captured interstellar objects from simulations is shown; it exhibits a maximum at about 0.6 km s⁻¹, which corresponds to an inbound velocity at large distance from the solar system of ~ −0.6 km s⁻¹. Most virtual interstellar objects associated with P/2019 LD₂ have inbound velocities at over 1 pc from the Sun close to the most probable value in Hands & Dehnen (2020) (see Fig. 7, left panel, green histogram). On the other hand, Fig. 7, right panel, shows the distribution of inbound velocities for the case of hyperbolic comet C/2018 F4 (PANSTARRS) as discussed in de la Fuente Marcos et al. (2019); for this object the dynamics and the observational data both strongly suggest that it is a former member of the Oort Cloud (Licandro et al. 2019). Figure 8 uses data from the same simulations plotted in Fig. 7 and further hints at an extrasolar origin for P/2019 LD₂; most clones are found well beyond the radius of the Hill sphere of the solar system 5 Myr ago. Although an origin outside the solar system for P/2019 LD₂ seems plausible and reasonably well supported by the available evidence, we also have to admit that the actual sequence of events that led to what is observed today could have been more complex. A former member of the scattered disk may have experienced a very close encounter with one of the giant planets after becoming part of the centaur dynamical class in the relatively recent past, less than 0.5 Myr ago; such an encounter may have produced a fragmentation event induced by the planetary tidal force that was eventually able to form the observed present-day P/2019 LD₂. The feasibility of such events was dramatically confirmed by comet Shoemaker-Levy 9 in 1992 (see, e.g., Nakano et al. 1993). However, the available orbit determination of P/2019 LD₂ is not robust enough to either confirm or reject this more complex scenario that seems a priori more likely, taking into account its low orbital inclination.

Fig. 7 Distributions of inbound velocities for virtual objects associated with comet P/2019 LD2 (ATLAS) (left) and with comet C/2018 F4 (PANSTARRS) (right). For P/2019 LD2 the results of integrations 0.5 Myr (black), 1 Myr (violet), 3 Myr (blue), and 5 Myr (green) into the past are shown; for C/2018 F4, the results of the 1 Myr (violet) and 3 Myr (blue) integrations discussed in de la Fuente Marcos et al. (2019) are displayed. The bins were computed using the Freedman and Diaconis rule implemented in NumPy (van der Walt et al. 2011). Counts were used to form a probability density such that the area under the histogram sums to one. The vertical red line indicates the value − 0.6 km s−1 discussed by Hands & Dehnen (2020).

Fig. 8 Values of the barycentric distance as a function of the velocity parameter 0.5 Myr (black), 1 Myr (violet), 3 Myr (blue), and 5 Myr (green) into the past (left panel) and future (right panel) for 500 control orbits of P/2019 LD2 (ATLAS). The velocity parameter is the difference between the barycentric and escape velocities at the computed barycentric distance in units of the escape velocity. Positive values of the velocity parameter identify control orbits that could be the result of capture (left panel) or lead to escape (right panel). The thick black line corresponds to the aphelion distance, a (1 + e), limiting casee = 1, that defines the domain of dynamically old comets with a−1 > 2.5 × 10−5 au−1 (see Królikowska & Dybczyński 2017); the thick red line indicates the radius of the Hill sphere of the solar system (e.g., Chebotarev 1965).

6 Conclusions

In this paper, we presented observations of Jupiter’s transient co-orbital comet P/2019 LD₂ (ATLAS) obtained on May 16, 2020 and May 17, using the OSIRIS camera-spectrograph at the 10.4 m GTC. We used the obtained images and spectra to characterize its overall level of cometary activity. We also carried out direct N-body simulations to investigate its orbital evolution. Our conclusions can be summarized as follows:

• (i)

  LD₂ shows a conspicuous coma and tail with a length of about 1′;

• (ii)

  There is no evidence of CN, C₂, or C₃ emission within the 3 σ level in the acquired spectrum. In particular, there is no sign of the CN (0–0) emission at 3880 Å that it is usually the strongest emission observed in comets. We obtain an upper limit to the CN gas production rate Q(CN) < (1.4 ± 0.7) × 10²⁴ mol s⁻¹. The non-detection of CN at 4.5 au for a small comet like P/2019 LD₂ has allowed us to place a rather conservative upper limit. A’Hearn et al. (1995) reported a production rate Q(CN) = 1.7 × 10²⁴ mol s⁻¹ for Jupiter family comet 74P/Smirnova-Chernykh at 3.56 au. Therefore, it is expected that at the heliocentric distances of LD₂, the CN production rate, if any, should be below the value 74P at r_h = 3.56 au. The acquired data does not allow us for a more stringent upper limit;

• (iii)

  The comet brightness in a 2.6′′ aperture diameteris r′ = 19.34 ± 0.02. The coma is redder than the Sun with colors (g′− r′) = 0.78 ± 0.03, (r′− i′) = 0.31 ± 0.03, and (i′− z′) = 0.26 ± 0.03;

• (iv)

  According to our model, the dust emission of LD₂ can be described by a Gaussian with a FWHM = 354 days, a maximum = 60 kg s⁻¹ reached on August 15, 2019 (t₀ = 275 days before the observations), which then decreases again, with a current (May 16, 2020) dust loss rate of 11 kg s⁻¹. This implies a total dust mass loss of 1.9 × 10⁹ kg since the start of the dust emission, and very little activity (0.6 kg s⁻¹) at the time of precovery Pan-STARRS 1 observations in May 2018;

• (v)

  The origin of the observed activity is most likely linked to a thermally driven process, associated with sublimation of crystalline water ice and clathrates, either by a seasonal effect or triggered by a collision with a smaller body;

• (vi)

  From the image photometry, we obtained a lower limit for the absolute magnitude H_g = 13.10 ± 0.03 mag thus an upper limit for the nucleus radius R_N between 5.0 and 8.0 km. With the Monte Carlo dust tail fitting code and considering the precovery Pan-STARRS 1 magnitude data of the nearly bare nucleus, a nuclear radius of ~3 km is derived, and this value is found to be compatible with the GTC dust tail brightness in the near-nucleus region. The derived particle speed of ~0.8 m s⁻¹ corresponds to the escape velocity from a R_N = 1.1 km object with the assumed nominal density of ρ_N = 1000 kg m⁻³. All results show that LD₂ is a kilometer-sized object, in the typical size-range of the JFCs;

• (vii)

  LD₂ is now an ephemeral co-orbital comet of Jupiter, following what looks like a short arc of a quasi-satellite cycle that started in 2017 and will end in 2028;

• (viii)

  LD₂ will experience a very close encounter with Jupiter at perhaps 0.016 au on January 23, 2063. If it survives the close approach, its probability of escaping the solar system during the next 0.5 Myr is 0.53 ± 0.03;

• (ix)

  The origin of LD₂ is still an open question. The probability of this comet having been captured from interstellar space during the last 0.5 Myr is 0.49 ± 0.02 (average and standard deviation), 0.67 ± 0.06 during the last 1 Myr, 0.83 ± 0.06 for 3 Myr, and of 0.91 ± 0.09 for 5 Myr, suggesting that LD₂ may be a temporarily captured interstellar comet. However, it cannot be discarded that a very close encounter with one of the giant planets of a former member of the scattered disk may have triggered a fragmentation event that was eventually able to produce the observed present-day LD₂.

Although the physical characterization of this object based on the data presented here can be regarded as robust, we note that reconstructing its past behavior as well as predicting its future dynamical evolution remains very challenging within the context of its current orbit determination. If observations acquired prior to 2018 are made public (they have already been found, as discussed by Kareta et al. 2021), the uncertainty associated with its past and future will be reduced considerably.

After the acceptance of this paper, an initial characterization of LD₂ was published by Bolin et al. (2021). They report an absolute magnitude H_V = 15.53 ± 0.05, which corresponds to a nucleus radius R_N = 2 km assuming a value for the albedo p_V = 0.07 (the one used in our dust model). Although their value is slightly smaller than our R_N = 3 km determination, it lends further support to our conclusion that LD₂ is a kilometer-sized object. The colors they report ((g′− r′) = 0.60 ± 0.03, (r′ − i′) = 0.18 ± 0.05, (i′ − z′) = 0.01 ± 0.07) are also slightly different. As this is an extended object the aperture used is important for comparison purposes. Using the same aperture (equivalent to 10 000 km at the comet distance) we obtain (g′ − r′) = 0.73 ± 0.03, (r′ − i′) = 0.37 ± 0.03, (i′ − z′) = 0.21 ± 0.07, still redder than Bolin’s values. Based on tail dimensions, they assume the dust tail as being populated by particles with a mean radius of only 400 μ. The inferred velocity of ~ 1 m s⁻¹ agrees with our findings. However, the total dust loss rate they infer is smaller than ours by about one order of magnitude, which we attribute to the different dust models used, mainly the consideration of a size distribution in our more realistic Monte Carlo dust tail model. They also did not detect C₂ in the spectrum of LD₂, providing an upper limit for its production rate that is 5.4 times higher than our computed upper limit for the CN production rate. The Q(C₂/CN) production rate is lower than 2 for the large majority of observed comets (A’Hearn et al. 1995) so our spectrum is likely more sensitive to gas production.

Acknowledgements

J.dL. acknowledges support from MINECO under the 2015 Severo Ochoa Program SEV-2015-0548. F.M. and L.L. acknowledge financial support from the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa" award to the Instituto de Astrofísica de Andalucía (SEV-2017-0709). F.M. also acknowledges financial support from the Spanish Plan Nacional de Astronomía y Astrofísica LEONIDAS project RTI2018-095330-B-100, and project P18-RT-1854 from Junta de Andalucía. M.DeP. acknowledges funding from the Preeminent Postdoctoral Program of the University of Central Florida and from the SRI/FSI project “Finding the recipe to cook a primitive small body in the Solar System”. N.P.-A. acknowledges support from SRI/FSI funds through the project “Diggin-up Ice Rocks in the Solar System”. This research was partially supported by MINECO under grant ESP2017-87813-R. RdlFM and CdlFM thank S.J. Aarseth for providing one of the codes used in this research and for comments on the evolution of hyperbolic comets, and A.I. Gómez de Castro for providing access to computing facilities. Part of the calculations and the data analysis were completed on the Brigit HPC server of the “Universidad Complutense de Madrid”, and we thank S. Cano Alsúa for his help during this stage. In preparation of this paper, we made use of the NASA Astrophysics Data System, the ASTRO-PH e-print server, and the MPC data server.

References

All Tables

All Figures

	Fig. 1 Imageof P/2019 LD2 (ATLAS) obtained on May 16, 2020. The image shown is the sum of four images, 60 s exposure time each, obtained using the r′ -band filter. The field is 150′′ × 60′′; north is up, east to the left. The object presents a conspicuous comet-like coma and tail almost aligned with the extended Sun-to-target radius vector (θ).
In the text

	Fig. 2 Brightness normalized radial profile of P/2019 LD2 (ATLAS), plotted as stars, compared to the profile of a field star, plotted as open circles. The profiles were obtained from one single 60 s exposure time r-band image taken on May 16, 2020. The comet profile is wider than that of the field stars because of the presence of a conspicuous coma.
In the text

	Fig. 3 Average of the two individual and flux calibrated spectra of comet P/2019 LD2 (ATLAS) obtained with the 10.4 m GTC (upper panel). The spectrum of the Sun was used to obtain the continuum shown in red that was then used to remove the dust reflected spectrum and obtain the emission spectrum (bottom panel). No CN emission band was detected at the 3σ level.
In the text

	Fig. 4 Time evolution of the modeled dust loss rate from 2019 LD2 relative to the GTC observation time (2020 May 16.19). The perihelion date and the peak emission time are indicated, as well as the time of the ATLAS discovery (June 10, 2019) and the precovery observations from Pan-STARRS 1 (May 21, 2018).
In the text

Fig. 5 Fitting of the dust model with the observations. Left panel: comparison of the observed (black contours) and modeled (red contours) tail brightness isophotes. The innermost contour corresponds to 8 × 10−14 solar disk units, and the brightness decreases in factors of two outward. One solar disk intensity unit corresponds to 2.4 × 106 erg s−1 cm−2 sr−1 Å−1. The images are rotated to the conventional orientation (north up, east to the left). The x- and y-axes are labeled in kilometers projected on the sky at the object distance. The blob observed to the right of the comet optocenter is a bright star. Right panel: comparison of observed (black line) and modeled (red line) tail brightness along the direction described by the blue dotted line in the left panel.

In the text

Fig. 6 Short-term evolution of comet P/2019 LD2 (ATLAS). The top three panels show the evolution of some representative parameters computed from the orbit determination in Table 3 (in black) and relevant control orbits with Cartesian vectors separated +3σ (in green), −3σ (in cyan), +9σ (in brown), and −9σ (in orange) from the nominal values in Table 4. The top panel displays the evolution of the distance from Jupiter to the comet; the red dashed line indicates the Hill radius of Jupiter at 0.338 au, and the red thick line indicates a distance of 10 Jovian radii. The second panel shows the variation over time of the value of the Tisserand parameter with respect to Jupiter TJ for the same set of reference orbits. The third panel shows the evolution of the resonant angle λr for the same sample of control orbits. The bottom panel displays the trajectory followed by P/2019 LD2 in a frame of reference centered at the Sun and rotating with Jupiter, projected on to the ecliptic plane. The diagram also includes the orbit of Jupiter, its position at (5.2, 0) au, and the Sun at (0, 0) au. This panel shows the evolution of the nominal orbital solution as in Table 3. The zero time instant in the top three panels corresponds to epoch JD 2459000.5 TDB. The blue ellipse in the bottom panel is described by Jupiter that has a small but non-zero value of the orbital eccentricity. The time span plotted in the bottom panel is the same as shown in the previous panels.

In the text

Fig. 7 Distributions of inbound velocities for virtual objects associated with comet P/2019 LD2 (ATLAS) (left) and with comet C/2018 F4 (PANSTARRS) (right). For P/2019 LD2 the results of integrations 0.5 Myr (black), 1 Myr (violet), 3 Myr (blue), and 5 Myr (green) into the past are shown; for C/2018 F4, the results of the 1 Myr (violet) and 3 Myr (blue) integrations discussed in de la Fuente Marcos et al. (2019) are displayed. The bins were computed using the Freedman and Diaconis rule implemented in NumPy (van der Walt et al. 2011). Counts were used to form a probability density such that the area under the histogram sums to one. The vertical red line indicates the value − 0.6 km s−1 discussed by Hands & Dehnen (2020).

In the text

Fig. 8 Values of the barycentric distance as a function of the velocity parameter 0.5 Myr (black), 1 Myr (violet), 3 Myr (blue), and 5 Myr (green) into the past (left panel) and future (right panel) for 500 control orbits of P/2019 LD2 (ATLAS). The velocity parameter is the difference between the barycentric and escape velocities at the computed barycentric distance in units of the escape velocity. Positive values of the velocity parameter identify control orbits that could be the result of capture (left panel) or lead to escape (right panel). The thick black line corresponds to the aphelion distance, a (1 + e), limiting casee = 1, that defines the domain of dynamically old comets with a−1 > 2.5 × 10−5 au−1 (see Królikowska & Dybczyński 2017); the thick red line indicates the radius of the Hill sphere of the solar system (e.g., Chebotarev 1965).

In the text