Two small satellites of Mars, Phobos and Deimos, are thought to have formed from debris ejected by a giant impact onto early Mars (Craddock 2011; Rosenblatt & Charnoz 2012; Citron et al. 2015; Rosenblatt et al. 2016; Hesselbrock & Minton 2017; Hyodo et al. 2017; Canup & Salmon 2018) . The small masses of the satellites, compared to expectations from giant impact simulations, has lead to proposals that the original satellite system was more massive, with most of its mass falling onto Mars over several Gyr following the impact (see Canup & Salmon 2018, for an alternate view). The early massive satellite hypothesis offers at least two competing explanations of Phobos's short (40 Myr; Lainey et al. 2007 ) lifetime against inward tidal evolution. Craddock (2011) and Rosenblatt et al. (2016) suggest that several satellites formed from the impact-generated debris disk and most have since tidally decayed onto Mars, with Phobos and Deimos being the last survivors due to their formation just interior and exterior to the martian synchronous orbit, respectively. Alternatively, Hesselbrock & Minton (2017) propose that Phobos is only the latest product of a repeating ring-satellite cycle at Mars, with each successive inner satellite being less massive than the preceding one. In the ring-satellite cycle model, satellites form from the outer edge of the ring, and then migrate outward through gravitational interaction with the ring. The ring loses mass to the planet at its inner edge, and once the ring is sufficiently depleted the satellite migrates inward due to tides. The inward migration of the satellite 2Ćuk et al.
continues until tidal stresses disrupt it to form a new ring (cf. Black & Mittal 2015; Hesselbrock & Minton 2019) , which viscously spreads, thus re-starting the cycle.
One difficulty for the giant impact scenario is the relatively large distance at which Deimos orbits Mars, beyond the martian synchronous orbit. In the ring-satellite cycle model of Hesselbrock & Minton (2017) , only satellites interior to the synchronous orbit would cycle, and therefore Deimos would have to have formed from the original impactgenerated circumplanetary disk, as in Craddock (2011) . Alternatively, Rosenblatt et al. (2016) have suggested that Deimos (or its precursors) may have been "pushed out" through a resonance with an ancient massive inner satellite. While attractive, this scenario may have difficulty in producing Deimos on its present orbit, as we would expect the orbit of a resonantly-pushed satellite to be quite excited (e.g. Malhotra 1993) , and tidal dissipation within Deimos may be too inefficient to damp its eccentricity over the age of the solar system (Szeto 1983) . While early work on Deimos's cratering record estimated that its surface dates back to the end of the Late Heavy Bombardment (Thomas & Veverka 1980) , that work did not consider planetocentric and, especially, sesquinary impacts. Nayak et al. (2016) show that the sesquinary impacts on Deimos are an important source of cratering, making determination of its surface age very difficult. This would be particularly true if Deimos had a more eccentric or inclined orbit in the past, which would increase relative energy of sesquinary impactors which have similar eccentricities and inclinations as Deimos, but whose orbits have precessed out of alignment with it. Therefore, if the orbit of Deimos was highly eccentric and inclined in the past, it could have suffered runaway sesquinary bombardment, which in turn could have damped Deimos's e and i, potentially reconciling an early period of strong dynamical excitation with its current near-circular/equatorial orbit.
The present-day orbit of Deimos can be used to constrain the martian satellites' dynamic past. The current mean inclination of Deimos relative to the Laplace plane is i f = 1.8 • , while its mean eccentricity is only e = 2.7 × 10 −4 (Jacobson & Lainey 2014) . Yoder (1982) found that the inward orbital migration of Phobos through the satellites' mutual 2:1 mean-motion resonance (MMR) excites the eccentricity of Deimos to e = 0.002, an order of magnitude higher than it is now. Yoder (1982) speculated that higher-than-expected tidal damping of eccentricity by Deimos (possibly due to resonantly excited rotational librations) could explain this discrepancy, but available data on the moments of inertia of Deimos do not support exceptionally large librational response (Thomas 1993) . Recently, Quillen et al. (2020) proposed that Deimos could have experienced enhanced dissipation due to chaotic rotation (c.f. Wisdom 1987) down to very low eccentricities. However, we think that the results of (Quillen et al. 2020) are an artifact of their accelerated simulations, in which the eccentricity damping timescale was artificially made shorter than the spin-down and wobble damping timescales. The problem of Deimos's eccentricity arises for any scenario in which Phobos starts its inward evolution beyond the location of 2:1 resonance at 4.3 R M . In such "ancient" Phobos model, including those of Craddock (2011), Rosenblatt et al. (2016) , and Canup & Salmon (2018) , Phobos should have crossed the 2:1 resonance with Deimos and excited its orbit about 2 Gyr ago (Yoder 1982) . On the other hand, the cyclic martian satellite hypothesis of Hesselbrock & Minton (2017) has Phobos (and previous-generations of inner satellites) migrating outward from the edge of the ring at 3.2 R M , and then reversing migration before ever reaching the 2:1 resonance with Deimos, consistent with the observed low eccentricity of Deimos.
The inclination of Deimos presents a different challenge, as the satellite's almost 2-degree tilt is quite remarkable. High orbital inclinations in otherwise dynamically cold satellite systems are often explained by convergent migration of two satellites into a second-order mean motion resonance (MMR), resulting in a capture into an inclination type resonance (Tittemore & Wisdom 1989; Malhotra & Dermott 1990 ). Deimos is not expected to migrate significantly due to tides (Yoder 1982) , and should have migrated only a very small amount due to radiation force (chiefly the Binary YORP effect;Ćuk & Burns 2005), even if we assume that Deimos's surface shape was unchanged over several Gyr. Therefore, we would need Phobos, or a more massive precursor inner satellite of Mars from an earlier cycle, to migrate outward through a second order resonance with Deimos. The innermost eligible resonance is the 3:1 MMR with Deimos at about 3.3 R M , which could have plausibly been crossed by an inner satellites of Mars (including Phobos) migrating outward due to torques from the martian ring (Hesselbrock & Minton 2017) . In the next Section, we present numerical simulations of this resonance, and show that this process produces an exceptionally good match for the present eccentricity and inclinations of Deimos. In further sections we will briefly present the analysis of the dynamical mechanism involved, and discuss implications for the history of the martian satellite system.
For integrations of the 3:1 MMR between a hypothetical past martian inner satellite and Deimos we use the numerical integrator simpl, which was extensively tested byĆuk et al. (2016) . simpl is a N-body mixed-variable symplectic The inner satellite is assumed to be migrating outward due to ring torques at a net rate that is 16% of its inward tidal migration rate. In the top panel, the gray line shows the semimajor axis at which a satellite would be in exact 1:3 commensurability with the inner satellite. The bottom panel shows the resonant argument 3λ − λI − Ω f − ΩL (see text). The resonance breaks due to libration growth caused by a secondary resonance starting at 0.7 Myr.
integrator that integrates the full equations of motion. We refer the reader toĆuk et al. (2016) for the details of the integrator, and we will here simply specify the parameters used in our simulations. We included the second-order oblateness moment of Mars J 2 (assumed to be same as current), the Sun (assuming a Keplerian orbit for Mars) and two satellites as point-masses. The two satellites in the model include Deimos, with its current mass, and an inner satellite, which we varied between 1 − 200× Phobos mass. We also varied the obliquity of Mars and the eccentricity of the martian heliocentric orbit. We assumed that both satellites had very low eccentricities and inclinations at the time of the 3:1 resonance crossing (see Section 4). We included both the tidal acceleration (which makes the inner satellite migrate inward), and a larger artificial tangential acceleration on the inner satellite that approximates ring torques, resulting in net outward migration. We otherwise ignored the ring's gravity, as a 100 M P ring would increase the martian obliquity moment J 2 by only about 0.1%. We used a timestep of 4 × 10 −5 yr (0.35 h) in all simulations shown here. We find that capture into an inclination-type subresonance of 3:1 MMR is a common outcome for a range of masses, whenever the net orbital migration rate of the inner satellite is slow enough. The resonance capture is not indefinite, but the resonance breaks at Deimos inclinations between zero and three degrees, with Deimos's final inclination varying both stochastically and systematically. Figure 1 shows the evolution of the satellites' semimajor axes and inclinations during the resonance capture for an inner satellite 18 times more massive than Phobos. The inclination of the inner satellite is not affected by this resonance, therefore implying that we are looking at i 2 type subresonance, in which 4Ćuk et al. only the inclination of Deimos is affected (Murray & Dermott 1999, Section 8.8 .2). However, we find that the relevant resonant argument is 3λ − λ I − Ω f − Ω L (Fig. 1, bottom panel) , where λ and Ω f are the mean longitude and longitude of the node of Deimos, and λ I is the longitude of the inner satellite. The resonant argument also depends on the orientation of Deimos's Laplace plane Ω L (measured with respect to the martian equator), which is unusual and will be discussed in the next Section. The dependence of the final inclination of Deimos after the resonance passage on the mass of the inner satellite is plotted in Fig. 2 . Each point plots the end-state of a numerical simulation. For many inner satellite mass values, we ran sets of simulations that differed only in the separation of the system from the resonance, enabling us to survey the inherent stochasticity of the process. Red diamonds plot simulations that used a martian obliquity of ǫ = 25 • , comparable with the present state, while black circles plot simulations with martian obliquity of ǫ = 40 • , closer to the long-term average (Laskar et al. 2004) . While those simulations used the martian heliocentric obliquity of e M = 0.09, blue triangles plot simulations with e M = 0.04 (with ǫ = 40 • ); we find no systematic dependence of the resonant outcomes on martian eccentricity. Plotted simulations generally used net outward migration rates that were about 16% of the (inward) tidal evolution rate. Fig. 2 shows that there is an increase of final inclination of Deimos with decreasing mass of the inner satellite in the 3-200 Phobos mass (M P ) range. Less massive inner satellites do not capture Deimos into resonance at this relative migration rate, and Deimos is left with a very low final inclination.
The simulations plotted in Fig. 2 visibly cluster along discrete curves, and these curves are offset for simulation sets that used different obliquities of Mars. We find that this pattern is caused by secondary resonances, which increase the libration amplitude of the i f i L resonance, eventually leading to breaking of the resonance lock. A secondary resonance is visible in the bottom panel of Fig. 1 , starting at 0.7 Myr into the simulation. These secondary resonances are very similar to those found by Malhotra & Dermott (1990) in the context of the Miranda-Umbriel 3:1 inclination-type MMR, and their origin will be discussed in the next Section.
The planes of all satellite orbits precess around their respective Laplace planes, with the orientation of the Laplace plane being determined by the sum of perturbations on the orbit. Close-in satellites like Phobos are chiefly perturbed by the planet's oblateness, and their Laplace plane is very close to the planet's equator. Distant satellites like the Moon are overwhelmingly perturbed by the Sun, so their Laplace plane is almost identical to the plane of the planet's heliocentric orbit. Deimos is intermediate between these two regimes, as its Laplace plane is not exactly in the plane of the martian equator, but is inclined by about a degree toward Mars's heliocentric orbit. Starting with the gravitational potential derived in Tamayo et al. (2013) (see also Tremaine et al. 2009 ), we can derive Deimos's Laplace plane tilt by requiring the terms containing cos Ω in the disturbing function to be zero (this is the requirement for uniform nodal precession):
where M Sun is the solar mass, and M , a M , R and J 2 are respectively the mass, semimajor axis, radius, and oblateness moment of Mars, and a is Deimos's semimajor axis. We also assumed circular orbits and i L << ǫ. Eq. 1 shows that the maximum i L happens at ǫ = 45 • , when i L ≃ 1.15 • . The inclination-type 3:1 MMR we are interested in has (apart from terms with no angular dependence) a resonant term in the disturbing function of the type
where f 57 is a function of masses and semimajor axes (Murray & Dermott 1999 , Table B .8). These terms can be derived using the geometry of the interaction between the inner satellite and Deimos, and this derivation assumes that the inclinations and nodes of both orbits are measured relative to the same reference plane. This last assumption is incorrect, as the Laplace plane of Deimos is tilted by a non-trivial amount. In order to obtain an accurate description of mutual positions and gravitational interactions, we will take the martian equator as the reference plane, and recognize that sin i sin Ω = sin i f sin Ω f + sin i L sin Ω L and sin i cos Ω = sin i f cos Ω f + sin i L cos Ω L , where the subscript f now refers to the "free" inclination and node of Deimos (Murray & Dermott 1999, Section 7.4 ). We can then transform the resonant argument in Eq. 2 so that we separate factors of sin i cos Ω and sin i sin Ω, express them in terms of free and forced components, and then reconstitute the new resonant terms. It can be shown that the original term produces three new resonant arguments, and the one we are interested in is
Therefore, due to the unusual Laplace plane of Deimos, we have a resonance that is derived from a pure inclination (i 2 -type) subresonance of the 3:1 MMR, but in many ways behaves as a mixed ii ′ -type resonance. Physically this makes sense, as the inner satellite is a perturber on an inclined orbit, despite this tilt being due to orientation of Deimos's Laplace plane, and not the inner satellite's own inclination. The above description of the resonance enables us to understand the secondary resonances that lead to Deimos exiting the resonance, as seen in Fig. 1 . The closest neighbor of our i f i L resonance is the pure free-inclination subresonance i 2 f with the resonant argument 3λ − λ I − 2Ω f , which is circulating with Deimos's orbital precession period of about 54 yr while the resonant argument of i f i L resonance is librating. The fact that our simulations do not include Mars's axial and orbital precession is not important, as in realityΩ L /Ω f < 10 −3 . In our simulations we find that the system can be captured into a secondary resonance when the libration frequency in the i f i L resonance is a simple multiple of the circulation frequency of the i 2 f argument. The libration period becomes shorter with increasing i f , and longer for larger amplitudes of libration of the resonant argument. Therefore, in order to maintain the secondary resonance as the free inclination grows, the amplitude of libration of the i f i L resonance has to grow over time. This widening of the libration amplitude is visible in Fig. 1 starting at 0.7 Myr, and eventually leads to breaking of the resonant lock.
In Fig. 1 , the libration frequency at which the secondary resonance happens has a period of about 108 years, making this a 2:1 secondary resonance. In Fig. 2 , we can see the signature of the 2:1 secondary resonance as the topmost lines our simulations cluster around. As the plot suggests, the same secondary resonance happens at higher free inclinations for smaller masses, as librations in i f i L resonance slow down with decreasing mass. The 2:1 secondary resonance also tends to happen at higher free inclination for smaller obliquities, as the line corresponding to the 2:1 secondary resonance for red diamonds (integrated at ǫ = 25 • ) is higher than that for black circles (ǫ = 40 • ). The line below that of the 2:1 secondary resonance corresponds to the 3:1 secondary resonance (i.e. libration periods of about 162 yr), and a few points hint at higher-order secondary resonances. The distribution of probabilities for capture into each of these resonances vary with the inner satellite mass, and while for some inner satellite masses (like 18M P used in Fig. 1) we get one secondary resonance to be dominant, for the inner satellite masses there is a diversity of outcomes. Net Migration Rate (Tidal=-1) Figure 3 . Effects of the inner satellite's net migration rate (normalized to tidal) on the final free inclination of Deimos. One set of simulations uses inner satellite mass of MI = 18MP and the martian obliquity of ǫ = 40 • (black circles), while the other set has MI = 25MP and ǫ = 25 • . Tidal migration rate has been calculated using tidal quality factor Q = 80 and tidal Love number k2 = 0.14 for Mars. High final inclinations of Deimos indicate resonance capture, while low values reflect resonance crossing without capture. While the transition between resonance "jump" and capture is dynamically significant, the variation of final inclinations on the left-hand side may be stochastic and reflect probabilistic capture into different secondary resonances, rather than sensitive dependence on the relative migration rate.
While we have demonstrated in Figs. 1 and 2 that capture into an inclination-type 3:1 MMR between an inner satellite and Deimos is plausible, resonance capture is known to sensitively depend on the rate of orbital evolution. To examine this dependence, we ran two sets of simulations with varying net migration rates, one featuring a M I = 18M P inner satellite and martian obliquity ǫ = 40 • (same as in Fig. 1) , and the other with a M I = 25M P inner satellite and martian obliquity of ǫ = 25 • . The resulting free inclinations of Deimos generated in these simulations are plotted in Fig. 3 . The maximum net migration rate at which resonance capture still occurs is at about 30-35% of the tidal migration rate (itself proportional to the inner satellite's mass). The relative critical migration rate is larger for the more massive inner satellite, implying a dependence on the inner satellite mass (the difference in martian obliquity in Fig. 3 is expected to have the opposite effect, favoring the M I = 18M P , ǫ = 40 • simulations). One direct implication is that the boundary between capture and non-capture into resonance seen on the left-hand side of Fig. 2 is not an intrinsic property of the inner satellite mass, but depends on the migration rate. We ran an additional simulation and confirmed that an inner satellite with M I = 3M P can capture Deimos into the i f i L resonance, if the net migration is about 8% of the tidal evolution rate. Therefore, resonance capture is still possible at very low masses of the inner satellite, but increasingly closely matched ring-torques and tidal acceleration are necessary in order for it to happen.
More broadly, what is the likelihood of this scenario in which the inner satellite migrates outward through the 3:1 resonance with Deimos at 3.3R M at a fraction of the tidal evolution rate? We would expect slow evolution to happen just before the ring torques become too weak to overcome the inward tidal acceleration, and the orbital migration reverses. Hesselbrock & Minton (2017) envision inner satellite formation from the outer edge of the ring at the Fluid Roche limit, which is ∼ 3.1 − 3.4R M for material with mass density in the range of that of Phobos and Deimos (∼ 1.5 − 1.9 g cm −3 ). Interestingly, Hu et al. (2020) recently suggested that the overall shape of Phobos is the best match to tidal forces at ∼ 3.3R M , rather than its present distance, suggesting formation or a prolonged period of residence at that distance, consistent with this model. Our model would favor a formation of the inner satellite at < 3.3R M , followed by some amount of outward orbital evolution from ring torques through the 3:1 MMR. While the exact numbers are model-dependent, Hesselbrock & Minton (2017) do expect the last three generations of the inner satellites (including Phobos) to only migrate out to about that distance. As the ring would be dissipating, the ring Fig. 1 ) of a 18MP inner satellite crossing the 3:1 MMR with Deimos at 16% of its tidal migration rate (martian obliquity of 40 • ). The left-hand side shows the crossing of the e 2 sub-resonance of the 3:1 MMR that closely follows the capture into the i f iL subresonance shown in Fig. 1 . Subsequently, we reduced the outward torque and the inner satellite migrated back through the 3:1 MMR with Deimos, with the second crossing shown on the right-hand side. The modern-day eccentricity of Deimos is e = 2.7 × 10 −4 (Jacobson & Lainey 2014) . torques would become weaker over time, and a very slow outward migration through the 3:1 MMR with Deimos at 3.3R M is certainly plausible in the context of the model of ring-satellite cycling. We also note that tidal dissipation within Mars may have been weaker in the distant past Samuel et al. (2019) , relaxing some of the constraints on the near-canceling of ring and tidal torques.
Apart from the inclination, the 3:1 MMR between the inner satellite and Deimos also affects Deimos's eccentricity. Since there is no equivalent of Laplace plane tilt for eccentricity, the strongest eccentricity subresonance of 3:1 MMR is the e 2 (Murray & Dermott 1999, Section 8.8.2). We do not see capture in that resonance in any of our simulations, but instead Deimos experiences a kick to its eccentricity. Figure 4 shows the eccentricity of Deimos at the moment of the 3:1 MMR crossing in the same simulation as shown in Fig. 1 , with the additional integration of the reverse passage of the inner satellite through the same resonance during its later inward migration. In general, the addition of eccentricity jumps generated in the two passages is random, so sometimes they add up and sometimes they cancel out. We conclude that the excitation of Deimos's eccentricity by resonance passages by a "Cycle 3" satellite (in the terminology of Hesselbrock & Minton 2017) is comparable and consistent with Deimos's current e = 2.4 × 10 −4 . Fig. 2 suggests that the best match to Deimos's inclination is obtained by a slow outward migration of a satellite with a mass of about M I = 20M P through the 3:1 MMR with Deimos, assuming that both satellites have had circular and planar orbits at that time. How does this situation fit into the broader understanding of the history of the martian system? This mass of the inner satellite is close to that predicted for the "Cycle 3" satellite by Hesselbrock & Minton (2017) (numbering goes backward in time; Phobos itself is Cycle 1). This is also expected to be the first generation of past martian inner satellites to contain only one satellite, which is consistent with undisturbed evolution through this relatively weak resonance with Deimos. As Deimos would have its current inclination by the time the Cycle 2 satellite and Phobos encounter the 3:1 resonance with Deimos, resonance capture is excluded, and only small kicks to e and i or Deimos are expected.
As the inner satellite would have recently formed from the ring, we would expect its orbit to be circular and planar. However, Deimos may have experienced resonances with past generations of inner satellites, making its low inclination at the time problematic (note that this is a problem for the eccentricity of Deimos in any scenario involving migrating massive inner satellites). While more modeling work is needed to understand this, we speculate that multiple inner satellites of Cycle 4 (as per Hesselbrock & Minton 2017) excited the eccentricity and/or inclination of Deimos enough that outside bombardment (heliocentric or planetocentric) triggered a runaway cascade of sesquinary impacts which damped the eccentricity and inclination of Deimos, after which re-accretion was possible Nayak et al. (2016) . According to Hesselbrock & Minton (2017) , re-accretion during the Cycle 4 would make Deimos about 3.5-4 Gyr old, which would be broadly consistent with much of its surface craters being heliocentric (Thomas & Veverka 1980) , although we would also expect a significant sesquinary component (Nayak et al. 2016) .
To conclude, we find that the present orbital inclination and eccentricity of Deimos strongly suggest a past 3:1 mean-motion resonance with an outward migrating past inner satellite of Mars about 20 times more massive than Phobos. Since the 3:1 MMR with Deimos is well within the synchronous orbit, this outward migration requires dynamical interaction with a substantial ring of Mars at the time (3-4 Gyr ago). Our results strongly support the cyclic martian satellite hypothesis of Hesselbrock & Minton (2017) . Broader implications include a much larger mass of circumplanetary material around early Mars than is currently in Phobos and Deimos, and a relatively young age for Phobos, requiring most of the craters on its surface to be made by planetocentric impactors.