Col-OSSOS: Color and Inclination are Correlated Throughout the Kuiper Belt

Both physical and dynamical properties must be considered to constrain the origins of the dynamically excited distant Solar System populations. We present high-precision (g-r) colors for 25 small (Hr>5) dynamically excited Trans-Neptunian Objects (TNOs) and centaurs acquired as part of the Colours of the Outer Solar System Origins Survey (Col-OSSOS). We combine our dataset with previously published measurements and consider a set of 229 colors of outer Solar System objects on dynamically excited orbits. The overall color distribution is bimodal and can be decomposed into two distinct classes, termed `gray' and `red', that each has a normal color distribution. The two color classes have different inclination distributions: red objects have lower inclinations than the gray ones. This trend holds for all dynamically excited TNO populations. Even in the worst-case scenario, biases in the discovery surveys cannot account for this trend: it is intrinsic to the TNO population. Considering that TNOs are the precursors of centaurs, and that their inclinations are roughly preserved as they become centaurs, our finding solves the conundrum of centaurs being the only outer Solar System population identified so far to exhibit this property (Tegler et al. 2016). The different orbital distributions of the gray and red dynamically excited TNOs provide strong evidence that their colors are due to different formation locations in a disk of planetesimals with a compositional gradient.


Introduction
Trans-Neptunian objects (TNOs) represent some of the most unaltered remnants of the planetary formation process. The surface of these objects can be used as tracers of their dynamical evolution which, in turn, helps us understand the current structure of the Kuiper Belt (e.g., Gladman 2005;Petit et al. 2011). Numerical models have been used to assess the evolution of TNOs from their formation locations onto their current orbits (e.g., Malhotra 1995;Gomes 2003;Levison et al. 2008;Nesvorny 2015aNesvorny , 2015b. Different formation locations might correspond to distinct surface compositions, though this remains disputed (Gil-Hutton 2002;Stern 2002;Trujillo & Brown 2002;Delsanti et al. 2004;Santos-Sanz et al. 2009). As such, both physical and dynamical properties must be considered to constrain the origins of TNOs. Unlike the majority of the known dwarf-planet-sized bodies, the smaller TNOs are too faint to be studied through optical and infrared spectroscopy. For these objects, one must instead rely on what broadband colors reveal.
The TNOs can be divided into two broad dynamical groups of objects: the dynamically quiescent cold classicals, and the dynamically excited hot population (Tegler & Romanishin 2000;Brown 2001). The latter comprises several subpopulations: (1) the hot classicals located at 40-48 au between Neptune's 3:2 and 2:1 mean-motion orbital resonances (MMRs; Levison & Stern 2001;Doressoundiram et al. 2002;Gladman et al. 2008), (2) the resonant objects locked in MMRs with Neptune (Malhotra 1995;Levison et al. 2008;Gladman et al. 2012), (3) the scattering population on dynamically unstable and highly eccentric orbits (Lykawka & Mukai 2007;Gladman et al. 2008;Gomes et al. 2008), and (4) the detached TNOs that experienced gravitational excitation in the past and are now on stable orbits unperturbed by Neptune (Gladman et al. 2002;Emel'yanenko et al. 2003;Gomes et al. 2005;Brasser & Schwamb 2015). Centaurs are thought to be former TNOs that evolved into inner orbits with semimajor axes and perihelia between the orbits of Jupiter and Neptune (e.g., . We consider them here as part of the larger group of dynamically excited TNOs. There are indications that colors and orbital properties of the dynamically excited TNOs are related. The largest dynamically excited TNOs uniformly exhibit neutral, ice-dominated surfaces (Fornasier et al. 2004; Barkume et al. 2008;Barucci et al. 2008Barucci et al. , 2011Brown 2012), possibly because they have sufficient gravity to retain their volatiles . In contrast, smaller dynamically excited TNOs have a bimodal color distribution, with classes termed gray and red that are divided around optical color (g-r)=0.78 (Tegler & Romanishin 1998Peixinho et al. 2003Peixinho et al. , 2012Peixinho et al. , 2015Tegler et al. , 2008Tegler et al. , 2016Fraser & Brown 2012;Lacerda et al. 2014;Wong & Brown 2017); this is fully quantified in Section 3. A correlation between colors and orbital inclinations was found for the population of classical TNOs (Tegler & Romanishin 2000;Hainaut & Delsanti 2002;Trujillo & Brown 2002;Doressoundiram et al. 2005;Peixinho et al. 2008). However, the data sets in these works were contaminated by the presence in their sample of peculiar compositional groups, including neutral-color Haumea family members Snodgrass et al. 2010) and red cold classical objects (Tegler & Romanishin 2000;Brown 2001). A careful removal of these objects diminishes the strength of the correlation between color and inclination, but the correlation remains statistically significant (Peixinho et al. 2015). More recently, Tegler et al. (2016) divided the observed population of centaurs into two groups by color and analyzed their orbital inclination distributions. They found that red centaurs have lower inclinations than the gray ones. As minor planets on trans-neptunian orbits are perturbed into the centaur region, they roughly preserve the inclinations they had before becoming centaurs (Volk & Malhotra 2013). As such, the different orbital distributions of the two color classes of centaurs should be reflective of a similar trend in the more distant TNOs. Such a trend might be the origin of the observed color versus inclination correlation of the hot classical population. Two color groups with different inclination distributions could indeed manifest as such a correlation. In short, different inclinations for the gray and red dynamically hot objects are anticipated but not yet observed in the Kuiper Belt.
We re-examine here the orbital inclination distribution of the two color classes of the dynamically excited TNOs. We use a large data set of high-quality color measurements that includes optical (g-r) colors obtained for 44 objects as part of the Colours of the Outer Solar System Origin Survey (Col-OSSOS; Schwamb et al. 2018), 25 of which are new measurements (Section 2). We further use robust color measurements extracted from several works reporting optical colors of TNOs, bringing our sample to 229 dynamically excited TNOs. We show in Section 3 that red dynamically excited TNOs have smaller orbital inclinations than their gray counterparts and that this trend appears to be present in every dynamical population. The centaurs are no longer the only population to exhibit this peculiarity. This has strong implications for the origins of the dynamically excited TNO populations (Section 4).

Sample Selection
Our analysis uses both newly obtained colors and color measurements from the literature for TNOs confirmed to be on dynamically excited orbits. We consider the subset of 44 dynamically excited objects from Outer Solar System Origin Survey (OSSOS; Bannister et al. 2016Bannister et al. , 2018 for which our team acquired optical and near-infrared colors as part of the Col-OSSOS survey (Schwamb et al. 2018). Measurements for 19 of these objects are previously published colors from Schwamb et al. (2018), and 25 are new color measurements reported here. Those include measurements for four objects that we previously published in Pike et al. (2017) and reprocessed here using the latest version of our photometric pipeline. We complement this Col-OSSOS data set with available highquality color measurements from the literature. These include optical colors from Peixinho et al. (2015), which compiles data from many authors, and Tegler et al. (2016), as well as Hubble Space Telescope (HST) data from Fraser & Brown (2012) and Fraser et al. (2015). Only colors observed with filters in the ∼500-850nm range were considered, i.e., Johnson-Cousins (V-R) and (R-I) and HST (F606w-F814w). One object from Peixinho et al. (2015), 2002GO9 (Crantor), has (R-I) color inconsistent with its (V-R) color and measured spectra (e.g., Alvarez-Candal et al. 2007). We consider only the (V-R) measurement reported in Peixinho et al. (2015) for that target.
Only small objects are considered for this analysis. The size transition from ice-rich surfaces on large TNOs to potentially volatile depleted surfaces on small TNOs is a matter of debate, but the transition is generally proposed to be at absolute magnitude H r ∼6 (Fraser & Brown 2012;Peixinho et al. 2012). Here, we find that including objects in the 5<H r <6 mag range does not affect the conclusions of our analysis (see Section 3). We therefore adopt a magnitude cut of H r > 5.0, corresponding to diameters of 440 km assuming a visible geometrical albedo of 0.09 (Lacerda et al. 2014).
We classified all of the TNOs in our sample using the procedure outlined in Gladman et al. (2008), which we briefly review here. Each object's orbit was fit using the Bernstein & Khushalani (2000) orbit-fitting routine. Starting from that best-fit orbit, a Monte-Carlo search was performed for the minimum and maximum semimajor axis orbits that do not produce orbit-fit residuals more than 1.5 times worse than those of the best-fit orbit. All three orbits (the best-fit, minimum-a, and maximum-a) were integrated forward in time for 10Myr and classified according to their evolution. A secure classification is one for which all three orbits behave the same way, and insecure classification is one where at least one orbit behaves differently. We have 202 secure classifications and 27 insecure classifications. Most of the insecure classifications are a result of an object being in or near an MMR rather than having a particularly large uncertainty in semimajor axis because all of the TNOs in our data set have reasonably long (>3 oppositions) observational arcs.
Objects belonging to compositionally distinct TNO populations, which we term as interlopers, were rejected from our sample. Specifically, we excluded all known Haumea family members, as well as gray-colored OSSOS/Col-OSSOS objects orbiting within the Haumea family cluster (confirmed members have 41.6<a<44.2 au, 0.08<e<0.19, 24.2<i<29°.1; Brown et al. 2007;Snodgrass et al. 2010). We further rejected cold classical objects, which are known for exhibiting different colors (Tegler & Romanishin 2000;Pike et al. 2017) and albedos (Brucker et al. 2009;Vilenius et al. 2014) from those of the dynamically excited TNOs. Gulbis et al. (2010) and Petit et al. (2011) found that the orbital inclination distribution of cold classicals can be best reproduced by a Gaussian of width 2 . 0 0.5 0.6  -+ or sin(i) times a Gaussian of width 2°.6, respectively. Adopting an inclination cut of i>5°ensures that contamination is minimal.
The same inclination cut was applied to all of the objects in our data set, independently of their dynamical classification, in order to perform a consistent analysis throughout all of the populations. We further verified that changing the inclination cut does not affect the conclusions of our analysis (Section 3). Centaurs with Tisserand parameter T J <3 were also excluded from our sample as their orbits are strongly coupled to Jupiter (e.g., Gladman et al. 2008) and thus have likely experienced larger orbital perturbations and changes in inclination (e.g., Volk & Malhotra 2013). The smallest (H>11) centaurs and scattering disk objects appear to have altered colors compared to the rest of the centaur and scattering population (Jewitt 2015), so they were also excluded. No other populations in our data set contain objects with H>11. Finally, we also rejected the two known objects with retrograde orbits because their origin population remains elusive (Gladman et al. 2009;Chen et al. 2016). Orbital elements and dynamical classifications for the resulting data set of 229 objects are reported in Appendix A.

Observations, Color Measurements, and Spectral Slopes
All of the new color measurements presented in this paper were acquired through the Col-OSSOS large program (PI: W. Fraser) on Gemini North between 2014 August and 2016 November. Col-OSSOS collects near-simultaneous g-, r-, and J-band photometry of a magnitude-limited (r<23.6) subset of the OSSOS sample. Optical measurements were acquired with the Gemini Multi-Object Spectrograph (GMOS; Hook et al. 2004), and the J-band sequence was obtained with the Near InfraRed Imager and Spectrometer (NIRI; Hodapp et al. 2003). We observed using a rgJgr sequence to account for any brightness variation due to light-curve effects.
A complete description of the data reduction and photometric extraction for Col-OSSOS is provided in Schwamb et al. (2018). We summarize it briefly here. Processing of the Col-OSSOS data was performed with the Gemini-IRAF package. We first used the bias images acquired as part of the GMOS calibration plan to remove the bias offset from the science frames. Master sky flats were produced from a set of twilight flats and used to flatten the science frames. Photometric measurements were achieved using the TRIPPy software package, which makes use of a pill-shaped aperture, built from the convolution product of a circular aperture with a line describing the direction and rate of motion of the target (Fraser et al. 2016). This approach minimizes the background contribution to the total flux in the photometric aperture.
Next, the Col-OSSOS data were calibrated to the Sloan Digital Sky Survey (SDSS) magnitude catalog (SDSS Collaboration et al. 2017) in cases where the field of view of the target was covered by SDSS; in other cases, we calibrated to the Pan-STARRS magnitude catalog (Magnier et al. 2013(Magnier et al. , 2016. We used in-frame catalog stars to determine the color conversion between the SDSS or Pan-STARRS and GMOS filter sets, using the relationships provided in Schwamb et al. (2018). We then fit the calibrated g and r magnitudes with a linear light-curve model assuming a constant (g-r) color across the observing sequence. The fit was weighted by the uncertainty of the individual photometric measurements calculated as the quadratic sum of the photometric uncertainty, the zero-point uncertainty, and the aperture correction uncertainty. The (g-r) color was then derived from our best fit, and the uncertainty on the color computed as the uncertainty on the fit. New colors obtained from our observations are reported in Table 1.
All of the color measurements from Col-OSSOS and the literature were converted to spectral slope (s), defined as the percentage increase in reflectance per 10 3 Åchange in wavelength normalized to 550nm, using the Synphot tool in the STSDAS software package. 13 Reported spectral slopes are the mean values from all the measurements of a target weighted by the inverse of the squared uncertainties. We consider only well-quantified measurements with a spectral slope uncertainty of Δs<12%/(10 3 Å), because larger error bars do not allow clear differentiation between the two color classes. The full list of measurements we consider is provided in Table 4 (Appendix A).

Red Dynamically Excited TNOs have Lower Orbital Inclinations
The whole data set is shown in Figure 1. The data suggest a relationship between color and inclination. In particular, there is a paucity of red TNOs with high orbital inclinations in our data set. In the following sections, we assess whether this trend is statistically robust. To do so, we divide our data set into two groups by color (Section 3.1) and test whether the two color groups have different inclination distributions (Section 3.2). An alternative approach would be to divide our data set into two inclination groups and compare the color distribution of these two groups. We take the first approach as our data set is clearly bimodal in color, as discussed in the next section. Effects of survey biases on our analysis are considered in Section 3.3.

Color Classification
Visually, the overall color distribution of all dynamically excited TNOs in our sample appears to be bimodal ( Figure 2). The test of normality rejects a single Gaussian distribution at the 99.8% confidence level (CL). We therefore test if it can be decomposed into multiple normal distributions. We calculated the Bayesian information criterion (BIC) for a set of N-component Gaussian mixture models (GMM), and confirmed that the best BIC is found for a two-component model. The GMM decomposes our sample into two normal distributions, with respective means s=11.5%/ (10 3 Å) and s 29.0% 10 3 = ( Å), and standard deviations s s = 5.8% 10 3 ( Å) and 8.4% 10 s 3 s = ( Å). This confirms a previously known property of the Kuiper Belt, which is that the dynamically excited red TNOs exhibit a broader range of colors than the gray ones. The origin of this feature remains unclear (e.g., Brown et al. 2011), and we do not investigate it further, as this goes beyond the scope of this paper. The two normal distributions intersect at s 20.6% 10 3 = ( Å), corresponding to optical colors (g-r)=0.78, (V-R)=0.56, and (F606w-F814w)=−0.11. Therefore, we classify all objects with s s 20.6% d + < 10 3 ( Å) as gray and all objects with s s 20.6% 10 3 d -> ( Å) as red, where δs is the uncertainty on the spectral slope. Objects at the boundary between the two color groups were not classified as they are ambiguous. Our classification is consistent with previous works that placed the bifurcation between gray and red TNOs at comparable spectral slope values ( Table 2).

Statistical Discrepancy between the Inclination
Distributions of Gray and Red TNOs

Statistical Tests
Four main tests were used to assess the orbital inclination versus spectral slope distribution of our data set. The onedimensional two-sample (2-s) Kolmogorov-Smirnov (KS; Chakravarti et al. 1967) and Anderson-Darling (AD; Stephens 1974) tests were used to determine whether the inclination distributions of the two TNO color groups defined in Section 3.1 come from a single parent distribution. These two tests are broadly similar and also complementary. The KS test is more sensitive to the center of the distributions, while the AD test gives more weight to the tails. The AD test is generally considered more powerful than the KS test. These two tests are nonparametric with no reference to a Gaussian. We therefore report 1 minus the p-value of the test (in percentage) as the level of confidence that the two distributions are different, rather than a number of σ deviation from a mean. We also used the student's T-test of means (Snedecor & Cochran 1989) to determine whether the mean inclinations of the gray and red TNOs are equal. Finally, we derived the Spearman coefficient (Zwillinger & Kokoska 2000) of our data set to search for a correlation between orbital inclination and spectral slope. Unlike the other three tests, the Spearman test does not require dividing our data set into two color groups. All of the results discussed in the following sections are summarized in Table 3. p-values were determined by computing each statistic for sets of inclination/spectral slope pairs constructed by independently randomly sampling the observed distribution of each parameter (Efron & Tibshirani 1993). We sampled with replacement so the same value of inclination or spectral slope could be sampled several times. The reported CL then indicates the probability that a random sample produces a larger test statistic than the observed population.

Full Data Set
We first consider the full data set as a whole. All statistical tests indicate a strong (CL>99.7%) discrepancy between the inclination distributions of the gray and red TNOs, as well as between their mean inclinations. The Spearman test indicates that the spectral slope anti-correlates with inclination at the same significance level ( Table 3). The color versus inclination trend is particularly obvious when comparing the ratios of grayto-red objects between low and high inclinations. From Figure 1, it appears that the inclination transition from where gray and red TNOs are roughly equally numbered to where red TNOs are almost nonexistent is located around i=21°. We therefore chose to compare the ratios of red-to-gray objects above and below this value. Only five TNOs-i.e., 11% of the objects-can be unambiguously classified as red at i>21° (  Figure 1). In contrast, red objects account for 42% of the lowinclination (5°<i<21°) population. To assess the likelihood of producing such an observation from a random population of objects, we bootstrap with replacement from the observed TNO samples of N g and N r objects, where N g and N r are the observed Notes. a OSSOS r′ and H r mag are at discovery, at an earlier observing geometry than the Col-OSSOS color measurement. b Object with previously published colors in Pike et al. (2017) and reprocessed here using the latest version of TRIPPy (Fraser et al. 2016). Δ = heliocentric range; r = distance at discovery; α = phase angle.
number of gray and red TNOs. We find that less than 0.1% of the simulations produce a higher or equal ratio of gray-to-red objects than in the observed sample of high-i (i>21°) TNOs. While the above tests confirm a trend in color versus inclination, they cannot determine whether this trend results from the two color classes having different inclination distributions or from a direct correlation between color and inclination. To assess the possibility of a correlation, we consider only objects with 5°<i<21°, where both color classes are thoroughly sampled by surveys ( Figure 1). Using the Spearman test, we find no evidence for a correlation between color and inclination for these objects ( Table 3). This indicates that the Spearman test finds a correlation for all (low-i and high-i) TNOs because the number of red TNOs drops off near 21°inclination. It follows that spectral slope and inclination of the dynamically excited TNOs are not directly correlated. Rather, red and gray objects have different inclination distributions. This opens interesting prospects for our understanding of the origin of TNO colors (see further discussion in Section 4).

Individual Populations
Here, we test whether the observed inclination distribution discrepancy between gray and red TNOs holds for the individual dynamical subpopulations of objects. We consider three distinct groups of dynamically excited TNOs. The first two groups, the hot classicals and the resonant objects, are each considered separately as they comprise enough objects to allow a statistically robust analysis. The remaining objects, including centaurs, scattering disk, and detached objects, are considered as a single group. We emphasize that considering those objects together does not mean we assume they have a common origin. We adopt this approach because these populations, especially the scattering disk and detached objects, are too sparsely sampled in our data set to be analyzed individually. Our goal here is to check whether the inclination discrepancy between gray and red TNOs holds for the remaining TNOs in our data set. Nevertheless, it is believed that the scattering disk directly feeds the centaur population . Since objects from the Kuiper Belt roughly preserve their inclinations as they become centaurs (Volk & Malhotra 2013), this seems the most sensible class merger. Figure 3 shows the inclination versus color distribution for the individual populations of TNOs.
Applying previous statistical tests to the individual groups, we find that the discrepancy between the mean inclinations of the gray and red TNOs holds with CL=96.1% on average for the hot classicals ( Table 3). The Spearman correlation between inclination and spectral slope is less significant (CL=93.2%). The KS and AD test return somewhat weaker results but this slightly depends on the magnitude and inclination cuts used to define our sample (see Appendix B). Our results therefore are in Figure 1. Orbital inclination vs. spectral slope of the dynamically excited TNOs and centaurs with i>5°. The vertical dashed line at s 20.6% 10 3 = ( Å) marks the limit between the two TNO color classes as determined by the Gaussian mixture model. The yellow star indicates the spectral slope value of the Sun. The reddest object in our sample is Centaur (5145)Pholus. The number of classified gray and red objects is indicated at the bottom of the corresponding panel. Unclassified objects are those for which the error bar on the color measurement intersects with the division between the two color classes. A smoothed density plot is used to highlight the density of the data points. Note the sparseness of red TNOs on highly inclined orbits. At 5°<i<21°, where both color classes are well sampled, inclinations and spectral slopes are not correlated; a correlation would be expected in the context of a collision origin for TNO colors (Section 4). Rather, TNOs appear to form two intrinsically different populations, each with its own orbital inclination and color distributions. Figure 2. Histogram of spectral slope (s) for our data set (229 objects in total). Application of a Gaussian mixture model decomposes our data set into two normal distributions (dashed lines) that intersect at s 20.6% 10 3 = ( Å). We adopt this value to classify our data set into the gray and red color classes. Table 2 shows that our value for this division is consistent with those of previous works. A similar trend is found in the resonant group and in the group of centaurs/scattering/detached objects at a comparable significance level ( Table 3). We also observe that out of the 102 resonant TNOs in our data set, where 49 are gray and 31 are red (the rest being unclassified), the 14 highest-inclination objects all belong to the gray class. When bootstrapping two random samples of 49 and 31 objects from the resonant population, we find that this occurs in less than 0.1% of the simulations, thus strongly supporting the conclusion that the gray and red resonant TNOs have different inclination distributions. Finding that resonant objects follow the same trend as other dynamical populations of TNOs is interesting. It implies that the initial orbital inclinations these objects had before being captured into orbital resonances with Neptune were not completely randomized by the capture process and the dynamical mechanisms, such as Kozai oscillation (e.g., Gomes et al. 2005), that take place in those resonances.
We stress that the group of centaurs/scattering/detached TNOs is strongly biased by the presence of centaurs in the sample. Excluding those objects completely removes the trend of gray and red TNOs having different orbital inclinations in that group. We note however that the scattering disk population does not contain any red objects at i>21°, similar to what is observed in the resonant population and, to a lesser extent, the hot classical (1 interloper) and the centaur (2 interlopers) populations ( Figure 3). Scattering disk objects might therefore follow the same trend as other TNO populations, but the current small number of scattering TNOs with well-measured colors prevents us from drawing any firm conclusion. Detached objects seem to constitute the only population that does not follow the trend, but their sample is also too small for conclusive analysis.

Summary of Statistical Tests
To summarize, we divided our data set into two color classes at 20.6%/(10 3 Å) and found a highly significant difference between the orbital inclination distributions of the gray and red TNOs. This discrepancy cannot be detected at the CL>99.7% significance level in the individual dynamical population of TNOs, but it does show up in each of these populations at lower statistical significance. This indicates that (1) this trend is common to every population of dynamically excited TNOs and (2) that our overall data set is not biased by one particular dynamical population. We emphasize that the results presented here slightly vary depending on the adopted magnitude and inclination cuts used to define the analyzed sample, as well as the spectral slope value of 20.6%/(10 3 Å) we adopt to differentiate between gray and red TNOs. We find however that a reasonable modification of those values does not affect the conclusions of our study in any significant way (see Appendix B).

Is the Observed Color-inclination Distribution an Artifact
of the Discovery Surveys?

Analytical Argument
The data set built for this work uses photometric measurements from a large variety of discovery surveys, each with its own detection biases and discovery efficiency. This raises the question as to whether these observing biases could be responsible for the observed inclination discrepancy between gray and red TNOs reported in this work. A discovery bias in favor of red TNOs comes from the combination of their higher optical albedos (Fraser et al. 2014;Lacerda et al. 2014) and redder colors with respect to gray objects. The observed ratio of red-to-gray objects depends on the filter used for the observations. Observations performed in broadband filters at longer wavelengths will discover a higher fraction of red TNOs than observations made at shorter wavelengths. This can introduce a color-inclination correlation if low-latitude discovery surveys (which are biased toward lower inclination TNOs) were performed at longer wavelengths than high-latitude surveys. To test whether this effect can be responsible for the observed distribution, let us consider the worst-case scenario, which would be if off-ecliptic TNO surveys were made in the V-band (550 nm), while on-ecliptic surveys were made in the R-band (660 nm).
The number of objects with albedos between a and a+da, between radii R and R+dR, and heliocentric distances between r and r+dr detected by the survey is is the Kuiper Belt radial distribution, h(a) is the TNO albedo distribution, and A is a constant. Assuming a power-law distribution in R, f R R q µ -( ) , and recasting the expression for R in terms of apparent magnitude m, it can be shown (Schwamb et al. 2018) that the number of objects with albedo a 0 aa 1 and magnitude m 0 mm 1 observed within the boundaries r 0 and r 1 of the Kuiper Belt is given by  where C is a constant, Δ is the geocentric distance to the object, and q 1 5 1.0 a = -» ( ) (Fraser et al. 2014) is the power-law slope of the apparent magnitude distribution. Here, we have assumed that the albedo, radial, and size distributions are independent. For simplicity, this expression does not reflect the apparent latitudinal and longitudinal variation of TNO density of the resonant populations. We are only interested in the ratio of gray-to-red objects observed within a given survey. Assuming that within each orbital class, separate color populations share the same spatial distribution, the latitudinal and longitudinal variations will be reflected in both color populations equally. Thus, for our derivation such variations can be ignored. In the case of a  Figure 1 for the individual dynamical populations of the dynamically excited TNOs. Top left: hot classicals; middle: resonant objects; right: combined group of centaurs, scattering disk, and detached objects. Bottom: the combined group is divided into centaurs (left), scattering disk objects (middle ), and detached objects (right). The absence of red TNOs on highly inclined orbits appears to be a common property of all dynamical populations. Centaurs are the only population for which different orbital inclinations of the two color classes have been shown so far (Tegler et al. 2016). Our data set reveals a similar trend to that of Tegler et al. (2016) for hot classicals, resonant objects, and centaurs, but with a lower statistical significance (see Table 3). The small number of scattering disk and detached objects in our sample does not permit us to statistically confirm the existence of a matching trend in those populations. However, we do observe a lack of red scattering disk objects on highly inclined orbits (i>21°). population composed of two color groups, where each group has its own albedo distribution, but similar size and radial distributions, and where the intrinsic number of objects in the two groups is given by A 2 =γ A 1 , where γ is a constant, the observed ratio of populations 1 and 2 is Here, we have assumed that the fraction of objects in the two groups is not a function of the discovery survey's observing depth. This is true as long as the size distribution can be modeled as a single power law ( f R R q µ -( ) ), i.e., as long as the limiting magnitude of the surveys is not sensitive to the break (or divot) magnitude (see, e.g., Shankman et al. 2013Shankman et al. , 2016Fraser et al. 2014;Lawler et al. 2018b). A fraction of objects in our data set, mainly centaurs and scattering disk objects, are fainter than the magnitude break measured for these objects (H=7.7 or 8.3; Lawler et al. 2018b). The presence of a break in the size distribution is explored numerically below in Section 3.3.2.
For simplicity, we consider the simple scenario where the two color groups each have a unique albedo, a 1 and a 2 . Then, the albedo distribution is given by the Dirac delta function h a a a n n d = -( ) ( ), and Equation (3) In the V-band (550 nm), gray and red TNOs have mean albedos of 6% and 12%, respectively (Fraser et al. 2014;Lacerda et al. 2014), implying a red-to-gray ratio of N N red gray = 5.66g. In this work, we derived mean spectral slopes of 11.4%/(10 3 Å) and 29.1%/(10 3 Å) for the gray and red color classes of TNOs, respectively. This implies albedos of 6.7% and 15.7% in the R-band (660 nm), and 8.41 Therefore, in the V spectral band the observed red-to-gray ratio is 67% that in the R spectral band.
In our data set, we measure a red-to-gray ratio of 0.821 for 5°<i<10°(23 red objects, 28 gray, and 17 unclassified). Assuming the worst-case scenario, where all on-ecliptic discovery surveys were performed in R and all off-ecliptic surveys were performed in V, the expected red-to-gray ratio at high inclination would be 0.550. Yet, the observed ratio for i>21°is 0.119 (5 red objects, 42 gray, and 6 unclassified). If considering g-(480 nm) and r-band (620 nm) filters instead of V and R, then the observed ratio in g is 56% that in r, and the expected ratio in our data set would be 0.450 at i>21°. Choosing different inclination cuts for the two subsamples slightly varies these numbers but does not affect the conclusion. Therefore, it appears that even the observationally heavily biased scenario considered here cannot account for the observed relationship between color and inclination.

Survey Simulation
We further test the effects of discovery biases on the color versus inclination distribution of TNOs using the OSSOS survey simulator (Bannister et al. 2016Lawler et al. 2018a). This simulator replicates the detection characteristics (pointings, detection and tracking efficiencies) of a given survey to determine which input model objects would be detected by the survey. We use this simulator to model a heavily biased survey composed of a single on-ecliptic block performed in the r-band and two high-latitude blocks, with mean inclinations of 12°and 22°, performed in g. All of the blocks are equally sensitive to the input objects in our simulation. Considering TNOs' mean (g-r)=0.78 color, the observational depth was set to m r =24.0 for the on-ecliptic block, and m g =24.75 for the off-ecliptic ones. We use as input populations the Canadian-France Ecliptic Plane Survey Each input object was cloned several times, allowing its orbital parameters a, e, and i to vary by 10% with respect to their nominal values, while randomizing its position angles (ω, Ω, and M). This cloning ensures that a sufficient number of simulated TNOs are randomly produced. We ran the simulator until a sample of 10,000 objects sampling the full parameter space was detected. Simulations were run independently for the on-ecliptic block and the high-latitude ones to ensure the same number of objects (5000) were detected by the blocks. Optical (g-r) colors were drawn from the double normal distribution measured in Section 3.1 and randomly attributed to the simulated objects. The input color ratio is irrelevant here; we are interested only in the relative color ratios between low-i and high-i objects in the output survey population. In order to determine the detectability of the simulated objects, an absolute magnitude distribution was assumed in the magnitude range of interest (5<H r <11). We used two different distributions: a single-slope distribution with an exponent of α=1.0 and a broken-slope distribution with slope exponents of α b =1.0 and α f =0.2 for the bright and faint end of the distribution, respectively (Fraser et al. 2014). The break magnitude of the broken-slope distribution was set to H 7.7 r = for the gray TNOs and H 6.9 r = for the red ones in order to account for their respective mean albedos of 6% and 12%.
The population of detected objects in our simulated survey reveals a trend of increasing color with decreasing inclination (Figure 4), in agreement with our survey setup. We first consider the population drawn from the single-slope distribution. The red-to-gray ratio of the objects found in the onecliptic block is 54% that in the off-ecliptic blocks for the L7 model and 62% for the K11 model, in good agreement with the predicted value of 56% derived in Section 3.3.1. However, a fraction of the high-i TNOs in our simulation are detected by the low-latitude block. Because of this, the color discrepancy between the low-i and high-i objects in our simulation is weaker than predicted analytically (Section 3.3.1). The red-togray ratio of the high-inclination (>21°) objects is 0.689 that of the low-inclination (i<10°) ones for the L7 model, and 0.854 for the K11 model. This strengthens our proposal that survey biases cannot account for the observed color versus inclination distribution. Finally, we find that the presence of a break in the magnitude distribution increases the red-to-gray ratio of the population detected in the high-latitude blocks relative to the on-ecliptic one; the on-ecliptic block ratio is 63% that in the off-ecliptic blocks for the L7 model, and 69% for the K11 model. This reinforces our conclusion that the observed color versus inclination correlation cannot be a result of the survey detection biases.

Real-life Surveys
We assess a subset of past surveys to further show that the correlation we report of color versus inclination distribution of the dynamically excited TNO populations is not the result of a discovery bias. Indeed, real surveys are far from being as biased as in the worst-case scenario we explored in Sections 3.3.1 and 3.3.2. We select the surveys with the most discoveries and an obtainable pointing history, including Chiang & Brown (1999), Gladman et al. (1998Gladman et al. ( , 2001 (Petit et al. 2011(Petit et al. , 2017, and OSSOS (Bannister et al. 2016. This list is by no means exhaustive, but we ensure it covers well the range of ecliptic latitudes and of filters sampled by past surveys, and it includes the surveys that have reported the most detections of TNOs with H r >5.0 (the DES, CFEPS, and OSSOS). We omit the two surveys of Schwamb et al. (2010) and Rabinowitz et al. (2012), which were up to large ecliptic latitudes, in red-band and in broadband red+blue filters, respectively, as we could not retrieve their precise pointing history. Figure 5 highlights that these surveys did not favor the detection of red TNOs on highly inclined orbits. On the contrary, most of the surveys performed far from the ecliptic plane used red-band filters, whereas the majority of surveys using blue-band and broadband blue+red filters were conducted close to the ecliptic. We should therefore expect objects detected on highly inclined orbits to be preferentially red compared to the low-inclination ones, which is contrary to the correlation we report. We also note that only a small fraction of the known TNOs were detected in the high-latitude fields, in agreement with their narrow orbital inclination distribution. This implies that even if a discovery bias exists in our data set, it probably only has a marginal effect on its color-inclination distribution.

Discussion
Considering the commonly accepted idea that centaurs are dynamically evolved TNOs, it would appear surprising if the inclination dichotomy between gray and red centaurs found by Tegler et al. (2016) did not occur in the Kuiper Belt. Indeed, how would such a signal appear after the scattering of centaurs from the Kuiper Belt? Here, we have shown that this dichotomy also exists in the Kuiper Belt, and that it appears to be a common feature to all dynamical classes of the dynamically hot TNOs. This removes the conundrum of centaurs being the only population observed so far to exhibit that property.
Finding different orbital properties for the gray and red dynamically hot TNOs opens up interesting prospects for understanding the origin of the color diversity of small TNOs. So far, two main interpretations have been put forward to explain this diversity. The first hypothesis suggests that all TNOs were originally similar, and evolutionary-collisional and resurfacing-processes altered them differently. The second hypothesis is that the two TNO color classes are composed of intrinsically different objects, with distinct compositions, that probably formed at different locations in the protoplanetary disk.
Orbital inclination is directly linked to the strength of the collisional environment. As such, previous reports of a color versus inclination correlation in the classical belt (Tegler & Romanishin 2000;Hainaut & Delsanti 2002;Trujillo & Brown 2002;Doressoundiram et al. 2005;Peixinho et al. 2008) were interpreted to be the result of collisional resurfacing processes (Gil-Hutton 2002;Stern 2002;Trujillo & Brown 2002;Moroz et al. 2003;Delsanti et al. 2004;Santos-Sanz et al. 2009). In that scenario, we would expect a direct correlation between color and inclination because the collisional kinetic energy directly correlates with inclination. Yet, we find no such correlation in our data set for orbital inclinations 5°<i<21°, where both color classes are well sampled by surveys (Section 3.2.2). As such, the strong Spearman correlation found for the complete set of TNOs most likely comes from the absence of red TNOs above the 21°i nclination. The lack of TNO color variation with rotation  and the lack of a correlation between color and orbital eccentricity (Thébault & Doressoundiram 2003;Thébault 2003) also appear to contradict a collisional origin for the inclination versus color relationship. Numerical simulations reveal that the kinetic energy TNOs receive from collisions always correlates better with eccentricity than inclination. Applying the same statistical tests as presented in Section 3.2 to the eccentricity versus color distribution in our data set reveals no correlation between these parameters.
It therefore appears that the dynamically excited TNOs formed in two intrinsically different populations, each with its . Inclination-dependent color ratio observed in our TNO sample (blue) compared to models of a heavily color-biased survey (Section 3.3.2; black). The survey model cannot reproduce the observed color ratio, implying it is intrinsic to the TNO population. Error bars on the TNO sample are derived from the number of color-unclassified TNOs (Section 3.1) in each inclination bin of 5°. The simulated survey observed a uniform color-ratio model on the ecliptic in r and at high latitudes in g. The indents in model color ratios at i≈12°, 22°correspond to the mean inclinations of our simulated survey highlatitude sky blocks. The dotted vertical lines indicate theoretical extreme colorratio cases observable for our simulation (Section 3.3.1). own orbital inclination and color distributions. Dynamical simulations have shown that the high-inclination TNOs could have formed closer to the Sun and been emplaced at their current location following a phase of planetary migration (Gomes 2003). Our study therefore suggests that gray TNOs accreted closer to the Sun compared to the red objects. During migration, their orbits were more excited by Neptune, which widened their orbital inclination distribution. The presence of gray class interlopers in the cold classical region, the so-called blue binaries , may however complicate this picture. The fact that these interlopers are all gray suggests that the gray class originally bordered the inner edge of the cold classical region (Schwamb et al. 2018). Future detailed modeling of planetary migration in a color-varying planetesimal disk will provide predictions which can be tested against our observed color-inclination distribution.
Different formation locations and emplacement efficiency as the origin of the two TNO color groups also appears compatible with the observation that gray TNOs are proportionally more numerous in the dynamically excited TNO populations. In our data set, gray TNOs represent 61%, 64%, 83%, and 87% of the total number of objects in the resonant, centaur, scattering, and detached populations, respectively, whereas they represent only 56% of the hot classicals. We however stress that this observation must be considered with extreme caution as it relies on a heavily biased data set built from various surveys with different detection biases. The Col-OSSOS survey will, upon completion, provide a complete fully debiased sample of TNO colors from which the intrinsic ratio of gray-to-red objects in each dynamical population will be derived (Schwamb et al. 2018).

Conclusion
We report new (g-r) colors for 25 small (H mag >5) dynamically excited TNOs observed in the Col-OSSOS survey. Combined with previously published color measurements from Col-OSSOS and other surveys, we consider the colors of a set of 229 dynamically excited TNOs and centaurs. We show that, when dividing our data set into two color classes at a spectral slope of s 20.6% 10 3 = ( Å), red objects have a lower orbital inclination distribution than their gray counterparts. This is true whether or not centaurs are included in our sample. This trend appears to be common to every dynamical population in the Kuiper Belt. The scattering and detached populations are too sparsely sampled to confirm a matching trend in those populations. Even in the worst-case scenario, observing biases in the discovery surveys cannot account for this correlation: it is intrinsic to the underlying TNO and centaur populations. Considering that TNOs are the precursors of centaurs, and that their inclinations are roughly preserved as they become centaurs (Volk & Malhotra 2013), our finding solves the conundrum of centaurs being the only outer solar system population identified so far to exhibit this property (Tegler et al. 2016). Finally, the observed inclination difference between gray and red TNOs tends to favor the hypothesis of different formation locations for these objects, and a compositional gradient in the original protodisk of planetesimals, rather than different evolutionary pathways.
The authors acknowledge the sacred nature of Maunakea, and appreciate the opportunity to observe from the mountain. This work is based on observations from the Large and Long Program GN-LP-1 (2014B through 2016B), obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), and Ministério da Ciência, Tecnologia e Inovação (Brazil). We thank the staff at Gemini North for their dedicated support of the Col-OSSOS program. Data was processed using the Gemini-IRAF package. STSDAS and PyRAF are products of the Space Telescope Science Institute, which is operated by AURA for NASA. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. The authors thank C. Shankman for sharing his Figure 5. On-sky pointing coordinates of the TNO discovery surveys listed in Section 3.3.3. We show the observing filter in which these surveys were performed, and schematically indicate their number of discoveries. Only the survey fields with positive detections are displayed for clarity. The dotted black curve and the blue curve indicate the location of the ecliptic and the galactic planes, respectively. Almost all of the surveys performed with blue-band and broadband blue+red filters were conducted close to the ecliptic plane. Most high-latitude observations come from the high-latitude component of CFEPS (Petit et al. 2017), which was performed with a red-band filter. These past observations imply that trans-Neptunian objects detected on highly inclined orbits should be preferentially red compared to the lowinclination ones. This is contrary to the correlation we report.  (Bertin & Arnouts 1996), STSDAS, TRIPPy (Fraser et al. 2016).

Appendix A Full Data Set
The complete set of 229 outer solar system objects on dynamically-excited orbits considered in this work is provided in Table 4, together with the orbital and physical properties of these objects.    We explore here the effects of sample selection and color classification on our statistical analysis. More specifically, we test whether our statistical results are sensitive to the magnitude and inclination cuts used to define our data set. We further test whether varying the spectral slope value used to divide our data set into the gray and red color classes make any difference. Large TNOs are characterized by different surfaces from those of small objects (see introduction, Section 1). The transition size between large and small TNO surfaces is usually proposed at H r ∼6 (Fraser & Brown 2012;Peixinho et al. 2012), i.e., at a diameter of 280 km assuming a visible albedo of 0.09 (Lacerda et al. 2014). In this work, our statistical analysis was performed including all objects with H r >5 (Section 2.1). We test here whether our analysis was biased by the presence in our data set of large 5<H r <6 TNOs by repeating our analysis considering only objects with H r >6. Our results are summarized at Table 5. They are statistically similar to those presented in Section 3.
We further test the effects of changing the inclination cut used for the selection of our data set. Our inclination cut of i=5°should ensure minimal contamination of our sample by cold classicals (Section 2.1). Let us assume, however, that our sample comprises a significant fraction of cold classical objects up to i=6°. Repeating our statistical analysis only for objects located at orbital inclinations of i>6°reveals similar results to that of Section 3 ( Table 5).
Finally, we test if our statistics are sensitive to the spectral slope value s lim used to divide our data set into the gray and red TNO classes. By adopting the two most extreme spectral slope values from Table 2 (17.4% and 23.8%/(10 3 Å)), we find that our results hold for the overall data set ( Table 6). Choosing the lowest spectral slope value weakens our results for the hot classical population as one previously unclassified high-i object falls into the red class. Choosing the highest value weakens our results for the centaur/scattering/detached objects as a significantly higher fraction of low-i objects fall into the gray class.

ORCID iDs
Michaël Marsset https:/ /orcid.org/0000-0001-8617-2425 Wesley C. Fraser https:/ /orcid.org/0000-0001-6680-6558 Rosemary E. Pike https://orcid.org/0000-0003-4797-5262 Michele T. Bannister https:/ /orcid.org/0000-0003-3257-4490 Note. Dynamical classification: cla=classical belt; sca=scattering disk; cen=centaur; x:y=resonators, where x and y indicate the specific mean-motion resonance. (i) Insecure classification. * 2007JK43 is nearly on a Uranus crossing orbit and has a semimajor axis that falls within the classical belt. It is classified as detached because no close encounter happens during the 10My integration of its orbit. Orbital elements are in the barycentric reference frame and were computed using the Bernstein & Khushalani (2000) orbit-fitting procedure. Dynamical classification uses the Gladman et al. (2008) classification scheme. See Bannister et al. (2016) for additional information. H r are the absolute magnitudes provided in "Ref." and converted to r-band using the Synphot/STSDAS tool. Peixinho et al. (2015) do not report the absolute magnitude for several objects in their data set. We retrieved these magnitudes from the database of Minor Bodies in the Outer Solar System (MBOSS) (Hainaut et al. 2012) when available, or from the original papers (Boehnhardt et al. 2002;Benecchi et al. 2009Benecchi et al. , 2011Sheppard 2010)  (This table is available in its entirety in machine-readable form.) Note. Same acronyms as those in Table 3. Note. Same acronyms as those in Table 3.