HydroSyMBA: A 1D Hydrocode Coupled with an N-body Symplectic Integrator

The numerical study of a massive disk of particles, such as Saturn's rings, is a technical challenge. High-resolution modeling, in which one considers the individual particles composing the disk, is strongly limited by numerical capabilities. As a result, studies limit themselves to a small patch of the rings and perform integration over a few orbits (e.g., Wisdom & Tremaine 1988; Salo 1992, 1995). These studies are fundamental, as they provide a better understanding of the basic physics affecting the evolution of the rings.

One of the most basic processes driving the evolution of a disk of material is the redistribution of angular momentum via collisions. Particles on close, nearly circular orbits move at slightly different velocities: inner particles move faster and can then overtake and collide with outer particles. The angular momentum of the pair of particles is redistributed such that the outer (inner) particle gains (loses) angular momentum. Since angular momentum is proportional to the square root of the particles' distance to the planet, the outer (inner) particle moves outward (inward). As a result, the disk spreads (e.g., Lynden-Bell & Pringle 1974).

Collisions (or any friction that dissipates energy) thus cause a global transfer of angular momentum from the inner to the outer regions of the disk. By analogy with the transfer of angular momentum in a nonideal fluid, the efficiency of the transfer is modeled by a viscosity, whose value is generally a function of distance and local mass of the disk. This allows one to model the disk not by a collection of discrete objects, but as a continuous medium or fluid in which angular momentum is redistributed at a rate specified by an appropriate viscosity. The latter can be modeled by performing local simulations such as those mentioned above, or through analytical treatments.

When two objects collide within the Roche limit, they typically cannot remain bound solely due to gravity as the differential gravitational pull from the primary body is stronger than the gravity between the two objects (Canup & Esposito 1995). This prevents large objects from growing, and only small particles can exist, as their integrity is preserved by the strength of the material they are made of. As the disk spreads due to its viscosity, it will bring material through the Roche limit. This material can then coalesce by gravitational instabilities which can cause the particles from the disk to clump into larger objects (Goldreich & Ward 1973). Being outside the Roche limit, these clumps will not be sheared apart by the planet and can form new moonlets that separate from the disk. As these moonlets evolve inside the "satellite region" they may collide with each other and progressively grow into larger objects through two-body collisions, potentially forming larger satellites. This growth mechanism has been observed both in numerical (Charnoz et al. 2010) and analytical studies (the "pyramidal regime" of Crida & Charnoz 2012).

A third physical process that needs to be considered is the gravitational interactions between the disk and the orbiting moons. For an inner disk, these contribute to transferring angular momentum from the disk to the moons. This results in a contraction of the disk, and an expansion of the satellites' orbits. This occurs through a type of resonance called Lindblad resonances, which occur when the ratio of the orbital frequency Ω(r) at distance r in the disk is a rational multiple of that of the satellite Ω_s: pΩ(r) = qΩ_s, where p and q are integers greater than 1. First-order Lindblad resonances, which occur when q = p + 1, are the strongest, such that limiting the computation to these generally gives a good approximation to the system's behavior.

A fourth process that affects the satellites is tidal dissipation. A satellite will cause a deformation of the planet it orbits via the gravitational attraction it exerts. Since the planet is not a perfect fluid, the deformation will not occur instantaneously in response to the satellite's gravity. As a result, the satellite is not directly above the tidal bulge it generates on the planet. The deformation of the planet changes its gravitational potential, and the lag between the satellite and the position of the bulge creates an asymmetry that gives rise to a torque on the planet, and an equal and opposite torque on the satellite. If the satellite's orbit is such that its angular orbital velocity is faster than the planet's angular rotation velocity, then the satellite leads the bulge and it receives a negative torque, which causes its orbit to contract. The planet receives an equal and opposite torque that accelerates its spin. If the satellite orbits more slowly than the planet spins, then it receives a positive torque and its orbit expands, while the planet's spin is decelerated.

Additionally, the planet may also raise tides on the satellite ("satellite tides"). These tend to decrease the satellite's orbital eccentricity and drive its spin state to synchronous rotation (see Section 3.5).

The new integrator is composed of two components: a disk component whose evolution is computed using a hydrocode, and a satellite component using an N-body symplectic integrator. The disk is discretized on a one-dimensional grid (in the radial direction) comprised of N equal-size bins, and at each timestep its surface density is evolved by computing: (1) the mass fluxes between each cell of the grid due to the viscous torque, and (2) the variation of the disk's angular momentum due to resonant torques (if satellites are present). The N-body integrator is a modified version of the symplectic integrator SyMBA (Duncan et al. 1988; Salmon & Canup 2012, 2017). At each timestep it evolves the orbits of discrete objects due to tides, close encounters with other bodies, and resonant torques from the disk.

We have designed the integrator such that SyMBA handles the initialization of the system and serves as the main code to drive the model. At each timestep, the integrator calculates the evolution of the disk for half a timestep dt. The latter is performed with a second-order Runge–Kutta integrator, which evolves the surface density of the disk due to viscous spreading (using Equation (5)) and resonant interactions with satellites (using Equations (13) and (14)). During this operation, it calculates accelerations on each body due to resonant torques. If tides are included, accelerations on each particle due to this process are calculated as well. These accelerations are used to apply an additional "kick" to the momentum of the particles, before the core N-body integration of SyMBA is performed. We repeat the same operation for the latter half of the timestep. Thus the additional accelerations are added symmetrically around the main SyMBA operator, so that the disk and tidal accelerations are treated as small "kicks" to each object's motion (McNeil et al. 2005; Canup & Ward 2006). One can express our algorithm using the following operator form:

$$\begin{eqnarray}&&{E}_{\mathrm{Disk}}\left(\displaystyle \frac{{dt}}{2}\right){E}_{\mathrm{Tides}}\left(\displaystyle \frac{{dt}}{2}\right){E}_{\mathrm{SyMBA}}\left({dt}\right){E}_{\mathrm{Tides}}\left(\displaystyle \frac{{dt}}{2}\right){E}_{\mathrm{Disk}}\left(\displaystyle \frac{{dt}}{2}\right).\end{eqnarray} \tag{ 1 }$$

The equations of mass and angular momentum conservation read (e.g., Pringle 1981)

$$\begin{eqnarray}&&r\displaystyle \frac{\partial \sigma }{\partial t}+\displaystyle \frac{\partial \left(r\sigma {v}_{r}\right)}{\partial r}=0\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \left(\sigma {r}^{2}{\rm{\Omega }}\right)}{\partial t}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial \left(\sigma {r}^{3}{\rm{\Omega }}{v}_{r}\right)}{\partial r}=\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left(\nu \sigma {r}^{3}\displaystyle \frac{\partial {\rm{\Omega }}}{\partial r}\right),\end{eqnarray} \tag{ 3 }$$

where σ, ν, and v_r are the disk's surface density, viscosity, and radial velocity, respectively. These quantities are functions of time t and distance r. Combining Equations (2) and (3), one gets the equation governing the radial evolution of the disk's surface density due to viscous spreading:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \sigma }{\partial t}=\displaystyle \frac{3}{r}\displaystyle \frac{\partial }{\partial r}\left(\sqrt{r}\displaystyle \frac{\partial \left(\nu \sigma \sqrt{r}\right)}{\partial r}\right).\end{eqnarray} \tag{ 4 }$$

The disk's surface density is evolved using a second-order Runge–Kutta integrator to solve Equation (4), using an approach similar to Salmon et al. (2010). At each timestep, the surface density in each bin 2 ≤ j ≤ (N − 1), centered at distance R_j, is evolved by combining the mass fluxes through each edge of the bin:

$$\begin{eqnarray}&&{\sigma }_{j}(t+{dt})={\sigma }_{j}(t)+\displaystyle \frac{{F}_{j+1}-{F}_{j}}{2\pi {R}_{j}{\rm{\Delta }}R}{dt},\end{eqnarray} \tag{ 5 }$$

where the mass flux F_j through the edge E_j = R_j − ΔR/2 of bin j is given by

$$\begin{eqnarray}&&{F}_{j}=6\pi \sqrt{{E}_{j}}\displaystyle \frac{\left({\nu }_{j}{\sigma }_{j}\sqrt{{R}_{j}}-{\nu }_{j-1}{\sigma }_{j-1}\sqrt{{R}_{j-1}}\right)}{{\rm{\Delta }}R}.\end{eqnarray} \tag{ 6 }$$

We then apply user-selected boundary conditions to set F₁ and ${F}_{N+1}$, e.g., "free-flow" $({{\boldsymbol{F}}}_{1}={{\boldsymbol{F}}}_{2},{{\boldsymbol{F}}}_{N+1}={{\boldsymbol{F}}}_{N})$, "stop-flow" $({{\boldsymbol{F}}}_{1}=0,{{\boldsymbol{F}}}_{N+1}=0)$, or "no-inflow" (${\boldsymbol{F}}$₁ = 0 if ${\boldsymbol{F}}$₂ > 0, ${\boldsymbol{F}}$_N+1 = 0 if ${\boldsymbol{F}}$_N < 0), which allows one to evolve the surface density in the innermost and outermost bins. It is possible to select different boundary conditions for the innermost and outermost bin. We have included various viscosity models in the integrator: viscosity from gravitational instabilities (Ward & Cameron 1978), multicomponent viscosity for particulate disks (Daisaka et al. 2001), two-phase protolunar disk (Thompson & Stevenson 1988), and α–prescription (Shakura & Sunyaev 1973).

Resonant interactions transfer angular momentum between the disk and orbiting bodies. The torque exerted by a satellite on a disk at an (m:m-1) inner Lindblad resonance is given by (Goldreich & Tremaine 1980):

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{m}=-\displaystyle \frac{{\pi }^{2}\sigma }{3\omega {{\rm{\Omega }}}_{{ps}}}{\left[r\displaystyle \frac{d{{\rm{\Phi }}}_{m}}{{dr}}+\displaystyle \frac{2{\rm{\Omega }}}{\omega -{{\rm{\Omega }}}_{{ps}}}{{\rm{\Phi }}}_{m}\right]}^{2},\end{eqnarray} \tag{ 7 }$$

where σ is the disk's surface density, ω is the orbital frequency at distance r, Ω_ps is the pattern speed of the resonance and Φ_m is the mth-order Fourier component of the satellite's gravitational potential. In our integrator we limit ourselves to first-order resonances, which are the strongest, and for which Ω_ps is equal to the satellite's orbital frequency Ω_S. The satellite's potential can then be written (Goldreich & Tremaine 1978):

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{m}=-\displaystyle \frac{{{GM}}_{S}}{{a}_{s}}{b}_{1/2}^{(m)}(\alpha ),\end{eqnarray} \tag{ 8 }$$

where M_S and a_S are the satellite's mass and semimajor axis, $\alpha =r/{a}_{S}={\left(1-1/m\right)}^{2/3}$, and ${b}_{1/2}^{(m)}(\alpha )$ is the Laplace coefficient of order 1/2 defined by

$$\begin{eqnarray}&&{b}_{i}^{(m)}(\alpha )=\displaystyle \frac{2}{\pi }{\int }_{0}^{\pi }\displaystyle \frac{\cos (m\theta )d\theta }{{\left(1-2\alpha \cos \theta +{\alpha }^{2}\right)}^{i}}\end{eqnarray} \tag{ 9 }$$

At a first-order resonance location we have Ω = mΩ_S/(m − 1) and the resonant torque can be written

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{m}=\displaystyle \frac{{\pi }^{2}}{3}\displaystyle \frac{{M}_{S}^{2}}{{M}_{P}}G\sigma {a}_{s}{c}_{m},\end{eqnarray} \tag{ 10 }$$

where M_P is the mass of the central planet, G is the gravitational constant, and ${c}_{m}={\alpha }^{3/2}{\left(\alpha \tfrac{{{db}}_{1/2}^{(m)}}{d\alpha }+2{{mb}}_{1/2}^{(m)}\right)}^{2}$.

We then use the following approximation (Goldreich & Tremaine 1980)

$$\begin{eqnarray}\begin{array}{rcl}\left(\alpha \displaystyle \frac{{{db}}_{1/2}^{(m)}}{d\alpha }+2{{mb}}_{1/2}^{(m)}\right) & \approx & \displaystyle \frac{2m}{\pi }\left[{K}_{1}\left(\displaystyle \frac{2}{3}\right)+2{K}_{0}\left(\displaystyle \frac{2}{3}\right)\right]\\ & \approx & \displaystyle \frac{2m}{\pi }2.51,\end{array}\end{eqnarray} \tag{ 11 }$$

where K₀ and K₁ are modified Bessel functions, so that ${c}_{m}\approx 2.55{m}^{2}\left(1-1/m\right)$ (see also Salmon & Canup 2012).

To implement the effect of resonant torques on the disk, we use a different approach than that of Charnoz et al. (2010). The disk's angular momentum L_d can be written as the sum of the angular momentum of each of its cells ${L}_{d}={\sum }_{j=1}^{N}{\sigma }_{j}{S}_{j}\sqrt{{{GM}}_{P}{R}_{j}}$, where S_j and R_j are the surface area and central position of bin number j. To apply the torque on the disk, we first determine the bin j in which the resonance falls. For an inner resonance, the torque will cause an inward flux of mass from bin j to bin $j-1$. The resulting change in the disk angular momentum ΔL_m due to the (m: m − 1) inner Lindblad resonance is

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{L}_{m} & = & {L}_{d}^{{\prime} }-{L}_{d}=\sqrt{{{GM}}_{P}}\left[\left({\sigma }_{j}^{{\prime} }-{\sigma }_{j}\right){S}_{j}\sqrt{{R}_{j}}\right.\\ & & +\left.\left({\sigma }_{j-1}^{{\prime} }-{\sigma }_{j-1}\right){S}_{j-1}\sqrt{{R}_{j-1}}\right],\end{array}\end{eqnarray} \tag{ 12 }$$

where primed symbols indicate resulting quantities after the torque is applied. Conservation of mass implies that $\left({\sigma }_{j}^{{\prime} }-{\sigma }_{j}\right){S}_{j}\,=-\left({\sigma }_{j-1}^{{\prime} }-{\sigma }_{j-1}\right){S}_{j-1}$. Finally, using Γ_m = dL/dt ≈ ΔL_m/Δt, the resulting surface density of bin j can be written as

$$\begin{eqnarray}&&{\sigma }_{j}^{{\prime} }=-\displaystyle \frac{{{\rm{\Gamma }}}_{m}{\rm{\Delta }}t}{\sqrt{{{GM}}_{P}}{S}_{j}\left(\sqrt{{R}_{j}}-\sqrt{{R}_{j-1}}\right)}+{\sigma }_{j}.\end{eqnarray} \tag{ 13 }$$

Similarly, we can write that

$$\begin{eqnarray}&&{\sigma }_{j-1}^{{\prime} }=-\displaystyle \frac{{{\rm{\Gamma }}}_{m}{\rm{\Delta }}t}{\sqrt{{{GM}}_{P}}{S}_{j-1}\left(\sqrt{{R}_{j-1}}-\sqrt{{R}_{j}}\right)}+{\sigma }_{j-1}.\end{eqnarray} \tag{ 14 }$$

We are assuming that the torque from a given resonance is entirely deposited locally. In reality, the resonant torque will generate a wave that will propagate away from the perturber. However, the wave is progressively damped due to the disk's viscosity. For example, in the case of a viscous disk such as Saturn's rings, most of the torque is deposited in the region ξ < 2, where ξ = ^−1/2ΔR_L/R_L and ΔR_L is the distance to the resonance located at R = R_L (Shu et al. 1985). The factor is a function of the local surface density and distance of the resonance. In Saturn's rings, is typically 10⁻⁸ (Shu 1984), such that most of the torque is deposited within ΔR_L < 2 × 10⁻⁴R_L.

To determine the evolution of the satellite's orbit due to the resonant torque, we follow the approach of Papaloizou & Larwood (2000; as in Salmon & Canup 2012). We define an orbital migration timescale t_mig = L_S/Γ_m, where L_S is the satellite's orbital angular momentum. This is used to calculate an acceleration on the satellite ${\boldsymbol{a}}$_mig = ${\boldsymbol{v}}$/t_mig, where ${\boldsymbol{v}}$ is the satellite's velocity. The disk torque on the satellite results in an additional "kick" that modifies the velocity following ${\boldsymbol{v}}$' = ${\boldsymbol{v}}$ + ${\boldsymbol{a}}$_migΔt.

Viscous spreading can bring disk material through the Roche limit, whose position a_R is provided to the integrator as an input parameter. When this happens, we trigger the formation of a new N-body particle. The mass of the particle is set to the larger of the two following quantities: (1) the entire mass in the cell containing the Roche limit, or (2) the mass of the fragment that would form from gravitational instabilities (Goldreich & Ward 1973)

$$\begin{eqnarray}&&{m}_{f}=\displaystyle \frac{16{\pi }^{4}{\xi }^{3}{\sigma }_{{\rm{R}}}^{3}{a}_{{\rm{R}}}^{6}}{{M}_{P}^{2}},\end{eqnarray} \tag{ 15 }$$

where σ_R is the disk surface density at the Roche limit, and ξ is a factor of order, but less than, unity. We use the largest of the two criteria to prevent the creation of an unmanageable number of small bodies. The corresponding mass is removed from the disk's cell containing the Roche limit (and sometimes inner neighboring cells when the second criterion is used), and the new body's semimajor axis is set so as to conserve the angular momentum of the system. Finally, we set the the particle's eccentricity to be the ratio of its escape velocity to the local orbital velocity (Lissauer & Stewart 1993; Salmon & Canup 2012), and its inclination to zero (expected for formation from an equatorial disk).

We have also implemented a (direct accretion) mode to replicate the "continuous regime" from Crida & Charnoz (2012) in which a satellite orbiting close to a disk could grow from direct accretion of disk material. If this option is enabled, when disk material is being brought to the Roche limit we look for the existence of a nearby satellite. The threshold we have selected is to require the disk's edge to be within two Hill radius of the satellite, as in Crida & Charnoz (2012). If there is no such moon, then we proceed as previously explained. If there is such a qualifying satellite, then we add the disk mass directly onto it and adjust its orbital parameters to conserve angular momentum.

The acceleration on a satellite of mass M_S due to tides raised by the satellite on a planet with mass M_P, radius R_P, spin vector ${\boldsymbol{s}}$_P, and Love number k_2P is given by (Mignard 1980; Touma & Wisdom 1994):

$$\begin{eqnarray}\begin{array}{rcl}{\left.\displaystyle \frac{{\partial }^{2}{\boldsymbol{r}}}{\partial {t}^{2}}\right|}_{P} & = & -\displaystyle \frac{3{k}_{2P}{{GMR}}_{P}^{5}}{{r}^{10}}\left(1+\displaystyle \frac{{M}_{S}}{{M}_{P}}\right){\rm{\Delta }}{t}_{P}\\ & & \times \left[2\left({\boldsymbol{r}}\cdot {\boldsymbol{v}}\right){\boldsymbol{r}}+{r}^{2}\left({\boldsymbol{r}}\times {{\boldsymbol{s}}}_{P}+{\boldsymbol{v}}\right)\right],\end{array}\end{eqnarray} \tag{ 16 }$$

where ${\boldsymbol{r}}=\left(x,y,z\right)$ and ${\boldsymbol{v}}=\left({v}_{x},{v}_{y},{v}_{z}\right)$ are the satellite's position and velocity in a reference frame centered on the planet. The lag time Δt_P is defined as the fixed time between the tide-raising potential and the rise of the tidal distortion. This parameter is planet-dependent and is supplied as an input to the integrator.

We evolve the spin of the planet due to the tidal torque exerted by the satellite. This torque can be written (Mignard 1980; Touma & Wisdom 1994)²

$$\begin{eqnarray}&&{\boldsymbol{T}}=-\displaystyle \frac{3{k}_{2P}{{GM}}_{P}^{2}{R}_{P}^{5}}{{r}^{8}}{\rm{\Delta }}{t}_{P}\left[\left({\boldsymbol{r}}\cdot {{\boldsymbol{s}}}_{P}\right){\boldsymbol{r}}-{r}^{2}{{\boldsymbol{s}}}_{P}+{\boldsymbol{r}}\times {\boldsymbol{v}}\right].\end{eqnarray} \tag{ 17 }$$

Over a time step dt, the tides kick the planet's spin angular momentum by ΔL_P = Tdt, where $T=\parallel {\boldsymbol{T}}\parallel$. Using ${\rm{\Delta }}{L}_{P}\,={\lambda }_{P}{M}_{P}{R}_{P}^{2}{\rm{\Delta }}{s}_{P}$, where λ_P is the planet's moment of inertia constant, the variation of the planet's spin due to tides over dt is

$$\begin{eqnarray}&&{\rm{\Delta }}{s}_{P}=\displaystyle \frac{{Tdt}}{{\lambda }_{P}{M}_{P}{R}_{P}^{2}}.\end{eqnarray} \tag{ 18 }$$

The parameter λ_P is supplied to the integrator as an input parameter. The planet's radius is usually set constant, but it can also be progressively adjusted to, e.g., account for the planet contracting (as in Salmon & Canup 2017), over a timescale supplied in input. In such cases, the planet's spin is also adjusted to conserve the rotational angular momentum.

Similarly, the acceleration on the satellite due to tides raised by the planet on the satellite can be written as (Mignard 1980; Touma & Wisdom 1994):

$$\begin{eqnarray}\begin{array}{rcl}{\left.\displaystyle \frac{{\partial }^{2}{\boldsymbol{r}}}{\partial {t}^{2}}\right|}_{S} & = & -\displaystyle \frac{3{{GM}}_{P}{k}_{2P}{R}_{P}^{5}}{{r}^{10}}\left(1+\displaystyle \frac{{M}_{S}}{{M}_{P}}\right){ \mathcal A }{\rm{\Delta }}{t}_{P}\\ & & \times \left[2\left({\boldsymbol{r}}\cdot {\boldsymbol{v}}\right){\boldsymbol{r}}+{r}^{2}\left({\boldsymbol{r}}\times {{\boldsymbol{s}}}_{S}+{\boldsymbol{v}}\right)\right],\end{array}\end{eqnarray} \tag{ 19 }$$

where ${\boldsymbol{s}}$_S is the satellite's spin vector. The factor ${ \mathcal A }$ reflects the strength of satellite versus planetary tides, with

$$\begin{eqnarray}&&{ \mathcal A }={\left(\displaystyle \frac{{M}_{S}}{{M}_{P}}\right)}^{-2}{\left(\displaystyle \frac{{R}_{S}}{{R}_{P}}\right)}^{5}\left(\displaystyle \frac{{k}_{2S}}{{k}_{2P}}\right)\left(\displaystyle \frac{{\rm{\Delta }}{t}_{S}}{{\rm{\Delta }}t}\right),\end{eqnarray} \tag{ 20 }$$

where k_2S, Δt_S, and R_S are the satellite's Love number, tidal time lag, and physical radius. When implementing the latter equation in our integrator, we assume that the satellite is rotating synchronously, so that ${s}_{S}\approx n=\sqrt{{{GM}}_{P}/{a}^{3}}$, where a is the satellite's semimajor axis. This is a reasonable assumption, given that despinning timescales are generally short compared to tidal evolution timescales (e.g., Salmon & Canup 2017).

The overall integration timestep dt is set initially by SyMBAs requirement to be no larger than 1/20th of the innermost orbital period to be considered. The disk's evolution is calculated for half a time step dt/2 before and after the N-body integration is performed. However, this timestep can be too big for the hydrocode, and can lead to instabilities and negative surface densities. We thus allow for subdividing the timestep for the evolution of the disk calculation. The diffusive stability constraint (DSC) requires the timestep to be no larger than

$$\begin{eqnarray}&&{{dt}}_{\mathrm{DSC}}=\displaystyle \frac{1}{2}\displaystyle \frac{{\rm{\Delta }}{R}^{2}}{{\nu }_{\mathrm{Max}}},\end{eqnarray} \tag{ 21 }$$

where ΔR is the grid's radial resolution, and ν_Max is the maximum value of the viscosity throughout the entire disk (Johnston & Liu 2004). Through our testing we found that limiting the timestep of the hydrocode to dt_hydro = 0.1 ∗ dt_DSC ensures stability of the integrator. In practice, at each iteration we calculate dt_hydro. If it is larger than dt/2, we set dt_hydro = dt/2. Otherwise, we perform an integration using dt_hydro and repeat the operation until the disk's evolution has been computed for a full dt/2.

4.1.1. Constant Viscosity

We perform a first series of tests in which we only consider the viscous spreading of a disk, in the absence of any orbiting bodies. In order to assess the accuracy of the code, we perform a first test using a constant and uniform viscosity model and compare the disk's evolution to the analytical result from Pringle (1981). The latter is a solution for an initial Dirac distribution centered at r = r₀, which we cannot exactly replicate due to the finite resolution of the grid. We approximate this distribution by setting the disk's surface density to 0 except in the bin i containing r = r₀, where we set σ = M_d/S_i. We set the primary to be the Earth, and the disk mass is equal to one Lunar mass $\left({M}_{{\rm{L}}}\right)$. We use a radial grid extending from 0 to $2.4\,{R}_{\oplus }$ with 240 equal-sized bins. The boundary conditions are set to F₁ = F₂ and F_N+1 = F_N, which allows material to flow freely through the grid's edges. We set the Roche limit to an artificially large value in order to deactivate the creation of moonlets.

The surface density as a function of time and distance is given by Pringle (1981):

$$\begin{eqnarray}&&\sigma \left(x,\tau \right)=\displaystyle \frac{{M}_{d}}{\pi {r}_{0}^{2}}{\tau }^{-1}{x}^{-1/4}{\exp }^{-\left(1+{x}^{2}\right)/\tau }{I}_{1/4}\left(2x/\tau \right),\end{eqnarray} \tag{ 22 }$$

where x = r/r₀, $\tau =12\nu t/{r}_{0}^{2}$, and I_1/4 is a modified Bessel function of the first kind. Figure 2 shows the surface density profile calculated by the integrator in black, and the analytical solution from Equation (22) in red dashes, at four times of evolution. Our results are in excellent agreement with the analytical solution.

Zoom In Zoom Out Reset image size

4.1.2. Variable Viscosity

We use a similar setup to the previous section, but this time start with a Gaussian surface density distribution between 4 and 6 R_⊕ (Figure 3, panel (a)). We use a radial grid extending from 1 to 11 R_⊕. We perform four simulations with 1000, 2000, 4000, and 8000 grid cells, which correspond to radial resolutions ΔR of 10⁻², 5 × 10⁻³, 2.5 × 10⁻³, and $1.25\,\times \,{10}^{-3}\,{R}_{\oplus }$, respectively. We use a gravitational instability viscosity model (Ward & Cameron 1978):

$$\begin{eqnarray}&&{\nu }_{\mathrm{WC}}=\displaystyle \frac{{\pi }^{2}{G}^{2}{\sigma }^{2}}{{{\rm{\Omega }}}^{3}}.\end{eqnarray} \tag{ 23 }$$

Zoom In Zoom Out Reset image size

We integrate the system for 5 × 10⁵T_K, where T_K is the orbital period at 1 R_⊕. This integration time corresponds to about 3 × 10⁶, 1.5 × 10⁷, 6 × 10⁷ and 2.4 × 10⁸ times the DSC timestep at t = 0 for the four grid resolutions. Figure 3 shows the disk's surface density at t = 0 in panel (a), and at t = 5 × 10⁵T_K in panel (b), for the N = 1000 case. Panels (c) and (d) show the angular momentum error L_err and mass error M_err, which we define as

$$\begin{eqnarray}&&{L}_{\mathrm{err}}=\displaystyle \frac{{L}_{d}-{L}_{d,0}+{L}_{{pl}}+{L}_{\mathrm{out}}}{{L}_{d,0}},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{M}_{\mathrm{err}}=\displaystyle \frac{{M}_{d}-{M}_{d,0}+{M}_{{pl}}+{M}_{\mathrm{out}}}{{M}_{d,0}},\end{eqnarray} \tag{ 25 }$$

where L_d,0 and M_d,0 are the disk's initial angular momentum and mass, and L_pl and L_out are the angular momentum of the mass falling on the planet M_pl and being removed through the outer edge M_out, respectively. The error in the angular momentum of the disk depends strongly on the resolution, but is small in all cases. At t = 5 × 10⁵T_K, the angular momentum error is 1.55 × 10⁻⁷, 3.88 × 10⁻⁸, 9.72 × 10⁻⁹, and 2.43 × 10⁻⁹ for the grids with 1000, 2000, 4000, and 8000 bins, respectively. The error decreases as N², i.e., increasing the resolution by a factor of 2 decreases the angular momentum error by a factor of 4.

We find that the integrator conserves mass to machine precision with all resolutions, except the N = 8000 cells case which is significantly worse than the other three coarser resolutions. We believe this is a consequence of rounding errors due to the integration timestep becoming too small in order to satisfy the DSC (which is a function of ΔR²). Our algorithm is currently written in double precision. While implementing quadruple precision might allow for finer grid resolutions (which might be desirable when satellites are included in order to accurately track deposition of the resonant torque in the disk), we believe that the accuracy of the integrator with slightly coarser resolutions is satisfactory for the most desired purposes.

In this second test we consider a slab of material extending from 2 to 7 R_⊕around an Earth-mass planet, with an initial mass of 1 M_L. We consider radial grids identical to the previous case, but do not consider the N = 8000 bins case. Our boundary conditions are again set to let material flow freely through the inner and outer edges, and we again use the viscosity model from Ward & Cameron (1978). The Roche limit is again set to an artificially large value in order to deactivate the creation of moonlets. We place a 0.1M_L satellite at an initial semimajor axis of 8 R_⊕. We turn on the tidal dissipation in the planet, using a spin period of 5 hr (typical of a post-giant impact Earth), and a tidal lag ${\rm{\Delta }}{t}_{P}=6.9\times {10}^{6}\,{\rm{s}}$ (about 1000 times the current tidal lag of the Earth), and integrate the system for 5 × 10⁵T_K. Figure 4 shows the disk at t = 0 and 1000 T_K in panels (a) and (b). In the latter, note how the disk's surface density is sculpted by the satellite's resonances, plotted as the vertical red lines. Panel (c) shows the evolution of the satellite's semimajor axis and position of its 2:1 resonance. The satellite moves outward rapidly until it reaches ∼17.46 R_⊕, at which point its 2:1 resonance lies just outside of the disk's outer edge at 11 R_⊕. It then continues to move outward slowly due to the primary's tidal torque. The final semimajor axis of the satellite is $17.691\,{R}_{\oplus }$, $17.698\,{R}_{\oplus }$, and $17.702\,{R}_{\oplus }$ with N = 1000, 2000, and 4000 bins.

Zoom In Zoom Out Reset image size

Panels (d) and (e) show conservation of the angular momentum and mass of the system (disk + satellite). In this case, we define the angular momentum error as ${L}_{\mathrm{err}}=\left({L}_{\mathrm{tot}}(t)-{L}_{\mathrm{tot},0}\right)/{L}_{\mathrm{tot},0}$, with ${L}_{\mathrm{tot}}={L}_{d}+{L}_{{pl}}+{L}_{\mathrm{out}}+{L}_{P}+{L}_{\mathrm{orb}}$, where L_P is the spin angular momentum of the planet and L_orb is the orbital angular momentum of the satellite. Our integrator again conserves mass down to machine precision with all three resolutions. Conservation of angular momentum is similar to the purely viscous case. For comparison, the widely used code Fargo conserves angular momentum in the heliocentric frame exactly by construction, but the total angular momentum in the barycentric frame only to ∼10⁻³ precision over ∼2500 T_K in the case of a Jupiter mass planet embedded in a protoplanetary disk (Crida et al. 2007).

We perform a final test in which we model the evolution of a protolunar disk, and the accumulation of the Moon from disk material. The Roche-interior disk contains $1.5\,{M}_{L}$ between 1.2 and $2.9\,{R}_{\oplus }$, with a power-law surface density profile σ(r) ∝ r⁻³. The radial grid contains 200 bins extending from 1 to $3\,{R}_{\oplus }$, corresponding to a resolution ${\rm{\Delta }}R=\,{10}^{-2}\,{R}_{\oplus }$. The outer disk contains 2000 bodies following a cumulative size-frequency distribution N(>R) ∝ R^−1.5, totaling $0.5\,{M}_{L}$ and extending from 3 to $6\,{R}_{\oplus }$. Figure 5(a) shows the initial disk. We use a viscosity model appropriate for a two-phase protolunar disk (Thompson & Stevenson 1988; Salmon & Canup 2012) $\nu =\min \left({\nu }_{\mathrm{WC}},{\nu }_{\mathrm{TS}}\right)$, with

$$\begin{eqnarray}&&{\nu }_{\mathrm{TS}}=\displaystyle \frac{{\sigma }_{\mathrm{SB}}{T}_{P}^{4}}{\sigma {{\rm{\Omega }}}^{2}},\end{eqnarray} \tag{ 26 }$$

where σ_SB is the Stefan–Boltzmann constant and T_P is the vapor disk's photospheric temperature. We use ${T}_{P}=2000\,{\rm{K}}$, appropriate for a vertically unstratified silicate vapor disk (Thompson & Stevenson 1988). Finally, we do not include tidal interactions.

Zoom In Zoom Out Reset image size

As found in previous works (e.g., Salmon & Canup 2012), the accretion of the material in the outer disk occurs in less than a year (Figure 5b), while delivery of inner-disk material occurs over decades. We find that as the growing Moon scatters inner moonlets, it gains significant angular momentum, which expands its semimajor axis well beyond what could be achieved solely via direct disk–satellite interactions ($\sim 4.6\,{R}_{\oplus }$). Capture into mean-motion resonances is also efficient, with the two bodies in Figure 5(d) in 2:1 resonance.