The Intrinsic Neptune Trojan Orbit Distribution: Implications for the Primordial Disk and Planet Migration
Abstract
The presentday orbit distribution of the Neptune Trojans is a powerful probe of the dynamical environment of the outer solar system during the late stages of planet migration. In this work, I conservatively debias the inclination, eccentricity, and libration amplitude distributions of the Neptune Trojans by reducing a priori unknown discovery and followup survey properties to nuisance parameters and using a likelihoodfree Bayesian rejection sampler for parameter estimation. Using this surveyagnostic approach, I confirm that the Neptune Trojans are a dynamically excited population: at 95% confidence, the Neptune Trojans’ inclination width must be . For comparison and motivation purposes, I also model the Jupiter Trojan orbit distributions in the same basis and produce new estimates of their parameters (Jupiter Trojan , , and ). The debiased inclination, libration amplitude, and eccentricity distributions of the Neptune Trojans are nominally very similar to those of the Jupiter Trojans. I use these new constraints to inform a suite of simulations of Neptune Trojan capture by an eccentric, rapidlymigrating Neptune from an initially dynamicallyhot disk. These simulations demonstrate that if migration and eccentricitydamping timescales were short ( Myr, Myr), the disk that Neptune migrated into must have been preheated (prior to Neptune’s appearance) to a width comparable to the Neptune Trojans’ extant width to produce a captured population with an inclination distribution width consistent with that of the observed population.
1. Introduction
A small sample of Neptune Trojans has been accumulated by a variety of surveys; however, inferences drawn from this sample about the intrinsic distributions of Neptune Trojan orbital properties have been limited and generally qualitative. The challenge inherent in extracting meaningful information from this sample is accurately determining the properties of the surveys that discovered them, and the properties of all surveys that, while sensitive to Neptune Trojans, did not discover any. In this work, I treat these unknown survey properties as nuisance parameters, and marginalize over them to extract as much useful information about the intrinsic orbital distributions of the Neptune Trojans as possible.
Of particular dynamical interest are the inclination, eccentricity, and libration amplitude distributions. These distributions encode information about the formation mechanism (insitu formation, chaotic capture, or other processes) and postformation evolution. Several Neptune Trojans have remarkably high inclinations (), even though surveys have byandlarge targeted fields near the Ecliptic where objects on inclined orbits spend relatively little time. Previous works have noted that this qualitatively indicates the existence of a large, poorlysampled highinclination Neptune Trojan population (Sheppard & Trujillo 2006, Sheppard & Trujillo 2010a).
In this work, I simultaneously consider the inclination, eccentricity, and libration amplitude distributions, generate synthetic populations of Neptune Trojans defined by these distributions, then pass these synthetic populations through “coverage functions:” simplified observational filters that are treated as independent functions of heliocentric ecliptic latitude , and heliocentric longitudinal separation from the Trojan libration centers (libration centers located roughly from Neptune), and inclination. The properties of these observational coverage functions are then marginalized over, effectively marginalizing the unknown properties of the surveys which discovered the Trojans. Like the Jupiter Trojans and other transNeptunian populations, the inclination distribution is modeled as a Brown’s distribution ( ). The libration amplitude and eccentricity distributions are both modeled as Rayleigh distributions, motivated by the distributions of the Jupiter Trojans. The inclination, libration amplitude, and eccentricity distributions are all truncated at upper limits derived from stability constraints, requiring appropriate corrections to their proposal volume and probability density functions. All statistical analysis is performed in the conceptually simple yet analytically powerful “Approximate Bayesian Computation” framework, described in section 3.
2. Sample
Table 1: Adopted Neptune Trojan Properties  
Name  
2001 QR  1.3  25.5  0.57  10.46 
2004 UP  1.4  10.8  0.73  10.24 
2005 TN  25.0  8.7  0.62  8.51 
2005 TO  5.2  9.2  1.62  9.12 
2006 RJ  8.2  6.3  7.99  0.58 
2007 VL  28.1  14.2  11.25  9.44 
2008 LC  27.6  16.4  2.80  0.84 
2011 HM  29.4  9.8  2.60  7.72 
: J2000 ecliptic inclination. : Halfpeak RMS libration amplitude and uncertainty. : Absolute value of J2000 heliocentric ecliptic longitude separation of object and nominal Trojan center; see Eqn. 7. : J2000 Heliocentric ecliptic latitude.
The Minor Planet Center (MPC) lists nine Neptune Trojans (six L4, three L5), but one of the L5 Trojans is unstable and likely a recently captured Centaur (Gladman et al. 2012, Horner et al. 2012). This object is therefore not considered to be reflective of the intrinsic inclination distribution of the (putatively primordial) Neptune Trojans. With the addition of the newly discovered L5 Trojan 2011 HM (Parker et al. 2013), the 8 known longterm stable Trojans have ecliptic inclinations ranging from to , and heliocentric ecliptic latitudes at discovery ranging in amplitude from to . Figure 1 illustrates these properties, and it is clear that Trojans have generally been discovered at latitudes significantly lower than their inclinations, even though an object spends roughly 50% of their time at latitudes greater than 70% of their inclination. Only the object 2006 RJ was higher than its median latitude at the time of discovery. This indicates that it is likely that most surveys that discovered Neptune Trojans targeted the ecliptic, and were therefore strongly biased toward detecting low inclination objects, and yet discovered a surplus of highinclination objects. The larger number of known L4 Neptune Trojans compared to L5 is likely an artifact of the L5 cloud being more poorly surveyed due to its current proximity to the Galactic plane.
These Neptune Trojans were discovered by a variety of surveys, performed at a variety of facilities and under varying conditions, and normally would not represent a sample from which estimating an intrinsic, debiased orbit distribution would be statistically advisable. However, using appropriate statistical care, we can make a conservative estimate of the range of plausible properties of the Neptune Trojan orbital distribution by marginalizing over the plausible volume of the unknown characteristics of all discovery surveys. This surveyagnostic approach can conceivably be applied to other populations, and since it is performed in a Bayesian framework, the outcomes can be meaningfully combined with results from large, monolithic, wellcharacterized surveys such as DES (eg., Gulbis et al. 2010), CFEPS (eg., Petit et al. 2011) and the ongoing Outer Solar System Origins Survey
The libration amplitudes listed in Table 1 were generated with the same technique as Parker et al. (2013). Each object’s motion was integrated with mercury6 (Chambers 1999) in the presence of the giant planets for 1 Myr. 100 clones of each object were integrated, with initial state vectors centered on the JPL Horizons solution, perturbed to populate the Cartesian uncertainty manifold generated by fitting all groundbased observations of each object with the fit_radec and abg_to_xyz routines developed in association with Bernstein & Kushalani (2000). Libration amplitudes for each clone were measured by assuming that libration is sinusoidal and deriving the sinusoidal halfamplitude from the RMS of the samples of the resonant angle over the entire 1 Myr integration:
(1) 
which produces the more appropriate halfamplitude for scaling a sinusoidal model than the usual peaktopeak definition of . The RMS produces a value that better reflects the mean libration behavior, while defining the amplitude from peak to peak is sensitive to large, singlecycle excursions of the resonant argument. As such, the RMSdefined amplitude is always smaller than the peaktopeak definition.
3. Approximate Bayesian Computation
For all parameter estimation in this work I utilize a likelihoodfree rejection sampler — specifically, the “Approximate Bayesian Computation” rejection (ABCr) scheme first presented in Pritchard et al. (1999). Approximate Bayesian Computation (ABC) is conceptually simple but statistically powerful, and has the primary advantage of not requiring the computation of any true likelihood value. I briefly outline this approach below, and refine its description later as merited by the specifics of each application. The literature on ABC is welldeveloped, with much more sophisticated methods available than those employed here; for a recent review see Marin et al. (2011).
In the context of this work, I utilize ABCr in order to approximate the posterior probability density functions of a set of parameters which define the properties of a model that describes an observed population of minor planets. Given a metric of similarity that reflects the relative similarity of some set of properties of an observed sample and a synthetic sample , where the synthetic set is generated by a model , and where the parameters of are drawn from prior distribution :

Propose a new set of model parameters , drawing from the prior distribution for each parameter.

Generate a synthetic sample of observations from model with parameters .

Generate the similarity metric for this synthetic sample compared to the observed sample.

If similarity metric is less than some cutoff value , keep as a successful trial.

Repeat steps 1—4 until sufficient number of trials have been successful.
The distribution of parameters of the retained sample from successful trials approximates the posterior distributions of those parameters in the Bayesian framework, given is sufficiently small and the number of trials is sufficiently large. In practice, I set dynamically by sampling many times and saving each sampled , then determining a value for which will retain less than some fraction of the total sample (in this work, this acceptance rate was set between 0.1%—0.01%, depending on the number of dimensions of the parameter space under consideration).
The similarity metric is usually defined by some distance on a set of summary statistics of the properties of and , but implementations vary on a casebycase basis. The use of summary statistics to define a metric of similarity is computationally advantageous for large datasets where other methods might be computationally prohibitive. If the summary statistics adopted are sufficient for the underlying distributions, then there is no information lost by using them. In this work, however, the data set is very small, and in the absence of a clear choice of sufficient summary statistics I instead utilize the twosample AndersonDarling test to determine a similarity metric (the twosample ADtest statistic). For multiple dimensions, the adopted metric is the sum of all onedimensional, twosample statistics over all dimensions, divided by the number of dimensions tested. Qualitatively, this is small when all distributions are similar to each other, and large when one or more of the distributions are discrepant. In this implementation, the critical values of or and samplesize effects are accounted for as the ABCr analysis effectively bootstraps these out. The adoption of a metric based on a collection of 1D 2sample ADtests is motivated by its use by the CFEPS survey (eg., Petit et al. 2011), but by utilizing the ABCr framework instead of simply adopting the lowest value inferred from each independent ADtest over all dimensions, this work reduces the risk of erroneously identifying a low value in any given application of the metric.
4. The Orbit Distributions of the Jupiter Trojans
The only other large population of Trojan objects in the solar system belongs to Jupiter. While evidently less populous than Neptune’s Trojan swarms (Sheppard & Trujillo 2006), the Jupiter Trojans are much better sampled, and are nearly complete to . To motivate the functional forms of inclination, libration amplitude, and eccentricity distributions I will adopt for the Neptune Trojans, I consider those of the Jupiter Trojans as templates.
To date, no populationwide orbit distribution models have been computed for the Jupiter Trojans. Here I consider only the relatively complete sample, and model the inclination, libration amplitude, and eccentricity distributions.
I find that the inclination distribution is very wellmodeled by a truncated Brown’s distribution (Brown 2001):
(2) 
while both the libration amplitudes are wellmodeled by truncated Rayleigh distributions defined by widths and and truncations and . For the Jupiter Trojans I hold the truncations of each distribution fixed at the maximum observed value of each parameter; , , and . In the later Neptune Trojan analysis, these thresholds are defined by dynamical constraints.
For consistency, I perform parameter estimation in the same ABCr framework that I later apply to the Neptune Trojans, though I only consider one distribution at a time. In the ABCr analysis, the parameters of interest (, , and ) are sampled uniformly in their distributions’ means. That is to say, for each trial, is selected uniformly, and then mapped to given that the model inclination distribution is a Brown’s distribution. Because the nominal functional form for the inclination, libration amplitude, and eccentricity distributions are truncated at some upper limit, their means are not trivial to compute. For example, the mean of a nontruncated Rayleigh distribution is simply given by , but for a truncated Rayleigh the mean is given by
(3) 
where is the truncation value (see Appendix A.1 for a derivation and more details). The functional form for a truncated Brown’s distribution is similar (see Appendix A.2 for a derivation and more details).
Proposing in the distribution mean has several advantages; chiefly, it translates an infinite prior volume () to a finite one, as the distribution mean is limited to a maximum value by the truncation of the distribution. For a truncated Rayleigh distribution, the mean asymptotically approaches as . A truncated Brown’s distribution has similar asymptotic behavior (see Appendices A.1 & A.2 for details).
For the Jupiter Trojans, the uniformmean prior is not particularly important because the sample is so large and welldefined. However, when I later consider the more poorlysampled Neptune Trojans, the utility of this prior becomes much more apparent.
For each distribution, ABCr trials were run, and the 1,000 best trials were retained to generate the posterior PDF for each parameter. Figure 2 illustrates the sample of observed osculating inclination, libration amplitude, and osculating eccentricity, compared to the accepted model distributions, and also illustrates the posterior PDFs for , , and defined by the accepted 1,000 ABCr trials. The observed Jupiter Trojan sample is very well modeled by the adopted functional forms for the orbit distribution; a KStest indicates peak values of 92%, 70%, and 98% for the inclinations, libration amplitudes, and eccentricities, respectively; indicating that there is no statistically significant evidence for discrepancy between the model distributions and the observed sample and qualitatively indicating good agreement between the two.
5. Technique for debiasing the Neptune Trojans
5.1. Survey coverage
Because of the lack of consistent information regarding their discovery surveys, a CanadaFrance Ecliptic Plane Survey (CFEPS)style surveysimulator approach (Jones et al. 2006) is not feasible for the Neptune Trojans. Instead, I adopt a set of basis functions informed by typical transNeptunian object survey strategies which represent the coverage of the sky by all surveys which did discover or plausibly could have discovered Neptune Trojans. The parameters of these “coverage functions” are allowed to vary, and the discovery circumstances of each known Neptune Trojan are used to quantify how well a given choice of parameters reproduces them, given a set of assumed orbit distribution parameters. Finally, the properties of these coverage functions are marginalized over, and conservative orbit distribution posterior probability density functions are determined.
To model the a priori unknown survey properties, I assume that the probability of detecting a given Trojan is essentially independent of its orbital phase save for effects of areal coverage in ecliptic latitude and longitude . I fold all the effects of survey depth (which couples with the luminosity function) and areal coverage into two nuisance “coverage functions” and . These represent the probability of any survey detecting a Neptune Trojan given the Trojan’s ecliptic latitude and longitude at the time of discovery. Nominally, these coverage functions will (1) only depend on ecliptic latitude and longitude, (2) likely be monotonic away from the centers of the Trojan clouds, and (3) have a functional range bracketed by the interval [0,1] and a domain over the entire sky.
In addition to modeling survey coverage in (), I also include the possibility of a direct bias against highinclination objects in the MPC database. The CFEPS survey found moderately higher inclination distribution widths for the 3:2 population than previous studies; because CFEPS carefully characterized their tracking biases in a selfconsistent way, they suggested that earlier samples of objects with sufficiently accurate orbits to be confidently characterized as members of the 3:2 population may have had an unacknowledged bias against high inclinations. That such biases might exist is not surprising, and the possible sources of such biases were explored in Jones et al. (2006). Because the surveys that contributed to this inclinationbiased sample were also sensitive to Neptune Trojans, it is likely that a similar inclination bias also exists in the observed sample of Neptune Trojans.
Coverage function priors from Plutino population
In order to determine the shape of the ecliptic latitude coverage function, I use the wellcharacterized 3:2 resonant population as a tracer. It is likely that similar ecliptic latitude and inclination biases apply to the 3:2 population as to the Neptune Trojan population; most surveys that discovered 3:2 objects were capable of discovering Neptune Trojans as well, as both populations share broad inclination distributions (and therefore both occupy a broad swath of ecliptic latitudes), similar onsky rates of motion, and probably share similar luminosity functions.
The intrinsic properties of the 3:2 population have been measured by the carefully calibrated CFEPS survey (Petit et al. 2011, Gladman et al. 2012). By considering the observed sample of (nonKozai) 3:2 population members (provided by Buie, private communication), and adopting the intrinsic inclination distribution measured by CFEPS, it is possible to derive the longitudinallyaveraged ecliptic coverage function that applied to the amalgam of all surveys that reported 3:2 discoveries to the Minor Planet center.
I begin with an ansatz of an appropriate functional form for the ecliptic latitude coverage function :
(4) 
This form has substantial flexibility provided by three free parameters, but the 3:2 observations constrain them to be strongly correlated. In general, its behavior for small is driven by index and halfwidth , while for large the function flattens to a constant probability floor set by . The function is truncated above , where is the highest ecliptic latitude (at discovery) of any observed 3:2 object, since beyond this point the population provides little constraint on the properties of the coverage function.
I then consider the functional form of a possible inclination bias. In this case, there is useful information in literature to inform our functional ansatz: namely, that CFEPS and other surveys all found that the form of the 3:2 population’s inclination distribution was wellmodeled by a Brown’s distribution — — but they found significantly different widths. In that case, I determine that the plausible form of the inclination bias is a Gaussian:
(5) 
so
(6) 
where . To map the bestfit CFEPS inclination width to the bestfit DES inclination width (Gulbis et al. 2010), a Gaussian inclination coverage function requires a width of . Gladman et al. (2012) notes that of 24 3:2 objects detected by CFEPS, four had inclinations over , while no DES 3:2 objects had inclinations so high; adopting a Gaussian inclination coverage function width , the probability of the observed sample retaining an otherwise observable object on an orbit is less than 40%. At (the highest 3:2 inclination observed by CFEPS), this drops to roughly 7%.
The goal of the following analysis is to utilize the 3:2 population as a “backlight” for the Neptune Trojans, illuminating the biases that likely apply to both populations. Specifically, I do this by using the 3:2 population to inform the functional form of the coverage functions and , and to derive priors for the parameters of these coverage functions. To do this, I used the sample of all objects in the MPC which are confidently identified as nonKozai 3:2 objects by the DES classification scheme (Marc Buie, private communication). The observables I model are their inclination and their ecliptic latitude at discovery . I model their inclination distribution as a Brown’s distribution with width , and adopt the functional forms for and outlined above, and use ABCr to back out the PDFs of , , , , and . Given a sample of synthetic inclinations drawn from a model , synthetic ecliptic latitudes are generated by , where is a uniform random variate over .
To define the prior on , I adopt the inclination distribution results for the Plutino population from the CFEPS analysis in Gladman et al. (2012); I approximate the posterior probability distribution for the Brown’s distribution as a triangular PDF with minimum value , modal value , and maximum value . This triangular pdf reproduces the quoted mode from Gladman et al. (2012), 2.5% of its integrated likelihood below their quoted 95% lower confidence limit and 2.5% of its integrated likelihood above their 95% upper confidence limit.
The parameter is drawn from a uniform prior over the 5th and 95th percentiles of (). The parameter is drawn from a uniform prior over . The parameter was drawn from a uniform prior over . The parameter was defined as , where was drawn from a uniform prior over — this effectively defines the prior on such that the value of is drawn from a uniform prior.
Given these priors, the model for generating and , and the observed samples of and , I perform ABCr trials with defined as the mean of and , with equal sample sizes for the observed and model samples. The 1,000 trials which produced the smallest are reserved, and the distribution of , , , , and which defined these trials defines the PDFs of each parameter.
Given the simple population model, the proposed functional forms for and work remarkably well. Figure 3 illustrates the quality with which the model inclinations and ecliptic latitudes at discovery reproduce the observed distributions. Figures 4 and 5 illustrate the marginalized PDFs of each parameter for the Plutino coverage function model that produced the excellent match illustrated in Figure 3. These PDFs (except for that of ) will become the priors for the coverage functions applied to the model Neptune Trojans in the next section.
5.2. Longitude coverage function
Using a similar procedure to determine the form of the longitude coverage function is less informative, as there may be a bias present for targeting the center of the Trojan clouds that is not present in the discovered and reported 3:2 population.
Since I am considering two Trojan clouds, but am making the implicit assumption that both clouds have identical properties (aside from some relative occupation normalization, which does not affect this model), the observable derived from longitude I consider is the absolute value of the separation from the mean Trojan centers (defined for zero libration amplitude), or
(7) 
I adopt a rollover function, normalized such that :
(8) 
This functional form also has two free parameters (a turnover width and a turnover location ). To consider coverage ranging from effectively uniform in longitude to sharply biased toward the center of the Trojan swarms, I adopt a uniform distribution as the prior on , allowing it to range from (uniform coverage) to (strongly favoring detection toward the core of the Trojan clouds). is drawn from a uniform prior over , roughly twice the range of all observed Neptune Trojan .
5.3. Neptune Trojan Orbit distribution forms
The functional forms adopted for the inclination and libration amplitude distributions are motivated by the distributions of the Jupiter Trojans and other transNeptunian populations. The intrinsic inclination distribution is taken to be a truncated Brown’s distribution (Brown 2001),
(9) 
Similarly, the intrinsic libration amplitude distribution is taken to be a truncated Rayleigh distribution, as was adopted for the distribution of libration amplitudes for the Jupiter Trojans),
(10) 
Finally, the intrinsic eccentricity distribution is also taken to be a truncated Rayleigh distribution, as was adopted for the eccentricity distribution of Jupiter Trojans.
(11) 
These distributions are clipped by regions of dynamical instability. For the inclination distribution, there are no stable orbits above (Zhou, Dvorak & Sun 2009), and in the following analysis the inclination threshold is set at . For the libration amplitude distribution, the truncation is placed at as larger libration amplitudes are unstable (Nesvorný & Dones, 2002). Eccentricities are conservatively limited to values below (Zhou, Dvorak & Sun 2010). Besides the physicallymotivated fixed truncation locations, all of these these distributions are defined with one free parameter each.
5.4. Pericenter bias
For a variety of reasons, the eccentricity distribution of the Neptune Trojans is not likely to be strongly biased by observational effects. The eccentricity of Neptune Trojans are limited by stability constraints to relatively low values, ; additionally, by definition, all Neptune Trojans have nearly identical semimajor axes. These two facts conspire to produce a relatively weak pericenter bias. I demonstrate this by considering the relative detection rate of the eccentricity extrema of the population using a simple Monte Carlo test: draw two large samples of magnitudes from a powerlaw luminosity function with slope and maximum value . One sample I preserve, and take it to represent the apparent magnitudes of objects on circular orbits . To the other sample I add a quantity where is drawn from a sample of Keplerian orbits with unit semimajor axis,uniformly distributed Mean Anomaly, and nonzero, finite eccentricity — this sample represents the apparent magnitudes of objects on eccentric orbits . I then define a detection probability function of similar form to real survey sensitivities,
(12) 
where for typical ground based surveys. Amalgams of multiple surves tend to have broader turnovers still. I select a value for and set . The potential for bias is estimated by determining the ratio of the detection rate of objects on eccentric orbits to the detection rate of objects on circular orbits.
The results of this simple analysis are illustrated in Fig. 6; for reasonable powerlaw slopes, the expected preference for eccentric orbits over circular ones never exceeds 20% over unity — and for the most eccentric known Neptune Trojan, the preference never exceeds 10%.
In addition, the Neptune Trojan luminosity function has been observed to break to very shallow slope at (Sheppard & Trujillo 2010a). If a survey reaches fainter than this break, it effectively “punches through” the entire Trojan population and is complete over all heliocentric distances occupied by them. The majority of discovered Neptune Trojans were found in surveys that reached substantially fainter than this (six from Sheppard & Trujillo 2006 & 2010b, one from Parker et al. 2013). For the purposes of the present analysis, then, I proceed assuming that the eccentricity bias in the detected sample is weak and can be neglected.
5.5. Priors
As with the Jupiter Trojans, for the main parameters of interest (, , and ) I propose uniformly in their distributions’ means. This becomes particularly important for the Neptune Trojans because of the inversion of an infinite prior volume to a finite one, allowing the consideration of very large widths in a meaningful way and without imposing an artificial cutoff to the proposal volume.
The latitude and inclination coverage function parameters are drawn randomly from the retained sample of 1,000 Plutino model trials. This preserves all correlations present in the posteriors of these parameters when informed by the observed Plutino population. The longitude coverage function’s inverse rollover width is drawn from a uniform prior over [, ], while is drawn from a uniform prior over , containing the range of all observed Neptune Trojan .
5.6. Sample Generation
Given the orbit distributions and priors described in the previous sections, synthetic samples of Neptune Trojans can be generated. The model I use to generate the synthetic properties and discovery circumstances for the Neptune Trojans is described in detail in Appendix B; in brief, it takes as input the orbit distribution model parameters , , , , , and , and the coverage function forms and parameters, and returns a sample of synthetic () that reflects both the model intrinsic distributions and the model biases simulated by the coverage functions.
6. Neptune Trojan Orbit Distribution Results
Table 2: Neptune Trojan Orbit Distribution Parameter Confidence Intervals  

Parameter  Mode  Lower 68%  Upper 68%  Lower 95%  Upper 95% 
10  9  16  7  26  
0.044  0.039  0.070  0.033  0.125  
With Inclination Bias  
21  17  46  13  91  
Without Inclination Bias  
17  14  26  11  46  
Parameter estimation was performed in the ABCr framework. For each case under consideration, ABCr trials were run, and the 10,000 best trials from each were retained to generate the marginalized posterior PDFs for each parameter under consideration. Two cases were considered: in case (1), I included the effect of the inclination coverage function informed by the Plutinos. In case (2), the inclination coverage function was ignored. This allows the effect of the putative inclination bias on the results to be quantified.
6.1. Case 1: With inclination bias
Figure 7 illustrates the allowed intrinsic distributions of Neptune Trojan orbit properties and observables, while Figure 8 illustrates the allowed apparent distributions of these same properties after the survey coverage functions are applied. The most drastic effect of these coverage functions is apparent in the distribution of latitude at discovery; while the intrinsic distribution is substantially different from the observed sample, the biased distribution matches the observed sample very well (much as it did for the 3:2 sample).
With the inclination coverage function drawn from the Plutino priors, the width of the Neptune Trojan inclination distribution is allowed to climb to very large values () within the 95% confidence interval, effectively allowing solutions approaching a uniform distribution. The lower limit on the 95% confidence interval is , and the mode of the posterior PDF (illustrated in Figure 9) is . This is admittedly a large range; but even so it is a formal definition of the uncertainty when adopting a Brown’s distribution for the inclination distribution, and can be utilized for studying the origins of the population.
While their overall distribution is reasonably wellreproduced by the models adopted here, the observed Neptune Trojan inclinations appear somewhat bimodal, perhaps suggesting structure more complex than the simple single Brown’s distribution used. Preliminary tests did not indicate any statistical evidence for a more complex twocomponent Brown’s distribution (similar to the superimposed hot and cold classical Kuiper Belt), but such a model cannot be ruled out with the current sample.
The libration amplitude distribution is somewhat better constrained at its upper end; the maximum allowed value of within the 95% confidence interval is , which is nine degrees shy of the truncation limit of the libration amplitudes, indicating a strong requirement for the presence of the gaussian component of the Rayleigh distribution. Widths as low as are allowed, and the peak of the distribution, is very similar to that of the Jupiter Trojans (). The posterior PDF for the accepted libration amplitude distribution parameter is illustrated in Figure 10.
The eccentricity distribution is a more moderate case; if a slightly lower eccentricity truncation were adopted (eg., 0.08, Nesvorny & Vokrouhlicky 2009) the distribution would be consistent with . With the adopted truncation of , the maximum distribution width allowed in the 95% confidence interval is . The lower limit is , and the peak falls at , similar again to the Jupiter Trojans (). The posterior PDF for the accepted eccentricity distribution parameter is illustrated in Figure 11.
6.2. Case 2: Without inclination bias
When the inclination coverage function is set to unity over all inclinations, and all other priors remain as before, the entire confidence interval for the inclination distribution width is shifted to lower widths. The upper limit of the 95% confidence interval becomes better constrained, at , significantly lower than the truncation inclination and indicating evidence for the presence of the gaussian component of the Brown’s distribution. The lower limit drops by a few degrees to , and the mode of the PDF (illustrated in the lower panel of Figure 9) decreases to , similar to the nominal width of the Plutinos (, Gladman et al. 2012) and the Jupiter Trojans ().
All other distributions were unaffected, and retained identical posterior PDFs.
7. Simulated Capture of Neptune Trojans
Table 3: Synthetic Captured Trojan Properties  
ID  Capture Eff.  Captured  Peak prob.  Peak prob.  
1  26  17.820  10  0.05  1  3.20  8  0.335  0.935  
2  26  17.820  10  0.1  1  2.52  63  0.544  0.972  
3  26  17.820  10  0.15  1  1.74  87  0.745  0.158  
4  26  17.820  10  0.15  1  8.80  22  0.984  0.425  
5  26  17.820  10  0.2  1  6.00  15  0.658  0.841  
6  26  17.820  15  0.15  1  1.00  25  0.739  0.837  
7  26  17.820  2  0.15  1  3.32  83  0.861  0.144  
8  26  17.820  5  0.15  1  2.52  63  0.926  0.461  
9  28  19.191  10  0.0086  1  2.67  2  —  —  —  — 
10  28  19.191  10  0.05  1  1.20  30  0.993  0.985  
11  28  19.191  10  0.10  1  4.00  1  —  —  —  — 
12  28  19.191  10  0.15  1  5.40  54  0.366  0.813  
13  28  19.191  10  0.15  1  3.20  8  0.912  0.608  
14  28  18.531  10  0.1  1  5.20  26  0.964  0.852  
15  28  18.531  10  0.2  1  2.13  16  0.273  0.884  
16  28  17.870  10  0.15  1  1.60  40  0.325  0.879  
17  28  19.191  15  0.15  1  2.00  10  0.896  0.705  
18  28  19.191  2  0.15  1  7.20  18  0.779  0.998  
19  28  19.191  5  0.15  1  7.48  56  0.940  0.599  
20  28  18.531  8  0.2  1  8.00  8  0.396  0.969  
21  28  19.191  10  0.15  1  5.60  28  0.378  0.381  
22  28  19.191  2  0.15  1  3.60  27  0.539  0.784  
23  28  19.191  5  0.15  1  7.60  19  0.619  0.839  
24  28  19.191  10  0.15  5  5.40  27  0.467  0.690  
25  28  19.191  10  0.15  1  6.67  5  0.616  0.764  
In order to determine the implications of the measured properties of the Neptune Trojan orbit distribution, I performed an exploratory suite of Neptune Trojan capture simulations. These simulations were similar to those of Lykawka et al. (2009), except that I include an initial eccentric epoch for Neptune and consider a preexcited disk. These simulations and their outcomes are detailed in the following section.
These simulations were performed with a version of mercury6 (Chambers 1999) that includes an artificial, userdefined force designed to drive semimajor axis migration and eccentricity damping (Wolff et al. 2012). The semimajor axis evolution followed an exponential form
(13) 
where is a given planet’s semimajor axis at time in the integration, and represent the planet’s initial and final semimajor axes, respectively, and is a characteristic timescale for migration.
Eccentricity damping was performed in a slightly more general way than has been implemented in previous works. I use an identical functional form as for semimajor axis migration,
(14) 
where is a given planet’s eccentricity at time in the integration, and represent the planet’s initial and final eccentricities, respectively, and is a characteristic timescale for eccentricity damping. This allows me to damp eccentricity asymptotically to some final, finite value. Because I am interested in the longterm behavior and stability of captured Trojans, I wish to replicate the final architecture of the giant planets’ orbits as closely as possible. As such, this ability to drive eccentricities to values near their current vaues (instead of to zero) is important.
The initial configuration of the giant planets was motivated by the Nice Model and by the recent work of Wolff et al. (2012) and Dawson & MurrayClay (2012), which argued that in order to preserve the cold classical Kuiper Belt while populating the hot classical Kuiper Belt, Neptune must have had a significant eccentricity () at the start of its final epoch of smooth outward migration. This initial eccentricity must have damped away relatively quickly, with a characteristic timescale less than a few times years. They argue that the characteristic timescale for semimajor axis must have been longer than the characteristic timescale for eccentricity damping.
In these simulations, Neptune is started with a moderate initial eccentricity. However, in order to avoid early scattering interactions with Uranus that would add an unwanted stochastic element to the outcome of these simulations, for many initial configurations Neptune’s eccentricity is not set as high as dictated by the limits of Dawson & MurrayClay (2012) given Neptune’s initial semimajor axis.
The disks that Neptune migrates into are composed of massless tracer particles. The disks are started with a Brown’s distribution for inclination, a Rayleigh distribution for eccentricity, and a powerlaw semimajor axis distribution with index (approximately producing a surface density distribution) defined over the range of 22 AU 30 AU. The width of the eccentricity distribution is scaled from of the width of the inclination distribution, such that . If any object is generated with a pericenter lower than Saturn’s initial semimajor axis, it is redrawn. Initial inclination widths between and were explored. The smallest disk considered was populated with 25,000 particles, while the largest was populated with 100,000 particles.
Two initial semimajor axes were considered (26 AU and 28 AU), and the initial semimajor axis of Uranus was set in each case to sweep through a desired range of Neptune:Uranus period ratios. Since Uranus and Neptune are not trapped in their 2:1 MMR today, it is unlikely that they ever passed through it. As such, due to their current proximity to the 2:1, the parameter volume from which the system could migrate through decreasing period ratios is very small, and I only consider initial configurations where the Neptune:Uranus period ratio increased (or remained constant) with time. At one extreme, Uranus is placed as close to its current location as possible without causing earlytime interactions with the initial state of Neptune. For the AU simulation, this places Uranus at its current semimajor axis. At the other extreme, Uranus was placed such that over the course of the simulation, the Neptune:Uranus period ratio did not evolve. For equal , this is trivially achieved by setting . Several intermediate cases were explored as well.
All giant planets other than Neptune experienced no eccentricity migration; their initial eccentricities were set at values representative of their current mean values. For Neptune, regardless of initial eccentricity, the eccentricity damping timescale was set to years.
The integrations were performed using mercury’s hybrid integrator, and run for ; at this point, the disk particles were downselected to only those with semimajor axes in proximity to Neptune’s (semimajor axis within one Neptune Hill radius). The integrations were then continued with only the giant planets and the downselected disk particles until 100 Myr total time had elapsed. For two small simulations, the integration was continued for 1 Gyr. After the prime integration was complete, a short 100,000 year integration was performed where each particle’s state vectors were saved very frequently. From this 100,000 year integration, the properties of the remaining objects were determined; if they were librating around the L4 or L5 points of Neptune with libration amplitudes less than and had eccentricities less than 0.12, they were identified as stable Neptune Trojans.
The number and properties of the objects identified as stable Neptune Trojans are compiled for each simulation and modeled with the same distributions used to model the observed sample in the preceding sections. Table 3 outlines the initial conditions of each simulation, the number of Trojans captured and the implied capture efficiency, and the properties of the inclination and libration amplitude distributions of the captured sample. Capture efficiencies are inferred from the structure of the primordial disk imposed; since Neptune may only sweep up Trojans from a small radial extent within the disk, these efficiencies could have strong spatial variance and simply using a different radial disk structure could produce substantially different capture efficiencies. These efficiencies (), are often sufficiently high that Trojan populations captured through the processes modeled here could be substantially reduced without being problematically low compared to the extant population size (Nesvorný & Vokrouhlický 2009).
Figure 12 illustrates the initial conditions and capture results of all 25 simulations. In general, it is clear that for most simulations, particles prefer to become captured onto Trojan orbits with inclinations similar to their primordial inclinations, and the mean primordial inclinations of particles that end up becoming Neptune Trojans generally closely match their mean inclinations as Neptune Trojans. This will be discussed further in the following section, and individual simulations will be highlighted.
7.1. Capture Results
In all the simulations run, no captured Trojan population was found to be significantly asymmetric between the L4 and L5 swarms. Figure 13 illustrates the number captured into each cloud by each simulation, and no points fall outside the 95% confidence region.
Inclinations of captured Trojan particles were found to change little for objects on initially inclined orbits. In general, most simulations resulted in the preferential capture of particles with higher initial inclinations, resulting in final widths higher than initial disk widths due in large part not to stirring or scattering of particles through the capture process, but rather due to this preference for capturing initially inclined objects. This produces a trend whereby the captured Trojan width is nearly identical to the initial width of just those particles that would go in to become Trojans; see Figure 14. This indicates that given the rates of Neptune migration explored, the planetesimal disk must have had a significant component of highinclination objects before Neptune’s arrival in order to produce the broad inclination widths that we see today; indeed, initial widths of are required to produce widths consistent with the lower limit of the 95% confidence interval for the Neptune Trojan inclination width.
In two integrations (15 and 20), the captured Trojan width is substantially higher than the initial disk width; in these cases, Neptune’s initial eccentricity was set to 0.2, and the Uranus:Neptune period ratio was allowed to sweep through a moderate range of values. Other simulations with the same initial giant planet architecture, but colder initial disks () did not produce any captured Trojans and are not listed in Table 3. However, as in most other simulations, in these two integrations the captured Trojans were not heated by the capture process; on average, they ended up with nearly identical final inclinations and initial inclinations. This indicates that the capture process did not heat the Trojans; instead, the capture process preferred preheated objects. As such, in order to produce large observed widths, a primordial disk must already be stirred to large widths to have sufficient highinclination objects to populate the high inclination Trojan orbits; if the disk is too cold, and few objects are of sufficiently high inclination to be captured, the mean Neptune Trojan capture probability plummets. In both of these simulations, the lowest initial inclination of a captured object was . Thus, while these integrations seem to suggest that very broad inclination distributions can be generated from initially modest disk widths, they still require that the initial disk have a significant population of highinclination objects in order to populate the highinclination Trojan orbits.
In several integrations with cold () primordial disks (eg., 18 and 22), several particles are visible in Figure 12 that attained captured inclinations substantially higher than their primordial values. However, the frequency of this occurrence in these simulations was low enough that the mean inclination was not substantially increased.
Longer migration timescales tended to produce larger captured inclination widths, as illustrated in Figure 15. Again, this increase in width seems to be largely due to an increasing preference for capturing inclined objects, instead of from stirring during the capture process. The increase in captured Trojan inclination width is relatively slow, increasing by only 20% over an order of magnitude change in migration timescale. Recent arguments that the migration of Neptune may have occurred on timescales much longer than those explored here (Morbidelli et al. 2014) may provide an alternative explanation. If the relatively slow increase in captured population width with increasing migration timescale changes slope beyond the regimes explored here, the long migration timescales suggested in Morbidelli et al. (2014) could produce broad captured Trojan inclination distributions from a cold initial disk. Naive extrapolation from the simulations performed here suggest that capture efficiencies in such a scenario would be low, and may require problematicallymassive initial disks if reduced by several orders; however, as previously noted, the capture efficiencies determined by the simulations presented here could be very sensitive to the assumptions of the primordial disk’s radial structure and extent. Investing in simulations of capture simulations with much longer migration timescales than those considered here would be worthwhile.
Integration Length
Clones of two integrations (3 and 12) were extended beyond the nominal 100 Myrs, and continued out to 1 Gyr (4 and 13) to explore the longterm evolution of the captured populations. In both cases, the captured population decayed by factors of a few, and the nominal inclination distribution widths decreased marginally but the change was not statistically significant. In both cases, the nominal libration amplitudes decreased; for integration 4, the libration amplitude width decrease was statistically significant, dropping from to . For integration 13, the nominal libration amplitude width decrease was smaller and not statistically significant given the captured Trojan sample sizes.
8. Summary
The Neptune Trojans have a high mean inclination; this work has demonstrated that if the population’s inclination distribution is wellmodeled by a truncated Brown’s distribution like other resonant minor planet populations, then the width of that distribution must be greater than under conservative assumptions. By simulating the capture of Neptune Trojans from a dynamicallyexcited disk by a migrating, eccentric Neptune, I have shown that the inclination widths generally do not change drastically between the initial disk and the captured Trojan population; what change there is does not come in large part from stirring during the capture process, but rather from preferential capture of inclined objects. This indicates that objects with excited orbits today were likely excited prior to Neptune’s arrival.
A mechanism that might be responsible for this excitation is not yet clear. Migration of planets driven by interactions with excited disks has not yet been explored in great detail. The evidence presented in this work indicates that if Neptune migrated quickly, the disk Neptune migrated into must have been heated prior to Neptune’s arrival; further work investigating precisely how Neptune’s migration would have proceeded when driven by such a disk is merited.
9. Acknowledgements
The computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University. I thank Matthew Holman, Rebekah Dawson, and Konstantin Batygin for their input regarding computational techniques and processes involved in planet migration and resonant capture, and Wesley Fraser for useful discussions on statistical techniques for minor planet population synthesis. Much of this manuscript was written in the 3 Little Figs cafe in Somerville, Massachusetts, and I thank the staff and owners for their coffee, lavender biscuits, and patience.
Appendix A Means of truncated distributions
a.1. Rayleigh distribution
The PDF of the Rayleigh distribution is given by
(A1) 
If this distribution is allowed to be continuous over all positive reals, then the mean is given by . However, in the case where there is a finite upper limit applied, the mean is not as straightforward. The mean of such a truncated distribution with known PDF and CDF (in the case where the PDF and CDF are defined for arguments greater than 0) is given by
(A2) 
In the case of a Rayleigh distribution, this evaluates to
Since as grows large, the truncated distribution approaches , the mean of this truncated Rayleigh distribution has the asymptotic value of as .
Eqn. A.1 is a transcendental equation, so in order to propose a as done in the text, Eqn. A.1 must be numerically solved for . In my implementation, this was accomplished with the SciPy implementation of Brent’s method.
As long as the truncation value is greater than , then the mode of the truncated Rayleigh distribution remains . The median of the truncated Rayleigh distribution is trivial to compute, but I include it for completeness here:
(A3) 
which, as approaches the limit of the untruncated distribution of .
a.2. Brown’s distribution
The Brown’s distribution (Brown 2001) is colloquially given by
(A4) 
For this work, we will assume that it applies only to the range (neglecting retrograde orbits). This means that, by nature, the Brown’s distribution does not have a simple PDF, as it is already truncated. Additionally, since integrals of the form
(A5) 
do not evaluate to a convenient form, determining the normalization of a PDF or CDF is not trivial. In order to determine the mean value, then, we must evaluate
(A6) 
To simplify this prospect, we use the Taylor expansion of to make these integrals more similar to those we computed for the Rayleigh distribution. For truncations up to , the Taylor series up to order 5 is sufficient for our purposes.
where
Given that as grows large, the truncated distribution approaches , the mean of a truncated Brown’s distribution has the asymptotic value of as . For the approximation given by Eqn. A.2 the astymptotic limit is . For truncations up to , these asymptotic values agree to within 0.05%.
As with the truncated Rayleigh distribution, this must be solved numerically for given a desired , and this was accomplished with the same numerical method described previously.
Appendix B Generating a sample of synthetic Neptune Trojan orbits
Given a proposed parametric model for orbit distributions, I generate a synthetic sample of Neptune Trojan orbits and observational circumstances in the following way. First, inclinations , libration amplitudes , and eccentricities are drawn from their respective proposed distributions (semimajor axes are assumed to be identical). Libration phases are uniformly sampled over , and the current resonant arguments are generated by assuming sinusoidal libration behavior: .
Since I am defining the distribution with respect to the centers of the L4 and L5 clouds, I set instead of , where is a small, positive angle derived from the empirical libration properties of the synthetic Neptune Trojans generated by the capture simulations in this work. This offset was found to depend on , , and (weakly) on . True majority of the variation in mean libration center was found to be wellrepresented by ; see Fig. 16. Longitudes of ascending node are uniformly sampled over .
Instead of proposing an argument of pericenter and mean anomaly, then solving Kepler’s equation, I use a smalleccentricity approximation to probabilistically determine an offset from the object’s mean longitude caused by the object’s eccentricity (the epicyclic component of motion). For eccentricities , a sample of approximate offsets between true anomaly and mean anomaly can be drawn from where is a uniform random variate drawn over , and is in Radians.
Now it is possible to construct the cartesian coordinates of the Trojans with respect to the libration center:
(B1)  
(B2)  
(B3) 
The two other observational parameters needed are the longitudinal angle between the Trojan and the mean Trojan libration center, , and the absolute value of the Heliocentric latitude . From the cartesian coordinates,
(B4)  
(B5) 
Now that , , and are available, they are passed through their respective coverage functions; these functions represent the probability that a given object will be observed given its , , or . “Synthetic observed” Trojans are randomly sampled given the product of their coverage functions. The biased, “synthetic observed” sample of () is returned to be compared directly to the real Neptune Trojan properties.
This approach implicitly assumes that libration is approximately sinusoidal, eccentricities are low, and the L4 and L5 clouds are identical.
Footnotes
 CFHT Large Program proposal: http://cfht.hawaii.edu/en/science/LP_13_16/OSSOS.pdf
References
 Bernstein, G. & Kushalani, B. 2000. AJ 120, 3323.
 Brown, M. 2001. AJ 121, 28042814.
 Chambers, J. 1999. MNRAS 203, 793799.
 Dawson, R. & MurrayClay, R. 2012. ApJ 750, 43.
 Gladman, B., and 10 coauthors. 2012. ApJ 144, 23.
 Gulbis, A., Elliot, J., Adams, E., Benecchi, S., Buie, M., Trilling, D., & Wasserman, L. 2010. AJ 140, 350369.
 Horner, J. & Lykawka, P. 2012. MNRAS in press.
 Jones, R. & 16 coauthors. 2006. Icarus 185, 508522.
 Lykawka, P., Horner, J., Jones, B., & Mukai, T. 2009. MNRAS 298, 17151729.
 Marin, J.M., Pudlo, P., Robert, C., & Ryder, R. 2011. arXiv:1101.0955v2
 Morbidelli, A., Gaspar, H., Nesvorny, D. 2014. Icarus 232, 8187.
 Nesvorný, D. & Vokrouhlický, D. 2009. AJ 137, 50035011.
 Parker, A. & 22 coauthors. 2013. AJ 145, 96.
 Petit. J.M. and 15 coauthors. 2011. AJ 142, 131.
 Pritchard, J., Seielstad, M., PerezLezaun, A., & Feldman, M. 1999. Molecular Biology and Evolution 16, 17918.
 Sheppard, S. and Trujillo, C. 2006. Science 313, 511.
 Sheppard, S. and Trujillo, C. 2010a. Science 329, 1304.
 Sheppard, S. and Trujillo, C. 2010b. ApJL 723, L233.
 Wolff, S., Dawson, R., & MurrayClay, R. 2012. ApJ 764, 171.
 Zhou, L.Y., Dvorak, R. & Sun, Y.S. 2009. MNRAS 398, 12171227.