Modeling the Nearly Isotropic Comet Population in Anticipation of Lsst Observations
Abstract
We run simulations to determine the expected distribution of orbital elements of nearly isotropic comets (NICs) in the outer solar system, assuming that these comets originate in the Oort Cloud at thousands of AU and are perturbed into the planetary region by the Galactic tide. We show that the Large Synoptic Survey Telescope (LSST) should detect and characterize the orbits of hundreds to thousands of NICs with perihelion distance outside 5 AU. Observing NICs in the outer solar system is our only way of directly detecting comets from the inner Oort Cloud, as these comets are dynamically excluded from the inner solar system by the giant planets. Thus the distribution of orbital elements constrains the spatial distribution of comets in the Oort cloud and the environment in which the solar system formed. Additionally, comet orbits can be characterized more precisely when they are seen far from the Sun as they have not been affected by nongravitational forces.
1 Introduction
The distribution of cometary orbital elements in the Oort cloud depends on the dynamical evolution of the solar system’s planetesimal disk and the environment in which the solar system formed. Unfortunately, the vast majority of Oort cloud comets are unobservable, usually being seen only when they are perturbed onto orbits with perihelion AU. Furthermore, the orbits of visible comets have generally been modified by poorly understood nongravitational forces (Yeomans et al. 2004). For both reasons it is difficult to infer the properties of the Oort cloud from the statistical distribution of comet orbits.
This paper describes the expected distribution of orbital elements of nearly isotropic comets (NICs). We define these to be comets that have been perturbed into the planetary region from the Oort cloud. Theoretical models of the distribution of NICs have been constructed by others e.g., Wiegert & Tremaine (1999), Levison et al. (2001) and Fouchard et al. (2013, 2014a, 2014b); the novel feature of the present study is that we focus on much larger heliocentric distances (up to 45 AU) in anticipation of deep wideangle sky surveys that are currently under development. We use a simple physical model that assumes a static spherical Oort cloud with comets uniformly distributed on a surface of constant energy in phase space. Models of the formation of the Oort cloud (e.g., Dones et al. 2004) suggest that this approximation is reasonable except perhaps at the smallest semimajor axis we examine, AU, where the cloud is somewhat flattened. Furthermore, we assume that these comets evolve solely under the influence of the Galactic tide and perturbations from the giant planets. We do not consider the stochastic effects of passing stars, as they have effects qualitatively similar to the effects of the Galactic tide (Heisler et al. 1987; Collins & Sari 2010); see Fouchard et al. (2013, 2014a, 2014b) for a detailed comparison of the influence of these two agents on the Oort Cloud. We follow the evolution of these comets through Nbody simulations for up to 4.5 Gyr and use the results of these orbit integrations to construct a simulated comet catalog.
This topic is of special interest now as the Large Synoptic Survey Telescope (LSST) will likely see many of these NICs in the outer solar system. LSST will survey 20,000 square degrees of sky (48% of the sphere) about 2,000 times over 10 years, down to an band magnitude of 24.5 (LSST Science Collaboration et al. 2009). LSST has a flux limit 3.2 magnitudes fainter than, and more than three times the area of, the current leader in finding faint distant solar system objects — the Palomar Distant Solar System Survey (Schwamb et al. 2010). It is expected to find tens of thousands of transNeptunian objects (LSST Science Collaboration et al. 2009); however we are not aware of predictions made specifically for objects originating in the Oort cloud.
Francis (2005) studied the longperiod ( years) comet population using the Lincoln NearEarth Asteroid Research (LINEAR) survey (Stokes et al. 2000). Most observed longperiod comets likely originated in the Oort cloud. He found a sample of 51 longperiod comets which were either discovered by LINEAR or would have been, had they not previously been discovered by another group. Fifteen of these had perihelion distances beyond 5 AU, but none beyond 10 AU. He used this sample to estimate properties of the Oort cloud. He found a“suggestive” discrepancy between the distribution of cometary perihelion distances in the observed sample and in theoretical models (Tsujii 1992; Wiegert & Tremaine 1999), but cautioned that the difference could be the result of a poor understanding of the rate at which comets brighten as they approach the Sun due to cometary activity. LSST will address this question by observing many comets at large heliocentric radii where they are inactive (see Section 2.2).
Hills (1981) proposed that the apparent inner edge of the Oort cloud at around 10,000 AU is not due to a lack of comets at smaller semimajor axes, but rather because the perihelion distances of those comets evolve slowly, so they are ejected or evolve to even smaller semimajor axes due to perturbations from the outer planets before they become visible from Earth. In contrast, comets with semimajor axis 10,000 AU have their perihelion distance changed by more than the radius of Saturn’s orbit in one orbital period, so they are able to jump across the dynamical barrier of the outer planets, and be seen in the inner solar system (Hills 1981). This barrier is not 100% leakproof, but as is demonstrated later in the paper, one expects the number density of comets with initial semimajor axes of 10,000 AU to decline by over two orders of magnitude interior to 10 AU. LSST should detect NICs at distances 15 AU and so will enable us to estimate the population of this inner Oort cloud directly, because we will be able to see NICs outside the region of phase space from which they are excluded by the giant planets. The properties of this cloud may contain information about the density and mass distribution in the Sun’s birth cluster (Brasser et al. 2012).
Observing NICs far from the Sun also probes in unique ways the parts of the Oort cloud that do send comets near Earth. Nongravitational forces due to outgassing when the comet comes near the Sun are the primary source of error in determining the original orbits of these comets (Yeomans et al. 2004). It is somewhat uncertain at what radius outgassing begins, but a reasonable estimate would be around 5 AU (see discussion in Section 2.2). Therefore, astrometric observations of comets beyond AU should allow much more precise determination of their original orbits (see discussion at the end of Section 2.1).
2 Simulation Description
We divide phase space into three regions, based on the perihelion distance of the cometary orbit. We define the “visibility region” as containing orbits with perihelion distance in the range AU. A “buffer region” includes orbits with . All other orbits are defined to be in the “Oort region”.
We simulated orbits with the Rebound package, developed by Rein & Liu (2012). We used their IAS15 integrator, a 15th order adaptivetimestep integrator that is sufficiently accurate to be symplectic to double precision (Rein & Spiegel 2015).
Our goal is to model the steadystate distribution of NICs with perihelia within 45 AU of the Sun that are produced from an Oort cloud with orbital elements uniformly distributed on an energy surface in phase space (so , where and are the cometary inclination and eccentricity). This approximation assumes that perturbations from the Galactic tide, passing stars or molecular clouds over the last four Gyr have been sufficient to isotropize comets both in position space (seen from the Sun) and velocity space at a fixed position. This has been shown to be roughly true for comets with semimajor axes greater than 2,000 AU (Duncan et al. 1987).
To initialize the simulation, we generated comets at random from this phasespace distribution for four discrete values of initial semimajor axis, 5,000, 10,000, 20,000, and 50,000 AU, with perihelion distances in the range 60 AU to 60 + AU. Then, using an analytic approximation to the torque from the Galactic tide (see appendix A), we determined an upper bound on the time (as a function of ) such that no comet from outside (60 + ) AU could enter the buffer or visibility regions within the next years. We chose , but it is straightforward to see that the results do not depend on this choice. We picked and years, for = 5,000, 10,000, 20,000 and 50,000 AU respectively. These numbers satisfy our upper bound.
We then evolved the comets under the influence of the Sun, the Galactic tide, and the four outer planets for . After had elapsed, we removed any comet with perihelion greater than 60 AU from the simulation. The remaining comets were allowed to evolve under the influence of the Galactic tide and gravity from the four giant planets and the Sun. At fixed intervals (taken to be 10 years), we recorded the position and velocity of any comet that was within 45 AU of the Sun in a catalog. This procedure gives us the same expected comet count and distribution of orbital elements as if we had allowed the system to evolve to steady state, and then catalogued the comets visible within 45 AU at an instant in time, and multiplied the flux by .
Comets are removed from the simulation if they collide with a planet, come within 0.1 AU of the Sun, move outside 200,000 AU, or are perturbed back into the Oort region ( 60 AU)
2.1 Orbital Elements
The treatment of orbital elements for highly eccentric orbits that pass through the orbits of massive planets is somewhat subtle. A comet that is having a close encounter with one of the giant planets will undergo large shortterm perturbations to its orbital elements that do not reflect changes to its orbit that will last longer than the duration of the encounter. Shortterm perturbations from such encounters are more serious for comets with large semimajor axes because the energy of the comet in the planetary potential well can equal or exceed the total orbital binding energy of the comet. At distance from a planet with mass , the fractional change in energy due to the potential energy of the planet is
(1) 
where is the semimajor axis of the comet prior to the close encounter. We stress that and are not the same quantity: is the current semimajor axis of the comet, whereas is the semimajor axis of the comet when the simulation was initialized.
To lessen the shortterm planetary perturbations to cometary orbital elements, we adopt the following procedure. For comets with AU, we simply report heliocentric orbital elements. These comets have large enough binding energy that they would have to pass close to a planet (within 2.0 AU for Jupiter) to obtain enough extra kinetic energy for to vary by more than 10% during the close passage (see Equation (1)). In order to prevent very close planetary approaches from contaminating our results, we discarded any observation in which the comet is currently close enough to a planet that the specific potential energy due to the planet is more than 10% of the specific binding energy in an orbit around the Sun with semimajor axis of 100 AU. This occurs for only of all catalog entries, or of all catalog entries with AU.
Comets in the visibility region with AU often have potential energies due to the planets which are comparable to their binding energies. For this reason, we report the barycentric orbital elements of the comet the last time it was in the range [90 AU, 110 AU]. These elements are wellbehaved, since they are calculated far outside the orbits of the giant planets where it is appropriate to represent the solar system as having all its mass located at the barycenter.
The ease with which LSST can determine orbital elements for slowly moving nearly unbound objects is also of interest to this study. To address this question, we searched the JPL small body database
2.2 Disrupted Comets
There is a substantial body of evidence suggesting that comets “fade” over time (e.g., Fernández 1981; Wiegert & Tremaine 1999). A number of processes have been proposed to explain this phenomenon (Weissman 1980): a comet could run out of volatile material, it could have its surface covered with a crust that prevents volatiles from escaping, or it could physically be broken apart by outgassing or tidal stress.
Fernández (2005) gives 3 AU as a likely cutoff to comet activity based on the sublimation temperature of water, but cautions that many comets experience some activity outside 3 AU due to the sublimation of more volatile elements. Comet 67P/ChuryumovGerasimenko first showed signs of activity when it was 4.3 AU from the Sun (Snodgrass et al. 2013). Comet HaleBopp showed substantial activity on its approach to perihelion as far out as 7.2 AU (Weaver et al. 1997), and at 27.2 AU after perihelion passage (Kramer et al. 2014). When we calculate numbers of visible NICs, we restrict ourselves to NICs further than 5 AU from the Sun. For this reason, we do not consider the effect of comet activity on magnitude, and just calculate the magnitude from the size and albedo of the bare nucleus, see Section 5.2. We believe that our assumption that there is negligible activity beyond 5 AU is reasonable, though not certain, given existing observations. In any event this assumption gives us a conservative estimate of the number of comets that a survey like LSST will discover.
Because comet activity does not affect brightness in our model, we are only sensitive to physical disruption of comets, not loss of volatiles. For this reason, throughout this paper, we refer to “disruption”, rather than “fading”.
In the results that follow, we remove a comet after it has made 10 apparitions in the catalog with radius AU (corresponding to a total exposure to the Sun at AU of about 100 years, since years). Comets are also removed if they ever travel within 0.1 AU of the Sun (even if they do not appear so close in the catalogue) or if they suffer a collision with one of the gas giants.
3 Comet Lifetimes in the Simulation Region
Figure 1 shows the fraction of NICs in our simulations with AU (period 200 years) appearing in the region with AU for more than years, as a function of . Different curves correspond to different values of the initial semi major axis . In this plot we terminate each orbit integration after 4.5 Gyr. The error bars are derived from the resampling procedure described in Section 4.1.2. In this and all subsequent plots, only comets with periods greater than 200 years (corresponding to semimajor axes greater than 34.2 AU) are counted.
Yabushita (1979) argued using a simple random walk model that the number of NICs surviving more than perihelion passages should scale as . Assuming for the sake of argument that NICs spend a fixed amount of time in the visibility region per perihelion passage, then the number of NICs having a given number of catalog entries is proportional to the number of NICs surviving for more than a given number of perihelion passages. We should therefore recover the same powerlaw as Yabushita (1979) if his model is a good approximation to the full physics captured by the simulation. This plot largely confirms the predictions of Yabushita (1979), but the exponent of the powerlaw seems to be slightly steeper than his value of .
Deviations from powerlaw behavior at short times occur because the visibility region is larger than the region of influence of the planets, so there is some delay before NICs that have entered the visibility region interact with the planets. Comets with smaller values of experience less torque due to the Galactic tide, so the delay is larger. They also have more binding energy that must be overcome prior to ejection. This explains the trend seen in Figure 1 that NICs with smaller take longer to be ejected.
The fact that some of our simulated particles survive for longer than 1 Gyr leads to concern about the physical validity of our assumption that there is a static Oort cloud. Likely, many of the NICs that will be observed with LSST exited the Oort cloud more than 1 Gyr in the past, when it may have had different physical properties.
Even if the properties of the Oort cloud have not changed over 4.5 Gyr, longlived comets may bias our simulations, as the following argument shows. Suppose, as seems reasonable from Figure 1, that the true distribution of time that a comet spends in the visibility region is given by a power law, i.e.
(2) 
for some in the interval . Then, if we terminate our integrations at some time , the average time that a comet spends in the visibility region during our simulation is
(3) 
In the limit that , this becomes
(4) 
In our simulations, (8, 1, 0, 0) particles with = (5,000, 10,000, 20,000, 50,000) AU survive with AU for the duration of the integration (4.5 Gyr). For AU and AU respectively, the longest lived particles survived for 1.1 and 1.4 Gyr. Now, consider the case of the simulation with AU. We drew particles with lifetimes from the distribution in Equation (2). By chance, we drew no particles with lifetime greater than Gyr. We will assume that we drew randomly from the distribution in Equation (2) subject to the constraint that we draw no particles with lifetimes greater than (this is not quite true, given that we did draw one particle with lifetime exactly ). We can renormalize the distribution in Equation (2) subject to the constraint that no comet have lifetime longer than , and calculate the mean. This is given by (assuming )
(5) 
Taking the ratio of Equations (4) to (5), and assuming that a comet makes a fixed number of appearances in the catalog per unit time spent in the visibility region, we find we have likely underestimated our total flux by something close to a factor of . This expression is 4.0 if we take , as suggested by Yabushita (1979), but declines to unity for . In principle, this bias can be reduced by simulating more comets but the computational resources necessary to reduce the bias substantially are prohibitive.
In the following sections we present simulation data for different values of . is defined relative to the time when the comet first enters within 45 AU except when Gyr, in which case it is defined relative to the start of the simulation. We believe that a value of a few Gyr is most observationally relevant, and most of the following discussion is based on such cutoff times, however given the limitations discussed above, we show results for shorter cutoff times as well.
4 Concentration of Nics Due to the Giant Planets
In this section we describe how the distribution of NICs is affected by the planets.
4.1 Expected Distribution of and with no Planets
We first calculate the expected distribution of orbital radius and perihelion distance in the absence of perturbations from the planets, constructing what we call the zeroplanet model. We assume a uniform distribution of cometary orbital elements in phase space at four fixed energies. These results provide a natural normalization to the plots in the following subsection, which show the distributions of and in our simulated catalog.
Distribution in Radius and Perihelion at a Snapshot in Time
Let there be comets distributed uniformly on the energy surface corresponding to semimajor axis . There is no need to distinguish between the initial semimajor axis and the current semimajor axis here, since the semimajor axis is not changed in the absence of planetary perturbations. Approximating the orbits as parabolic in the visibility region, we find that the radial velocity at radius of an orbit with perihelion distance is
(6) 
Since we assume a uniform distribution of comets on the energy surface, the probability density of the squared eccentricity is
(7) 
Then, since , the number of comets per unit perihelion distance is given by
(8) 
where the last equality holds because we are interested in comets with . A comet on a nearparabolic orbit with perihelion spends a fraction of its time in the radial interval between and , where
(9) 
and is the period of the orbit. Then, using Equations (8) and (9), we can solve for the number of comets in a radial interval :
(10) 
The number of comets with perihelion in the range to expected to be present out to a maximum value of is given by
(11) 
This evaluates to
(12) 
As mentioned in Section 2, because of the way we have set up our simulation, and the fact that we sample every years, we would expect our catalog to contain a number of comets with orbital elements equal to
(13) 
if we had not included any planets in the simulation.
Concentration Factors
In this subsection we compare our catalog to the one which would be produced had we not included planets in our simulations. In Figure 2, we compare our simulated comet catalog with the zeroplanet model (Equations (10) and (13)). We have plotted a “concentration factor” — the ratio of comets appearing in our catalog within a given radial interval to the number calculated from Equations (10) and (13), assuming the same density of comets in the Oort cloud as was used to initialize our simulations. As stated previously, only comets with periods greater than 200 years are counted. Each panel corresponds to a different value of the initial cometary semimajor axis . Different colored lines correspond to different values of as shown in the legend.
These data represent the results from following comets where 10,000, 16,000, 40,000, and 420,000 for 5,000, 10,000, and 20,000 AU, and 50,000 AU respectively. Note that fewer than 10% of these ever enter the region AU (955, 1548, 3949, and 29215 respectively). Most comets do not evolve to AU within the entry period and are therefore discarded at the end of the entry period. Some of those that do come within AU never reach AU (if the angular momentum is nearly perpendicular to the torque).
We would like to know the true distribution of orbital elements in the limit that we simulate a very large number of comets. To estimate our random error, we employed resampling. For each point on the curve, we drew comets with replacement from our simulated comets. The points are the mean of the resampled distribution, and the error bars correspond to the 16 and 84 percentiles of the resampled distribution (if the distribution were normally distributed, these would be 1sigma error bars). The error bars on nearby points are highly correlated in most of our plots because the same comet contributes to several bins in the course of its evolution, hence the low pointtopoint scatter relative to the error bars.
Although these data reflect the orbits of thousands of comets, the statistical errors are large in many cases, since often the majority of the contribution to a particular bin comes from only one or two longlived comets, (see the discussion in §3).
Figure 3 is similar to Figure 2, except that instead of plotting the number of appearances in our catalog in bins of heliocentric radius , we plot appearances in bins of perihelion . The distribution of provides more information about the NIC orbits: a sharp feature in the perihelion distribution is smoothed out when one looks at the distribution of heliocentric radius, since the same orbit can be observed at a range of values of .
In Figures 2 and 3, and in most of the subsequent plots, we have horizontally offset the blue, green, yellow and red curves curves by , , , and respectively in order to make the curves distinguishable in regions where they overlap.
The following qualitative features of Figures 2 and 3 are straightforward to explain:

The flux of comets with AU shows a sharp dropoff (in both and ) interior to 10 AU. This is because comets originating at small semimajor axes are subjected to weak Galactic tides and change their perihelia slowly. The majority of kicks given to a comet with AU are large enough to either unbind a comet infalling from outside a few thousand AU, or to reduce to the point that the timescale to change the perihelion distance is much longer than an orbital period. We therefore expect to see a jump in the number of comets appearing inside 5 AU at values of exceeding that at which a comet can go from perihelion AU (largely unaffected by Jupiter and Saturn) to perihelion AU in one orbit. This occurs at approximately 30,000 AU. Therefore, in order for a comet starting from inside AU to appear in the inner solar system, it needs to have a lucky orientation with respect to Jupiter and Saturn
^{5} . This lucky orientation can either yield multiple small energy kicks on subsequent perihelion passages, or yield a kick that increases the semimajor axis so that the comet receives a larger torque from the Galaxy (Kaib & Quinn 2009).Thus we conclude that comets with AU are mostly ejected by interactions with the outer planets before they reach small heliocentric radii. Hills (1981) arrives at a similar result considering the effects of passing stars discretely rather than as a smooth Galactic tidal potential. Collins & Sari (2010) provide a discussion of when it is appropriate to treat the influence of the Galaxy as being due to a smooth tidal field, and when it should be modeled as discrete stellar encounters.
The lower concentration factors for comets with smaller values of do not mean that we will see fewer NICs for a given number of Oort cloud comets at that energy. This is because the fraction of the total comets that are in the visibility region at a given time in the zeroplanet model scales with (see Equation (10)).

The difference between the green curves and the orange and red curves grows with increasing . This is because the kicks from the planets are smaller at large , so the comets survive longer. The exception to this is the curve for Myr and 5,000 AU. In this case, comets have mostly had insufficient time to be torqued to small values of . It should be noted however, that the time to reach a given perihelion distance is not completely determined by the initial semimajor axis because a comet could be scattered to larger by an early encounter with Neptune, and subsequently evolve more rapidly.
There is little difference between the results for Gyr, and Gyr for comets with AU. As discussed previously, this is likely an artifact of our simulations including insufficiently many comets with these values of to capture the tail of the lifetime distribution (see Section 3).

We do not expect to drop exactly to unity as soon as is larger than the extent of the planetary perturbations, because NICs which have interacted with the planets may be systematically carried away from the planetary region by the tide at a different rate than they were carried in (due to a change in orientation or semimajor axis).
5 Distribution of Orbital Elements for Visible Comets
The above analysis shows the degree to which the giant planets concentrate NICs in the outer solar system and exclude them from inside the orbit of Jupiter. In this section, we use an estimate of the size distribution of NICs, the relationship between magnitude, size and heliocentric distance, and the concentration effect due to interactions with the planets to calculate the number of NICs expected to be seen in an allsky survey as a function of the limiting magnitude .
5.1 Size Distribution
Comets have so far been observed primarily within a few AU of the Sun, where their brightness is influenced strongly by their activity. At the larger distances that we focus on here, comets are believed to be generally inactive (see discussion in Section 2.2), so their brightness is determined solely by their size, distance, and albedo. Let be the apparent magnitude of a comet 1 AU from the Sun and 1 AU from the observer, seen from zero phase angle. Based on a sample of longperiod comets from about to , Sosa & Fernández (2011) derive a relation for active comets between the radius (in kilometers) and :
(14) 
where = 2.072 and . We caution that our use of this formula requires substantial extrapolation: the largest comet used to determine the formula has a radius of 1.8 km, more than an order of magnitude smaller than the smallest comets detectable at 30 AU in a survey with the limiting magnitude of LSST (see Section 5.2). Sosa & Fernández (2011) note that the relation in Equation (14) predicts a radius (13 km) for comet HaleBopp that is somewhat below other estimates (mostly falling in the 20–35 km range). This discrepancy suggests that Equation (14) may underpredict the radii of large comets, in which case our estimates of the observable comet population will be conservative. Note that by using Equation (14) we are assuming that longperiod comets are mostly the same population as NICs (or at least have the same size distribution).
Hughes (2001) finds that the number of longperiod comets with brightness passing through perihelion in the inner solar system per year per AU of perihelion distance is given by
(15) 
with and . We can then transform variables to using Equation (14). We find that
(16) 
This distribution holds down to km, however, for simplicity we extrapolate it down to km — the smallest comet visible at 5 AU in our model (see next section). This size distribution leads to a weak divergence in total mass at the large end of the spectrum. Nevertheless, we assume that the powerlaw behavior holds up to several tens of kilometers. The size distribution in Equation (16) is steeper than the relation () estimated in Hughes (2001) because he uses a different relation between and . It is also substantially steeper than the relation () found in Snodgrass et al. (2011) for the Jupiterfamily comets. If the size distribution is shallower than we have estimated at large radii, then our estimates of the observable comet population will be conservative.
5.2 Visibility Model
In this section we describe our model for determining how likely a given simulated comet is to be visible. We assume that comets have an band geometric albedo of 0.04 as suggested in Lamy et al. (2004). We find that the magnitude of an inactive comet is
(17) 
where we have used as the apparent band magnitude of the Sun. Equation (5.2) is only valid for comets far from the Sun, since we have assumed that , that the phase angle is zero, and most importantly, that we are seeing the bare nucleus of an inactive comet. Lamy et al. (2004) state that magnitude drops off at about a rate of 0.04 magnitudes/degree of phase angle, meaning that error due to this nonzero phase angle is limited to at most 0.23 magnitudes for a comet at 10 AU. Similarly, at 10 AU, the largest possible error arising from the approximation that the Suncomet distance is equal to the Earthcomet distance also corresponds to a magnitude error of .
5.3 Weighting of Observations
Using Equation (5.2) we can solve for , the radius in kilometers of the smallest comet visible at distance in a survey with limiting magnitude , assuming :
(18) 
The comet size is not explicitly tracked in our simulations. We assume that the sizes of comets in our simulation are drawn from the distribution in Equation (16), with a lower cutoff radius of km — the smallest comet visible at 5 AU in our model. To account for the fact that not all comets are visible at all orbital radii, we weight a simulated comet appearance at radius by the fraction of comets that would be visible at the observed value of given the assumed size distribution in the simulation. An appearance at high will receive a low weight, since most comets would not be visible so far away. We assign weight 1 to observations at AU, as all comets in our assumed size distribution would be visible at 5 AU. Then at general , we assign weight
(19) 
We quote numbers of comets with a given set of orbital elements per comets larger than one kilometer in a spherical distribution at the assumed initial semimajor axis. To achieve this normalization, we multiply our counts by
(20) 
where is the fraction of the phase space of orbits with semimajor axis that consists of orbits with perihelia between 60 and 65 AU, given by
(21) 
and is the number of comets we initialize between 60 and 65 AU.
5.4 Distribution of Visible Comets
In this section we present results from our simulations showing how many NICs are visible over the whole sky at a given snapshot of time as a function of and , for observations taken between and . In all cases, we assume a limiting band magnitude of 24.5, equivalent to the oneexposure limit for LSST (LSST Science Collaboration et al. 2009). The number of distant NICs expected to be discovered by LSST differs from the results presented here for two reasons. First, LSST is expected to operate for 10 years, so it should see more than just the comets visible in a snapshot, particularly in the case of the closer comets where is less than 10 years. Second, LSST will only survey 48% of the sky, so will only see about half of the comets that would be visible in an equivalent allsky survey. Comets move slowly enough that trailing losses will be insignificant given the 30 second exposure time. Using Equation (8) from Ivezić et al. (2008), we estimate a comet at 10 AU will have a limiting magnitude only 0.06 magnitudes brighter due to trailing losses.
Figures 4 and 5 show the number of NICs expected to be visible outside AU per unit of and respectively, per comets with greater than 1 km at the labeled initial semimajor axis in the Oort cloud. The shapes of the curves are substantially different for different values of , particularly in the region between 5 and 10 AU, where the statistics are the best. In the 5,000 AU case, the expected count increases by a factor of 5 from 5 AU to 10 AU for Gyr. In the 50,000 AU case it decreases by a factor of . Observations of comets in this regime will therefore allow us to observationally constrain the distribution of . The peak of moves smoothly from around 15 AU for AU to less than 5 AU for AU.
As shown in Appendix B, we expect and to decline as and respectively in the zeroplanet model. Deviations from this behavior are due to variation in the concentration factor as shown in Figures 2 and 3.
We also examined what happened if we broke up the sample into two groups depending on the current semimajor axis of the comet. Figure 6 is identical to Figure 5 except that we have only considered those comets that have semimajor axes greater than 300 AU. The error bars are smaller, because the comets with the most appearances tend to diffuse to smaller values of , leaving a population with less spread in number of appearances. For this reason, this subset of comets, although smaller in number, has more power to discriminate between Oort cloud models.
In Figure 7 we plot for only the comets in Figure 5, but not Figure 6, i.e., those comets whose orbits have AU. We see that declines more sharply with than in the whole sample of longperiod comets. This is because it is difficult for a comet to attain AU at large perihelion, because the kicks are too small. We also note that the overall numbers are larger by a factor of a few for comets with AU.
5.5 Distribution in Semimajor Axis
Figure 8 shows the semimajor axis distribution (number of appearances per unit logarithmic interval in semimajor axis) for all the NICs in a given perihelion bin (denoted by the color of the curve) and initial semimajor axis (panel). The error bars are generally larger for the points at small semimajor axis, implying that the statistics in these bins are dominated by a few comets.
We find it more illuminating to make this plot in the coordinate , which is proportional to the energy. Doing so allows us to test the predictions of random walk models such as those in Yabushita (1979) and Everhart (1976). In their simplest form, one can imagine a comet starting with energy . Every perihelion passage, it takes a step of size towards higher or lower energy. It is removed if (it becomes unbound), or if , where is the critical energy level to be reclassified as a shortperiod comet. We ignore the possibility that a short period orbit could be perturbed back to the longperiod regime. Comets are injected near . In the limit that , one can show that a steadystate distribution of comet energies normalized by the rate of perihelion passage is given by a linear equation of the form
(22) 
where is a constant depending on the comet injection rate and the size of the kick.
We see some support for Equation (22) in Figure 9. We have only considered comets with AU ( AU). This is a small enough range that we would expect the curves to be nearly flat if Equation (22) were correct (since is much more negative than the energies shown in Figure 9). The curves are generally flat for perihelion distances 5 AU 30 AU. Inside 5 AU, the error bars are too large to say much. Outside 30 AU, there is a definite trend favoring energies closer to 0, as would be expected since the diffusion rate at these perihelion distances is so slow that the steadystate energy distribution cannot be achieved.
We see evidence of a weak “Oort spike” when AU  an excess of comets with AU interpreted as the result of the initial conditions of cometary orbits. This is seen only for comets with AU, since for larger values of , the spike is drowned out by the large numbers of comets which remain at more negative energies for many orbits. As a quantitative measure of the observed Oort spike, Fernández & Sosa (2012) find that 30% of comets discovered since the year 1900 with AU have (10,000 AU). We find that we cannot come close to reproducing this unless we assume some disruption. We applied a disruption law in which comets were removed after spending years within 5 AU of the sun. For comets with 50,000 AU and AU, the fraction seen with AU is {, , }, for = {100, 1,000, 10,000 years}. This confirms the wellknown result that one cannot reproduce the Oort spike without some model for comet disruption or fading that removes comets after only a few appearances (e.g., Wiegert & Tremaine 1999).
We also notice a broader “spike” of comets around the original energy for comets with or 10,000 AU, and perihelion distances outside 30 AU. A typical kick in energy at 30 AU might be on the order of AU (Duncan et al. 1987), so comets with small enough to not be rapidly carried away by the tide will show a broad peak around their original energy at large perihelion distance.
5.6 Orbital Inclination
Figure 10 shows the orbital inclination distribution (in ecliptic coordinates) of the visible NICs. The longest period NICs are preferentially retrograde. In a random walk model, the density of comets varies inversely with the step size. The energy kick is 23 times larger for a prograde comet (Duncan et al. 1987), translating into an expected prograde fraction of to .
We find that our error bars on the prograde fraction are large for the entire comet population, so for the following analysis we split the population into two groups depending on the semimajor axes of the comets. Once again we find that elements of comets with large values of have less statistical noise. We obtain prograde fractions among the visible comets of , and for 5,000, 10,000, 20,000, and 50,000 AU respectively, for the comets with AU. For comets with AU, we find prograde fractions of , and . These data are consistent with the random walk model.
While the simulation data agree with the random walk model, they contradict the observations. There is only a slight preference in the observational data for retrograde comets with high perihelion. 64 out of 110 comets (58%) with period greater than 200 years and perihelion greater than 5 AU in the database at http://ssd.jpl.nasa.gov/dat/ELEMENTS.COMET are retrograde.
5.7 Size Distribution
A bigger telescope enables us to see rare large comets because it can search more volume. It is impossible to say exactly what size distribution of comets to expect in the observed sample, however we can make an estimate based on extrapolation of the size distribution from Fernández & Sosa (2012). It is instructive to first consider the zeroplanet model with a fixed powerlaw for the size distribution.
In Appendix B, we calculate the number of NICs greater than a certain size expected to be seen in the zeroplanet model. We find in Equation (B) that the number visible with size greater than is given by assuming comets greater than 1 km at 10,000 AU, but the results depend sensitively on the assumed number density and semimajor axis. Additionally, due to the concentration at larger orbital radii, we actually expect to see objects much larger than one would infer from the zeroplanet model used to derive Equation (B).
Figure 11 shows the expected number of observations of NICs larger than a given size in the distance range 5 AU 45 AU. This calculation assumes comets greater than 1 km at the given initial semimajor axis in the Oort cloud and slope of the size distribution given in Equation (16). Exclusion of comets with orbital radii less than 5 AU should have a negligible effect on the counts except at the smallest sizes, as the concentration factors are lower there, and it is a negligible fraction of the visible volume for the larger comets. Exclusion of comets with orbital radii greater than 45 AU has no effect for the size of comets that we have plotted, since none of them would be visible past 45 AU. We see that assuming the size distribution from Equation (16) holds to these sizes, we observe on the order of 10 NICs greater than 80 km in size if the majority of the NICs are coming from 5,000 AU, and a few if the majority of the NICs are coming from 10,000 AU.
6 Conclusions
We simulated the evolution of NICs originating in the Oort cloud as they interact with the giant planets. We used these simulations to create a catalog of simulated comet positions and velocities. We observe different distributions of orbital elements including perihelion distance, semimajor axis and inclination depending on the semimajor axes at which the NICs originate. Observations by LSST will therefore let us determine the absolute numbers of comets in the Oort cloud as a function of semimajor axis, and test Oort’s standard model for the origin of comets. The distribution of NICs outside the orbits of Jupiter and Saturn will provide direct evidence for the presence or absence of the hypothesized “inner Oort cloud” corresponding to semimajor axes between and AU.
One surprising result is that we expect at least tens of percent of the comets observed by LSST in the outer solar system to have been interacting with the giant planets for more than 1 Gyr. This result makes the interpretation of the comet population detected by LSST more difficult and interesting, since the population and spatial distribution of comets in the Oort cloud almost certainly evolves on timescales of a few Gyr.
We will also get a better measurement of the Oort spike — the excess of comets in nearly parabolic orbits — as we will be able to measure highprecision orbits in a regime where comets are likely unaffected by nongravitational forces. This will put constraints on models of comet fading, as well as the original semimajor axes of comets.
We will furthermore be able to constrain the size distribution of Oort cloud comets out conservatively to several tens of kilometers, and perhaps even to larger bodies, depending on the size distribution and number density of comets. Our results are based on a relatively steep slope for the size distribution of comets, (Equation 16) and may substantially underestimate the total number of comets that will be discovered at large distances by LSST.
We thank the referee, Ramon Brasser, for helpful and constructive comments and advice.
Appendix A Maximum Torque on a Radial Orbit
We calculate the torque on a radial orbit from the Galactic tide. This calculation allows us to exclude comets from the Oort cloud that will definitely not enter the buffer region during the entry period as described in Section 2. The results are also used in Section 4.1.2. While our method is approximate, we have checked that our prediction is conservative in the sense that no comet can enter the buffer region which we predicted could not enter the buffer region. For simplicity, we made the following approximations.
In calculating the torque, we assume that the orbit is completely radial, pointing in the direction of apocenter. This is welljustified by the characteristic ratios of . We ignore the components of the Galactic tidal field in the plane of the Galaxy, as they are smaller than the outofplane () component by about a factor of 10.
Heisler & Tremaine (1986) give the tidal force (per unit mass) as
(A.1) 
where , and the Oort and constants are taken to be 14.6 km s kpc and km s kpc respectively (Binney & Tremaine 2008). This leads to an instantaneous torque given by
(A.2) 
where is the Galactic latitude of the comet (which is constant in the radial orbit approximation). The orbit averaged torque is calculated by noting that .
Appendix B Comets Visible Assuming a Uniform Distribution of Orbital Elements on an Energy Surface and No Planets
In this appendix, we calculate the number of comets observable at a snapshot in time as a function of heliocentric radius in the zeroplanet model. We approximate comet orbits as parabolic in the observation region. If we have a population of comets with perihelion , of which pass perihelion per year that are large enough to be visible at , then the density of comets visible in a radial interval at radius , is
(B.1) 
where is the radial velocity (Equation (6)). Using Equations (16) and (18), we find that in the zeroplanet model,
(B.2) 
where is the size (in km) of a body that is marginally detectable at 1 AU in our model (see Equation (18)). Therefore, using Equations (6), (B.1) and (B.2) we can write the density of comets as an integral over :
(B.3) 
We also calculate how many comets we see as a function of perihelion given a magnitude cut at a given snapshot in time. We find that
(B.4) 
where is the amount of time that a comet with radius in kilometers and perihelion in AU is visible during one perihelion passage. From Equation (6), we see that
(B.5) 
where is the maximum heliocentric distance at which the comet is visible. Putting this together, we find that
(B.6) 
Finally, we ask how many comets we expect to see above a certain size in the zeroplanet model. In this calculation, we take the size distribution in Equation (16), but normalize the counts to those expected if there were to comets with semimajor axis and radii greater than 1 km. A cloud of comets with orbital elements uniformly distributed on a surface of constant energy in phase space will have comets passing perihelion per year per AU perihelion, for . Therefore, the correct normalization is
(B.7) 
In this model, we can expect to see comets larger than where is given by
(B.8) 
Setting this equal to , and solving for , we find that we expect to see one comet larger than 14.4 km if = 10,000 AU, but this estimate is clearly quite sensitive to the assumed number of comets and that value of .
Footnotes
 affiliation: Department of Astrophysical Sciences, Princeton University, Ivy Lane, Princeton, NJ 08544; ksilsbee@astro.princeton.edu
 affiliation: Institute for Advanced Study, 1 Einstein Drive Princeton, NJ 08540; tremaine@ias.edu
 Because the boundary between the buffer and Oort region at 60 AU corresponds to a perihelion distance twice the semimajor axis of Neptune’s orbit, we expect the planets to have a negligible effect on the orbits of comets in the Oort region. Therefore, it is reasonable to assume that a comet with an orbit aligned such that the Galactic tide pulls it from the buffer region into the Oort region will not return for a long time.
 http://ssd.jpl.nasa.gov/?horizons
 Or a kick from a star that passes unusually close to the Sun, a rare event not included in our model.
References
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton, NJ: Princeton University Press)
 Brasser, R., Duncan, M. J., Levison, H. F., Schwamb, M. E., & Brown, M. E. 2012, Icar, 217, 1
 Collins, B. F., & Sari, R. 2010, AJ, 140, 1306
 Dones, L., Weissman, P. R., Levison, H. F., & Duncan, M. J. 2004, in Astronomical Society of the Pacific Conference Series, Vol. 323, Star Formation in the Interstellar Medium: In Honor of David Hollenbach, ed. D. Johnstone, F. C. Adams, D. N. C. Lin, D. A. Neufeeld, & E. C. Ostriker, 371
 Duncan, M., Quinn, T., & Tremaine, S. 1987, AJ, 94, 1330
 Everhart, E. 1976, NASA Special Publication, 393, 445
 Fernández, J. A. 1981, A&A, 96, 26
 —. 2005, Comets: Nature, Dynamics, Origin, and their Cosmogonical Relevance (Dordrecht: Springer)
 Fernández, J. A., & Sosa, A. 2012, MNRAS, 423, 1674
 Fouchard, M., Rickman, H., Froeschlé, C., & Valsecchi, G. B. 2013, Icarus, 222, 20
 —. 2014a, Icarus, 231, 110
 —. 2014b, Icarus, 231, 99
 Francis, P. J. 2005, ApJ, 635, 1348
 Heisler, J., & Tremaine, S. 1986, Icar, 65, 13
 Heisler, J., Tremaine, S., & Alcock, C. 1987, Icar, 70, 269
 Hills, J. G. 1981, AJ, 86, 1730
 Hughes, D. W. 2001, MNRAS, 326, 515
 Ivezić, Ž., Tyson, J. A., Abel, B., et al. 2008, ArXiv eprints, arXiv:0805.2366
 Kaib, N. A., & Quinn, T. 2009, Science, 325, 1234
 Kramer, E. A., Fernandez, Y. R., Lisse, C. M., Kelley, M. S. P., & Woodney, L. M. 2014, Icar, 236, 136
 Lamy, P. L., Toth, I., Fernandez, Y. R., & Weaver, H. A. 2004, in Comets II, ed. M. Festou, H. Keller, & H. Weaver (Tucson, AZ: University of Arizona Press), 223–264
 Levison, H. F., Dones, L., & Duncan, M. J. 2001, AJ, 121, 2253
 LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, ArXiv eprints, arXiv:0912.0201
 Marsden, B. G., & Williams, G. V. 1993, Catalogue of Cometary Orbits 1993. Eighth edition. (Cambridge, MA: Minor Planet Center)
 Rein, H., & Liu, S.F. 2012, A&A, 537, A128
 Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424
 Schwamb, M. E., Brown, M. E., Rabinowitz, D. L., & Ragozzine, D. 2010, ApJ, 720, 1691
 Snodgrass, C., Fitzsimmons, A., Lowry, S. C., & Weissman, P. 2011, MNRAS, 414, 458
 Snodgrass, C., Tubiana, C., Bramich, D. M., et al. 2013, A&A, 557, A33
 Sosa, A., & Fernández, J. A. 2011, MNRAS, 416, 767
 Stokes, G. H., Evans, J. B., Viggh, H. E. M., Shelly, F. C., & Pearce, E. C. 2000, Icar, 148, 21
 Tsujii, T. 1992, Celestial Mechanics and Dynamical Astronomy, 54, 271
 Weaver, H. A., Feldman, P. D., A’Hearn, M. F., & Arpigny, C. 1997, Science, 275, 1900
 Weissman, P. R. 1980, A&A, 85, 191
 Wiegert, P., & Tremaine, S. 1999, Icar, 137, 84
 Yabushita, S. 1979, MNRAS, 187, 445
 Yeomans, D. K., Chodas, P. W., Sitarski, G., Szutowicz, S., & Królikowska, M. 2004, in Comets II, ed. M. Festou, H. Keller, & H. Weaver (Tucson, AZ: University of Arizona Press), 137–151