Tidal Evolution of the Eccentric Moon around Dwarf Planet (225088) Gonggong

We define γ_X,G ≡ M_X/M_G as the secondary-to-primary mass ratio of the system formed after giant impacts: satellite systems with 10⁻³ ≲ γ_X,G ≲ 10⁻¹ can be directly formed via giant impacts as intact fragments (Canup 2005; Arakawa et al. 2019). We also define γ_i,t ≡ M_imp/M_tar as the impactor-to-target mass ratio, where M_tar and M_imp are the masses of the target and impactor, respectively. Arakawa et al. (2019) performed giant impact simulations with three settings for the targets and impactors: γ_i,t = 0.5 with an ice composition (Set A), γ_i,t = 0.2 with an ice composition (Set B), and γ_i,t = 0.5 with a basalt composition (Set C). The total mass of the target and the impactor is M_tar + M_imp = 6 × 10²¹ kg in the simulations.

Figure 1(a) shows the distribution of the semilatus rectum from the results of the giant impact simulations (Arakawa et al. 2019). The semilatus rectum, p_orb, and the periapsis distance, q_orb, are defined as follows (Murray & Dermott 1999):

$$\begin{eqnarray}&&{p}_{\mathrm{orb}}\equiv a\left(1-{e}^{2}\right)={q}_{\mathrm{orb}}\left(1+e\right),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{q}_{\mathrm{orb}}\equiv a\left(1-e\right),\end{eqnarray} \tag{ 2 }$$

where a is the semimajor axis and e is the eccentricity. Arakawa et al. (2019) performed more than 400 runs of giant impact simulations of 1000 km sized TNOs, and revealed that the initial periapsis distance of satellites formed via giant impacts, q_ini, is in the range

$$\begin{eqnarray}&&3\lesssim \displaystyle \frac{{q}_{\mathrm{ini}}}{{R}_{{\rm{G}}}}\lesssim 4,\end{eqnarray} \tag{ 3 }$$

where R_G is the radius of the primary. As the initial eccentricity of satellites formed via giant impacts, e_ini, is widely distributed in the range 0–1 (Arakawa et al. 2019), the initial semilatus rectum, p_ini, should be in the range

$$\begin{eqnarray}&&3\lesssim \displaystyle \frac{{p}_{\mathrm{ini}}}{{R}_{{\rm{G}}}}\lesssim 8.\end{eqnarray} \tag{ 4 }$$

We confirmed that p_ini/R_G is in the range 3–8 for all outcomes resulting in satellite formation, as predicted in Equation (4).

Zoom In Zoom Out Reset image size

Figure 1(b) shows the distribution of the initial spin period of secondaries after giant impacts, ${P}_{{\rm{X}},\mathrm{ini}}=2\pi /{\dot{\theta }}_{{\rm{X}},\mathrm{ini}}$, where ${\dot{\theta }}_{{\rm{X}},\mathrm{ini}}$ is the initial spin angular velocity of secondaries. Although it is difficult to analytically estimate P_X,ini, the numerical results show that the initial spin period of the secondaries is in the range 3 hr ≲ P_X,ini ≲ 18 hr. The lower limit of P_X,ini is given by P_X,ini > 2π/Ω_cr, where

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{cr}}\equiv \sqrt{\displaystyle \frac{4\pi { \mathcal G }\rho }{3}},\end{eqnarray} \tag{ 5 }$$

is the critical spin angular velocity for rotational instability (e.g., Kokubo & Genda 2010); where ${ \mathcal G }$ is the gravitational constant, and ρ is the bulk density.

In this section, we discuss the initial spins and orbit of the Gonggong–Xiangliu system from the perspective of the conservation of angular momentum. The total angular momentum of the system, L_tot, is

$$\begin{eqnarray}&&{L}_{\mathrm{tot}}={I}_{{\rm{G}}}{\dot{\theta }}_{{\rm{G}}}+{I}_{{\rm{X}}}{\dot{\theta }}_{{\rm{X}}}+\displaystyle \frac{{M}_{{\rm{G}}}{M}_{{\rm{X}}}}{{M}_{\mathrm{tot}}}\sqrt{{ \mathcal G }{M}_{\mathrm{tot}}{p}_{\mathrm{orb}}},\end{eqnarray} \tag{ 6 }$$

and L_tot should be constant over time during the tidal evolution of the system. Here, M_tot ≡ M_G + M_X is the total mass of the satellite system, and I_i and ${\dot{\theta }}_{i}$ are the moment of inertia and the spin angular velocity, respectively. Hereafter, the subscript i denotes the primary (i = G) or the secondary (i = X). For the case of undifferentiated bodies, the density distribution of the planetary interior is homogeneous, and the moment of inertia is approximately given by

$$\begin{eqnarray}&&{I}_{i}=\displaystyle \frac{2}{5}{M}_{i}{R}_{i}^{2}.\end{eqnarray} \tag{ 7 }$$

As the spin angular momentum of the secondary, Xiangliu, is negligibly small, we obtain the following equation:

$$\begin{eqnarray}&&\begin{array}{l}{I}_{{\rm{G}}}{\dot{\theta }}_{{\rm{G}},\mathrm{ini}}+\displaystyle \frac{{M}_{{\rm{G}}}{M}_{{\rm{X}}}}{{M}_{\mathrm{tot}}}\sqrt{{ \mathcal G }{M}_{\mathrm{tot}}{p}_{\mathrm{ini}}}\\ \ \ ={I}_{{\rm{G}}}{\dot{\theta }}_{{\rm{G}},\mathrm{obs}}+\displaystyle \frac{{M}_{{\rm{G}}}{M}_{{\rm{X}}}}{{M}_{\mathrm{tot}}}\sqrt{{ \mathcal G }{M}_{\mathrm{tot}}{p}_{\mathrm{obs}}},\end{array}\end{eqnarray} \tag{ 8 }$$

where ${\dot{\theta }}_{{\rm{G}},\mathrm{ini}}$ and ${\dot{\theta }}_{{\rm{G}},\mathrm{obs}}$ are the initial and observed (current) spin angular velocities of the primary, Gonggong; and p_ini and p_obs are the initial and observed semilatus recta. Figure 2 shows the initial semilatus rectum, p_ini, as a function of the initial spin period of Gonggong, P_G,ini, the observed spin period of Gonggong, P_G,obs, and the secondary-to-primary mass ratio, γ_X,G.

Zoom In Zoom Out Reset image size

We note that the current spin period of Gonggong has not yet been uniquely determined; the observation by Pál et al. (2016) suggested P_G,obs = 22.4 hr or 44.8 hr, although the most likely light-curve solution is the double-peaked solution with a slight asymmetry. This corresponds to a spin period of P_G,obs =44.8 hr. In contrast, Kiss et al. (2019) suggested that the most plausible solution for the system would be a single-peaked light curve in the visible range (i.e., P_G,obs = 22.4 hr) caused by surface features rather than by a distorted shape, and the primary's equator and the secondary's orbit are coplanar.

We also calculated the initial location of the corotation radius, a_corot, which is given by

$$\begin{eqnarray}&&{a}_{\mathrm{corot}}={\left(\displaystyle \frac{\sqrt{{ \mathcal G }{M}_{\mathrm{tot}}}}{2\pi }{P}_{{\rm{G}},\mathrm{ini}}\right)}^{2/3}.\end{eqnarray} \tag{ 9 }$$

For an initial eccentricity of e_ini = 0, the tidal evolution of the satellite results in inward migration when the initial semimajor axis is a_ini < a_corot. Then, the initial semilatus rectum is expected to be p_ini ≳ a_corot to explain the current semimajor axis of the system, a_obs, which is an order of magnitude larger than a_corot.

Figure 2 indicates that the tidal evolution of the satellite results in outward migration when the secondary-to-primary mass ratio is approximately γ_X,G = 10⁻² or larger. In contrast, the tidal evolution of the satellite would result in inward migration when γ_X,G ≪ 10⁻² and p_ini/R_G ≲ 4 (i.e., e_ini ≃ 0). Therefore, for the case of γ_X,G ≪ 10⁻², the initial orbit of the Gonggong–Xiangliu system formed via a giant impact must have a finite eccentricity; otherwise, the satellite would migrate inward and become tidally disrupted. The minimum value of e_ini (and p_ini) is expected to depend on γ_X,G and P_G,ini; we confirm this prediction in Sections 5 and 6.

We summarize the initial parameters for the calculations of tidal evolution in Table 2. In Section 5, we parameterize the radius of Xiangliu, R_X, the initial eccentricity, e_ini, and the initial temperature, T_ini. For simplicity, we set ρ_G = ρ_X, where ρ_G and ρ_X are the bulk densities of Gonggong and Xiangliu. We assume that the initial temperature of Xiangliu is equal to that of Gonggong (i.e., T_G,ini = T_X,ini = T_ini) in Section 5. We also investigated the dependence of the outcomes of tidal evolution on the current spin period of Gonggong, P_G,obs, the initial spin period of Xiangliu, P_X,ini, the different settings of the initial temperature of Xiangliu (we set T_X,ini = 120 K regardless of T_G,ini), and the value of the reference viscosity (see Table 3), η_ref, in Section 6.

Table 2. Initial Condition for Calculations of Tidal Evolution

Symbol	Parameter	Value	Unit	Comment
qini	Initial periapsis distance	2.1 × 103	km	3 ≲ qini/RG ≲ 4 (Arakawa et al. 2019)
eini	Initial eccentricity	0.1, 0.2, ..., 0.8	⋯	0 < eini < 1 (Arakawa et al. 2019)
PG,ini	Initial spin period of Gonggong	⋯	hr	${P}_{{\rm{G}},\mathrm{ini}}=2\pi /{\dot{\theta }}_{{\rm{G}},\mathrm{ini}}$ is given by Equation (8)
PX,ini	Initial spin period of Xiangliu	12 or Porb,ini	hr	3 hr < PX,ini < 18 hr (see Figure 1(b))
RG	Radius of Gonggong	600	km	Kiss et al. (2019)
RX	Radius of Xiangliu	20, 40, ..., 120	km	RX ≳ 18 km (see Kiss et al. 2019)
TG,ini	Initial temperature of Gonggong	120, 140, ..., 240	K	TG,ini is below the melting point of ice
TX,ini	Initial temperature of Xiangliu	T G,ini or 120	K	See Section 6.3

Note. Parameters for standard runs are indicated by bold face.

Download table as:  ASCIITypeset image

Table 3. Material Properties Used in This Study

Symbol	Property	Value	Unit	References
c	Specific heat capacity	$8.8\left(T/{\rm{K}}\right)$	J kg−1 K−1	Hammond et al. (2016)
kth	Thermal conductivity	$0.48+488{\left(T/{\rm{K}}\right)}^{-1}$	W m−1 K−1	Hammond et al. (2016)
αexp	Thermal expansion coefficient	1 × 10−4	K−1	Kamata et al. (2019)
Ea	Activation energy	60	kJ mol−1	Kamata et al. (2019)
Tref	Reference temperature	273	K	Kamata et al. (2019)
ηref	Reference viscosity at Tref	1014 or 1010	Pa s	See Section 6.4
μ	Shear modulus	3.33	GPa	Robuchon & Nimmo (2011)
α	Andrade exponent	0.33	⋯	Rambaux et al. (2010)

Note. Parameters for standard runs are indicated by bold face.

Download table as:  ASCIITypeset image

In this study, we assume that the primary's equator and the secondary's orbit are coplanar, and the orbit of the satellite system is prograde (i.e., the spin–orbit angle is close to zero). This assumption is also consistent with the numerical simulations of Arakawa et al. (2019); however, their simulations do not consider the pre-impact spin of the target and impactor. Observations by Kiss et al. (2019) are consistent with both prograde and retrograde orbit. If the Gonggong–Xiangliu system is in retrograde orbit, the semimajor axis will decrease with increasing time. Thus the determination of the spin–orbit angle by future observations could provide critical constraints on the origin of the satellite system.

We also note that the observed coplanar orbit is for the present-day system, and this does not indicate that the system was born in a coplanar orbit. Although the giant impact simulations preferred a coplanar orbit, we should discuss the effect of nonzero initial inclination/obliquity on the tidal evolution in future studies. Several studies (e.g., Correia 2020; Renaud et al. 2021) pointed out that obliquity evolution can affect how bodies fall into and out of higher-order spin–orbit resonances.

Tides are raised on both the primary and secondary. The tidal lag caused by friction leads to angular momentum exchange, which also leads to spin and orbital evolution. In this study, we use the tidal evolution equations following Boué & Efroimsky (2019). The orbit-averaged variations of the spin rates of the primary and secondary, ${\dot{\theta }}_{{\rm{G}}}$ and ${\dot{\theta }}_{{\rm{X}}}$, semimajor axis a, and eccentricity e are given by

$$\begin{eqnarray}&&\displaystyle \frac{1}{n}\displaystyle \frac{d{\dot{\theta }}_{{\rm{G}}}}{{dt}}=-\displaystyle \frac{3n}{2}\displaystyle \frac{{M}_{{\rm{G}}}{M}_{{\rm{X}}}}{{M}_{{\rm{tot}}}}\displaystyle \frac{{M}_{{\rm{X}}}}{{M}_{{\rm{G}}}}\displaystyle \frac{{R}_{{\rm{G}}}^{2}}{{I}_{{\rm{G}}}}{\left(\displaystyle \frac{{R}_{{\rm{G}}}}{a}\right)}^{3}\sum _{q}{{ \mathcal A }}_{{\rm{G}},q},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{n}\displaystyle \frac{d{\dot{\theta }}_{{\rm{X}}}}{{dt}}=-\displaystyle \frac{3n}{2}\displaystyle \frac{{M}_{{\rm{G}}}{M}_{{\rm{X}}}}{{M}_{{\rm{tot}}}}\displaystyle \frac{{M}_{{\rm{G}}}}{{M}_{{\rm{X}}}}\displaystyle \frac{{R}_{{\rm{X}}}^{2}}{{I}_{{\rm{X}}}}{\left(\displaystyle \frac{{R}_{{\rm{X}}}}{a}\right)}^{3}\sum _{q}{{ \mathcal A }}_{{\rm{X}},q},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{a}\displaystyle \frac{{da}}{{dt}}=\displaystyle \frac{1}{a}\left[{\left(\displaystyle \frac{{da}}{{dt}}\right)}_{{\rm{G}}}+{\left(\displaystyle \frac{{da}}{{dt}}\right)}_{{\rm{X}}}\right],\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{r}\begin{array}{c}\displaystyle \frac{1}{a}{\left(\displaystyle \frac{{da}}{{dt}}\right)}_{{\rm{G}}}=2n\displaystyle \frac{{M}_{{\rm{X}}}}{{M}_{{\rm{G}}}}{\left(\displaystyle \frac{{R}_{{\rm{G}}}}{a}\right)}^{5}\\ \,\times \,\displaystyle \sum _{q}\left[\displaystyle \frac{\left.3(2+q\right)}{4}{{ \mathcal A }}_{{\rm{G}},q}+\displaystyle \frac{q}{4}{{ \mathcal B }}_{{\rm{G}},q}\right],\end{array}\end{array}\end{array}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}\displaystyle \frac{1}{a}{\left(\displaystyle \frac{{da}}{{dt}}\right)}_{{\rm{X}}}=2n\displaystyle \frac{{M}_{{\rm{G}}}}{{M}_{{\rm{X}}}}{\left(\displaystyle \frac{{R}_{{\rm{X}}}}{a}\right)}^{5}\\ \,\times \,\displaystyle \sum _{q}\left[\displaystyle \frac{\left.3(2+q\right)}{4}{{ \mathcal A }}_{{\rm{X}},q}+\displaystyle \frac{q}{4}{{ \mathcal B }}_{{\rm{X}},q}\right],\end{array}\end{array}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{e}\displaystyle \frac{{de}}{{dt}}=\displaystyle \frac{1}{e}\left[{\left(\displaystyle \frac{{de}}{{dt}}\right)}_{{\rm{G}}}+{\left(\displaystyle \frac{{de}}{{dt}}\right)}_{{\rm{X}}}\right],\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\begin{array}{c}\frac{1}{e}{\left(\frac{{de}}{{dt}}\right)}_{{\rm{G}}}=\frac{n}{{e}^{2}}\frac{{M}_{{\rm{X}}}}{{M}_{{\rm{G}}}}{\left(\frac{{R}_{{\rm{G}}}}{a}\right)}^{5}\\ \times \displaystyle \sum _{q}\,\left[\frac{3(2+q)(1-{e}^{2})-6\sqrt{1-{e}^{2}}}{4}{{ \mathcal A }}_{{\rm{G}},q}+\frac{q}{4}(1-{e}^{2}){{ \mathcal B }}_{{\rm{G}},q}\right],\end{array}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\begin{array}{c}\frac{1}{e}{\left(\frac{{de}}{{dt}}\right)}_{{\rm{X}}}=\frac{n}{{e}^{2}}\frac{{M}_{{\rm{G}}}}{{M}_{{\rm{X}}}}{\left(\frac{{R}_{{\rm{X}}}}{a}\right)}^{5}\\ \times \displaystyle \sum _{q}\,\left[\frac{3(2+q)(1-{e}^{2})-6\sqrt{1-{e}^{2}}}{4}{{ \mathcal A }}_{{\rm{X}},q}+\frac{q}{4}(1-{e}^{2}){{ \mathcal B }}_{{\rm{X}},q}\right],\end{array}\end{eqnarray} \tag{ 17 }$$

where $n=\sqrt{{ \mathcal G }{M}_{\mathrm{tot}}/{a}^{3}}$ denotes the mean motion. The initial orbital period, P_orb,ini, is ${P}_{{\rm{orb}},{\rm{ini}}}=2\pi /\sqrt{{ \mathcal G }{M}_{{\rm{tot}}}/{a}_{{\rm{ini}}}^{3}}$. Here, we assume that the orbital inclinations on both the primary's and secondary's equators are negligibly small. We note that Kiss et al. (2019) suggested that the orbit of the secondary would be coplanar with the equator of the primary. The coefficients ${{ \mathcal A }}_{{\rm{G}},q}$, ${{ \mathcal A }}_{{\rm{X}},q}$, ${{ \mathcal B }}_{{\rm{G}},q}$, and ${{ \mathcal B }}_{{\rm{X}},q}$, are given by

$$\begin{eqnarray}&&{{ \mathcal A }}_{{\rm{G}},q}={\left[{G}_{\mathrm{2,0},q}(e)\right]}^{2}\ \mathrm{Im}\left[{\tilde{k}}_{2,{\rm{G}}}\left((2+q)n-2{\dot{\theta }}_{{\rm{G}}}\right)\right],\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{{ \mathcal A }}_{{\rm{X}},q}={\left[{G}_{\mathrm{2,0},q}(e)\right]}^{2}\ \mathrm{Im}\left[{\tilde{k}}_{2,{\rm{X}}}\left((2+q)n-2{\dot{\theta }}_{{\rm{X}}}\right)\right],\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{{ \mathcal B }}_{{\rm{G}},q}={\left[{G}_{\mathrm{2,1},q}(e)\right]}^{2}\ \mathrm{Im}\left[{\tilde{k}}_{2,{\rm{G}}}\left({qn}\right)\right],\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{{ \mathcal B }}_{{\rm{X}},q}={\left[{G}_{\mathrm{2,1},q}(e)\right]}^{2}\ \mathrm{Im}\left[{\tilde{k}}_{2,{\rm{X}}}\left({qn}\right)\right],\end{eqnarray} \tag{ 21 }$$

where G_2,0,q and G_2,1,q are the eccentricity functions (see Appendix A), and ${\tilde{k}}_{2,{\rm{G}}}$ and ${\tilde{k}}_{2,{\rm{X}}}$ are the Love numbers of Gonggong and Xiangliu, which depend on the tidal frequency (see Appendix B). We consider the terms G_2,p,q (p = 0 or 1) when the following two conditions, ∣q∣ ≤ 200 and ${[{G}_{2,p,q}]}^{2}\geqslant {10}^{-20}$, are satisfied.

In our numerical integration, we calculate ${\dot{\theta }}_{{\rm{G}}}/n$ and ${\dot{\theta }}_{{\rm{X}}}/n$ instead of ${\dot{\theta }}_{{\rm{G}}}$ and ${\dot{\theta }}_{{\rm{X}}}$ (e.g., Cheng et al. 2014). The relation between ${dn}/{dt}$ and da/dt is

$$\begin{eqnarray}&&\displaystyle \frac{1}{n}\displaystyle \frac{{dn}}{{dt}}=-\displaystyle \frac{3}{2}\displaystyle \frac{1}{a}\displaystyle \frac{{da}}{{dt}},\end{eqnarray} \tag{ 22 }$$

and we obtain the following equations:

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\left(\displaystyle \frac{{\dot{\theta }}_{{\rm{G}}}}{n}\right)=\displaystyle \frac{1}{n}\displaystyle \frac{d{\dot{\theta }}_{{\rm{G}}}}{{dt}}+\displaystyle \frac{3{\dot{\theta }}_{{\rm{G}}}}{2n}\displaystyle \frac{1}{a}\displaystyle \frac{{da}}{{dt}},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\left(\displaystyle \frac{{\dot{\theta }}_{{\rm{X}}}}{n}\right)=\displaystyle \frac{1}{n}\displaystyle \frac{d{\dot{\theta }}_{{\rm{X}}}}{{dt}}+\displaystyle \frac{3{\dot{\theta }}_{{\rm{X}}}}{2n}\displaystyle \frac{1}{a}\displaystyle \frac{{da}}{{dt}}.\end{eqnarray} \tag{ 24 }$$

In this study, we use the simplifying assumption of undifferentiated homogeneous bodies. We do not consider the internal temperature structure but calculate the effective temperature of the interior for simplicity ⁵ (e.g., Ojakangas & Stevenson 1986; Shoji et al. 2013). Then, the temperature evolution of the body i is given by the following equation:

$$\begin{eqnarray}&&\displaystyle \frac{{{dT}}_{i}}{{dt}}=\displaystyle \frac{{Q}_{\mathrm{tot},i}}{{M}_{i}{c}_{i}},\end{eqnarray} \tag{ 25 }$$

where c_i is the specific heat capacity of the primary and secondary. The total heat generation rate within bodies, Q_tot,i , is given by

$$\begin{eqnarray}&&{Q}_{\mathrm{tot},i}={Q}_{\mathrm{tide},i}+{Q}_{\mathrm{con},i}+{Q}_{\mathrm{dec},i},\end{eqnarray} \tag{ 26 }$$

where Q_tide,i , Q_con,i , and Q_dec,i are the tidal heating, the conduction/convection cooling, and the decay heating terms, respectively. We note that the main heat source is not tidal heating but decay heating. We discuss the balance between conduction/convection cooling and decay heating for the primary in Appendix C.

3.2.1. Tidal Heating

The tidal heating rate, Q_tide,i , is given by the sum of two terms:

$$\begin{eqnarray}&&{Q}_{\mathrm{tide},i}={Q}_{\mathrm{spin},i}+{Q}_{\mathrm{orb},i},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{Q}_{\mathrm{spin},i}=-{I}_{i}{\dot{\theta }}_{i}\displaystyle \frac{d{\dot{\theta }}_{i}}{{dt}},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{Q}_{\mathrm{orb},i}=-\displaystyle \frac{{{GM}}_{{\rm{p}}}{M}_{{\rm{s}}}}{2{a}^{2}}{\left(\displaystyle \frac{{da}}{{dt}}\right)}_{i}.\end{eqnarray} \tag{ 29 }$$

However, the contribution of the tidal heating is negligibly small when we consider the tidal evolution of 100 km sized satellites around 1000 km sized dwarf planets (e.g., Arakawa et al. 2019).

3.2.2. Conduction/Convection Cooling

The term for cooling by conduction/convection, Q_con,i , is given by

$$\begin{eqnarray}&&{Q}_{\mathrm{con},i}=\min \left({Q}_{\mathrm{cond},i},{Q}_{\mathrm{conv},i}\right),\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{Q}_{{\rm{cond}},i}=-4\pi {R}_{i}^{2}{\rho }_{i}{c}_{i}{\kappa }_{i}\displaystyle \frac{{T}_{i}-{T}_{{\rm{surf}}}}{{R}_{i}},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{Q}_{{\rm{conv}},i}={{Nu}}_{i}{Q}_{{\rm{cond}},i},\end{eqnarray} \tag{ 32 }$$

where κ_i is the thermal diffusivity, which is given by

$$\begin{eqnarray}&&{\kappa }_{i}=\displaystyle \frac{{k}_{\mathrm{th},i}}{{\rho }_{i}{c}_{i}},\end{eqnarray} \tag{ 33 }$$

and k_th,i is the thermal conductivity. Reese et al. (1999) proposed that the Nusselt number, Nu_i , is given by the following scaling equation (see also Solomatov & Moresi 2000):

$$\begin{eqnarray}&&{{Nu}}_{i}=2.51{{\rm{\Theta }}}_{i}^{-1.2}{{Ra}}_{i}^{0.2},\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{i}=\displaystyle \frac{{E}_{{\rm{a}}}\left({T}_{i}-{T}_{\mathrm{surf}}\right)}{{R}_{\mathrm{gas}}{{T}_{i}}^{2}},\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{{Ra}}_{i}=\displaystyle \frac{{\alpha }_{\exp }{\rho }_{i}{g}_{i}\left({T}_{i}-{T}_{{\rm{surf}}}\right)}{{\kappa }_{i}{\eta }_{i}}{R}_{i}^{3},\end{eqnarray} \tag{ 36 }$$

where Θ_i and Ra_i are the Frank–Kamenetskii parameter and the Rayleigh number, respectively. The surface gravity, g_i , is given by ${g}_{i}={ \mathcal G }{M}_{i}/{{R}_{i}}^{2}$, and E_a is the activation energy, R_gas is the gas constant, α_exp is the thermal expansion coefficient, and η_i is the viscosity. The temperature dependence of the viscosity is given as follows (e.g., Goldsby & Kohlstedt 2001):

$$\begin{eqnarray}&&{\eta }_{i}={\eta }_{\mathrm{ref}}\exp \left[\displaystyle \frac{{E}_{{\rm{a}}}}{{R}_{\mathrm{gas}}{T}_{\mathrm{ref}}}\left(\displaystyle \frac{{T}_{\mathrm{ref}}}{{T}_{i}}-1\right)\right],\end{eqnarray} \tag{ 37 }$$

where η_ref is the reference viscosity and T_ref is the reference temperature (Table 3). We set the surface temperature T_surf = 40 K as assumed in previous studies on the thermal evolution of dwarf planets (e.g., Robuchon & Nimmo 2011; Kamata et al. 2019). The material parameters used in this study are summarized in Table 3.

We stress that the temperature evolution of bodies significantly impacts their spin/orbital evolution. As the viscosity is a function of the temperature, the complex Love number also depends on the temperature (see Appendix B). Then the complex Love number influences the spin/orbital evolution, which is calculated from Equations (10)–(17).

Although it is beyond the scope of this study, ices in undifferentiated icy bodies would not be pure water ice, and the partial melting of ice mixtures may have a great impact on the effective viscosity (e.g., Henning et al. 2009). For the case of ice mixtures containing 1% ammonia with respect to water, the effective viscosity is orders of magnitude lower than that of pure water ice when the temperature is above the ammonia–water eutectic temperature (176 K; e.g., Arakawa & Maeno 1994; Neveu & Rhoden 2017). We could mimic this effect by changing the reference viscosity, which is similar to the impact of crystal grain size on the reference viscosity (see Section 6.4).

3.2.3. Decay Heating

For 100–1000 km sized planetary bodies, the radioactive decay of long-lived isotopes is the principal heat source (e.g., Robuchon & Nimmo 2011). We assume that the elemental abundances of long-lived isotopes of the rocky parts of TNOs are equal to those of carbonaceous chondrites (Lodders 2003; Robuchon & Nimmo 2011). We consider four species as heat sources, namely ²³⁸U, ²³⁵U, ²³²Th, and ⁴⁰K. The half-life, t_HL,j , and the initial heat generation rate per unit mass of rock, H_0,j , of the element j are shown in Table 4.

Table 4. Radioactive Species and Decay Data (Robuchon & Nimmo 2011)

Element j	Half-life tHL,j (Myr)	H0,j (W kg−1 of Rock)
238U	4468	1.88 × 10−12
235U	703.81	3.07 × 10−12
232Th	14,030	1.02 × 10−12
40K	1277	2.15 × 10−11

Download table as:  ASCIITypeset image

Then the decay heating term is given by

$$\begin{eqnarray}&&{Q}_{\mathrm{dec},i}={f}_{\mathrm{rock}}{M}_{i}\sum _{j}{2}^{-t/{t}_{\mathrm{HL},j}}{H}_{0,j},\end{eqnarray} \tag{ 38 }$$

where f_rock is the mass fraction of rock in bodies. We set f_rock = 0.4 in this study. Kiss et al. (2019) determined that the bulk density of Gonggong is 1.75 × 10³ kg m⁻³; our assumption of f_rock = 0.4 seems within an acceptable range if TNOs are made of a mixture of ices, rocks, and organic materials. In future studies, we will investigate the effect of the difference in f_rock on the thermal and orbital evolution.

Figure 3 shows the summary of the final state for standard runs of our simulation. In our standard runs, we set P_G,obs =22.4 hr, P_X,ini = 12 hr, T_X,ini = T_G,ini = T_ini, and η_ref =10¹⁴ Pa s. We change the radius of Xiangliu, R_X, the initial eccentricity, e_ini, and the initial temperature, T_ini. The different markers represent the different final states of the satellite system classified according to the spin state of Xiangliu, semimajor axis, and eccentricity.

Zoom In Zoom Out Reset image size

Figure 3(a) shows the outcomes of the tidal evolution for the case of R_X = 20 km. In this setting, Xiangliu migrates inward and finally collides with Gonggong when the initial eccentricity is e_ini ≤ 0.2 (Type Z). The critical value of e_ini required for the outward migration of Xiangliu is clearly dependent on its radius; smaller satellites have a larger critical value of e_ini. This trend is consistent with our analytic arguments in Section 2.2.

.Figure 3(b) shows the outcomes of the tidal evolution for the case of R_X = 40 km. In this setting, some runs with large e_ini and high T_ini result in extremely eccentric systems with e > 0.9 within 4.5 Gyr (Types Bx/Cx). The parameter space resulting in Types Bx/Cx is wider for larger values of R_X. In addition, the runs with e_ini ≳ 0.7 tend to become extremely eccentric when R_X ≥ 40 km.

Figure 3(f) shows the outcomes of the tidal evolution for the case of R_X = 120 km. With this setting, runs with small e_ini and high T_ini result in Type A, i.e., Xiangliu is trapped in 1:1 spin–orbit resonance at t = 4.5 Gyr. In general, runs with small e_ini and high T_ini tend to turn into Type A after tidal evolution unless their initial eccentricity is too small to satisfy the condition for the outward migration of Xiangliu (see Figures 3(c) and (d)).

We did not find any clear trends for the boundary between Types B and C. In other words, it is difficult to predict whether the secondary is trapped in a higher-order spin–orbit resonance or not trapped in resonances. The fractions of Type B and Type C were 65/336 = 19% and 52/336 = 15%, respectively, according to all our runs with the standard setting (see Table 5). We note that we did not consider the realistic probability distributions of R_X, e_ini, and T_ini for satellite systems around TNOs.

Table 5 shows the statistics of the final state for standard runs of our simulation. Figure 3 shows that the fraction of each type is clearly dependent on the radius of the secondary, R_X. The fractions of Type A, Type Bx, and Type Cx increase with increasing R_X. In contrast, the fractions of Type B, Type C, and Type Z decrease with increasing R_X.

Our results also indicate that other satellites around 1000 km sized TNOs may not be in 1:1 spin–orbit resonance. In the Haumea–Hi'iaka system, the spin period of the secondary, Hi'iaka, is approximately 120 times shorter than its orbital period (Hastings et al. 2016). Although the fast spin of Hi'iaka could be explained by a scenario in which the current satellite system was formed via the catastrophic disruption of the first-generation moon and subsequent reaccretion of fragments (e.g., Schlichting & Sari 2009; Ćuk et al. 2013), we could also explain the fast spin of Hi'iaka without the disruption event. We note, however, that the radius of Hi'iaka is 150 km (Ragozzine & Brown 2009), and the fraction of the nonresonant case (Type C) is small in our standard simulations for the Gonggong–Xiangliu system with R_X ≳ 100 km. Moreover, the spin period of the primary, Haumea, is P_G,obs = 3.9 hr (Rabinowitz et al. 2006), which is an order of magnitude shorter than that of Gonggong. In addition, the nonspherical and highly elongated shape of Haumea (e.g., Ortiz et al. 2017) may also play an important role in their tidal evolution. We will apply our tidal evolution model to the Haumea–Hi'iaka system to unveil its spin–orbit evolution in future studies.

Section 4.1 shows some cases in which the secondary is not in spin–orbit resonance after a 4.5 Gyr orbital evolution (i.e., Type C). Here we discuss the condition for capture into spin–orbit resonance.

Appendix B shows that the imaginary part of the Love number, $\mathrm{Im}[{\tilde{k}}_{2,{\rm{X}}}({\omega }_{{\rm{X}},q})]$, is a minimum/maximum when ${\omega }_{{\rm{X}},q}\simeq \pm {\left({\mu }_{\mathrm{eff},{\rm{X}}}{\tau }_{{\rm{A}},{\rm{X}}}\right)}^{-1}$, where ${\omega }_{{\rm{X}},q}=(2+q)n-2{\dot{\theta }}_{{\rm{X}}}$ is the tidal frequency. When the secondary is in the spin–orbit resonance of order q, the tidal frequency satisfies $-{\left({\mu }_{\mathrm{eff},{\rm{X}}}{\tau }_{{\rm{A}},{\rm{X}}}\right)}^{-1}\lt {\omega }_{{\rm{X}},q}\lt {\left({\mu }_{\mathrm{eff},{\rm{X}}}{\tau }_{{\rm{A}},{\rm{X}}}\right)}^{-1}$, and the spin angular velocity of the secondary, ${\dot{\theta }}_{{\rm{X}}}$, is given by

$$\begin{eqnarray}&&{\dot{\theta }}_{{\rm{X}}}=\displaystyle \frac{2+q}{2}n+\displaystyle \frac{x}{2{\mu }_{\mathrm{eff},{\rm{X}}}{\tau }_{{\rm{A}},{\rm{X}}}},\end{eqnarray} \tag{ 39 }$$

where ∣x∣ < 1 is a dimensionless parameter. For ∣x∣ ≪ 1, the imaginary part of the Love number is given by

$$\begin{eqnarray}&&\mathrm{Im}\left[{\tilde{k}}_{2,{\rm{X}}}\left((2+q)n-2{\dot{\theta }}_{{\rm{X}}}\right)\right]\simeq \displaystyle \frac{3}{2}x.\end{eqnarray} \tag{ 40 }$$

Here, we assume that the following relations are satisfied for integers $q^{\prime} \ne q$ (see Section 3.1 for the definitions of ${{ \mathcal A }}_{{\rm{X}},q}$ and ${{ \mathcal B }}_{{\rm{X}},q}$):

$$\begin{eqnarray}&&{{ \mathcal A }}_{{\rm{X}},q}\gg {{ \mathcal A }}_{{\rm{X}},q^{\prime} },{{ \mathcal B }}_{{\rm{X}},q^{\prime} },{{ \mathcal B }}_{{\rm{X}},q}.\end{eqnarray} \tag{ 41 }$$

Under this assumption, $d{\dot{\theta }}_{{\rm{X}}}/{dt}$, (da/dt)_X, and $d({\dot{\theta }}_{{\rm{X}}}/n)/{dt}$ near the spin–orbit resonance are approximately given by

$$\begin{eqnarray}\displaystyle \begin{array}{ccc}\frac{1}{n}\frac{d{\dot{\theta }}_{{\rm{X}}}}{{dt}} & \simeq & -\frac{9}{4}\frac{{M}_{{\rm{G}}}{M}_{{\rm{X}}}}{{M}_{\mathrm{tot}}}\frac{{M}_{{\rm{G}}}}{{M}_{{\rm{X}}}}\frac{{R}_{{\rm{X}}}^{2}}{{I}_{{\rm{X}}}}{\left(\frac{{R}_{{\rm{X}}}}{a}\right)}^{3}{\left[{G}_{2,0,q}(e)\right]}^{2}xn,\\ & \simeq & \frac{45}{8}\frac{{M}_{{\rm{G}}}}{{M}_{{\rm{X}}}}{\left(\frac{{R}_{{\rm{X}}}}{a}\right)}^{3}{\left[{G}_{2,0,q}(e)\right]}^{2}{xn},\end{array}\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{a}{\left(\displaystyle \frac{{da}}{{dt}}\right)}_{{\rm{X}}}\simeq \displaystyle \frac{9(2+q)}{4}\displaystyle \frac{{M}_{{\rm{G}}}}{{M}_{{\rm{X}}}}{\left(\displaystyle \frac{{R}_{{\rm{X}}}}{a}\right)}^{5}{\left[{G}_{\mathrm{2,0},q}(e)\right]}^{2}{xn},\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d}{{dt}}\left(\displaystyle \frac{{\dot{\theta }}_{{\rm{X}}}}{n}\right) & \simeq & -\displaystyle \frac{45}{8}\left[{\left(\displaystyle \frac{a}{{R}_{{\rm{X}}}}\right)}^{2}-\displaystyle \frac{6}{5}{\left(\displaystyle \frac{2+q}{2}\right)}^{2}\right]\\ & & \times \,\displaystyle \frac{{M}_{{\rm{G}}}}{{M}_{{\rm{X}}}}{\left(\displaystyle \frac{{R}_{{\rm{X}}}}{a}\right)}^{5}{\left[{G}_{\mathrm{2,0},q}(e)\right]}^{2}{xn}+\displaystyle \frac{3(2+q)}{4}\displaystyle \frac{1}{a}{\left(\displaystyle \frac{{da}}{{dt}}\right)}_{{\rm{G}}}.\end{array}\end{eqnarray} \tag{ 44 }$$

Assuming that (da/dt)_G > 0, the required conditions for a stable equilibrium around ${\dot{\theta }}_{{\rm{X}}}\simeq (2+q)n/2\,+\,x/({\mu }_{\mathrm{eff},{\rm{X}}}{\tau }_{{\rm{A}},{\rm{X}}})$ with ∣x∣ ≪ 1 are

$$\begin{eqnarray}&&{\left(\displaystyle \frac{a}{{R}_{{\rm{X}}}}\right)}^{2}\gg \displaystyle \frac{6}{5}{\left(\displaystyle \frac{2+q}{2}\right)}^{2},\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&{\left[{G}_{\mathrm{2,0},q}(e)\right]}^{2}\gg \displaystyle \frac{2(2+q)}{15n}{\left(\displaystyle \frac{a}{{R}_{{\rm{G}}}}\right)}^{3}\displaystyle \frac{1}{a}{\left(\displaystyle \frac{{da}}{{dt}}\right)}_{{\rm{G}}}.\end{eqnarray} \tag{ 46 }$$

On the other hand, for the case of ${\left({da}/{dt}\right)}_{{\rm{G}}}\lt 0$, the required conditions for a stable equilibrium are

$$\begin{eqnarray}&&{\left(\displaystyle \frac{a}{{R}_{{\rm{X}}}}\right)}^{2}\ll \displaystyle \frac{6}{5}{\left(\displaystyle \frac{2+q}{2}\right)}^{2},\end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray}&&{\left[{G}_{\mathrm{2,0},q}(e)\right]}^{2}\gg -\displaystyle \frac{4}{9(2+q)n}\displaystyle \frac{{M}_{{\rm{X}}}}{{M}_{{\rm{G}}}}{\left(\displaystyle \frac{a}{{R}_{{\rm{X}}}}\right)}^{5}\displaystyle \frac{1}{a}{\left(\displaystyle \frac{{da}}{{dt}}\right)}_{{\rm{G}}}.\end{eqnarray} \tag{ 48 }$$

Except for the case of q = 0, ${\left[{G}_{\mathrm{2,0},q}(e)\right]}^{2}\to 0$ when e → 0 (see Appendix A). Thus, a finite eccentricity is needed for capture into higher-order spin–orbit resonances with q ≠ 0.

Figure 4 shows the distributions of the final eccentricity and semimajor axis, e_fin and a_fin. The observed values for the Gonggong–Xiangliu system are e_fin = 0.3 and a_fin/R_G =24,000 km/600 km = 40; the green stars indicate the observed values.

Zoom In Zoom Out Reset image size

We found that the final eccentricity is e_fin ≲ 0.05 for Type A. Xiangliu has a moderate eccentricity of e_obs = 0.3 (Kiss et al. 2019), and therefore it would not be in 1:1 spin–orbit resonance. The final eccentricity of satellite systems classified as Type B and Type C is in the wide range 0 < e_fin < 0.9. We note, however, that the fraction of systems whose e_fin is comparable to e_obs depends strongly on the radius of Xiangliu. The fraction of systems whose final eccentricity is in the range 0.2 ≤ e_fin ≤ 0.5 is as follows: 10/56 = 18% for the case of R_X = 20 km, 5/56 = 8.9% for 40 km, 4/56 = 7.1% for 60 km, 3/56 = 5.4% for 80 km, 4/56 = 7.1% for 100 km, and 0/56 = 0% for 120 km. Thus, the radius of Xiangliu is likely not larger than 100 km.

We also found a clear relation between a_fin and e_fin. The dashed lines in Figure 4 are given by

$$\begin{eqnarray}&&{p}_{{\rm{f}}{\rm{i}}{\rm{n}}}={a}_{{\rm{f}}{\rm{i}}{\rm{n}}}\left(1-{e}_{{\rm{f}}{\rm{i}}{\rm{n}}}^{2}\right)={\rm{c}}{\rm{o}}{\rm{n}}{\rm{s}}{\rm{t}}\end{eqnarray} \tag{ 49 }$$

for each R_s. We explain the reason why a_fin and e_fin are given by Equation (49) in Section 5.2. Figure 4 shows that the radius of Xiangliu might be close to 100 km from the point of view of p_fin. We note that we should consider the dependence of p_fin on many properties of the system, such as the viscoelastic properties of icy bodies and their thermal histories (see Sections 6.3 and 6.4).

Our results suggest that not only Gonggong–Xiangliu but also the Quaoar–Weywot system would not be in 1:1 spin–orbit resonance. Fraser et al. (2013) reported that the eccentricity of the Quaoar–Weywot system is e = 0.13–0.16. This eccentricity indicates that the system might be classified as Type B or Type C, because systems classified as Type A have a small eccentricity of e ≲ 0.05 in general.

Figure 4 shows that the final semilatus rectum, p_fin, hardly depends on the final eccentricity, e_fin. In this section, we derive an analytic solution for the evolution of the semilatus rectum.

The time derivative of p_orb is given by the following equation:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{1}{{p}_{\mathrm{orb}}}\displaystyle \frac{{{dp}}_{\mathrm{orb}}}{{dt}} & = & \displaystyle \frac{1}{a}\displaystyle \frac{{da}}{{dt}}-\displaystyle \frac{2{e}^{2}}{1-{e}^{2}}\displaystyle \frac{1}{e}\displaystyle \frac{{de}}{{dt}},\\ & = & \ 3n\displaystyle \frac{{M}_{{\rm{X}}}}{{M}_{{\rm{G}}}}{\left(\displaystyle \frac{{R}_{{\rm{G}}}}{a}\right)}^{5}{\left(1-{e}^{2}\right)}^{-1/2}\displaystyle \sum _{q=-\infty }^{+\infty }{{ \mathcal A }}_{{\rm{G}},q}\\ & & +\ 3n\displaystyle \frac{{M}_{{\rm{G}}}}{{M}_{{\rm{X}}}}{\left(\displaystyle \frac{{R}_{{\rm{X}}}}{a}\right)}^{5}{\left(1-{e}^{2}\right)}^{-1/2}\displaystyle \sum _{q=-\infty }^{+\infty }{{ \mathcal A }}_{{\rm{X}},q}.\end{array}\end{eqnarray} \tag{ 50 }$$

Here, we assume that the contribution of the tidal torques caused by the secondary is negligible:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{1}{{p}_{\mathrm{orb}}}\displaystyle \frac{{{dp}}_{\mathrm{orb}}}{{dt}} & \simeq & 3n\displaystyle \frac{{M}_{{\rm{X}}}}{{M}_{{\rm{G}}}}{\left(\displaystyle \frac{{R}_{{\rm{G}}}}{a}\right)}^{5}\\ & & \times \,{\left(1-{e}^{2}\right)}^{-1/2}\displaystyle \sum _{q=-\infty }^{+\infty }{{ \mathcal A }}_{{\rm{G}},q}.\end{array}\end{eqnarray} \tag{ 51 }$$

Assuming that the spin period of the primary is sufficiently shorter than the orbital period (i.e., ${\dot{\theta }}_{{\rm{G}}}\gg n$), we can rewrite the term ${\sum }_{q}{{ \mathcal A }}_{{\rm{G}},q}$ as follows:

$$\begin{eqnarray}&&\sum _{q=-\infty }^{+\infty }{{ \mathcal A }}_{{\rm{G}},q}\simeq \mathrm{Im}\left[{\tilde{k}}_{2,{\rm{G}}}\left(-2{\dot{\theta }}_{{\rm{G}}}\right)\right]\sum _{q=-\infty }^{+\infty }{\left[{G}_{\mathrm{2,0},q}(e)\right]}^{2}.\end{eqnarray} \tag{ 52 }$$

We found that ${\sum }_{q}{\left[{G}_{\mathrm{2,0},q}(e)\right]}^{2}$ is approximately given by

$$\begin{eqnarray}&&\sum _{q=-\infty }^{+\infty }{\left[{G}_{\mathrm{2,0},q}(e)\right]}^{2}\simeq {\left(1-{e}^{2}\right)}^{-6},\end{eqnarray} \tag{ 53 }$$

in the range 0 ≤ e ≤ 0.8 (see Figure A2). Therefore, the time evolution of p_orb is approximately given by the following equation:

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}\displaystyle \frac{{{dp}}_{{\rm{orb}}}}{{dt}} & \simeq & 3\sqrt{{ \mathcal G }{M}_{{\rm{tot}}}}\displaystyle \frac{{M}_{{\rm{X}}}}{{M}_{{\rm{G}}}}{R}_{{\rm{G}}}^{5}{p}_{{\rm{orb}}}^{-11/2}\\ & & \cdot {\rm{Im}}\left[{\tilde{k}}_{2,{\rm{G}}}\left(-2{\dot{\theta }}_{{\rm{G}}}\right)\right],\end{array}\end{array}\end{eqnarray} \tag{ 54 }$$

and we obtain p_orb = p_orb(t) as

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{p}_{{\rm{orb}}}(t) & \simeq & \left[\displaystyle \frac{39}{2}\sqrt{{ \mathcal G }{M}_{{\rm{tot}}}}{R}_{{\rm{G}}}^{2}\right.\\ & & \cdot {\left.{\displaystyle \int }_{0}^{t}{{dt}}^{{\rm{{\prime} }}}\,{\rm{Im}}\left[{\tilde{k}}_{2,{\rm{G}}}\left(-2{\dot{\theta }}_{{\rm{G}}}\right)\right]\right]}^{2/13}{R}_{{\rm{X}}}^{6/13},\end{array}\end{array}\end{eqnarray} \tag{ 55 }$$

which does not include the eccentricity, e, in the equation.

Figure 5 shows the final semilatus rectum, p_fin, as a function of the radius of the secondary, R_X. We only used the data classified as Type A/B/C in the figure. We plotted the average value and twice the standard error of p_fin for each R_X. Equation (55) shows that the final semilatus rectum is proportional to ${R}_{{\rm{X}}}^{6/13};$ our numerical results are consistent with the theoretical prediction.

Zoom In Zoom Out Reset image size

Figure 6 shows the distribution of the final eccentricity and spin period of the secondary, e_fin and P_X,fin, for standard runs of our simulation. The spin period of Xiangliu at t = 4.5 Gyr is in the range 10 hr ≲ P_X,fin ≲ 10³ hr, and the range of P_X,fin depends on R_X. In the case of R_X = 20 km, the final spin period of Xiangliu is in the range 10 hr ≲ P_X,fin ≲ 10² hr and is always shorter than the orbital period, P_orb,fin. In contrast, for R_X = 100 and 120 km, the final spin period of Xiangliu is in the range 10² hr ≲ P_X,fin ≲ 10³ hr, and P_X,fin is significantly longer than the initial spin period, P_X,ini = 12 hr. Thus, the determination of the spin period of Xiangliu by future observations and light-curve analyses is the key to estimating its radius.

Zoom In Zoom Out Reset image size

The final spin period of Xiangliu is close to the initial spin period, P_X,fin ≃ P_X,ini, for the case of Type C. The spin period of satellites formed via a giant impact should provide abundant information regarding the condition of the colliding TNOs, including the mass ratio, the spin state before the impact, and the (un)differentiated state. Therefore, we want to determine the spin period of satellites from astronomical observations for Gonggong–Xiangliu and for other satellite systems around large TNOs.

Figure 7 shows the final eccentricity of the system, e_fin, for standard runs of our simulation. We note that we do not plot the results classified into Type Z.

Zoom In Zoom Out Reset image size

Figure 7 shows that satellite systems with an initial eccentricity of 0.2 ≤ e_ini ≤ 0.5 tend to a final eccentricity of 0.2 ≤ e_fin ≤ 0.5. When the initial eccentricity is e_ini ≤ 0.1, the final eccentricity at t = 4.5 Gyr is typically e_fin ≤ 0.1 or the system becomes Type Z. When the initial eccentricity is e_ini ≥ 0.6, on the other hand, the final eccentricity at t = 4.5 Gyr is typically e_fin ≥ 0.5 for the case of R_X ≤ 60 km, and e_fin ≥ 0.9 or e_fin ≤ 0.1 for the case of R_X ≥ 80 km.

Therefore, we conclude that the initial eccentricity of the Gonggong–Xiangliu system after the moon-forming giant impact may be in the range 0.2 ≲ e_ini ≲ 0.5. This nonzero e_ini is also consistent with the fact that e_ini ≳ 0.2 is needed to avoid the inward migration of Xiangliu (see Figure 3).

Debris disks can be formed by a giant impact (Canup 2005), collisional disruption of a pre-existing satellite (Hyodo & Charnoz 2017), or tidal disruption of a passing body (Hyodo et al. 2017). If a single large satellite is formed via accretion of debris disk materials around Gonggong, the initial eccentricity becomes small: e ≲ 0.1 (e.g., Ida et al. 1997; Kokubo et al. 2000). In this case, we cannot reproduce the observed eccentricity of the system. Thus, the Gonggong–Xiangliu system was more likely to be born as an "intact moon" rather than "disk-origin moon," and their formation process is similar to that of the Pluto–Charon system (e.g., Canup 2005; Sekine et al. 2017; Arakawa et al. 2019).

In our standard models, we set the current spin period of Gonggong as P_G,obs = 22.4 hr. The spin period of Gonggong was obtained from the light-curve analysis by Pál et al. (2016) and has two possible values: P_G,obs = 22.4 hr or 44.8 hr. In Section 2.2, we discuss the impact of P_G,obs on the location of the initial corotation radius. We showed that a large e_ini is necessary to avoid the inward migration of Xiangliu when P_G,obs = 44.8 hr; however, a nonzero value of e_ini is required even for the case of P_G,obs = 22.4 hr. As we cannot exclude the possibility of P_G,obs = 44.8 hr from the light-curve analyses, we discuss the case in this section.

Panels (a) in Figures 8–11 show the results for the case of P_G,obs = 44.8 hr. Figure 8(a) shows that the satellite system results in Type Z for an initial eccentricity of e_ini ≤ 0.4. The critical value of e_ini to avoid the inward migration of Xiangliu is significantly larger than that for our standard runs (see Figure 3(c)), which is consistent with our prediction from the analytical estimate of the corotation radius in Section 2.2.

Figure 9(a) shows the distribution of e_fin and a_fin. The final eccentricity of satellite systems classified as Type B and Type C is in the wide range 0 < e_fin < 0.9, as in the case shown in Figure 4(c). The final semilatus rectum was slightly lower than the observed value for the Gonggong–Xiangliu system. Figure 10(a) shows the distribution of e_fin and P_X,fin. The general trend is similar to that of the standard shown in Figure 6(c).

Figure 11(a) is the color map for e_fin. In the case of P_G,obs = 44.8 hr, satellite systems with an initial eccentricity of e_ini ≃ 0.6 tend to a final eccentricity of 0.2 ≤ e_fin ≤ 0.5. Systems with large e_ini tend to result in a large e_fin; this trend is similar to that observed in the standard runs (see Figure 7(c)); although the suitable range of e_ini to reproduce the observed value for the Gonggong–Xiangliu system is different.

In our standard models, we set the initial spin period of Xiangliu as P_X,ini = 12 hr. This assumption is based on numerical simulations of moon-forming giant impacts of 1000 km sized TNOs (see Figure 1(b)). In this case, P_X,ini ≪ P_orb,ini and a large number of runs are trapped in higher-order spin–orbit resonances (see Figure 3). In this section, we instead set P_X,ini = P_orb,ini. Although the assumption of P_X,ini = P_orb,ini seems unrealistic, we discuss the effects of the initially synchronized Xiangliu on the tidal evolution of the system to better understand our results.

Panels (b) in Figures 8–11 show the results for the case of P_X,ini = P_orb,ini. Figure 8(b) shows that the final state of the system is not Type A (i.e., 1:1 spin–orbit resonance) in some cases even if the initial spin period of the secondary is P_X,ini = P_orb,ini.

Figure 9(b) shows the distribution of e_fin and a_fin. For the case of P_X,ini = P_orb,ini, the final eccentricity is in the wide range, even for Type A (0 < e_fin < 0.7). We note, however, that the assumption of P_X,ini = P_orb,ini is unrealistic if they are formed via giant impacts (see Figure 1(b)). Figure 10(b) shows the distribution of e_fin and P_X,fin. In this case, the final spin period of Xiangliu lies in the range 10² hr ≲ P_X,fin ≲ 10³ hr. In addition, all satellite systems with e_fin ≤ 0.6 are classified as Type A.

Figure 11(b) is the color map for e_fin. In the case of P_X,ini = P_orb,ini, satellite systems with an initial eccentricity of 0.2 ≤ e_ini ≤ 0.4 tend to a final eccentricity of 0.2 ≤ e_fin ≤ 0.5. This result is similar to that of our standard runs (see Figure 7).

In our standard models, we assumed that the initial temperatures of Gonggong and Xiangliu are the same: T_G,ini = T_X,ini = T_ini. Although this assumption is somewhat natural, we can imagine a case in which the initial temperature of Xiangliu is lower than that of Gonggong. The motivation for this setting is as follows.

Giant impact simulations of 1000 km sized TNOs (e.g., Canup 2005; Sekine et al. 2017; Arakawa et al. 2019) revealed that satellites are directly formed as intact fragments of the impactor ("intact moons"). As the mass of the impactor is several times smaller than that of the target body, the cooling rate of the impactor would be faster than that of the target, and the temperature of the impactor might be lower than that of the target upon collision. Although we need to assess the effect of the impact heating on the initial temperatures of the primary and secondary, the setting in which the initial temperature of Xiangliu is lower than that of Gonggong seems to be reasonable.

Panels (c) in Figures 8–11 show the results for the case of T_X,ini = 120 K. Figure 8(c) shows that no systems classified as Type A are formed in this case.

Figure 9(c) shows the distribution of e_fin and a_fin, and Figure 10(c) shows the distribution of e_fin and P_X,fin. These distributions are generally similar to those for our standard runs (Figures 4(c) and 6(c)); however, no systems result in 1:1 spin–orbit resonance.

Figure 11(c) is the color map for e_fin. In the case of T_X,ini = 120 K, satellite systems with an initial eccentricity of 0.2 ≤ e_ini ≤ 0.4 tend to a final eccentricity of 0.2 ≤ e_fin ≤ 0.5. This result is also similar to those of our standard runs (see Figure 7(c)).

In our standard models, we set the reference viscosity of icy bodies as η_ref = 10¹⁴ Pa s. This value is widely used in studies of differentiated icy bodies, including of the thermal evolution of Pluto (e.g., Kamata et al. 2019), the orbital evolution of Pluto–Charon (e.g., Renaud et al. 2021), and the tidal heating of Europa and Titan (e.g., Tobie et al. 2005).

However, we note that the reference viscosity of icy bodies depends strongly on the grain size of ice crystals (e.g., Goldsby & Kohlstedt 2001; Kubo et al. 2006), and η_ref for ice mixed with a small volume fraction of dust grains may be orders of magnitude lower than the canonical value for differentiated ice bodies (e.g., Kubo et al. 2009). When the ice viscosity is controlled by diffusion creep, the reference viscosity, η_ref, at the reference temperature, T_ref = 273 K, is given by (Kalousová et al. 2016)

$$\begin{eqnarray}&&{\eta }_{\mathrm{ref}}=\displaystyle \frac{{{d}_{\mathrm{ice}}}^{2}}{A}\exp \left(\displaystyle \frac{{E}_{{\rm{a}}}}{{R}_{\mathrm{gas}}{T}_{\mathrm{ref}}}\right),\end{eqnarray} \tag{ 56 }$$

where d_ice is the grain size of the ice crystals, A = 3.3 ×10⁻¹⁰ Pa⁻¹ m² s⁻¹ is the pre-exponential constant, R_gas =8.31 J K⁻¹ mol⁻¹ is the gas constant, and E_a is the activation energy. Assuming E_a = 60 kJ mol⁻¹ (Kamata et al. 2019), we obtain the following relation between d_ice and η_ref:

$$\begin{eqnarray}&&{\eta }_{\mathrm{ref}}=9.3\times {10}^{8}\ {\left(\displaystyle \frac{{d}_{\mathrm{ice}}}{1\,\mu {\rm{m}}}\right)}^{2}\ \mathrm{Pa}\ {\rm{s}}.\end{eqnarray} \tag{ 57 }$$

When undifferentiated icy bodies are mixtures of micron-sized dust grains and matrix ices, the grain size of ice crystals would be maintained in the micron size owing to the Zener pinning effect (Kubo et al. 2009). Thus, the reference viscosity of η_ref ∼ 10¹⁰ Pa s may be somewhat reasonable for undifferentiated icy bodies.

Panels (d) in Figures 8–11 show the results for the case of η_ref = 10¹⁰ Pa s. Figure 8(d) shows that runs with a small e_ini and high T_ini tend to result in Type A.

Figure 9(d) shows the distribution of e_fin and a_fin, and Figure 10(d) shows the distribution of e_fin and P_X,fin. These distributions are similar to that of our standard runs. The final semilatus rectum, p_fin, for the case of η_ref = 10¹⁰ Pa s is not significantly different from that for the standard case of η_ref = 10¹⁴ Pa s. This is because the temperature of Gonggong, T_G, also depends on η_ref; and the viscosity of Gonggong, η_G = η(T_G), is not strongly dependent on η_ref between the two settings, η_ref = 10¹⁰ and 10¹⁴ Pa s. Figure D8 shows the thermal and orbital evolutions of Gonggong and Xiangliu for the case of η_ref = 10¹⁰ Pa s, T_ini = 140 K, and e_ini = 0.4.

Figure 11(d) is the color map for e_fin. In the case of η_ref = 10¹⁰ Pa s, satellite systems with an initial temperature of T_ini ≥ 160 K tend to result in Type A or Type Bx/Cx (i.e., a final eccentricity of e_fin ≤ 0.1 or e_fin ≥ 0.9). Therefore, the initial temperature of T_ini ≤ 140 K might be suitable for explaining the formation of moderately eccentric satellite systems when η_ref = 10¹⁰ Pa s.

Viscoelastic models are widely used to describe the tidal response of small icy bodies. The Andrade model is an empirical model based on experimental results for the viscous flow of metals and ice (e.g., Andrade 1910; Glen 1955). The model is widely used to compute tidal dissipation within icy satellites. For example, Shoji et al. (2013) calculated the tidal heating of Enceladus using the Andrade model. Neveu & Rhoden (2019) also performed numerical simulations of coupling thermal, geophysical, and orbital evolution using the model.

The Andrade viscoelastic model contains an anelastic part, and the effect of elasticity decreases in the high-frequency region. The complex shear modulus, $\tilde{\mu }(\omega )$, is expressed in terms of a creep function as follows:

$$\begin{eqnarray}&&\tilde{\mu }(\omega )=\displaystyle \frac{1}{\tilde{J}(\omega )},\end{eqnarray} \tag{ B1 }$$

where $\tilde{J}(\omega )$ is the creep function and ω is the angular velocity of the forcing cycle. The creep function is given by (e.g., Efroimsky 2012a, 2012b)

$$\begin{eqnarray}\begin{array}{rcl}\tilde{J}(\omega ) & = & \left[\displaystyle \frac{1}{\mu }+\displaystyle \frac{1}{\mu {\left({\tau }_{{\rm{A}}}\omega \right)}^{\alpha }}\cos \left(\displaystyle \frac{\alpha \pi }{2}\right){\rm{\Gamma }}\left(\alpha +1\right)\right]\\ & & -\,i\left[\displaystyle \frac{1}{\eta \omega }+\displaystyle \frac{1}{\mu {\left({\tau }_{{\rm{A}}}\omega \right)}^{\alpha }}\sin \left(\displaystyle \frac{\alpha \pi }{2}\right){\rm{\Gamma }}\left(\alpha +1\right)\right],\end{array}\end{eqnarray} \tag{ B2 }$$

where μ is the shear modulus, η is the viscosity, and 0 < α < 1 is the Andrade exponent. We assume α = 0.33 for ice (e.g., Rambaux et al. 2010). The relaxation time of the Andrade model is given by τ_A = η/μ.

We note that $\tilde{\mu }$ depends on a dimensionless parameter, τ_A ω, and can be divided into two regions. For the case of τ_A ω ≫ 1, the complex shear modulus is approximately given by

$$\begin{eqnarray}&&\tilde{\mu }(\omega )\simeq \mu \left[1+i{\left({\tau }_{{\rm{A}}}\omega \right)}^{-\alpha }\sin \left(\displaystyle \frac{\alpha \pi }{2}\right){\rm{\Gamma }}\left(\alpha +1\right)\right]\end{eqnarray} \tag{ B3 }$$

and $| \tilde{\mu }(\omega )| /\mu \simeq 1$, respectively. On the other hand, when τ_A ω ≪ 1, the complex shear modulus is approximately given by

$$\begin{eqnarray}&&\tilde{\mu }(\omega )\simeq \eta \omega \left[{\left({\tau }_{{\rm{A}}}\omega \right)}^{1-\alpha }\cos \left(\displaystyle \frac{\alpha \pi }{2}\right){\rm{\Gamma }}\left(\alpha +1\right)+i\right]\end{eqnarray} \tag{ B4 }$$

and $| \tilde{\mu }(\omega )| /\mu \simeq {\tau }_{{\rm{A}}}\omega$, respectively.

For a homogeneous spherical body, the tidal potential Love number, ${\tilde{k}}_{2}(\omega )$, is given by

$$\begin{eqnarray}&&{\tilde{k}}_{2}(\omega )=\displaystyle \frac{3\rho {gR}}{2\rho {gR}+19\tilde{\mu }(\omega )}=\displaystyle \frac{3/2}{1+{\mu }_{\mathrm{eff}}\left(\tilde{\mu }/\mu \right)},\end{eqnarray} \tag{ B5 }$$

where $g={ \mathcal G }M/{R}^{2}$ is the surface gravity, ${ \mathcal G }$ is the gravitational constant, and ρ and R are the density and radius of the body, respectively. We introduce the (dimensionless) effective rigidity, μ_eff, as follows (e.g., Goldreich & Sari 2009; Cheng et al. 2014):

$$\begin{eqnarray}&&{\mu }_{\mathrm{eff}}\equiv \displaystyle \frac{19\mu }{2\rho {gR}},\end{eqnarray} \tag{ B6 }$$

and we found that μ_eff ∼ 10² for 1000 km sized TNOs. The imaginary part of the complex Love number is negative when ω is positive. When the angular velocity of the forcing cycle is negative, the complex Love number is given by

$$\begin{eqnarray}&&{\tilde{k}}_{2}(-\omega )=\overline{{\tilde{k}}_{2}(\omega )},\end{eqnarray} \tag{ B7 }$$

where $\overline{{\tilde{k}}_{2}(\omega )}$ is the complex conjugate of ${\tilde{k}}_{2}(\omega )$.

In the case of μ_eff ≫ 1, the tidal potential Love number is classified into three regions: (i) high-frequency region, (ii) intermediate region, and (iii) low-frequency regions (see also Efroimsky 2012b). We found that $\left|\mathrm{Im}\left[{\tilde{k}}_{2}(\omega )\right]\right|$ reaches a maximum at ${\tau }_{{\rm{A}}}\omega \simeq \pm {\mu }_{{\rm{eff}}}^{-1};$ the maximum value is $\left|\mathrm{Im}\left[{\tilde{k}}_{2}(\omega )\right]\right|\simeq 3/4$.

B.2.1. High-frequency Region

When τ_A ω ≫ 1, the absolute value and the imaginary part of the complex Love number, $\left|{\tilde{k}}_{2}(\omega )\right|$ and $\mathrm{Im}\left[{\tilde{k}}_{2}(\omega )\right]$, are approximately given by

$$\begin{eqnarray}&&\left|{\tilde{k}}_{2}(\omega )\right|\simeq \displaystyle \frac{3}{2{\mu }_{\mathrm{eff}}},\end{eqnarray} \tag{ B8 }$$

$$\begin{eqnarray}&&\mathrm{Im}\left[{\tilde{k}}_{2}(\omega )\right]\simeq -\displaystyle \frac{3}{2{\mu }_{\mathrm{eff}}}{\left({\tau }_{{\rm{A}}}\omega \right)}^{-\alpha }\sin \left(\displaystyle \frac{\alpha \pi }{2}\right){\rm{\Gamma }}\left(\alpha +1\right).\end{eqnarray} \tag{ B9 }$$

The tidal quality factor, ${ \mathcal Q }(\omega )\equiv \left|\mathrm{Im}\left[{\tilde{k}}_{2}(\omega )\right]/{\tilde{k}}_{2}(\omega )\right|$, is given by

$$\begin{eqnarray}&&{ \mathcal Q }(\omega )\simeq {\left({\tau }_{{\rm{A}}}\omega \right)}^{-\alpha }\sin \left(\displaystyle \frac{\alpha \pi }{2}\right){\rm{\Gamma }}\left(\alpha +1\right).\end{eqnarray} \tag{ B10 }$$

B.2.2. Intermediate Region

When ${{\mu }_{\mathrm{eff}}}^{-1}\ll {\tau }_{{\rm{A}}}\omega \ll 1$, $\left|{\tilde{k}}_{2}(\omega )\right|$, $\mathrm{Im}\left[{\tilde{k}}_{2}(\omega )\right]$, and ${ \mathcal Q }(\omega )$ are approximately given by

$$\begin{eqnarray}&&\left|{\tilde{k}}_{2}(\omega )\right|\simeq \displaystyle \frac{3}{2}{\left({\mu }_{\mathrm{eff}}{\tau }_{{\rm{A}}}\omega \right)}^{-1},\end{eqnarray} \tag{ B11 }$$

$$\begin{eqnarray}&&\mathrm{Im}\left[{\tilde{k}}_{2}(\omega )\right]\simeq -\displaystyle \frac{3}{2}{\left({\mu }_{\mathrm{eff}}{\tau }_{{\rm{A}}}\omega \right)}^{-1},\end{eqnarray} \tag{ B12 }$$

and

$$\begin{eqnarray}&&{ \mathcal Q }(\omega )\simeq 1.\end{eqnarray} \tag{ B13 }$$

B.2.3. Low-frequency Region

When ${\tau }_{{\rm{A}}}\omega \ll {\mu }_{{\rm{eff}}}^{-1}$, $\left|{\tilde{k}}_{2}(\omega )\right|$, $\mathrm{Im}\left[{\tilde{k}}_{2}(\omega )\right]$, and ${ \mathcal Q }(\omega )$ are approximately given by

$$\begin{eqnarray}&&\left|{\tilde{k}}_{2}(\omega )\right|\simeq \displaystyle \frac{3}{2},\end{eqnarray} \tag{ B14 }$$

$$\begin{eqnarray}&&\mathrm{Im}\left[{\tilde{k}}_{2}(\omega )\right]\simeq -\displaystyle \frac{3}{2}{\mu }_{\mathrm{eff}}{\tau }_{{\rm{A}}}\omega ,\end{eqnarray} \tag{ B15 }$$

and

$$\begin{eqnarray}&&{ \mathcal Q }(\omega )\simeq {\mu }_{\mathrm{eff}}{\tau }_{{\rm{A}}}\omega .\end{eqnarray} \tag{ B16 }$$

Figure D1 shows a typical tidal evolution pathway of a satellite system resulting in Type A. The initial condition of this model is R_X = 60 km, T_ini = 200 K, and e_ini = 0.4 (see Figure 3(c)). The final semimajor axis and eccentricity are a_fin/R_G = 26.0 and e_fin = 5.6 × 10⁻³, respectively; and the spin of the secondary is in 1:1 spin–orbit resonance at t = 4.5 Gyr.

The thermal evolution of both the primary and secondary hardly depends on their spin/orbital evolution. This is because their thermal evolution pathways are determined by the balance between decay heating and conduction/convection cooling.

Figure D2 shows the first 10 Myr of the tidal evolution of this case. Both the primary and secondary are initially captured into spin–orbit resonances. At t = 2.872 Myr, the eccentricity becomes e = 4.15 × 10⁻², and the secondary is released from 5:2 spin–orbit resonance. Then, the primary is also released from 2:1 spin–orbit resonance at t = 2.875 Myr. The secondary is recaptured into 2:1 spin–orbit resonance at t = 2.897 Myr, and finally trapped in 1:1 resonance at t = 3.54 Myr.

Figure D3 shows a typical tidal evolution pathway of a satellite system resulting in Type B. The initial condition of this model is R_X = 60 km, T_ini = 140 K, and e_ini = 0.4 (see Figure 3(c)). The final semimajor axis and eccentricity are a_fin/R_G = 29.7 and e_fin = 0.37, respectively, and the spin of the secondary is in 5:2 spin–orbit resonance at t = 4.5 Gyr.

The temperature of the 1000 km sized primary hardly depends on the initial temperature approximately 1 Gyr after the start of calculations (see Figures D1(d), D3(d), and D5(d)). This shows that the temperature of the primary reaches the equilibrium determined by the balance of decay heating and conduction/convection cooling at t ≲ 1 Gyr.

Figure D4 shows the first 300 Myr of the tidal evolution of this case. Both the primary and secondary are initially captured into spin–orbit resonances. The primary is released from 2:1 spin–orbit resonance at t = 113 Myr, at which time the eccentricity begins increasing. Conversely, 5:2 spin–orbit resonance of the secondary is stable in this case.

Figure D5 shows a typical tidal evolution pathway of a satellite system resulting in Type C. The initial condition of this model is R_X = 60 km, T_ini = 120 K, and e_ini = 0.4 (see Figure 3(c)). The final semimajor axis and eccentricity are a_fin/R_G = 28.0 and e_fin = 0.32, respectively, and the spin of the secondary is not in spin–orbit resonances at t = 4.5 Gyr.

Figure D6 shows the first 300 Myr of the tidal evolution of this case. The secondary is captured into 3:1 spin–orbit resonance at first, and then the primary is captured into 2:1 spin–orbit resonance. Both the primary and secondary are released from the spin–orbit resonance at t = 193.8 Myr, and the breaking of the spin–orbit resonance of the primary is prior to that of the secondary. The eccentricity starts to increase at this time, similar to the case shown in Figure D4. In Section 4.2, we show the conditions required to maintain spin–orbit resonances.

Figure D7 shows a typical tidal evolution pathway of a satellite system resulting in Type Z. The initial condition of this model is R_X = 60 km, T_ini = 200 K, and e_ini = 0.1 (see Figure 3(c)). The final semimajor axis and eccentricity are a_fin/R_G = 1.1 and e_fin = 1.9 × 10⁻⁵, and the spin of the secondary is in 1:1 spin–orbit resonance at t = 56 kyr. In this case, the secondary collides with the primary at the end of the simulation. Or, in reality, the secondary is tidally disrupted and second-generation satellites/rings might be formed after the disruption event.

Figure D8 shows a tidal evolution pathway of a satellite system for the case of η_ref = 10¹⁰ Pa s. The initial condition of this model is R_X = 60 km, T_ini = 140 K, and e_ini = 0.4 (see Figure 8(d)). The final semilatus rectum is not very different from that of the standard case of η_ref = 10¹⁴ Pa s (see Figures D1(d), D3(d), and D5(d)), because the temperature of Gonggong, T_G, also depends on η_ref. In the case of η_ref = 10¹⁰ Pa s, T_G is lower than that for the case of η_ref = 10¹⁴ Pa s; and the viscosity of Gonggong, η_G = η(T_G), is not very different between the two settings, η_ref = 10¹⁰ Pa s and η_ref = 10¹⁴ Pa s.