Hsing Wen Lin
email: hsingwel@umich.edu

Ying-Tung Chen

Brett Gladman

Ruth Murray-Clay

Wing-Huen Ip

Stephen D J Gwyn

Jean-Marc Petit

Typeset using L A T E X twocolumn style in AASTeX63

Keywords:

The minor planets on orbits that are dynamically stable in Neptune's 1:1 resonance on Gyr timescales were likely em:laced by Neptune's outward migration. We explore the intrinsic libration amplitude, eccentricity, and inclination distribution of Neptune's stable Trojans, using the detections and survey efficiency of the Outer Solar System Origins Survey (OSSOS) and Pan-STARRS1. We find that the libration amplitude of the stable Neptunian Trojan population can be well modeled as a Rayleigh distribution with a libration amplitude width σ A φ of 15 • . When taken as a whole, the Neptune Trojan population can be acceptably modeled with a Rayleigh eccentricity distribution of width σ e of 0.045 and a typical sin(i) × Gaussian inclination distribution with a width σ i of 14 ± 2 • ; however, these distributions are only marginally acceptable. This is likely because, even after accounting for survey detection biases, the known large (H r < 8) and small (H r ≥ 8) Neptune Trojans appear to have markedly different eccentricities and inclinations. We propose that like the classical Kuiper belt, the stable intrinsic Neptunian Trojan population have dynamically 'hot' and dynamically 'cold' components to its eccentricity/inclination distribution, with σ e-cold ∼ 0.02/σ i-cold ∼ 6 • and σ e-hot ∼ 0.05/σ i-hot ∼ 18 • . In this scenario, the 'cold' L4 Neptunian Trojan population lacks the H r ≥ 8.0 members and has 13 +11 -6 'cold' Trojans with H r < 8.0. On the other hand, the 'hot' L4 Neptunian Trojan population has 136 +84 -75 Trojans with H r < 10 -a population 2.4 times greater than that of the L4 Jovian Trojans in the same luminosity range.

INTRODUCTION

The Neptunian Trojans (NTs) are minor planets that co-orbit with Neptune at semi-major axes ∼ 30.1 au. These objects librate in 1:1 resonance. Like Nep-tune's other n:1 resonators, the 1:1 resonances also contain symmetric and asymmetric libration islands (Beauge 1994;Malhotra 1996;Murray-Clay & Chiang 2005;Gladman et al. 2012;Pike et al. 2015;Volk et al. 2018;Chen et al. 2019). Asymmetric 1:1 librators are termed Trojans and librate around Neptune's L4 (leading) and L5 (trailing) Lagrange points; the extent of the stable region around each point depends on both orbital inclination and eccentricity. Zhou et al. (2009Zhou et al. ( , 2011) ) showed that NTs can be dynamically stable for billion years even at very high orbital inclinations (> 25 • ).

The symmetric librators have horseshoe coorbital motion, encompassing both the L4 and L5 Lagrange points, but they are generally not long-term stable (e.g. Brasser et al. 2004a; de la Fuente Marcos & de la Fuente Marcos 2012). Additionally, some known unstable NTs librate around their Lagrange point for timescales of only Myr or shorter, suggesting that they are temporarily captured and not ancient (Guan et al. 2012;Horner et al. 2012). The existence of temporary NTs is expected (Horner & Lykawka 2012). Numerical simulations show that several percent of the Centaur population can be sticking to the 1:1 co-orbital resonances of Neptune and Uranus at any given time (Alexandersen et al. 2013).

In contrast, stable Trojans, i.e. asymmetric librators with Gyr-long dynamical lifetimes, are part of a population that dates back to events early in Solar System history (Lykawka et al. 2011;Parker 2015;Gomes & Nesvorný 2016). Many numerical studies suggest that during Neptune's outward migration, the initial NT population could be captured into its current orbits from a minor planet population that was previously excited to a wide range in eccentricities (e), inclinations (i), and libration amplitudes (A), which may explain the observed population [START_REF] Nesvorný | [END_REF]Lykawka et al. 2009;Lykawka & Horner 2010;Parker 2015;Chen & Zheng 2016;Gomes & Nesvorný 2016). Furthermore, analysis of the long-term behaviour of these captured NTs revealed that they are not expected to change their e, i, A significantly over Gyr timescales. In other words, once an object is captured as a NT, it can hold some memory of its primordial e/i at the time of capture. The same is true for NTs that may have formed in-situ (Lykawka et al. 2009(Lykawka et al. , 2011)). Thus, these facts imply that the e-/i-distributions of currently known NTs can probe the conditions of the primordial population of NTs and place insightful constraints on cosmogonic models of the outer solar system.

Thirteen of the currently known 23 NTs are dynamically stable, maintaining Trojan behavior for more than 1 Gyr forward integrations (see Ta-ble 1; c.f. Lykawka et al. 2011;Parker et al. 2013;Alexandersen et al. 2016;Gerdes et al. 2016;Lin et al. 2016;Wu et al. 2019). Using stable NTs to identify the present-day NT e-/i-distributions offers constraints on these migration models and the possible origin of NTs. For example, Lykawka et al. (2009Lykawka et al. ( , 2011) ) predicted that stable NTs possess similarly wide e-distributions for captured objects and colder e-distributions for objects formed locally.

The sensitivity of Solar System surveys to the NT's e, i, and A distribution is a function of their sky coverage (Lin et al. 2016(Lin et al. , 2019)). For example, a low-inclination orbit will spend its orbital period entirely within a lowecliptic-latitude sky, which is the predominant survey coverage for general minor planet surveys. More inclined orbits spend a smaller fraction of their orbital period within low-ecliptic-latitude fields, and thus geometrically lower the likelihood of their population's detection. Additionally, only Neptune's L4 point has been well targeted by surveys, due to the overlap of the L5 point on the Galactic plane in the era of digital sky surveys. Merely three NTs are known from L5, while twenty are known from L4 (Table 1).

Past NT surveys have suggested that the current NT population has a broad inclination distribution. Sheppard & Trujillo (2006) surveyed near L4 and found a stable high-inclination NT; based on this discovery, they concluded that the population has a thickdisk distribution. Parker (2015) noted that the eight then-known stable NTs mostly had high orbital inclinations, despite their detection in a variety of low-eclipticlatitude surveys, such as Sheppard & Trujillo (2010a,b) and Parker et al. (2013). A statistical method has been applied in Parker (2015) to de-bias the observed distributions of orbital inclinations, eccentricities, and libration amplitudes. Typically, the intrinsic inclination distribution of outer minor planet populations is modeled using a sin(i) × Gaussian distribution (Brown 2001). Parker (2015) confirmed the broad distribution, finding a width of σ i > 11 • , with statistically inconclusive hints of bimodality. The broad inclination distribution would indicate a primordial NT cloud at least as thick as the Jovian Trojan cloud, the other major Trojan population in the Solar System.

Motivation in the understanding of the NT inclination distribution became apparent following the NT observations of the Pan-STARRS1 survey (PS1; Chambers et al. 2016) and the initial results of the Dark Energy Survey (Lin et al. 2019). Both PS1 and the Dark Energy Survey have extensive sky coverage with large areas at high ecliptic latitudes, making them thoroughly sensitive to NTs on highly inclined orbits. In the PS1 analysis of Lin et al. (2016), six stable L4 NTs were detected -yet only one had an inclination larger than 20 degrees. As a raw observational result, this seems to suggest a colder NT inclination distribution than that found by Parker (2015). Lin et al. (2016) found that the acceptable range of inclination widths, 7 • < σ i < 27 • , still generally agreed with Parker (2015)'s result of σ i > 11 • . On the other hand, the first-pass analysis of part of the Dark Energy Survey (DES) by Lin et al. (2019) detected five NTs in high-latitude fields, and two of the five NTs have > 30 • inclinations. After removing the observation bias of DES using a DES survey simulator (Hamilton 2019), Lin et al. (2019) found that this result indicates a high inclination widths (σ i ∼ 26 • ) NT population. A full search of that survey for Neptunian Trojans is ongoing (Bernardinelli et al. 2020).

Additional NT discoveries and non-detections from surveys with well-characterized detection biases can improve our understanding of the intrinsic NT inclination distribution. The Outer Solar System Origins Survey (OSSOS; Bannister et al. (2016Bannister et al. ( , 2018))) discovered four NTs. As OSSOS has highly quantified detection efficiencies, its survey biases in orbital parameter space can be thoroughly modeled, using a survey simulator (Lawler et al. 2018). The smaller survey of Alexandersen et al. (A16;2016) discovered an additional stable NT, 2012 UV 177 . A useful non-detection constraint comes from the sky covered by the related Canada-France Ecliptic Plane Survey and its highlatitude (HiLat) extension (CFEPS; Petit et al. 2011Petit et al. , 2017)). We refer to the combination of these three wellcharacterized surveys as OSSOS+. Their sample of five NTs, together with their quantified survey characteristics, can provide a valuable constraint on the NT orbital distribution.

In this study, we combine OSSOS+ and PS1 to assess the intrinsic eccentricity and inclination distribution of the Neptunian Trojans. Section 2 introduces the OS-SOS+ and PS1 survey coverage and NT sample, with an analysis of the dynamical stability of the four of five OSSOS+ NTs in Section 3. We use a survey simulator for both OSSOS+ and the PS1 survey to investigate the orbital distribution of the intrinsic population (Section 4), and discuss the evidence for two size-dependent components in the NT inclination distribution. We discuss possible formation mechanisms for two-component (Section 5). We conclude in Section 6.

THE STABLE NEPTUNIAN TROJAN SAMPLE

FROM THE OSSOS+ AND PS1 SURVEYS

NTs from OSSOS+

The survey coverage of the well-characterized surveys that we term OSSOS+ neatly samples both sides and the middle of the on-sky distribution of the L4 Neptunian Trojans (Fig. 1). OSSOS had eight survey blocks of sky on or within 12 • of the invariable plane, half of which were bracketing the margins of the L4 zone. A16 had two survey regions, both within the L4 zone, with one on-plane and one at higher ∼ 14 • ecliptic latitudes, intended to constrain the L4 distribution. Both surveys are deep, to 24th-25th magnitude in r; roughly consistent with H r = 10.5 at 30 au, and H r is the r-band absolute magnitude. CFEPS had fields distributed on or near the ecliptic across a wide range of longitudes, with its high ecliptic latitude (HiLat) extension sampling fields from latitudes of 12 • to 85 • . CFEPS/HiLat are shallower at 22nd-24th magnitude in g and r. Together, the deep and mostly low-latitude coverage of OS-SOS+ is ideally placed to preferentially constrain the low-inclination L4 NTs. However, as these are all surveys with well-quantified detection efficiencies, they also constrain the more rarely detected high-inclination NTs.

There are five NTs from OSSOS+, all stable (Sec. 3.2) and from Neptune's L4 Lagrange region (Table 1). OS-SOS discovered four: 2015 VV 165 , 2015 VW 165 , 2015 VX 165 , all in the on-plane 15BC block, and 2015 RW 277 , in the on-plane 15BS block. These were among 843 outer Solar System discoveries by this wide-field imaging survey in 2013-2017, using the Canada-France-Hawaii Telescope (CFHT) (Bannister et al. 2018). The NT 2012 UV 177 was discovered by A16 in their higher-latitude survey region, also with CFHT. No NTs were detected in CFEPS.

NTs from Pan-STARRS1

PS1 has two surveys in the period between May 2010 and 2014 that we consider: the shallower all-sky 3π steradians survey, which has a limiting magnitude of r PS1 -band limiting magnitude ∼ 21.5, and the slightly deeper w PS1 -band Solar System Survey with limiting magnitude ∼ 22.5 (H r ∼ 8 at 30 au). Their sky coverage is across both L4 and L5 points (Fig. 1), though like many other surveys, PS1's analysis does not provide moving-object detections in the Galactic plane, and thus is sensitive only to large-amplitude librators and not the core of the L5 region.

In 2010-2014, PS1 detected seven NTs using the PS1 Outer Solar System Pipeline, including re-detection of two that were previously known, 2001 QR 322 and 2006 RJ 103 (Table 1). The detection algorithm was described in Holman et al. (2018). Two of the seven, 2006 RJ 103 and the brightest detection, 2013 KY 18 , were detected in the 3π survey, with the other five NTs found in the Solar System Survey. We exclude the only L5 NT in this PS1 sample, 2013 KY 18 , from our population analysis (Sec. 4), as its dynamical lifetime of about a million years indicates it is only temporarily captured (Lin et al. 2016).

DYNAMICAL CHARACTERIZATION OF THE OSSOS NEPTUNIAN TROJANS

Orbits on 10 Myr timescale

Table 1 lists the best-fit barycentric orbits for the four OSSOS NTs, and for the other currently known NTs, based on all available astrometry listed at the Minor Planet Center on April 2018. The best-fit barycentric orbit and the 1-sigma level uncertainties of the orbital parameters were obtained by the orbit fit routine of Bernstein & Khushalani (2000).

The libration characteristics of all the known NTs were determined from 10 Myr integrations of the bestfit orbit and 250 clones generated from sampling the covariance matrix of the best-fit orbit. All integrations were performed using SWIFT (Levison & Duncan 1994), with the Sun and the four giant planets included as the only massive bodies in the system, and the mass of the terrestrial planets added into the Sun. We list in Table 1 the libration center (φ center ) and amplitude (A φ ) for the resonant angle φ 1:1 = λ N -λ T , where λ N and λ T are the mean longitude of Neptune and the Trojan, respectively. In order to be consistent with the results of Alexandersen et al. (2016), Lin et al. (2016) and Parker (2015), the A φ is presented as a half-peak amplitude, which is √ 2 × σ φ , where the σ φ is the standard deviation of the resonant angle φ 1:1 .

Remarkably, the five NTs from the OSSOS+ NT sample mostly have high orbital inclinations. Four have inclination 17 • , three of which were detected in OSSOS survey regions centered on the invariable plane. Only one, 2015 VW 165 , has a low inclination of ∼ 5 • . Considering the predominantly low ecliptic latitudes of the two surveys, the generally high inclinations of these NTs are surprising. We consider the implications for a wide NT inclination distribution in Section 4.

Dynamical Stability over 1 Gyr

The long-term dynamnical stability of each of the known NTs is indicated in Table 1; the stability of most of these NT has been determined in previous works. To understand the long-term dynamical stability of the four OSSOS NTs, we integrated 1000 clones of each object for (2016,2018). ¶: Also detected by PS1.

1 Gyr. The clones were sampled from the six-parameter covariance matrix of the best-fit orbit and integrated with the same set-up as for the 10 Myr computations (Sec. 3.1). We calculated the mean lifetime τ 0 (τ 0 = 1/λ, λ is the exponential decay constant) of the clones by fitting an single exponential decay function to the number of clones remaining in the 1:1 resonance, as a function of the integration time. We consider that once a clone leaves the resonance, it is lost (not remaining) regardless of a potential return into the resonance. The results are summarized in Table 2, and we consider each NT in term below.

Highly stable: 2015 VW165 and 2015 VX165

The clones of 2015 VW 165 are extremely stable: none of its thousand clones were lost during the 1 Gyr integration. The clones of 2015 VX 165 also show good stability with τ 0 ∼ 2.81 Gyr; less than a third of its thousand clones were lost during the 1 Gyr integration. The fact that most of the allowable orbit-fit parameter space for these two NTs is long-term stable supports the 2.81 a τ0 and τ1 are the mean lifetimes of two-phase exponential decay. b None of the 1000 clones were lost during the 1 Gyr integration. conclusion that they were likely captured into the NT population early in the Solar System's history. is only 0.65 Gyr. 2015 RW 277 's clones show diverse behaviour (Figure 2). About half of 2015 RW 277 's clones have τ 0 > 1 Gyr; the other half have a shorter lifetime of τ 1 ∼ 130 Myr. Their decay is only well fit by two exponential decay functions, each with a different mean lifetime, τ 0 = 1.4 Gyr, and τ 1 = 0.11 Gyr. The lifetime of 2015 RW 277 's clones are correlated with their orbital elements, so the diversity of the decay rates is likely due to this object's observational uncertainties; their dynamical stability will be worth to re-assessing if additional observations are acquired.

Meta-stability has been seen previously in the NTs: for instance, 2001 QR 322 and 2008 LC 18 showed similar behaviors and lifetime distributions (Brasser et (Guan et al. 2012;Horner & Wyn Evans 2006;Horner et al. 2012;Alexandersen et al. 2013). We do not consider it plausible that these meta-stable NTs are captured objects on unusually long Gyr timescales; Alexandersen et al. (2013) found a mean lifetime of just 78 kyr and a max-imum of 18.2 Myr in their simulations of temporary NTs. Given this, it is likely that future refinements to their orbits will reveal that they are, in-fact, long-term stable. Thus, we consider 2015 VV 165 and 2015 RW 277 as in the stable population, rather than as temporary captures.

SURVEY CONSTRAINTS ON THE INTRINSIC NEPTUNE TROJAN ORBITAL DISTRIBUTION AND POPULATION

To place constraints on the intrinsic population of NTs, we use a survey simulator to apply the surveys' quantified detection biases (to the degree available) to an NT population model. We adjust the model's parameters until, using the two-sample Anderson-Darling test (AD test, which is more sensitive to the tails of the distribution than the Kolmogorov-Smirnov test), the model produces simulated detections that match the observed number and orbital element distribution of our observed NT sample.

The OSSOS survey simulator (Bannister et al. 2016;Lawler et al. 2018) has been used in previous works to model orbital distributions and population estimates (e.g. Kavelaars et al. 2009;Petit et al. 2011;Volk et al. 2016Volk et al. , 2018;;Shankman et al. 2017). The quantified detection efficiencies have been measured for all the OS-SOS+ surveys and incorporated into the survey simulator.

For the PS1 sample, a quantified survey efficiency is not yet available. PS1 is a highly multiplexed survey, and the chance of detecting an object increases with the number of times that it could have been observed. Therefore, we adopted the observing selection function from Lin et al. (2016) for the PS1 Solar System Survey, which based on the number of visits and worked well to estimate the NT detectability. It can be summarized as the total sum of probability mass functions of binomial distribution:

f ef f (n) = n i=10 n i 0.35 i × (1 -0.35) n-i . (1)
Here, n is the total number of exposures in a specific survey region, i is the minimal number of detections required for finding an object in Lin et al. (2016), which is 10, and 0.35 is the 50% of detectability of r=22.5 object for Solar system survey, or 21.5 objects for 3π survey times the filling factor of 70%. Thus, the detectability of PS1 is a function of number of exposures. For the Solar system survey, the typical number of exposures is between 10 to 60 and vary with sky positions. For the 3π survey, because there were so many repeated visits to various areas of the sky, typically larger than 100, the equation 1 would eventually close to 1 for a magnitude 21.5 object. However, since the galactic plane complicates the detection efficiency due to a much more crowded stellar field, and the assumption of 100% efficiency presumably would not be sufficient. Therefore, We define a effective 3π survey coverage field of the sky around L4 (0 • < RA < 60 • and -30 • < Decl. < 45 • , see Figure 1).

We generated models of the NT orbital distribution with a fixed semi-major axis (a=30.1 au), a Rayleigh distribution of eccentricity, a Rayleigh distribution of the libration amplitude and a sin(i) × Gaussian distribution of inclination. Each of these distributions has just a single free parameter, to avoid the over-fitting of a small number of observations. The Rayleigh distribution is given by:

p(x; σ) ∝ xe -1 2 (x/σ) 2 , x ≥ 0. ( 2
)
Here σ is the width of the Rayleigh distribution, and we use it to be the width of the eccentricity distribution, σ e , or the width of the libration amplitude distribution, σ A φ . The sin(i) × Gaussian distribution is an Rayleigh distribution with sin(i) instead x in equation 2, which is:

p(sin(i); σ i ) ∝ sin(i)e -1 2 (sin(i)/σi) 2 , sin(i) ≥ 0. (3)
Similar models were used in Alexandersen et al. (2016) and Lin et al. (2016).

The cumulative absolute magnitude distribution function is given by a single power-law distribution:

Σ(H) ∝ 10 αH , (4)
Here H is the absolute magnitude. We tested two different absolute magnitude distributions applied to the NT models, with power-law index α = 0.9, the bright end slope of dynamical excited TNOs (Fraser et al. 2014;Lawler et al. 2018), and α = 0.65, the average slope of TNOs (Bernstein et al. 2004). We expect the real orbital/absolute magnitude distributions to be more complicated (Volk et al. 2016;Lawler et al. 2018). However, based on the results of Figure 3 and later, since there is little difference between the different models, the single parameter/slope is fine for such a small number of detections.

To constrain acceptable models for the intrinsic orbital element distribution (the libration amplitudes, the eccentricity or the inclination), we ran the survey simulations and calculated the AD probabilities of the relevant orbital element distributions for the simulated and real observations, to test which model is most consistent with the OSSOS+ and PS1 NT detections. We simulate the two surveys entirely independently to permit cross-checks. For each model, we simulated 2000 NT detections. The AD statistics were calculated by selecting 1000 sub-samples from the 2000 simulated NT detections, with each sub-sample containing the same number of NTs in the observed samples (i.e. 6 observed NTs for PS1 and 5 for OSSOS+). Finally, we calculated the AD statistics of the sub-samples' relevant orbital element distribution, versus those of the whole simulated samples, to understand the null distributions of each NT model. We estimated the rejectability of each model by comparing the AD statistics of the observed samples with the null distributions.

Libration Amplitude Width

We test the acceptable range of σ A φ for both OSSOS+ and PS1. The inclination width of the NT model for OS-SOS+ and PS1 survey simulations were set as 22 • and 6 • , respectively (see Section 4.3). The eccentricity width were set as 0.07 and 0.27, respectively (see Section 4.2). With the inclination width and eccentricity width fixed, we generated models with σ A φ set from 5 • to 40 • with 1 • steps. We widen the libration amplitude from Parker (2015)'s best suggested σ A φ = 10 • , as most of the NTs in our sample have A φ > 10 • (Table 1).

We find, as shown in Figure 3, if the libration amplitude distribution of NTs follow a Rayleigh distribution, the libration amplitude width, σ A φ , must be greater than 5 • . Moreover, PS1 and OSSOS+ have results consistent with each other. Figure 3 shows that σ A φ = 15 • has the lowest rejectability to both samples, though σ A φ = 10 • as suggested by Parker (2015) is not rejectable. We generated models with the eccentricity width (σ e ) of the Rayleigh distribution varied from 0.01 to 0.09, and constrain acceptable models in the same way as with the simulations of σ A φ (Sec. 4.1), but while varying the eccentricity. Fig 4 shows that the rejectability of various eccentricity models is startlingly different between the two surveys. The result suggests that if the eccentricity distribution follows Rayleigh distribution, the PS1 sample favors a generally smaller and narrower eccentricity distribution than the distribution of the OSSOS+ sample. Using a smaller slope (α = 0.65) of absolute magnitude distribution function does not significantly affect the acceptable range of eccentricity width.

Eccentricity Width

Inclination Width

With the distribution of all orbital elements fixed to their least rejectable values except inclination, we generated 14 different models with the inclination width (σ i ) of the sin(i) × Gaussian distribution varied from 4 to 30 • with a 2 • step. The reference plane is the invariable plane of the Solar System (Souami & Souchay 2012). We constrain acceptable models in the same way as with the simulations of σ A φ (Section 4.1) and σ e (Section 4.2), but while varying the inclination.

Fig 5 shows that the rejectability of various inclination models is startlingly different between the two surveys. The PS1 NT sample suggests that if the stable NT population has a sin(i) × Gaussian inclination distribution, the acceptable σ i range should be 16 • , while the OS-SOS+ sample requires a higher inclination distribution with σ i 12 • . Similar to the σ e test, the two different H distributions that we chose do not affect the result.

For the mutually acceptable range of σ i , we consider the overlapping region where both surveys are not rejectable at the 95% level. This intrinsic inclination distribution has a sin(i) × Gaussian distribution with σ i of 12-16 • . This result is very similar to the measured inclination widths for other resonant populations: Jovian Trojans (σ i ∼ 14 • , Parker ( 2015 However, the acceptable σ i and σ e intervals of the two surveys are very different, and there is no easily evident detection bias in PS1 against finding highinclination Trojans, considering PS1 has detected a high-inclination/retrogade object (471325) 2011 KT 19 (Chen et al. 2016). Moreover, PS1 and OSSOS+ have different limiting magnitude (H r ∼ 8 and 10.5 at 30au, respectively) it is worthwhile to consider if the data are better represented by an alternative, size-Dependent model.

A Bimodal and Size-Dependent NT Eccentricity/Inclination Distribution

Bimodal and size-dependent inclination distributions are known elsewhere in the trans-Neptunian populations.

The classical Kuiper belt has a wellknown bimodal inclination distribution (Brown 2001;Levison & Stern 2001). The inclinations of classical KBOs can be fitted by two sin(i) × Gaussian distributions with different σ i (Brown 2001;Gulbis et al. 2010). At least two ('cold' and 'hot';Elliot et al. 2005) or more ('kernel', 'stirred', and 'hot';Petit et al. 2011) populations exist in the classical Kuiper belt. The luminosity function of each is distinct (Fraser et al. 2014).

We test if the NT population could have a bimodal and size-dependent inclination distribution. In this scenario, the physically larger (H r < 8) NTs mostly occupy dynamically cold, low eccentricity and low inclination orbital phase space (i < 10 • ), and the physically smaller (H r > 8) NTs exist in dynamically hotter, high eccentricity and high inclination phase space. To test this scenario, we mixed a 'cold' component and an 'hot' component in our NT model. We found that to match the real detections of PS1 and OSSOS+, the 'cold' component needs a shallower H distribution with α ∼ 0.2, a cutoff at H r ∼ 8, a eccentricity distribution with σ e ∼ 0.02, and an inclination distribution with σ i ∼ 6 • , which is similar to the 2:1 resonant population (∼ 7 • in Gladman et al. (2012) and σ i ∼ 6 • in Chen et al. ( 2019)). We set the the cutoff at H r ∼ 8 based on the fact of only one H r < 8 NT has inclination > 10 • . Moreover, this cutoff avoids OSSOS+ detecting too many synthetic low inclination NTs.

The 'hot' component does not need any cutoff in H distribution. We set a 'divot' absolute magnitude distribution, which have a bright-end slope and transitioning to a faint-end slope with a differential number contrast (see Lawler et al. (2018) for detail description), σ e ∼ 0.05 and larger σ i ∼ 18 • , similar to the 5:1 resonant population (σ i ∼ 22 • , Pike et al. (2015)). We note that the selection of the two power laws 'divot' or single power law absolute magnitude distribution has very little effect on the results of orbital distribution. The reason we choose 'divot' here is based on the fact of that NT population lack the intermediate-sized members (Sheppard & Trujillo 2010a), and it is reasonable for the total population estimation (see section 4.5).

However, does the more complex two-component model represent the underlying NT population better than the simple single component model? We therefore evaluate the model preference using the Bayesian model comparison (Kass & Raftery 1995). To do so, we calculate the Bayes factor, K, between two models:

K = P r(D|M 1 ) P r(D|M 2 ) = P r(θ 1 |M 1 )P r(D|θ 1 , M 1 )dθ 1 P r(θ 2 |M 2 )P r(D|θ 2 , M 2 )dθ 2 = P r(M 1 |D) P r(M 2 |D)
P r(M 2 ) P r(M 1 ) .

(

) 5
Here D is observed data, M 1 and M 2 are the plausibility of the two different models, and θ 1 , θ 2 are the model parameters of M 1 and M 2 , respectively. If we propose the model priors are equal, P r(M 2 ) = P r(M 1 ), the Bayes factor is the ratio of the posterior probabilities of M 1 and M 2 .

To estimate the posterior probabilities, we use Approximate Bayesian computation (ABC) rejection algorithm. In this schema, a set of parameters θ that define the properties of model M are first sampled from a prior distribution. From this sampled parameter set θ, a data set D sim is simulated under the model M specified by θ. A similarity metric ρ(D sim , D obs ) represents the similarity between simulated data set D sim and observed data D obs . If the simulated data set D sim is too different from the observed data D obs the sampled parameters θ can be rejected. Therefore, we can set a cut off value ǫ, and the acceptance criterion in ABC rejection algorithm is:

ρ(D sim , D obs ) ≤ ǫ. (6)
We preform the ABC rejection and the Bayesian model comparison by the following steps:

1. Pick θ p (σ e -cold, σ e -hot, σ i -cold σ i -hot for two components model, and σ e , σ i for single component model) from the prior distribution (Table 3) to build Neptune Trojan model M p .

2. Run survey simulation to generate 2000 simulated detections D simp based on the above model.

3. Calculate similarity metric ρ(D simp , D obs ). Use two-sample AD statistics to determine a similarity ρ. For two dimensions (σ e and σ i ), the similarity metric ρ is the sum of the two one-dimensional AD statistics divided by two.

4. The cut off value ǫ set dynamically. Similar to the bootstrapped AD statistics performed in section 4, randomly drew 6 PS1 and 5 OSSOS+ simulated detections for hundred times to determine the null distribution. Set a dynamical ǫ at top 5% level of null distribution. Therefore, if a trial has ρ smaller than ǫ, keep θ i as a successful trial.

5. Repeat steps 1-4 until a sufficient number of trials have been successful.

6. Finally, the ratio of the acceptance rates of two models is approximately the Bayes factor.

We noted that the Bayes factors could be sensitive to the prior distribution of parameters. We use wider prior distributions on the two-component model than the single-component model (Table 3), and such prior

σ i(-hot) uniform 10 • -18 • 14 • -26 •
distributions selection is a benefit for the acceptance rate of the single-component model; nevertheless, we find a Bayes factor K ∼ 7, which is substantial evidence that the two-component model is the preferred model for modeling NT population (Kass & Raftery 1995).

Fig 6 shows the cumulative distributions of the synthetic and real detections with our nominal values of σ e φ , σ i φ and σ A φ and σ i . The libration amplitude Rayleigh distribution with a width σ A φ = 12 • agrees well with both the PS1 and OSSOS+ detections. The inclination distributions of the PS1 and OSSOS+ detections are well-matched with a bimodal sin(i) × Gaussian inclination distribution with width σ i = 6 • and 18 • , respectively. The bimodal Rayleigh eccentricity distribution with width σ e = 0.02 and 0.05 provides a fine-match with the PS1 and ossos+ detections.

Fig 7 shows the cumulative distributions of synthetic detections generated from the single-component model. While the single component model is non-rejectable, it shows the difficultly to reproduce several observational results: the expected eccentricity distribution is either too high (PS1) or too low (OSSOS+) comparing with the observations, and while the inclination distribution seems well match with OSSOS+, it is too high for the PS1. Moreover, the single component model also expects to have more faint end detections (H r > 9) from OSSOS+, which is not the case. We summarize the parameters of our NT models in Table 4.

The Population of Stable L4 Neptune Trojans

With the bimodal inclination distribution from Sec. 4.4, we use the survey simulator to estimate the intrinsic number of both 'cold' and 'hot' stable L4 NTs. This is the total number of objects in the model population necessary to match the 6 and 5 detections of PS1 and OSSOS+, respectively. The absolute magnitude distribution of 'cold' component is the same as Sec. 4.4 with α = 0.2 and a cutoff at H r = 8. For the absolute magnitude distribution of 'hot' component, we use the parameters in Lawler et al. (2018) and set N (< H) ∝ 10 α b Hr with a bright end slope α b = 0.9 from the largest observed NTs, H = 6.9, down to H = 8.3. We set a 'divot' (Shankman et al. 2013(Shankman et al. , 2016) ) at H = 8.3, and after that we have a faint-end slope α f = 0.5 to satisfy the fact of no NT detected between 8.1 <H r < 8.6. To match the number of H r < 8 L4 NTs that were detected in our surveys (6 NTs, all from PS1), we require the L4 island to have a total 23 +20 -11 NTs for H r < 8.0, with 13 +11 -6 and 10 +9 -5 NTs from a cold and hot component, respectively. To match the total number of L4 NTs that were detected in both of the surveys, which have H r < 10 (roughly equal to a diameter of 50 km for an albedo of 0.05), the total L4 NT population should be 149 +95 -81 . Because the L4 NTs have no cold members H r ≥ 8.0, 136 +84 -75 of the 149 NTs with H r < 10 belong to hot component. This result is in good agreement with the result of Lin et al. (2019) of 162 ± 73, the more approximate estimate of 250 Sheppard & Trujillo (2010a), and the upper limit of 300 L4 NTs estimated by Gladman et al. (2012).

On the other hand, if we consider the singlecomponent Neptunian Trojan model, the total popula- tion estimation of L4 NT is 180 +115 -98 , which is close to the estimation using the two-component model.

Moreover, the MPC (Minor Planet Center) lists 63 L4 Jovian Trojans with H V < 10.2 (2020 Feb. 7), which is equal to H r < 10 if assuming the V-r color of Trojan is 0.2 (Jewitt 2018;Smith et al. 2002). Thus, our result indicates that the L4 NT population is about 2.4 times greater than that of the L4 Jovian Trojans, for H r < 10.

Non-Detection of L5 NTs and Total Population

The OSSOS+ survey did not detect any L5 NTs. However, Figure 1 shows that two OSSOS survey regions and the CFEPS survey overlaps with the L5 NT region. Could the OSSOS+ non-detection of L5 Trojans mean the L5 population is smaller than the L4 population?

Assuming that the populations at L4 and L5 are equal and symmetric, the expected number of detected L5 NTs is only 8% of the total detections of simulated NTs. Thus, the non-detection of L5 NTs by OSSOS+ is expected: there is approximately a 67% chance to have this non-detection result. We conclude that the assumption of same-sized populations in the L4 and L5 camps is not rejectable by the current NT sample. Such a result consists of the similar populations in the L4 and L5 clouds suggested by Sheppard & Trujillo (2010b). Note that most of the simulated L5 NT detections belong to the cold component; without cold NT component in L5 island, the chances of a non-detection result in OSSOS+ are much higher.

If the L4 and L5 camps are symmetric, the total population estimation using two-component model will be: number of cold Trojans = 26 +22 -12 for H r < 8.0, and number of hot Trojans is 272 +168

-150 for H r < 10. By using single-component model, we estimate the total number of NT population is 360 +230 -196 for H r < 10. Both results with the estimation of Lin et al. (2019).

DISCUSSION

Since the two-component model is our preferred model for modeling the NT population, in the section, we discuss the possible scenarios for the formation of the two components.

In Section 4.1, we found that the PS1 and OS-SOS+ NTs have consistent libration amplitude distri-butions. Since the current libration amplitude distribution of NTs correlates to how they were captured (Fleming & Hamilton 2000;Lykawka et al. 2009;Gomes & Nesvorný 2016), such a result suggests that if the NT population has two components (see Section 4.4), they were likely captured at the same stage of Neptune's migration. Thus, we propose the following three possible scenarios:

1. Formed as two components: To form a 'hot' component of NTs with σ i ∼ 20 • , the NTs need to be captured from an initially widely dispersed disk (Parker 2015). On the other hand, the 'cold' component has only σ i ∼ 6 • and can be captured from a thin disk [START_REF] Nesvorný | [END_REF]Lykawka & Horner 2010;Parker 2015;Chen & Zheng 2016). If Neptune swept through a disk with both a thin and thick component, that might capture a population with two inclination distributions but one libration amplitude distribution. We know that the current classical Kuiper belt has two such overlapping populations (though it is generally thought that the hot and cold components formed at separate locations in the original disk; see, e.g. [START_REF] Morbidelli | The Trans-neptunian Solar System[END_REF]). Was such an overlapping population in place in the early, closer part of the disk?

A more likely two-component formation scenario is that, like Lykawka et al. (2009Lykawka et al. ( , 2011) ) discuss in detail, the origin and evolution of Trojans formed locally with Neptune (referred to in those papers as "pre-formed") and were captured from trans-Neptunian disks. So, there remains the possibility that a small fraction of local NTs survived to this date. Unsurprisingly, in general local NTs display low e, i (< 0.1; < 5 -10 • ), while their captured counterparts display wide ranges of e, i. Depending on the initial conditions or model details, the surviving fraction of local NTs may range from virtually zero to tens of % (Chen & Zheng (2016) found similar results). Thus, it's possible that local and captured NTs survived with similar fractions and that they may be akin to cold/hot components in the cloud. This scenario has a observable consequence: we would expect the same cold/hot NT components can be also observed in L5 cloud. This scenario might also imply that the color distributions could differ between the cold and hot components, since the corresponding origins of the two components could have different color distributions.

2. Formed as one component, evolved into two: This scenario originated in Lin et al. (2019). The NT population could form with an intermediate-width inclination distribution and then evolve into two components. Almeida et al. (2009) found that the collision rate between Trojans and Plutinos is much higher than Plutino-Plutino or Trojan-Trojan collisions, and it is more effective for the low inclination objects. Almeida et al. (2009) also suggested that this finding could explain the existence of size-and color-inclination dependencies in the Plutino population. If this is true, the same size-and color-inclination dependencies should also present in NT population. Especially since the Plutino population has a lower inclination distribution than the NT population, the high inclination NTs have higher chance to avoid collisions. This scenario may also explain why the colors of NTs differ from the color bimodality of the other trans-Neptunian object populations: the collisions remove the ultra-red matter (Gil-Hutton 2002;Grundy 2009) from the surfaces of NTs. However, would the Trojan-Plutino collision rate be high enough? Can it eliminate the small-sized NTs to produce the cutoff after H r > 8? This is also questionable, and beyond the scope of this paper.

3. A collisional NT family in low eccentricity/inclination orbital space: The only known collisional family in Kuiper Belt, the Haumea family, has a shallow H-distribution slope and lacks small family members (Pike et al. 2020). These facts suggest that the Haumea family formed in a graze-and-merge scenario rather than a catastrophic collision. Similar to the Haumea family, the cold component of NTs also has a shallow H-distribution slope and lacks members smaller than H r > 8. Moreover, all of the cold component candidates (inclination < 10 • and H r < 8) have the same color (Jewitt 2018). Could the cold component NTs belong to a collisional NT family in low eccentricity/inclination orbital space? Unlike scenario 1, there would likely be an asymmetry between the L4 and L5 NT populations, since there may not exist another collisional family in the L5 population. Such a consequence can be tested by future L5 NT surveys.

SUMMARY

We present the orbital properties of four newly discovered NTs by the near-ecliptic survey OSSOS. Three of them have orbital inclination 17 • , as expected for the dynamically hot inclination distribution of the NT population. Our numerical integrations for the four new OS-SOS NTs showed that they are long-term dynamically stable in the 1:1 resonance, with two showing metastability within their orbit fit uncertainty ranges.

We explored the intrinsic libration amplitude, eccentricity, and inclination distributions of the stable NT population, using both the OSSOS+ surveys and the Pan-STARRS1 survey, via a survey simulator. Combined with an NT found earlier by Alexandersen et al. (2016), there are five NTs discovered by OSSOS+, and six stable NTs from the PS1 survey. The libration amplitude distribution can be described as a Rayleigh distribution with libration amplitude width σ A φ > 5 • . The best matched σ A φ is 15 • for both PS1 and OSSOS+. Using a Rayleigh eccentricity distribution model, the acceptable eccentricity width (σ e ) for both surveys is ∼ 0.045. For a sin(i) × Gaussian inclination distribution model, the acceptable inclination width (σ i ) for both surveys is 12-16 • .

Considering the two surveys have very different magnitude limits and latitude coverage: OSSOS+ is much deeper than PS1 and focuses on the ecliptic. The overlapping acceptable region for the eccentricity and inclination distributions derived from each survey are small and near the rejectable level, so we also consider an alternative scenario.

We propose size-dependent and bimodal eccentricity and inclination distributions for the stable NT population to explain the detections of NTs in the surveys we considered. One group, dynamically 'cold' NTs, has a shallow H distribution with slope ∼ 0.2, and only contain large NTs (H r < 8) on low eccentricity and low inclination orbits with σ e ∼ 0.02 and σ i ∼ 6 • , respectively. The other group, dynamically 'hot' NTs, have a wider range of sizes and eccentricity and inclination width of σ e ∼ 0.05 and σ i ∼ 18 • . We perform a Bayesian model comparison to find the preferred model between this more complex two-component model and the simpler single-component model. The result shows substantial evidence that the two-component model is the preferred model to describe Neptunian Trojan population.

With the two-component NT population model, we found that there are 13 +11 -6 'cold' NTs with H r < 8.0, and 136 +84 -75 'hot' NTs with H r < 10.0 in the L4 island. On the other hand, the population of L4 NTs is 180 +115 -98 if we use the single-component model. This result suggests that the NT population is 2.5 to 3 times larger than that of the Jovian Trojans within the same size range.

Although OSSOS has completed its observing, PS1 and now PS2 continue surveying and may discover more bright NTs. The Dark Energy Survey with CTIO has detected NTs (Gerdes et al. 2016;Lin et al. 2019;Bernardinelli et al. 2020;Khain et al. 2020), including those that are too small for PS1 to detect. Future faint NT detections will be tremendously enhanced by the discoveries of the Vera C. Rubin Observatory (LSST), particularly with the proposed North Ecliptic Spur survey (Olsen et al. 2018;Schwamb et al. 2018). These new bright-and faint-end NT samples will test the sizedependent bimodal NT eccentricity/inclination distribution that we propose.

Figure 1 .

 1 Figure 1. Sky coverage of the surveys used in this analysis relative to the Neptunian Trojan population model developed in Sec 4 (inclination width σi = 15 • , epoch = 2015 November 1 (grey dots) and epoch = 2010 November 1 (yellow dots) ). The L4 (left) and L5 (right) clumps are both shown relative to the ecliptic (thick orange line) and galactic (blue line) planes. Note that CFEPS ran most of a decade earlier than OSSOS, so they sample different parts of a moving population, despite the footprints partly overlapping. The PS1 coverage is approximately indicated in two ways: the Solar System survey to r ∼ 22.5 encompasses the region within ±12 • of the ecliptic plane (between two thin orange lines), while the 3π survey to r ∼ 21.5 covered all the area north of Decl. -30 • . The pale blue shading area shows the effective coverage (avoid galactic plane) of PS1 3π survey on the L4 clump.

  Figure2. The retention over a Gyr in the 1:1 resonance of a thousand clones of the Neptune Trojan 2015 RW277's orbit, sampled from within its covariance matrix of orbital uncertainties. The best-fit two-phase exponential decay curve shows in red dashed line, which is the superposition of two exponential decay curves, the orange dash-dotted line with mean lifetimes τ0 = 1.4 Gyr and the green dotted line with τ1 = 0.11 Gyr. This metastable behaviour has been seen previously in the NTs, by 2001 QR322.

Figure 3 .

 3 Figure 3. The bootstrapped AD rejectability of various values of the libration amplitude width σA φ . Values below the red line are rejectable at greater than 95% confidence.

Figure 5 .

 5 Figure 5. The bootstrapped AD rejectability of various values of the inclination width σi of the Neptune Trojan orbital distribution. Values below the red line are rejectable at greater than 95% confidence. The non-rejectable overlap area between the OSSOS+ and PS1 constraints on the inclination width is 12-16 • . Varying the size distribution α makes little difference to the results.

  )), plutinos (σ i ∼ 12 • in Volk et al. (2016) and ∼ 14 • in Alexandersen et al. (2016)), 5:2 resonators (σ i ∼ 11 • in Volk et al. (2016) and 14 • in Gladman et al. (2012)) and the 5:3 resonators (σ i ∼ 11 • , Gladman et al. (2012)).

Figure 6 .

 6 Figure 6. The cumulative distributions of the eccentricity (top), inclination (2nd row), libration amplitude σA φ (3rd row), and absolute magnitude Hr (bottom) for the PS1 (left) and OSSOS+ (right) NT detections. The solid lines are the distribution of synthetic detections generated by the survey simulator. The dash lines are the intrinsic distributions of NT model. The dots show the real defections. Since the eccentricity and inclination distributions are absolute magnitude dependence, the top and 2nd row left panels show the intrinsic distributions of the NT model with H ≤ 8.

  i = 60 • . b (α b , α f , c, H b ), α b : bright end power law index, α f : faint end power law index, c: contrast between bright and faint end, H b : break point. See Lawler et al. (2018) for details.

Figure 7 .

 7 Figure 7. Same as Figure 6, but modeling by a single component model.

Table 1 .

 1 Barycentric orbits and resonant dynamics of the known Neptune Trojans

	Designation	a (au)	e	i ( • )	φcenter ( • )	A φ ( • )	Hr	L	Stability Discovery Survey
	2001 QR322	30.232 ± 0.001	0.02849 ± 0.00002 1.323 ± < 0.001	68.0 +0.1 -0.1	26.4 +0.2 -0.2	7.7	L4 M *	Deep Ecliptic Survey ¶
	(385571) Otrera	30.144 ± 0.003	0.0270 ± 0.0002	1.431 ± < 0.001	61.3 +0.2 -0.1	10.8 +0.8 -0.7	8.6	L4 S *	ST06
	2005 TN53	30.125 ± 0.002	0.06584 ± 0.00003 25.001 ± < 0.001 59.16 +0.09 -0.09	8.1 +0.6 -0.6	8.8	L4 S *	ST06
	(385695) Clete	30.132 ± 0.001	0.05280 ± 0.00005 5.252 ± < 0.001	61.10 +0.08 -0.06	8.5 +0.4 -0.3	8.1	L4 S *	ST06
	2006 RJ103	30.0393 ± 0.0009 0.03002 ± 0.00002 8.163 ± < 0.001	60.380 +0.009 -0.009	5.01 +0.07 -0.07	7.3	L4 S *	SDSS ¶
	(527604) 2007 VL305	30.027 ± 0.002	0.06331 ± 0.00003 28.117 ± < 0.001 61.00 +0.05 -0.05	14.3 +0.2 -0.2	7.7	L4 S *	SDSS
	2008 LC18	29.957 ± 0.003	0.08308 ± 0.00006 27.539 ± < 0.001 297.6 +0.4 -0.4	17.2 +1.2 -1.2	8.0	L5 M *	ST10
	2010 TS191	30.008 ± 0.001	0.04586 ± 0.00003 6.563 ± < 0.001	61.83 +0.02 -0.02	11.38 +0.07 -0.07	7.8	L4 M *	PS1
	2010 TT191	30.087 ± 0.004	0.07029 ± 0.00008 4.276 ± < 0.001	65.3 +0.1 -0.1	19.5 +0.3 -0.3	7.6	L4 M *	PS1
	2011 MH102	30.111 ± 0.005	0.0806 ± 0.0005	29.377 ± 0.001	299.7 +0.4 -0.4	9.9 +2.0 -1.8	7.9	L5 S *	P13
	(530664) 2011 SO277	30.162 ± 0.002	0.01187 ± 0.00003 9.639 ± < 0.001	63.9 +0.2 -0.2	18.9 +0.5 -0.4	7.4	L4 S *	PS1
	(530930) 2011 WG157	30.030 ± 0.003	0.02791 ± 0.00007 22.299 ± < 0.001 61.38 +0.02 -0.01	15.69 +0.05 -0.04	6.9	L4 S *	PS1 (3π)
	2012 UD185	30.201 ± 0.002	0.04406 ± 0.00004 28.299 ± < 0.001 ?	?	7.4	L4 ?	PS1 (IfA)
	2012 UV177	30.024 ± 0.004	0.0723 ± 0.0008	20.833 ± < 0.001 60.7 +0.2 -0.2	9.8 +0.9 -0.7	9.0	L4 S +, *	A16
	2013 KY18	30.149 ± 0.005	0.123 ± 0.002	6.658 ± < 0.001	293.4 +1.3 -2.0	20.6 +4.1 -2.6	6.6	L5 T *	PS1
	2013 VX30	30.0876 ± 0.0006 0.08374 ± 0.00002 31.259 ± < 0.001 59.09 +0.02 -0.02	5.1 +0.2 -0.2	8.1	L4 S ♣	Dark Energy Survey
	2014 QO441	30.1019 ± 0.0007 0.10528 ± 0.00004 18.831 ± < 0.001 61.71 +0.05 -0.04	10.3 +0.2 -0.2	8.1	L4 S *	Dark Energy Survey
	2014 QP441	30.074 ± 0.003	0.0650 ± 0.0008	19.403 ± < 0.001 59.46 +0.07 -0.04	2.0 +0.8 -0.5	9.1	L4 S *	Dark Energy Survey
	2014 UU240	30.057 ± 0.001	0.0484 ± 0.0001	35.744 ± < 0.001 57.1 +0.3 -0.1	2.6 +0.3 -0.4	8.1	L4 S ♣	Dark Energy Survey
	2015 RW277	30.013 ± 0.005	0.077 ± 0.001	30.826 ± 0.001	65.8 +0.9 -1.1	27.8 +2.6 -4.1	9.92 L4 M †	OSSOS
	2015 VV165	30.120 ± 0.001	0.0850 ± 0.0001	16.855 ± < 0.001 64.5 +0.9 -0.6	20.5 +2.0 -1.3	8.83 L4 M †	OSSOS
	2015 VW165	30.102 ± 0.001	0.0511 ± 0.0002	4.998 ± < 0.001	63.3 +0.6 -0.4	16.3 +1.7 -1.0	8.03 L4 S †	OSSOS
	2015 VX165	30.073 ± 0.001	0.07522 ± 0.00003 17.140 ± 0.001	61.63 +0.02 -0.06	12.6 +0.1 -0.2	8.85 L4 S

† OSSOS Notes -a, e and i were fitted from all available observations, with an epoch of the first day of observation. φcenter and A φ were found by 10 Myr integration (Sec. 3.1). Hr of the non-OSSOS NTs were computed from the Minor Planet Center's reported HV , assuming that V-r = 0.2. L -Lagrange point of 1:1 resonance. Stability -S: dynamically stable for ≥ 1 Gyr; M: metastable for most of a Gyr; T: transiently in 1:1 resonance with lifetime on the order of Myr. Reference: † This work . ♣ Lin et al. (2019). + Alexandersen et al. (2016). * Wu et al. (2019). Discovery surveys -Deep Ecliptic Survey: Elliot et al. (2005). ST06: Sheppard & Trujillo (2006). ST10: Sheppard & Trujillo (2010b). P13: Parker et al. (2013). PS1: Lin et al. (2016). PS1 (IfA): Only recently listed at the MPC in MPS 990734, MPS 990735. A16: Alexandersen et al. (2016). Dark Energy Survey: Gerdes et al. (2016); Lin et al. (2019). OSSOS: Bannister et al.

Table 2 .

 2 Dynamical lifetimes of the OSSOS Neptunian Trojans

	Trojans	τ0 (Gyr) a τ1 (Gyr) a
	2015 RW277	1.4	0.11
	2015 VV165	0.65	-
	2015 VW165	> 4.5 b	-
	2015 VX165

  Lin et al.

		0.0 0.2				PS1, α = 0.9 PS1, α = 0.65 OSSOS+, α = 0.9 OSSOS+, α = 0.65
	Rejectability	0.4 0.6
		0.8
		0.01 1.0	0.02	0.03	0.04 eccentricity Width σe 0.05 0.06	0.07	0.08
	Figure 4. The bootstrapped AD rejectability of various
	values of the eccentricity width σe of the Neptune Trojan
	orbital distribution. Values below the red line are rejectable
	at greater than 95% confidence. The non-rejectable overlap
	area between the OSSOS+ and PS1 constraints on the ec-
	centricity width is about 0.045. Varying the size distribution
	α makes little difference to the results.

Table 3 .

 3 Prior distributionsParameter Distributions single-component two-component

	σe-cold	uniform	-	0.015 -0.045
	σ e(-hot)	uniform	0.04 -0.06	0.045 -0.07
	σi-cold	uniform	-	4 • -14 •

Table 4 .

 4 Neptunian Trojan Population Models

			single-component model	two-component model
	Parameter Distribution		cold	hot
	a (au)	constant	30.1	30.1	30.1
	e	Rayleigh	σe = 0.045	σe = 0.02	σe = 0.05
	i ( •

ACKNOWLEDGMENTS

The authors thank Matthew Holman, who kindly let us use the pointing records of the PS1 survey in this study. The authors also acknowledge the anonymous referees' useful suggestions for improving the manuscript.

HWL acknowledges the support of NASA Grant: NNX17AF21G. KV acknowledges support from NASA (grants NNX15AH59G and 80NSSC19K0785) and NSF (grant AST-1824869). M.T.B. appreciates support from UK STFC grant ST/P0003094/1 and travel support provided by STFC for UK participation in LSST through grant ST/N002512/1. WHI acknowledges the support from MOST Grant: MOST 104-2119-008-024 (TANGO II), MOE under the Aim for Top University Program NCU, Macau Technical Fund: 017/2014/A1 and 039/2013/A2. This research was supported by funding from the National Research Council of Canada and the National Science and Engineering Research Council of Canada.

Based on observations obtained with MegaPrime/MegaCam, a joint project of the Canada-France-Hawaii Telescope (CFHT) and CEA/DAPNIA, at CFHT which is operated by the National Research Council (NRC) of Canada, the Institute National des Sciences de l'universe of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii, with this project receiving additional access due to contributions from the Institute of Astronomy and Astrophysics, Academia Sinica, National Tsing Hua University, and Ministry of Science and Technology, Taiwan.

The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, and the University of Maryland.

Facilities: CFHT (MegaPrime), PS1

Software: orbfit (Bernstein & Khushalani 2000), SWIFT (Levison & Duncan 1994), OSSOS Survey Simulator (Bannister et al. 2016;Lawler et al. 2018)