Co-accretion + Giant Impact Origin of the Uranus System: Post-impact Evolution

We assume that the impulse to the planet from the giant impact destabilizes a prior satellite system, leading to mutually disruptive collisions that produce an outer debris disk with mass ∼10⁻⁴ M_U. To estimate the mass of the inner c-disk needed to cause rapid nodal regression out to distances consistent with low-inclination Oberon, Morbidelli et al. (2012) used a Laplace–Lagrange ring code and an N-body code, mimicking the secular effect of an inner c-disk by treating it as a moon of mass M_c orbiting at 3R_U (hereafter we will refer to this body as the c-moon). We adopt a similar approach, but model the outer debris disk with an N-body simulation (Duncan et al. 1998) that includes inelastic collisions (adopting normal and tangential coefficients of restitution) and mergers when the rebound velocity is below a mutual escape velocity (as in Salmon & Canup 2012). The goal is to assess when nodal randomization driven by the inner c-disk overcomes the tendency for gravitational interactions among fragments in the outer disk to maintain an inclined disk.

Table 1 lists the initial parameters for our simulations. All begin with an outer disk of 5000 equal-mass particles arranged in a radially flat surface density profile, with a total outer disk mass of 10⁻⁴ M_U and with semimajor axes extending from 4–40R_U. The initial N-body particles are large, >50 km in radius, which sets the granularity of the treatment of self-gravity in our simulations. The actual initial size of outer debris fragments is uncertain. Gravitational interactions are treated with an N² algorithm with mutual interactions included for bodies closer than about 6 mutual Hill radii (Duncan et al. 1998). All particles initially have the same eccentricity (0.1), longitude of ascending node (either 0°, 10°, or 20°), and orbital inclinations (either 30°, 45° or 60°). The different initial inclinations correspond to different assumptions for the angle between Uranus' pre-impact and post-impact spin axes, δ, which itself is a function of the assumed value for θ₀, Uranus' pre-impact obliquity, and the azimuthal positions of the planet's spin axes before and after the impact. The angle δ must be less than 90° so that once the disk realigns with the planet's new equatorial plane, disk material will orbit in the same direction as the planet's spin as needed to yield the current satellite system. For a randomly oriented giant impact, the impact orientation needed to satisfy the δ < 90° condition becomes more probable as θ₀ is increased, with 30%–60% of orientations yielding δ < 90° for 10° ≤ θ₀ ≤ 70° (Morbidelli et al. 2012). We identify a second, more stringent requirement that δ must be less than ∼60° to account for Oberon's orbital distance if one assumes the prior satellite system orbited within 10² R_U (see Section 2.3).

Table 1.  N-body Simulations of Nodal Randomization in the Outer Debris Disk

Run	Mc‐moon	Ωd	id	$\sigma \left({{\rm{\Omega }}}_{d,\gt 30{R}_{{\rm{U}}}}\right)$	$\sigma \left({{\rm{\Omega }}}_{d,\gt 30{R}_{{\rm{U}}}}\right)$	〈ed 〉	〈id 〉
				at 100 yr	at 1000 yr
	$\left({M}_{{\rm{U}}}\right)$	(deg)	(deg)	(deg)	(deg)		(deg)
1	0.003	0	30	1.6	4.7	0.088	23.8
2	0.003	10	30	1.6	7.0	0.091	23.7
3	0.003	20	30	1.7	5.8	0.086	23.7
4	0.003	0	45	1.9	2.4	0.106	33.5
5	0.003	10	45	4.8	2.5	0.098	32.5
6	0.003	20	45	7.2	5.3	0.108	32.7
7	0.003	0	60	6.1	5.1	0.119	46.8
8	0.003	10	60	2.5	1.7	0.142	46.2
9	0.003	20	60	1.2	10.2	0.141	46.7
10	0.010	0	30	9.2	46.5	0.098	20.3
11	0.010	10	30	8.7	48.3	0.101	20.6
12	0.010	20	30	8.8	48.1	0.098	20.9
13	0.010	0	45	8.3	50.3	0.126	34.0
14	0.010	10	45	8.5	49.1	0.119	33.4
15	0.010	20	45	6.2	46.1	0.129	33.8
16	0.010	0	60	5.2	35.9	0.147	48.6
17	0.010	10	60	5.6	43.0	0.149	48.9
18	0.010	20	60	3.2	43.9	0.141	49.8
19	0.030	0	30	42.4	52.1	0.167	20.1
20	0.030	10	30	44.5	51.2	0.175	19.7
21	0.030	20	30	46.6	50.8	0.172	20.2
22	0.030	0	45	37.7	50.1	0.198	33.2
23	0.030	10	45	42.1	51.4	0.196	33.4
24	0.030	20	45	42.4	51.4	0.202	33.6
25	0.030	0	60	27.9	53.4	0.209	56.0
26	0.030	10	60	29.6	52.0	0.235	56.4
27	0.030	20	60	26.6	52.7	0.225	55.3

Note. Initial parameters and results from our numerical simulations of nodal randomization in the outer debris disk. M_c‐moon is the mass of the inner c-moon used to mimic the c-disk. We set its initial longitude of ascending node to 0 in all runs. Ω_d and i_d are the initial longitude of ascending node and inclination of disk particles. $\sigma \left({{\rm{\Omega }}}_{d,\gt 30{R}_{{\rm{U}}}}\right)$ is the standard deviation of the longitude of ascending node of disk particles with semimajor axis beyond 30R_U. We use this quantity to measure how much nodal randomization has occurred in the distant regions of the disk; randomization out to ≥30R_U is needed to successfully explain Oberon for i_d ≥ 30°. 〈e_d 〉 and 〈i_d 〉 are the mean eccentricity and inclination of disk particles at t = 1000 yr. All particles have an initial mass of 2 × 10⁻⁸ M_U.

Download table as:  ASCIITypeset image

The large c-moon used to mimic the effect of the c-disk in our simulations has a mass of 0.003, 0.01, or 0.03M_U, and is placed initially at 3R_U on a circular, non-inclined orbit. We set the initial longitude of ascending node of the c-moon to 0. We use normal and tangential coefficients of restitution of _n = 0.01 and _t = 1, similar to other studies of satellite accretion (Ida et al. 1997; Salmon & Canup 2012; Canup & Salmon 2018). We ignore precessional forcing by the Sun. For the current large Uranian moons, solar effects are minimal and the Laplacian plane is coincident with the planet's equatorial plane out to Oberon's distance (e.g., Dobrovolskis 1991). Because the simulations here also include the effect of a massive inner c-disk, the distance at which solar forcing becomes important will be larger still, and so we neglect it here.

Uranus' physical radius at the time of a late giant impact would have been somewhat larger than its current radius (with plausibly R_p ∼ 1.2–1.5R_U; e.g., Bodenheimer & Pollack 1986; Fortney et al. 2007), and its early rotation rate slower by conservation of spin angular momentum. A slower rotation rate would decrease the planet's J₂, since J₂ ∝ ω² (e.g., Bertotti & Farinella 1990), where ω is the spin frequency of the planet. The spin angular momentum of a planet of mass M_U and radius R_P is ${L}_{\mathrm{spin}}={{kM}}_{{\rm{U}}}{R}_{P}^{2}\omega$, where k is the moment of inertia constant. Conservation of spin angular momentum gives ${R}_{{\rm{U}}}^{2}{\omega }_{{\rm{U}}}={R}_{P}^{2}\omega$, where ω_U ∼ 10⁻⁴ s⁻¹ is the current spin of Uranus and we assume for simplicity an early Uranus with a moment of inertia constant comparable to that of the current planet. For an early Uranus with a radius of R_P = 1.3R_U, this yields a post-impact ω ∼ 6.2 × 10⁻⁵ s⁻¹. Accordingly, we set J₂ ∼ 1.3 × 10⁻³ in our simulations, about a factor of 2.6 smaller than the current J₂ of Uranus.

Figure 2 shows the node, eccentricity, and inclination of disk particles at 0, 100, and 1000 yr, for cases with a smaller c-moon (M_c = 0.003M_U; Run 1) and a larger c-moon (M_c = 0.03M_U; Run 26). In the first, the nodes of disk particles efficiently randomize out to distances of about 20R_U. Beyond that distance, the particles retain a common node (i.e., they are precessing coherently), indicating a warped outer disk structure that would yield an inclined outer satellite. In the case with a larger c-moon (Figure 2, right), the nodes of the particles are efficiently randomized across the entire disk after 1000 yr.

Zoom In Zoom Out Reset image size

The planet's J₂ and the c-moon cause outer particle nodes to regress with a regression rate that decreases rapidly with orbital distance while producing no (direct) secular change in particle semimajor axes (a). However, as the initially tilted ring begins to regress, inner regions regress more rapidly, and this differential nodal regression produces a temporary mass distribution that is akin to a leading spiral wave pattern. Self-gravity across this structure produces a negative torque on the outer disk regions in our simulations that affects debris orbital elements (including a). The effect is short-lived, existing only prior to substantial dispersion of the nodes.

Columns 5–8 in Table 1 show several properties of the resulting outer disk particles. To quantify the efficiency of nodal randomization, we measure the standard deviation of the longitude of ascending node of disk particles located beyond 30R_U after 100 and 1000 yr of evolution (this quantity is initially 0 as all disk particles have the same node). We find, as expected, larger randomization for runs with a larger c-moon. For M_c = 0.003M_U, nodes in the 30–40R_U region are minimally dispersed by <10° and retain coherency, while for M_c ≥ 0.01M_U, nodal dispersion reaches about 50° across this region within 10³ yr. Overall, our results agree with those of Morbidelli et al. (2012) on the c-disk mass required to realign the outer regions of the debris disk to Uranus' post-impact equatorial plane.

Consider the Uranus system just after a giant impact by an approximate Earth mass projectile has produced the planet's 98° obliquity and a massive c-disk. The needed c-disk would be compact, with a radius of a few R_U, and it would likely be composed predominantly of water and be (approximately) aligned with the planet's post-impact equatorial plane (Slattery et al. 1992; Reinhardt et al. 2020; Rufu & Canup 2020). The energy of the impact will have heated the planet's surface to temperatures ∼10⁴ K, and vaporized the c-disk. Outer debris produced by disruptive collisions among the prior satellites is expected to contain roughly half ice and half rock, consistent with expected compositions of moons accreted during Uranus' late gas co-accretion (e.g., Canup & Ward 2006). Ice may sublimate due to Uranus' luminosity, even at large distances; e.g., a rotating particle at distance r, with Bond Albedo A_b ∼ 0.1, and emissivity _ir ∼ 1, will be heated to a temperature ${T}_{\mathrm{par}}\sim {[{({R}_{{\rm{U}}}/r)}^{2}(1-{A}_{b})/(4{\epsilon }_{\mathrm{ir}})]}^{1/4}{T}_{{\rm{U}}}$, which for an effective temperature for Uranus of T_U > 2000 K implies T_par > 200 K for (r/R_U) < 50. Resulting water vapor thermal velocities would be less than the local escape velocity. Thus we expect the outer disk may initially contain water vapor and (primarily) rocky debris.

With time Uranus cools, losing the heat delivered by the giant impact. To (crudely) estimate its cooling timescale, consider a case in which the impact energy, $\sim \tfrac{1}{2}{M}_{\oplus }{\left(20\,\mathrm{km}\,{{\rm{s}}}^{-1}\right)}^{2}\sim 1.2\times {10}^{40}$ erg, is deposited in an outer layer of the planet that is heated by ΔT = 10⁴ K, which implies an outer layer mass ∼0.15M_U. For comparison, an SPH simulation of the impact of a 1 Earth mass object into a Uranus-like planet with an impact velocity of 1.1 times the mutual escape velocity and a 45° impact angle heats the outer ∼10% of the planet's mass by ≥10⁴ K (R. M. Canup 2017, personal communication; see also Rufu & Canup 2020). For a well-mixed layer, the time for Uranus to cool to temperature T_U is

$$\begin{eqnarray}&&t\sim \displaystyle \frac{0.15{M}_{{\rm{U}}}C}{4\pi {R}_{{\rm{U}}}^{2}{\sigma }_{\mathrm{SB}}{T}_{{\rm{U}}}^{3}},\end{eqnarray} \tag{ 1 }$$

where C ∼ 10⁸ erg K⁻¹ g⁻¹ is specific heat and σ_SB is the Stefan–Boltzmann constant. As Uranus cools, the ice condensation distance at which T ∼ 200 K moves inward to smaller orbital radii. An opaque vapor disk passively heated by the planet has a temperature $T\sim 0.3{\left(3{R}_{{\rm{U}}}/r\right)}^{3/4}{T}_{{\rm{U}}}$ (e.g., Ruden & Pollack 1991). Combining this with T_U(t) from Equation (1) provides a simple estimate of the time to ice condensation (i.e., T ≤ 200 K) as a function of r. Beyond 10R_U, ice may condense after 10²–10³ yr, which is less than or comparable to the local satellite accretion timescale (see below). However, for r ≤ 5R_U Uranus' luminosity would maintain a water-dominated c-disk as a vapor for 10⁴–10⁵ yr. Deeper energy deposition in the planet and/or less efficient mixing would yield slower cooling than these estimates. However in general, one expects satellite accretion in the outer debris disk would occur before the inner c-disk cools and begins to condense. Accordingly, we first model accretion in the outer debris disk (Section 2.3), and then separately consider the later viscous evolution of the inner c-disk as it starts to condense, spread, and spawn moonlets (Section 3).

We simulate the reaccumulation of the outer debris disk to identify conditions needed to yield a system broadly similar to today's Ariel, Umbriel, Titania, and Oberon in terms of satellite number, masses, and orbital distribution. We perform a suite of N-body accretion simulations, with initial disk parameters indicated in Table 2. The combined mass of Ariel, Umbriel, Titania, and Oberon is M_AUTO = 1.044 × 10⁻⁴ M_U; we use a somewhat larger initial disk mass (1.15 × 10⁻⁴ M_U) to allow for some material loss. We consider a vapor-free disk that is roughly half rock and half ice (with a particle bulk density of 1.5 g cm⁻³), which is plausible if Uranus cools efficiently, and neglect the gravitational potential of the c-disk, since we assume that precession forced by the planet's J₂ and the c-disk has already randomized the outer debris disk nodes.

Table 2.  N-body Simulations of Outer Debris Disk Accretion

Run	N	Rout	q	〈ed 〉	〈id 〉	N	Nsats	Msats	amin	amax	〈esats〉	〈isats〉
		$\left({R}_{{\rm{U}}}\right)$			(deg)			$\left({M}_{\mathrm{AUTO}}\right)$	(RU)	$\left({a}_{\mathrm{Oberon}}\right)$		(deg)
1	5000	23	0	0.1 ± 0.03	45 ± 3	0.01	4	0.94	4.52	0.41	0.018 ± 0.015	1.7 ± 0.35
2	5000	23	1	0.1 ± 0.03	45 ± 3	0.01	5	0.93	2.97	0.41	0.012 ± 0.006	1.4 ± 0.09
3	5000	23	0	0.2 ± 0.03	60 ± 3	0.01	2	0.91	3.45	0.21	0.025 ± 0.004	2.0 ± 0.22
4	5000	23	1	0.2 ± 0.03	60 ± 3	0.01	4	0.82	2.54	0.20	0.003 ± 0.003	2.0 ± 0.13
5	5000	23	0	0.1 ± 0.03	30 ± 3	0.01	4	0.93	6.70	0.65	0.016 ± 0.011	0.6 ± 0.41
6	5000	23	1	0.1 ± 0.03	30 ± 3	0.01	6	1.10	3.82	0.71	0.042 ± 0.020	1.7 ± 0.86
7	5000	23	0	0.1 ± 0.03	45 ± 3	0.1	3	0.91	5.59	0.43	0.006 ± 0.003	1.1 ± 0.27
8	5000	23	1	0.1 ± 0.03	45 ± 3	0.1	5	0.93	3.59	0.39	0.015 ± 0.006	1.3 ± 0.42
9	5000	23	0	0.2 ± 0.03	60 ± 3	0.1	4	1.03	1.92	0.20	0.028 ± 0.008	1.8 ± 0.52
10	5000	23	1	0.2 ± 0.03	60 ± 3	0.1	4	0.86	2.41	0.23	0.007 ± 0.003	2.0 ± 0.04
11	5000	23	0	0.1 ± 0.03	30 ± 3	0.1	4	0.96	7.19	0.68	0.007 ± 0.004	0.8 ± 0.21
12	5000	23	1	0.1 ± 0.03	30 ± 3	0.1	4	1.03	4.06	0.63	0.027 ± 0.012	0.7 ± 0.21
13	1000	23	0	0.1 ± 0.03	10 ± 3	0.1	4	0.99	8.89	0.94	0.021 ± 0.018	0.3 ± 0.12
14	1000	23	1	0.1 ± 0.03	10 ± 3	0.1	4	0.91	8.66	0.97	0.030 ± 0.019	0.9 ± 0.47
15	1500	30	0	0.1 ± 0.03	45 ± 3	0.1	4	0.96	6.67	0.65	0.018 ± 0.006	2.3 ± 0.25
16	1500	30	1	0.1 ± 0.03	45 ± 3	0.1	6	1.04	3.47	0.59	0.014 ± 0.002	2.2 ± 0.48
17	1500	30	0	0.2 ± 0.03	60 ± 3	0.1	5	1.01	2.98	0.31	0.007 ± 0.002	4.2 ± 0.28
18	1500	30	1	0.2 ± 0.03	60 ± 3	0.1	4	0.73	3.52	0.28	0.007 ± 0.006	4.1 ± 0.24
19	1500	30	0	0.1 ± 0.03	30 ± 3	0.1	4	0.88	11.75	1.07	0.014 ± 0.008	1.4 ± 0.59
20	1500	30	1	0.1 ± 0.03	30 ± 3	0.1	4	1.00	5.62	0.79	0.050 ± 0.019	1.5 ± 0.51
21	1500	30	0	0.1 ± 0.03	10 ± 3	0.1	5	1.08	7.52	1.35	0.069 ± 0.029	1.8 ± 0.86
22	1500	30	1	0.1 ± 0.03	10 ± 3	0.1	4	0.87	10.28	1.13	0.030 ± 0.008	0.7 ± 0.51

Note. Results from N-body simulations of satellite accretion in the outer debris disk. All disks have an initial disk mass M_d = 1.15 × 10⁻⁴ M_U, and contain N particles distributed with a surface density profile σ ∝ a^−q . R_out is the semimajor axis of the outermost particles, and 〈e_d 〉 and 〈i_d 〉 are the mean eccentricities and inclination of the disk particles at the beginning of the simulation. _N is the normal coefficient of restitution used when treating collisions between particles. N_sats is the number of large satellites at t = 3 × 10⁴ yr, where we define large as having a mass greater than 10⁻⁵ M_U. M_sats is the total mass of large satellites in units of the combined mass of Ariel, Umbriel, Titania, and Oberon M_AUTO. a_min and a_max are the semimajor axes of the innermost and outermost final large satellites in units of R_U and the current semimajor axis of Oberon a_Oberon, respectively. 〈e_sats〉 and 〈i_sats〉 are the mean eccentricities and inclination of the final large satellites.

Download table as:  ASCIITypeset image

Disk particles are assigned average starting eccentricities of 0.1 or 0.2, and initial inclinations of 30°, 45°, or 60° with random longitudes of ascending node, so that the disk is a thick torus centered on the planet's equatorial plane. We consider two surface density profiles (σ(r) ∝ r^−q with q = 0, 1), and two values for the normal coefficient of restitution (_N = 0.01, 0.1). We continue the simulations for 3 × 10⁴ yr.

Figure 3 (left) shows the evolution of the system from Run 14. Satellites grow as particles collide and merge, and this process is most rapid in the inner region of the disk because collision rates depend on orbital frequency, Ω. Consider an outer debris disk with mass 10⁻⁴ M_U distributed uniformly out to ∼ 25R_U, with surface density σ ∼ 700 g cm⁻². As the disk collisionally evolves, the balance between gravitational stirring and collisional damping will yield an equilibrium dispersion velocity, u, that is comparable to the escape velocity for the object size that contains most of the swarm's mass. The accretion timescale for a radius R body (or alternatively, the time spent at radius R) is approximately ${\tau }_{\mathrm{acc}}\sim {(\rho R/\sigma )({f}_{g}{\rm{\Omega }})}^{-1}$, where ρ is the body density and f_g is a gravitational focusing factor that is a function of ${({v}_{\mathrm{esc}}/u)}^{2}$, where v_esc is the escape velocity of the growing body and u is the dispersion velocity of the accreted material. If accreting objects are similar in size, or if the largest objects contain most of the system mass (as is true in the end stages of accretion), then u ≥ v_esc and f_g is of order unity, implying growth of a R ∼ 750 km, ice–rock satellite at 15R_U in ∼500 yr, consistent with our simulation results. Earlier growth could be much faster if a larger body is accreting smaller material and most of the swarm mass is contained in the smaller material, because (v_esc/u) and f_g can then be large. Our simulations are limited in their ability to resolve such effects by their numerical resolution.

Zoom In Zoom Out Reset image size

After 3 × 10⁴ yr, the simulation in Figure 3 (left) obtains a satellite system broadly similar to Ariel, Umbriel, Titania, and Oberon. However there is no low-mass Miranda analog, perhaps pointing to a different origin for this moon (see below). The outer debris disk accretion process is completed within a few 10⁴ yr. Properties of the final large satellites, defined as those having a mass greater than 10⁻⁵ M_U, are shown in the right portion of Table 2. The average final number of large satellites and total satellite system mass across the simulations are 〈N_sats〉 = 4.24 ± 0.9 and 〈M_sats/M_AUTO〉 = 0.95 ± 0.09, in good agreement with the current system. Disks with an initial outer edge equal to Oberon's current distance of 23R_U (Runs 1–14) produce systems that are on average too radially compact, with an outermost moon well interior to Oberon's current distance ($\langle {a}_{\max }/{a}_{\mathrm{Oberon}}\rangle =0.50$). Increasing the initial outer edge to 30R_U (Table 2, Runs 15–22) improves this result, with $\langle {a}_{\max }/{a}_{\mathrm{Oberon}}\rangle =0.77$. For both (R_out/R_U) = 23 and 30, the final systems are more compact as the initial i is increased.

The latter effect is simply understood. As an initially thick, high-i torus collisionally damps to form a flattened disk and satellites, the components of its particles' orbital angular momentum in the equatorial plane will tend to cancel out, while the components perpendicular to this plane will be approximately conserved. The latter would be ${h}_{\mathrm{perp}}\approx \sqrt{{{GM}}_{{\rm{U}}}{a}_{o}}\cos i$ for initial debris with inclination i, small eccentricity, and semimajor axis a_o . Conservation of this quantity as material collisionally damps to low-i orbits implies contraction to a distance ${a}_{f}\sim {a}_{o}{(\cos i)}^{2}$. Requiring a final maximum semimajor axis consistent with Oberon's distance, i.e., a_f ≈ 25R_U, then implies that an initial disk with a substantially larger maximum semimajor axis, ${R}_{\mathrm{out}}\sim 25{R}_{{\rm{U}}}/{(\cos i)}^{2}$, is needed, a condition most closely met by Runs 5–6, 11–14, and 19–22, whose results are (generally) the most consistent with the current Uranian satellites.

This constraint has implications for the needed giant impact configuration. Just after the giant impact and prior to substantial collisional inclination damping, the outer debris disk inclination relative to Uranus' new, post-impact equatorial plane would have been ≈ δ, the angle between Uranus' pre-impact and post-impact spin axes. The initial R_out value for this debris would have been comparable to the outer radius of the preexisting regular satellite system formed by co-accretion. Requiring a_f ≈ 25R_U to account for Oberon then constrains δ, with $\delta ={\cos }^{-1}[{(25{R}_{{\rm{U}}}/{R}_{\mathrm{out}})}^{1/2}]$. It seems reasonable to assume that a preexisting Uranian satellite system formed via co-accretion would have had a broadly comparable radial scale to that of the Jovian and Saturnian regular satellites, with an outermost large satellite interior to 100 planetary radii. For R_out ≤ 100R_U, the maximum allowable value for δ is then 60° if one requires a_f ≈ 25R_U. Together with the requirement that the giant impact leave Uranus with its current 98° obliquity, this means that Uranus' obliquity prior to the giant impact in the Morbidelli et al. (2012) scenario must have been substantial, with θ₀ ∼ 40° or greater. This may not be implausible, given Neptune's 30° obliquity and the possibility for multiple large impacts and/or spin–orbit resonant effects (Rogoszinski & Hamilton 2020), but it is more restrictive than the arguments advanced in Morbidelli et al. (2012), which considered only the δ < 90° constraint cited in Section 2.1 that allows for smaller values of θ₀.

The Table 2 simulations produced on average an innermost large satellite well interior to Ariel, with $\langle {a}_{\min }/{R}_{{\rm{U}}}\rangle =5.4\pm 2.8$ versus (a_Ariel/R_U) = 7.5. Our first suite of simulations retained all objects that avoided direct collision with the planet. However, objects that strayed within a few R_U may instead have been lost due to gas drag by the inner water vapor c-disk. A c-disk with mass 10⁻² M_U that extends from the planet's surface to 3R_U has a surface density σ ∼ 6 × 10⁶ g cm^–2. The lifetime of a satellite with radius R and density ρ orbiting within such a disk is

$$\begin{eqnarray}&&{\tau }_{\mathrm{gd}}\sim \displaystyle \frac{8}{3{C}_{D}}\displaystyle \frac{\rho R}{\sigma }{\left(\displaystyle \frac{r{\rm{\Omega }}}{c}\right)}^{3}{{\rm{\Omega }}}^{-1},\end{eqnarray} \tag{ 2 }$$

where C_D ∼ O(1) is a drag constant, and c/(rΩ) is the vapor scale height of the disk, where c is the vapor sound speed, which is few × 10⁵ cm s⁻¹ soon after the impact (Rufu & Canup 2020). Near 3R_U, the loss timescale may then be only τ_gd ∼ 10 yr (1/C_D ) (6 × 10⁶ g cm⁻²/σ)(R/500 km)(0.1/[c/rΩ])³.

As such, we repeated a subset of the runs with the condition that any object that strayed within 3R_U was removed; results are shown in Table 3 and in Figure 3 (right). The average number and total mass of large satellites in these runs are 〈N_sats〉 = 3.47 ± 1.2 and 〈M_sats/M_AUTO〉 = 0.79 ± 0.28, while the average semimajor axis of the innermost final large satellite is increased to $\langle {a}_{\min }/{R}_{{\rm{U}}}\rangle =5.69\pm 1.58$, in somewhat better agreement with Ariel. However, we still do not see Miranda analogs: the final innermost moons in these simulations all have masses more than 15 times that of Miranda.

Table 3.  N-body Simulations of Outer Debris Disk Accretion with Removal Inside 3R_U

Run	Nsats	Msats	amin	amax	〈esats〉	〈isats〉
		$\left({M}_{\mathrm{AUTO}}\right)$	(RU)	$\left({a}_{\mathrm{Oberon}}\right)$		(deg)
1	4	1.02	4.88	0.48	0.017 ± 0.004	1.7 ± 1.02
2	4	0.90	4.44	0.42	0.040 ± 0.027	1.6 ± 1.14
3	2	0.47	4.33	0.23	0.012 ± 0.003	3.9 ± 0.10
4	2	0.33	4.32	0.22	0.015 ± 0.007	3.1 ± 0.02
6	5	0.96	5.46	0.74	0.012 ± 0.005	1.2 ± 1.03
9	2	0.44	4.23	0.23	0.036 ± 0.022	2.6 ± 0.72
10	2	0.30	4.35	0.22	0.010 ± 0.002	3.9 ± 0.33
13	4	1.02	8.25	0.98	0.037 ± 0.018	1.0 ± 0.54
14	4	1.10	4.81	0.98	0.047 ± 0.037	1.9 ± 1.23
15	4	0.93	7.55	0.70	0.025 ± 0.008	2.7 ± 0.54
16	4	0.85	5.46	0.53	0.017 ± 0.015	2.7 ± 1.04
17	2	0.62	4.89	0.26	0.026 ± 0.005	5.0 ± 0.44
18	2	0.44	4.87	0.26	0.015 ± 0.007	4.2 ± 0.56
19	4	0.97	7.95	0.87	0.028 ± 0.019	1.5 ± 0.87
20	5	0.98	5.85	0.94	0.019 ± 0.010	1.8 ± 0.69
21	4	1.03	9.26	1.33	0.033 ± 0.018	1.0 ± 0.45
22	5	1.06	5.80	1.22	0.057 ± 0.033	2.5 ± 2.05

Note. Results from a subset of the runs from Table 2 where we have this time removed any particle that passes within 3R_U to mimic removal via gas drag from the vapor-rich c-disk.

Download table as:  ASCIITypeset image

Inner Miranda is distinct not just because of its much smaller mass compared with the outer large moons, but also because of its apparent composition. The large moons have similar densities between 1.52 and 1.66 g cm⁻³ (Jacobson 2014), implying similar bulk compositions with ≥50% rock by mass, consistent with material expected in a co-accretion disk. In contrast, Miranda's density is 1.17 g cm⁻³. Miranda has Enceladus-like tectonic features, attributed to upwelling of partially melted ices due to past tidal heating (e.g., Pappalardo & Schubert 2013; Beddingfield et al. 2015), and such endogenic activity seems inconsistent with preservation of large-scale porosity in its interior. Thus, Miranda's density seems to imply a much higher ice content than in the large moons, which is difficult to reconcile with an origin from the reaccumulation of material produced during earlier co-accretion. Instead, it is possible Miranda originated from material in the ice-rich c-disk, perhaps consistent with the model of Hesselbrock & Minton (2019); we return to this issue in Section 3.5.

With these caveats, we conclude that reassembly of a disrupted prior satellite system formed by co-accretion could plausibly produce a satellite system resembling the current four large Uranian moons, so long as the pre-impact Uranian obliquity was substantial (θ₀ ≥ 40°). For the remainder of the paper we focus on what we find to be the more constraining final phase of the co-accretion + giant impact model, in which the inner massive c-disk cools and viscously evolves.

We simulate the c-disk using a model developed in the context of the Moon's accretion after a giant impact (Salmon & Canup 2012). Our model represents material within the Roche limit by a uniform surface density disk that is described analytically. Material outside the Roche limit is described by an N-body code (Duncan et al. 1998) modified to include tidal accretion criteria relevant near the Roche limit (Canup & Esposito 1995, 1996). The Roche-interior disk's total mass (M_d ) and its outer edge (R_d ≤ a_R ) evolve with time due to a radiation-limited viscosity per above and interactions with outer moons. Material that spreads inward onto the planet is removed. As material spreads outward past the Roche limit, we remove it from the continuum disk portion of the model and add it to the N-body code in the form of new moonlets just exterior to the Roche limit. We include the strongest resonant interactions (i.e., the 2:1, 3:2, etc.) for all moons close enough to the disk to have one or more of their strong resonances fall within the disk, i.e., for all moons with a semimajor axis a < 1.6R_d , where R_d is the outer radius of the Roche-interior disk. Objects passing close enough to the planet to tidally disrupt (Sridhar & Tremaine 1992) are removed from the N-body code and their mass and angular momentum are added to that of the Roche-interior disk. We consider tidal time delay values Δt = 2.7 or 270 s, corresponding respectively to $\left(Q/{k}_{2}\right)\approx {10}^{4}$ or ≈10² at a distance of 3.8R_U.

We perform simulations (Table 4) with three different c-disk masses: M_d = 3 × 10⁻³, 10⁻², and 3 × 10⁻² M_U. The smallest of these appears insufficient to realign the outer disk to distances consistent with Oberon's orbit (Sections 2.1 and 2.3), but is included here for comparison's sake. For each disk mass, we perform three simulations with different random values for the longitudes of ascending node and mean anomalies of spawned moonlets.

Table 4. Simulations of the c-disk

Run	Mc‐disk	(Q/k2)	Msats > async
	$\left({M}_{{\rm{U}}}\right)$		$\left({M}_{\mathrm{AUTO}}\right)$
1	0.003	104	7.45
2	0.003	104	6.56
3	0.003	104	6.99
4	0.010	104	15.19
5	0.010	104	15.78
6	0.010	104	15.09
7	0.030	104	35.20
8	0.030	104	40.59
9	0.030	104	24.07
10	0.003	100	6.88
11	0.003	100	7.03
12	0.003	100	7.50
13	0.010	100	14.53
14	0.010	100	13.63
15	0.010	100	17.53
16	0.030	100	1.00
17	0.030	100	1.00
18	0.030	100	1.00

Note. Initial parameters and results from simulations of the evolution of a massive c-disk in the presence of analogs for Ariel, Umbriel, Titania, and Oberon (Figure 4). M_c‐disk is the mass of the initial Roche-interior c-disk, (Q/k₂) are Uranus' tidal parameters (see text), and M_sats > a_sync is the total mass of moons beyond a_sync at the end of the simulation, in units of the combined mass of the current Uranian moons M_AUTO. In most cases (Runs 1–15), massive moons spawned from the c-disk expand beyond a_sync and disrupt/absorb the outer moons. Only in the case of a very massive c-disk and very strong tidal dissipation in Uranus (Runs 16–18) is the outer moon system retained.

Download table as:  ASCIITypeset image

We assume a compact c-disk that initially lies entirely within the Roche limit for material with density ≈1 g cm⁻³. A more extended c-disk is certainly plausible (and perhaps probable), but the compact case appears the most likely to produce a successful outcome in which c-disk material remains interior to synchronous orbit. We consider a non-fully contracted Uranus with a radius of R_P = 1.3R_U and rotation period of ∼28 hr, so that a_sync is shifted outward to 4.5R_U. Finally, we place the four largest moons at their current positions, assuming that they accreted at these locations during the prior debris disk accretion phase (Figure 4).

Zoom In Zoom Out Reset image size

Figure 5 shows the evolution of a system with a c-disk of mass M_d = 0.003M_U. The disk spawns an initial moonlet (orange dot) whose orbit rapidly expands because the resonant interactions with the disk are stronger than the tidal torque from the planet. The moonlet first ceases its outward migration at ∼ 4.25R_U where its 2:1 resonance lies at the Roche limit. Shortly thereafter, the disk spawns additional moonlets (blue dot). These mostly collide with the first moonlet leading to its growth in mass. However, as some of them get scattered or trapped into mean-motion resonances, they transfer additional angular momentum to the first moonlet whose orbit continues to expand until it crosses the synchronous orbit. Once this first moonlet has receded away sufficiently, a second moonlet can grow and evolve outward through the same process. These accretion dynamics are similar to the continuous, discrete, and pyramidal regimes described in the analytical work of Crida & Charnoz (2012).

Zoom In Zoom Out Reset image size

After ∼1.5 × 10⁵ yr, the first moonlet spawned from the c-disk absorbs Ariel. After ∼4.2 × 10⁵ yr, it absorbs Umbriel. After 6 × 10⁵ yr, Titania and Oberon have been mostly unaffected, but two very massive interior moons have been brought beyond the synchronous orbit. These would be long lived, and yield a total satellite system mass ∼7 times greater than the current system, an unsuccessful result.

Figure 6 shows the evolution of a system with a high-mass c-disk with M_d = 3 × 10⁻² M_U. The disk again spawns a first moonlet (orange dot) whose orbit rapidly expands because the positive resonant toque is much stronger than the negative tidal torque from the planet. However, in this case the first spawned moonlet is massive enough to temporarily confine the c-disk inside the Roche limit, so that for a time R_d < a_R . This delays the time until the next moonlet is spawned by the time required for the c-disk to viscously spread back to the Roche limit. After ∼2 × 10⁵ yr, a second moonlet is spawned at the Roche limit and is immediately absorbed by the first one. This repeats twice, until the original spawned moonlet becomes so massive that it perturbs the orbit of Ariel, which is eventually absorbed at ∼4.7 × 10⁵ yr. The dynamical exchange results in the moonlet's new semimajor axis expanding to outside synchronous orbit. After ∼4.2 × 10⁵ yr, the outer three moons go through an instability due to the inner moonlet crossing their mean-motion resonances. This results in Umbriel being ejected from the system.

Zoom In Zoom Out Reset image size

In this case Titania and Oberon survive at ∼6 × 10⁵ yr, but their orbits have been strongly affected, in particular their eccentricities. Here again, a moon whose mass is ∼ 10M_AUTO and is comprised primarily of c-disk material has been brought beyond the synchronous orbit, at odds with the Uranus system.

In the previous case we found that a more massive disk can spawn a moon massive enough to confine the disk inside the Roche limit. This limits the number of moonlets produced, resulting in reduced scattering and angular momentum transfer, which keeps the outer system stable for more than ∼4 × 10⁵ yr. However, with weak tidal dissipation in Uranus the large spawned moonlet remains near a_sync and does not tidally evolve inward before mutual interactions drive its orbit beyond a_sync. As such, we performed a second suite of simulations with the same setup, but with $\left(Q/{k}_{2}\right)={10}^{2}$, which produces 100 times faster tidal evolution.

Figure 7 shows the evolution of a system with a high-mass c-disk (M_d = 3 × 10⁻² M_U) with $\left(Q/{k}_{2}\right)={10}^{2}$. A massive moonlet is spawned early on and confines the disk inside the Roche limit while its orbit again rapidly expands to ∼ 4.25R_U. But now the tidal torque is sufficient to cause the moonlet to spiral inward before the inner disk can viscously expand back out to a_R and spawn additional moonlets. As the moonlet's orbit contracts, it drives the c-disk toward the planet as well, allowing the outer four satellites to be retained.

Zoom In Zoom Out Reset image size

As expected based on prior works (Crida & Charnoz 2012; Salmon & Canup 2012; Hesselbrock & Minton 2017, 2019; Canup & Salmon 2018), the most massive moon spawned by a Roche-interior c-disk initially forms near the Roche limit and rapidly recoils via disk torques to a distance ∼1.6r_d , with r_d < a_R for the initial moon spawned from a sufficiently massive c-disk. Subsequently the disk's edge viscously spreads back out to a_R on timescale τ_ν , which tends to drives the moon outward on this timescale too. The competing effect is tidal interaction with the planet, which so long as the moon remains interior to a_sync causes the moon's orbit to lose angular momentum and spiral inward on timescale τ_tidal. For the protolunar disk, τ_ν ≪ τ_tidal and disk accretion produces a moon whose semimajor axis prior to any tidal evolution is ∼2.2 times the Roche limit (Salmon & Canup 2012). This would here imply a massive moon that recoils out to ∼6R_U and survives, an unsuccessful result. However, in a water-dominated c-disk, the viscous spreading timescale is orders-of-magnitude longer, and it is possible for very low (Q/k₂) to instead have τ_tidal < τ_ν . In this case, inward tidal evolution decreases the maximum distance obtained by a spawned moon, keeping it within a_sync. Very strong tidal dissipation in early Uranus then provides for a potentially successful outcome, as seen in Figure 7.

In actuality, an inwardly decaying spawned moon will tidally disrupt and resupply the Roche-interior disk, an effect not included in our simulations. On longer timescales the c-disk will continue to viscously spread, lose mass, and spawn ever-smaller inner moons that may be increasingly likely to stay within a_sync (Hesselbrock & Minton 2017, 2019). Hesselbrock & Minton (2019) proposed that Miranda (as well as smaller more interior Uranian moons) is a spawned moonlet from a Roche-interior disk with a very low mass, ∼few × 10⁻⁶ M_U. Whether Miranda could originate from the initially vastly more massive c-disk of the Morbidelli et al. (2012) model is an intriguing question. If Miranda did originate from c-disk material, while the 4 large outer moons instead were re-assembled from a prior co-accreted system, this would offer a natural explanation for why Miranda is ice-rich while the outer large moons are ≥50% rock.

Could Miranda have accreted from the outer portions of a massive c-disk if a low-mass component of the disk initially extended to ∼5R_U? For this to be viable, low-mass Miranda (with mass 10⁻⁶ M_U) at ∼5R_U must avoid being accreted by the much more massive moon spawned at the Roche limit as the Roche-interior c-disk viscously expanded. The simulations shown in Figures 5–7 demonstrate that the first, most massive spawned moonlet has a mass between ∼10⁻⁴ and 10⁻³ M_U for the Morbidelli et al. (2012) c-disk, and that this moonlet rapidly expands outward to 4.25R_U due to disk torques. At this distance, the spawned moonlet would likely dynamically destabilize or accrete Miranda because they would be separated by less than several mutual Hill radii. Broadly analogous results were seen in the Canup & Salmon (2018) simulations of the evolution of an impact-generated disk around Mars: massive Roche-interior disks produced massive spawned moons that accreted and destroyed outer smaller Deimos analogs unless the initial Roche-interior disk mass was below a critical value. Thus, it seems unlikely that Miranda could accrete from the initial outer portions of a massive c-disk and survive.

The second possibility is that Miranda was a spawned moon from the c-disk at a much later stage in its evolution when viscous spreading had reduced its mass by orders-of-magnitude to a few × 10⁻⁶ to 10⁻⁵ M_U. As the c-disk evolves, it eventually cools sufficiently to completely condense. Once this occurs, the disk's viscous timescale becomes inversely proportional to the square of the disk surface density (Ward & Cameron 1978; Salmon & Canup 2012) and it spreads more and more slowly as its mass progressively decreases. Our code is too computationally expensive to model this protracted spreading evolution. Future simulations using a more dynamically simplified, but computationally efficient model such as that of Hesselbrock & Minton (2019) will be needed to assess whether an initial c-disk massive enough for the Morbidelli et al. (2012) model may, much later in its evolution, spawn a stable Miranda-like moon.

We considered a pure water c-disk because impact simulations indicate that this would be the dominant component for a differentiated rock–ice impactor (Slattery et al. 1992; Reinhardt et al. 2020; Rufu & Canup 2020). The c-disk may also initially contain a minority rock component originating from the impactor's core, which would condense before the water. Its evolution should be less important to the viability of a co-accretion + giant impact model than the evolution of the water because the total mass in rock will likely be much less and because the rock Roche limit is closer to the planet (<2R_U), so that rocky moonlets will be more likely than their icy counterparts to remain interior to synchronous orbit and be lost. The c-disk might also contain some gas (H, He) from the outer layers of Uranus that could affect the late inner disk evolution once its water vapor has fully condensed, perhaps via gas drag, before the H-rich component was removed (e.g., via photoevaporation).