C/2014 UN₂₇₁ (Bernardinelli-Bernstein): The Nearly Spherical Cow of Comets

The DES astrometry is mapped to Gaia DR2 (Gaia Collaboration et al. 2018) using the astrometric model presented in Bernstein et al. (2017). ⁵³ All distortions due to the telescope, instrument, and detectors are known to ≈1 mas rms, and the color-dependent effects (differential chromatic refraction in the atmosphere and lateral color distortions) are corrected using the object's mean g − i color. The position uncertainties for the DES exposures (2014–2018) include the shot noise from each detection, as well as an anisotropic contribution from the atmospheric turbulence (Bernardinelli et al. 2020). For the VISTA (2010) and CFHT (2014) images, we retrieve detrended images from the archives, remeasure the positions using SExtractor windowed centroiding, and produce a polynomial astrometric solution in the vicinity of BB by referencing nearby stars from the Gaia DR2 catalog. Astrometry for the PanSTARRS1 (PS1; 2014–2019) exposures is extracted using Gaussian fits to the comet and field stars, with the latter referenced to DR2. One PS1 exposure with a barely detectable signal-to-noise ratio (S/N) is a >3σ outlier from the orbit fit and excluded from further consideration.

We determine the object's orbit using the method of Bernstein & Khushalani (2000), and we do not include nongravitational forces in the orbit fit, as these have not been detected for BB yet. The orbital elements and derived uncertainties are presented in Table 1 and yield χ²/dof = 116.5/96. Considering only DES observations yields a consistent orbit with ≈1.5× larger uncertainties and χ²/dof = 66.4/58. The semimajor axis and inclination of the incoming orbit are 20,200 au and 955, respectively, fully characteristic of Oort cloud membership. Perihelion of 10.95 au will be reached on 2031 January 21. Forward integration of the orbit for ≈50 yr shows that the semimajor axis will be increased by 40% after this perihelion.

Table 1. Osculating Barycentric Orbital Elements

Epoch	a (au)	e	i	Ω	ω	tperi (JD)
1950	20, 200 ± 130	0.999458(4)	954663	1900029	3262793(3)	2,462,887.94(4)
2100	28, 070 ± 170	0.999610(2)	954606	1900093	3262438(4)	2,462,887.87(9)

Note. Elements are given at epochs before and after the current passage through the realm of the giant planets, assuming only gravitational forces. Uncertainties in the last digit are given in parentheses where they are sufficiently large.

Download table as:  ASCIITypeset image

It is of substantial interest to determine whether BB has been appreciably warmed on previous perihelion passages. We study the past dynamics of BB using a numerical procedure similar to that of Królikowska & Dybczyński (2018), and also by analytic approximations. We will consider perturbations by the Galactic tidal tensor ${ \mathcal G }$, assumed (as in Heisler & Tremaine 1986) to be diagonal in the Galactic frame rotating with the Sun, where $\hat{{\boldsymbol{x}}}$ points to the Galactic center and $\hat{{\boldsymbol{z}}}$ points to the north Galactic pole. This translates to a contribution to the Hamiltonian of the system in the form (Fouchard 2004)

$$\begin{eqnarray}&&{{ \mathcal H }}_{{\rm{G}}}=\displaystyle \frac{1}{2}{\boldsymbol{x}}\cdot { \mathcal G }\cdot {\boldsymbol{x}}={{ \mathcal G }}_{1}\displaystyle \frac{{x}^{2}}{2}+{{ \mathcal G }}_{2}\displaystyle \frac{{y}^{2}}{2}+{{ \mathcal G }}_{3}\displaystyle \frac{{z}^{2}}{2}.\end{eqnarray} \tag{ 1 }$$

We adopt the nominal Oort constants (Oort 1927) A = 15.1 and B = −13.4 km s⁻¹ kpc⁻¹ from Li et al. (2019) and a local stellar density ρ₀ = 0.15 M_☉ pc⁻³ from Vokrouhlický et al. (2019), so we have ${{ \mathcal G }}_{1}\equiv -(A-B)(3A+B)=-9.49\,\times {10}^{-16}$, ${{ \mathcal G }}_{2}\equiv {(A-B)}^{2}\,=\,8.48\times {10}^{-16}$, and ${{ \mathcal G }}_{3}\equiv 4\pi \mu {\rho }_{0}\,-2({B}^{2}-{A}^{2})=8.59\times {10}^{-15}\,{\mathrm{yr}}^{-2}$. The angular velocity of the Sun is Ω₀ ≡ B − A = −28.5 km s⁻¹ kpc⁻¹ = −2.91 × 10⁻⁸ yr⁻¹.

The numerical approach is to integrate, backward in time, clones of the orbit solution sampled from the state vector covariance matrix. We use the WHFast (Wisdom & Holman 1991; Rein & Tamayo 2015) integrator of REBOUND (Rein & Liu 2012) and include the giant planets as active perturbers, as well as the effects of the Galactic tide using REBOUNDx (Tamayo et al. 2020). Figure 1 shows the histogram of the previous orbit perihelion distance and time from numerical integration of the sampled orbits, which are near 18.2 au and 3.41 Myr, respectively.

Zoom In Zoom Out Reset image size

Analytically, we calculate the change in angular momentum Δ L imparted by the tidal torque over a full orbit in the limit where the orbit is fully radial, e → 1. Near this limit, L² = 2kq, where k is the barycentric gravitational constant 1.0014 GM_⊙, and q is the perihelion. We define $\hat{{\boldsymbol{e}}}$ as the unit vector toward perihelion (inverse of aphelion direction). The Born approximation then yields ${\rm{\Delta }}{\boldsymbol{L}}=5\pi \sqrt{\tfrac{{a}^{7}}{k}}\hat{{\boldsymbol{e}}}\times ({ \mathcal G }\cdot \hat{{\boldsymbol{e}}}).$ This yields a previous perihelion of 18.3 au for the nominal orbit, in agreement with the numerical integration.

The ascending node of the previous passage has r_h ≈ 20 au in the numerical integrations, suggesting the possibility that an encounter with Uranus had significantly altered the orbit, and some earlier perihelion had been <18 au. The chances of a sufficiently close Uranus encounter are very low, however. If BB encountered Uranus at impact parameter b with relative speed v at heliocentric distance r_h , then the impulse was I = 2GM_U /bv, and the maximum change in specific angular moment would be ∣Δ L ∣ ≤ Ir_h . Since BB is in a near-parabolic orbit that is nearly perpendicular to the ecliptic, we have

$$\begin{eqnarray}&&\displaystyle \frac{| {\rm{\Delta }}{\boldsymbol{L}}| }{| {\boldsymbol{L}}| }\leqslant \displaystyle \frac{{{Ir}}_{h}}{\sqrt{2{{GM}}_{\odot }q}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\leqslant \,\displaystyle \frac{2{{GM}}_{U}{r}_{h}}{b}\sqrt{\displaystyle \frac{{r}_{h}}{2{{GM}}_{\odot }}}\displaystyle \frac{1}{\sqrt{2{{GM}}_{\odot }q}}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{{M}_{U}}{{M}_{\odot }}\displaystyle \frac{{r}_{h}}{b}\sqrt{\displaystyle \frac{{r}_{h}}{q}}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\approx 5\times {10}^{-5}\displaystyle \frac{20\,\mathrm{au}}{b}.\end{eqnarray} \tag{ 5 }$$

For this perturbation of the angular momentum (and perihelion) to be as large as 10% requires b < 0.01 au. The chances that Uranus is this close to BB's previous node are <0.01% given the 120 au circumference of Uranus's orbit.

We may also use the impulse approximation to assess the angular momentum imparted by passages of stars close to the Sun during the previous orbit. For a star with mass M_⋆ with a closest approach to the Sun at point b and velocity v , while the comet is at position r , the angular momentum imparted is

$$\begin{eqnarray}&&{\rm{\Delta }}{\boldsymbol{L}}=\displaystyle \frac{2{{GM}}_{\star }}{v}{\boldsymbol{r}}\times \left[{\boldsymbol{b}}\left(\displaystyle \frac{1}{| {{\boldsymbol{b}}}^{{\prime} }{| }^{2}}-\displaystyle \frac{1}{| {\boldsymbol{b}}{| }^{2}}\right)+\hat{{\boldsymbol{v}}}\displaystyle \frac{{\boldsymbol{r}}\cdot \hat{{\boldsymbol{v}}}}{| {{\boldsymbol{b}}}^{{\prime} }{| }^{2}}\right],\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\boldsymbol{b}}}^{{\prime} }={\boldsymbol{b}}-{\boldsymbol{r}}+\hat{{\boldsymbol{v}}}({\boldsymbol{r}}\cdot \hat{{\boldsymbol{v}}}).\end{eqnarray} \tag{ 7 }$$

We use the list of reliable close stellar encounters (b < 1 pc) derived from the Gaia DR2 catalog by Bailer-Jones et al. (2018), restricted to those with perihelion times −4 Myr < t_ph < 0 Myr relative to present. We updated the stellar parameters for each star to the values and uncertainties given in Gaia EDR3 (Gaia Collaboration et al. 2021). This removes a very strong perturber from this list (DR2 955098506408767360) and some others, leaving eight potential encounters over BB's previous orbit (see also Bobylev & Bajkova 2020). We assume linear motion for the perturbing stars and sample from the Gaia uncertainties. The net effect of these encounters is to slightly decrease L , i.e., to raise the previous perihelion by an amount that is well below the effect of the Galactic tide. The left panel of Figure 1 shows the histograms of previous perihelion distances derived from the analytic approximations for the Galactic tide alone (orange; sampling from BB's orbital uncertainties) and for the Galactic tide plus stellar encounters (green).

The conclusion, which is robust to the details of the tidal model or these eight stars' dynamics, is that the previous passage of BB was further from the Sun than the current one. Indeed, under the tidal model, the perihelion has been getting smaller with each successive passage for many orbits into the past. We conclude that BB is a "new" comet in the sense that there is no evidence for a previous approach closer than 18 au to the Sun since ejection into the Oort cloud. Indeed, the observations of BB may be the first ever obtained for an object that has not crossed within Uranus's orbit since the solar system formation epoch. It remains true, however, that our knowledge of stellar encounters is incomplete, and it is possible that some yet-unknown star's passage could have lifted BB's perihelion from a lower value to its present one.

We measure the flux in each DES image using scene-modeling photometry, similar to Brout et al. (2019). We define a target region around each detection of 272 × 272 pixels (at 0264 pixel⁻¹) and simultaneously fit a model for the object's flux and the background sources to all DES images from the same filter in this region of the sky. The background is modeled as a grid of point sources that is present in all images, while the object is modeled as a point source present in only the detection image. Each point source is convolved with the point-spread function (PSF; see Jarvis et al. 2021 for a detailed description of the DES PSF model) of each pixel location in each exposure. This procedure also allows us to measure fluxes in exposures in which there is no detection of the object but the orbit indicates its presence. Thus, as seen in Table 2, there are DES images having photometry but no useful astrometric data. The resultant fluxes and errors are rigorously correct for an unresolved image, essentially using the central few arcseconds' signal, and thus insensitive to any coma that does not have a central concentration. Flux calibration for all DES exposures is determined to mmag precision, as described in the Dark Energy Survey Collaboration (2021).

Table 2. Observational Data for C/2014 UN₂₇₁

UTC	R.A.	Decl.	Error	Source	Mag	Band	Helio. Dist.	Geo. dist.
							(au)	(au)
2010-11-15.238223	01:19:36.7532	−30:42:27.260	0105	VISTA	22.47 ± 0.26	z	34.07	33.52
2014-08-14.599780	01:42:23.9960	−35:36:37.600	0078	PS1	23.10 ± 0.20*	i	29.28	28.74
2014-08-28.601410	01:41:25.1473	−36:02:11.803	004	CFHT	22.39 ± 0.09	r	29.23	28.59
2014-08-28.606693	01:41:25.1180	−36:02:12.397	009	CFHT	23.05 ± 0.08	g	29.23	28.59
2014-08-28.609338	01:41:25.1015	−36:02:12.632	006	CFHT	22.14 ± 0.17	i	29.23	28.59
2014-08-28.612016	01:41:25.0891	−36:02:13.023	004	CFHT	22.69 ± 0.12	r	29.23	28.59
2014-10-20.294348	01:34:35.0187	−37:14:46.144	0068	DES	22.62 ± 0.10	r	29.05	28.36
2014-11-04.122035	01:32:21.0290	−37:23:21.138	0050	DES	22.09 ± 0.12	z	28.99	28.41
2014-11-14.258286	⋯	⋯	⋯	DES	21.93 ± 0.46	Y	28.96	28.46
2014-11-15.251906	01:30:47.0392	−37:25:39.196	0047	DES	22.32 ± 0.16	z	28.95	28.47
2014-11-18.235595	01:30:23.4255	−37:25:40.010	0040	DES	22.52 ± 0.05	r	28.94	28.49
2014-11-18.238347	01:30:23.4067	−37:25:39.984	0041	DES	22.22 ± 0.07	i	28.94	28.49
2014-11-27.221480	01:29:17.5832	−37:24:12.514	0113	DES	⋯	g	28.91	28.55
2014-12-11.115676	01:27:54.8295	−37:17:49.638	0084	DES	22.43 ± 0.09	i	28.86	28.67
2014-12-11.144452	⋯	⋯	⋯	DES	22.57 ± 0.62	Y	28.86	28.67
2015-01-09.107010	01:26:35.3245	−36:51:38.750	0130	DES	22.65 ± 0.30	z	28.76	28.94
2015-08-11.612980	01:47:05.4300	−37:26:12.860	0102	PS1	23.33 ± 0.22*	i	27.99	27.49
2015-08-17.357470	01:46:47.5049	−37:37:39.814	0053	DES	22.64 ± 0.08	r	27.97	27.42
2015-08-17.358826	01:46:47.5170	−37:37:39.886	0080	DES	22.99 ± 0.09	g	27.97	27.42
2015-08-24.344520	01:46:18.6053	−37:51:28.614	0060	DES	23.12 ± 0.09	g	27.94	27.35
2015-08-24.347311	⋯	⋯	⋯	DES	22.53 ± 0.07	r	27.94	27.35
2015-08-24.348700	⋯	⋯	⋯	DES	22.38 ± 0.08	i	27.94	27.35
2015-09-01.306625	⋯	⋯	⋯	DES	22.46 ± 0.54	Y	27.92	27.27
2015-09-02.377866	01:45:30.1964	−38:08:52.304	0085	DES	21.94 ± 0.28	Y	27.91	27.26
2015-09-13.390509	01:44:15.8625	−38:28:51.622	0073	DES	22.47 ± 0.10	i	27.87	27.17
2015-10-06.274354	01:40:59.7742	−39:03:19.836	0147	DES	22.18 ± 0.17	z	27.79	27.08
2015-11-20.228513	01:33:55.0781	−39:29:32.567	0067	DES	22.18 ± 0.14	z	27.63	27.21
2015-11-20.235561	01:33:55.0098	−39:29:32.478	0058	DES	22.50 ± 0.11	i	27.63	27.21
2016-01-11.093863	01:29:57.7382	−38:50:26.056	0044	DES	22.45 ± 0.06	r	27.44	27.64
2016-01-11.095236	01:29:57.7285	−38:50:25.964	0055	DES	23.08 ± 0.09	g	27.44	27.64
2016-08-09.585360	01:52:24.3560	−39:29:37.720	0090	PS1	23.53 ± 0.18*	w	26.68	26.21
2016-08-09.595282	01:52:24.3230	−39:29:39.130	0129	PS1	23.45 ± 0.22*	w	26.68	26.21
2016-10-01.297615	⋯	⋯	⋯	DES	22.10 ± 0.07	i	26.49	25.80
2016-10-01.298994	01:46:40.1803	−41:11:31.779	0043	DES	22.16 ± 0.05	r	26.49	25.80
2016-10-01.300363	01:46:40.1727	−41:11:32.002	0069	DES	22.63 ± 0.05	g	26.49	25.80
2016-10-03.314076	⋯	⋯	⋯	DES	22.43 ± 0.10	i	26.48	25.80
2016-10-03.315455	01:46:20.1770	−41:14:24.527	0038	DES	22.36 ± 0.05	r	26.48	25.80
2016-10-03.316819	01:46:20.1667	−41:14:24.592	0050	DES	22.81 ± 0.07	g	26.48	25.80
2017-08-13.571449	01:58:27.6990	−41:56:20.770	0114	PS1	22.52 ± 0.18*	i	25.34	24.86
2017-08-14.585122	01:58:24.7020	−41:58:43.760	0054	PS1	21.70 ± 0.09*	i	25.34	24.85
2017-09-08.554246	01:56:09.8760*	−42:55:16.230*	0102	PS1	22.00 ± 0.20*	i	25.24	24.62
2017-09-30.421878	01:52:49.6350	−43:36:32.720	0112	PS1	21.96 ± 0.15*	i	25.16	24.51
2017-10-15.272083	01:50:05.7460	−43:57:07.295	0056	DES	22.29 ± 0.05	g	25.11	24.49
2017-10-15.337769	01:50:04.9880	−43:57:11.797	0063	DES	21.42 ± 0.08	z	25.11	24.49
2017-10-15.339148	01:50:04.9697	−43:57:11.890	0031	DES	21.69 ± 0.04	r	25.11	24.49
2017-10-29.358251	01:47:22.7320	−44:09:45.240	0145	PS1	21.54 ± 0.20*	i	25.06	24.52
2017-11-06.330497	01:45:51.8610	−44:13:42.320	0131	PS1	22.08 ± 0.21*	i	25.03	24.54
2017-12-11.250030	01:40:25.1950	−44:04:24.940	0127	PS1	22.01 ± 0.15*	w	24.90	24.76
2017-12-15.180418	01:39:59.8371	−44:00:56.398	0029	DES	22.05 ± 0.04	r	24.89	24.79
2017-12-15.181814	01:39:59.8244	−44:00:56.379	0049	DES	22.58 ± 0.06	g	24.89	24.79
2017-12-25.149228	⋯	⋯	⋯	DES	21.87 ± 0.29	Y	24.85	24.87
2018-09-10.355826	⋯	⋯	⋯	DES	22.06 ± 0.03	g	23.90	23.31
2018-10-04	⋯	⋯	⋯	TESS	20.29 ± 0.15	T	23.82	23.21
2018-10-21.243368	01:55:54.0391	−46:47:14.935	0075	DES	21.60 ± 0.17	Y	23.75	23.20
2018-10-27.178589	01:54:38.2037	−46:52:42.842	0041	DES	21.59 ± 0.07	z	23.73	23.21
2018-11-08.235135	01:52:05.9144	−46:59:29.069	0022	DES	21.71 ± 0.03	r	23.69	23.25
2018-11-08.236509	01:52:05.8988	−46:59:29.102	0026	DES	21.51 ± 0.04	i	23.69	23.25
2018-11-08.237892	01:52:05.8841	−46:59:29.051	0032	DES	22.23 ± 0.04	g	23.69	23.25
2019-08-19.598297	02:14:55.2390	−47:31:52.370	0109	PS1	⋯	i	22.64	22.18
2019-08-29.569869	02:14:08.3560	−47:59:21.660	0188	PS1	⋯	w	22.61	22.10
2019-08-29.585234	02:14:08.2450	−47:59:24.350	0149	PS1	⋯	w	22.61	22.10
2019-08-29.592922	02:14:08.2250	−47:59:25.260	0144	PS1	⋯	w	22.61	22.10
2020-09-21	⋯	⋯	⋯	TESS	18.24 ± 0.15	T	21.18	20.67

Note. Data marked with an asterisk are considered unreliable and not used in the analyses. Magnitudes for the 2019 PS1 observations are aperture-dependent, so they are not tabulated here.

Download table as:  ASCIITypeset image

For the VISTA and CFHT detections, we acquire detrended images from their respective archives ⁵⁴ and use MAG_AUTO measurements from SExtractor (Bertin & Arnouts 1996). Each exposure is placed on the DES magnitude system by choosing a zero-point to match the magnitudes found in the DES coadd catalogs in the corresponding filter for matching objects in the field. Bandpass differences between the VISTA z and CFHT gri filters and their DES counterparts lead to color corrections that are well below the measurement errors on these points and are ignored.

The PS1 photometry given in Table 2 is derived by fitting Gaussians to the comet images and stars of known magnitude (the w-band images use r-band magnitudes of the standards) and scaling the Gaussian fits. This photometry is less reliable and has a lower S/N than the DES data, so we will not make use of it in characterizing the nucleus. The PS1 measurements taken in 2019 are, however, valuable for characterizing the development of a coma between the end of DES in 2018 and the 2021 recoveries. We extract aperture photometry for these images around the predicted positions of BB to form the curves of growth shown in Section 4. The magnitude zero-points of the i and w images are determined by comparison of 6'' diameter aperture photometry of bright stars to their i- and r-band magnitudes in the DES catalogs. The color terms between the PS1 and DES bands are again well below the measurement errors (Equation (B6) of the Dark Energy Survey Collaboration 2021).

TESS observations consist of "sectors," 24° × 96° regions of the sky observed nearly continuously for approximately 4 weeks (Ricker et al. 2015). In its survey of the southern sky, TESS observed BB in three sectors, one in late 2018 and two in late 2020. For all three sectors, we identify and cut out an approximately 15 × 05 region of the TESS full-frame images along the path of the comet with the tesscut tool (Brasseur et al. 2019). We then apply a difference imaging scheme aimed toward removing background stars by, for each frame and pixel, subtracting the mean flux observed in that pixel in all cadences observed between 5 and 10 hr from the time of the frame of interest.

At each frame, we then measure the flux of the target in an aperture of 5 × 5 of TESS's 21'' pixels. We apply the same method to nearby stars on the detector with low (<1%) levels of photometric variability and well-characterized TESS magnitudes to transform our measured fluxes to magnitudes. The scatter of the residuals for stars on this scale is 0.15 mag, likely due to crowding of faint stars and intrapixel sensitivity variations on the TESS detector (Vorobiev et al. 2019). These issues should be less dramatic for BB due to its motion across the detector; nonetheless, we apply this 0.15 mag uncertainty conservatively on the individual magnitudes. The brightening of 2.01 ± 0.04 mag between the two TESS epochs is more reliably determined than the magnitude at either epoch. For an object of solar color, we should find r − T = r_⊙ − T_⊙ = 4.61 − 4.26 = 0.35 mag (Stassun et al. 2018; Willmer 2018).

We fit all of the valid photometry from VISTA, CFHT, and DES to a model in which there is a fixed absolute magnitude H_b in each band b and an achromatic light-curve fluctuation ${\rm{\Delta }}H=A\sin \phi .$ These data are plotted in Figure 2. The illumination phase is between 14 and 25 for all observations here, so we ignore the phase terms in converting observed magnitudes to H. We have insufficient data to determine the light-curve phases ϕ_i for each exposure, so we consider each observation to have a random, independent ϕ_i ∈ [0, 2π]. The posterior probability of the light-curve amplitude and the "true" H_b , given observations of H_i in band b_i with uncertainty σ_i for each observation i, is

$$\begin{eqnarray}&&\begin{array}{l}p\left(A,\{{H}_{b}\}| \{{H}_{i},{\sigma }_{i}\}\right)\\ \quad \propto \displaystyle \prod _{i\in \mathrm{obs}}\int {\rm{d}}\phi \,\exp \left[-{\left({H}_{i}-{H}_{{b}_{i}}-A\sin \phi \right)}^{2}/2{\sigma }_{i}^{2}\right].\end{array}\end{eqnarray} \tag{ 8 }$$

Figure 3 (left) plots the posterior probability for the light-curve semiamplitude A. All of the bands are consistent and combine to yield A = 0.20 ± 0.03 mag. This is a strong detection of variability in excess of the measurement errors. From Figure 2, it is clear that most of this variability is in short-term variation, not a long-term trend, consistent with a nuclear body with 10%–20% departures from sphericity. Ridden-Harper et al. (2021) reported a nondetection of variation in the TESS photometric time series, though no upper limit is reported. The TESS photometry has the majority of its flux coming from the coma, which will suppress the amplitude of any nuclear light curve, and hence does not exclude the possibility of ±0.2 mag of nuclear variation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Unresolved outbursts could also change the measured "point-source" fluxes. We can conclude, though, that the ±20% fluctuations are upper bounds for both fluctuations from outbursts and for nuclear variability. Thus, BB's nucleus is close either to sphericity or to pole-on rotation.

The fitting process yields estimates of the mean absolute magnitude of the comet of H = {8.51 ± 0.04, 7.96 ± 0.03, 7.91 ± 0.05, 7.68 ± 0.06, 7.79 ± 0.14} in the grizY bands, respectively. Figure 3 (right) plots the implied reflectivity in each band, normalized to unity in the r band, showing a color only slightly redder than neutral and perhaps even slightly blue in r − i. Jewitt (2015) presented colors for various outer solar system bodies quantified by a fit to a model of linear reflectance versus wavelength with a slope of S% per 100 nm when normalized to unit reflectivity at 550 nm (V band). Fitting the above H values (omitting Y) to such a model yields S = 4.9 ± 1.4, with χ²/dof = 4.4/2. The solid line in the right panel of Figure 3 shows the S = 5 law.

An alternative method of deriving colors for BB is to find pairs of exposures taken in different bands within 5–10 minutes of each other, so that light-curve variations are unimportant. This results in color estimates of g − r = 0.49 ± 0.01, r − i = 0.22 ± 0.02, i − z = 0.32 ± 0.09, and g − z = 0.87 ± 0.04. These agree with the mean-H method, except that the pair-based r − i color is significantly redder. Fitting the 24 closely timed observations from 10 distinct nights to the linear-reflectance model yields S = 6.1 ± 1.1, with χ²/dof = 15.1/13, in good agreement with the value derived from the H values.

Jewitt (2015) reported that potential relatives of BB, namely LPC nuclei and Damocloids, have typical S values of 10 and 15, respectively. These Oort bodies are significantly bluer than the TNO populations. Comet BB shares this deviation from the TNO colors, in fact appearing a bit more neutral than the few other well-measured Oort cloud migrants.

If the flux measurements in these ≤2018 exposures were significantly contaminated with comae rather than being predominately nuclear, we might expect the measured H to increase as the comet approaches the Sun. In Figure 4, we present the H_r averaged over all photometric observations from a given season. Exposures from band b are shifted to the r band using the H_r − H_b from the posterior maximization above. An rms error of $0.20/\sqrt{2}$ mag is added in quadrature to each measurement error to include noise from random sampling of a sinusoidal light curve. While the year-to-year means of the DES observations are formally inconsistent with a constant magnitude, the potential year-scale variation is small (≈0.1 mag) and shows no long-term trend. Indeed, the 2010 VISTA observation is consistent with a constant magnitude as well, so there is no evidence for a brightening of BB's absolute magnitude as it moves from r_h = 34.1 to 23.7 au for apertures of ≈1'' size.

Zoom In Zoom Out Reset image size

Under the assumption that the H_r = 7.96 derived above is entirely from a spherical nucleus with geometric albedo in the r band of p_r and density ρ, the diameter, mass, and escape velocity of BB are

$$\begin{eqnarray}&&D={\left(\displaystyle \frac{{p}_{r}}{0.04}\right)}^{-1/2}\,155\,\mathrm{km},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&m={\left(\displaystyle \frac{{p}_{r}}{0.04}\right)}^{-3/2}\left(\displaystyle \frac{\rho }{1\,{\rm{g}}\,{\mathrm{cm}}^{-3}}\right)2.0\times {10}^{18}\,\mathrm{kg},\end{eqnarray} \tag{ 10 }$$

and

$$\begin{eqnarray}&&{v}_{\mathrm{esc}}={\left(\displaystyle \frac{{p}_{r}}{0.04}\right)}^{-1/2}{\left(\displaystyle \frac{\rho }{1\,{\rm{g}}\,{\mathrm{cm}}^{-3}}\right)}^{1/2}\,58\,{\rm{m}}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 11 }$$

At the nominal assumed albedo, this makes BB a factor of 2.5 larger in diameter than C/1995 O1 (Hale-Bopp) (Fernández 2002), another LPC that is the largest of any comet in the past century (Lamy et al. 2004) and had H_r ≈ 9.7 at incoming r_h = 6.4 au (Szabó et al. 2012). Interestingly, this places BB at a comparable size to the 160–180 km diameter estimated for the active Centaur (2060) Chiron using occultation data (Bus et al. 1996; Sickafoose et al. 2020). It will be of interest to compare the behavior of these two objects of comparable size but very different dynamical pathways to active states at r_h of 10–20 au.

Of course, BB could be smaller than this if its albedo is above the canonical p_r = 0.04, e.g., if it has extensive clean ice patches on its surface.

For a single species with molecular mass m_mol sublimating into vacuum from a surface, the mass-loss rate per unit area A is

$$\begin{eqnarray}&&\dot{m}\equiv \displaystyle \frac{\dot{M}}{A}={P}_{\mathrm{sat}}\sqrt{\displaystyle \frac{{m}_{\mathrm{mol}}}{2\pi {kT}}}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\propto {e}^{-{\rm{\Delta }}H/{RT}}\sqrt{\displaystyle \frac{{m}_{\mathrm{mol}}}{2\pi {kT}}}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\Rightarrow \qquad \mathrm{log}\dot{M}\sqrt{T}=\mathrm{const}-\displaystyle \frac{{\rm{\Delta }}H}{{RT}},\end{eqnarray} \tag{ 16 }$$

where P_sat is the saturation vapor pressure, and we use the Clausius–Clapeyron formula to express its dependence on temperature T in Equation (15). Here ΔH is the enthalpy of sublimation (which we assume varies little with T), and R is the ideal gas constant.

Under radiative equilibrium with negligible heat conduction with the cometary interior and negligible heat loss to sublimation, a section of the surface attains a temperature T with

$$\begin{eqnarray}&&\epsilon \sigma {T}^{4}=\displaystyle \frac{{L}_{\odot }}{4\pi {r}_{h}^{2}}(1-p)\langle \cos \theta \rangle \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\Rightarrow \ T={\left[\displaystyle \frac{(1-p)\langle \cos \theta \rangle }{\epsilon }\right]}^{1/4}{\left(\displaystyle \frac{{L}_{\odot }}{4\pi \sigma {\left(1\,\mathrm{au}\right)}^{2}}\right)}^{1/4}{\left(\displaystyle \frac{{r}_{h}}{1\,\mathrm{au}}\right)}^{-1/2}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&=\eta \times \left(397\,{\rm{K}}\right){\left(\displaystyle \frac{{r}_{h}}{1\,\mathrm{au}}\right)}^{-1/2}.\end{eqnarray} \tag{ 19 }$$

In these equations, σ is the Stefan–Boltzmann constant, p is the Bond albedo of the surface (nominally 0.04), and is the infrared emissivity (nominally 0.9). With θ as the angle between illumination and the normal, the average $\langle \cos \theta \rangle$ over the thermal timescale for the warmest part of the comet (which will dominate the sublimation rate at low T) will be between 1 (at the subsolar point for a short thermal time constant or a pole-on rotator) and 1/π for the equator of an orthogonal rotator with a long time constant. We bundle all of these physical/geometric constants in the first term of Equation (18) into a factor η, which is nominally close to unity but could be as low as ≈0.75. ⁵⁵

The third part of the simple model is to relate the coma brightness A f ρ to the sublimation rate. If the scattering is dominated by solid particles with albedo p_d , radius a_d , near-spherical geometric cross section, density ρ_d , production rate ${\dot{M}}_{d}$, and velocity v_d , then

$$\begin{eqnarray}&&Af\rho =\displaystyle \frac{3{p}_{d}{\dot{M}}_{d}}{8\pi {a}_{d}{\rho }_{d}{v}_{d}}.\end{eqnarray} \tag{ 20 }$$

If we assume that p_d , a_d , ρ_d , and the dust-to-gas ratio $\chi ={\dot{M}}_{d}/\dot{M}$ are independent of heliocentric distance over the 20–30 au range, we obtain a scaling

$$\begin{eqnarray}&&Af\rho \propto \dot{M}/{v}_{d}.\end{eqnarray} \tag{ 21 }$$

Note that this proportionality does not require the geometric scattering limit to hold, only that the scattering per unit mass of dust is time-invariant, i.e., no grain destruction or breakdown. Combining this with Equation (16) yields

$$\begin{eqnarray}&&\mathrm{log}\left(Af\rho {v}_{d}{T}^{1/2}\right)=\mathrm{const}-\displaystyle \frac{{\rm{\Delta }}H}{{RT}}.\end{eqnarray} \tag{ 22 }$$

One working assumption for v_d is that it will scale with the thermal velocity, i.e., ∝ T^1/2. A stronger dependence would be expected if the dust velocity is driven by radiation pressure: ${v}_{d}\propto {r}_{h}^{-2}\propto {T}^{4}.$

The left panel of Figure 7 plots A f ρ v_d T^1/2 on the (logarithmic) y-axis versus 1/T on the x-axis, assuming radiative equilibrium for T and v_d ∝ v_th. Per Equation (22), the data should follow a line with slope −ΔH/R on this plot for sublimation of a single species. The data are seen to be consistent with this model, with a value of ΔH that depends on whether we take the fast- or slow-rotator limit for η in the radiative equilibrium formula.

Zoom In Zoom Out Reset image size

The right panel of Figure 7 plots the result of a quantitative fit of Equation (22) to the A f ρ measurements as a function of ΔH, marginalizing over the unknown constant. In the nominal (red) case, the inferred enthalpy of sublimation is fully consistent with the lab-measured values for NH₃ and/or CO₂ and clearly inconsistent with the more volatile species N₂, CH₄, and CO. As is well known, H₂O is ruled out as a volatile at these distances as well. The values of ΔH are taken from Luna et al. (2014) and Feistel & Wagner (2007). The dashed green curve changes the scaling of v_d versus T from thermal to radiation-pressure laws and does not change the conclusion. Moving to the fast-rotator value η = 0.75 yields proportionately lower ΔH values, still better associated with CO₂ and NH₃.

The coma growth rate is thus strongly suggestive of activity powered by CO₂ and/or NH₃ sublimation at 20 au < r_h < 25 au, with slow and/or pole-on rotation somewhat favored to yield higher peak surface temperatures. Across the range of potential surface temperatures, the mass-loss rate per square meter from pure CO₂ is >100× larger than that from a pure NH₃ ice surface, so CO₂ would have to be far less abundant than NH₃ in order for the latter to be driving the activity rate in pure-ice forms.

At r_h > 25 au, the coma is too weak to be well measured in the data. It is possible that more volatile species (N₂, CH₄, or CO) could have dominated sublimation at these times. These latter species, however, have orders-of-magnitude higher vapor pressure (and specific sublimation rates) than CO₂ and NH₃ at T = 78 K, the η = 1 subsolar temperature for r_h = 25 au. Any surface abundance in their pure-ice forms on the surface of BB must be very low relative to the CO₂/NH₃ that appears to dominate sublimation at r_h ≤ 25 au. A scenario that economically ties the dynamical and thermodynamic results above is that these more volatile species were heavily sublimated from BB during its previous perihelion passage to ≈18 au, leaving behind a crust that is largely depleted in these most volatile species. Or perhaps this depletion is a remnant of the thermal environment at the original location of formation of BB.

The dynamics of dust production on BB are made more interesting by the fact that the escape velocity of ≈60 m s⁻¹ is above or comparable to the dust velocities estimated for other comets, e.g., <50 m s⁻¹ for Comet Boattini (Hui et al. 2019) or ≈4 m s⁻¹ for C/2017 K2 (Jewitt et al. 2019). The low sublimation rates and pressures for BB at r_h > 22 au may limit the size of particles that can be raised and attain escape by the gas, and/or the coma may be dominated by small grains that are driven to escape by radiation pressure. The coma is observed to be asymmetric in its 2021 observations, the 2018 DES data (Figure 5), and potentially the TESS observations (Farnham 2021), implicating radiation pressure or tidal escape mechanisms for the dust. High-quality imaging of the coma as soon as possible would be of great interest in constraining the dust size and dynamics.

If a fraction f_active of the surface of BB is sublimating CO₂ at the rates described in Equation (14), then the gas mass-loss rate from BB rises from 400f_active to 7 × 10⁴ f_active kg s⁻¹ when closing from 26 to 20 au. This assumes η ≈ 1 (slow rotator) and the vapor pressure cited by Fray & Schmitt (2009). A nonporous, pure-CO₂ ice patch at the subsolar point would be eroded by 0.1 mm yr⁻¹ at 26 au, increasing to 2 cm yr⁻¹ at its current r_h ≈ 20 au.

Adopting the simplistic model of Equation (20) and a uniform streaming velocity v_d = f_v v_th, where v_th is the rms 1D thermal velocity of a CO₂ molecule, we can solve for the dust particle radius as

$$\begin{eqnarray}&&{a}_{d}=\sqrt{\displaystyle \frac{9}{8\pi }}\displaystyle \frac{{p}_{d}\chi {f}_{\mathrm{active}}}{{f}_{v}}\displaystyle \frac{{P}_{\mathrm{sat}}{a}_{\mathrm{BB}}^{2}{m}_{\mathrm{mol}}}{{\rho }_{d}[Af\rho ]{kT}}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\approx \left(\displaystyle \frac{p\chi }{0.04}\right)\left(\displaystyle \frac{0.1}{{f}_{v}}\right)\left(\displaystyle \frac{{f}_{\mathrm{active}}}{0.1}\right)\times 10\,\mu {\rm{m}}.\end{eqnarray} \tag{ 24 }$$

We take a_BB as the nominal radius of the comet, m_mol as the mass of a CO₂ molecule, albedos of p = 0.04 for both the comet and its debris, and the observed values of A f ρ to obtain Equation (24). The nominal values of χ, f_active, and f_v are ill-supported order-of-magnitude estimates. The implied dust particle radius may even be overestimated, as the assumption of f_active = 0.1 of the surface being CO₂ ice near the subsolar temperature may be an overestimate. It nonetheless suggests that the coma has small particles, indeed small enough that the assumption of a geometric scattering cross section needs to be relaxed. Small dust particles might be expected because of the anemic gas flow density and the importance of radiation pressure in overcoming BB's gravity. Clearly, a more detailed model of the dust dynamics would be of interest.

Opportunities to study incoming comets at r_h > 20 au have been rare enough that there is no definition of "typical" behavior. It is already clear that the behavior is diverse. It is good to keep in mind that selection biases will favor the discovery at large distances of comets that are unusually active or, like BB, unusually large.

Comet C/2017 K2 (PanSTARRS), hereafter K2, has a similar aphelion and inclination to BB, but Królikowska & Dybczyński (2018) reported it to have a 97% chance of having passed within r_h < 10 au on its previous passage. It was discovered at r_h = 16 au, but pre-discovery images at r_h = 24 au also display a coma (Hui et al. 2018).

Comet K2 is much smaller than BB, with an estimated nuclear radius of <9 km (Jewitt et al. 2017). These authors inferred K2's coma to be comprised of millimeter-scale particles moving at speeds of ≈4 m s⁻¹—which would probably not escape BB—and estimate that K2 has been expelling such particles at a relatively steady and isotropic rate since r_h ≈ 30 au. This is in stark contrast to BB's exponentially growing scattering cross section in the 20–30 au range.

Comet C/2010 U3 (Boattini) was discovered inbound at r_H = 18.4 au and found on earlier images at r_h = 25.8 and 24.6 au with a visible coma (Hui et al. 2019). Like K2, Boattini is thought to have had its previous perihelion at <10 au, yet its coma behavior is quite distinct from K2's, showing intermittent outbursts and a tail inferred to be composed of much smaller particles than K2's.

Comets K2, Boattini, and BB display very diverse activity patterns at r_h > 15 au. Comet K2 has a steady, large-particle coma, Boattini exhibits outbursts and a tail, and BB undergoes exponential growth in cross section over this period that is consistent with simple sublimation thermodynamics of carbon dioxide and/or ammonia. There is also diversity in the apparent role of CO in distant activity. CO emission has been detected from K2 at 7 au (incoming; Yang et al. 2021), and CO was detected in the coma of the large LPC Hale-Bopp and inferred to drive its activity at distances as large as 26 au (Gunnarsson et al. 2003; Szabó et al. 2008). Womack et al. (2017) reviewed the data on CO emission at distances of 4–11 au for Hale-Bopp, Chiron, and 29P/Schwassmann-Wachmann 1 (another active Centaur) and concluded that CO drives activity in all three objects, albeit with very diverse temporal behavior. This makes the current agreement of BB's coma growth with CO₂ ice sublimation, rather than pure CO ice sublimation, of interest; is there a skin-deep lack of CO, such that large quantities of buried CO could be liberated as fresh material is uncovered nearer perihelion? Searches for CO sublimation in the coming years are clearly desirable.

No generalized pattern is apparent yet for the behavior of distant LPCs, aside from the existence of activity in some form out to ≈30 au. Even this ubiquity can be explained by selection effects on a more diverse population.