The Feasibility and Benefits of In Situ Exploration of 'Oumuamua-like Objects

The interstellar asteroid 1I/2017 U1 ('Oumuamua), presented astronomers with their first glimpse of an object entering the inner solar system with a hyperbolic trajectory. The detection of this remarkable body was made on 2017 October 19 by the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) project (Chambers et al. 2016), and its extrasolar origin was rapidly confirmed, leading to a 2017 October 25 discovery announcement through the Minor Planet Center (IAU Minor Planet Center 2017a). Its securely determined perihelion distance, p = 0.255 au, and eccentricity, e = 1.2, imply a pre-encounter velocity ${v}_{\infty }=26\,\mathrm{km}\,{{\rm{s}}}^{-1}$. No known or proposed outer solar system objects are remotely capable of imparting such an impulse; 'Oumuamua is a singular outsider.

In the weeks following 'Oumuamua's discovery, observations by many groups provided significant additional detail concerning its properties. The near-daily flux of articles was large enough that a review of the current state of knowledge is now useful.

Moderate-resolution spectra of 'Oumuamua were obtained by Masiero (2017) and Ye et al. (2017) immediately after the Minor Plant Center's circulation of the discovery announcement. These measurements showed no sign of narrow emission or absorption features and pointed to a spectral gradient S' that is skewed toward the red. Ye et al. (2017) measured a gradient S' = 0.1 ± 0.06(100 nm)⁻¹ at 650 nm, based on a 400–900 nm bandpass. This value for S' is consistent with that displayed by cometary nuclei, Trojan asteroids, and inactive Centaurs, but is less red than active Centaurs and all classes of Kuiper Belt objects, which typically exhibit ultrared gradients S' = 0.23 ± 0.02(100 nm)⁻¹ in the optical region (Jewitt 2015). Bannister et al. (2017) obtained g', r', and J photometry, which extended 'Oumuamua's color assessment beyond 1 μm and slightly into the near-infrared. Their results indicate a spectral slope S' = 0.22 ± 0.15(100 nm)⁻¹, consistent within the errors to the results obtained in the visible region of the spectrum. Additionally, Ye et al. (2017) searched radar meteor data taken during the relevant arrival time and found that no incoming meteors appeared to have been associated with 'Oumuamua (which passed only 0.161 au from Earth at close approach).

The arrival direction of 'Oumuamua's incoming trajectory was located close to the solar apex (the direction of the Sun's peculiar motion with respect to the Local Standard of Rest), providing strong circumstantial evidence that it is a member of a large population of planetesimals ejected from planetary systems throughout the Milky Way's disk during their formation phases (Laughlin & Batygin 2017; Trilling et al. 2017). Its kinematics suggest that it belongs to a population with kinematics similar to the Population I stars of the solar neighborhood (Mamajek 2017). Indeed, Wright (2017) reiterated that known (or unknown) planets in our own solar system have no chance of scattering 'Oumuamua from an Oort cloud origin into its hyperbolic trajectory.

Several authors have discussed the possibility that specific source stars of origin for 'Oumuamua might be identifiable. Gaidos et al. (2017) draw on velocity similarities to suggest that 'Oumuamua was formed ∼45 Myr ago in a protoplanetary disk of one of the 50–85 pc distant stars of the Carina/Columba association. Zuluaga et al. (2017) presented a general method to assign probabilities of specific nearby stars being parent systems for interstellar asteroids, and they presented their most likely source candidate for 'Oumuamua as the double system HD 200325. Dybczyński & Królikowska (2018) carried out a detailed dynamical assessment of more than 200,000 well-characterized nearby stars, with the aim of identifying potential source systems, and found no particularly compelling candidates. Intriguingly, Dybczyński & Królikowska do note that roughly 820,000 years ago, 'Oumuamua encountered the nearby planet-bearing star Gliese 876 (Rivera et al. 2010) with a ${v}_{\infty }=5\,\mathrm{km}\,{{\rm{s}}}^{-1}$ relative velocity and an impact parameter b ∼ 2.2 pc. Dybczyński & Królikowska (2018) conclude that this encounter likely occurred at too large a distance to admit a physical connection when the various uncertainties are taken into account. Zhang (2018) emphasized the general difficulties associated with back-tracing interstellar objects (ISO) to their source systems, noting in particular that points of origin are obscured within tens of millions of years by the chaotic nature of orbits in the galactic potential.

Remarks have also given consideration to the idea that 'Oumuamua was originally ejected from the solar system and has re-encountered the Sun. This possibility has an exceptionally low probability. The phase mixing (Binney & Tremaine 2008) implied by 'Oumuamua's large heliocentric ${v}_{\infty }$ implies that its trajectory was fully independent from that of the Sun's prior to its celebrated recent encounter. Adopting a σ = π au² cross section (which includes the Sun's gravitational focusing), a $\langle v\rangle \,=25\,\mathrm{km}\,{{\rm{s}}}^{-1}$ velocity, and solar excursions of ΔR = ±350 pc and ΔZ = ±350 pc centered on an R = 8.5 kpc mean galactocentric radius, one calculates an expected wait time $\tau =(2\pi R{\rm{\Delta }}R{\rm{\Delta }}Z)/(\sigma \langle v\rangle )\sim {10}^{24}$ years between departure from the Sun's vicinity and its subsequent return.

Although the physical origin of 'Oumuamua is speculative, a number of observations provide clues to its composition and layout. At periastron, 'Oumuamua passed within 0.25 au of the Sun, and despite achieving solar irradiation levels, I > 20 kW m⁻², deep images of the object produced no sign of a coma (IAU Minor Planet Center 2017b), demonstrating that it has a nonvolatile near-surface composition. Naively, this seems surprising, since estimates of the ratio of cometary to asteroidal material in the Oort cloud range from 200:1 to 10,000:1 (Meech et al. 2017). Strong color-independent time variability in the photometric light curve observations presented by Knight et al. (2017), Jewitt et al. (2017), Meech et al. (2017), and Bolin et al. (2018) suggests that the object is highly elongated and moderately rapidly rotating, with an aspect ratio of up to 10:1, a mean radius of r_o = 102 ± 4 m, and a P_rot ∼ 7.3 hr rotation period. This rate of spin is sufficient to overwhelm the self-gravity of a rubble pile with ρ ≲ 6 g cm⁻³, implying that 'Oumuamua has nonnegligible tensile strength. Fraser et al. (2018) and Drahus et al. (2017) independently analyzed 'Oumuamua's photometric light curve measurements and concluded that it is exhibiting nonprincipal axis rotation. Like the Saturnian satellite Hyperion (Wisdom et al. 1984), 'Oumuamua may be tumbling chaotically.

Domokos et al. (2017) propose (drawing on an Eikonal model) that 'Oumuamua's highly elongated shape is the expected consequence of abrasion incurred by extended passage through a sea of tiny particles, and Gaidos (2017) notes that the large aspect ratio could potentially arise if 'Oumuamua is a contact binary composed of two less-extremely elongated ellipsoidal components with heterogeneous surfaces.

'Oumuamua's unusual observed properties generated discussion of a full range of theories regarding its provenance. Ensuing attempts to detect technologically generated radio signals from the asteroid using eight hours of time on the Robert C. Byrd Green Bank Telescope were, however, unsuccessful (Enriquez et al. 2018).

Planet formation theory suggests that giant planet migration may be responsible for the ejection of very large quantities of planetesimals into interstellar space (Tsiganis et al. 2005; Levison et al. 2008). If 'Oumuamua is a planetesimal ejected from another system, then the observation of one such comet during Pan-STARRS's operational lifetime to date implies that an average of approximately $1\,{M}_{\oplus }$ of material is ejected from every star (Laughlin & Batygin 2017; Trilling et al. 2017). Raymond et al. (2018b) point out that this estimate assumes a monosized population of ISOs, and that if the true planetesimal mass distribution is top heavy, similar estimates yield an implausibly high mass of ejected material for each star. They demonstrate that this may be reconciled if 'Oumuamua is a member of a subpopulation of 0.1%–1% of planetesimals that are tidally disrupted by gas giant planets on their pathway to ejection.

The advent of 'Oumuamua, which was detected in spite of a number of observational biases, suggests that a large population of similar objects may pervade the galaxy. Extrapolating from a sample size of one and using estimates based on cross section, Laughlin & Batygin (2017) calculate that the Galaxy contains N ∼ 2 × 10²⁶ 'Oumuamua-like planetesimals. This result concurs with estimates by Jewitt et al. (2017), Meech et al. (2017), and Trilling et al. (2017), who estimate space densities of order n ∼ 0.1 au⁻³, a value implying that on average, of order one such object can always be found traversing an R = 1 au sphere centered on the Sun. Most recently, Do et al. (2018) carried out a detailed assessment of Pan-STARRS's ability to detect 'Oumuamua-like objects. Their evaluation folds in detectability as a function of direction of approach and considers possible power-law indices for the unknown size distribution. Do et al. (2018) find an interstellar number density n = 0.2 of similar objects per au³. Clearly, additional detections will be both welcome and informative.

Jewitt et al. (2017) draw analogies with the composition and behavior of comets to argue that 'Oumuamua's lack of coma is not inconsistent with an ice-rich bulk composition. Long-duration exposure to the radiation environment of the interstellar medium is capable of mantling a comet-like object with a crust of volatile-depleted "regolith" of very low thermal conductivity (Cooper et al. 2003). Fitzsimmons et al. (2018) present a more detailed thermal model that supports the discussion in Jewitt et al. (2017). Raymond et al. (2018a) posit that 'Oumuamua is a fragment of a tidally disrupted comet-like, icy planetesimal that underwent enough close passages to its host star to become extinct. Additionally, Jewitt et al. (2017) raise the interesting possibility that heat transfer into 'Oumuamua may cause it to develop an observable coma as it departs the solar system. A related phenomenon has been observed for select comets, including comet P/Halley, which produced a postperihelion outburst at 14 au that has been attributed to release of trapped gases during a phase transition at depth from amorphous to crystalline ice (Prialnik & Bar-Nun 1992). It is unlikely, but nonetheless possible, that HST observations of 'Oumuamua, completed in 2017 November and December and at the beginning of 2018,¹ may show signs of coma, thereby placing important constraints on both the surface and interior properties.

Ćuk (2018) and Jackson et al. (2017), however, find no particular surprise in observing an ISO with an asteroidal bulk composition. These authors argue that 'Oumuamua was ejected by a binary system, in which case it (and many bodies similar to it) may naturally be composed of nonvolatile material. Hansen & Zuckerman (2017) point to the observed rate of pollution of white dwarfs by asteroidal debris, and they argue that accompanying ejection of such material from host systems must also occur. Moreover, they note that large post-main-sequence stellar luminosities will lead to volatile loss for many ejected bodies. As a point of reference, the Sun is expected to achieve a peak luminosity of $5200\,{L}_{\odot }$ (Sackmann et al. 1993) and will lose a significant burden of its minor-body retinue in the dynamical reshuffling that will befall the solar system as the Sun loses mass during its transition into a white dwarf (Duncan & Lissauer 1998). Rafikov (2018) also considers the hypothesis that 'Oumuamua's lack of coma implies a refractory composition and demonstrates that tidal disruption of asteroidal bodies that make close approaches to white dwarf stars can potentially explain its origin.

Given the range of possible origins, it would be useful to investigate a body such as 'Oumuamua in close proximity, especially in connection with the use of a kinetic impactor to excavate a debris plume that could be examined spectroscopically. Projects of this type certainly have precedent. Over the past three decades, a rich heritage of in situ investigations of comets, asteroids, and Kuiper Belt objects has been established, and in particular, a number of missions have been targeted to solar system comets. The first missions of this type included the NASA International Sun-Earth Explorer (ISEE-3), which was directed through the tail of Comet Giacobini-Zinner in 1985 (Tsurutani et al. 1986), and the Vega 1, Vega 2, Sakigake, Suisei, and Giotto missions, which all investigated Comet Halley in 1986 (LePage 2011).

NASA's Deep Space 1 Mission encountered Comet Borrelly in 2001 (Nelson et al. 2004). In 2004, NASA's Stardust Mission flew to within 240 km of the nucleus of comet Wild 2 (Tsou et al. 2004), collected dust particles by both trapping them in aerogel and cratering them on foil, and returned the payload to Earth in 2006 (Brownlee et al. 2006; Flynn et al. 2006; Sandford et al. 2006). During the near-contemporaneous Deep Impact Mission, a purpose-built kinetic impactor was directed to strike the comet Tempel 1, producing a debris plume that was observed both from Earth and from the Mission's main spacecraft. The impact excavated contents from beneath the comet's radiation-hardened surface, allowing an assessment of pristine planetesimal-forming material (A'Hearn et al. 2005).

More recently, ESA's Rosetta Mission (Glassmeier et al. 2007) deployed a small lander on the surface of comet 67/Churyumov–Gerasimenko, while the main Rosetta craft orbited at a distance of d ∼ 20 km and extensively imaged the comet's surface. At mission's end in 2016, the main craft was directed to strike the comet, returning extremely close-up images as it did so (El-Maarry et al. 2015). Finally, although not directed to visit a comet, NASA's Osiris-REx mission is flying to the asteroid Bennu, where it will sample the surface and return the contents to Earth (Lauretta & OSIRIS-Rex Team 2012).

The inference that extrasolar planetesimals frequently infiltrate the inner solar system raises the possibility of designing a current-technology space mission to intercept and investigate one of these objects at close range. Such a mission would face a number of challenges, including (1) the large heliocentric velocities of objects on hyperbolic trajectories, (2) the lack of substantial time following the discovery of the target object for mission planning and execution, and (3) uncertainty in targeting during final approach. It is worth noting that when 'Oumuamua was detected and announced in late 2017 October, it had already passed its periastron location (which occurred on 2017 September 9), and indeed, was already more than 1 au from the Sun.

Hein et al. (2017) reported a preliminary investigation of mission designs aimed at intercepting 'Oumuamua itself. At present (and especially going forward), such missions are very challenging to mount, due to the large Δv's required to catch up with the rapidly departing object. Hein et al. (2017) considered trajectories that include direct transfers, as well as missions that employ a combination of Jovian gravity assist and Solar Oberth maneuvers (see, e.g., Bond & Allman 1996), and they speculated on the potential use of advanced technologies such as solar sails, laser propulsion, and magnetospheric "chipsat" acceleration.

In this paper, we adopt a complementary approach and assess the near-term prospects for intercepting future 'Oumuamua-like arrivals using conventional chemical propulsion. We envision a mission with scientific charge that is similar to the Deep Impact Mission, in which a kinetic impactor strikes an impinging ISO at high velocity, enabling a companion flyby probe to carry out a spectroscopic examination of its surface and, crucially, its subsurface contents. The feasibility of such missions hinges primarily on whether an incoming body can be detected with sufficient lead time to launch an interceptor that can meet an ISO traveling on a hyperbolic trajectory. An evaluation of this problem thus involves coupling detectability analyses (Moro-Martín et al. 2009; Cook et al. 2016) to trajectory evaluations for a population of incoming bodies.

The plan for this paper is as follows. In Section 2, we outline a thermal model for 'Oumuamua during its passage through the inner solar system, allowing us to place quantitative limits on the composition and the thermal conductivity of its surface. This analysis delineates that the primary scientific goal for any mission to such a body would be to determine whether the regolith-sheathed ice model is correct; to augment efforts by other authors, the model that we consider takes the elongated shape, the full mass, and the chaotic tumbling of the body into account. In Section 3, we assess the likelihood of arrival directions of future ISOs as a function of arrival direction (based on a heliocentric viewpoint), and we predict detection rates as a function of observing season (on Earth) and limiting survey magnitude. In Section 4 we evaluate the expected fuel (or specifically Δv) requirements for short-notice direct-transfer interception missions launched from Earth to investigate arriving ISOs. We pay particular attention to the range of trajectories that could have been used to intercept 'Oumuamua (and find that they could easily have been mounted if the ISO had been detected in time). In Section 5, we enumerate the potential scientific benefits of high-relative-velocity intercept missions and conclude.

'Oumuamua's unusually elongated shape and its lack of a coma were both unexpected, yet these two features may also be connected. The highly elongated figure of the object, coupled with its tumbling (that is, nonprincipal axis rotational) motion, imply that during the multimonth journey through the inner solar system, the surface was repeatedly and randomly illuminated from all directions. Depending on the viewing angle, 'Oumuamua's apparently spindle-like geometry allows the ratio, ξ, of the projected area, A_p(θ, ϕ), to the surface area, S_A, to instantaneously be substantially smaller or larger than the value ξ = 1/4 appropriate to a sphere.

Moreover, as discussed in the previous section, it is reasonable to expect that 'Oumuamua is covered with a porous regolith of material with low thermal conductivity that has been significantly processed during billions of years of exposure to both the galactic cosmic-ray flux and to penetrating ionizing radiation. Given the observed lack of coma and adopting an ellipsoidal three-dimensional shape, we can estimate the combinations of regolith thickness and thermal diffusivity required to protect a frozen, volatile-rich interior from sublimation. Our estimate proceeds by first outlining the calculation of ξ for an adopted body shape and a choice of orientation (or ranges of orientation), and then solving for the rate of heat transfer through the regolith given (1) the solar illumination during the known orbital trajectory, (2) an appropriately time-averaged form factor, ξ, and (3) a parameterized grid of assumptions for the regolith depth and thermal conductivity.

2.1. Geometry of a Projected Ellipsoid

A highly elongated object displays a relatively small projected area relative to its surface area when viewed along its long axis, and a similarly large one when viewed along its short axis. This may have served to limit the solar radiation energy that 'Oumuamua's surface received during its closest approach to the Sun. For quantitative heat transfer calculations, we need an analytic framework to estimate the value of the geometric factor, ξ, as a function of orientation.

In general, a two-dimensional ellipse of the form²

$$\begin{eqnarray}&&{a}_{11}{x}^{2}+2{a}_{12}{xy}+{a}_{22}{y}^{2}+c=0,\end{eqnarray} \tag{ 1 }$$

where x and y denote the spatial coordinates and a₁₁, a₁₂, a₂₂, and c are general coefficients, has an area, A,

$$\begin{eqnarray}&&A=\displaystyle \frac{\pi }{\sqrt{{\lambda }_{+}{\lambda }_{-}}},\end{eqnarray} \tag{ 2 }$$

where λ₊ and λ₋ are the eigenvalues of the general quadratic form for the ellipse equation:

$$\begin{eqnarray}&&{\lambda }_{+,-}=\displaystyle \frac{1}{2c}[-({a}_{11}+{a}_{22})\pm \sqrt{{({a}_{11}-{a}_{22})}^{2}+4{a}_{12}^{2}}].\end{eqnarray} \tag{ 3 }$$

If the tumbling timescale (which is of order of the rotation period) is short in comparison to the encounter timescale, ξ is simply the projected surface area of an ellipsoid, averaged over all viewing angles. Cauchy proved that for any convex surface, this is always $\tfrac{1}{4}$. Simply stated, for the sphere, ξ = 1/4, and because of the symmetry of the sphere, this result holds for each infinitesimal surface area element on the sphere. Therefore, any infinitesimal element of area has ξ = 1/4. Since the surface of any convex body consists of a number of these infinitesimal surface area elements, Cauchy's theorem follows by summing all of these contributions (Meltzer 1949).

If, however, 'Oumuamua's orientation was predominantly head-on or face-on during its periastron passage (when the received flux was changing rapidly), this could serve to significantly increase or decrease the amount of solar radiation that it was exposed to.

Adopting the notation of Stark (1977), who derived a similar quantity for the isophotes of elliptical galaxies, we consider an ellipsoid centered in a three-dimensional coordinate system (x, y, z) with coefficients t, u, and a_v of the form

$$\begin{eqnarray}&&{({tx})}^{2}+{({uy})}^{2}+{z}^{2}={a}_{v}^{2}.\end{eqnarray} \tag{ 4 }$$

A rotation through the first two Euler angles is sufficient to transform to any viewing angle, so consider a rotation of (θ, ϕ):

$$\begin{eqnarray}&&\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=\left[\begin{array}{ccc}\cos \phi & -\cos \theta \sin \phi & \sin \theta \sin \phi \\ \sin \phi & \cos \theta \cos \phi & -\sin \theta \cos \phi \\ 0 & \sin \theta & \cos \theta \end{array}\right]\left[\begin{array}{c}x^{\prime} \\ y^{\prime} \\ z^{\prime} \end{array}\right].\end{eqnarray} \tag{ 5 }$$

Then, in the rotated coordinate system, (x', y', z'), Equation (4) is

$$\begin{eqnarray}&&{{fz}}^{{\prime} 2}+{{gz}}^{{\prime} 2}+h={a}_{v}^{2},\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}f & = & f(\phi ,\theta ,t,u)={t}^{2}{\sin }^{2}\theta {\sin }^{2}\phi \\ & & +{u}^{2}{\sin }^{2}\theta {\cos }^{2}\phi +{\cos }^{2}\theta ,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}g & = & g(x^{\prime} ,y^{\prime} ,\phi ,\theta ,t,u)=\sin \theta \sin 2\phi ({t}^{2}-{u}^{2})x^{\prime} \\ & & +\sin 2\theta (1-{t}^{2}{\sin }^{2}\phi -{u}^{2}{\cos }^{2}\phi )y^{\prime} ,\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}h & = & h(x^{\prime} ,y^{\prime} ,\phi ,\theta ,t,u)=({t}^{2}{\cos }^{2}\phi +{u}^{2}{\sin }^{2}\phi ){x}^{{\prime} 2}\\ & & +\cos \theta \sin 2\phi ({u}^{2}-{t}^{2})x^{\prime} y^{\prime} +({t}^{2}{\cos }^{2}\theta {\sin }^{2}\phi \\ & & +{u}^{2}{\cos }^{2}\theta {\cos }^{2}\phi +{\sin }^{2}\theta ){y}^{{\prime} 2}.\end{eqnarray} \tag{ 9 }$$

Therefore, at a given viewing angle, (θ, ϕ), the two-dimensional projection of the ellipse is defined by the maxima of the rotated coordinates orthogonal to the line of sight along the ellipsoid. For a rotation through the first Euler angle ϕ, with line of sight defined by x', the projected ellipsoid has a semimajor axis along y', defined by the maxima of y', where $\tfrac{\partial y^{\prime} }{\partial x}=0$, and semiminor axis along the z' axis.

The area of this projected ellipse, A_p, is then

$$\begin{eqnarray}&&{A}_{p}=\displaystyle \frac{\pi }{{tu}}{\left(\displaystyle \frac{{u}^{4}{\sin }^{2}\phi +{t}^{4}{\cos }^{2}\phi }{{t}^{2}{\cos }^{2}\phi +{u}^{2}{\sin }^{2}\phi }\right)}^{\tfrac{1}{2}}.\end{eqnarray} \tag{ 10 }$$

For a triaxial ellipsoid of aspect ratio 10:1:1, the value reported by Meech et al. (2017), if viewed head-on, the form factor is minimized with ξ ∼ 0.03, whereas if the body is viewed along the short axis, one finds a maximum value, ξ ∼ 0.31.

2.2. 'Oumuamua's Thermal Response to Solar Heating

Assuming that the timescale for randomization of the principal axis orientation from tumbling, τ_ran ∼ τ_rot, is shorter than the time over which the distance from the Sun to 'Oumuamua changes significantly, we can adopt the approximation ξ = 1/4 that follows from Cauchy's theorem and calculate the solar flux, Φ, that 'Oumuamua received at any point during its trajectory:

$$\begin{eqnarray}&&{\rm{\Phi }}(t)=\displaystyle \frac{{L}_{\odot }}{16\pi d{(t)}^{2}}(1-\alpha ),\end{eqnarray} \tag{ 11 }$$

where d is the instantaneous solar distance, L_⊙ is the solar luminosity, and α is the surface Bond albedo.

To model the evolution of the temperature profile of 'Oumuamua during its encounter with the Sun, we solve the one-dimensional heat conduction equation in cylindrical coordinates

$$\begin{eqnarray}&&\displaystyle \frac{{dT}}{{dt}}=\displaystyle \frac{\kappa }{\rho {c}_{P}}\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left(r\displaystyle \frac{\partial T}{\partial r}\right)\end{eqnarray} \tag{ 12 }$$

for the temperature T, where κ is the thermal conductivity, ρ is the density, c_P is the specific heat, and r is the radial depth.

We solve the above equation on a 200-zone grid to model the entirety of the 1800 cm semiminor axis of the best-fit ellipsoid from Meech et al. (2017). The model has an outer shell of low-conductivity "regolith" and an inner core of H₂O ice. For computational efficiency, we use ice zones that are substantially larger than the regolith zones.

To enforce global energy conservation, we adopt boundary treatments for (1) the outermost regolith surface zone, (2) the regolith–ice transition interface, and (3) the innermost ice zone. For the surface regolith zone, we use a Newton–Raphson technique to find a surface temperature T that balances absorption of solar flux, heating of the surface zone material, reradiation from the surface, and diffusion of heat energy from the surface zone into the interior:

$$\begin{eqnarray}&&{\rm{\Phi }}(t)-\epsilon \sigma {T}^{4}-\kappa {\displaystyle \frac{{dT}}{{dr}}}_{r=0}-{c}_{P}\rho {\rm{\Delta }}r{\rm{\Delta }}T=0.\end{eqnarray} \tag{ 13 }$$

In Equation (13), Φ(t) is the energy flux received in the outer zone, σT⁴ is the energy radiated, κ dT/dr_r=0 is the energy diffused into the adjoining zone, and c_PρΔrΔT is the internal energy change of the zone resulting from the change in temperature. We use  ∼ 0.95 for the bolometric emissivity, and σ is the Stefan–Boltzmann constant. At the ice–regolith interface and the inner boundary, we also use Newton–Raphson iteration to solve similar equations, albeit without the radiative energy terms. For the regolith-to-ice transition condition, we consider both the flux diffused from the innermost regolith zone to the first ice zone, and from the first ice zone to the second ice zone.

Figure 1 shows the results of this thermal modeling applied to 'Oumuamua during a solar encounter. The temperature profile is plotted at depths of 0.0, 3.0, and 6.0 cm and at the center over the 200 days prior to and after periastron. The full simulation was run over a multiyear period spanning 2013 February 09 to 2019 December 09. We use a 200-zone grid representing the asteroid coated in a 6 cm thick layer of porous regolith material. We use values of density ρ_r = 1.0 g cm⁻³ and ρ_ice = 0.92 g cm⁻³, specific heat capacities c_Pr = 0.55 J g⁻¹ K⁻¹ and c_Pice = 2.0 J g⁻¹ K⁻¹, and thermal conductivities κ_r = 9.0 × 10⁻⁵ W K⁻¹ cm⁻¹ and κ_ice = 2.25 × 10⁻² W K⁻¹ cm⁻¹ for the regolith and ice, respectively. We track the energy flux received, radiated, and contained in the simulation, and we chart the reradiation efficiency in the top panel of Figure 1. We require that the total energy is conserved to a factor of better than ∼1% of the sum of the total energy received minus radiated. With this choice of depth and conductivities, a visible coma is absent at all times.

Zoom In Zoom Out Reset image size

It is worthwhile to discuss the apparent differences between our results and those of Jewitt et al. (2017) and Fitzsimmons et al. (2018). Jewitt et al. (2017) make an order-of-magnitude estimation based on the heat equation, t_c ∼ d²ρc_P/κ, to estimate the thermal skin depth d (at which the surface temperature is a factor of ∼1/e times the surface temperature T ≲ 560 K) of the regolith coating, where t_c is the timescale to conduct heat. Using thermal properties similar to those assumed in this simulation and t_c = 2 × 10⁷ s, they estimate that 'Oumuamua has a skin depth d ∼ 50 cm and conclude that, to prevent a visible coma, the depth of the thermal coating must be D ≳ 50 cm. While this provides a first estimate for the coating's thickness, it does not account for the geometric factor ξ from the smaller projected surface area, the efficiency of thermal reradiation over the course of the flyby, and the potentially sharp difference in thermal conductivities at the (assumed discontinuous) ice–regolith interface. We find that with ξ ∼ 0.25 and a radiation efficiency that approaches unity close to and after periastron, the surface temperature heats to T ≲ 550 K. We also see that the presence of a large reservoir of high-conductivity ice beneath the regolith permits a much steeper gradient in regolith temperature due to the ice–regolith boundary condition. Additionally, a slight decrease in temperature diffusion stems from the adoption of cylindrical coordinates.

Fitzsimmons et al. (2018) present high-resolution simulations of the thermal regolith coating down to 5 m in depth. This analysis includes the thermal radiation and uses a geometric factor of ξ = 1/π, a regolith density of ρ_r = 1.0 g cm⁻³, specific heat capacity of c_Pr = 0.55 J g⁻¹ K⁻¹, and thermal conductivity of κ_r = 1.0 × 10⁻⁵ W K⁻¹ cm⁻¹. In their simulations, the surface temperature heats to T ≲ 580 K, and the temperature is reduced by a factor of 1/e at a depth of ∼10–20 cm. Using the order-of-magnitude argument from Jewitt et al. (2017) applied to the thermal properties assumed in the Fitzsimmons et al. (2018) results, one would expect the skin depth to be d ∼ 15 cm, which is consistent with their simulation. The authors note that the skin depth would scale with different thermal properties and albedos than those chosen in their simulation. We verify that our simulation converges to the Fitzsimmons et al. (2018) results with the same thermal and geometric properties and an increased amount of regolith coating.

The major discrepancy between our analysis and that of Fitzsimmons et al. (2018) is the smaller effective skin depth of the regolith, which can be predominantly explained by the additional modeling of the ice. With the thermal properties assumed in Figure 1, the order-of-magnitude estimate of the regolith skin depth is d ∼ 40 cm, which is inconsistent with our simulation. The timescale for heat diffusion over a constant distance scales as t_c ∼ ρc_P/κ. This implies that with the values used by Fitzsimmons et al. (2018), the timescale for heat diffusion of ice is t_cice ∼ t_cr/100. With these choices of thermal properties, heat diffused across the ice–regolith boundary effectively heats up the entire interior of the icy asteroid, maintaining a substantially colder inner regolith temperature, very close to the temperature of the core. With the thermal properties used in Figure 1, the timescale for heat diffusion of ice is t_cice ∼ t_cr/10. With the extra factor of 10 in the diffusion timescale, a temperature gradient is now maintained in the ice zones, evident in the difference in the core and 6 cm temperature profiles in Figure 1, but the first regolith zone is still substantially colder than it would be without the ice. Therefore, equating the skin depth of regolith to the depth of coating at which the icy interior will sublimate is potentially misleading, as the regolith will exhibit a much higher temperature gradient due to the ice–regolith boundary condition. We have verified that the results of these simulations have converged and that the discrepancy is not due to the resolution. In summary, we conclude that the discrepancies can be attributed to the geometric factor, thermal radiation, and the modeling of the icy interior. We note that if ξ > 0.25 or the asteroid was smaller and had less ice to heat up, more regolith coating would be required to prevent development of a coma.

We present similar simulations for a range of depths and thermal conductivities of the regolith coating, and the results are presented in Figure 2. For our grid, we choose a range of thermal conductivities κ_r ∈ (1.0 × 10⁻⁵, 1.0 × 10⁻⁴) W K⁻¹ cm⁻¹ that allow for a factor of 10 difference in the timescale for diffusion across the regolith. At the minimum value of κ_r, the temperature gradient in the ice is very small, leading to a much colder inner regolith zone, while at the maximum κ_r, a temperature gradient is maintained in the ice. We see that, generally, a coma is only produced very close to periastron, with a very small amount of insular coating. We conclude that future ISOs may similarly fail to produce comas in response to solar irradiation. As a consequence, the most definitive way to determine the bulk composition of ISOs may be through excavation of subsurface material using a kinetic impactor.

Zoom In Zoom Out Reset image size

If estimates of high space densities for 'Oumuamua-like objects prove correct, discoveries of similar objects will certainly be forthcoming, and it is useful to assess the expected distributions of arrival directions, speeds, and impact parameters with respect to both the Sun and the Earth.

If we assume that the population of ISOs in the galaxy has a kinematic distribution similar to Population I stars, we can adopt the phase space distribution function for the stellar population in the solar neighborhood (Binney & Tremaine 2008) as a proxy for that of the free-streaming interstellar asteroids. Given that the number density of 'Oumuamua-like interstellar asteroids may be as much as ∼10¹⁶ higher than the stellar number density, the collisionless fluid approximation is even more applicable to this population.

We assume that the interstellar asteroids in the solar system's vicinity follow a Schwarzschild velocity distribution, where the dispersions σ₁, σ₂, and σ₃ come from the Hipparcos catalog (Binney & Tremaine 2008). The coordinate system with respect to a right-handed galactic coordinate is given by the vertex deviation for the primary eigenvector of the ellipsoidal distribution function for a given population of stars. The number density of interstellar asteroids, $n({\rm{\Delta }}{\boldsymbol{v}})$, in a volume of velocity space ${\rm{\Delta }}{\boldsymbol{v}}=({\rm{\Delta }}{v}_{1},{\rm{\Delta }}{v}_{2},{\rm{\Delta }}{v}_{3})$ is

$$\begin{eqnarray}&&n({\rm{\Delta }}{\boldsymbol{v}})=\displaystyle \frac{n}{8}\displaystyle \prod _{i=1}^{3}\mathrm{erf}(\displaystyle \frac{{v}_{i}}{\sqrt{2}{\sigma }_{i}}){| }_{{v}_{i}-{\rm{\Delta }}{v}_{i}}^{{v}_{i}+{\rm{\Delta }}{v}_{i}}.\end{eqnarray} \tag{ 14 }$$

In the vicinity of the Sun, gravitational focusing increases the flux of arriving bodies, especially for periastron distances r < 1 au. To quantify this effect, consider a test particle on a hyperbolic trajectory

$$\begin{eqnarray}&&\displaystyle \frac{1}{r}=\displaystyle \frac{{{GM}}_{\odot }}{{b}^{2}{v}_{\infty }^{2}}(1+e\cos (\theta -\varpi )),\end{eqnarray} \tag{ 15 }$$

where (r, θ) are the polar coordinates with respect to the Sun's position, ${v}_{\infty }$ is the hyperbolic excess velocity, e is the eccentricity, G is the gravitational constant, M_⊙ is a solar mass, ϖ is a phase, and b is the impact parameter. In the absence of the Sun's gravitational field, bodies that enter the 1 au sphere will have b ≤ 1 au. Gravitational focusing augments the effective cross section by a factor ζ = b²/1 au², which is a function of ${v}_{\infty }$ for hyperbolic asteroids. To calculate ζ, we find the maximum $b({v}_{\infty })$ with a periastron distance of 1 au. Extremizing Equation (15) gives the perihelion distance r_p:

$$\begin{eqnarray}&&{r}_{p}=\displaystyle \frac{{b}^{2}{v}_{\infty }^{2}}{{{GM}}_{\odot }(1+e)}.\end{eqnarray} \tag{ 16 }$$

Substituting and canceling give an expression for the impact parameter b. For the case of the solar system, with numerical values for r_p = 1 au, $M=1\,{M}_{\odot }$, the augmentation for the cross section is

$$\begin{eqnarray}&&\zeta \approx 1+{\left(\displaystyle \frac{42\mathrm{km}/s}{{v}_{\infty }+\phi }\right)}^{2},\end{eqnarray} \tag{ 17 }$$

where ϕ is an artificial gravitational softening parameter set to 2 km s⁻¹ in order to avoid the limit of ${v}_{\infty }=0$. We have verified that varying this parameter does not significantly alter the results of this analysis.

We want to connect the velocity distribution in the Sun's galactic neighborhood to incoming ISO trajectories. Every point on the sky isomorphically defines a unit vector with respect to the Sun. Likewise, any one of these unit vectors can be mapped to a three-dimensional galactic velocity vector, unique up to a scaling of the associated speed. We construct a three-dimensional grid of galactic velocities, map each point in the grid to its unit vector and corresponding position on the sky, and calculate the associated probability with Equations (14) and (17). We bin the resulting grid by position to produce a full-sky probability distribution, shown in Figure 3. This map does not present the probability that an ISO will perforate the 1 au sphere at that location in the sky. Rather, it presents the probability that one will approach the solar system with a velocity vector parallel to that direction, with an arbitrary impact parameter.

Zoom In Zoom Out Reset image size

Given a trajectory, the apparent magnitude, m, of an object with absolute magnitude H can be calculated using

$$\begin{eqnarray}&&m=H+2.5{\mathrm{log}}_{10}\left(\displaystyle \frac{{d}_{\mathrm{BS}}^{2}{d}_{\mathrm{BO}}^{2}}{p({\rm{\Theta }}){d}_{\mathrm{OS}}^{4}}\right),\end{eqnarray} \tag{ 18 }$$

where d_BS, d_BO, and d_OS are the distances between the body and Sun, body and observer, and observer and Sun, respectively, Θ is the phase angle defined by

$$\begin{eqnarray}&&\cos {\rm{\Theta }}=\displaystyle \frac{{d}_{\mathrm{BO}}^{2}+{d}_{\mathrm{BS}}^{2}-{d}_{\mathrm{OS}}^{2}}{2{d}_{\mathrm{BO}}{d}_{\mathrm{BS}}},\end{eqnarray} \tag{ 19 }$$

and p(Θ) is the phase integral, approximated by

$$\begin{eqnarray}&&p({\rm{\Theta }})=\displaystyle \frac{2}{3}\left(\left(1-\displaystyle \frac{{\rm{\Theta }}}{\pi }\right)\cos {\rm{\Theta }}+\displaystyle \frac{1}{\pi }\sin {\rm{\Theta }}\right).\end{eqnarray} \tag{ 20 }$$

For given choices of b and ${v}_{\infty }$ of an 'Oumuamua-like asteroid assumed to have H = 22.4, we calculate the two-dimensional trajectory and evaluate the visibility at all points on the 1 au unit sphere, which (in the approximation of a circular orbit) represents every possible location of Earth. On a 100 × 100 grid of b and ${v}_{\infty }$, we calculate the fractional surface area on the sphere for which the object is detectable by the Large-scale Synoptic Survey Telescope (LSST; m ≤ 24) at any point before periastron, and we show contours of the calculated recovery fraction in Figure 4.

Zoom In Zoom Out Reset image size

In order to assess the future detectability and arrival rates of ISOs, we perform two independent Monte Carlo integrations over all possible incoming trajectories. The first integration is computationally efficient because we reduce the problem to a two-dimensional integral via the method described in the previous paragraph. We randomly draw velocities from the Schwarzschild distribution function (without the increased gravitational focusing) and impact parameters from (0, 10) au given by a p(b) ∼ b². Given these two parameter choices, we calculate the percentage of objects that are visible from at least one location on the 1 au unit sphere. We do this for 10,000 draws, and we see that the solution converges to ∼1.6% of incoming ISOs in the Sun-centered 10 au sphere being detectable.

In order to estimate wait times and seasonal variability in detectability, we perform a full six-dimensional Monte Carlo integration. We randomly draw a position in the sky given by the distribution in Figure 3, ${v}_{\infty }$ from the Schwarzschild distribution, b from (0, 5) au given by a p(b) ∼ b² distribution, an angle $\psi \in (0,2\pi )$ representing a direction for b, and a starting epoch placing Earth at a particular point on its orbit. We position the asteroid in space by first moving to the point 5 au along the unit vector corresponding to the point on the sky selected, and then translating within the plane normal to this vector by b in a direction defined by ψ. (Note that in this setup, the maximum starting distance from the Sun is $5\sqrt{2}$ au.) We boost the starting velocity to account for the change in gravitational potential energy at the initial position relative to that at infinity. We then use the N-body package Rebound (Rein & Liu 2012) to integrate both the Earth and the asteroid for 3 years, and we evaluate the magnitude of the ISO using Equations (18)–(20) along its trajectory. With a sample size of one and no observational constraints on the intrinsic magnitude distribution of the galactic population of ISOs, we assume that every ISO has an absolute magnitude of 22.4. We note, however, that due to Equation (18), the results of the simulation scale linearly with the absolute magnitude of the ISO. We perform the integration over 3.65 million trajectories.

This method samples the entire distribution of asteroids in the 10 au cube surrounding the Sun and accurately pinpoints when and where and for how long these objects will be detectable. As a result of our Monte Carlo integration, we find that 2.6% of the ISOs considered during the course of the Monte Carlo simulation will reach m ≤ 24, and 0.017% will reach m ≤ 19. Additionally, we require that the angle between the Sun, Earth, and ISO be greater than 120° to be detectable by LSST. With this constraint, we find that 0.9% of the ISOs will reach m ≤ 24, and 0.013% will reach m ≤ 19. A number density of 0.2 au⁻³ (Do et al. 2018) implies that  ∼200 such objects are in the 10 au cube at any point in time. Adopting the Pan-STARRS's search volume estimates of Do et al. (2018), these statistics can be connected to an expected rate for observers on Earth with a given detection capability. An average crossing time t_cross ∼ 1 year implies an occurrence rate of objects with m ≤ 24 of two per year. The two histograms in Figure 5 show the number of ISOs that reach m ≤ 24 and m ≤ 19 per month, under the constraint that the Sun, Earth, and ISO angle be less than 120°. These thresholds are roughly the limiting magnitudes of LSST and Pan-STARRS, although it is important to note that since these telescopes are not capable of seeing the entire night sky, these numbers should be roughly halved to convert to detections. From the plot, LSST will have very little seasonal variability in detections, while Pan-STARRS is more likely to detect these objects in the spring. This makes sense physically, since at this time, the Earth is moving in the direction of the solar apex and is thus encountering more of these objects.

Zoom In Zoom Out Reset image size

The kinetic energy imparted via a high-impact (v ≳ 40 km s⁻¹) collision with an impinging ISO would excavate a substantial plume of subsurface debris. The content of this material could be examined spectroscopically by a companion flyby spacecraft, permitting an ISO's true composition (and the associated consequences for the planet formation process) to be better assessed. This general strategy, which was adopted by the Deep Impact mission (A'Hearn et al. 2005), requires far less propulsive energy than velocity-matching or sample-return missions and requires only that the ISO trajectory is determined in sufficient time for an interceptor to be marshaled.

To simplify our analysis, we assume that launches onto interception trajectories are mounted with a single impulsive velocity change, ΔV, from a phase space point belonging to a circular, heliocentric, and a = 1 au orbit. Energetically, this task is very similar to the launch of a body at rest (in the corotating frame) from Earth's L1 Lagrange Point (Murray & Dermott 1999). In the calculations below, we assume that Earth's gravity is negligible and that the orbital trajectories of both the ISO and the spacecraft are Keplerian conic sections. This approach amounts to assuming that early stages of the rocket are used to lift the mission stage up out of Earth's gravitational potential well. We have numerically checked that the trajectories of our resulting optimal orbits are not sensitively dependent on the initial conditions, by verifying that the same impulsive thrust sent from L1 and the Earth, under the influence of Earth's gravitational field, results in endpoints that are closer than 0.1 au apart.

With these constraints, we can formulate the problem in the following manner. Given two initial position and velocity vectors at time t₀, ${{\boldsymbol{r}}}_{0}^{{\prime} }$, ${{\boldsymbol{r}}}_{0}$, ${{\boldsymbol{v}}}_{0}^{{\prime} }$, and ${{\boldsymbol{v}}}_{0}$, for an interceptor and target, respectively, we ask for the optimal flight time Δt and the impulsive change to the initial velocity, ${\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} }$, that simultaneously minimize the energy $E=\tfrac{1}{2}| {\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} }{| }^{2}$ while ensuring that a collision occurs.

For the specific case of two elliptical orbits, Leeghim (2013) developed a robust optimization algorithm to the optimal trajectory problem using Lagrange multipliers. Here we generalize this approach for elliptical or hyperbolic trajectories for both target and interceptor.

Following Bate et al. (1971) and Chobotov (1991), it is useful to adopt the so-called universal variable, defined as

$$\begin{eqnarray}\chi =\left\{\begin{array}{ll}\sqrt{a}E & e\lt 1\\ \sqrt{-a}F & e\gt 1\end{array}\right.,\end{eqnarray} \tag{ 21 }$$

given an orbit's semimajor axis, a, and eccentric anomaly, E, or hyperbolic eccentric anomaly, F. In this formalism, the semimajor axis for the hyperbola has a negative value.

Given a value of χ at a time t₀ + Δt, Kepler's equation can be written as

$$\begin{eqnarray}&&\sqrt{\mu }{\rm{\Delta }}t=\displaystyle \frac{{{\boldsymbol{r}}}_{0}\cdot {{\boldsymbol{v}}}_{0}}{\sqrt{\mu }}{\chi }^{2}C(\alpha {\chi }^{2})+(1-\alpha {r}_{0}){\chi }^{3}S(\alpha {\chi }^{2})+{r}_{0}\chi ,\end{eqnarray} \tag{ 22 }$$

where α = 1/a, μ = GM, and C(z) and S(z) are Stumpff functions, defined as

$$\begin{eqnarray}S(z)=\left\{\begin{array}{ll}\displaystyle \frac{\sqrt{z}-\sin \sqrt{z}}{{\sqrt{z}}^{3}} & z\gt 0\\ \displaystyle \frac{\sinh \sqrt{-z}-\sqrt{-z}}{{\sqrt{-z}}^{3}} & z\lt 0\end{array}\right.\end{eqnarray} \tag{ 23 }$$

and

$$\begin{eqnarray}C(z)=\left\{\begin{array}{ll}\displaystyle \frac{1-\cos \sqrt{z}}{\sqrt{z}} & z\gt 0\\ \displaystyle \frac{\cosh \sqrt{-z}-1}{-z} & z\lt 0\end{array}\right..\end{eqnarray} \tag{ 24 }$$

In this formalism, the four Lagrange coefficients,

$$\begin{eqnarray}&&f=1-\displaystyle \frac{{\chi }^{2}}{{r}_{0}}C(\alpha {\chi }^{2}),\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&g={\rm{\Delta }}t-\displaystyle \frac{1}{\sqrt{\mu }}{\chi }^{3}S(\alpha {\chi }^{2}),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\dot{f}=\displaystyle \frac{\sqrt{\mu }}{{{rr}}_{0}}[\alpha {\chi }^{3}S(\alpha {\chi }^{2})-\chi ],\end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray}&&\dot{g}=1-\displaystyle \frac{{\chi }^{2}}{r}C(\alpha {\chi }^{2}),\end{eqnarray} \tag{ 28 }$$

uniquely define the position and velocity vector of an orbit after a time Δt via

$$\begin{eqnarray}&&{\boldsymbol{r}}=f{{\boldsymbol{r}}}_{0}+g{{\boldsymbol{v}}}_{0}\end{eqnarray} \tag{ 29 }$$

and

$$\begin{eqnarray}&&{\boldsymbol{v}}=\dot{f}{{\boldsymbol{r}}}_{0}+\dot{g}{{\boldsymbol{v}}}_{0}.\end{eqnarray} \tag{ 30 }$$

We define two functions, η₁ and η₂, that express the orbital motion of the target and interceptor in terms of their respective universal variables:

$$\begin{eqnarray}{\eta }_{1}(\chi ,{\rm{\Delta }}t) & = & \displaystyle \frac{{{\boldsymbol{r}}}_{0}\cdot {{\boldsymbol{v}}}_{0}}{\sqrt{\mu }}{\chi }^{2}C(\alpha {\chi }^{2})\\ & & +(1-\alpha {r}_{0}){\chi }^{3}S(\alpha {\chi }^{2})+{r}_{0}\chi -\sqrt{\mu }{\rm{\Delta }}t\end{eqnarray} \tag{ 31 }$$

and

$$\begin{eqnarray}&&{\eta }_{2}(\chi ^{\prime} ,{\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} },{\rm{\Delta }}t)=\displaystyle \frac{{{\boldsymbol{r}}}_{0}^{{\prime} }\cdot ({{\boldsymbol{v}}}_{0}^{{\prime} }+{\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} })}{\sqrt{\mu }}{\chi }^{{\prime} 2}C(\alpha {\chi }^{{\prime} 2})\\ &&\quad +(1-\alpha {r}_{0}^{{\prime} }){\chi }^{{\prime} 3}S(\alpha {\chi }^{{\prime} 2})+{r}_{0}^{\prime} \chi ^{\prime} -\sqrt{\mu }{\rm{\Delta }}t.\end{eqnarray} \tag{ 32 }$$

The problem amounts to finding values ${\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} }$ and Δt that satisfy

$$\begin{eqnarray}&&{\boldsymbol{r}}^{\prime} -{\boldsymbol{r}}=0,\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{\eta }_{1}(\chi ,{\rm{\Delta }}t)=0,\end{eqnarray} \tag{ 34 }$$

and

$$\begin{eqnarray}&&{\eta }_{2}(\chi ^{\prime} ,{\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} },{\rm{\Delta }}t)=0,\end{eqnarray} \tag{ 35 }$$

subject to the constraint that

$$\begin{eqnarray}&&J=\displaystyle \frac{1}{2}| {\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} }{| }^{2}\end{eqnarray} \tag{ 36 }$$

is minimized.

Numerical solutions for this system of equations are straightforward to find by the following iterative process:

• 1.  
  Choose a trial flight time Δt.
• 2.  
  Solve for roots of the transcendental Equation (31) for the universal variable of the target χ corresponding to Δt.
• 3.  
  Solve the system of four transcendental equations defined by Equations (32) and (33) for all three components of ${\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} }$ and $\chi ^{\prime}$.
• 4.  
  Repeat process until a Δt that minimizes Equation (36) is found.

For the purpose of this paper, where a reasonable range of possible flight durations are known a priori, the foregoing iterative procedure is sufficient. A more general algorithm is, however, useful, so we present one in the Appendix.

4.1. A Specific Example: 'Oumuamua

'Oumuamua was discovered on 2017 October 19, well after experiencing its periastron passage on 2017 September 9 and, indeed, when its heliocentric distance exceeded r = 1 au. Given the late discovery date and the subsequent span of inaction, it would be extremely difficult to reach. It is interesting, however, to use the above-described procedures to evaluate the mission trajectories that could have been employed over a range of launch dates assuming a sufficiently early detection.

Figure 6 was produced using the procedure described above and charts the optimal values of ΔV(t) (from Earth's orbit, treating Earth as a test particle) as a function of launch dates ranging from 2017 July 7 to 2017 August 30. Each point on the chart is colored by the time of flight for the optimal trajectory.

Zoom In Zoom Out Reset image size

It is clear that, given sufficient warning, it would not have been energetically difficult to intercept 'Oumuamua. Missions launched during the first several weeks of 2017 July require Δv ≈ 4 km s⁻¹. Figure 7 shows the trajectory of the Earth, the interceptor, and 'Oumuamua for the 2017 July 25 minimum-ΔV mission (which had a flight time of 83.38 days). We have verified these matched-conic trajectories with full integrations of the solar system using the Rebound code (Rein & Liu 2012) and assuming departure trajectories from the Earth–Sun L1 point. Moreover, it is evident that, as a first-order rule of thumb, the minimum-energy trajectory will have the date of impact be when the target reaches closest approach to Earth.

Zoom In Zoom Out Reset image size

4.2. Missions to Future Interstellar Asteroids

We now turn to an evaluation of prospects for feasible interception missions to future ISOs, drawing on the integrations that comprise the Monte Carlo simulation described in Section 3. The specific example of 'Oumuamua suggests that optimal interception missions reach impact close to the moment when the ISO either enters or exits the 1 au sphere centered at the Sun. The ΔV required for such an interception mission can roughly be estimated as twice the ΔV required to traverse the distance from the Earth to the asteroid at these dates, given a warning time, defined as the difference in the date of detection and the date of impact. Therefore, we only consider ISO trajectories in the Monte Carlo simulation that (1) are visible using LSST (given both magnitude and angle constraints), (2) have a perihelion distance that is closer than 1 au, and (3) are visible prior to the proposed impact dates. This population is ∼1/8 the total number that are visible by LSST, as shown in Figure 5. Figure 8 shows histograms of the required ΔV for impactor missions prior to and after periastron for each ISO trajectory in this population. The y axis is weighted by the yearly occurrence rates for each ΔV bin. We find that, for the population that satisfies constraints (1) to (3), the date of detection is always prior to the date that the ISO enters the 1 au sphere, and that the median ΔV required for the interception missions is ∼12 km s⁻¹.

Zoom In Zoom Out Reset image size

An order-of-magnitude estimate follows quickly from this estimate of the expected distributions of ΔV for such a mission. The SpaceX Falcon Heavy quotes a payload capability to Mars of 16,800 kg, which we conservatively use for the payload constraint to L1. The Deep Impact mission to Tempel I had an impactor weighing ∼400 kg and a scientific package weighing ∼600 kg (A'Hearn et al. 2005). Due to the uncertainty of the position of the ISO, it seems appropriate to use ∼16 impactors with a total weight of 400 kg. The mission program is greatly assisted by the expected  40 km s⁻¹ velocity of impact with the hyperbolic ISO. Assuming that the remainder of the payload consists of fuel and oxidants, to account for the oxidants and efficiency of the rocket, we allow ∼1200 kg of fuel (with specific energy similar to compressed hydrogen) to produce the ΔV. Equating the kinetic energy to the energy produced by the fuel, we calculate that a maximum ΔV ∼ 15 km s⁻¹ should be attainable, to impart the same amount of kinetic energy (per impact) as the Deep Impact Tempel I interception did. This capability is above the median ∼12 km s⁻¹ of impact missions before periastron in Figure 8. Given the population of favorable ISO trajectories and order-of-magnitude assessment of ΔV capabilities, we predict wait times of order of 10 years between favorable mission opportunities with the detection capabilities of LSST.

'Oumuamua departed as quickly as it arrived, leaving both insights and unanswered questions in its wake.

The developing opinion prior to 'Oumuamua's detection (Moro-Martín et al. 2009; Jura 2011) was that incoming planetesimals with genuinely hyperbolic orbital elements would be a rare occurrence, and upper limits on the number densities from surveys including the Catalina Sky Survey (Drake et al. 2009) and the Pan-STARRS1 projects (Chambers et al. 2016) have become progressively more restrictive (Engelhardt et al. 2017). Moreover, it was expected that when such objects did appear, they would behave like comets. Our analysis suggests that for 'Oumuamua, the combination of a rapid approach, a conductive and extremely cold interior, chaotic tumbling, and an elongated figure provided the subsurface ice with protection from the briefly intense surface heating that transpired during its hyperbolic encounter with the Sun.

It is not yet clear if 'Oumuamua's size and shape are particularly indicative of the usual properties of bodies from its source population, and it seems likely that a resource such as Pan-STARRS will detect kindred objects only every few years. LSST, however, with its 2019 first light, is likely to increase the discovery rate to several per year, based on the work presented here and in Trilling et al. (2017). Additional ISOs, if and when they are found, will set 'Oumuamua into a much more definitive context.

Over the past several decades, targeted missions to comets, asteroids, and outer solar system bodies have provided a raft of important insights into the sequence of events that molded the solar system. By extension, there is significant scientific motivation for investigating interstellar asteroids in close proximity. Given the large velocities involved, both orbit-matching and sample-return missions to a passing ISO are likely to be extremely challenging. Our analysis indicates, however, that when LSST is on line, the wait time between opportunities to intercept an ISO with a kinetic impactor is of order of 10 years or less.

This material is based upon work supported by the National Aeronautics and Space Administration through the NASA Astrobiology Institute under Cooperative Agreement Notice NNH13ZDA017C issued through the Science Mission Directorate. We acknowledge support from the NASA Astrobiology Institute through a cooperative agreement between NASA Ames Research Center and Yale University. D.S. thanks the Gruber Foundation and Patricia Gruber for the Gruber Fellowship and the Stephen B. Butler Scholarship Fund for the Stephen B. Butler Fellowship, both of which also supported the work reported here. We thank Darin Ragozzine, Scott Sandford, Edward Wright, and J.J. Zanazzi for particularly useful conversations and Joel Ong for insightful comments. We also thank the anonymous referee, whose constructive report led to modifications that significantly improved the scientific content of this paper.

Using Lagrange multipliers, Leeghim (2013) demonstrated that for the case of two elliptical orbits, the minimal energy trajectory satisfies

$$\begin{eqnarray}&&h={{\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} }}^{T}{L}_{0}\left[\left(\displaystyle \frac{\partial {\boldsymbol{r}}}{\partial \chi }-\displaystyle \frac{\partial {\boldsymbol{r}}^{\prime} }{\partial \chi ^{\prime} }\right)+\displaystyle \frac{r}{\sqrt{\mu }}({{\boldsymbol{v}}}_{0}-{{\boldsymbol{v}}}_{0}^{{\prime} }-{\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} })\right]=0,\end{eqnarray} \tag{ 37 }$$

where

$$\begin{eqnarray}&&{L}_{0}={\left(\displaystyle \frac{\partial {\boldsymbol{r}}^{\prime} }{\partial {\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} }}-\displaystyle \frac{1}{r^{\prime} }\displaystyle \frac{\partial {\boldsymbol{r}}^{\prime} }{\partial \chi ^{\prime} }\displaystyle \frac{\partial {\eta }_{2}}{\partial {\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} }}\right)}^{-1}\end{eqnarray} \tag{ 38 }$$

is a 3 × 3 matrix, and all partial derivatives for the elliptical case are given in Appendix A of Leeghim (2013).³ Equations (37) and (38) hold in the general case of hyperbolic or elliptical orbits, but the partial derivatives are different if either or both orbits are hyperbolic. The hyperbolic Lagrange coefficients are

$$\begin{eqnarray}&&f=1-\displaystyle \frac{(-a)}{{r}_{0}}[\cosh \left(\displaystyle \frac{\chi }{\sqrt{-a}}\right)-1]\end{eqnarray} \tag{ 39 }$$

and

$$\begin{eqnarray}&&g={\rm{\Delta }}t-\displaystyle \frac{(-a)}{\sqrt{\mu }}[\sqrt{-a}\sinh \left(\displaystyle \frac{\chi }{\sqrt{-a}}\right)-\chi ].\end{eqnarray} \tag{ 40 }$$

In order to solve for $\tfrac{\partial {\boldsymbol{r}}}{\partial \chi }$ and $\tfrac{\partial {\boldsymbol{r}}^{\prime} }{\partial \chi ^{\prime} }$, the partial derivatives in the hyperbolic case are

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial \chi }=-\displaystyle \frac{\sqrt{-a}}{{r}_{0}}\sinh \left(\displaystyle \frac{\chi }{\sqrt{-a}}\right)\end{eqnarray} \tag{ 41 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{\partial g}{\partial \chi }=-\displaystyle \frac{(-a)}{\sqrt{\mu }}(\cosh \left(\displaystyle \frac{\chi }{\sqrt{-a}}\right)-1).\end{eqnarray} \tag{ 42 }$$

To solve for the matrix $\tfrac{\partial {\boldsymbol{r}}^{\prime} }{\partial {\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} }}$,

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial a}=-\displaystyle \frac{1}{{r}_{0}^{\prime} }\left[1-\cosh \left(\displaystyle \frac{\chi ^{\prime} }{\sqrt{-a}}\right)+\displaystyle \frac{\chi ^{\prime} }{2\sqrt{-a}}\sinh \left(\displaystyle \frac{\chi ^{\prime} }{\sqrt{-a}}\right)\right],\end{eqnarray} \tag{ 43 }$$

and

$$\begin{eqnarray}\displaystyle \frac{\partial g}{\partial a} & = & \displaystyle \frac{-1}{\sqrt{\mu }}\left[\chi ^{\prime} -\displaystyle \frac{3\sqrt{-a}}{2}\sinh \left(\displaystyle \frac{\chi ^{\prime} }{\sqrt{-a}}\right)\right.\\ & & +\left.\displaystyle \frac{\chi ^{\prime} }{2}\cosh \left(\displaystyle \frac{\chi ^{\prime} }{\sqrt{-a}}\right)\right].\end{eqnarray} \tag{ 44 }$$

Finally, the vector $\tfrac{\partial {\eta }_{2}^{h}}{\partial {\rm{\Delta }}{v}_{0}^{{\prime} }}$ for the hyperbolic case is

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\eta }_{2}^{h}}{\partial {\rm{\Delta }}{v}_{0}^{{\prime} }}=\displaystyle \frac{\partial {\eta }_{2}^{h}}{\partial a}\displaystyle \frac{\partial a}{\partial {\rm{\Delta }}{v}_{0}^{{\prime} }}+\displaystyle \frac{-a}{\sqrt{\mu }}\left[\cosh \left(\displaystyle \frac{\chi ^{\prime} }{\sqrt{-a}}\right)-1\right]{{{\boldsymbol{r}}}_{0}^{{\prime} }}^{T},\end{eqnarray} \tag{ 45 }$$

where

$$\begin{eqnarray}&&\displaystyle \frac{\partial a}{\partial {\rm{\Delta }}{v}_{0}^{{\prime} }}=\displaystyle \frac{2{a}^{2}}{\mu }{({{\boldsymbol{v}}}_{0}^{{\prime} }+{\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} })}^{T},\end{eqnarray} \tag{ 46 }$$

and the scalar is

$$\begin{eqnarray}\displaystyle \frac{\partial {\eta }_{2}^{h}}{\partial a} & = & \chi ^{\prime} -\displaystyle \frac{3}{2}\sqrt{-a}\sinh \left(\displaystyle \frac{\chi ^{\prime} }{\sqrt{-a}}\right)+\displaystyle \frac{\chi ^{\prime} }{2}\cosh \left(\displaystyle \frac{\chi ^{\prime} }{\sqrt{-a}}\right)\\ & & +\displaystyle \frac{{{\boldsymbol{r}}}_{0}^{{\prime} }\cdot ({{\boldsymbol{v}}}_{0}^{{\prime} }+{\rm{\Delta }}{v}_{0}^{{\prime} })}{\sqrt{\mu }}\left[1-\cosh \left(\displaystyle \frac{\chi ^{\prime} }{\sqrt{-a}}\right)\right.\\ & & \left.+\displaystyle \frac{\chi ^{\prime} }{2\sqrt{-a}}\sinh \left(\displaystyle \frac{\chi ^{\prime} }{\sqrt{-a}}\right)\right]+\displaystyle \frac{{r}_{0}}{2}\left[\displaystyle \frac{-1}{\sqrt{-a}}\sinh \left(\displaystyle \frac{\chi ^{\prime} }{\sqrt{-a}}\right)\right.\\ & & \left.+\displaystyle \frac{\chi ^{\prime} }{-a}\cosh \left(\displaystyle \frac{\chi ^{\prime} }{\sqrt{-a}}\right)\right].\end{eqnarray} \tag{ 47 }$$

Given these definitions, a robust algorithm for finding the optimal trajectory for elliptical or hyperbolic orbits is as follows:

• 1.  
  Choose a trial flight time Δt.
• 2.  
  Solve for roots of the transcendental Equation (31) for the universal variable of the target χ corresponding to Δt.
• 3.  
  Solve the system of four transcendental equations defined by Equations (32) and (33) for all three components of ${\rm{\Delta }}{{\boldsymbol{v}}}_{0}^{{\prime} }$ and $\chi ^{\prime}$.
• 4.  
  Evaluate Equation (37).
• 5.  
  Repeat process until Equation (37) is satisfied.