Authors: Robin Canup, Julien Salmon
Categories: Research Articles, Astronomy, Space Sciences, Research Article, SciAdv r-articles
Source: Science Advances
Phobos and Deimos are likely remnants of a low-mass disk generated by the impact of a Vesta to Ceres-sized body with Mars.
Unlike the Moon, whose mass is about 1% of the mass of the Earth, the combined mass of Phobos and Deimos is only MPD = 2 × 10^−8^
MM, where MM = 6.4 × 10^26^ g is Mars’ mass. How Phobos and Deimos formed remains uncertain. Understanding their origin would provide constraints on Mars’ formation and early dynamical environment, and the moons are the focus of the Japan Aerospace Exploration Agency’s planned Martian Moons eXploration (MMX) mission.
The synchronous orbit at which the orbital period around Mars equals the martian day is located at async ≈ 6 Mars radii (RM) from Mars’ center. Interior (exterior) to this distance, a moon’s orbit spirals inward (outward) because of tides raised by the moon on Mars. Phobos and Deimos’s cratered surfaces suggest that they are ancient objects formed ~4 billion years ago (1). Integrating back in time from Phobos’s current orbit at 2.8 RM implies that it formed at a distance of ~5 to 6 RM, whereas Deimos’s orbit has expanded only slightly over its history to its current orbit at 6.9 RM (2). Thus, both moons likely originated in the region between ~5 and 7 RM, with orbits near but on opposite sides of async. The moons’ reflectance spectra resemble those of primitive asteroids, inspiring the idea that they are asteroids that were captured intact into Mars orbit (3). However, intact capture is difficult to reconcile with the moons’ nearly circular and coplanar orbits, which, instead, suggest that they accreted from an equatorial disk around Mars (3, 4).
A straightforward way to produce a disk is by an impact with Mars, in analogy to a Moon-forming impact with Earth (4–6). Mars’ 25-hour day is consistent with an oblique collision by an impactor of mass Mimp ~ few × 10^−2^
MM (7), whereas Mars’ Borealis basin implies an impactor with Mimp ~10^−3^ to few × 10^−2^
MM (8, 9). Impacts of this scale would produce disks orders of magnitude more massive than Phobos and Deimos (see the Supplementary Materials) (10–12). It has been proposed that Phobos and Deimos accumulated from the outer edge of such a disk, while more massive inner moons spiraled inward and were lost (13, 14). An alternative idea is that Deimos formed at the disk’s edge, but Phobos formed billions of years later from debris left over from tidal disruption of previous inner moons (15). In either case, an impact origin requires an initial disk that extends to ~6 to 7 RM (to account for Deimos’s position and the lack of more distant martian moons), and that tiny moons near async avoided being swept up by transient massive inner moons. Here, we use state-of-the-art models to determine the implied properties of a Phobos-Deimos–forming impact and whether these conditions can be met. We first simulate the accretion of moons from a circum-martian disk using a hybrid N-body model (16, 17) to constrain the disk properties needed to form Phobos-Deimos type systems. Then, we perform hydrodynamical simulations of impacts into Mars to identify those capable of producing disks with these properties.
The Roche limit is located at a distance aR ≈ 2.7 RM for rocky material orbiting Mars. Exterior to aR, collisions lead to the growth of moons, and we track this process by a full N-body simulation (18). Interior to aR, debris undergoes mutual collisions but remains dispersed because of the planet’s differential gravity. Our model describes this inner material as a continuous disk of uniform surface density, whose mass, Min, and outer edge position, rd, evolve via gravitational interactions with outer moons (19) and a collisional viscosity that causes the disk to spread radially (20). Inner disk material is removed when it spreads inward onto the planet or outward past the Roche limit. Material that spreads past aR is added to the N-body portion of the code in the form of new moonlets with mass of ~ 10^−9^MM ~ 0.05 MPD. Outcomes do not strongly depend on this mass so long as it is much smaller than the mass of moons that accrete near the Roche limit (16), as is the case in our simulations (for example, see Fig. 1B). Resonant interactions between outer moons and the inner disk produce a positive torque on moon orbits, which causes them to expand, and a negative torque on the disk, which opposes its tendency to spread outward. Our model considers the strongest disk resonances (the 1, 2, etc.) that produce direct torques only on moons with semi-major axes a < 1.6 rd. However, moons with a > 1.6 rd often experience indirect disk torques due to their interactions with interior moons whose resonances lie in the disk (16, 17). We include inward (outward) migration of moon orbits due to Mars tides for all moons interior (exterior) to async using a constant time-lag model in which the tidal acceleration is proportional to (k2Δt), where k2 is the planet’s Love number (17, 21). The planet’s tidal time lag Δt is related to its tidal dissipation factor Q by Q ~ 1/(2Δt|ω − n|), where 2|ω − n| is the dominant tidal frequency, n is the satellite’s mean motion, and ω is the planet’s angular spin rate. Bodies with nonzero strength disrupt at smaller distances than the classical Roche limit. On the basis of the disruption periapse estimated for low-viscosity bodies in the study of Sridhar and Tremaine (22), we assume that a moon that evolves to within ~2 RM tidally disrupts and we add its mass and angular momentum to the inner disk.
Fig. 1 Accretion of a Phobos-Deimos type system from an impact-generated disk with initial mass M
disk= 10^−5^ MM.The Roche interior disk’s mass and radial extent are indicated by thick horizontal bar; black circles show masses and semi-major axes of simulated moons, with horizontal lines indicating radial variation due to orbital eccentricity. (A) The initial outer disk is described by 1500 particles with a size frequency distribution ∝ d^−3^ (d is particle diameter) and a surface density profile ∝ r^−5^, as suggested by SPH simulation results (see the Supplementary Materials). (B and C) Moons with up to ~10 times the mass of Phobos and Deimos accrete in the mid-region of the outer disk, but strong tidal interaction with Mars causes them to spiral inward and be lost. (D) After 10^7^ years, two small moons with similar properties to those inferred for Phobos and Deimos (red circles) remain on orbits straddling async.
Figure 1 shows the evolution of a disk with initial mass Mdisk = 10^−5^
MM. Strong tidal interaction with the planet is assumed, with (k2Δt) equivalent to martian tidal parameters (Q/k2) ≈ 30 for a moon orbiting at 5 RM. The outer disk first accretes into a distribution of moons in 10^4^ to 10^5^ years. The most massive inner moons remain interior to ~5 RM and spiral inward over 10^5^ to 10^6^ years because of tides. After 10^7^ years, the system has two Phobos-Deimos class objects that accreted from material initially near the disk’s outer edge on low-eccentricity and low-inclination orbits on opposite sides of async. On longer time scales, the inner disk will continue to viscously spread, lose mass, and spawn ever-smaller inner moons (15). Eventually, this material may be lost through a combination of moon inward tidal evolution, tidal disruption, and solar radiation forces (see the Supplementary Materials).
Figure 2 shows the percentage of our disk accretion simulations that leave final satellites orbiting beyond async (solid curve), as well as the total mass of these satellites (symbols), as a function of the initial disk’s mass and the planet’s tidal parameters. Across a range of disk conditions, we find that the survival of small satellite(s) near async (where we define “small” as having a total mass of < 10^−7^MM ~ 5 MPD, per dashed line in Fig. 2) requires Mdisk ≤ 3 × 10^−5^
MM and (Q/k2) < 80 (see the Supplementary Materials). The latter is plausible for early Mars with k2 ~ unity (23).
Fig. 2 Results of disk accretion simulations.Left axis and Total mass of final moons orbiting beyond the synchronous orbit across ~90 accretion simulations as a function of initial disk mass and assumed martian tidal parameters. Survival of small moons (defined as having a total mass < 10^−7^ M
M, which is ~5 MPD, represented by dashed line) requires Mdisk≤ 3 × 10^−5^ MMand (Q/k2) <80. Right axis and solid Percentage of all accretion simulations with final moons orbiting beyond asyncas a function of the initial disk mass. For cases with Mdisk≥ 5 × 10^−5^ MM, either no moons survive beyond asyncor those that do survive are more massive than Phobos-Deimos by a factor of 10 or more (see also the Supplementary Materials). The latter are moons that initially accreted interior to asyncbut were subsequently driven outward beyond asyncby mutual interactions and disk torques.
The finding that there is an upper limit on the initial disk mass and effective (Q/k2) to avoid sweep-up of Phobos and Deimos analogs can be understood by comparing the time scale for orbital expansion of an inner moon of mass m and semimajor axis a due to disk torques, τres ~ MM^2^[(a − rd)/a]^3^ [1.7 a^2^/σΩm]^−1^ (19), where σ is the inner disk’s surface density, to the time scale for the moon’s orbital contraction due to tides, τtides ~ (1/3)(Q/k2)(MM/m)(a/RM)^5^ Ω^−1^, where Ω is orbital frequency. For a given a, the ratio (τtides/τres) is independent of m but varies linearly with (Q/k2)σ. This ratio must be sufficiently small so that tides counteract the tendency for direct and indirect disk torques to drive inner moons outward where they destabilize small bodies near async.
For disk masses >3 × 10^−5^MM, we find no cases that leave small moons near async (see the Supplementary Materials). For these disks, massive moons that accrete from the Roche-exterior disk can have “feeding zones” that extend to the region near async, allowing them to accrete material initially present there. Additional large moons spawned from the Roche-interior disk also produce more significant radial excursions of inner moons than assumed in the previous work by Rosenblatt et al. (14) due to large, scattering-induced eccentricities and/or capture into mean-motion resonances that can drive moons outward beyond 1.6 rd (see the Supplementary Materials).
We next simulate impacts into Mars using smoothed particle hydrodynamics (SPH; see the Supplementary Materials). We find that an Mimp = 0.03 MM impactor, as advocated in previous works (5, 10, 12, 14), produces overly massive disks with 10^−4^ < (Md/MM) ≤ few × 10^−3^ (see the Supplementary Materials). Comparable results are seen in the study of Hyodo et al. (12), in which a 0.03-MM impactor produces Mdisk ~ 9 × 10^−4^
MM, with ~800 MPD orbiting beyond 4 RM. Such a large impactor may have produced Mars’s rotation and a temporary satellite system (for example, see fig. S7), but its massive disk appears inconsistent with forming Phobos and Deimos.
Instead, we find that a much smaller impactor with Mimp ≤ 3 × 10^−3^
MM is needed to produce an appropriately low-mass disk (see Fig. 3 and the Supplementary Materials). Standard SPH uses approximately equal-mass particles. The 10^6^-particle simulation in Fig. 3 (A to C) then describes the disk material with only tens of particles, which is insufficient for determining the disk’s radial mass distribution. To address this, we implement a particle splitting algorithm into our SPH code [see the Supplementary Materials and the study of Kitsionas and Whitworth (24)]. We first use a 10^6^-particle simulation with standard SPH to identify the region on the precollision planet that is ejected above the planet’s surface. We then split each “parent” particle within and neighboring this region into 13 lower-mass “child” particles, whose vector velocities, total mass, and thermal energies are the same as those of the parent. We also split all particles in the precollision impactor. We then repeat the simulation using these objects. The final simulation resolves the impact ejecta with an order-of-magnitude finer mass resolution (Fig. 3, D to F), so that the disk material is described by several hundred particles whose individual masses are comparable to MPD. This provides a factor of ~10 to 10^2^ greater resolution of the disk compared to previous SPH simulations (10–12, 14).
Fig. 3 Simulations of the impact of a Vesta-mass body with Mars, with M
imp= 0.5 × 10^−3^ MM, vimp= 1.5 vesc(7 km s^−1^), and a 45° impact angle.Color scales with temperature in kelvin; distances shown in units of 10^3^ km. Top 10^6^-particle simulation with standard SPH. Bottom Same impact modeled with SPH + particle splitting and an order-of-magnitude higher resolution of the ejected material. After 10 hours, the disk mass is 8.5 × 10^−6^ MM, with ~2 MPDhaving equivalent circular orbits at and beyond ~5 RM, consistent with subsequent accumulation of Phobos and Deimos in the 5- to 7-RMregion. Disk material orbiting beyond the Roche limit is 85% martian in origin and 12% vapor by mass, with temperatures between 1800 and 2000 K.
Figure 4 shows results of SPH simulations with particle splitting. Initially, disk particles have high-eccentricity orbits, leading to high-velocity, disruptive collisions (12) that dissipate energy and circularize debris orbits. Once collisions have sufficiently damped relative velocities, moon accumulation can begin. To estimate the disk’s distribution at the start of moon accretion, we compute the equivalent radius, aeq, of a circular orbit with the same angular momentum as each initial SPH disk particle, that is, with aeq = a(1 − e^2^), where a and e are the initial semi-major axis and eccentricity. The maximum value for aeq provides an estimate of the disk’s radial extent, with aeq,max being the distance beyond which the disk mass is less than the mass of a single SPH particle.
Fig. 4 Properties of disks simulated with SPH + particle splitting of the ejected material.Shape indicates impactor mass, with (M
imp/MM) = 3 × 10^−3^ (circles), 10^−3^ (triangles), and 0.5 × 10^−3^ (squares). Color scales with impact velocity, with blue, yellow, orange, and red corresponding to (vimp/vesc) = 1.1, 1.5, 2, and 3, respectively. Dashed lines indicate upper limit on disk mass needed to preserve small moons near asyncbased on our accretion simulations. Plots show disk mass versus (A) scaled impact parameter (equal to the sin of the impact angle, with a 90° angle corresponding to a grazing impact), (B) percentage of Roche-exterior disk originating from Mars, and (C) aeq,max. The cumulative radial disk mass profiles for disks with masses <3 × 10^−5^ MMare shown in (D).
To produce favorable disk conditions, we find an upper limit of Mimp ~ 3 × 10^−3^
MM because, for this impactor mass, most disks are too massive, and a lower limit of Mimp ~ 0.5 × 10^−3^
MM (for example, see Fig. 3) because for this impactor mass, some disks appear too compact to yield Deimos (see fig. S2 and the Supplementary Materials). This implies that a Phobos-Deimos–forming impactor would be between the mass of Vesta and twice the mass of Ceres.
The successful Fig. 4 impacts occur across a range of impact angles (30° to 60°) and impact speeds (5 to 15 km s^−1^). The successful cases have <aeq,max> = 6.1 ± 0.9 RM, consistent with the disk scale needed to explain Phobos-Deimos. However, we caution that even with particle splitting, the outermost disk “edge” is described by only a few SPH particles, and thus its position remains crudely resolved. Outer disk material (with aeq > aR) from which Phobos and Deimos would accrete is predominantly martian in all cases (77 ± 12%). The latter result differs from canonical lunar-forming impacts, which produce much more massive disks derived primarily from portions of the impactor that are gravitationally torqued into orbit (6). For the much smaller impacts here, disk material is produced by pressure gradients on a portion of the leading, downrange ejecta whose velocities are at low angles relative to the planet’s surface, and this material is largely derived from the target (see the Supplementary Materials).
A Phobos-Deimos–forming impact would have produced one of Mars’s largest basins, which include Borealis, Utopia, and Hellas. Current uncertainties in the impactor size/energy associated with each basin are substantial. Phobos-Deimos–forming impacts in Fig. 4 have energies between 5 × 10^34^ and 2 × 10^36^ erg, with the higher end falling within the range estimated for a putative Borealis-forming collision [3 × 10^35^ erg to 6 × 10^36^ erg; (8, 9)]. Our smallest impactor has a 570-km diameter, within the range of projectile diameters estimated for Utopia (~400 to 800 km) and Hellas (~300 to 700 km) based on two different crater scaling relationships (25). The lower end of the basin-forming impactor size/energy estimates is consistent with Phobos-Deimos originating from the Borealis event, with subsequent Utopia and Hellas impacts producing less massive and more compact disks that did not yield long-lived moons. Initially, eccentric debris from these later events would have bombarded Phobos and Deimos. The higher end of the basin impactor size/energy estimates implies that a Borealis-forming impact would have produced a disk too massive to yield Phobos-Deimos. We find that these disks often produce no long-lived moons because massive moons accrete outer disk material but remain interior to async and are rapidly tidally lost (for example, see fig. S8). In this case, the later Utopia- or Hellas-forming impacts could have produced Phobos-Deimos. Improved models of martian basin formation may allow for these histories to be better constrained.
In a Phobos-Deimos–forming impact, disk material is heated sufficiently to dehydrate OH-bearing minerals (26). In Moon-forming impacts, H2O vapor is gravitationally bound to Earth and may later condense to be at least partially accreted by the Moon (27, 28). However, Mars’s smaller mass implies that H2O vapor at 1800 to 2000 K is vulnerable to rapid escape once ejecta expands beyond ~5 to 6 RM (29), which occurs after only a few hours (Fig. 3). Although detailed modeling is required, this suggests that a Phobos and Deimos formed via impact would have dry endogenic compositions. Because their source materials are predominantly martian in origin, their refractory elemental compositions would be Mars-like.
Phobos and Deimos’s compositions remain uncertain. Reflectance spectra for both moons are similar to those of primitive, D-type asteroids, or alternatively to space-weathered, iron-bearing silicates (30). Phobos’s thermal emission spectra most closely resemble silicates rather than chondritic materials (31). Recent work argues that the moons’ spectra are consistent with surfaces composed of submicron-sized grains that condensed from vapor at temperatures <2200 K in the outer portions of an impact-generated disk (32), potentially in agreement with the conditions found here. The MMX mission will assess the moon compositions through remote characterization of their subsurfaces by a neutron and gamma-ray spectrometer and the eventual return of samples. These data will be crucial to determining whether the moons formed by impact or through an alternative process.
SPH is a Lagrangian hydrodynamic method that describes colliding objects by a large number of spherically symmetric particles, each of which represents a fixed mass of a given composition. The three-dimensional spatial distribution of each particle was defined with a density weighting function, known as the kernel, and a characteristic radius, known as the smoothing length, h. The functional form of the kernel did not change during a simulation, but the smoothing length of each particle was varied so as to maintain overlap with a desired number of other particles (typically a few tens), which allowed low-density regions to be smoothly resolved, although with coarse spatial resolution. In the version of SPH used here [a descendant of that of W. Benz, as in the study of Canup (6)], the evolution of each particle’s position, velocity, internal energy, and density were evolved because of gravity, pressure forces, and shock dissipation. Material strength was neglected. A tree code was used for the gravity and nearest neighbor calculations, and the code was parallelized.
The equation of state relates a particle’s specific internal energy and local density to pressure at each time step. Our simulations used the semianalytic equation of state known as ANEOS (33, 34). We considered dunite/forsterite mantles [with equation of state parameters provided in table S1 of the study of Canup (35)] and iron ANEOS cores, as in the study of Canup (6). We considered a differentiated target Mars containing 30% iron and 70% dunite by mass. We considered impactors that are either differentiated (with the same iron-to-mantle mass ratio as the planet) or undifferentiated (100% dunite).
Each simulation was continued for ~10 hours by which time the short-term effects of the collision on the planet’s shape and any substantial secondary impacts had ended (for example, see Fig. 2). We used an iterative procedure (6) to determine whether material is in the planet, in bound orbit around the planet, or escaping. For each bound particle that is outside the planet, we computed an equivalent circular orbit semimajor axis, aeq, defined by setting GMMaeq equal to the particle’s specific angular momentum normal to the equatorial plane of the planet. The equivalent circular orbit is representative of that to which the mass represented by a particle would settle after undergoing mutual collisions, which rapidly damp orbital eccentricities and inclinations but transport angular momentum much more slowly. Initially, ejecta may disperse over an expanding volume, but mutual collisions would become increasingly probable as bound material reapproaches the planet, particularly for small expected ejecta sizes (12). We defined those particles with aeq greater than the equatorial radius of the planet as being “in the disk” and those that are energetically unbound as escaping. At 10 hours, the calculation of the disk mass by this method is within 10 to 40% of that calculated by summing over the SPH particles whose instantaneous periapses are above the planet’s surface.
Standard SPH simulations do not provide sufficient resolution to accurately determine the properties of a low-mass disk with 10^−5^
MM~, because the disk is represented by only of order 10 SPH particles even when 10^6^ particles are used to describe the colliding objects. We designed a particle splitting technique in which we “split” a selection of particles to increase resolution of the disk by an order of magnitude. Details on this technique are given in the Supplementary Materials.
We thank two reviewers for their helpful comments and suggestions. Funding: R.C. and J.S. acknowledge support from NASA’s Emerging Worlds and Solar System Exploration Research Virtual Institute programs and Southwest Research Institute for software development. **Author ** R.C. designed the study and wrote the paper. R.C. was responsible for the modeling, design, and data processing of the SPH simulations. J.S. was responsible for the modeling, design, and data processing of the disk accretion simulations. **Competing ** The authors declare that they have no competing interests. **Data and materials ** All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.