Pre-discovery Activity of New Interstellar Comet 2I/Borisov beyond 5 au

The ZTF is a wide-field optical survey utilizing the 1.2 m Palomar Oschin Schmidt and a dedicated camera with a 55 deg² field of view (Bellm et al. 2019; Graham et al. 2019; Masci et al. 2019). Apart from two baseline surveys and a number of mini-surveys, ZTF executes a mini-survey ("Twilight Survey") that observes regions with solar elongation down to 35° (Ye et al. 2020), which is particularly suitable for the search of pre-discovery detections of 2I. The first phase of the ZTF Twilight Survey was executed from 2018 November to 2019 June, fortuitously covering the period before 2I reached solar conjunction in 2019 July.

Using the ZChecker comet monitoring package (Kelley et al. 2019) and JPL orbit solution #12 (the most recent JPL solution at the time that the search was conducted), we identified a total of 202 images from 2018 October 1 to the discovery date of 2I (2019 August 30) that potentially contained the comet. All images use 30 s exposure times, and were taken in a wide range of observing circumstances, with 5σ limiting magnitudes between 16 and 21 and seeing between 1'' and 5''. The comet reached solar conjunction in 2019 July, and the last three sets of images acquired before that were taken in the course of the Twilight Survey on 2019 April 29, May 2, and 5 respectively, when 2I was at r_H = 5.20–5.09 au. Prior to these Twilight Survey data, the last image was taken on 2019 April 16, when the comet was at r_H = 5. 45 au. The April 29 and May 5 observations were not optimal because only a few images were acquired, there was a high sky background, and/or the comet was near image edges. The May 2 observations were particularly fortuitous as the comet was located in an overlapping strip between two fields, and thus had 2× more images than we typically have (eight images total, instead of the typical four images per field). Therefore, we focused on the 2019 May 2 data for our initial search of a pre-discovery detection.

The astrometric observations available up to this point did not provide meaningful constraints on the cometary non-gravitational perturbations (e.g., Yeomans et al. 2004), which lead to different orbit solutions. To capture the ephemeris variations using different model assumptions, we fit the available astrometry as of 2019 October 3 using a gravity-only model, and two non-gravitational models that either assumed sublimation of H₂O (Marsden et al. 1973) or CO (with sublimation rate ∝r_H⁻²)¹³ to be the primary driver of comet activity. We then combined all eight images from May 2 following the apparent motion of 2I. The three solutions are tabulated in Table 1 (together with the final orbit, after we have successfully identified the comet in the ZTF pre-discovery data, as we discuss below), and the combined image as well as the uncertainty ellipses of the three solutions are shown as Figure 1. Because of the short arc and the low elongation of many of the astrometric positions, we note that the estimated non-gravitational parameters could be unreliable. In particular, the parameters for the H₂O-driven model appear to be too large to be credible, and are possibly caused by astrometric biases at such a small solar elongation.

Zoom In Zoom Out Reset image size

Table 1.  Orbital Elements of Various Solutions in Figure 1, as Well as the Best Solution (JPL 37), which Uses Precovery Data and a CO-driven Non-gravitational Forces

Parameter	Gravity-only	Nongrav-H2O	Nongrav-CO	Final-CO
	(without precovery)	(without precovery)	(without precovery)	(with precovery)
Epoch (TDB)	2019 Sep 16.0	2019 Sep 16.0	2019 Sep 16.0	2019 Dec 20.0
Perihelion time (TDB)	2019 Dec 8.65 ± 0.07 d	2019 Dec 13 ± 1 day	2019 Dec 9.1 ± 0.4 day	2019 Dec 8.551 ± 0.001 day
Perihelion distance q (au)	2.003 ± 0.003	1.85 ± 0.05	1.98 ± 0.02	2.00664 ± 0.00004
Eccentricity e	3.34 ± 0.01	2.7 ± 0.2	3.27 ± 0.07	3.3576 ± 0.0003
Inclination i	4408 ± 003	456 ± 05	443 ± 02	440515 ± 0.0004
Long. Ascending Node Ω	30812 ± 003	3067 ± 04	3080 ± 02	3081488 ± 00004
Argument of Perihelion ω	20921 ± 006	213° ± 1°	2097 ± 05	2091227 ± 00009
Radial accel. A1 (au day−2)	⋯	(3.5 ± 0.7) × 10−4	(1.4 ± 1.5) × 10−5	(−4.4 ± 3.2) × 10−8
Transverse accel. A2 ($\mathrm{au}\,{\mathrm{day}}^{-2}$)	⋯	(−2.7 ± 0.3) × 10−4	(−6.2 ± 5.6) × 10−6	(1.3 ± 0.3) × 10−7

Download table as:  ASCIITypeset image

We identified a possible detection at a signal-to-noise ratio of ∼10 on the combined May 2 image (Figure 2). The detection is within the 3σ uncertainty ellipses of gravity-only and CO non-gravitational solutions, and about 5'' southeast of the predicted Minor Planet Center (MPC) position. Visually, the object is about ∼5'' in diameter and is largely circular in shape, with no apparent sign of a tail. The object is barely visible on individual frames (Figure 2), with a motion consistent with that of 2I.

Zoom In Zoom Out Reset image size

We then examined other images in the time period in question, using the same shift-and-stack technique as outlined above. The detection on 2019 May 2, if real, would have greatly reduced the ephemeris uncertainty from a few arcminutes to a few arcseconds back to 2019 January, and would enable more pre-discovery data to be found. Since the comet is extremely faint, in order to eliminate contamination due to variable sky conditions and passing background stars, we only use frames that have (1) a 5σ limit of r ≳ 19.5, (2) an average full-width-half-maximum (FWHM) of <3 pixels, and (3) no background stars that are within 10'' from the predicted position of the comet.

All the pre-discovery detections and non-detections are summarized in Table 2. We were able to trace the object back to 2018 December 13. Apart from the 2019 May 2 data, the object is not visible in individual frames, and a clear detection usually requires stacking images from multiple nights. By including these astrometric positions¹⁴ with the post-discovery astrometric measurements of 2I and considering the non-gravitational effect, we were able to achieve a satisfactory orbital fit with residuals of order 1'', which is slightly above the average compared to typical ZTF astrometry (better than 05). This is due to a weak systematic bias in the data, which is possibly caused by the fact that most astrometric data were taken at low solar elongation and therefore at high airmass, introducing some differential color refraction bias (see also the discussion in Section 4.4). Nevertheless, we identify the object observed from 2018 December 13 to 2019 May 5 as 2I.

Table 2.  Summary of All the Pre-discovery Observations

Images date	Median date (UT)	Survey	rH (au)	Δ(au)	Images used	FWHM	Res.	rPS1 mag
2018 Oct 31–2018 Nov 21	2018 Nov 8.82	ZTF	8.55	7.90	g: 12; r: 16	14–30	⋯	m3σ > 22.69
2018 Dec 13–2018 Dec 22	2018 Dec 19.15	ZTF	7.75	6.99	g: 6; r: 6	19–27	07	21.19 ± 0.15
2019 Jan 17	2019 Jan 17.30	PS1	7.18	6.58	i: 4	09–11	⋯	⋯
2019 Feb 24–2019 Feb 26	2019 Feb 25.18	ZTF	6.42	6.26	g: 3; r: 4	16–36	06	20.97 ± 0.18
2019 Mar 1	2019 Mar. 1.10	CSS	6.35	6.24	Clear: 4	∼3''	11	≳21
2019 Mar 16–2019 Mar 18	2019 Mar 17.18	ZTF	6.02	6.14	g: 4; r: 5	21–29	15	20.46 ± 0.16
2019 Apr 9–2019 Apr 13	2019 Apr 12.16	ZTF	5.53	5.97	g: 3; r: 3	16–30	14	20.42 ± 0.09
2019 May 2	2019 May 2.16	ZTF	5.15	5.79	r: 8	20–21	11	20.12 ± 0.11
2019 May 5	2019 May 5.15	ZTF	5.09	5.76	r: 4	18–25	07	20.16 ± 0.21

Note. The heliocentric and the geocentric distances of the comet (r_H and Δ) are given at the median time of the image epochs. FWHM is the average full-width-half-maximum of the image.

Download table as:  ASCIITypeset image

A point worth addressing is the non-detection in November 2018. The 3σ uncertainty of the position is ∼4'' (4 pixels) along the major axis, and the motion rate difference between different orbit solutions is less than one-fourth of a pixel over the entire time span, but no comet is detected in the stack (Figure 3). A subsequent light-curve analysis (see Section 4.1) shows that 2I should be ∼1 mag above than the 3σ limit of the stack.

Zoom In Zoom Out Reset image size

The Pan-STARRS survey (Chambers et al. 2016) is a wide-field asteroid survey comprised of two identical 1.8 m telescopes, the Pan-STARRS1 (PS1) and Pan-STARRS2 (PS2). The survey has an image archive extending back to 2010. Using the ZTF precovery data along with the available astrometry from the Minor Planet Center up to 2019 October 1, we generate an ephemeris covering the period from 2018 January 1 (when 2I was well below any practical ground-based telescope sensitivity) until the discovery date of 2019 August 30. Four 45 s i-band images taken by the PS1 system on 2019 January 17 were identified, as summarized in Table 2. The normal processing for Pan-STARRS applies masking to remove areas not optimal for photometry because of non-uniform charge transfer efficiency, and for this study, the exposures were reprocessed without this mask. They were then visually inspected carefully over a large region centered at the expected ephemeris. The predicted ephemeris and uncertainty, however, place 2I firmly into a 70''-wide chip gap (Figure 4).

Zoom In Zoom Out Reset image size

While no detections were found in the Pan-STARRS pre-discovery images, the fact that the expected location of 2I is contained within a chip gap favors the ZTF precovery positions being correct, as even a weakly active comet should have been visible otherwise given the G ∼ 23 limiting magnitude and near arcsecond seeing. The chip gap only occupies a small fraction (∼10% based on the astrometry up to 2019 October 1) of the uncertainty ellipse, therefore there was a much larger chance of being proven wrong.

However, we note that PS1 was not operational between 2018 August 23 and December 12 due to a dome shutter failure, and between 2019 February 10 and March 27 due to loss of power after a winter storm on Haleakala where the telescope is installed.

The CSS (e.g., Christensen et al. 2018) operates three telescopes dedicated to the discovery and follow-up of NEOs: the 0.7 m Catalina Schmidt at Mt. Bigelow, and 1.0 m and 1.5 m telescopes at Mt. Lemmon, Arizona.

Following a strategy similar to that for the Pan-STARRS data search, we searched the archival data of all three telescopes dating back to 2018 January 1. We identified various sets of images covering the predicted ephemeris of 2I: multiple sets of images taken by the 0.7 m Catalina Schmidt, but all dated no later than 2018 December 15, and a set of four images taken by the 1.5 m telescope on 2019 March 1.

The comet was too faint (≳1.5 mag beyond the limit of the image) for the 2018 December and earlier images, while the 2019 March images reached a limiting magnitude that could be compatible with a faint detection of 2I. However, the position corresponding to the ZTF detections fell on a region heavily contaminated with background field sources. In three of the four frames, the object would have overlapped the point-source-function (PSF) of a field star. The remaining frame was also marginally affected by an even brighter nearby star, but it might show a slight enhancement that is within 1.5–2 pixels from the predicted position of 2I (Figure 5). A deeper stack of historical images obtained by Catalina with the same telescope reveals no background source at that position, down to a limiting magnitude much fainter than the individual frame. Unfortunately, the enhancement was extremely faint, and might be compatible with a noise feature enhanced by the tails of the bright star's PSF, making it difficult to draw a solid conclusion.

Zoom In Zoom Out Reset image size

The best-fit light curve as shown in Figure 6 revealed a shallow slope, K₁ = 5.3 ± 0.2, taking a linear phase coefficient β = 0.04 mag deg⁻¹. A shallow slope means that the activity stays largely constant as the comet approaches the Sun. Whipple (1978) shows that shallow brightening is common on the "dynamically new" solar system comets (with orbital periods P ≫ 10⁴ yr), which have ${K}_{1}\in (5,8)$), but uncommon on short-period comets and other long-period comets, which have K₁ ≳ 10, though his data set is dominated by small heliocentric distances with r_H < 5 au. In this respect, 2I is analogous to dynamically new comets in the solar system.

The K₁ = 5.3 slope seemingly deviates from the data beyond ∼8 au, preceded by what appears like a steep brightening phase, with ${K}_{1}\approx 29$, though additional pre-discovery observations are needed to verify it. If such a steep phase is real, it might indicate the onset of sublimation of cometary volatiles. The most compatible major volatile would be CO₂, which has an onset "knee" distance r_H = 13 au. Other cometary species, such as CH₃CN, HCN, and CH₃OH, have similar turn-on distances (Meech & Svoren 2004), but they have low abundances in solar system comets (Cochran et al. 2015). CO, another cometary volatile commonly found on dynamically new comets in solar system that has an onset distance r_H = 120 au, is not as compatible. However, it has been suggested for the case of 'Oumuamua that the volatiles may be buried beneath the surface and are only activated when the comet is much closer to the Sun, due to the time lag for the heat wave to penetrate to the depth of the ice (Fitzsimmons et al. 2018; Seligman & Laughlin 2018). Therefore, it is too early to exclude CO as the main driver of 2I's activity.

An alternative explanation of the activity is the exothermic crystallization of amorphous water ice, a mechanism that may be responsible for the activity of distant comets (Prialnik & Bar-Nun 1990). Amorphous ice forms below an environmental temperature of ∼130 K and is capable of trapping gas as they form, a phenomenon that has been observed in laboratory experiments (Bar-Nun et al. 1987; Jenniskens & Blake 1994), though the presence of amorphous ice is yet to be directly observed on comets. Depending on the illumination of the cometary nucleus, crystallization of amorphous ice on the surface can start around 6–12 au (Jewitt et al. 2017), consistent with the observed turn-on distance of 2I.

To gain deeper insight into the light curve, we tested a sublimation model (Meech & Svoren 2004) that computes the amount of gas sublimating from an icy surface exposed to solar heating to explore the activity. The total brightness within a fixed aperture combines radiation scattered from both the nucleus and the dust dragged from the nucleus in the escaping gas flow, assuming a dust-to-gas mass ratio of 1. We used a nucleus radius of 0.5 km (Jewitt & Luu 2019), assuming an albedo of 0.04 for the nucleus and a linear phase function of 0.04 mag deg⁻¹ for the nucleus and 0.02 mag deg⁻¹ for the coma typical of other comets (Meech & Jewitt 1987; Krasnopolsky et al. 1987), a nucleus density of 400 kg m⁻³ similar to that seen for comets 9P/Tempel 1, 103P/Hartley 2 (Thomas 2009), and 67P/Churyumov–Gerasimenko (Pätzold et al. 2016), a grain density of 800 kg m⁻³ (Fulle et al. 2016), and large (10–100 μm) grains (Fitzsimmons et al. 2019). Unsurprisingly, our model confirmed that the activity of 2I must be driven by a species more volatile than water, otherwise it would have been well below the detection limit of any of the surveys at ∼5 au. We also found that the differences between the shape of the sublimation curves for CO and CO₂ near 8 au is minimal, so it is impossible to distinguish between these volatiles without further pre-discovery observations.

The non-detection in 2018 November data can be used to constrain the size of the nucleus. The effective cross-section area for scattering can be calculated using

$$\begin{eqnarray}&&{C}_{{\rm{e}}}=\displaystyle \frac{\pi {\left(1\mathrm{au}\right)}^{2}}{p}{10}^{-0.4({M}_{r}-{m}_{\odot ,r})},\end{eqnarray} \tag{ 4 }$$

where p = 0.04 is the assumed optical geometric albedo of the nucleus (Lamy et al. 2004), ${m}_{\odot ,r}=-26.9$ is the apparent r-band magnitude of the Sun (Willmer 2018), and the absolute r-band magnitude (technically, the upper limit) M_r is defined as

$$\begin{eqnarray}&&{M}_{r}={m}_{r}-5\mathrm{log}{r}_{{\rm{H}}}{\rm{\Delta }}-{\rm{\Phi }}\alpha ,\end{eqnarray} \tag{ 5 }$$

with the variables following the same definitions for Equation (3). By inserting all the numbers, we have ${C}_{{\rm{e}}}\,\lt 140\,{\mathrm{km}}^{2}$ for the non-detection in 2018 November, at 8.5 au. This upper bound indicates that the radius of the nucleus is no larger than ∼7 km.

The size of the active area on the nucleus can be estimated with knowledge of the mass-loss rate of the comet and the mass flux of the activity-driving volatiles. The mass-loss rate of 2I can be estimated using the cross-section area of the dust and the speed of the dust flow. Using Equation (4), we derive the cross-section area of the coma to be ∼280–$360\,{\mathrm{km}}^{2}$, from March to October 2019, in which r_H decreases from 6.4 to 2.6 au. The mass-loss rate $\dot{M}$ can be calculated by

$$\begin{eqnarray}&&\dot{M}=\displaystyle \frac{4\rho \bar{a}{C}_{{\rm{e}}}}{3\tau },\end{eqnarray} \tag{ 6 }$$

where $\rho =800\,\mathrm{kg}\,{{\rm{m}}}^{-3}$ is the bulk density of the dust (Fulle et al. 2016), which is admittedly not yet constrained for interstellar comets, but we do not have reason to believe it is much different from solar system comets, and therefore have assumed it to be analogous to the latter, $\bar{a}={10}^{-4}\,{\rm{m}}$, the characteristic size of the dust, which is similarly assumed based on the observation of dynamically new solar system comets (e.g., Ye & Hui 2014; Jewitt et al. 2019), and τ is the timescale on which a dust particle moves out of the aperture, which can be estimated by $\tau =l/v\approx {10}^{5}\,{\rm{s}}$, where l is the linear length of the aperture at the comet, and $v\sim 1\,{\rm{m}}\,{{\rm{s}}}^{-1}$ is the ejection speed of the dust, taking the dust speed constrained by Guzik et al. (2019) and assuming a classic ∝a^−0.5 dependence. By inserting all the numbers, we obtain $\dot{M}\approx 10\,\mathrm{kg}\,{{\rm{s}}}^{-1}$, with the uncertainty around an order of magnitude mainly owing to the parameter $\bar{a}$ (which can vary by an order of magnitude among solar system comets; see Fulle 2004).

We then solve the energy balance equation for CO and CO₂ ice, which are likely to be responsible for 2I's activity following the discussion in Section 4.1, at the sub-solar point on the nucleus:

$$\begin{eqnarray}&&\displaystyle \frac{(1-{A}_{B}){L}_{\odot }}{4\pi {r}_{{\rm{H}}}^{2}}=\epsilon \sigma {T}^{4}+Z(T)L(T),\end{eqnarray} \tag{ 7 }$$

where A_B = 0.01 is the Bond albedo of the nucleus measured for 9P/Tempel 1 and 103P/Hartley 2 (Li et al. 2013a, 2013b), ${L}_{\odot }$ is the luminosity of the Sun,  = 0.9 is the infrared emissivity of the nucleus, σ is the Boltzmann constant, L(T) is the latent heat of the sublimation of the ice, and Z(T) is the mass flux. We solve L(T) using the model by Cowan & A'Hearn (1979),¹⁶ and obtain Z(T) ≈ 3 × 10⁻⁶–$2\times {10}^{-5}\,\mathrm{kg}\,{{\rm{m}}}^{-2}\,{{\rm{s}}}^{-1}$ from 6.4 to 2.6 au for CO, and 1 × 10⁻⁶–$9\times {10}^{-6}\,\mathrm{kg}\,{{\rm{m}}}^{-2}\,{{\rm{s}}}^{-1}$ for the same r_H span for CO₂. The active area required to support the mass-loss rate would then be

$$\begin{eqnarray}&&A=\displaystyle \frac{\dot{M}}{({M}_{{\rm{d}}}/{M}_{{\rm{g}}})Z(T)},\end{eqnarray} \tag{ 8 }$$

where ${M}_{{\rm{d}}}/{M}_{{\rm{g}}}$ is the dust-to-gas mass ratio, which is again unknown for interstellar comets. If we take M_d/M_g = 1 based on the measurement of long-period comet C/1995 O1 (Hale-Bopp) at similar heliocentric distance (Weiler et al. 2003), we have A = 0.5–$10\,{\mathrm{km}}^{2}$ over the time span between 2018 December and 2019 September, in line with the number derived by McKay et al. (2020) based on their observation in 2019 October ($1.7\,{\mathrm{km}}^{2}$). Taking the upper limit of the size of the nucleus derived in Section 4.2, the active fraction of 2I is >0.1% of the nucleus. This is in line with the known solar system comets, which have fractional active areas from a few 0.1% to >100% (Tancredi et al. 2006), though we note that most of these measured comets are short-period comets with H₂O as the driving species. However, we caution that the uncertainty in A is about 1–2 orders of magnitude when we consider the uncertainties in $\dot{M}$ and M_d/M_g, therefore the derived A and active fraction are highly uncertain.

The inclusion of the precovery data in the orbit estimation process provides more stringent constraints on the trajectory of 2I, but also introduces challenges to correctly modeling the dynamics. A gravity-only model of the orbit struggles to fit the data of 2019 March and earlier. In particular, the 2018 December position is rejected as an outlier (using the outlier rejection algorithm by Carpino et al. 2003). At this stage it is not entirely clear whether this behavior is caused by systematic errors in the bulk of the astrometric data, which were taken at low solar elongation and therefore at high airmass, or by non-gravitational accelerations. We note that non-gravitational accelerations were detected in the motion of 'Oumuamua, despite the lack of visible outgassing (Micheli et al. 2018).

Table 1 reports JPL solution 37, which fits all the precovery observations and uses non-gravitational forces assuming CO as the primary driver (Section 2.1), more consistent than H₂O with the photometric data. The non-gravitational model for CO₂ is not available at this point; but since both CO and CO₂ are more volatile than water and 2I was in the regime of both volatiles when it was observed, we believe that the two models should behave similarly over the fit span.

We also tested the rotation jet model (Chesley & Yeomans 2005), which computes the non-gravitational perturbations from a discrete number of jets, whose acceleration is averaged over a nucleus rotation. For the driver of the activity we again used CO. We considered two jets, a nearly polar one at $10^\circ$ of colatitude and a midlatitude one on the southern hemisphere at 135° of colatitude. Then, we scanned a raster for the spin pole's R.A. and decl., estimating the strengths of the two jets from the fit to the astrometry.

As shown in Figure 7, we find two minima for the pole's R.A. and decl.: (340°, +30°) and (205°, −55°). Because of the larger number of parameters, the jet model provides a better fit to the data than the Marsden et al. (1973) model. However, its reliability will need to be validated by its ability to make accurate predictions. Past experience has shown that the jet model can provide more accurate comet trajectory estimates (Farnocchia et al. 2016) and therefore the jet model solutions are worth consideration, especially as the observed arc extends into late 2019 and 2020.

Zoom In Zoom Out Reset image size