HydroSyMBA: a 1D hydrocode coupled with an N-body symplectic integrator

TL;DR

HydroSyMBA is a novel numerical model that combines a 1D hydrocode with an N-body integrator to simulate the co-evolution of circumplanetary disks and satellites, enabling detailed studies of satellite formation and orbital dynamics.

Contribution

It introduces a self-consistent framework integrating a hydrocode with an N-body integrator, including features like moonlet spawning and tidal effects, for modeling satellite-disk systems.

Findings

Accurately models the evolution of circumplanetary disks and satellite formation.

Incorporates resonant torques and tidal dissipation effects.

Enables simulation of diverse planetary systems and satellite origins.

Abstract

The numerical modeling of co-existing circumplanetary disks/rings and satellites is particularly challenging because each part of the system requires a very different approach. Disks are generally well represented by a fluid-like dense medium, whose evolution can be calculated by a hydrocode. On the other hand, the orbital evolution of satellites is generally performed using N-body integrators. We have developed a new numerical model that combines a 1-dimensional hydrocode with the N-body integrator SyMBA. The disk evolves due to its viscosity, and resonant torques from satellites. The latter is applied to the satellites as an additional "kick" to their accelerations. The integrator also includes the ability to spawn new moonlets at the disk's outer edge if the latter expands beyond a material-dependent Roche limit, as well as the effects of tidal dissipation in the planet and/or the satellite on the satellite orbits. The resulting integrator allows one to accurately model the evolution of an inner circumplanetary disk, and the formation of satellites by accumulation of disk material, all within a single self-consistent framework. Potential applications include the formation of Earth's Moon, the evolution of the inner Saturn system, the martian and uranian moons, and compact exoplanetary systems.

Equations

$$E_{\mathrm{Disk}}\left(\frac{dt}{2}\right)E_{\mathrm{Tides}}\left(\frac{dt}{2}\right)E_{\mathrm{SyMBA}}\left(dt\right)E_{\mathrm{Tides}}\left(\frac{dt}{2}\right)E_{\mathrm{Disk}}\left(\frac{dt}{2}\right).$$

$$r\frac{\partial \sigma }{\partial t}+\frac{\partial \left(r\sigma v_r\right)}{\partial r}=0$$

$$\frac{\partial \left(\sigma r^2\Omega \right)}{\partial t}+\frac{1}{r}\frac{\partial \left(\sigma r^3\Omega v_r\right)}{\partial r}=\frac{1}{r}\frac{\partial }{\partial r}\left(\nu \sigma r^3\frac{\partial \Omega }{\partial r}\right)$$

$$\frac{\partial \sigma }{\partial t}=\frac{3}{r}\frac{\partial }{\partial r}\left(\sqrt{r}\frac{\partial \left(\nu \sigma \sqrt{r}\right)}{\partial r}\right)$$

$${\sigma }_j(t+dt)={\sigma }_j(t)+\frac{F_{j+1}-F_j}{2\pi R_j\Delta R}dt.$$

$$F_j=6\pi \sqrt{E_j}\frac{\left({\nu }_j{\sigma }_j\sqrt{R_j}-{\nu }_{j-1}{\sigma }_{j-1}\sqrt{R_{j-1}}\right)}{\Delta R}.$$

$${\Gamma }_m=-\frac{{\pi }^2\sigma }{3\omega {\Omega }_{ps}}{\left[r\frac{d{\Phi }_m}{dr}+\frac{2\Omega }{\omega -{\Omega }_{ps}}{\Phi }_m\right]}^2,$$

$${\Phi }_m=-\frac{G{M}_S}{a_s}{b}_{1/2}^{(m)}(\alpha ),$$

$$b_i^{(m)}(\alpha )=\frac{2}{\pi }{\int }_0^{\pi }\frac{\cos (m\theta )d\theta }{{\left(1-2\alpha \cos \theta +{\alpha }^2\right)}^i}$$

$${\Gamma }_m=\frac{{\pi }^2}{3}\frac{M_S^2}{M_P}G\sigma a_s{c}_m,$$

$$\left(\alpha \frac{d{b}_{1/2}^{(m)}}{d\alpha }+2m{b}_{1/2}^{(m)}\right)\approx \frac{2m}{\pi }\left[K_1\left(\frac{2}{3}\right)+2K_0\left(\frac{2}{3}\right)\right]\approx \frac{2m}{\pi }2.51,$$

$$\Delta L_m=L_d^{\prime }-L_d=\sqrt{G{M}_P}\left[\left({\sigma }_j^{\prime }-{\sigma }_j\right)S_j\sqrt{R_j}+\left({\sigma }_{j-1}^{\prime }-{\sigma }_{j-1}\right)S_{j-1}\sqrt{R_{j-1}}\right],$$

$${\sigma }_j^{\prime }=-\frac{{\Gamma }_m\Delta t}{\sqrt{G{M}_P}{S}_j\left(\sqrt{R_j}-\sqrt{R_{j-1}}\right)}+{\sigma }_j.$$

$${\sigma }_{j-1}^{\prime }=-\frac{{\Gamma }_m\Delta t}{\sqrt{G{M}_P}{S}_{j-1}\left(\sqrt{R_{j-1}}-\sqrt{R_j}\right)}+{\sigma }_{j-1}.$$

$$m_f=\frac{16{\pi }^4{\xi }^3{\sigma }_R^3{a}_R^6}{M_P^2},$$

$${\frac{{\partial }^2\mathbf{r}}{\partial t^2}|}_P=-\frac{3k_{2P}GM{R}_P^5}{r^{10}}\left(1+\frac{M_S}{M_P}\right)\Delta t_P\left[2\left(\mathbf{r}\cdot \mathbf{v}\right)\mathbf{r}+r^2\left(\mathbf{r}\times {\mathbf{s}}_{\mathbf{P}}+\mathbf{v}\right)\right],$$

$$\mathbf{T}=-\frac{3k_{2P}G{M}_P^2{R}_P^5}{r^8}\Delta t_P\left[\left(\mathbf{r}\cdot {\mathbf{s}}_{\mathbf{P}}\right)\mathbf{r}-r^2{\mathbf{s}}_{\mathbf{P}}+\mathbf{r}\times \mathbf{v}\right].$$

$$\Delta s_P=\frac{Tdt}{{\lambda }_P{M}_P{R}_P^2}.$$

$${\frac{{\partial }^2\mathbf{r}}{\partial t^2}|}_S=-\frac{3G{M}_P{k}_{2P}{R}_P^5}{r^{10}}\left(1+\frac{M_S}{M_P}\right)\mathcal{A}\Delta t_P\left[2\left(\mathbf{r}\cdot \mathbf{v}\right)\mathbf{r}+r^2\left(\mathbf{r}\times {\mathbf{s}}_{\mathbf{S}}+\mathbf{v}\right)\right],$$

$$\mathcal{A}={\left(\frac{M_S}{M_P}\right)}^{-2}{\left(\frac{R_S}{R_P}\right)}^5\left(\frac{k_{2S}}{k_{2P}}\right)\left(\frac{\Delta t_S}{\Delta t}\right),$$

$$d{t}_{DSC}=\frac{1}{2}\frac{\Delta R^2}{{\nu }_{Max}},$$

$$\sigma \left(x,\tau \right)=\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),$$

$${\nu }_{WC}=\frac{{\pi }^2{G}^2{\sigma }^2}{{\Omega }^3}$$

$$L_{err}=\frac{L_d-L_{d,0}+L_{pl}+L_{out}}{L_{d,0}}$$

$$M_{err}=\frac{M_d-M_{d,0}+M_{pl}+M_{out}}{M_{d,0}}$$

$${\nu }_{TS}=\frac{{\sigma }_{SB}{T}_P^4}{\sigma {\Omega }^2},$$

Full text

1 Introduction

The formation and the evolution of planetary satellites is a fundamental topic in planetary science. The history of these systems is closely tied to that of the planet they orbit, such that understanding satellites is an indirect way of studying planet formation in general, and the history of our Solar System in particular. The irregular satellites, such as Triton around Neptune, are likely objects that formed in the outer Solar System, were later scattered inward, passed close to the planet and got captured via gravitational interactions. The regular satellites likely formed from a disk around the planet itself.

Planets form inside the proto-planetary disk, a giant disk of gas and dust orbiting the Sun. If the planet grows massive enough, some gas and dust from the main nebula can become bound to the planet and orbit it, forming a subnebula. We believe that the Galilean Satellites formed by accretion of gas and dust in Jupiter’s subnebula (e.g. Canup & Ward, 2006). The Earth’s Moon likely formed by accretion from a disk, itself the product of a giant impact suffered by Earth toward the end of its formation (Cameron & Ward, 1976; Salmon & Canup, 2012). Some of Saturn’s satellites may be the children of Saturn’s rings, as they may have assembled by accumulation of material from rings that were much more massive in the past (Canup, 2010; Charnoz et al., 2010, 2011; Salmon & Canup, 2017). Phobos and Deimos, satellites of Mars, are often viewed as captured asteroids, but their equatorial, quasi-circular and low-inclination orbits point to formation from an ancient circum-Mars ring (Rosenblatt et al., 2016; Hesselbrock & Minton, 2017; Hyodo et al., 2017; Canup & Salmon, 2018). Uranus and Neptune’s regular satellites are also compatible with formation from a ring (Crida & Charnoz, 2012).

Of particular importance for these systems is the dynamical transition that separates an inner disk, or ring, from the region of growing satellites: the Roche limit. This distance is a function of the densities of the central body and the accreting disk material. Inside this distance, the differential gravitational force from the primary body across two colliding bodies is strong and generally prevents the bodies from “sticking” to each other by mutual gravity. As a result, only small objects, which are held together by the strength of the material they are made from rather than from their gravity, are generally present. Outside of the Roche limit, two colliding bodies can remain bound by their mutual gravity, allowing larger objects to grow, eventually forming satellites (Figure 1). A prime example is Saturn’s system. Saturn’s rings are composed primarily of water ice, setting the Roche limit at $\sim 140000$ km. Inside this distance lie the rings, composed of particles of order $1\text{\,}\mathrm{m}$. In the transition region around this distance lies the F ring, which has a spiral structure and within which the appearance and destruction of clumps of material are observed. Beyond this distance we find the many Saturnian satellites.

While circumplanetary disks and satellites appear to share a common history, these systems are inherently difficult to study in a single, self-consistent model, given the different numerical approaches necessary to correctly model the rings and the satellites. Modeling the entirety of a cold particulate disk, such as Saturn’s rings, by a collection of individual particles would require one to include billions of small objects, which would rapidly limit computational efficiency. However, if the disk is dense enough, then the evolution of the disk is dominated by collective effects, so that it becomes possible to model it as a continuous medium (Spahn & Schmidt, 2006), which can be represented on a discrete grid with a few $10^{2}$ to $10^{4}$ radial or annular cells , greatly improving numerical efficiency. Hydrocodes are generally the tool of choice to study the evolution of circumplanetary disks.

To model the growth of individual satellites, it is necessary to accurately integrate their orbital evolution, which is dominated by few-body gravitational interactions, with the planet and possibly other satellites. This is generally done using N-body integrators, which model each satellite individually and evolve its orbit due to the gravitational potential of the planet and that of other orbiting bodies during close encounters. Most state-of-the-art integrators can efficiently integrate up to 10000 bodies over $\sim 10^{7}$ orbits.

Due to the large difference in modeling needs, most numerical integrators use simplified approaches to model part of the system, which generally prohibits including important physical processes. For instance, the integrator used in Charnoz et al. (2010) did not explicitly model the interactions between growing satellites, such that important interactions such as capture in mean-motion resonances could not effectively occur. On the other hand, the integrator used in Salmon & Canup (2012) did not resolve the radial structure of the disk, such that the resonant torque from an outer moonlet was not deposited at the actual position of the resonance, but instead was applied at the disk’s outer edge, which would strongly affect the outward spreading of the disk in that region.

In an attempt to remedy the limitations of current integrators, we have developed a new code that combines a hydrocode component to model a continuous inner disk, and an N-body integrator to fully track the orbital evolution of discrete outer bodies. While Fargo-SyMBA (Morbidelli & Nesvorny, 2012) uses a similar approach, its numerical complexity and the absence of treatment of accretion processes at the Roche limit make it unsuitable for studying the formation of satellites from a spreading, Roche-interior circumplanetary disk. The disk section in our new model is broadly inspired from that in Charnoz et al. (2010), although we have largely rewritten the disk-satellite interaction portion to address issues with conservation of the system’s angular momentum. The N-body part is handled by a modified version of the symplectic integrator SyMBA (Duncan et al., 1998), in which the particle-particle interactions remain unchanged, but we have added additional accelerations on the bodies due to interactions with the disk and tidal dissipation.

The resulting code allows one to model the full system of rings/disk and satellites with a tailored approach for each part, and to include the most important dynamical processes. While in the following we focus our presentation on systems where the primary is a planet, the latter could also, for example, be a star in systems with very close-in orbiting exoplanets (e.g. van Lieshout et al., 2018), many of which have been identified in the last decade, in particular via the Kepler mission. In Section 2 we describe the details of the integrator. In Section 3 we present some test results, and in Section 4 we discuss some potential future developments.

2 Relevant physical processes

2.1 Viscous spreading of a disk

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. Salo, 1992, 1995; Wisdom & Tremaine, 1988). 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 non-ideal 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.

2.2 Accretion outside the Roche limit

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)

2.3 Disk-satellite interactions

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 $\Omega(r)$ at distance $r$ in the disk is a rational multiple of that of the satellite $\Omega_{s}$: $p\Omega(r)=q\Omega_{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.

2.4 Tides

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 which 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 that 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).

3 Numerical implementation.

3.1 Basic architecture of the integrator

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 1-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 Eq. 5) and resonant interactions with satellites (using Eq. 13 and Eq. 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:

[TABLE]

3.2 Viscous spreading

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

[TABLE]

where $\sigma$, $\nu$ 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:

[TABLE]

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\leq j\leq(N-1)$, centered at distance $R_{j}$, is evolved by combining the mass fluxes through each edge of the bin:

[TABLE]

where the mass flux $F_{j}$ through the edge $E_{j}=R_{j}-\Delta R/2$ of bin j is given by

[TABLE]

We then apply user-selected boundary conditions to set $F_{1}$ and $F_{N+1}$, e.g. “free-flow” $(\bf F_{1}=F_{2},F_{N+1}=F_{N})$, “stop-flow” $(\bf F_{1}=0,F_{N+1}=0)$, or “no-inflow” $(\bf F_{1}=0$ if $\bf F_{2}>0,F_{N+1}=0$ if $\bf 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), multi-component viscosity for particulate disks (Daisaka et al., 2001), two-phase protolunar disk (Thompson & Stevenson, 1988), and $\alpha-$prescription (Shakura & Sunyaev, 1973).

3.3 Disk-satellite interactions

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):

[TABLE]

where $\sigma$ is the disk’s surface density, $\omega$ is the orbital frequency at distance $r$, $\Omega_{ps}$ is the patern speed of the resonance and $\Phi_{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 $\Omega_{ps}$ is equal to the satellite’s orbital frequency $\Omega_{S}$. The satellite’s potential can then be written (Goldreich & Tremaine, 1978):

[TABLE]

where $M_{S}$ and $a_{S}$ are the satellite’s mass and semi-major 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:

[TABLE]

At a first-order resonance location we have $\Omega=m\Omega_{S}/(m-1)$ and the resonant torque can be written

[TABLE]

where $M_{P}$ is the mass of the central planet, $G$ is the gravitational constant and $c_{m}=\alpha^{3/2}\left(\alpha\frac{db_{1/2}^{(m)}}{d\alpha}+2mb_{1/2}^{(m)}\right)^{2}$.

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

[TABLE]

where $K_{0}$ and $K_{1}$ are modified Bessel functions, so that $c_{m}\approx 2.55m^{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 $\Delta L_{m}$ due to the $(m:m-1)$ inner Lindblad resonance is

[TABLE]

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 $\Gamma_{m}=dL/dt\approx\Delta L_{m}/\Delta t$, the resulting surface density of bin $j$ can be written as

[TABLE]

Similarly, we can write that

[TABLE]

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 $\xi<2$, where $\xi=\epsilon^{-1/2}\Delta R_{L}/R_{L}$ and $\Delta R_{L}$ is the distance to the resonance located at $R=R_{L}$ (Shu et al., 1985). The factor $\epsilon$ is a function of the local surface density and distance of the resonance. In Saturn’s rings, $\epsilon$ is typically $10^{-8}$ (Shu, 1984), such that most of the torque is deposited within $\Delta R_{L}<2\times 10^{-4}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}/\Gamma_{m}$, where $L_{S}$ is the satellite’s orbital angular momentum. This is used to calculate an acceleration on the satellite $\mathbf{a}_{mig}=\mathbf{v}/t_{mig}$, where $\mathbf{v}$ is the satellite’s velocity. The disk torque on the satellite results in an additional “kick” that modifies the velocity following $\mathbf{v}^{\prime}=\mathbf{v}+\mathbf{a}_{mig}\Delta t$.

3.4 Accretion outside the Roche limit

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)

[TABLE]

where $\sigma_{R}$ is the disk surface density at the Roche limit, and $\xi$ 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 semi-major 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.

3.5 Tides

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 $\mathbf{s_{P}}$ and Love number $k_{2P}$ is given by (Mignard, 1980; Touma & Wisdom, 1994):

[TABLE]

where $\mathbf{r}=\left(x,y,z\right)$ and $\mathbf{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 $\Delta 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)111There is a typo in Equation (111) of Touma & Wisdom (1994): the $r_{em}^{10}$ on the denominator should be $r_{em}^{8}$.

[TABLE]

Over a time step $dt$, the tides kick the planet’s spin angular momentum by $\Delta L_{P}=Tdt$, where $T=\left\|\mathbf{T}\right\|$. Using $\Delta L_{P}=\lambda_{P}M_{P}R_{P}^{2}\Delta s_{P}$, where $\lambda_{P}$ is the planet’s moment of inertia constant, the variation of the planet’s spin due to tides over $dt$ is

[TABLE]

The parameter $\lambda_{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):

[TABLE]

where $\mathbf{s_{S}}$ is the satellite’s spin vector. The factor $\mathcal{A}$ reflects the strength of satellite versus planetary tides, with

[TABLE]

where $k_{2S}$, $\Delta 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 semi-major axis. This is a reasonable assumption, given that despinning timescales are generally short compared to tidal evolution timescales (e.g. Salmon & Canup, 2017)

3.6 Control of the time step

The overall integration timestep $dt$ is set initially by SyMBA’s 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

[TABLE]

where $\Delta R$ is the grid’s radial resolution, and $\nu_{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 Tests

4.1 Viscous spreading

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_{0}$, 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_{0}$, where we set $\sigma=M_{d}/S_{i}$. We set the primary to be the Earth, and the disk mass is equal to one Lunar mass $\left(M_{L}\right)$. We use a radial grid extending from 0 to $2.4\text{\,}\mathrm{R}_{\oplus}$ with 240 equal-sized bins. The boundary conditions are set to $F_{1}=F_{2}$ 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 de-activate the creation of moonlets.

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

[TABLE]

where $x=r/r0$, $\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.

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 $6R_{\oplus}$ (Figure 3, panel (a)). We use a radial grid extending from 1 to $11R_{\oplus}$. We perform 4 simulations with 1000, 2000, 4000 and 8000 grid cells, which correspond to radial resolutions $\Delta R$ of $10^{-2}$, $5\times 10^{-3}$, $2.5\times 10^{-3}$, and $1.25\text{\times}{10}^{-3}\text{\,}\mathrm{R}_{\oplus}$, respectively. We use a gravitational instability viscosity model (Ward & Cameron, 1978):

[TABLE]

We integrate the system for $5\times 10^{5}T_{K}$, where $T_{K}$ is the orbital period at $1R_{\oplus}$. This integration time corresponds to about $3\times 10^{6}$, $1.5\times 10^{7}$, $6\times 10^{7}$ and $2.4\times 10^{8}$ times the DSC timestep at t=0 for the 4 grid resolutions. Figure 3 shows the disk’s surface density at $t=0$ in panel (a), and at $t=5\times 10^{5}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:

[TABLE]

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\times 10^{5}T_{K}$, the angular momentum error is $1.55\times 10^{-7}$, $3.88\times 10^{-8}$, $9.72\times 10^{-9}$ and $2.43\times 10^{-9}$ for the grids with 1000, 2000, 4000 and 8000 bins, respectively. The error decreases as $N^{2}$, 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 $\Delta R^{2}$). 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.

4.2 Disk and satellite

In this second test we consider a slab of material extending from 2 to $7R_{\oplus}$ around an Earth-mass planet, with an initial mass of $1M_{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 de-activate the creation of moonlets. We place a $0.1M_{L}$ satellite at an initial semi-major axis of $8R_{\oplus}$. We turn on the tidal dissipation in the planet, using a spin period of 5 hours (typical of a post-giant impact Earth), and a tidal lag $\Delta t_{P}=$6.9\text{\times}{10}^{6}\text{,}\mathrm{s}$$ (about 1000 times the current tidal lag of the Earth), and integrate the system for $5\times 10^{5}T_{K}$. Figure 4 shows the disk at t=0 and $1000T_{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 semi-major axis and position of its 2:1 resonance. The satellite moves outward rapidly until it reaches $\sim 17.46R_{\oplus}$, at which point its 2:1 resonance lies just outside of the disk’s outer edge at $11R_{\oplus}$. It then continues to move outward slowly due to the primary’s tidal torque. The final semi-major axis of the satellite is $17.691\text{\,}\mathrm{R}_{\oplus}$, $17.698\text{\,}\mathrm{R}_{\oplus}$ and $17.702\text{\,}\mathrm{R}_{\oplus}$ with N=1000, 2000 and 4000 bins.

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_{err}=\left(L_{tot}(t)-L_{tot,0}\right)/L_{tot,0}$, with $L_{tot}=L_{d}+L_{pl}+L_{out}+L_{P}+L_{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 3 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 $\sim 10^{-3}$ precision over $\sim 2500T_{K}$ in the case of a Jupiter mass planet embedded in a protoplanetary disk (Crida et al., 2007).

4.3 Lunar accretion

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\text{\,}\mathrm{M}_{\mathrm{L}}$ between 1.2 and $2.9\text{\,}\mathrm{R}_{\oplus}$, with a power-law surface density profile $\sigma(r)\propto r^{-3}$. The radial grid contains 200 bins extending from 1 to $3\text{\,}\mathrm{R}_{\oplus}$, corresponding to a resolution $\Delta R=${10}^{-2}\text{,}\mathrm{R}_{\oplus}$$. The outer disk contains 2000 bodies following a cumulative size-frequency distribution $N(>R)\propto R^{-1.5}$, totaling $0.5\text{\,}\mathrm{M}_{\mathrm{L}}$ and extending from 3 to $6\text{\,}\mathrm{R}_{\oplus}$. Figure 5a 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_{WC},\nu_{TS}\right)$, with

[TABLE]

where $\sigma_{SB}$ is the Stefan-Boltzmann constant and $T_{P}$ is the vapor disk’s photospheric temperature. We use $T_{P}=$2000\text{,}\mathrm{K}$$, appropriate for a vertically unstratified silicate vapor disk (Thompson & Stevenson, 1988). Finally, we do not include tidal interactions.

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 semi-major axis well beyond what could be achieved solely via direct disk-satellite interactions ($\sim$4.6\text{,}\mathrm{r}_{\oplus}$$). Capture into mean-motion resonances is also efficient, with the two bodies in Figure 5d in 2:1 resonance.

5 Summary

We have developed a new numerical integrator that allows one to model the evolution of a continuous disk and individual bodies, within a single, self-consistent environment. The evolution of the disk is computed using a second-order Runge-Kutta integration scheme, while the orbital evolution of the particles is handled by the symplectic N-body integrator SyMBA. The disk evolves due to its viscosity (with various viscosity models available), and with satellite interactions through first-order Lindblad resonances. The evolution of the satellites includes particle-particle interactions (including collision and mergers), back reaction from resonant interactions with the disk, and tidal dissipation in the primary and/or the satellite itself.

The current integrator represents a significant step forward in the modeling on circumplanetary disks and satellites. It is suitable to study a variety of systems, e.g.: a well-mixed non-stratified protolunar disk; the rings and satellites of Saturn; the formation of Phobos and Deimos from a circum-Mars disk; the formation of some of Uranus’s satellites (Morbidelli et al., 2012); the formation of the moons of exoplanets (Nakajima et al., 2014). We have also identified some additional features that would be useful to incorporate

• •

Tidal disruption In the current integrator, an object scattered inside the Roche limit remains intact. However, it may instead be tidally disrupted. Sridhar & Tremaine (1992) indicate that an object of density $\rho$ will be entirely disrupted in a single pass if its pericenter $q<1.05\left(M_{P}/\rho\right)^{(1/3)}$. This can be a significant process as shown in simulations of lunar accretion in Salmon & Canup (2012). Alternatively, a body would be absorbed by the inner disk if as it passed through the disk it encountered a mass comparable to or greater than its own.

• •

Layered disk Models of lunar accretion consider disks in which the vapor and melt phases are well-mixed and co-evolve (Salmon & Canup, 2012, 2014). However, recent works (Ward, 2012) argue that the condensed phase of the protolunar disk rapidly settles in the midplane, resulting in a disk of magma surrounded by a vapor atmosphere. The latter can be represented by two additional disks, above and below the condensed disk. The evolution of each phase needs then to be modeled separately as they each will respond to different physical processes, or to similar processes but with different characteristic parameters.

• •

Irregular grid It would be valuable to implement grids with variable resolution, as it would allow one to model more precisely certain regions of the disk (e.g. the outer region where most resonant interactions occur) at a reduced computational cost.

Acknowledgements

This work has been funded by Southwest Research Institute Internal Research program, and NASA’s Solar System Exploration Research Virtual Institute grant NNA14AB03A. We thank an anonymous referee for their careful reading of the paper and their valuable comments. Southwest Research Institute’s Internal Research program provided primary funding for the code’s development. While the code ir proprietary to SwRI, we are happy to make it available to other scientists in the context of mutually collaborative projects. We invite interested readers to contact us directly for such requests.