Observations of Disintegrating Long-period Comet C/2019 Y4 (ATLAS): A Sibling of C/1844 Y1 (Great Comet)

We took measurements of comet C/2019 Y4 in the NEXT and LDT images using an aperture of fixed linear radius $\varrho ={10}^{4}$ km projected at the distance of the comet from the optocenter. The equivalent apparent angular size of the aperture is always large enough such that the slight trailing of the comet in the NEXT data from 2020 January to February would not be a concern. Figure 3(a) shows our multiband lightcurve measurements as functions of time, in terms of time from the epoch of perihelion passage of C/2019 Y4 (t_p= TDB 2020 May 31.0). The comet apparently brightened on its way to perihelion in a continuous manner until $t-{t}_{{\rm{p}}}\gtrsim -70$ days, after which the downtrend in brightness was seen.

Zoom In Zoom Out Reset image size

We also show the color of the comet in terms of g − r and r − i color indices respectively in the left and right panels of Figure 4. Interestingly, while the r − i color index of the comet remained constant and Sun-like (i.e., ${\left(r-i\right)}_{\odot }=+0.12\,\pm 0.02;$ Willmer 2018) given the measurement uncertainties, the color across the g and r bands seems to indicate that the comet had a bluing trend from a reddish color (g − r ≈ 0.6, in comparison to the Sun's ${\left(g-r\right)}_{\odot }=+0.46\pm 0.04;$ Willmer 2018) since 2020 January, reached a dip at an epoch of ∼60 days preperihelion, when the comet appeared even bluer than the Sun (g − r ≈ 0.2), and began to be reddening afterwards. Accordingly, we argue that the bluing dip was caused by gas emission from a massive amount of previously buried fresh volatiles suddenly exposed to sunlight, indicative of a disintegration event in 2020 mid-March.

Zoom In Zoom Out Reset image size

We evaluated the intrinsic lightcurve of the comet by correcting the varying observing geometry and computed its absolute magnitude from the apparent magnitude using

$$\begin{eqnarray}&&{m}_{\lambda }\left(1,1,0\right)={m}_{\lambda }\left({r}_{{\rm{H}}},{\Delta },\alpha \right)-5{\rm{log}}\left({r}_{{\rm{H}}}{\Delta }\right)+2.5{\rm{log}}\phi \left(\alpha \right),\end{eqnarray} \tag{ 1 }$$

in which λ is the magnitude bandpass, r_H and Δ are the heliocentric and topocentric distances, respectively, both expressed in astronomical units, and $\phi \left(\alpha \right)$ is the phase function of the comet, approximated by the empirical Halley-Marcus phase function (Marcus 2007; Schleicher & Bair 2011). The resulting intrinsic lightcurve of C/2019 Y4 is plotted in Figure 3(b), from which we clearly notice that the lightcurve trend appears broadly the same as the one in Figure 3(a). Indeed the comet continuously brightened until ≲70 days preperihelion, thereafter followed by a conspicuous fading process in an intrinsic manner, which we think also acts as a piece of evidence that a disintegration event has occurred to comet C/2019 Y4.

The change in the intrinsic brightness is mostly related to the variation in the effective scattering geometric cross-section of the dust particles through

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

where C_e is the cross-section, ${m}_{\odot ,r}=-26.93$ is the apparent r-band magnitude of the Sun at the mean Earth–Sun distance r₀ $\approx$  1.5 × 10⁸ km (Willmer 2018), and p_r = 0.1 is the assumed value for the r-band geometric albedo of cometary dust (Zubko et al. 2017), as the true value remains unconstrained. The reason why we only focus on the r-band data here is that these images have less contamination from gaseous emission than the g-band ones do, and they have higher sensitivity than the i-band images do. By employing Equation (2), we estimated the change in the effective cross-section within the fixed photometric aperture during the brightening part ($-140\lesssim t-{t}_{{\rm{p}}}\lesssim -80$ days) of the lightcurve to be ${\rm{\Delta }}{C}_{{\rm{e}}}=\left(9.9\pm 0.5\right)\times {10}^{2}$ km², corresponding to an average growth rate in the effective cross-section of $\left\langle {\dot{C}}_{{\rm{e}}}\right\rangle =\left(2.0\pm 0.1\right)\times {10}^{2}$ m² s⁻¹. Assuming that the increased cross-section consists of dust grains with mean radius $\bar{{\mathfrak{a}}}$ and bulk density ρ_d, the average net mass-loss rate within the photometric aperture is then given by

$$\begin{eqnarray}&&\left\langle {\dot{M}}_{{\rm{d}}}\right\rangle =\displaystyle \frac{4}{3}{\rho }_{{\rm{d}}}\bar{{\mathfrak{a}}}\left\langle {\dot{C}}_{{\rm{e}}}\right\rangle .\end{eqnarray} \tag{ 3 }$$

The product of $\bar{{\mathfrak{a}}}$ and ρ_d is inversely proportional to the β parameter (0.03 ≲ β  ≲  0.1, see Section 2.3). Substitution into Equation (3) gives us $\left\langle {\dot{M}}_{{\rm{d}}}\right\rangle \approx 4\pm 2\,\mathrm{kg}$ s⁻¹ for C/2019 Y4 during the observed brightening period.

An approach to constrain the nucleus size of C/2019 Y4 is to estimate the minimum active surface that would be needed to supply the mass-loss rate during the brightening process, provided that the activity was all driven by sublimation of water (H₂O) ice. The corresponding lower bound to the nucleus size can then be estimated from

$$\begin{eqnarray}&&{R}_{{\rm{n}}}\gtrsim \sqrt{\displaystyle \frac{\left\langle {\dot{M}}_{{\rm{d}}}\right\rangle }{\pi {f}_{{\rm{s}}}}}.\end{eqnarray} \tag{ 4 }$$

The equilibrium sublimation mass flux of H₂O gas would be 1.2 × 10⁻⁶  ≲ f_s  ≲  1.2 × 10⁻⁴ kg s⁻¹ m⁻² at the range of the heliocentric distances during the time span (1.8 ≲ r_H  ≲  2.7 au). Inserting numbers in, we obtain R_n  ≳  60 m for the nucleus size of the comet.

An upper limit to the nucleus size of C/2019 Y4 could have been derived from our detection of the nongravitational acceleration (Section 3.1). However, strictly speaking, the nongravitational effect is only applicable to the barycenter of the unresolved fragments, rather than an intact nucleus. The equivalent bulk density of the barycenter should be much lower than those of typical cometary nuclei (e.g., 533 ± 6 kg m⁻³ for 67P/Churyumov–Gerasimenko; Pätzold et al. 2016), to a degree we cannot firmly constrain. We thus posit that applying this approach by assuming a typical bulk density for cometary nuclei is no longer valid.

C/2019 Y4 has been unambiguously cometary since our earliest observation in 2020 January. Its central optocenter had been strongly condensed until 2020 April 5, after which it started to become even more diffuse and elongated (see Figure 1).⁵ The morphological change, altogether with ongoing observations from NEXT in which we see multiple optocenters in the coma, strongly indicates that the cometary nucleus has split into multiple pieces of fragments (Ye & Zhang 2020). This feature has been confirmed by Steele et al. (2020) and Lin et al. (2020). A detailed study about the disintegration and follow-up observations will be presented in another paper in preparation.

Physical properties of cometary dust can be revealed by studying the morphology of cometary dust tails (e.g., Fulle 2004). To better understand comet C/2019 Y4, here we adopted the classical syndyne–synchrone computation (e.g., Finson & Probstein 1968). A syndyne line is loci of dust grains that are subject to the β parameter in common, which is the ratio between the solar radiation pressure force and the gravitational force of the Sun, and is inversely proportional to the product of the dust bulk density and the grain radius, ${\rho }_{{\rm{d}}}{\mathfrak{a}}$, but are freed from the nucleus at various release epochs. Dust grains that are driven by different values of β and are released from the nucleus at a common epoch constitute a synchrone line. In the syndyne–synchrone approximation, the ejection velocity of dust grains is ignored.

We focused on the r-band images taken since 2020 March from NEXT, in which the dust tail of C/2019 Y4 was recorded the most clearly and the optocenter appeared untrailed. In Figure 5, we plot four examples of the syndyne–synchrone grid computation. Here for better visualization, only sparse syndyne–synchrone grids are presented. A much denser grid was used to visually determine the β parameter and the dust release epochs about which the dust tail appeared to be symmetrical. For data from some single observed epoch, we find it difficult to judge whether the symmetry of the dust tail appeared more like a syndyne or synchrone. However, with results from multiple observed epochs, we realized that the dust tail of the comet can be better approximated by syndyne lines with 0.03 ≲  β  ≲  0.1, because the range of the β parameter remained roughly the same, whereas the dust release epoch would keep changing, were the tail closer to a synchrone line. This finding is consistent with the fact that the comet has been active protractedly since the discovery. Assuming a typical value of ρ_d = 0.5 g cm⁻³ for the bulk density of cometary dust of C/2019 Y4, we find the grain radius to be $10\lesssim \bar{{\mathfrak{a}}}\lesssim 40$  μm, fully within the known dust-size range of other long-period comets (Fulle 2004).

Zoom In Zoom Out Reset image size

As the syndyne–synchrone computation does not unveil the ejection speed of the observed dust grains of the comet, we estimate this quantity, denoted as V_ej, by means of measuring the apparent length of the sunward extent to the dust coma ℓ. The two quantities are connected by the following relationship

$$\begin{eqnarray}&&{V}_{{\rm{ej}}}=\displaystyle \frac{\sqrt{2\beta {\mu }_{\odot }{\Delta }\tan \ell \sin \alpha }}{{r}_{{\rm{H}}}},\end{eqnarray} \tag{ 5 }$$

where μ_⊙ = 3.96 × 10⁻¹⁴ au³ s⁻² is the heliocentric gravitational constant. We found ℓ ≈ 10'' on 2020 March 1, ∼15'' on March 13, and ∼25'' on March 31. Substituting, Equation (5) yields V_ej ≈ 30 m s⁻¹ at the beginning of 2020 March, ∼50 m s⁻¹ around halfway, and further increased to ∼80 m s⁻¹ at the end of the month for dust grains of β ∼ 0.1. This finding is similar to what has been identified for other long-period comets (e.g., Moreno et al. 2014).

We performed astrometry of comet C/2019 Y4 in the r-band images from LDT and NEXT by exploiting AstroMagic⁶ and codes developed by D. Tholen with the Gaia DR2 catalog (Gaia Collaboration et al. 2018). In this step, we realized that the tailward bias of the astrometric measurements was readily conspicuous in most of the observed images before early 2020 April, after which the comet apparently lost the central condensation visibly and so no astrometry was measured.⁷ To make sure that our orbit determination will not be skewed by the tailward bias in the astrometry, we performed least-square linear fits to the centroids in the J2000 equatorial east–west and decl. directions, respectively, as functions of the aperture size (Figure 6). The zero-aperture astrometry was then obtained, with its uncertainties propagated from the centroiding errors. In addition to our astrometric measurements of C/2019 Y4, we also included the astrometric observations of the comet from station T12 (Tholen NEO Follow-Up at the University of Hawai'i 2.24 m telescope) available from the MPC Observations Database,⁸ which have been corrected for the tailward bias as well (D. Tholen, from whom we obtained the corresponding measurement errors through private communication). The other astrometric observations available from the MPC Observations Database had to be discarded, because we do not think that they are zero-aperture astrometry of the comet, and no astrometric uncertainty is available either.

Zoom In Zoom Out Reset image size

We employed the orbit determination code FindOrb⁹ developed by B. Gray, which incorporates gravitational perturbation from the eight major planets, Pluto, the Moon, and the most massive 16 main-belt asteroids. The code applies post-Newtonian corrections, and uses the planetary and lunar ephemerides DE431 (Folkner et al. 2014). Initially we attempted to fit a purely gravitational orbit to the astrometric observations, however, the resulting astrometric residuals exhibit an obvious systematic trend in observations beyond 3σ from both the beginning and the end of the observed arc. The mean rms residual of the fit is 0266 from 104 observations in total. However, after we included the radial and transverse nongravitational parameters A₁ and A₂, first introduced by Marsden et al. (1973), in a nongravitational force model in which the nongravitational acceleration of the nucleus is assumed to be proportional to the mass flux of hemispherical H₂O-ice sublimation (M.-T. Hui & D. Farnocchia 2020, in preparation), the trend can be completely removed and the mean rms residual of the fit shrinks to 0126, only roughly the half of the one in the gravity-only solution. Adding the normal component of the nongravitational parameter A₃ does not improve the orbital fit. Nor is the obtained A₃ statistically significant. Furthermore, only in very few cases has A₃ been determined meaningfully (Yeomans et al. 2004), suggesting that A₃ plays a less significant role in comparison to A₁ and A₂. Therefore, we opted not to include A₃ but only A₁ and A₂. Our best-fit nongravitational solution to the orbit of C/2019 Y4 as well as the associated details are summarized in Table 2, where we can see that the radial nongravitational parameter of the comet, ${A}_{1}=\left(+2.25\pm 0.13\right)\times {10}^{-7}$ au day⁻², is particularly enormous among the whole comet population, but is by no means unseen among disintegrated comets, e.g., ${A}_{1}\,=\left(+1.21\pm 0.12\right)\times {10}^{-6}$ au day⁻² for C/2015 D1 (SOHO) by Hui et al. (2015) and ${A}_{1}\,=\left(+1.74\pm 0.11\right)\times {10}^{-7}$ au day⁻² for C/2017 E4 (Lovejoy) by JPL Horizons. The obtained transverse nongravitational parameter, ${A}_{2}=\left(-3.1\pm 1.0\right)\times {10}^{-8}$ au day⁻², is far less significant than its radial counterpart A₁ by almost an order of magnitude, and yet is nevertheless typical in the context of disintegrated comets, e.g., ${A}_{2}=\left(-1.55\pm 0.09\right)\times {10}^{-8}$ au day⁻² for C/1999 S4 (LINEAR) and ${A}_{2}=\left(+6.2\pm 0.8\right)\times {10}^{-8}$ au day⁻² for C/2010 X1 (Elenin), both computed by JPL Horizons.

Our result is similar to the nongravitational solutions to the orbit of C/2019 Y4 by the MPC¹⁰ (A₁ = +2.6 × 10⁻⁷ au day⁻¹ and A₂ = −2.9 × 10⁻⁸ au day⁻¹, no uncertainties given) and by JPL Horizons (${A}_{1}=\left(+2.86\pm 0.17\right)\times {10}^{-7}$ au day⁻² and ${A}_{2}=\left(-0.9\pm 1.2\right)\times {10}^{-8}$ au day⁻²), despite that most of their used astrometric observations are likely uncorrected for the tailward bias, and the old nongravitational force model by Marsden et al. (1973) was adopted. We think that in this specific case, although the tailward bias is obviously present in their used data from individual nights, it has been fortuitously negated by the varying observing geometry (Table 1), as well as by the enormous nongravitational effect.

In order to investigate the dynamical relationship between C/1844 Y1 and C/2019 Y4, we also derived a gravity-only orbit for the latter, using a shorter observed arc during which we found no statistically significant nongravitational effect, because we prefer that the enormity of the nongravitational effect is unlikely to be characteristic of the complete orbit, but reflects the ongoing disintegration in the current apparition. For this purpose, we only wanted to include data points from the period when the nongravitational force has not yet played an important role. The earliest astrometric data from T12 were always included. We tested with both gravity-only and nongravitational force models, and checked astrometric residuals of the astrometry, the significance of the radial nongravitational parameter A₁, and the mean rms residuals of the fits. What we found is that if any astrometric observations from 2020 March 13 and thereafter are included, the nongravitational solution to the orbit of C/2019 Y4 improves the fit considerably in comparison to the gravity-only version. Thus, by discarding all of the astrometric observations starting from 2020 March 13, we obtained the final version of the gravity-only solution to the orbit of the comet. The information is summarized in Table 2 as well.

For C/1844 Y1, neither the MPC nor JPL Horizons give the uncertainty information of its orbital elements. Thus, we extracted the topocentric astrometry from Bond (1850). Only observations with both R.A. and decl. measurements available from the same epochs were used and fed into FindOrb. The equatorial coordinates were precessed from epoch 1845.0 to 2000.0. Based upon the mean residual of the preliminary orbital solution, we downweighted all of the observations by an equal uncertainty of 15''. Ten (out of 80 in total) observations with astrometric residuals ≳3σ were rejected as outliers. Our best-fit orbital solution for C/1844 Y1, which we found to be in agreement with the published one by the MPC and JPL Horizons at the 1σ level, is presented in Table 2, together with our solutions for C/2019 Y4.

The similarity between the orbits of C/1844 Y1 and C/2019 Y4 obviously hints at a possible genetic relationship between the two comets that they are likely two components that split from a common progenitor. Our primary goal is to investigate when the split event between the comet pair most likely took place and how large the separation speed was.

We adopted a simplistic two-body dynamical model for the split event as follows. At some epoch t_frg, the progenitor of the comet pair experienced a split event, during which two major components—C/1844 Y1 and C/2019 Y4—were produced. The gravitational interaction between the pair was neglected, so was gravitational perturbation from the major planets in the solar system, as this effect is generally relatively unimportant (Sekanina & Kracht 2016). To make sure that this choice is valid to both C/1844 Y1 and C/2019 Y4, we created 1,000 clones for either of the comets based on the obtained best-fit gravity-only orbital elements and the covariance matrices. Then we utilized MERCURY6 (Chambers 1999) to integrate the clones together with the nominal orbits backward to 7 kyr ago, well past the previous perihelia that took place ∼5 kyr ago. What we found is that neither of the two comets had close approaches to any of the major planets since their last perihelion passages, validating the choice of neglecting planetary perturbation as an approximation. Our task is essentially equivalent to identifying how the orbits of C/1844 Y1 and C/2019 Y4 intersect and when the pair were both at the intersection point.

The major difference in the orbital elements of comets C/1844 Y1 and C/2019 Y4 lies in the epochs of their perihelion moments, whereas all other elements are not distinct (Table 2). Therefore for the simplistic dynamical model, we concentrated on the time interval of the perihelion epochs of the comet pair. As the orbit of C/2019 Y4 is more accurate, we used its elements, including periapsis distance q, eccentricity e, inclination i, longitude of ascending node Ω, and argument of periapsis ω at epoch TDB 1700 January 01.0, referenced to the solar system barycentric ecliptic J2000.0, for the sake that the barycentric orbital elements of both comets remain largely constant between two consecutive periapsis returns. In this reference system, the differences between all of the orbital elements but the periapsis epoch are even within the corresponding 1σ orbital element uncertainties of C/1844 Y1 (Table 3). We focus on how the separation velocity between the two comets and the split epoch should be can produce a difference in periapsis epochs of $\left|{\rm{\Delta }}{t}_{{\rm{p}}}\right|\approx 175.5\,\mathrm{yr}$ in the following.

Separation speeds between fragments of split comets are found to be in a range of 0.1 ≲ V_sep  ≲  10 m s⁻¹ (Boehnhardt 2004, and citations therein), of which only the in-plane component can alter the periapsis epoch, but not other orbital elements (Sekanina & Kracht 2016). We investigated the influence of the separation velocity in terms of radial, transverse, and normal (RTN) components upon the difference between the periapsis epochs of two fragments in the simplistic two-body dynamical model. The RTN coordinate system has its origin at the primary fragment, with the radial axis pointing away from the solar system barycenter, the normal axis directed along the total angular momentum of the primary fragment, and the transverse axis constructed to form a right-handed orthogonal system. At some split epoch t_frg, the state vector of the primary fragment was computed from the orbital elements. The total velocity of the secondary fragment was updated by adding the separation velocity V_sep to the velocity of the primary fragment. Thereby a new state vector was obtained, which was then converted into the orbital elements for the secondary fragment.

We searched for conditions that should be satisfied for the split between C/1844 Y1 and C/2019 Y4 from the latest possible orbital revolution, such that the observed gap between the periapsis moments of the comet pair can be achieved. Given the periodicity of their orbits, there is an indefinite number of desired conditions from more than a single orbital revolution ago. The reason why we consider only the latest possible orbital revolution for the split event is as follows. First, comet C/2019 Y4 is fragmenting in the current apparition, indicative of the nucleus as a loosely bound aggregate as typical cometary nuclei (e.g., Weissman et al. 2004, and citations therein). Second, given that the number of the observed split events of long-period comets was ∼30 among a total number of ∼10³ long-period comets discovered in the past ∼150 yr (identified using the JPL Small-Body Database Search Engine), we can estimate a lower limit to the splitting rate as ∼2% per century for each long-period comet, which is comparable to the one for short-period comets (∼3% per century per comet; Boehnhardt 2004). Assuming that all long-period comets behave alike, within one orbit (∼5 kyr) around the Sun, their progenitor (or any of its descendants) would experience at least one split event. Therefore, limiting our search for the split between comets C/1844 Y1 and C/2019 Y4 only within a time frame not too much more than one orbital revolution in the past is a reasonable confinement.

In Figure 7, we plot the change in the periapsis epoch Δt_p as a function of the RTN components of the separation velocity and the split epoch. We can learn that only if the split event that produced C/1844 Y1 and C/2019 Y4 occurred around the previous perihelion passage, which was ∼5 kyr ago, and the in-plane component of the separation velocity between the pair was ≳1 m s⁻¹, can the observed difference between the periapsis epochs of the pair be caused by separation speeds within the known range of split comets. Another conclusion we can draw is that the out-of-plane component of the separation velocity between the pair alone cannot bring about the observed periapsis epoch difference, which is similar to what Sekanina & Kracht (2016) found for another comet pair C/1988 F1 (Levy) and C/1988 J1 (Shoemaker–Holt).

Zoom In Zoom Out Reset image size

To ensure that the planetary perturbation will not drastically alter the conclusion we drew based on the simplistic two-body model, we also employed our code, which includes planetary perturbation and has been utilized to analyze the split event of active asteroid P/2016 J1 (Pan-STARRS; Hui et al. 2017), to find a best-fit nonlinear least-squared solution to the split parameters between C/1844 Y1 and C/2019 Y4. We treated C/2019 Y4 as the major component, and used its nominal orbit. The major difference here from the application in Hui et al. (2017) is that, instead of fitting a list of topocentric astrometric observations, we fitted the heliocentric orbital elements of C/1844 Y1 (Table 2) with the associated covariance matrix from Section 3.1, and minimized the following quantity

$$\begin{eqnarray}&&{\chi }^{2}\left({t}_{\mathrm{frg}},{V}_{{\rm{R}}},{V}_{{\rm{T}}},{V}_{{\rm{N}}}\right)={\rm{\Delta }}{{\boldsymbol{E}}}^{{\rm{T}}}{\boldsymbol{W}}{\rm{\Delta }}{\boldsymbol{E}},\end{eqnarray} \tag{ 6 }$$

where ΔE is the orbital element residual vector and W is the weight matrix determinable from the covariance matrix of the orbital elements ${\boldsymbol{E}}=\left(q,e,i,{\rm{\Omega }},\omega ,{t}_{{\rm{p}}}\right)$ of C/1844 Y1 (e.g., Milani & Gronchi 2010). Using different initial guesses, we soon realized that unless adopting an exhaustive and extensive search, which will be extremely time consuming, the code would converge to different solutions, indicative of the existence of multiple local minima. We present two of the best-fitted solutions we obtained in Table 4. Although there is no definite solution to the split parameters between C/1844 Y1 and C/2019 Y4 from the N-body dynamical model, we think that the conclusion based on the simplistic two-body model remains valid with inclusion of planetary perturbation, because the best-fit solutions all have in-plane components of the separation velocity ≳1 m s⁻¹ and the split epochs around the previous perihelion of C/2019 Y4, at heliocentric distance r_H ≲ 10 au.

Although the fragmentation mechanism that led to the disruption between the pair remains unclear, we can safely rule out the possibility of tidal disruption, because the perihelion distance is larger than the Roche radius of the Sun for comets (a few solar radii) by at least an order of magnitude, and there was no close approach to any of the major planets in the time frame we investigated either. We postulate that possible fragmentation mechanisms include rotational instability due to anisotropic mass loss, excessive internal thermal stress, or gas pressure (Boehnhardt 2004, and and citations therein) overwhelming differential stress due to dynamic sublimation pressure (Steckloff et al. 2015). A detailed discussion of the potential physical mechanisms that caused the observed ongoing disintegration event of C/2019 Y4 will be presented in another paper in preparation.