On the Impact Origin of Phobos and Deimos. II. True Polar Wander and Disk Evolution

The origin of the two small Martian moons, Phobos and Deimos, is still debated. It has been believed that they were captured asteroids, due to their spectral properties resembling those of D-type asteroids (e.g., Burns 1978; Murchie et al. 1991). However, the captured scenario has been confronted with the challenge of explaining their almost circular equatorial orbits around Mars (inclinations of 108 and 179 from Mars’ equator for Phobos and Deimos, respectively). In contrast, accretion within debris disks formed by a giant impact may naturally explain their orbital configurations (Craddock 1994, 2011; Ida et al. 1997; Hyodo et al. 2015; Hyodo & Ohtsuki 2015), but Crida & Charnoz (2012) showed that Phobos and Deimos cannot be formed directly from the spreading of a ring interior to the Roche radius. Then, Rosenblatt et al. (2016) have shown that these small moons can be formed by accretion within a thin debris disk that extends outside the Roche limit ( where is the Mars’ radius) and is sculpted by an outward migration of a large inner moon that is formed by the spreading of a thick disk lying inside the Roche limit. Furthermore, Hesselbrock & Minton (2017) have shown that the tidal disruption of such a large inner moon during the tidal decay creates a new generation of rings/disks around Mars followed by the spreading and accretion of smaller moons. They showed that such a process could have occurred repeatedly over the past 4.3 billion years, suggesting that Phobos is the last generation of moon we observe today. Note that currently the only successful scenario to form Deimos is the accretion within an extended outer disk proposed by Rosenblatt et al. (2016).

Such a Martian-moon-forming disk can be created by a giant impact that can also form the asymmetric northern lowland that is the Borealis basin: impactor mass of (about 1/3 Martian radius) and an impact velocity of ∼6 km s⁻¹ (Marinova et al. 2008; Citron et al. 2015; Rosenblatt et al. 2016; Hyodo et al. 2017b). However, two natural questions arise. First, without or a slow pre-impact spin of Mars (compared with the impact spin angular momentum), a giant impact spins up Mars, thus a debris disk is generated symmetrically around the equatorial plane. In this case, the impact point (that is the Borealis basin) is expected to be located near the equator, so why is the Borealis basin currently located on the northern hemisphere and not on the equatorial plane? Second, during protoplanet formation through successive accretion of planetesimals, a protoplanet may naturally have a rotation (Ohtsuki & Ida 1998). Thus, if Martian pre-impact spin is comparable and not aligned to the angular momentum delivered by the giant impact, the resultant debris disk is expected to be inclined with respect to the equatorial plane of Mars. Also, the location of the Borealis basin may not be in the polar region. So, why do Phobos and Deimos orbit almost on the Martian equatorial plane, and why is the Borealis basin currently located on the northern hemisphere? In addition, just after the giant impact, the orbits of the debris are eccentric (Section 3.1). Thus, we need a dynamical path to form a circular equatorial disk from which Phobos and Deimos can accrete (Rosenblatt et al. 2016).

In Section 2, using analytical arguments, we show that the Borealis basin-induced reorientation of the planet (true polar wander (TPW)) can explain the dynamical path of the Borealis basin settling in its current location in the northern hemisphere. In Section 3, we detail the dynamical evolution of the initial disk that is formed by a Borealis basin-forming impact and show that an initial “inclined” (with respect to Mars’ equatorial plane) and eccentric disk will settle into a thin equatorial circular disk, forming equatorial orbiting Phobos and Deimos. In Section 4, we summarize our results.

The Borealis impact basin is expected to be located near the equator of the planet, if one assumes that the giant collision that caused it gave Mars most of its present spin. Even if Mars had a significant pre-impact spin, it is unlikely that the Borealis basin could form directly beneath the polar region. Under this assumption, one needs to explain the current polar position of the Borealis basin. Here, we check for TPW as a mechanism for providing the required motion from the equator to polar areas.

We use the equilibrium theory proposed by Gold (1955) and adapted by Willemann (1984) and Matsuyama et al. (2006) for assessing the Mars TPW induced by the Tharsis bulge. This theory consists of computing the effect of surface mass excess load and the effect of the rotational bulge excess mass on the inertia tensor of the planet. The basic calculations of this theory are the diagonalization of the perturbed inertia tensor to provide the position of the new pole at the surface of the planet after TPW motion (see also Matsuyama et al. 2006). In this theory, the latitudinal shift δ of the pole at the surface of the planet is given by the following relationship:

where is the co-latitude of the center of the axisymmetric excess mass before TPW, is the ratio of the mass excess over the rotational bulge load, and α is a combination of load and tidal Love numbers of the planet (see the Appendix). Although this theory has been applied for Mars with surface mass excess (Tharsis), it can also be applied to mass deficit such as that created by impact basins (Willemann 1984). In that case, the center of the impact basin is expected to move from its initial location toward the pole of the hemisphere. Thus, it is likely that the center of the Borealis was initially located in the northern hemisphere in order to move toward its current location. Then, has a negative value and is about −1.102 for the Borealis basin, depending on the tidal Love numbers of Mars, i.e., on the lithospheric thickness (see the Appendix). Figure 1 displays the expected position of the center of the Borealis basin after TPW (final position) versus its initial position.

Zoom In Zoom Out Reset image size

However, we need to know the Borealis position after Borealis-TPW to find which pre-Borealis-TPW positions are suitable from Figure 1. As Tharsis TPW is the most recent large TPW event in Mars history, we need to compute the position of Borealis before Tharsis TPW. The Tharsis-TPW displacement is toward the equator with amplitude of 189 (see Matsuyama & Manga 2010; Bouley et al. 2016, for most recent estimations), meaning that Tharsis formed at higher latitudes by about 20° than it is today. This Tharsis-TPW has also moved the Borealis position. However, the present longitude is 208°E (Andrews-Hanna et al. 2008), which is off the Tharsis central meridian (259.5°E, Matsuyama & Manga 2010) by about 50°. Simple spherical trigonometry considerations show that the TPW amplitude decreases when the distance to its central meridian increases. As the difference between Borealis and Tharsis central meridians is about 50°, the Borealis center will not move by 189 in latitude but only by 5°. It corresponds to a pre-Tharsis TPW co-latitude of 18° (latitude of 72°) for Borealis center (see the Appendix).

If one assumes that no significant TPW events occur between Borealis and Tharsis-TPW ones, this 18° co-latitude of Borealis center before the Tharsis TPW in turn suggests two possible initial (pre-Borealis) latitude ranges between 45°–50° and 5°–10° (co-latitudes of 45°–50° and of 85°–80°, respectively) for Mars’ lithospheric thickness at time of TPW between 50 and 200 km (see Figure 1). Intermediate latitudes, between 10°N and 45°N are possible for lithosphere thinner than 50 km. The 5°–10° range of initial latitudes indicates large TPW displacement between −62° and −67° (three times larger than Tharsis-TPW), so a Borealis initial position near the equator in agreement with a collision having given Mars most if its spin. The 45°–50° range implies a more modest TPW between −22° and −27° (see Equation (1)) comparable to the absolute value of the Tharsis-TPW. The 45°–50° range also implies that the giant collision would not have given to Mars most of its spin, meaning in turn that Mars would have significant pre-impact spin. This spin is, however, difficult to compute because it also depends on the angle and the kinetic energy of the collision.

The linear theory of Matsuyama & Manga (2010) assumes elastic rheology for the lithosphere and the planet and so for instantaneous TPW displacement. Chan et al. (2014) have introduced a more realistic rheology allowing for assessing the timescale of the TPW displacement. They found a typical timescale of 20 Myr for the Tharsis TPW to reach the displacement predicted by the equilibrium theory (see Figure 4 in their paper), depending on the viscosity of the Martian mantle. At the time of Borealis formation, the viscosity of the mantle is expected to be lower than at time of Tharsis formation, hence implying an even shorter timescale for the Borealis-TPW displacement to take place. Thus, the Borealis center can easily reach its final position before the Tharsis-TPW takes place.

Although the two solutions for Borealis-TPW are mathematically possible, the solution with lower TPW displacement seems more plausible than the larger one. Large TPW is indeed expected to produce planetary scale stresses in the lithosphere, which should be seen as a global tectonic pattern. But such a global pattern is not observed at Mars surface today (see Bouley et al. 2016). Although these authors argue that the lithospheric stresses should be relaxed with time, the tectonic pattern should be still recorded in the lithosphere as the secular cooling of the planet stiffens the lithosphere with time. Actually, even modest TPW generates stresses in the lithosphere. However, as these stresses are more modest, they may less likely produce tectonic pattern than in the case of huge TPW. Therefore, a more modest TPW is more in agreement with the absence of TPW-induced global tectonic pattern. Such an issue deserves, however, more thorough investigations that will be performed in future works. In turn, the dichotomy might not be initially lying at the pre-Tharsis equator, in contrary to the hypothesis in Bouley et al. (2016). Nevertheless, large impacts occurring between Borealis and Tharsis events might have produced TPW of a couple of degrees. Considering a few large impacts (Bottke & Andrews-Hanna 2017), we indeed found that Utopia and Hellas could have produced such TPW, which, in turn, could have modified the position of Borealis center after its own TPW with regard to the position derived from the Tharsis-TPW alone (the estimation of the initial position of Borealis). Hence, thorough investigations are needed to fully draw Mars TPW history to check whether our Borealis TPW calculations are still in agreement with Tharsis TPW predictions within the error bar of the Tharsis center determination (Matsuyama & Manga 2010) and the error bar of the estimation of Tharsis TPW displacement from surface observations (e.g., Perron et al. 2007; Kite et al. 2009). This will be performed in future works.

The application of TPW linear theory suggests two possible solutions for Borealis-TPW: on one hand, large TPW and so near-equatorial position for Borealis impact center, and on the other hand, smaller TPW with latitude impact center. The latter solution is, however, favored for physical considerations, and in turn implies a non-equatorial debris disk. The evolution of this non-equatorial disk is studied in the next section.

The results of TPW in Section 2 imply that the giant impact is more likely to occur not on the equatorial plane, but at higher latitude. This indicates that the mean inclination of the debris just after the impact is expected to be inclined with respect to the equator, depending on the pre-impact spin state. In addition, the orbits of the initial debris are expected to be eccentric (see Section 3.1). However, to form Phobos and Deimos in the equatorial plane in a framework of the giant impact hypothesis, a dense inner equatorial and a light outer equatorial disk are required to form (Rosenblatt et al. 2016). Inelastic collision can decrease eccentricity through energy damping, but it conserves the angular momentum. The inner dense and outer light surface density profile of the Martian-moon-forming disk (Rosenblatt et al. 2016, see their Figure 1) was generated as a consequence of the angular momentum conservation from the same a − e distribution of the debris as that shown in Figure 2 by calculating an equivalent circular orbit of radius for every disk particles () as done also in the case of the Moon-forming impact and its disk (e.g., Canup 2004). However, the detailed collisional evolution and its timescale are unclear. In this section, we address a dynamical path that can potentially bring the impact-generated inclined and eccentric debris to the circular equatorial disk that can form Phobos and Deimos.

Zoom In Zoom Out Reset image size

3.1. Orbital Elements of Disk Particles Just After the Impact

Figure 2 shows the distribution of eccentricity of disk particles just after the impact, which is taken from the Martian-moon-forming impact simulation in Hyodo et al. (2017b). We find that the disk particles initially have a wide distribution of large eccentricity (up to almost one) around Mars. The distribution of semimajor axes and eccentricities of disk particles just after the impact can be analytically derived as follows. For particles that end up orbiting around Mars (called “disk particles”), the near impact point is expected to be their pericenters (because Keplerian orbits are closed trajectories, and if their pericenter is inside Mars, they would be accreted very rapidly). Therefore, pericenter distances of the disk particles can be written as where is the small distance from the surface of Mars, a is the semimajor axis and e is the eccentricity. Thus, we can derive the relationship between a and e as

Thus, the distribution of a and e of disk particles just after the impact is expected to be along the line obtained by Equation (2) (Figure 2), and this means that all particles share almost the same pericenter distances.

3.2. Orbital Evolution of Disk Particles

Just after the impact, the disk particles are almost aligned (phase alignment) at their longitude of the pericenters (Figure 3) with wide distributions of eccentricity and semimajor axis (Figure 2). In addition, the disk may be initially inclined from the Martian equatorial plane depending on the giant impact condition. Thus, under these circumstances, collision velocities between nearby particles are significantly small about their shear motion, thus collisional damping is not effective unless collision happens at pericenter between particles whose orbital elements are significantly different. In this subsection, we will discuss two extreme cases of the expected dynamical evolutions of the system where the system forms a torus-like structure and where the system damps very quickly to form a thin inclined nearly circular disk. Note that direct N-body simulations including fragmentation is necessary to understand more details about the disk evolution. However, this is beyond the scope of this work and we will leave this matter in the future work.

Zoom In Zoom Out Reset image size

3.2.1. The Case for Forming a Torus-like Structure

Here, we will discuss the case when particles experience orbital precession and the system is randomized before particle–particle collisions are effective. As was also discussed in Hyodo et al. (2017b), just after the impact, the particles have large eccentricities when they start to be influenced by Martian oblateness (mainly J₂ term), i.e., as their orbits (the argument of pericenter ω and the longitude of ascending node Ω) start to precess around Mars. Using value for Mars (http://pds-geosciences.wustl.edu/mro/mro-m-rss-5-sdp-v1/mrors_1xxx/document/shadr.pdf), we can calculate the timescales of precession as and , where and , respectively (Figure 4). Note that the actual value of J₂ may differ from the current value, and it is important to calculate the precession rate, but this is the beyond the scope of our paper. The timescales depend on the eccentricity, semimajor axis, and inclination from the equatorial plane of Mars but for parameters of interest here ( and ), timescales range from 1 to 100 years depending on the inclination of the disk with respect to the Martian equatorial plane (Figure 4). The ratio between these two precession timescales is that only depends on i, and these timescales are comparable for the nominal case of i = 45°.

Zoom In Zoom Out Reset image size

After the formation of a torus-like structure (see also Hyodo et al. 2017a, 2017b), particle–particle collision may collapse the system into equatorial plane. Here, we analytically estimate the collisional timescale between particles after they form a torus-like structure. In the particle-in-a-box approximation, the collision timescale can be written as

where is the number density of particles, is the collision cross section and is the relative velocity. As also discussed in Hyodo et al. (2017a), after the system forms a torus-like structure, orbits of particles are randomly oriented so that they cross each other and their relative velocity becomes about their local Keplerian velocity. Without gravitational focusing, the cross section can be written as , where is the particle size. Here, we can neglect the gravitational focusing term because when particles are close to or within the Roche limit, gravitational attraction between particles becomes negligible (Ohtsuki 1993; Hyodo & Ohtsuki 2014, 2015). Following the argument of Hyodo et al. (2017a), the volume number density of particles can be written as

where is the number of particles as a function of a, e, and i. and are the probability of finding a particle at radial distance r and an angle ψ from the equatorial plane of the planet and are written as

and

respectively.

Here, to estimate the collision timescale, we first integrate over a and e assuming the disk mean inclination of either i = 10°, 45° and 80°. We use the data obtained from SPH simulations from Hyodo et al. (2017b) (Figure 2) to know the distribution of a. Then, the distribution of e is obtained using Equation (2). Hyodo et al. (2017b) has considered that the particle size that formed during the impact ejection is regulated by their local shear velocity and material surface tension, and they estimated that the typical particle size is m. We assume particle density of m⁻³ to calculate the number of particles as , where is the mass of one SPH particle whose orbital elements are a and e, and is the mass of particle. The collision timescale decreases when r is close to the particles’ pericenter and/or ψ is closer to i. This is because near pericenter, apocenter, and at maximum elevation, the radial and vertical orbital velocity is zero so that the particle residence time is at its maximum there. Thus, in this study, we assume that collision takes place at and . Then, assuming the relative velocity of ) for orbital elements of e = 0.8, we estimate the collision timescale of less than year for i = 10°, 45°, and 80°. Note that in the above argument, we assume arbitrary choices of the location of collision and relative velocity, but we confirm that the collision timescale is always much less than a year, even with any other choices. The estimated collision timescales close to the pericenter are much smaller than the orbital period of particles. In contrast, the collision timescales at around apocenter distance are larger than the orbital period of particles.

Thus, as soon as the disk particles come back to their pericenter distances, they experience collisions. Therefore, after the formation of a torus-like structure, the system is expected to collapse through inelastic collisions to form a thin equatorial disk on a timescale comparable to their orbital period. In addition, the eccentric orbits of the debris are circularized at the same time when the inclination is damped, forming the inner dense and outer light radial profile of the disk in Rosenblatt et al. (2016) under the assumption that every particle converges to the circular orbit corresponding to its angular momentum. However, we may need more detailed investigations that consider collisional fragmentation. We leave this matter for future works.

3.2.2. The Case for Forming a Flat Inclined Low-eccentricity Disk

Even under the phase alignment, particles may experience high velocity collision at their pericenters because they share almost the same pericenter distances and have a wide distribution of large eccentricity (see Figures 2 and 3). In this subsection, we discuss the case when collisional damping at their pericenter is very efficient, thus a flat inclined low-eccentricity disk is formed. If collision occurs between particles that share the same pericenter distance but have different orbital elements ( where ) under the phase alignment as seen in our case (see Figures 2 and 3), collision velocity can be estimated as

Thus, if collision occurs between particles whose eccentricities are e = 0.9 and e = 0, respectively, we get km s⁻¹. Then, their initial collision timescale is expected to be their synodic period when there is no dense circular inner rings as

where is the difference in semimajor axis of two colliding particles. In contrast, after dense circular rings are formed due to collision damping, the collision timescale between the dense circular rings and eccentric particles becomes about their Keplerian period (∼few days for particles whose ). Such high collision velocity may quickly damp the system and be energetic enough to form sized particles (Hyodo et al. 2017b). As also discussed in Section 3.2.1, the initial a − e distribution is expected to be circularized through inelastic collisions while conserving the angular momentum of the disk, eventually forming the inner dense and outer light radial distribution of Rosenblatt et al. (2016).

The timescale of the circularization (∼tens of days; see Equation (8)) is much shorter than the precession timescale ( years, see Section 3.2.1). Therefore, the system forms an inclined low-eccentricity inner dense and outer light disk, while the disk keeps its average inclination from the equatorial plane of Mars resulted directly from the giant impact rather than the nodal precession forms a torus-like structure with large eccentricity.

After the formation of a flat inclined circular disk, neglecting the effect of self gravity (we will discuss this effect in the next paragraph), the evolution is expected to be either of the following two extreme cases. Case (a): if collision timescale is shorter than the differential precession timescale to form a torus-like structure, inelastic collision occurs between nearby particles whose longitudes of node are slightly different as a result of the differential nodal precession and their mean inclination decreases. This is because the nodal differential precession induced by J₂ term always tries to make the system symmetric to the equatorial plane of Mars. Thus, the disk takes gradual inside-out evolution by lowering its mean inclination to settle into the Martian equatorial plane with differential precession timescale (see also Figure 4 at small e). The other extreme case is case (b): if collision timescale is larger than the precession timescale, the system first forms a torus-like structure with the differential precession timescale. Then, collision will collapse the system into the equatorial plane. In our case, the collision timescale in the outer disk (in the case of a torus structure with its inclination i = 45° and with the radial width of where Phobos and Deimos are expected to form with a disk mass of , and particle size of 1 m and its density 2500 kg m⁻³) is estimated by using the same argument above (Equation (3)), where , , and m³. Thus, , m² and m s⁻¹. Thus, years. This timescale is actually the maximum timescale of collision because we assume the scale height of and that the actual particle size is expected to be smaller due to collisional cascade. So, collision timescale is much smaller than the precession timescale. In the inner disk, the collision timescale is shorter than that of the outer disk due to its larger number density and larger relative velocity. Thus, in reality, case (a) occurs.

Lastly, using N-body simulations that include self gravity between all particles, we investigate the evolution of the system by nodal precession. Under some condition, the self gravity of the disk can prevent differential precession, and the disk precesses rigidly without changing its inclination (Batygin 2012; Morbidelli et al. 2012). Our N-body code is the same as that used in Hyodo et al. (2015, 2017a) and Hyodo & Charnoz (2017), and we include the effect of J₂ (see more details Hyodo et al. 2017a). In the case of the canonical Martian-moon-forming disk, the inner disk () is massive ( kg) and the outer disk () has a mass of Phobos ( kg) (Rosenblatt et al. 2016; Hyodo et al. 2017b). Orbits of the inner disk are expected to precess much quicker than those of the outer disk (we have confirmed by N-body simulations that the self gravity within the inner disk does not prevent the differential precession and the case (a) discussed above is expected to occur on a timescale of ∼few 10 years (see Figure 4 at e = 0)). Thus, the inner part may quickly form a thin equatorial massive inner disk, while the outer disk remains inclined. The orbital period of particles within the inner massive disk is much shorter than the precession timescale of the outer disk (Figure 4). Thus, we can approximate the effect of the massive inner disk on the nodal precession of the outer disk as a secular perturbation of the inner massive satellite (see also Morbidelli et al. 2012). In addition, the inner massive disk quickly spreads and massive inner moons are formed (Crida & Charnoz 2012; Rosenblatt et al. 2016). Such inner massive disk or moon(s) increase the net effect of J₂ and the nodal precession can be accelerated. Thus, following the procedure in Morbidelli et al. (2012), we perform N-body simulation including a single massive satellite (mass of 10²⁰ kg) on circular orbit () on the equatorial plane of Mars. In contrast, the outer disk is represented by a swarm of 1000 equal-mass particles initially on circular orbits with their mean inclinations of 45° from the equatorial plane of Mars. In addition, we investigate the outer disk evolution without the inner massive satellite (which means without the effect of the inner disk).

Figure 5 shows time evolutions of the longitude of nodes of disk particles. The inner satellite enhances the net effect of J₂ and the differential precession is slightly accelerated compared with the case of no inner massive disk (see also Morbidelli et al. 2012). Our N-body simulations show that in the case of Phobos and Deimos forming disk ( kg), the self gravity is not effective to induce the rigid precession of the disk. In the framework of the giant impact hypothesis, the inner large moon is formed from the Roche-interior massive disk on a timescale of ∼100 years and migrates outward up to on a timescale of ∼1000 years (Rosenblatt et al. 2016). As discussed above, in the outer disk, the collision timescale is shorter than the nodal precession timescale, thus the disk inclination decreases with the precession timescale (case of (a)). Together with the results of N-body simulations (Figure 5), the nodal precession timescale in the outer disk is a few tens to a few hundreds years (see Figure 4, , and ), depending on the initial inclination. Therefore, it is likely that the outer disk eventually settles into a thin equatorial near the circular disk before the large inner moon migrates outward so that the equatorial Phobos and Deimos can accrete (Rosenblatt et al. 2016).

Zoom In Zoom Out Reset image size

The origin of Martian moons Phobos and Deimos is intensely debated. Recent works have shown that they might have accreted within a debris disk produced by the giant impact that formed the Borealis basin (Rosenblatt et al. 2016; Hesselbrock & Minton 2017). If so, two dynamical questions naturally arise. First, why is the Borealis basin not on the Martian equatorial plane, but close to the north pole? The Borealis basin-forming impact can produce almost all of the current spin period of Mars. Thus, if there is no pre-impact Martian spin or a slow spin, the Borealis basin is expected to form around the equatorial plane and not on the northern hemisphere. Second, why do Phobos and Deimos orbit almost on the equatorial plane of Mars? If there is a pre-impact spin on Mars comparable to that given by the Borealis basin-forming impact, and if the pre-impact spin axis is not aligned to the angular momentum vector of the impactor, the resultant disk orbital plane is expected to be different from the equatorial plane of Mars, thus non-equatorial Phobos and Deimos may form.

To answer the first question—why the Borealis basin is located close to the northern pole and not near the equatorial plane—we investigated the planetary reorientation, due to the mass deficit at the Borealis basin (Section 2). We found that Borealis-induced TPW, using the equilibrium theory, can provide the required reorientation of the planet to move the center of the Borealis to its current latitudinal position from its initial location, which is between 5°N and 50°N (more likely between 45°N and 50°N), for a lithospheric thickness between 50 and 200 km (intermediate latitudes, between 5°N and 45°N are possible, but for lithospheric thicknesses thinner than 50 km). However, this estimation could be refined by taking into account the post-Borealis TPW expected with large impact basins like Hellas or Utopia and with Tharsis, but this is out of the scope of this paper.

Based on the canonical Martian-moon-forming impact (an impact energy of J and an impact angle of 45° (Rosenblatt et al. 2016; Hyodo et al. 2017b)) that also formed the Borealis basin (Marinova et al. 2008), and using analytical arguments, we investigated the detailed post-impact disk. Just after the impact, disk particles have large eccentricities (Section 3.1) and almost the same pericenter distances. In addition, their orbits are almost aligned (“phase alignment”) at their longitude of the pericenters (Section 3.2) and the initial disk is expected to be inclined with respect to the equatorial plane of Mars. Thus, collisional damping is not efficient because collision velocity is only the order of their shear velocities except their pericenter distances, where the maximum collision velocity is ∼1 km s⁻¹. In this paper, we considered that the debris disk may experience either of the following two dynamical paths before forming a thin circular equatorial disk: formation of a torus-like structure (Section 3.2.1) or formation of a thin inclined (with respect to Mars’ equatorial plane) circular disk (Section 3.2.2). Comparing the timescales of these two dynamical evolutions, we found that the formation of a thin inclined circular disk is expected to occur preferentially, due to fast collisional damping at particle’s orbit pericenter. Then, due to the differential nodal precession and particle–particle inelastic collisions, the inclined disk is expected to experience an inside-out evolution to gradually lower the mean inclination of the disk to eventually settle into the equatorial plane within 1000 years (Section 3.2.2). Then, the thin equatorial circular disk is expected to form Phobos and Deimos near the equatorial plane (Rosenblatt et al. 2016). Thus, the above arguments are likely to be a dynamical pathway that can answer our second question: why do Phobos and Deimos orbit almost on the equatorial plane of Mars?

As discussed above, the results and analytical arguments presented in this work have strengthened the giant impact origin of Phobos and Deimos. Together with the expected material properties of the building blocks of Phobos and Deimos (such as particle sizes; material provenance: Mars material or impactor material; and thermodynamic properties) that have also been investigated in the framework of the giant impact hypothesis (Ronnet et al. 2016; Hyodo et al. 2017b), our results would finally be tested by a future sample return mission such as JAXA’s Martian Moons eXploration (MMX) mission.

This work was supported by JSPS Grants-in-Aid for JSPS Fellows (JP17J01269) and by the JSPS-MAEDI bilateral joint research project (SAKURA program). P.R. was financially supported by the Belgian PRODEX program managed by the European Space Agency in collaboration with the Belgian Federal Science Policy Office. H.G. was financially supported by JSPS Grants-in-Aid for Scientific Research (JP15K13562, JP17H02990, JP17H06457). S.C. was supported by the UnivEarthS Labex programme at Sorbonne Paris Cite (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02). Numerical computations were partly performed on the S-CAPAD platform, IPGP, France. S.C. acknowledges the support of the CNRS/INFINITI program. Finally, we thank Dr. Aurélien Crida for comments that greatly improved the presentation of the manuscript.

The equilibrium theory considers axisymmetric loads, thus allowing the computation of the polar wander in latitude alone. Although the Borealis basin is slightly elliptical, we will assume an axisymmetric load for direct application of the computation developed in the equilibrium theory. In this theory, planet tends to reorient in response to the modification of the planet inertia tensor induced by mass excess or deficit (Matsuyama et al. 2006). However, the rotational bulge mass excess counteracts this effect, and the efficiency of the TPW depends on the ratio between mass excess/deficit and rotational bulge load that corresponds to the following coefficient:

where a, g, and M, are the equatorial radius (3400 km), the gravity (3.711 m s⁻²), the mass ( kg), and the rotation rate of Mars ( rad s⁻¹), respectively, and is the tidal fluid Love number for a planet without a lithosphere (i.e., 1.1867, see Table 1 in Matsuyama et al. 2006). is the second degree zonal term of the spherical expansion of the surface density of the axisymmetric mass deficit (therefore it has a negative value; see Gold 1955). It is given as

where θ is the co-latitude, is the co-latitudinal extent of the axisymmetric deficit load, and its surface density. If we assume this surface density is constant and that the shape of the cavity can be approximated by a spherical cap, it becomes

where ρ is the volume density of the excavated material (3000 kg m⁻³), and h is the height of the spherical cap, corresponding to the crustal dichotomy thickness (26 km, Neumann et al. 2004).

The α coefficient in Equation (1) is computed as in Equation (26) of Matsuyama et al. (2006):

where k_f^T and k_f^L are the fluid Tidal and load Love numbers for a planet with a lithosphere of a given thickness. The values of the Love numbers are extracted from Figure 3 of Matsuyama & Manga (2010) and given in Table 1 as a function of lithospheric thickness. We considered a thin lithosphere as expected in early Mars history.

Considering the spherical triangle formed by the post-Tharsis TPW geographical pole, the pre-Tharsis TPW geographical pole (or paleopole) and the center of Borealis, basic spherical trigonometric relationship yields

with and are the initial and final co-latitudes of Borealis center, respectively. δ is the Tharsis-TPW displacement and is the longitude shift between Tharsis and Borealis central meridians. Given the Borealis final co-latitude, (latitude of 67°, Andrews-Hanna et al. 2008), the Tharsis TPW displacement δ = 189 (Matsuyama & Manga 2010) and the longitude shift of 50° (see Section 2), it yields .