Dynamical Implantation of Blue Binaries in the Cold Classical Kuiper Belt

The Kuiper Belt holds important clues about the solar system formation. Neptune's early migration (Fernández & Ip 1981; Malhotra 1993), in particular, is thought to be responsible for the complex orbital distribution of Kuiper Belt objects (KBOs). Observations and dynamical modeling of KBOs (see Morbidelli & Nesvorný 2020 and Gladman & Volk 2021 for recent reviews) help us to understand: (i) the mass and radial structure of the original outer disk of planetesimals, and (ii) the nature and timescale of Neptune's migration (e.g., Hahn & Malhotra 2005; Levison et al. 2008; Nesvorný & Vokrouhlický 2016). For example, it is inferred that Neptune relatively slowly migrated (one e-fold ≳10 Myr; Nesvorný 2015a) on a modestly eccentric orbit (e ∼ 0.03–0.1; Dawson & Murray-Clay 2012; Nesvorný 2021). Neptune's orbit could have been excited by a dynamical instability (Tsiganis et al. 2005; Nesvorný & Morbidelli 2012).

KBO colors and binarity provide crucial constraints. The color distribution of small KBOs is bimodal with less-red (LR) and very-red (VR) objects observed, in different proportions, in all dynamical categories (e.g., Wong & Brown 2017). Complex organics and surface weathering are thought to be responsible for different colors. The VR bodies are believed to have formed beyond an iceline of some volatile such as methanol or ammonia (Brown et al. 2011). The exact location of the LR-to-VR transition, r*, is unknown. The dynamically cold KBOs at r = 42–45 au are predominantly VR, which indicates r* < 42 au. If r* < 30 au, however, the dynamically hot populations such as Plutinos, hot classicals (HCs) and scattered-disk objects would presumably have a much larger proportion of VR objects (Nesvorný et al. 2020; Ali-Dib et al. 2021). This indicates 30 < r* < 42 au.

Observations reveal an unusually high fraction of binaries among cold classicals (CCs; ∼30%; Noll et al. 2020). The CC binaries have nearly equal-size components (R₂/R₁ > 0.5, where R₁ and R₂ are the primary and secondary radii) and wide separations (a_b/(R₁ + R₂) ∼ 10–1000, where a_b is the binary semimajor axis). They presumably formed during the formation of KBOs themselves, by the streaming instability (Youdin & Goodman 2005) and subsequent gravitational collapse (Nesvorný et al. 2010, 2019). Their survival in different KBO categories reflects how much these populations were affected by collisions and dynamical processes (Petit & Mousis 2004; Parker & Kavelaars 2010). For example, only a very small fraction of wide binaries is expected to survive the dynamical implantation from r < 30 au (most become unbound during encounters to Neptune; Parker & Kavelaars 2010; Nesvorný & Vokrouhlický 2019).

Fraser et al. (2017, 2021) reported detection of several "blue" binaries in the CC population (i.e., LR colored; spectral slope s < 17% as defined in Fraser et al. 2017; Table 1). This is surprising because if these binaries formed—like most CCs—beyond 42 au, they should be VR—like most CCs. They could not have formed in the massive disk below 30 au, which is the source of most dynamically hot KBOs, presumably because they would not, especially the wide ones, survive (Parker & Kavelaars 2010). It has been suggested (Fraser et al. 2017, 2021) that the blue binaries formed at r ∼ 38 au, from where that were pushed—by migrating Neptune's 2:1 resonance—to r > 42 au. As this implantation mechanism does not involve gravitational scattering by Neptune, binaries are retained. This interpretation would indicate the LR-to-VR transition at r* ≃ 38–42 au.

Here we reexamine the implantation of blue binaries from r < 42 au to r > 42 au. The best dynamical models from Nesvorný et al. (2020) are used to determine the implantation probability to an orbit with the semimajor axis 42 < a < 47 au, inclination i < 5°, and stable perihelion distance q > 36 au. The dynamical survival of binaries during planetary encounters is evaluated as a function of the mutual separation of binary components. We find that the widest blue binaries—2001 QW322 and 2003 UN284—would not survive if they formed below 30 au (≲1% survival probability). Their survival becomes progressively more likely with the increasing starting radius: ∼10%–25% for 30 < r < 36 au, ∼40%–50% for 36 < r < 38 au, and ∼100% for r > 38 au. We discuss the implication of these, and other results obtained here, for blue binaries, the LR-to-VR transition in the original disk, and the interpretation of KBO colors in general.

We make use of two Kuiper Belt simulations published in Nesvorný et al. (2020) they are referred to as 10/30 and 30/100 here (Table 2). See that work for the description of the integration method, planet migration, initial orbital distribution of disk planetesimals, and comparison of the results with the orbital structure of the Kuiper Belt. A shared property of the selected runs is that Neptune migrates outward by scattering planetesimals. Planetesimals are initially distributed in a disk extending from just beyond the initial orbit of Neptune at 22 au to >50 au. The radial distribution of planetesimals is assumed to have a drop at ≃30 au, with the massive inner disk and low-mass outer extension, or a smooth profile with the surface density exponentially decreasing with distance. Either of these initial distributions can match existing observational constraints (Nesvorný et al. 2020). The simulations were performed with the symplectic N-body integrator known as Swift (Levison & Duncan 1994).

All encounters of planetesimals with planets were recorded. This was done by monitoring the distance of each planetesimal from Jupiter, Saturn, Uranus, and Neptune, and recording every instance when the distance dropped below 0.5 R_Hill,j , where R_Hill,j are the Hill radii of planets (j = 5 to 8 from Jupiter to Neptune, the inner planets were not included). We verified that the results do not change when more distant encounters are accounted for. The disk planetesimals that ended up on orbits with with 42 < a < 47 au, q > 36 au, and i < 5° at the end of the simulations (t = 4.5 Gyr) were selected and used for the analysis presented here. See Nesvorný & Vokrouhlický (2019) for the binary survival during implantation into the dynamically hot populations.

Each selected planetesimal was assumed to be a binary object. We considered a range of binary separations, 1 < a_b/R_b < 2000, where (Agnor & Hamilton 2006), initially circular orbits (binary orbit eccentricity e_b = 0), and a random distribution of binary inclinations (i_b). Ten clones with different binary inclinations were assigned to each selected planetesimal to increase the statistic. Each binary was evolved through all recorded planetary encounters. We used the Bulirsch–Stoer (B-S) N-body integrator that we adapted from Numerical Recipes (Press et al. 1992). The center of mass of the binary was first integrated backward from the time of the closest approach to 3 R_Hill,j . It was then replaced by the actual binary and integrated forward through the encounter until the planetocentric distance of the binary exceeded 3R_Hill,j . The final binary orbit was used as an initial orbit for the next encounter and the algorithm was repeated over all encounters.

The B-S code monitored collisions between binary components. If a collision occurred, the integration was stopped and the impact speed and angle were recorded. A fraction of binaries became unbound. For the surviving binaries, we recorded the final values of a_b, e_b and i_b, which were then used to evaluate the overall change of orbits. The statistics of surviving binaries was used to compute the binary survival probability, p_b, as a function of a_b/R_b.

Figure 1 shows the implantation probability onto an orbit with 42 < a < 47 au, i < 5° and q > 36 au from different parts of the original planetesimal disk. The implantation probability is computed as p_j = n_j /N_j , where N_j is the number of planetesimals starting with r_j−1 < r < r_j in our simulations, and n_j is the number of planetesimals starting with r_j−1 < r < r_j and ending with 42 < a < 47 au, i < 5° and q > 36 au. Probability p_j does not account for the availability of planetesimals in different zones. We used eight zones in total with r₀ = 24 au r_j = 30 + 2(j − 1) au for j = 1,...,7, and r₈ = 47 au. Note that the first zone (j = 1) corresponds to the full radial span of the inner massive disk, here assumed to be 6 au wide, and the eighth zone (j = 8) is the orbital radius of the cold classical disk; all other zones are 2 au wide.

Zoom In Zoom Out Reset image size

The implantation probability increases with the initial radial distance, from p₂ ∼ 10⁻⁴ for r ≃ 31 au to p₇ ≃ (2–3) × 10⁻³ for r ≃ 41 au. It is apparently easier to implant bodies in the cold classical belt if they start radially closer to the belt. For the eighth zone, where bodies already start in the CC region, p₈ is the probability of survival on a CC orbit. We find, consistently with the previous results (Nesvorný 2015b), that 30%–40% of planetesimals starting in the cold classical belt survive (p₈ ≃ 0.3–0.4). The implantation probabilities from the massive disk below 30 au are low: p₁ ≃ 8 × 10⁻⁵ for 10/30 and p₁ ≃ 2 × 10⁻⁵ for 30/100. The probability profile less steeply raises with the initial orbital radius if Neptune's migration is faster, and is practically flat between 30 and 42 au for the 10/30 simulation.

We now multiply the implantation probability by the number of planetesimals initially available in each zone, η_j . We consider two cases from Nesvorný et al. (2020): the (1) truncated power-law profile (the surface density Σ ∝ 1/r, truncated at 30 au), and (2) exponential profile (, where r₀ denotes the inner edge radius at 24 au, and Δr au is one e-fold, no outer truncation here). For the truncated power-law profile, the change of the surface density at r = 30 au is parametrized by the contrast parameter, c, which is simply the ratio of surface densities on either side of 30 au. Here we use Δr = 2.5 au and c = 1000. Both these cases were shown to correctly reproduce the orbital and color distribution of KBOs (Nesvorný et al. 2020).

We define n_CC,j = p_j η_j , where n_CC,j represents the actual number of planetesimals implanted onto the CC orbits from different zones, and normalize n_CC,j such that n_CC,8 = 10⁴. This means that we expect 10⁴ planetesimals to survive in the CC region ( ∼ 1/3 of the original population; Figure 1). For comparison, Kavelaars et al. (2021) estimated from the Outer Solar System Origins Survey (Bannister et al. 2018) observations that there are ≃5000–15,000 CCs with diameter D > 100 km (the range given here reflects the uncertain albedo of CCs, p_V = 0.05–0.2).

With the truncated power-law profile, the 30–42 au region does not contribute much to the population of CCs, for either migration run considered here (Figure 2). In this case, we would expect ∼100 bodies to be implanted from r = 30–42 au, only ∼1% of today's CC population. For the inner part of the planetesimal disk, r < 30 au, we find n_CC,1 ∼ 2500 for 10/30 and n_CC,1 ∼ 800 for 30/100. This would represent ∼8%–25% of the CC population and imply that nearly all blue CC binaries would need to be implanted from r < 30 au. This is difficult to reconcile with observations (Fraser et al. 2017, 2021), because the wide binaries such as 2001 QW322 and 2003 UN284 would not dynamically survive (Section 4; Nesvorný & Vokrouhlický 2019). The population of blue cold classicals would be dominated by singles in this case—the opposite to what is observed.

Zoom In Zoom Out Reset image size

The results with the exponential disk profile are somewhat more plausible. Here, the contribution of the r < 30 au and 30 < r < 42 au source regions would be similar: 15%–30% in 10/30 and ∼6%–10% in 30/100 (percentages given relative to the present CC population; Kavelaars et al. 2021). This would also produce roughly the right proportion of LR/VR CCs assuming that 30 < r* ≲ 40 au for 10/30 or 35 < r* ≲ 40 au for 30/100 (Nesvorný et al. 2020). In the next section, we discuss the dynamical survival of binaries starting in different zones.

Figure 3 shows how the survival probability changes with the starting heliocentric distance and binary separation. We find that the binary survival depends on a_b/R_b, where , and not on a_b, R₁ and R₂ individually. This is a consequence of the binary dissociation condition described in Agnor & Hamilton (2006). A binary with the total mass m_b = m₁ + m₂ can become unbound when the planetocentric Hill radius of the binary, , where q is the distance of the closest approach and m_pl is the planet mass, becomes smaller that the binary separation; that is r_{Hill, b} < a_b. This condition yields

where R_pl and ρ_pl are the planet radius and density. Here we assumed that the primary and secondary components of binaries have the same density, ρ. For exactly equal-size binaries with R₁ = R₂, R_b = (R₁ + R₂)/2^2/3.

Zoom In Zoom Out Reset image size

The encounter statistics changes with the starting heliocentric distance and this projects into the expectation for binary survival. The survival probability is also a strong function of the binary separation (Figure 3). ⁴ For the ease of presentation, we consider three types of binaries: tight with a_b/R_b ∼ 10, wide with a_b/R_b ∼ 100, and extreme with a_b/R_b ∼ 1000. For reference, one tight, five wide, and two extreme blue CC binaries are currently known (Table 1).

The tight binaries with small separations are strongly bound together and have relatively high survival probabilities, p_b > 0.9 (independently of their starting orbital radius), except for the inner disk in 10/30, for which we obtain p_b,1 ≃ 0.85. The tight binaries are affected only during extremely close encounters to planets which do not happen too often. The wide binaries with a_b/R_b ≃ 100 are less likely to survive. Here the survival probability increases from p_b,1 ≃ 0.3 for r < 30 au to ≃0.8–1 for bodies starting near the inner edge of the cold belt (r = 36–42 au). The extreme binaries are not expected to survive (e.g., p_b ≃ 0.01 for r ≲ 30 au) unless they start with r > 37 au. The implications of these results for blue binaries depends on the color transition radius r* (Figure 4).

Zoom In Zoom Out Reset image size

In the model, we assign color to each simulated object depending on whether it started at r < r* ("blue" or LR) or r > r* (VR). The color transition at r* is assumed to be a sharp boundary between LR and VR. This can be a consequence of the sublimation-driven surface depletion in some organic molecules, such as methanol or ammonia (Brown et al. 2011). In this case, the color transition at r* would correspond to the sublimation radius (see Section 5 for a discussion).

The case with r* ∼ 30 au can be ruled out because the blue CC binaries would have to start with r < 30 au (to have blue colors), and the great majority of wide/extreme binaries would not dynamically survive their implantation in the CC population. The color transition somewhere between 30 and 42 au could be more plausible. The main constrains are: (i) to implant enough blue binaries in the cold disk, and (ii) not implant too many blue singles. Indeed, the existing observations reveal only two blue CC singles (s < 17%; Fraser et al. 2017, 2021). To satisfy both conditions, a very large binary fraction in the source region, r < r*, is needed (Fraser et al. 2017). In addition, a large fraction of initial binaries starting with r < r* must survive implantation; the ones that are dissolved are a source of blue singles.

None of the options considered here seems to work to satisfy the constraints: the blue singles always win over the wide/extreme blue binaries. This is illustrated in Figure 4. Here we assumed that all planetesimals are born as binaries with a fixed separation. In three different cases, we set a_b/R_b = 10, 100, and 1000. This is not meant to represent the real situation. In reality, single planetesimals must have formed as well, and many initial binaries in the inner massive disk should have been dissociated by impacts (Nesvorný & Vokrouhlický 2019; see Section 5 below), producing singles. In addition, for any reasonable assumption about their formation, the planetesimal binaries should have formed with a range of binary separations.

For wide binaries with a_b/R_b = 100, the predicted blue-binary to blue-single number ratio, R_b/s, is smaller than one—in direct contradiction to observations (Fraser et al. 2017, 2021; R_b/s ≃ 4, Table 1)—unless r* > 42 au (Figure 4). This would imply that the blue CC binaries had to form in situ at r = 42–47 au. We discuss this possibility in Section 5. The extreme binaries with a_b/R_b = 1000 are easier to dissolve, and yield much lower values of R_b/s. The tight binaries with a_b/R_b = 10 show R_b/s > 1 for any choice of the color transition radius. This is simply because they survive and dominate the final statistics over singles (recall our assumption of the 100% blue binary fraction in the source).

To explain the observations that indicate the preference for blue wide/extreme binaries over tight ones and singles (Fraser et al. 2021), it could be imagined that all planetesimals start as tight binaries and these binaries become wider. This is not the case. In fact, only a very few tight binaries became wide/extreme our simulations—this formation channel is simply not efficient enough. A similar result was obtained in Nesvorný & Vokrouhlický (2019) for binaries implanted in the dynamically hot Kuiper Belt populations. In a recent paper, Stone & Kaib (2021) suggested that a small subset of tighter binaries can give rise to a small population of very wide binaries, but their simulations show that only ∼1% of initial binaries with 30 < a_b/R_b < 100 can evolve to have a_b/R_b > 300, with most (∼90%) surviving with a_b/R_b < 100. This would not help to explain the observed statistics of blue CC binaries.

We end at an impasse: no simple combination of binary properties, radial distribution of planetesimals, and color transition radius can match the constraints. Something else must be going on. We first discuss the influence of observational biases. An important bias that favors detection of binaries over singles is that a binary is composed of two objects of some size and will appear brighter than a single object of the same size (e.g., Noll et al. 2020). For example, an exactly equal-mass binary is ≃0.75 mag brighter than a single object. A magnitude-limited survey is therefore expected to yield a higher-than-intrinsic binary fraction.

A similar bias occurs when bodies are monitored for color measurements where the observational target must be bright enough to measure accurate colors. For the size distribution of large CCs (Fraser et al. 2014), this could contribute by a multiplication factor of several in terms of R_b/s, and is not large enough to resolve the issue highlighted in Figure 4. In addition, the magnitude bias does not disfavor the detection of tight binaries (and unresolved tight binaries). We therefore find that this bias cannot be entirely responsible for the problem in question. ⁵

In this work, we used the dynamical models from Nesvorný et al. (2020). While these models were previously shown to reproduce the orbital structure of the Kuiper Belt, dynamical modeling is a source of additional uncertainties. For example, Volk & Malhotra (2019) showed that the mean orbital inclination of Plutionos and HCs can be obtained in a model where Neptune migrates on a very-low-eccentricity orbit. This migration mode could potentially have different implications for the implantation of blue binaries (although it is not clear, in detail, how this would work). To start with, however, it would be important to demonstrate that the very-low-eccentricity migration of Neptune can match other constraints as well (e.g., the inclination distribution of KBOs; see discussion in Nesvorný 2021).

Levison & Morbidelli (2003) proposed that CCs were pushed out from ∼35 au to >42 au by the 2:1 resonance with migrating Neptune. To maintain low eccentricities in this migration model, the resonant population was required to be massive (the total mass of several Earth masses). As bodies in the 2:1 resonance do not experience encounters with Neptune, this would favor binary survival. Moreover, if r* > 35 au, the implanted binaries would be blue, potentially resolving some of the tensions discussed above. On the down side, with several Earth masses in the 2:1 resonance, the collisional grinding in the resonant population would presumably be intense, casting doubt on the collisional survival of wide binaries. It is also not clear whether a very massive 2:1 population is plausible (modern dynamical models do not place much mass beyond 30 au; e.g., Levison et al. 2008; Volk & Malhotra 2019; Nesvorný et al. 2020).

Additional uncertainty is related to the initial inclination distribution of planetesimals. In Nesvorný et al. (2020), we adopted the Rayleigh distribution to set up the initial inclinations. For r < 30 au, the mean inclination was set to 〈i〉 = 0.06 (35). For r > 30 au, the mean inclination gradually decreased with the orbital radius such that 〈i〉 ≃ 17 for r = 45 au, a value comparable with the free inclination of cold classicals (Fraser et al. 2021). It is not known what the real inclination distribution was or how it varied with the orbital radius. It is possible, for example, that the initial inclination distribution was wider (e.g., because bodies starting in the massive planetesimal disk have their inclinations excited by large objects that formed in the disk; Nesvorný & Vokrouhlický 2016) and more strongly varied with the orbital radius. If so, we imagine that the implantation probability from r ≲ 35 au to 42 < a < 47 au, i < 5° and q > 36 au would be reduced—relative to our nominal model—because fewer bodies could reach orbits with very low inclinations. This would then presumably shift the balance toward having more CCs implanted from larger orbital radii, where the binary survival is better, increase R_b/s, and potentially improve the match to observations. This is left for future work.

It is also conceivable that the wide/extreme blue binaries were not implanted in the CC population but instead formed in situ at r = 42–47 au. ⁶ If that would be the case, however, we would need to explain why their LR colors differ from the predominantly VR colors of CCs. The time of formation could potentially be responsible for the color difference. KBOs presumably formed by the streaming instability (Youdin & Goodman 2005), or a related process, in a protoplanetary gas disk. During the early stages, when the protoplanetary disk was massive and the solar luminosity was high (Baraffe et al. 1998), the surface layers of a flared disk at 45 au were heated by solar radiation and the radiation was efficiently diffused to the disk midplane. The midplane temperature at r = 45 au is therefore expected to be relatively high (Bitsch et al. 2015). As the disk becomes less massive over time and the solar luminosity decreases, the temperature is expected to drop. For example, in some disk models, the midplane temperature at 45 au changes from above 30 K during the early stages to below 20 K during the late stages (e.g., Bitsch et al. 2015). This would cause carbon monoxide (CO), which is the second most important gas molecule in astronomical environments, to freeze.

The CO iceline near 45 au would have two important implications for the formation of CC planetesimals. First, the surface density of solids would be boosted near 45 au, including the contribution from the cold finger effect (Drazkowska & Alibert 2017). This could potentially produce, if grains stick and grow to become large enough (Birnstiel et al. 2016), favorable conditions for the streaming instability (Carrera et al. 2015; Yang et al. 2017; Li & Youdin 2021). Second, for temperatures below 30 K, CO can be destroyed on grains (chemistry mediated by cosmic ray ionization), producing CO₂, CH₃OH (methanol), and other hydrocarbons (e.g., Bosman et al. 2018), all of which have relatively high sublimation temperatures and can remain on a surface. Methanol, in particular, has been the only ice unambiguously detected on the surface of Arrokoth (Grundy et al. 2020)—providing evidence for hydrogenation of CO, possibly already on the surface of Arrokoth, but before the gas disk dispersal. ⁷ Methanol is known to retain high albedo and redder colors after irradiation (Brunetto et al. 2006), possibly the main reason behind the higher albedos and redder colors of CCs (Tegler & Romanishin 2000; Brucker et al. 2009).

Blue CC binaries could have formed in situ at 45 au during the earlier protoplanetary disk stages when the ambient disk temperature at 45 au was higher, CO remained in the gas phase and/or the production of CH₃OH was inefficient for T > 30 K (Bosman et al. 2018). This would lead to different chemistry and different composition of planetesimals formed under these conditions. The early-formed planetesimals would subsequently be exposed to physical conditions in a late-stage protoplanetary disk (e.g., very low temperatures), and could accrete a thin layer of methanol-rich materials (see above). This layer would have to be later disturbed (e.g., by small impacts) to expose the interior materials. Our hypothesis could potentially explain the color affinity of blue CC binaries to KBOs in the dynamically hot populations, which formed below 30 au and also experienced higher temperatures during their formation.

The model can incorporate, after the gas disk dispersal, the implantation of bodies from r < 42 au to the CC population: the majority of implanted bodies would be tight binaries and singles (Section 4). A large fraction of the implanted population—those bodies that formed earlier and/or closer to the Sun—would have LR colors. The two blue singles reported in Fraser et al. (2021) could have formed in situ or been implanted. They represent ∼2% of the total CC sample with known colors (Fraser et al. 2021). This fraction can be used to place an upper bound on the implantation process. For example, we expect ∼10% of the CC population to have LR colors in the 30/100 case (for any r* > 30 au, in a size-limited sample; Figure 2). This seems excessive but note that the LR/VR ratio is influenced by the magnitude (see above) and albedo biases (LR objects generally have lower albedos than VR objects), both of which disfavor the detection of blue singles in the CC population (relatively to red binaries). It is therefore possible that the LR/VR ratio in the CC population is higher than the magnitude-limited color surveys currently indicate.

Fraser et al. (2017, 2021) observationally characterized colors and binarity in the cold classical population of the Kuiper Belt. They reported detection of eight CC binaries with LR colors (spectral slope s < 17%), most of which have very wide separations (7 of 8 have a_b/R_b ≳ 100), and two blue singles (possibly unresolved binaries). The LR (or "blue") binaries were suggested to have been implanted in the CC population from r < 42 au. Here we studied the implantation process and found:

• 1.  
  The implantation probability increases with the increasing orbital radius of the source. For r < 30 au, it is ≃10⁻⁴ for fast Neptune's migration (run 10/30) and ≃2.5 × 10⁻⁵ for slow migration (30/100). For r ≳ 38 au, the implantation probability raises to ≳10⁻³.
• 2.  
  Our dynamical model indicates that ∼10% of cold classicals should have LR colors as these bodies were implanted onto CC orbits (42 < a < 47 au, i < 5° and q > 36 au) from below the color transition in the original planetesimal disk (this assumes r* > 30 au; Nesvorný et al. 2020).
• 3.  
  The blue binaries can be implanted from a wide range of initial heliocentric distances, but the survival of the widest binaries implies formation at r > 30 au.
• 4.  
  For any reasonable choice of parameters (Neptune's migration history, initial disk profile, etc.), the implantation model predicts a predominance of blue singles, rather than blue binaries, which contradicts existing observations.

We propose that wide blue binaries formed in situ at r = 42–47 au and their color reflects early formation and higher temperature of the young protoplanetary disk. The predominantly VR colors of cold classicals are linked to the production of methanol and other hydrocarbons during the late stages of the disk, when the temperature at 45 au dropped to ≃20 K and carbon monoxide was hydrogenated. Methanol, which is known to retain high albedo and redder colors after irradiation, was detected on the surface of Arrokoth.

The work of D.N. was supported by the NASA Emerging Worlds program. The work of D.V. was supported by the Czech Science Foundation (grant 21-11058S).