OSSOS. XXIII. 2013 VZ70 and the Temporary Coorbitals of the Giant Planet

We present the discovery of 2013 VZ70, the first known horseshoe coorbital companion of Saturn. Observed by the Outer Solar System Origins Survey (OSSOS) for 4.5 years, the orbit of 2013 VZ70 is determined to high precision, revealing that it currently is in `horseshoe' libration with the planet. This coorbital motion will last at least thousands of years but ends ~10 kyr from now; 2013 VZ70 is thus another example of the already-known `transient coorbital' populations of the giant planets, with this being the first known prograde example for Saturn (temporary retrograde coorbitals are known for Jupiter and Saturn). We present a theoretical steady state model of the scattering population of trans-Neptunian origin in the giant planet region (2--34 au), including the temporary coorbital populations of the four giant planets. We expose this model to observational biases using survey simulations in order to compare the model to the real detections made by a set of well-characterized outer Solar System surveys. While the observed number of coorbitals relative to the scattering population is higher than predicted, we show that the number of observed transient coorbitals of each giant planet relative to each other is consistent with a transneptunian source.

the mean longitude, P denotes the planet, Ω is the longitude of the ascending node, ω the argument of pericenter and M the mean anomaly. The resonant angle φ 11 must librate rather than circulate (ie. φ 11 must occupy a bounded range) in order for an object to be considered to be in coorbital resonance. Like other n : 1 resonances, the 1:1 mean-motion resonance includes multiple libration islands; objects in these islands are called leading Trojans (mean φ 11 = +60 • ), trailing Trojans ( φ 11 = 300 • = −60 • ), quasi-satellites ( φ 11 = 0 • ) or horseshoe coorbitals ( φ 11 = 180 • ). The motion of Trojans librate around one of the L4 or L5 Lagrangian points, while the path of horseshoe coorbitals encompass all of the L3, L4 and L5 Lagrangian points; quasi-satellites appear to orbit the planet (while not actually being bound to it). Quasi-satellites and horseshoe coorbitals are almost always unstable and thus temporary (eg. Mikkola et al. 2006;Ćuk et al. 2012;Jedicke et al. 2018) with the exception of Saturn's moons Epimetheus and Janus, which are horseshoe coorbitals of each other (Fountain & Larson 1978). Greenstreet et al. (2020) and Li et al. (2018) discuss the existence of high inclination (i > 90 • ) objects temporarily trapped in a 1:-1 retrograde "coorbital" resonance with Jupiter and Saturn, respectively, although these are not coorbitals in the traditional sense described above; since they orbit the Sun in the opposite direction than the planet, retrograde coorbitals are not protected from close approaches with the planet the way that prograde coorbitals are, nor do the resonant island librations (ie. Trojan, horseshoe, quasi-satellite motion) behave in the traditional sense in the retrograde configuration.
For many planets, the coorbital phase space is unstable due to perturbations from neighboring planets (eg. Nesvorný & Dones 2002;Dvorak et al. 2010). Innanen & Mikkola (1989) first suggested, at a time when only the Jovian Trojans were known, that populations of objects in stable 1:1 resonance with each of the other giant planets may exist; their analysis showed that the exact Lagrangian points are unstable for Saturn, but that Trojans farther from the resonance center (featuring larger libration amplitudes) could be stable for at least 10 Myr. These results were confirmed by Holman & Wisdom (1993). Using longer timescales than previous studies, de la Barre et al. (1996) specifically studied the stability of Saturnian Trojans and found that Saturnian Trojans could only be long-term (> 428 Myr) stable with very specific conditions: very small eccentricity (<0.028), φ 11 libration amplitude greater than 80 • , ω libration about a point 45 • ahead of Saturn's ω, and constraints on the timing of the maximum eccentricity relative to the timing of Jupiter's maximum eccentricity, so that Jupiter and the Trojans do not approach close enough to dislodge the Trojan from Saturn's 1:1 resonance. Nesvorný & Dones (2002) showed that while Neptunian Trojans may have only been depleted by a factor of 2 over the age of the Solar System, the Saturnian Trojans would have been depleted by a factor of 100. Studying the cause of the instability of Saturnian Trojans, Marzari & Scholl (2000) and Hou et al. (2014) found that the instability is caused by interactions between mean motion and secular resonances. Huang et al. (2019) investigated the stability of retrograde Saturnian coorbitals and found that they are always unstable due to an overlap with the ν 5 and ν 6 secular resonances.
Given these destabilizing factors, causing any primordial population to have been mostly depleted and allowing only small niches to be long term stable, it is not surprising that no long-term stable Saturnian Trojans have been discovered to date.
Only Mars, Jupiter and Neptune have known populations of long-term (>Gyr) stable Trojans (which thus might be primordial) (Wolf 1906;Bowell et al. 1990;Levison et al. 1997;Marzari et al. 2003;Scholl et al. 2005). These long-term stable Trojan populations are important for understanding planet formation processes. As a few examples: Polishook et al. (2017) suggested that the Martian Trojans are likely to be impact ejecta from Mars, and used the mass of the current Trojan cloud to constrain how much Mars' orbit could have evolved during the phase of collisions. Morbidelli et al. (2005) showed that in order to reproduce the wide inclination-distribution of the Jovian Trojans, the Trojans must have been captured from an excited disk during a migration phase rather than having formed in place together with Jupiter. Nesvorný et al. (2013) demonstrated that a sudden displacement of Jupiter's semi-major axis, can explain the asymmetry seen between the L4 and L5 clouds and use the mass of the Jovian Trojan clouds to estimate the mass of the primordial planetesimal disk. Gomes & Nesvorný (2016) used the observed mass of Neptunian Trojans to infer that Neptune migrated slightly past its current location and then back, destabilizing the cloud, as we would otherwise observe a more massive cloud. Parker (2015) demonstrated that if Neptune's migration and eccentricity-damping was fast, the disk that it migrated into and captured Trojans from must already have been dynamically excited prior to Neptune's arrival in order to reproduced our observed orbital distribution.
While only three planets are known to have long-term stable Trojans, scattering objects (scattering TNOs, Centaurs and even some objects originating in the asteroid belt 1 ) can become temporary coorbitals, transiently captured into unstable resonance (Alexandersen et al. 2013;Greenstreet et al. 2020). All Solar System planets except Mercury, Mars and Jupiter now have known populations of temporary coorbitals on prograde (i < 90 • ) orbits (Wiegert et al. 1998;Mikkola et al. 2004;Karlsson 2004;Horner & Lykawka 2012;Alexandersen et al. 2013;Greenstreet et al. 2020). Temporary "sticking" like this also occurs in other resonances (eg. Duncan & Levison 1997;Tsiganis et al. 2000;Alvarez-Candal & Roig 2005;Lykawka & Mukai 2007;Yu et al. 2018;Volk et al. 2018). While long-term stable Trojans inform us of conditions in the time of planet formation and migration, the temporarily captured coorbitals inform us about properties of the scattering population. For example, Alexandersen et al. (2013) confirmed the Shankman et al. (2013) finding that the size distribution of the scattering population must have a transition in order to explain the observed ratio of small nearby scattering object (including Uranian coorbitals) and larger more distant ones (including Neptunian coorbitals) Horner & Wyn Evans (2006) integrated the Centaurs known at the time, demonstrating that Centaurs do indeed get captured into temporary coorbital resonance with the giant planets, claiming that Jupiter should have by far the most temporary coorbitals, followed by Saturn and hardly any for Uranus and Neptune. Alexandersen et al. (2013) pointed out that using the known centaurs as the starting sample is biased towards having more objects nearer the Sun, and thus more captures for the inner giant planets, resulting in a disagreement with the sample of at-the-time known temporary coorbitals; they instead used a model that started with scattering TNOs that scatter inwards to become Centaurs and temporary coorbitals, to demonstrate that a TNO origin can explain the distribution of the temporary coorbitals of Neptune and Uranus.
1 For the rest of this work, "scattering objects" will be considered synonymous with scattering TNOs and Centaurs of TNO origin, ignoring Centaurs originating from the asteroid belt, unless asteroidal origin is explicitly mentioned.
In this paper we describe the discovery of the first known Saturnian horseshoe coorbital, 2013 VZ 70 , and demonstrate its temporary nature (Section 2). Furthermore we expand upon the analysis of Alexandersen et al. (2013) to analyze the populations of temporary coorbitals of all four giant planets, in an attempt to demonstrate the likely origin of 2013 VZ 70 and similar objects. We use numerical integrations to construct a steady-state distribution model of the scattering Trans-Neptunian Objects (TNOs) and temporary coorbitals of the giant planets (Section 3). Lastly we use survey simulations, exposing our model to the survey biases of a well-understood set of surveys, in order to compare our theoretical predictions to real detections of this population (Section 4).  (Folkner et al. 2014).
This orbit is very close to that of Saturn, although the two bodies are seperated by ∼ 180 • on the sky. From dynamical integrations we found that 2013 VZ 70 is in fact in the 1:1 mean motion resonance with Saturn, in a horseshoe configuration (see Figure 1). However, the best fit clone only remains resonant for about 11 kyr before leaving the resonance and re-joining the scattering population.
We investigated whether the best fit orbit could be near a stability boundary by generating orbit clones from appropriate resampling of our astrometry, allowing us to test whether any orbit consistent with the astrometry featured long-term stability. Each clone was produced by resampling all the astrometry (using a normal distribution with standard deviation equal to the mean residual of the to Saturn's mean motion (red dot). The start and end points are marked. The time between integration outputs (blue dots) is ∼ 0.33 years. Note that when the object appears near Saturn in this planar, meanmotion subtracted projection, it is not actually close to the planet due to the vertical motion caused by the orbital inclination and the fact that Saturn's true location does oscillate around the marked mean location.
The small cycles are caused by the eccentricity of the object's orbit, causing one little loop-and-shift motion for every orbit around the Sun. Each local minimum in distance from the center of the plot corresponds to a perihelion passage, and each local maximum corresponds to an aphelion passage; the small loops occur when the object's distance is close to Saturn's mean heliocentric distance, while the large shifts occur when the object is at a distance substantially different from Saturn, thus moving faster or slower around the Sun than Saturn does. It is thus clear from this figure that the horseshoe libration period at this time is roughly 30 orbital periods, although the libration period does vary slightly while always remaining near ∼ 1 kyr.
best fit, 0. 146) and fitting a new orbit. This process was repeated 10,000 times using Find_Orb.
The distribution of orbits generated by this process explicitly shows how the uncertainty of some of the orbital parameters are strongly coupled, as can be seen in Figure 2. From these 10,000 clones, we identified the most extreme orbits (largest and smallest value of each parameter) and integrated these 8 clones (labelled in Figure 2) as well as the best fit orbit. These dynamical integrations were done using Rebound (Rein & Liu 2012) with the WHFast (Wisdom & Holman 1991;Kinoshita et al. 1991;Rein & Tamayo 2015) symplectic integrator. The eight major planets and Pluto 2 were included as massive perturbers and an integration step size of 5% of Mercury's initial orbital period ( 0.012 yr 4.39 days) was used, while output was saved approximately 3 times per year. As can be seen in Figures 3 and 4, the future evolution of all of the clones involve an initial period in coorbital resonance, but all clones leave the resonance between 6 kyr and 26 kyr from now. The large range of resonance exit times is due to the highly chaotic nature of the orbit. We used a second set of numerical integrations, where clones were displaced infinitesimally (10 −13 -10 −12 AU, or 1.5-15 cm) relative to each of the above clones, to estimate the object's current timescale for chaotic divergence (the Lyapunov time scale); we found this to be 410±60 yr 3 . 2013 VZ 70 is thus definitely in coorbital resonance now, but the chaotic nature of the orbit means that the duration of this temporary resonance capture will likely not be constrained further, even with additional observations.

DERIVING THE STEADY STATE ORBITAL DISTRIBUTION
We proceed to investigate the potential origin of temporary coorbitals like 2013 VZ 70 . We model the source of the giant planet temporary coorbitals and investigate their detectability in characterised surveys. We produced a steady-state distribution of scattering objects in the a < 34 au region from orbital integrations similar to those used in Alexandersen et al. (2013); the details of those integrations can be found in the supplementary material of that paper. We primarily outline the deviations from those used in the previous paper below.
2 Pluto was primarily included as a test to ensure that the system was set up correctly, not because we expect the mass of Pluto to have any influence on the outcome of the integration. However, since Pluto's mass is known, there was no reason to not include it. Pluto was confirmed to be resonating in the 3:2 resonance with Neptune in our integrations, as expected. 3 A Jupyter notebook demonstrating how the Lyapunov time scale was calculated is available at DOI: 10.11570/21.0008.

Figure 2.
Orbital elements of 10,000 fits to resampled astrometry. Each clone was generated using Find_Orb's MCMC feature, using a Gaussian noise equal to the mean residual of the best fit orbit (0. 146).
The clones are color coded by semi-major axis (which is also the x-axis of the left most row) to give an additional rough indicator of its correlation with the other orbital elements. The best fit orbit is marked with a black dot, and orbits with either the smallest or largest value of one of the parameters are marked with a letter ((a-h, not to be confused with any orbital elements). This figure demonstrates how the uncertainties on the different parameters are related; most parameters are strongly coupled, while i and Ω have only weak or no coupling with other parameters.
To perform the dynamical integrations, we used the N-body code SWIFT-RMVS4 (provided by Hal Levison, based on the original SWIFT (Levison & Duncan 1994)) with a base time step of 25 days and an output interval of 50 years for the orbital elements of the planets and any particle which with the resonance. The "stair step" patterns occurs at the time when φ 11 is close to 0 • /360 • , which is the time the coorbital is closest to the planet; the close approach causes the switch from a slightly-larger-thanthe-planet's semi-major axis to a slightly-smaller-than-the-planet's (and vice versa) that ensures that the planet/coorbital never overtake each other. The close approaches also imparts small changes in the other orbital elements, seen as the "stair step" pattern.
at the moment had a < 34 au. The gravitational influences of the four giant planets and the Sun were included. The system starts with 8500 particles, derived from the 34 au < a < 200 au scattering Since 1 Gyr is substantially longer than the dynamical lifetime of Centaurs and scattering TNOs, we thoroughly sample the a <34 au phase space, despite the limited number of initial particles. The output is combined along the time-axis to produce a distribution of approximately 300 million sets of orbital elements. As in Alexandersen et al. (2013) we confirm that the distribution in the first 100 Myr is similar enough to the distribution in the following 900 Myr (because the a <34 au region is populated very quickly despite starting off empty) that we can treat the distribution as a whole as being in steady state. We also ran a similar simulation with particles drawn from the modified version of the Kaib et al. (2011) model also used in Alexandersen et al. (2013) and Shankman et al. (2013); this modified version was generated with an assumption of the primordial planetesimal disk being more dynamically excited. As in Alexandersen et al. (2013), we find that using the standard Kaib et al. (2011) model and the modified version as our starting condition makes little quantitative difference on the end results. For the rest of this work we will therefore only be referring to the results from the simulations using the standard Kaib et al. (2011) model.
The method for determining coorbital behavior in the particle histories is also very similar to that used in Alexandersen et al. (2013), with some small modifications. To diagnose whether particles are coorbital, the orbital histories (at 50 year output intervals) were scanned using a running window 30 kyr long for Uranus/Neptune and 5 kyr long for Jupiter/Saturn; this window size was chosen to be several times longer than the typical Trojan libration period at the given planet (∼ 1 kyr for Saturn as seen in Figure 4). A particle was classified as a coorbital if, within the running window, both its average semimajor axis was less than 0.2 AU from the average semimajor axis of a given planet and no individual semimajor axis value deviated more than R H from that of the planet. Here Our results are in good agreement with those for Uranus and Neptune in Alexandersen et al. (2013). Table 1 contains the fraction of the steady-state population in coorbital motion with each of the giant planets at any given time, as well as the distribution of coorbitals between horseshoe, Trojan and quasi-satellite orbits. The coorbital fraction for Uranus and Neptune are slightly higher than in Alexandersen et al. (2013), despite very similar methodology; however, these results agree within their expected accuracy. The capture fraction decreases from Neptune through to Jupiter (with almost scattering objects is beyond Neptune and that the dynamical timescales (orbital period, libration timescale) are longer farther from the Sun. An interesting result is that while Neptune's and Uranus' coorbitals are roughly equally distributed between horseshoe and Trojan coorbitals, Saturn seems to very preferentially capture scattering objects into horseshoe orbits, and Jupiter has a much larger fraction of quasi-satellites than any of the other planets. Table 1  it is noteworthy that the median number of orbital periods for coorbital captures for all four planets are within a factor of 2.5 of each other. However, going Neptune to Jupiter, particles are increasingly unlikely to have multiple captures, presumably due to the increasing ability of the planet to scatter the object to large semi-major axis; this results in particles on average spending both more total time and more total orbital periods in coorbital motion with Uranus and Neptune than Saturn and Jupiter.

COMPARING THEORY AND OBSERVATIONS
In order to compare our dynamical model to the real detections, we run the model through the OSSOS Survey Simulator (Bannister et al. 2018;Lawler et al. 2018a). The survey simulator generates one object at a time (with an orbit drawn from the dynamical model and an H-magnitude drawn from a parametric model discussed later) and assesses whether the object would have been discovered by the input surveys. We used all of the characterized surveys with sufficient characterization available Table 1. Steady state fractions of the a < 34 AU, q > 2 AU scattering objects that are in temporary co-orbital resonance with the giant planets. For reference, "Horseshoe" coorbitals librate about φ 11 = 180 • , "Trojan" librate about φ 11 = 60 • and 300 • , and "Quasi-satellites" librate about φ 11 = 0 • . Also listed are the mean, median and maximum duration of such captures seen in our simulations, and the median lifetime divided by the orbital period of the associated planet. Lastly the mean, median and maximum number of captures experienced by a particle that is trapped by the planet at least once.  (Bannister et al. 2018). These surveys combined will be referred to as OSSOS++.

Orbital distribution
For our survey simulations, it is preferable to have orbital distribution functions rather than an orbital distribution composed of a fixed number of discrete particles. This allows for the simulator to be run for as long as necessary, without producing duplicate identical particles. We have set up independent distributions for the coorbitals and the scattering objects, as described below, both inspired by the distribution seen in the integrations discussed in Section 3. Our model files and scripts for use with the OSSOS Survey Simulator (Lawler et al. 2018a) are provided at DOI: 10.11570/21.0008 for anybody curious to use this model distribution.

Scattering objects
We use the output from the Section 3 integrations, taking every particle's orbit at every time step and binning them using bin sizes of 0.5 au, 0.02 and 2 • .0 in a, e and i space. The survey simulator reads this binned table, randomly selects a bin weighted by the number of particles that went into the bin, and then randomly assigns a, e and i from a uniform distribution within the bin. Ω, ω and M are all assigned randomly from a uniform distribution from 0 • to 360 • , since the orientations of scattering objects' orbits are random. This process allows us to draw essentially infinite unique particles that follow a distribution consistent with the steady-state distribution from Section 3.

Coorbital objects
We cannot simply bin the coorbital distributions as we did for the scattering distribution. The numbers of coorbitals in the Section 3 integrations are low (particularly for Jupiter), and the number of dimensions we would need to bin is higher since the resonant angle φ 11 is also important for the coorbital distribution. Instead, we opted to use parametric distributions, fitted to the distributions seen in Section 3.
In this simplified parametric model, the semi-major axis of the coorbital is always set equal to that of the planet, since the few tenths of au variability do not influence detectability by sky surveys as much as the details of the eccentricity and inclination distribution. The eccentricity is modelled with a normal distribution, centred at 0 with a width w e , multiplied by sin 2 (e), truncated to [e min , e max ] 4 : This functional form has little physical motivation and was merely chosen as it in the end provides a good fit to the distribution seen in our integrations. The inclination is modelled as a normal distribution, with centre at 0 • and a width w i , multiplied by sin(i), truncated at i max : This is simply a Normal distribution modified to account for the spherical coordinate system. Lastly, Ω and M are chosen randomly from a uniform distribution [0 • , 360 • ) while ω is calculated from φ 11 , the value of which depends on the type of coorbital. The different types of coorbitals are generated using the ratios in Table 1. The details of the selection of a φ 11 value is similar to that used in   Table 2.

Absolute magnitude distributions
The Solar System absolute magnitude (H) distribution of the TNOs is not well constrained for objects fainter than about H r ≈ 8.0, although it is clear that there is a transition from a steep to shallower slope somewhere in 7.5 ≤ H r ≤ 9 (Sheppard & Trujillo 2010;Shankman et al. 2013;Fraser et al. 2014;Alexandersen et al. 2016;Lawler et al. 2018b). The scattering objects provide a clue to the small-end distribution, as many of these reach distances closer to the Sun, allowing us to more easily detect smaller objects. Lawler et al. (2018b) carefully analysed the size distribution of the scattering objects in OSSOS++; since our sample is a subset of their sample (we only use objects with a < 34 au), we will directly apply the two magnitude distributions favored by Lawler et al.  Table 3. Details of the sample (29 objects) used in this work. MPC name denotes the Minor Planet Center designation for the TNO, while O++ name is the internal designation used within the OSSOS++ surveys.

Population estimate
Cls is the classification of the object, where coorbitals of Saturn, Uranus and Neptune are indicated with the initial of the planet (S, U, N), with subscripted H, 4 or 5 for horseshoe coorbitals, leading Trojans and trailing Trojans, respectively; finally C indicates a non-coorbital Centaur/scattering object. Mag is the magnitude at discovery in the filter F, while H is the absolute magnitude in that same filter. The J2000 barycentric distance, semi-major axis, eccentricity and inclination are shown in d, a, e, i, respectively; for both d and i the uncertainty is 1 on the last digit or smaller, and has therefore been omitted.   We predict a population estimate for the scattering objects with a < 34 au of trans-Neptunian origin based on our model and the real detections. The OSSOS++ surveys discovered a total of 29 scattering objects with a < 34 au (including 4 temporary coorbitals), listed in Table 3. For the purposses of this work, 2013 VZ 70 is included in this sample, despite its uncharacterized status, as discussed in subsection 4.5. We thus ran the survey simulation with our scattering model (see Section 4.1.1) as input until it detected 29 objects, recorded how many objects had been drawn from the model, and repeated 1000 times to measure the uncertainty for the population estimate. For the divot and knee H distributions respectively, we predict the existence of (2.1±0.2)×10 7 and (4.9±1.0)×10 6 scattering TNOs with a < 34 au and H r < 19. Given the size and orbit distribution, most of these are small objects beyond 30 au and thus far beyond the detectability both of the surveys we consider here and similar-depth future surveys like the upcoming Legacy Survey of Space and Time (LSST) on the Vera Rubin Observatory.

Expected versus detected numbers
Using the population estimate of the a < 34 au scattering objects as measured in Section 4.3, we predict the number of temporary coorbitals of TNO origin that OSSOS++ should have detected. This is done by running the survey simulator for each planet's coorbital population separately (using the coorbital model defined in Section 4.1.2), inputting a fixed number of coorbital particles (equal to the total scattering object population estimate found in Section 4.3 multiplied by the coorbital fraction for the given planet as found in Section 3) and recording the number of detections, repeating 1000 times to sample the distribution. We find that for both the divot and knee distribution and for each planet, the most common (expected) value of temporary coorbital detections is zero. However, the probabilities of getting zero detected temporary coorbitals for Saturn, Uranus and Neptune are 84%, 74% and 71% (divot) or 91%, 73% and 59% (knee). The probability of getting zero detections for all three planets in these surveys is thus less than 50%. In other words, more often than not, we would expect OSSOS++ to detect at least one giant planet coorbital beyond Jupiter (the case of Jupiter is discussed in the next paragraph). From the distribution of simulated detections, we find that the detection of four coorbitals (as in the real surveys) is unlikely, at a probability of 0.8%, but not completely implausible. We expand on this below.
For Jupiter, the chance of zero detections is > 99.99% due to the rate of motion cuts imposed on/by the moving object detection algorithms of the OSSOS++ surveys; only the most eccentric Jovian coorbitals would have been detectable at aphelion (and only in a few fields). It is thus entirely reasonable that OSSOS++ found no Jovian coorbitals, neither temporary nor long term stable; these surveys were simply not sensitive to objects at those distances. We note that there is one known temporary retrograde (i > 90 • ) coorbital of Jupiter, 2015 BZ 509 (514107) Ka'epaoka'awela (Wiegert et al. 2017), whose origin is, according to Greenstreet et al. (2020), most likely the main asteroid belt and not the trans-neptunian/scattering object population. Our simulations produce no retrograde coorbitals of any of the planets, supporting that Ka'epaoka'awela likely originates from the asteroid belt and not the transneptunian region. of the detected a < 34 au scattering objects should be coorbitals, whereas the four real coorbitals make up 14% of detections (4 of 29); the observed fraction of the a < 34 au scattering objects that are in temporary coorbital resonance is thus ∼ 5 times higher than expected. Before the OSSOS survey, which was by far the most sensitive survey of the ensemble and discovered over 80% of the OSSOS++ TNOs and scattering objects, 60% were coorbital (3 of 5), so it would appear that the initially high fraction of coorbitals detected in the earlier surveys in our set was a fluke, and that the ratio is approaching the theoretical value predicted above as the observed sample increases. We thus do not feel it justified to hypothesize additional sources for the temporary coorbital population at this time. While we cannot rule out that the population of temporary coorbitals, particularly for Jupiter, is supplemented from other sources such as the asteroid belt and primordial Jovian Trojans, Greenstreet et al. (2020) finds that for Jovian temporary coorbitals, the asteroid belt is only the dominant source for retrograde (i > 90 • ) coorbitals, which they estimate comprise 1% of the temporary coorbital population. It is unlikely that the asteroid belt is a dominant source for the outer planets if it isn't for Jupiter. The contribution of the asteroid belt to the steady state temporary coorbital distribution of the giant planets is thus insignificant, and we are likely not missing any important source population in producing our population/detection estimates.

Caveat
While the orbit of 2013 VZ 70 was well determined by the OSSOS observations, it is not part of the characterized OSSOS dataset. 2013 VZ 70 was discovered in images taken in a "failed" observing sequence from 2013B (failed due to poor image quality and the sequence not being completed), which was thus not used for the characterised (ie. well understood) part of the OSSOS survey. This failed sequence, which should have been 30 high-quality images of 10 fields (half of the OSSOS "H" block), only obtained low-quality (limiting m r ≈ 23.5) images of 6 fields. A TNO search of these images was conducted (discovering 2013 VZ 70 ) to facilitate follow-up observations (color and light curve measurements), but this shallow search was never characterised due to the expectation that everything would be rediscovered in an eventual high-quality discovery sequence. A high-quality observing sequence of the full set of H-block fields was successfully observed in 2014B, with limiting magnitude m r = 24.67, which was used for the characterised search. However, as a year had passed, 2013 VZ 70 had already left the field due to its large rate of motion; unlike all other objects discovered in the failed 2013B sequence, 2013 VZ 70 was thus not re-discovered in the characterised discovery images. As such, 2013 VZ 70 is not part of the characterised sample of the survey, as that sample only includes objects discovered in specific images on specific nights through a carefully characterised process. However, because the failed discovery sequence points at the same area of the sky as parts of the characterised survey and it is a small minority of the total observed fields, it would make hardly any difference on the discovery biases whether these particular images are included in the characterization or not. From our simulations in Section 4 we can see that only about 8% of simulated detections of theoretical Saturnian Coorbitals were discovered in the OSSOS H block; this block is thus not in a crucial location for discovering Saturnian coorbitals in any way. The fact that the only Saturnian coorbital to have been discovered in OSSOS++ was among the very small minority of those surveys' total discoveries that were not characterised thus appears to be a low-probability event. We can therefore treat 2013 VZ 70 as effectively being part of the characterised survey for the purposes of this work, with the warning that this approach should not be used for other non-characterised objects from these surveys; most other objects are non-characterised for other reasons, mostly for being fainter than the well-measured part of the detection efficiency function. That being said, ignoring 2013 VZ 70 from the sample on grounds of being uncharacterized would brings the coorbital to total scattering ratio down to 11% (3 out of 28), closer to the 3% predicted in the previous section. Uranus' coorbitals should be roughly equally distributed between horseshoe and Trojan coorbitals, Saturn very preferentially captures scattering objects into horseshoe orbits, and Jupiter should have a much larger fraction of its temporary coorbitals be quasi-satellites than any of the other planets.
Accounting for observing biases in a set of well-characterized surveys (CFEPS (Petit et al. 2011), HiLat (Petit et al. 2015), the Alexandersen et al. (2016) survey, and OSSOS (Bannister et al. 2018)) we find that the fraction of a < 34 au scattering objects that are in temporary coorbital motion is higher in the real observations (13.7%) than in simulated observations (2.9%). However, for the distribution of the temporary coorbitals among the giant planets, we find that our predictions (∼ 0:1:2:2 for J:S:U:N) are consistent with the observations (0:1:1:2).