<figure><img src="image_1.png" /></figure>

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

Ryuki Hyodo¹,² ID, Pascal Rosenblatt³,<fnref n="4" />, Hidenori Genda¹ ID, and Sébastien Charnoz² ID

¹Earth-Life Science Institute/Tokyo Institute of Technology, 2-12-1 Tokyo, Japan; hyodo@elsi.jp

²Institut de Physique du Globe/Université Paris Diderot/CNRS, F-75005 Paris, France

³Royal Observatory of Belgium, B-1180 Brussels, Belgium

Received 2017 June 24; revised 2017 October 31; accepted 2017 November 6; published 2017 December 19

## Abstract

Phobos and Deimos are the two small Martian moons, orbiting almost on the equatorial plane of Mars. Recent works have shown that they can accrete within an impact-generated inner dense and outer light disk, and that the same impact potentially forms the Borealis basin, a large northern hemisphere basin on the current Mars. However, there is no a priori reason for the impact to take place close to the north pole (Borealis present location), nor to generate a debris disk in the equatorial plane of Mars (in which Phobos and Deimos orbit). In this paper, we investigate these remaining issues on the giant impact origin of the Martian moons. First, we show that the mass deficit created by the Borealis impact basin induces a global reorientation of the planet to realign its main moment of inertia with the rotation pole (True Polar Wander). This moves the location of the Borealis basin toward its current location. Next, using analytical arguments, we investigate the detailed dynamical evolution of the eccentric inclined disk from the equatorial plane of Mars that is formed by the Martian-moon-forming impact. We find that, as a result of precession of disk particles due to the Martian dynamical flattening $J_2$ term of its gravity field and particle-particle inelastic collisions, eccentricity and inclination are damped and an inner dense and outer light equatorial circular disk is eventually formed. Our results strengthen the giant impact origin of Phobos and Deimos that can finally be tested by a future sample return mission such as JAXA's Martian Moons eXploration mission.

**Key words:** planets and satellites: formation – planets and satellites: individual (Phobos and Deimos)

## 1. Introduction

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 1°08 and 1°79 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 (~3R<sub>Mars</sub> where R<sub>Mars</sub> 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 ~0.03 M<sub>Mars</sub> (about 1/3 Martian radius) and an impact velocity of ~6 km s<sup>-1</sup> (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

<footnote n="4">⁴ Current Address: ACRI-ST, Nice, France.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Predictions of true polar wander (TPW) due to Borealis impact basin as a function of lithospheric thickness as labeled on the plot. The horizontal dotted line represents the pre-Tharsis TPW (or post-Borealis TPW) co-latitude of the Borealis basin (see the text).</figcaption></figure>

circular disk, forming equatorial orbiting Phobos and Deimos. In Section 4, we summarize our results.

## 2. True Polar Wander

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 $\delta$ of the pole at the surface of the planet is given by the following relationship:

$$ \delta = \frac{1}{2} \tan^{-1} \left[ \frac{Q' \alpha \sin(2\theta_L)}{1 - Q' \alpha \cos(2\theta_L)} \right], \quad (1) $$

where $\theta_L$ is the co-latitude of the center of the axisymmetric excess mass before TPW, $Q'$ is the ratio of the mass excess over the rotational bulge load, and $\alpha$ 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, $Q'$ 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.

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 18°9 (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 18°9 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 ~45° 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.

### 3. Dynamical Evolution of the Disk

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 $a_{eq}$ for every disk particles ($a_{eq} = a(1 - e^2)$) 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

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Distribution of disk particles on semimajor axis to eccentricity plane in the case of the Borealis basin-forming canonical impact (impact angle of 45° with respect to the Martian surface and $N_{SPH} = 3 \times 10^6$ at $T = 5$ hr, where $N_{SPH}$ is the number of SPH particles and $T$ is the time after the start of simulation). Data are taken from Hyodo et al. (2017b). Here, we plot only particles whose pericenter distances are larger than the radius of Mars (disk particles). The solid black lines from top to bottom represent Equation (2) with $R_{Mars} = 3260$ km and $\Delta r = 0, 0.5$ and $1.0R_{imp}$, respectively, taken from the simulation (Hyodo et al. 2017b), where $R_{imp} \sim 1000$ km is the radius of the impactor.</figcaption></figure>

eccentric debris to the circular equatorial disk that can form Phobos and Deimos.

#### 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 $r_{peri} \sim R_{Mars} + \Delta r = a(1 - e)$ where $\Delta r$ 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

$$e = 1 - \frac{R_{Mars} + \Delta r}{a} \qquad (2)$$

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.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Distribution of the longitude of the pericenter in the case of the Borealis basin-forming canonical impact ($N_{SPH} = 3 \times 10^6$ at $T = 5$ hr obtained from Hyodo et al. (2017b), see also their Figure 1).</figcaption></figure>

## 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.

### 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_2$ term), i.e., as their orbits (the argument of pericenter $\omega$ and the longitude of ascending node $\Omega$) start to precess around Mars. Using $J_2 = 1.96 \times 10^{-3}$ 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 $T_{\omega} = 2\pi/\dot{\omega}$ and $T_{\Omega} = 2\pi/\dot{\Omega}$, where $\dot{\omega} = \frac{3n}{(1-e^2)^2} \left(\frac{R_{Mars}}{a}\right)^2 \left(1 - \frac{5}{4}\sin^2(i)\right) J_2$ and $\dot{\Omega} = -\frac{3n \cos(i)}{2(1-e^2)^2} \left(\frac{R_{Mars}}{a}\right)^2 J_2$, respectively (Figure 4). Note that the actual value of $J_2$ 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 ($e \sim 0.5-0.9$ and $a \sim 2-10 R_{Mars}$), 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 $|T_{\omega}/T_{\Omega}| = \left| \frac{\cos(i)}{2(1 - \frac{5}{4} \sin(i)^2)} \right|$ that only depends on $i$, and these timescales are comparable for the nominal case of $i = 45^{\circ}$.

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

$$T_{\text{col}} = \frac{1}{n_{\text{p}} \sigma_{\text{col}} v_{\text{rel}}}, \quad (3)$$

where $n_{\text{p}}$ is the number density of particles, $\sigma_{\text{col}}$ is the collision cross section and $v_{\text{rel}}$ 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 $\sigma_{\text{col}} = 4\pi r_{\text{p}}^2$, where $r_{\text{p}}$ 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

$$n_{\text{p}}(r, \psi) = \int da \int de \int di \frac{N(a, e, i) P(r|a, e) P(\psi|i)}{2\pi r^2}, \quad (4)$$

where $N(a, e, i)$ is the number of particles as a function of $a, e$, and $i$. $P(r|a, e)$ and $P(\psi|i)$ are the probability of finding a particle at radial distance $r$ and an angle $\psi$ from the equatorial plane of the planet and are written as

$$P(r|a, e) = \frac{r}{\pi a \sqrt{a^2 e^2 - (a - r)^2}}, \quad (5)$$

and

$$P(\psi|i) = \frac{|\cos(\psi)|}{\pi \sqrt{\sin(i)^2 - \sin(\psi)^2}}, \quad (6)$$

respectively.

Here, to estimate the collisional timescale, we first integrate $n_{\text{p}}$ over $a$ and $e$ assuming the disk mean inclination of either $i = 10^{\circ}, 45^{\circ}$ and $80^{\circ}$. 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 $r_{\text{p}} \sim 1.5$ m. We assume particle density of $\rho_{\text{p}} = 2500 \text{ kg m}^{-3}$ to calculate the number of particles as $N(a, e, i) = m_{\text{SPH}}/m_{\text{p}}$, where $m_{\text{SPH}}$ is the mass of one SPH particle whose orbital
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Precession timescale $T_{\text{pre}}$ as a function of semimajor axis and eccentricity at different disk inclination from the equatorial plane of Mars (from top to bottom panels, $i = 10^{\circ}$, $45^{\circ}$, and $80^{\circ}$ are shown, respectively). Left panel shows that for argument of pericenter ($T_\omega$) and right panel shows that of longitude of the ascending node ($T_\Omega$). White line shows the analytical $a-e$ distribution of the initial disk particles (Equation (2)). When $T_{\text{pre}} < 1$ years, it is plotted with black, and when $T_{\text{pre}} > 10^9$ years, it is plotted with white.</figcaption></figure>

elements are *a* and *e*, and $m_p = 4\pi\rho_p r_p^3/3$ 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 $r = 1.1r_{\text{per}}$ and $\psi = 0.9i$. Then, assuming the relative velocity of $v_{\text{rel}} \sim v_{\text{peri}} \sin(45^{\circ})$ for orbital elements of $e = 0.8$, we estimate the collision timescale of less than $1 \times 10^{-5}$ year for $i = 10^{\circ}$, $45^{\circ}$, and $80^{\circ}$. 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 $r_{peri}$ but have different orbital elements ($e_{1,2}$ where $e_1 > e_2$) under the phase alignment as seen in our case (see Figures 2 and 3), collision velocity can be estimated as

$$v_{\text{col,alig}} = \sqrt{\frac{GM}{r_{\text{peri}}}} (\sqrt{1+e_1} - \sqrt{1+e_2}) = 3.2 \text{ km s}^{-1} \\ \times \left(\frac{r}{3800 \text{ km}}\right)^{-1/2} \left(\frac{M}{6 \times 10^{23} \text{ kg}}\right)^{1/2} \\ \times (\sqrt{1+e_1} - \sqrt{1+e_2}). \quad (7)$$

Thus, if collision occurs between particles whose eccentricities are $e=0.9$ and $e=0$, respectively, we get $v_{\text{col,alig}} \sim 1.2 \text{ km s}^{-1}$. Then, their initial collision timescale is expected to be their synodic period when there is no dense circular inner rings as

$$T_{\text{syn}} = 2\pi a / \left(\frac{3}{2} \Delta a \Omega\right) \sim 15 \text{ days} \times \left(\frac{M}{6 \times 10^{23} \text{ kg}}\right)^{-1/2} \\ \times \left(\frac{\Delta a}{3400 \text{ km}}\right)^{-1} \left(\frac{a}{3.4 \times 10^4 \text{ km}}\right)^{5/2}, \quad (8)$$

where $\Delta a$ 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 $a = 10R_{\text{Mars}}$). Such high collision velocity may quickly damp the system and be energetic enough to form ~100 µm 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 (~1–100 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_2$ 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^\circ$ and $a=5.5R_{\text{Mars}}$ with the radial width of $\Delta a = 3R_{\text{Mars}}$ where Phobos and Deimos are expected to form with a disk mass of $M_{\text{Phobos}}$, and particle size of 1 m and its density 2500 kg m$^{-3}$) is estimated by using the same argument above (Equation (3)), where $n_p = N/V$, $N = M_{\text{Phobos}}/m_p = 10^{16} \text{ kg}/10^4 \text{ kg} = 1 \times 10^{12}$, and $V = (2\pi a) \Delta a (2a \tan(i)) \sim 4.4 \times 10^{22} \text{ m}^3$. Thus, $n_p \sim 2.2 \times 10^{-11}$, $\sigma_{\text{col}} = \pi (1.0 \text{ m})^2 \sim 3 \text{ m}^2$ and $v_{\text{rel}} \sim v_{\text{Kep}} \sin(i) = \sqrt{GM/a} \sin(i) = 1000 \text{ m s}^{-1}$. Thus, $T_{\text{col}} \sim 0.5$ years. This timescale is actually the maximum timescale of collision because we assume the scale height of $2a \tan(i)$ 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_2$ (see more details Hyodo et al. 2017a). In the case of the canonical Martian-moon-forming disk, the inner disk ($a < 4R_{\text{Mars}}$) is massive (~$10^{20}$ kg) and the outer disk (4–7$R_{\text{Mars}}$) has a mass of Phobos (~$10^{16}$ 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
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Longitude of the nodes of disk particles under the effect of $J_2$ at different epoch obtained from our *N*-body simulations. Initially, particles on circular orbits are distributed between 4–7$R_{Mars}$, where Phobos and Deimos are expected to form (Rosenblatt et al. 2016). The left panel shows the case of no inner massive satellite and the right panel shows the case with a massive inner satellite. A large dot in the right panel shows a massive satellite. Black dots show outer disk particles at $T=0$ year, green dots show those of 250 years, and blue dots show those after 500 years.</figcaption></figure>

*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_2$ 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^{20}$ kg) on circular orbit ($a = 2.5R_{Mars}$) 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_2$ 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 (~$10^{16}$ 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 $4R_{Mars}$ 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, $a = 4–7R_{Mars}$, and $e \sim 0$), 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).

## 4. Summary and Discussion

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 $3 \times 10^{29}$ 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 $\sim 1$ km s$^{-1}$. 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.

<table><caption>**Table 1**
Love Numbers of Mars and Parameter Values for Varying Lithospheric Thickness Values</caption><thead><tr><th>Lithospheric Thickness (km)</th><th>Load Fluid Love Number $k_f^L$</th><th>Tidal Fluid Love Number $k_f^T$</th><th>α Parameter</th></tr></thead><tbody><tr><td>50</td><td>-0.9150</td><td>1.1100</td><td>1.2708</td></tr><tr><td>100</td><td>-0.8529</td><td>1.0370</td><td>1.1686</td></tr><tr><td>200</td><td>-0.6906</td><td>0.8475</td><td>1.0822</td></tr></tbody></table>

## Appendix
The Equilibrium Theory

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 $Q'$ coefficient:

$$Q' = \frac{\frac{4\pi a^3 g}{5M} L'_{20}}{\frac{-1}{3\sqrt{5}} a^2 \Omega^2 k_f^{T*}}, \quad (9)$$

where $a$, $g$, and $M$, are the equatorial radius (3400 km), the gravity ($3.711 \text{ m s}^{-2}$), the mass ($6.4 \times 10^{23}$ kg), and the rotation rate of Mars ($7.08 \times 10^{-5} \text{ rad s}^{-1}$), respectively, and $k_f^{T*}$ is the tidal fluid Love number for a planet without a lithosphere (i.e., 1.1867, see Table 1 in Matsuyama et al. 2006). $L'_{20}$ 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

$$L'_{20} = 2\pi \int_{0}^{\theta_0} \frac{3 \cos^2 \theta - 1}{2} \rho_s(\theta) d\theta, \quad (10)$$

where $\theta$ is the co-latitude, $\theta_0$ is the co-latitudinal extent of the axisymmetric deficit load, and $\rho_s$ 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

$$L'_{20} = -\pi \cos(\theta_0) \sin(\theta_0)^2 \rho h / 2, \quad (11)$$

where $\rho$ is the volume density of the excavated material ($3000 \text{ kg m}^{-3}$), and $h$ is the height of the spherical cap, corresponding to the crustal dichotomy thickness (26 km, Neumann et al. 2004).

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

$$\alpha = \frac{1 + k_f^L}{1 - \frac{k_f^T}{k_f^{T*}}}, \quad (12)$$

where $k_f^L$ and $k_f^T$ 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

$$
\cos(\theta_L) = \cos(\theta_f)\cos(\delta) + \sin(\theta_f)\sin(\delta)\cos(\Delta L), \quad (13)
$$

with $\theta_L$ and $\theta_f$ are the initial and final co-latitudes of Borealis center, respectively. $\delta$ is the Tharsis-TPW displacement and $\Delta L$ is the longitude shift between Tharsis and Borealis central meridians. Given the Borealis final co-latitude, $\theta_f = 23^\circ$ (latitude of 67°, Andrews-Hanna et al. 2008), the Tharsis TPW displacement $\delta = 18^\circ 9$ (Matsuyama & Manga 2010) and the longitude shift of 50° (see Section 2), it yields $\theta_L = 18^\circ$.

ORCID iDs

Ryuki Hyodo https://orcid.org/0000-0003-4590-0988

Hidenori Genda https://orcid.org/0000-0001-6702-0872

Sébastien Charnoz https://orcid.org/0000-0002-7442-491X

References

Andrews-Hanna, J. C., Zuber, M. T., & Banerdt, B. W. 2008, Natur, 453, 1212

Batygin, K. 2012, Natur, 491, 418

Bottke, W. F., & Andrews-Hanna, J. C. 2017, NatGe, 10, 344

Bouley, S., Baratoux, D., Matsuyama, I., et al. 2016, Natur, 531, 344

Burns, J. A. 1978, VA, 22, 193

Canup, R. M. 2004, Icar, 168, 433

Chan, N.-H., Mitrovica, J. X., Daradich, A., et al. 2014, JGRE, 119, 169

Citron, R. I., Genda, H., & Ida, S. 2015, Icar, 252, 334

Craddock, R. A. 1994, LPI, 25, 293

Craddock, R. A. 2011, Icar, 211, 1150

Crida, A., & Charnoz, S. 2012, Sci, 338, 1196

Gold, T. 1955, Natur, 175, 526

Hesselbrock, A., & Minton, D. A. 2017, NatGe, 10, 266

Hyodo, R., & Charnoz, S. 2017, AJ, 154, 34

Hyodo, R., Charnoz, S., Ohtsuki, K., & Genda, H. 2017a, Icar, 282, 195

Hyodo, R., Genda, H., Rosenblatt, P., & Charnoz, S. 2017b, ApJ, 845, 125

Hyodo, R., & Ohtsuki, K. 2014, ApJ, 787, 56

Hyodo, R., & Ohtsuki, K. 2015, NatGe, 8, 686

Hyodo, R., Ohtsuki, K., & Takeda, T. 2015, ApJ, 799, 40

Ida, S., Canup, R. M., & Stewart, G. R. 1997, Natur, 389, 353

Kite, E. S., Matsuyama, I., Manga, M., Perron, J. T., & Mitrovica, J. X. 2009, EPSL, 280, 254

Marinova, M. M., Aharonson, O., & Asphaug, E. 2008, Natur, 453, 1216

Matsuyama, I., & Manga, M. 2010, Geophys. Res., 115, E12020

Matsuyama, I., Mitrovica, J. X., Manga, M., Perron, J. T., & Richards, M. A. 2006, JGR, 111, E02003

Morbidelli, A., Tsiganis, K., Batygin, K., Crida, A., & Gomes, R. 2012, Icar, 219, 737

Murchie, S. L., Britt, D. T., Head, J. W., et al. 1991, JGR, 96, 5925

Neumann, G. A., Zuber, M. T., Wieczorek, M. A., et al. 2004, JGRE, 109, E08002

Ohtsuki, K. 1993, Icar, 106, 228

Ohtsuki, K., & Ida, S. 1998, Icar, 131, 393

Perron, J. T., Mitrovica, J. X., Manga, M., Matsuyama, I., & Richards, M. A. 2007, Natur, 447, 840

Ronnet, T., Vernazza, P., Mousis, O., et al. 2016, ApJ, 828, 109

Rosenblatt, P., Charnoz, S., Dunseath, K. M., et al. 2016, NatGe, 9, 581

Willemann, R. J. 1984, Icar, 60, 701