Giant planets and brown dwarfs on wide orbits: a code comparison project

TL;DR

This study compares various hydrodynamical codes in simulating the migration and growth of gas clumps in circumstellar discs, highlighting the impact of numerical methods and the dominance of physical uncertainties.

Contribution

It provides a systematic comparison of multiple simulation codes for clump evolution, revealing the effects of numerical choices and the greater influence of physical modeling uncertainties.

Findings

Quantitative agreement in migration rates varies from 10% to 50%.

Artificial viscosity and sink prescriptions significantly affect results.

Discrepancies among population synthesis models exceed those among hydrodynamical codes.

Abstract

Gas clumps formed within massive gravitationally unstable circumstellar discs are potential seeds of gas giant planets, brown dwarfs and companion stars. Simulations show that competition between three processes -- migration, gas accretion and tidal disruption -- establishes what grows from a given seed. Here we investigate the robustness of numerical modelling of clump migration and accretion with the codes PHANTOM, GADGET, SPHINX, SEREN, GIZMO-MFM, SPHNG and FARGO. The test problem comprises a clump embedded in a massive disc at an initial separation of 120 AU. There is a general qualitative agreement between the codes, but the quantitative agreement in the planet migration rate ranges from $\sim 10$% to $\sim 50$%, depending on the numerical setup. We find that the artificial viscosity treatment and the sink particle prescription may account for much of the differences between the codes. In order to understand the wider implications of our work, we also attempt to reproduce the planet evolution tracks from our hydrodynamical simulations with prescriptions from three previous population synthesis studies. We find that the disagreement amongst the population synthesis models is far greater than that between our hydrodynamical simulations. The results of our code comparison project are therefore encouraging in that uncertainties in the given problem are probably dominated by the physics not yet included in the codes rather than by how hydrodynamics is modelled in them.

Figures

Tables

Table 1. Table 1: The parameters distinguishing the Runs presented in this paper. R a subscript 𝑅 a R_{\rm a} , h s subscript ℎ 𝑠 h_{s} , and M p0 subscript 𝑀 p0 M_{\rm p0} are the sink accretion radius, the gravitational softening parameter, and the mass of the planet, respectively. All the other parameters and initial conditions are the same for all four Runs.

Run	$R_a$ (AU)	$h_s$ (AU)	$M_{\mathrm{p0}}$ ($M_J$)
Run 1	0.0	1	2
Run 2	0.0	2	2
Run 3	0.5	0.01	2
Run 4	1.0	0.01	12

Equations

$$T_{\mathrm{irr}}=T_0{\left(\frac{R}{R_0}\right)}^{-1/2},$$

$$u_{\mathrm{irr}}=\frac{k_B{T}_{\mathrm{irr}}}{\mu (\gamma -1)},$$

$$\frac{\mathrm{du}}{\mathrm{dt}}=-\frac{u-u_{\mathrm{irr}}}{t_{\mathrm{cool}}},$$

$${\Omega }_{\mathrm{K}}(R)={\left(\frac{G{M}_{\ast }}{R^3}\right)}^{1/2}$$

$${\Sigma }_{\mathrm{in}}(R)=\frac{M_{\mathrm{d}}}{2\pi R(R_{\mathrm{out}}-R_{\mathrm{in}})}$$

$$\tau =\frac{R}{|\dot{R}|},$$

$${\tau }_{\mathrm{iso}}={\left(2.7+1.1\lambda \right)}^{-1}\frac{M_{\star }}{M_p}\frac{M_{\star }}{{\Sigma }_{\mathrm{p}}{r}_p^2}{\left(\frac{c_s}{r_p{\Omega }_p}\right)}^2\frac{1}{{\Omega }_p}.$$

$${\tau }_{\mathrm{sg}}\approx \frac{5.6}{(3.8-\lambda )}\gamma Q_p\frac{h_p^3}{q}{\left(\frac{0.1}{h_p}\right)}^2\frac{2\pi }{{\Omega }_p},$$

$$\dot{R}(t)=\frac{R(t)-R_0}{t},$$

$$\dot{\bar{R}}(t)=\frac{1}{T_{\mathrm{p}}}{\int }_{t-T_{\mathrm{p}}/2}^{t+T_{\mathrm{p}}/2}\dot{R}(t^{\prime })\mathrm{d}{t}^{\prime },$$

$$R_{\mathrm{H}}=R{\left(\frac{M_{\mathrm{H}}}{3M_{\star }}\right)}^{1/3},$$

$$C_{\mathrm{p}}=\frac{3}{4}\frac{H}{R_{\mathrm{H}}}+\frac{50\alpha H^2}{R^2}\frac{M_{\ast }}{M_{\mathrm{p}}}\le 1.$$

$$\alpha =\frac{1}{10}{\alpha }_{\mathrm{v}}\frac{h_{\mathrm{sml}}}{H},$$

$$t_{acc}=\frac{M_{\mathrm{p}}}{\dot{M}},$$

$${\tau }_{\mathrm{acc}}=\frac{t_{\mathrm{acc}}}{T_{\mathrm{p}}},$$

$${\dot{M}}_{\mathrm{acc}}=b\frac{M_p}{M_{\ast }}{\Omega }_{\mathrm{p}}\rho R^3,$$

$$\Sigma =\frac{c_{\mathrm{s}}{\Omega }_{\mathrm{p}}}{\pi GQ},$$

$${\tau }_{\mathrm{acc}}=b^{-1}Q.$$

$$\dot{M}=4\Sigma \Omega R_H^2$$

$${\tau }_{acc}=\frac{M_{\mathrm{p}}}{\dot{M}{T}_{\mathrm{p}}}=\frac{3}{8}\frac{R_{\mathrm{H}}}{H}Q.$$

$$t_{\mathrm{I}}={\left(\frac{M_{\mathrm{p}}}{M_{\ast }}\Omega \right)}^{-1}\frac{H}{R}.$$

$$t_{\mathrm{II}}=t_{\mathrm{visc}}=\frac{1}{\alpha \Omega }{\left(\frac{H}{R}\right)}^{-2}.$$

$$M_{\mathrm{t}}=2M_{\ast }{\left(\frac{H}{R}\right)}^3$$

$$t_{\mathrm{I}}=f_{\mathrm{mig}}\frac{M_{\ast }^2}{M_{\mathrm{p}}{M}_{\mathrm{d}}}\frac{H^2}{R^2}{\Omega }^{-1}$$

$$t_{\mathrm{II}}=t_{\mathrm{visc}}\left(1+\frac{M_{\mathrm{p}}}{M_{\mathrm{d}}}\right).$$

$${\dot{M}}_{\mathrm{p}}={\dot{M}}_{\mathrm{acc}}{\left[1+e^{-(C_{\mathrm{p}}-1)/\Delta C}\right]}^{-1},$$

$$N_{\mathrm{present}}=N_{\mathrm{birth}}\times P_{\mathrm{surv}},$$

Full text

††pubyear: 2017††pagerange: Giant planets and brown dwarfs on wide orbits: a code comparison project.–Giant planets and brown dwarfs on wide orbits: a code comparison project.

1 Introduction

Secondary star formation via gravitational instability (GI) of massive circumstellar discs has now been observed by ALMA (Tobin et al., 2016) and may be a viable explanation for the high frequency and the host metallicity correlations of stellar binaries with separations less than tens of AU (Moe et al., 2018). Modern star formation simulations (Bate, 2018) and observations of young discs (Tychoniec et al., 2018) also indicate that massive large gas discs could be abundant.

The conditions for disc fragmentation (Gammie, 2001; Rafikov, 2005) are similar to those for forming first hydrostatic cores in star formation (Larson, 1969), implying that the masses of gas clumps born in the discs must be initially similar to those of the opacity-limited fragments, e.g., $\sim 5-10{\,{\rm M}_{\rm J}}$ (Low & Lynden-Bell, 1976; Rees, 1976; Masunaga et al., 1998), although both smaller and larger initial clump masses were considered in the literature (Boley et al., 2010; Kratter et al., 2010; Forgan & Rice, 2013a). Due to these uncertainties and due to strong clump evolution after formation via inward migration (Mayer et al., 2004; Vorobyov & Basu, 2005; Machida et al., 2011; Baruteau et al., 2011), gas accretion (Zhu et al., 2012a; Stamatellos, 2015; Mercer & Stamatellos, 2017) and tidal disruption (Boley et al., 2010; Nayakshin, 2010), it is difficult to predict when and how often disc fragmentation leads to the formation of planets (Kuiper, 1951), brown dwarfs (Stamatellos & Whitworth, 2008, 2009) or secondary stellar companions (Kratter et al., 2010).

The scale of uncertainty in this problem is immense and affects our understanding of even the most basic questions, especially in the theory of planet formation. Direct imaging surveys show that the occurrence rate of wide separation (tens of AU or more) planetary mass companions to FGK stars, and also brown dwarfs, is just a few % (Biller et al., 2013; Chauvin et al., 2015; Reggiani et al., 2016; Vigan et al., 2017). This is much smaller than $\gtrsim 50$% observed planet occurrence rate at separations less than a fraction of AU from the star (see chapter 2 in Winn & Fabrycky, 2015). One interpretation of this result is that gravitational disc instability rarely makes planetary-mass objects (Kratter et al., 2010; Forgan & Rice, 2013a; Rice et al., 2015; Vigan et al., 2017). On the other hand, if radial migration and tidal disruption transmogrify planetary mass gas clumps into short period planets, including sub-Neptune mass planets (Boley et al., 2010; Nayakshin & Fletcher, 2015), then the rate at which GI fragmentation forms planetary-mass clumps could be much higher; the resulting planets are simply not where they were born.

Furthermore, there is now observational support that at least some initially widely separated objects end up at sub-AU separations from the star, presumably due to disc migration. The frequency of appearance of planets more massive than $\sim 4$ Jupiter masses and brown dwarf companions to stars do not correlate with the host star metallicity (Raghavan et al., 2010; Troup et al., 2016; Nayakshin, 2017a; Santos et al., 2017), indicating that these objects probably did not form by Core Accretion (which predicts an opposite correlation, see Mordasini et al., 2012). Additionally, the properties and statistics of very strong episodic flaring of young protostars, known as FU Ori outbursts (Hartmann & Kenyon, 1996), are consistent with stars tidally disrupting and devouring (Vorobyov & Basu, 2006; Takami et al., 2018) up to a dozen gas clumps per lifetime.

There are many physical uncertainties in the physics of the problem, e.g., disc opacity (Meru & Bate, 2010), initial conditions for disc fragmentation (Vorobyov & Basu, 2010; Zhu et al., 2012a), treatment of gas cooling close to and inside the Hill sphere of the planet (Nayakshin & Cha, 2013; Stamatellos, 2015; Mercer & Stamatellos, 2017), dust growth and dynamics inside the clump, which may stronly affect clump cooling and heating balance (Helled & Bodenheimer, 2011; Nayakshin, 2016), etc.

However, in addition to this, different simulation codes use different numerical algorithms to model the same processes, and it is not clear if applying these codes to the same problem will yield identical results. The goal of our paper is evaluate how the simulation results differ between some commonly used numerical codes. To focus on this issue alone, we set up a physically simple test problem of a gas giant planet embedded in a massive gas disc at an initial separation of 120 AU. The disc cooling is treated with the widely used idealised $\beta$-cooling prescription (Gammie, 2001; Rice et al., 2005).

To disentangle various effects, we perform four comparison runs. The initial planet mass is set to $M_{\rm p0}=2{\,{\rm M}_{\rm J}}$ in three of the runs and to $M_{\rm p0}=12{\,{\rm M}_{\rm J}}$ in the fourth. As explained above, gas accretion onto the gas clumps is an integral part of the problem. Therefore, in two of the $M_{\rm p0}=2{\,{\rm M}_{\rm J}}$ runs we turn off gas accretion onto the planet, setting instead a relatively large gravitational softening length parameter to reduce the amount of gas flowing into the gravitational potential well of the planet (as was also done by Baruteau et al., 2011). In the other two comparison runs, a sink particle prescription is used to absorb the gas accumulating at the planet location.

The paper is structured as follows. In §2 we describe the physical setup and initial conditions of the problem, and describe the contributing codes. In §3 we present main results of our paper. A comparison of the results to population synthesis prescriptions is made in §4, and in §5 we discuss observational implications of this work.

2 Problem and Numerical detail

2.1 Contributing Codes

There are five 3D SPH codes that we compare here: PHANTOM (Price et al., 2017), GADGET (Springel, 2005), SPHINX (Dehnen & Aly, 2012), SEREN (Hubber et al., 2011a; Hubber et al., 2011b), and SPHNG (Benz, 1990). The Meshless Finite Mass code GIZMO (Hopkins, 2015) builds on SPH methods and adds a kernel discretization of the volume, coupled to a high-order matrix gradient estimator. The GIZMO-MFM numerical scheme has a higher order consistency and appears to overcome some of the numerical viscosity issues in SPH, and has been recently shown to reproduce the expected convergence of the critical cooling timescale for fragmentation (see Deng et al., 2017), which has been hard to achieve with SPH methods previously (e.g., Meru & Bate, 2010). Finally, FARGO is a 2D fixed cylindrical grid finite differencing code (Masset, 2000) which has been widely used for studies of planet migration and has shown consistency with analytical solutions in the linear regime applicable to much lower mass planets (e.g., Baruteau & Masset, 2008) than studied here.

2.2 Problem choice

The potential formation of gas giant planets via gravitational instability of protoplanetary discs (e.g., Kratter & Lodato, 2016) motivates our study. To this end, all of our runs use a massive gas disc with initial mass $M_{\rm init}=0.2M_{\odot}$ as an initial condition for all of our runs. The disc is in circular rotation around a star with mass $M_{*}=1M_{\odot}$. At fragmentation, the disc Toomre (1964) $Q$-parameter is $Q\la 1$ (e.g., Boley et al., 2010). Such discs generate spiral density arms. Interactions of the planets with the arms give stochastic velocity kicks to the planets (e.g., Baruteau et al., 2011). In addition, fragmenting discs usually hatch more than one gas clump. Clump-clump interactions also lead to angular momentum exchange between the clumps (Cha & Nayakshin, 2011) and even mergers (Hall et al., 2017). These processes are stochastic and make numerical simulations of planet migration with different codes susceptible to small numerical detail.

To avoid this stochasticity, we simplified the task at hand by choosing the parameters of the problem such that the Toomre parameter of the disc is slightly larger than expected at fragmentation, i.e., $Q\gtrsim 2$ everywhere, which makes the disc gravitationally stable. We then inject a planet into the disc and follow its evolution numerically. It is clearly desirable to extend the code comparison in the future in the regime in which the disc is free to fragment and form more clumps.

An ideal gas equation of state is used in this paper with the adiabatic index $\gamma=7/5$, as appropriate for diatomic gas. The star irradiates the disc and sets the minimum irradiation temperature, which is a function of radius $R$:

[TABLE]

where $T_{0}=20$ K and $R_{0}=100$ au. The irradiation temperature corresponds to the specific internal energy,

[TABLE]

where $\mu=2.45m_{p}$ is the mean molecular weight of the gas.

The radiative cooling of the disc is modelled with the $\beta$-cooling prescription widely used in the literature to model marginally stable self-gravitating discs (e.g., Rice et al., 2005). The irradiation from the central star is additionally present as a heating term, so that the specific internal energy of the gas, $u$, evolves according to

[TABLE]

where $t_{\rm cool}=\beta\Omega_{\rm K}(R)^{-1}$,

[TABLE]

We use $\beta=10$ for the runs presented below. This value of $\beta$ is comparable to the critical fragmentation $\beta$ for $\gamma=7/5$ as found by Rice et al. (2005), although recent simulations with GIZMO-MFM suggest that disc fragmentation may occur at lower $\beta$ (e.g., Deng et al., 2017) for this code. However, the inclusion of external irradiation will also likely lead to fragmentation happening for lower values of $\beta$ (Rice et al., 2011)

2.3 Initial conditions

We first describe the initial conditions for the SPH codes and GIZMO-MFM. The star is treated as a sink particle that accretes any SPH particles that enter inside the sink radius, $R_{\rm sink}=3$ au. The gravitational softening of the star is set at $h_{\rm g}=0.25$ au. The disc is initially set up with the surface density profile

[TABLE]

where $R_{\rm in}=10$ au and $R_{\rm out}=300$ au are the inner and the outer initial disc radii, respectively. The disc is relaxed for about 10 orbits at the outer edge before the planet is inserted. This is done to allow the disc to settle into a vertical hydrostatic balance and to damp out radial disc oscillations. During the disc relaxation procedure, a small fraction ($\sim 3$%) of the SPH particles are accreted onto the central star. This is inevitable due to artificial viscosity of the disc increasing in regions of lower particle number, which is usually near the inner disc boundary.

These initial conditions, after the relaxation procedure was applied, are presented in fig. 1. The top panel shows the gas column density multiplied by radius, e.g., $\Sigma(R)\times(R/100$ AU) and the vertically averaged gas temperature profile $T(R)$. Both of these are compared to the respective column density and temperature profiles before the relaxation (blue dashed curves). We see that both the inner and the outer regions of the disc are depleted by the relaxation process, but that the region between $R\sim 30$ AU and $R\sim 200$ AU has a smooth $\Sigma\propto 1/R$ profile. The gas temperature profile is very close to equation 1, except for radii $R\lesssim 20$ AU where the artificial viscosity heating is not negligible. Since our relaxed disc has a strong roll-over at radii smaller than $R\sim 30$ AU, we expect that the planet migration process in this numerical setting will be strongly affected at radii of about $40$ AU. The bottom panel of Fig. 1 presents the disc aspect ratio $H/R$ normalised to 0.1 and the Toomre $Q$-parameter. As stated in the Introduction, the disc is everywhere stable to self-gravity and does not fragment.

For FARGO, the initial conditions were obtained in the same physical setup but the code was relaxed for 50 orbits at the outer edge.

In the simulations presented below, time is counted from the relaxed initial condition shown in Fig. 1. We inject the planet instantaneously at $t=0$ on a prograde circular orbit centred on the star at the initial separation of $R=120$ AU. No change is made to the initial velocity of either the gas or the star. Note that for $M_{\rm p0}=2{\,{\rm M}_{\rm J}}$, the planet mass is only 1% of the disc mass and just 0.2% of the total mass of the system, so this approach is justifiable. While for the $M_{\rm p0}=12{\,{\rm M}_{\rm J}}$ simulation the error is larger, we prefer this approach because keeping the planet orbit fixed while increasing its mass slowly (a common approach in studies of low mass planet migration) would lead to undesirable modifications of the disc structure for our problem. For example, as found by Malik et al. (2015), the gap opening criterion should include a gap-opening time scale. If the planet migrates across the gap sooner than the gap could be excavated, no gap is opened. However, keeping the planet on a fixed orbit implies an infinite source of angular momentum and therefore may result in the planet opening a gap in the disc where none should be present.

2.4 Approach to code algorithm differences

Numerical hydrodynamics codes, whether particle or grid based, employ different numerical algorithms to integrate equations of motions, various time-stepping criteria, and approximate techniques to resolve contact discontinuities such as shocks and singularities arising in the gravitational potential and forces near point masses (Bodenheimer et al., 2007). For example, by default GADGET uses the Monaghan (1997) form of the artificial viscosity with the Balsara-Monaghan switch to reduce artificial viscosity in shear flows (Monaghan, 1992; Balsara, 1995), and the spline kernel for SPH (for details see Springel, 2005). More modern formulations of artificial viscosity exist and different SPH kernels are adopted by some of the other codes (see §3.4). It is possible to modify GADGET to use the same approaches. However, it is not possible in practice to modify all of the codes to employ exactly the same numerical algorithms due to significantly different intrinsic code designs. Additionally, such code alterations would defeat the purpose of our code comparison project as the codes actually being compared would then be different from their current community-used versions.

Therefore, we attempted no code modification in this project with only a few exceptions that relate to the most salient physics of the problem. For each test problem presented below, all of the codes use the same gravitational softening parameters and the accretion radii for the two sinks in the problem, as detailed further below. The sections below discuss the implementation of sink particle accretion, gravitational softening and artificial viscosity in the codes used in this paper.

2.4.1 GADGET

Our implementation of GADGET is very similar to the code description given in the instrument paper by Springel (2005), with a few changes detailed below. GADGET uses the spline kernel (Monaghan & Lattanzio, 1985) for both the SPH density field and computing the gravitational softening around all particles, including the sink particles. We use 40 particles for the neighbour search. The artificial viscosity of SPH is that given by the Monaghan-Balsara formulation (Gingold & Monaghan, 1982; Balsara, 1995), modified by the viscosity limiter prescription (see eq. 11 in Springel, 2005) to alleviate unwanted angular momentum transport in the presence of shear flows. We follow the default GADGET settings in this paper, keeping the artificial viscosity coefficient $\alpha_{\rm v}$ set to 1 for all times, and $\beta_{\rm v}=2\alpha_{\rm v}$.

The sink particles are implemented in a very simple way. Any SPH particle that is separated from the sink by a distance smaller than the accretion radius $R_{\rm a}$ is accreted by the sink. The linear momentum and mass of the accreted particle are added to that of the sink. Some authors consider more complicated gas accretion criteria. For example, Bate et al. (1995a) consider the expected pressure of the gas within the sink region and the binding energy of the gas with respect to the sink. However, there is much physical uncertainty in picking these additional gas accretion criteria. The sink radius defines the region of space where we have insufficient information (usually, no information at all) about the gas properties. The interactions of that missing gas with the SPH particle in question could change the properties of the latter in ways that cannot be computed. For example, an SPH particle on a hyperbolic trajectory around the sink is formally not bound to the sink and thus would not be accreted if one accretes only particles with negative binding energies (Bate et al., 1995a). However, the same particle may be accreted if the particle were to interact with the missing gas within the sink radius, shock due to this interaction, and then lose the excess energy through radiation.

For further discussion of these issues and tests of our GADGET implementation of the sink particle prescription, see Cuadra et al. (2006) and Humphries & Nayakshin (2018). Nayakshin (2017b) found that the sink radius prescription tends to over-estimate the gas accretion rate onto a planet embedded in a massive gas disc for simulation parameters comparable to those used here (see Fig. A1 in Nayakshin, 2017b). Gas accretion rates measured in this paper should be thus taken as upper limits to the corresponding astrophysical problem.

2.4.2 PHANTOM

Cullen & Dehnen (2010) introduced an artificial viscosity switch which utilizes the derivative of the velocity divergence to detect shocks. Due to the switch, the artificial viscosity coefficient $\alpha_{\rm v}$ is varied between a minimum value, $\alpha_{\rm min}$, far from the shock, and the maximum, $\alpha_{\rm max}=1$, reached close to the shock. We use this method for PHANTOM in this paper, as described in detail in §§2.2.7-2.2.9 in Price et al. (2017). We fix the artificial viscosity coefficient $\beta_{\rm v}$ at 4 for our comparisons runs (see Price & Federrath, 2010). An exception to this is §3.4 where we explore how results depend on the choices of the artificial viscosity prescription for PHANTOM.

Gravitational softening in PHANTOM is different for interactions between sinks and interactions between sinks and SPH particles. The sink-sink softening is set to 0 by default. The sink-gas gravitational softening length is the maximum between the fixed softening length of the sink and the gas particles adaptive softening length. Gravity for SPH particles is softened by the SPH kernel function (see §2.12.2 in Price et al., 2017).

Compared to GADGET, the PHANTOM default sink particle implementation also sets constraints on the binding energy and relative angular momentum of the SPH particle to be accreted. In this paper we disable these additional checks and use the same approach as specified in §2.4.1.

2.4.3 SPHINX

SPHINX is an SPH code based on a conservative formulation (as derived from a variational principle, e.g. Price 2012) with individual artificial dissipation strengths $\alpha_{\rm v}$ adapted using the Cullen & Dehnen (2010) switch with $\beta_{\rm v}=2\alpha_{\rm v}$. The details of the artificial viscous force differ slightly (by an amount $O(h^{2})$) from traditional implementations to accommodate the one-sweep SPH algorithm, which avoids separate sweeps over all particle neighbours for the density and force computations. For the runs here, we use the Wendland (1995) $C^{2}$ smoothing kernel, which scales as $w\propto(h-r)^{3}(h+3r)$ for $r<h$ with smoothing length $h$, adjusted to obtain $N_{h}=4\pi\rho h^{3}/3m=80$ at each time step. Gravity is computed using a $C^{\infty}$ softening kernel with density $\propto(r^{2}+h_{\rm s}^{2})^{-7/2}$, which results in a smaller force bias than traditional Plummer softening (Dehnen, 2001). Individual softening lengths $\epsilon$ are scaled to the smoothing lengths $h$ such that the estimates for the gas and gravitating mass densities are mutually consistent (have the same bias). SPHINX uses an oct-tree for neighbour search (and gas-selfgravity which is computed using the fast multipole method Dehnen 2000) and the leap-frog (2nd order symplectic) time integrator. Star and planets are represented by sink particles, whose gravity is computed by direct summation. Any gas particle within one sink radius is accreted by a sink particle, whereby its mass, linear and angular momentum, as well as energy is absorbed by the sink particle (which carries a spin and internal energy for this book keeping).

2.4.4 SPHNG

SPHNG is based on the version developed by (Benz, 1990) and first presented by (Bate et al., 1995a). It uses variable individual smoothing lengths $h_{j}$ and adjusts them so that the number of nearest neighbours for any particle is $50\pm 20$. It also uses individual particle time-steps to simulate dense regions with sufficient precision while avoiding over-simulation of less dense regions, and integrates the particles using a second order Runge-Kutta scheme. The standard artificial viscosity (Monaghan, 1992), with $\alpha_{\rm v}=1.0$ and $\beta_{\rm v}=2.0$, and standard spline kernel are used. A binary tree is used to calculate neighbour lists and to determine gravitational forces between gas particles, with the gravitational force softened by the SPH kernel function (Price et al., 2017). The gravitational force between the gas particles and the sink particles is, however, done using a direct calculation, which is softened by replacing the $1/r^{2}$ gravitational force dependence with $1/(r^{2}+h_{\rm s}^{2})$. If accretion onto the sink particles is allowed, then particles are only accreted if they are bound and if the specific angular momentum of the particle is less than that required for them to form a circular orbit at the accretion radius (Bate et al., 1995a).

2.4.5 GIZMO

The GIZMO code is a multi-method code which inherits the tree-based gravity algorithm from GADGET3 (see Springel, 2005, for GADGET2 code description) and couples it with different Lagrangian hydrodynamical solvers. For this paper we employ the Meshless Finite Mass (MFM) hydro method in GIZMO which solves the inviscid fluid equations by partitioning the computational domain using volume elements associated with a particle distribution, and computing fluxes through the volume ‘overlap’ by means of a Riemann solver as in finite volume Godunov-type methods (Hopkins, 2015). Volume elements are constructed via convolution integrals with kernel functions analogous to those adopted in SPH. Owing to the use of a Riemann solver (here we use the HLLC solver and the minmod slope limiter), GIZMO-MFM employs no explicit artificial viscosity. This numerical method appears significantly less dissipative than SPH for differentially rotating flows, better conserving angular momentum and vorticity (Hopkins, 2015; Deng et al., 2017). The kernel for the volume partitioning, the gravitational softening and the sink particle implementation are all identical to those of GADGET used in this paper (§2.4.1).

2.4.6 SEREN

The SPH code seren was developed for star and planet formation simulations by Hubber et al. (2011b); Hubber et al. (2011a). The code uses an octal tree to compute gravity and find neighbours, multiple particle timesteps, and a 2nd order Runge-Kutta integration scheme. To simulate the effect of physical viscosity in discs, seren uses a time-dependent artificial viscosity (Morris & Monaghan, 1997) with parameters $\alpha_{\min}=0.1$, $\alpha_{\max}=1$ and $\beta_{\rm v}=2\alpha_{\rm v}$, so as to reduce artificial shear viscosity away from shocks (this scheme is the predecessor of the Cullen & Dehnen, 2010, method). Sink particles, which interact with the rest of the computational domain only through their gravity, are used to represent the central star and the planet (Bate et al., 1995b). Gas particles accrete onto a sink when they are within the sink radius and bound to the sink (see Hubber et al., 2011a). Once gas particles are accreted, their mass and linear angular momentum is added to sink. The gravitational force between gas particles and a sink is found through a direct calculation and softened according to $1/(r^{2}+h_{\rm s}^{2})$ to avoid unphysically large gravity forces.

2.4.7 FARGO

FARGO is a 2D grid based, staggered-mesh code (Masset, 2000; Baruteau & Masset, 2008) that has been used extensively to study planet migration (Masset, 2002; Masset & Casoli, 2010; Baruteau et al., 2011). For the runs presented here, we use a cylindrical grid with 508 and 1536 cells in the radial and azimuthal directions, respectively. The radial grid is logarithmic with the inner and outer boundary conditions set at 10 and 300 AU, respectively. Von Neumann–Richtmyer artificial bulk viscosity is used to treat contact discontinuities (Stone & Norman, 1992).

For this paper, FARGO also uses a fixed gravitational softening parameter $h_{s}$ as for all the other codes, which is a break with the common practice of scaling $h_{s}$ with the local disc scaleheight or the star-planet separation (e.g., Baruteau et al., 2011), but allows for a more uniform comparison between the codes. Specifically, the softening parameter used in grid simulations is typically set to $h_{s}=\epsilon H$ where $\epsilon\sim O(0.1)$ (Müller et al., 2012). In this case the gravitational softening would be a function of position as $H\propto R$ for our simulations (see fig. 1). The consequences of this for numerics are not immediately obvious, but we note that for $\epsilon=0.1$ and $H\sim 0.1R$, the adaptive softening is equivalent to $h_{\rm s}=0.4-1.2$ AU in the radial range 40-120 AU, which is not too dissimilar from the 1 AU and 2 AU fixed smoothing employed in Runs 1 and 2 (see below). For a relatively large fixed value $\epsilon=0.7$ we find that the FARGO migration timescales increase by $\approx 50\%$ compared to those presented in this paper.

2.5 The comparison runs

It is possible to resolve the pre-collapse gas giant planets (clumps) in modern computer simulations directly (e.g., Boley et al., 2010; Galvagni et al., 2012; Zhu et al., 2012b; Nayakshin, 2017b; Hall et al., 2017). However, while the clumps can be resolved and modelled from the point of view of hydrodynamics, other physics, e.g., a proper equation of state including molecular hydrogen internal degrees of freedom, dust dynamics and radiative transfer, are not yet implemented in most of the codes available to us here. Any simplified radiative transfer scheme applied to the clumps would necessarily over-simplify their internal physics (their cooling balance is significantly different from that of the disc; e.g., see Vazan & Helled, 2012) and would thus be riddled with its own uncertainties. A more prudent approach for us to follow here is to model the planet as a sink particle, just as the star, albeit with its own gas accretion (sink) radius.

Table 1 shows the parameters that distinguish the four different comparison runs that are presented below. In Runs 1-3, the initial planet mass is set to $M_{\rm p0}=2{\,{\rm M}_{\rm J}}$, whereas Run 4 starts with $M_{\rm p0}=12{\,{\rm M}_{\rm J}}$. In Run 1 and 2, gas accretion onto the sink is completely turned off by setting the accretion radius to zero. This is done to try to isolate the effects of planet migration versus gas accretion onto the planet. This is especially important since FARGO is a grid based code in which implementation of gas accretion is drastically different from the sink particle method of SPH codes. Therefore, Runs 1 and 2 can be simulated with SPH codes and FARGO, whereas Runs 3 and 4 are done with SPH only.

Turning off gas accretion onto a planet does not come free of numerical cost. Gas that gets bound to the planet may eventually get very close to the planet. A very high gas density around the planet is numerically challenging as the SPH particle time step becomes too short for the code to execute effectively. Therefore, to avoid that, in Runs 1 and 2 the planet softening radius, $h_{\rm s}$ is increased to 1 and 2 AU, respectively, from the much smaller value used in Run 3. For a similar reason Run 4 uses a larger accretion radius than Run 3.

The initial SPH particle number is $N=10^{6}$ for all of the runs presented here.

2.6 Analytical expectations

Tanaka et al. (2002) derived an analytical expression for type I migration of a low mass planet in an isothermal disc. The migration timescale, defined as

[TABLE]

where $\dot{R}$ is the rate of change of planet-star separation due to gravitational torques from the disc, is given by

[TABLE]

Here $\lambda$ is the exponent of the surface density power law, $\Sigma\propto R^{-\lambda}$, $\Sigma_{\rm p}$ is the surface density at the planet location, $M_{\star}$ and $M_{p}$ are the star and planet masses, respectively, $r_{p}$ is the planet-star separation, $c_{s}$ is the gas sound speed at the planet and $\Omega_{p}$ is the planet Keplerian angular velocity. For the initial parameters of our disc and $M_{\rm p}=2{\,{\rm M}_{\rm J}}$, we obtain a migration time scale of $\tau_{\rm iso}=14.6\times 10^{3}$ yr. Even though our discs are not isothermal, the results of Tanaka et al. (2002) are widely used, and serve as a useful comparison for us.

Baruteau et al. (2011) used the 2D code FARGO to study planet migration in very massive self-gravitating discs, for which the Toomre parameter $Q$ self-regulates to a value between $\sim 1.5$ and $\sim 3$ over a broad range of radii. These authors also offered an analytical expression for the migration time scale:

[TABLE]

where $q=M_{\rm p}/M_{\star}$ is the mass ratio; $Q_{p}$ is the Toomre parameter and $h_{p}=H/R$ at the planet position. For the initial parameters of our Runs 1-3, eq. 8 yields $\tau_{\rm sg}=5.0\times 10^{3}$ yr at a separation $R=120$ AU.

3 Results

3.1 At a glance

Fig. 2 shows the planet separation against time for Runs 1 & 2. To recap, gas accretion onto the planet is off, and instead a relatively large gravitational softening parameter is used. Despite this, some gas accumulates deep inside the Hill sphere, and differently so for different codes. This appears to be the primary reason why Runs 1 and 2 stalled for SEREN at around 5000 yr. Fig. 3 shows the results of Runs 3 & 4 (left and right panels, respectively) in which gas accretion onto the planet (sink particle) is allowed. The sink mass versus time is shown in the lower panels.

A cursory look at figs. 2 & 3 shows that there is a general qualitative agreement between the different codes. For example, in Runs 1-3 the planet manages to migrate to separations of $40-60$ AU for most of the codes, whereas in Run 4, in which the planet is much more massive, the planet stalls further out due to it opening a deep gap in the disc. At the same time, there are significant quantitative disagreements between the codes. All of the codes show that the planet develops orbital eccentricity, but the actual value of the eccentricity is different, varying between $\sim 0.01$ to the maximum of $\sim 0.1$.

3.2 Analysis of Runs 1-3

3.2.1 Migration rates

We now analyse Runs 1-3 in which the planet initial mass is $M_{\rm p0}=2{\,{\rm M}_{\rm J}}$. To aid quantitative analysis, we determine migration time scale, $\tau$, from the simulations. A straight-forward use of eq. 6 to calculate $\tau$ from the simulation data is ill advised due to planets having non-zero eccentricity: the instantaneously defined migration time varies significantly over a fraction of the planet orbital timescale. Some sort of time averaging of $\tau$ over times at least as long as an orbital period is thus needed.

To do so, we first define a time-dependent migration rate as the final difference $\Delta R/\Delta t$, where the separation and time differences are counted from the initial values:

[TABLE]

where $R_{0}=R(t=0)=120$ AU, $t>0$ is time, and $R(t)$ is the planet-star separation at that time. To remedy the oscillatory behaviour in $\dot{R}$ due to finite orbital eccentricity, we define an orbit-averaged quantity

[TABLE]

where $T_{\rm p}$ is the planet orbital period at location $R(t)$. We use this definition to define the planet migration rate after $t=4,000$ yr for all of the codes, which we label $\tau_{4}$. We then also define the migration time scale $\tau_{7}$, following the procedure outlined above, but uding the data between 4,000 and 7,000 yrs. Comparison of $\tau_{4}$ and $\tau_{7}$ tells us how the migration rate varies as the planet gets closer to the star. Due to a non zero planet orbital eccentricity a finer time-resolved analysis of the migration rate does not appear well justified.

Fig. 4 compares the migration time scales $\tau_{4}$ and $\tau_{7}$ (left and right panels, respectively) for all the codes for Runs 1-3, which are shown with the coloured symbols. The dashed and solid horizontal lines show the migration timescales given by eq. 7 and 8, respectively. These analytical estimates of $\tau$ are computed using the initial disc properties (see fig. 1).

Taking the full range of $\tau_{4}$ and $\tau_{7}$ values, we see that they vary by a factor of $3-4$ between the different codes for Runs 1 & 2, and by a smaller factor of $\sim 2$ for Run 3. For $\tau_{4}$, the mean of the migration time scales are closer to the Tanaka et al. (2002) expression, but for $\tau_{7}$ the mean lies between the analytical estimates of Tanaka et al. (2002) and Baruteau et al. (2011). The range in the migration time scales is similar to the factor of $\sim 3$ difference between these two analytical results. We also note that $\tau_{4}$ is longer than $\tau_{7}$ for most of the runs, implying that migration of the planet accelerates somewhat as the planet gets closer to the star (as long as it remains in the Type I). The same trend is predicted by the formulae shown in equations 7 and 8. We conclude from fig. 4 that there is a qualitative agreement not only between the different codes but also with the theory.

Comparing Runs 1 and 2, we note that the migration timescales vary by $\sim 10$% for most of the codes whereas $h_{\rm s}$ changes by a factor of two. However, for SEREN the difference the two runs is larger, and is in the opposite sense compared with most of the other codes. This is likely due to a non linear interplay of how gravitational softening affects gravitational torques vs planet accretion. To make further progress we must consider the role of planet accretion in greater detail, as the evolving planet mass certainly affects the migration rate.

3.2.2 Hill mass versus sink mass

The sink particle mass may not always properly reflect the mass of the planet. To quantify this, we define an effective Hill mass of the planet, $M_{\rm H}$, as the sink mass plus the mass of the gas within $R_{\rm H}/2$ of the planet. The choice of $R_{\rm H}/2$ is motivated by results of Nayakshin (2017b), who finds that gas bound to the planet is usually located within half the Hill radius; material between $R_{\rm H}/2$ and $R_{\rm H}$ is much more likely to be lost as the planet migrates inwards.

We should also note that the Hill radius definition needs to include the mass of the gas envelope around the sink itself, that is,

[TABLE]

where we use $M_{\rm H}$ rather than the sink mass, $M_{\rm p}$. When the Hill mass is dominated by the sink mass, $M_{\rm H}$ can be safely replaced by $M_{\rm p}$, and the calculation of $M_{\rm H}$ from the particle data is trivial. In general, however, the mass of the gas surrounding the sink is not negligible, so we iterate over $R_{\rm H}$ and $M_{\rm H}$ to find self-consistent values for these two quantities that obey eq. 11.

Fig. 5 shows the Hill mass and the sink mass for Runs 1-3 calculated for the different codes. For Run 3, where gas accretion onto the sink is allowed, we see that for all the codes $M_{\rm H}\approx M_{\rm p}$. In other words, the gas mass within the Hill sphere is negligible compared with the sink mass. As the sink mass grows rapidly by gas accretion, this also means that once gas enters the Hill sphere it accretes onto the sink rapidly, so there is never a dynamically significant gas envelope around the sink. This is expected since we use a relatively large value of $R_{\rm a}=0.5$ AU for Run 3. Nayakshin (2017b) found that the accretion rate onto the sink is roughly proportional to the sink radius (see Appendix in that paper) and that sink radii larger than $\sim 0.1$ AU over-estimate the rate of gas accretion onto the sink when compared with a simulation in which the clump was directly resolved111However, it is not clear what is the appropriate value of $R_{\rm a}$ to use in general as it also depends on the numerical resolution, e.g., the number of SPH particles used. Using too low a value of $R_{\rm a}$ may lead to an under-estimate of the accretion rate as the sink region may become unresolved due to a finite SPH particle resolution..

Fig. 5 shows that in Runs 1 & 2 the mass of gas surrounding the sink particle within $R_{\rm H}/2$ is comparable to the sink mass by the end of the runs, in stark contrast to Run 3. For PHANTOM in particular, at $t=10,000$ yr, the Hill mass is dominated by the envelope.

In a qualitative agreement between the codes, $M_{\rm H}$ is always larger in Run 3 than in Runs 1 and 2. This demonstrates that the gas envelope around the planet particle, which builds up in Runs 1 and 2 but not in Run 3, has a detrimental effect on further gas accretion onto the planet. This is likely due to the extra pressure of the envelope, which makes it more difficult for the gas entering the Hill sphere to remain there. However, the exact trend going from Run 1 to Run 2 in the Hill mass is not the same for the different codes. While for GADGET and GIZMO-MFM a larger gravitational softening results in a lower mass gas envelope, this is not the case for PHANTOM and SPHINX. Therefore, gas accretion onto the planet (or the planet envelope) remains a significant source of uncertainty even in the simulations where gas accretion is turned off. An exception to this could be problems where gas accretion onto the planet is physically insignificant, such as when the planet mass is very sub-Jovian or the gas cooling time is very long (as in the $\beta\gg 1$ regime in Nayakshin, 2017b).

Let us now compare the uncertainties in the planet accretion rate versus that in migration. The left panel of fig. 3 shows that there is more disagreement in the planet mass versus time plot between the different codes for Run 3 than in the planet migration tracks. The mass of gas accreted by the planet varies from a minimum of $\sim 4{\,{\rm M}_{\rm J}}$ to a maximum of $\sim 12{\,{\rm M}_{\rm J}}$, whereas the planet migration timescales vary by less than a factor of 2. We believe that this smaller disagreement in planet migration rates may be somewhat fortuitous. As the planet mass increases, the analytic formulae in the linear type I regime (e.g., eq. 7) predict that the migration rate should increase linearly with planet mass. However, as the planet starts to open a gap, it starts to transition into a slower type II regime. The migration rate therefore depends on the planet mass somewhat less strongly than can be expected based on the theoretical type I predictions.

3.3 Run 3 and Run 4

3.3.1 Gap opening

Runs 3 and 4 both use the sink particle prescription but differ in the initial sink mass, $2{\,{\rm M}_{\rm J}}$ and $12{\,{\rm M}_{\rm J}}$, respectively. These two simulations cover the parameter space in which a growing planet goes from migrating in type I (no gap in the disc) to type II (a deep gap opened). In the outer massive disc, both planet migration rates and gas accretion rates onto the planet are far larger in the Type I regime than in the Type II regime (e.g., Zhu et al., 2012b; Nayakshin, 2017b). The time and radial location where the switch between migration regimes occurs is thus of a significant importance.

Fig. 6 shows with different coloured lines the planet mass versus separation tracks for Runs 3 (left panel) and 4 (right panel) for all the eligible codes. The planets start at the lower right corner and move towards the upper left corner in this diagram.

There are also four black curves in the figure that show theoretical predictions from Crida et al. (2006) for when a deep gap in the disc should be opened. According to these predictions, the planet opens a gap when the parameter $C_{\rm p}$ is smaller than unity:

[TABLE]

Here $\alpha$ is the physical viscosity parameter of the gas disc (Shakura & Sunyaev, 1973). We do not set a physical viscosity parameter in the runs presented here (PHANTOM offers a facility for this but most other SPH codes do not). However, artificial viscosity in numerical schemes can mimic certain effects of a physical viscosity. Price et al. (2017) show that for the PHANTOM viscosity implementation, artificial viscosity parameter $\alpha_{\rm v}$, set to unity for all SPH codes here (but see §3.4), results in effective Shakura & Sunyaev (1973) viscosity parameter

[TABLE]

where $h_{\rm sml}$ is the SPH smoothing lengh and $H$ is the local disc vertical height scale (see Murray, 1996). At the separation where our planets open gaps, we have $h_{\rm sml}/H\approx 0.4$, and hence the effective disc viscosity of these codes is about $\alpha=0.03$.

Additionally, self-gravitating protoplanetary discs generate physical viscosity that saturates at a maximum value of $\alpha\sim 0.06$ (Gammie, 2001; Rice et al., 2005) for marginally stable discs. The value of the $Q$-parameter for our disc is significantly greater than the critical $\sim 1.5$ and we thus expect that the effective $\alpha$ from the disc self-gravity is much smaller than the maximum value.

Fig. 6 show the planet gap-opening mass as a function of separation for our initial discs, defined as the planet mass for which $C_{\rm p}=1$. The solid curve sets $\alpha=0.03$, whereas for the dashed and the dotted curves $\alpha=0.05$ and $\alpha=0.1$, respectively. Since planet migration effectively stalls (at least on the time scales of our simulations) when the planet switches to the type II migration regime, the radial location of this switch can be identified in the figure as the point where the planet track turns from being mainly horizontal to being more vertical. For Run 3, the left panel of fig. 6 shows that the location at which the migration type switches is approximately consistent with the Crida et al. (2006) prediction for $\alpha=0.1$, although the actual value of the separation and planet mass at that point are somewhat different for the codes. However, the estimated effective disc viscosity for the codes is $\alpha=0.03$, and the respective (solid) curve in fig. 6 yields significantly smaller masses. The only exception to that is GIZMO-MFM whose meshless finite mass scheme was shown to provide smaller artificial viscosity (Deng et al., 2017).

The results of Run 4 are largely consistent with this picture. We see that the gap opening value of planet mass and separation lie close to the $\alpha=0.1$ theoretical curve, with GIZMO-MFM transiting into type II migration somewhat earlier once again. One exception to this is PHANTOM, for which the planet seems to cross the migration type dividing line rather uneventfully.

The fact that our simulated gas clumps open gaps at higher masses and later in time than predicted by the Crida et al. (2006) analysis confirms the findings of Malik et al. (2015) who showed that in massive circumstellar discs, gap opening is more difficult than for less massive discs. As shown by Malik et al. (2015), if planets migrates through the horse-shoe region faster than the gap can be excavated by planet toques, the gap remains closed even if $C_{\rm p}$ falls below unity.

Finally, although our code migration comparison project is not designed to study the longer term planet evolution that occurs in the Type II regime, we can see from fig. 6 that there is a significant disagreement in the planet evolution once it crosses over into the Type II regime. While qualitatively we see that planets tend to stall in Type II, as expected, some codes predict that the planets continue to migrate in while others (PHANTOM in the left hand panel) start to migrate outward. This may indicate that the secular evolution of the planets in the Type II migration regime is even more model dependent than the Type II which we mainly aim to study here.

3.3.2 Gas accretion time scales

As emphasized by previous authors, there is a competition between the process of gas accretion onto the planet and its inward migration (e.g., Zhu et al., 2012b; Nayakshin, 2017b). This competition plays a significant role in shaping of the outcome of disc fragmentation. It is hence convenient to define, in addition to the migration time scale, an accretion time scale for the planet, $t_{\rm acc}$,

[TABLE]

where $\dot{M}$ is the gas accretion rate onto the planet. The corresponding dimensionless quantity $\tau_{\rm acc}$,

[TABLE]

where $T_{\rm p}=2\pi/\Omega_{\rm p}$ is the orbital period at the planet location, will be useful as well.

Bate et al. (2003) studied planet migration and accretion in isothermal discs and found that the following equation describes the gas accretion rate onto the planet well in the Type I migration regime,

[TABLE]

where $b\approx 2.3$ empirically and $\rho$ is the disc midplane density. By writing $\rho=\Sigma/(2H)$ and expressing

[TABLE]

where $Q$ is the Toomre parameter at the planet location, we can re-arrange the Bate et al. (2003) result as

[TABLE]

Zhu et al. (2012a) used a 2D code to study clump migration and accretion, and provided a 2D estimate for the rate of gas accretion onto the planet,

[TABLE]

Expressing $\Sigma$ through eq. 17 again, we obtain the corresponding gas accretion time scale

[TABLE]

Since for our planets $R_{\rm H}\sim H$ within a factor of two or so, eq. 20 is actually not very different from eq. 18.

Fig. 7 shows dimensionless accretion time scales for Runs 3 and 4. The black curves show the analytic estimates obtained with eqs. 18 and 20, respectively. For eq. 18, we show three curves which use $b=2.3$ (as in Bate et al., 2003), and then also $b=$ 1, and $1/3$. We can see that both analytic prescriptions predict much faster accretion rates onto the planet than actually measured in the simulations. This is most likely due to the analytic estimates assuming an isothermal equation of state and therefore the maximum efficiency for gas capture onto the planet. In the runs presented here, the gas is not isothermal and heats up due to adiabatic compression in the Hill sphere. The cooling rate $\beta$-parameter is $\beta=10$, which is relatively large. Nayakshin (2017b), see also Humphries & Nayakshin (2018), found that gas accretion onto planets is significantly suppressed for $\beta\gtrsim$ a few. The isothermal gas accretion rate estimates from Bate et al. (2003) and Zhu et al. (2012b) physically corresponds to the $\beta\ll 1$ regime investigated in Nayakshin (2017b), for which much higher accretion rates were indeed obtained. It appears that $b\approx 1/3$ in eq. 18 fits the gas accretion rates in the Type I migration regime best.

Fig. 7 also demonstrates that the accretion time increases strongly when the planet switches to the type II migration regime. This has also been seen in previous simulations (e.g., Bate et al., 2003) and is to be expected as the planet clears its immediate neighbourhood of gas, chocking its own growth.

The initial dips in the accretion time for both panels in fig. 7 are caused by our artificial initial conditions, in which a massive planet is injected in the disc. The gas within the Hill sphere of the planet then finds itself strongly bound to it and accretes onto the planet on a time scale shorter than the local dynamical time, $1/\Omega$. This initial transient is followed by a more self-consistent evolution in which the gas in the Hill sphere of the planet ”knows about its existence”.

3.4 Importance of artificial viscosity prescription

Artificial viscosity is used in SPH and grid based codes to treat flow discontinuities such as shocks (Monaghan, 1992; Bodenheimer et al., 2007). The codes we test here differ in their implementation of the artificial viscosity. Some part of the differences in the results of Runs 1-4 (discussed in §3) may be due to these numerical technique differences. Varying the viscosity prescriptions for all of the codes would make the presentation of this paper overly long. Instead we pick one code, PHANTOM, and investigate how different artificial viscosity choices affect the results for just Run 3.

All modern SPH codes employ artificial viscosity prescriptions that include a term linear in $\Delta v$, the velocity difference between two interacting SPH particles, and a term quadratic in $\Delta v$ (Springel, 2005; Price et al., 2017). That is, the first term enters artificial viscosity with a dimensionless coefficient $\alpha_{\rm v}$, and the second with coefficient $\beta_{\rm v}$. In some codes, e.g., GADGET, these coefficients are fixed whereas in others such as PHANTOM they are allowed to vary in time during simulations. Cullen & Dehnen (2010) in particular presented a method in which $\alpha_{\rm v}$ depends on the time derivative of the particle velocity divergence. The latter is used as a shock indicator and helps to eliminate artificial viscosity away from shocks, reducing unwanted numerical dissipation in dynamically quiet regions. Additionally, there are different suggestions on the appropriate values for the coefficient $\beta_{\rm v}$ to use, and in fact this may depend on the problem studied (Price et al., 2017).

Fig. 8 shows how the planet separation (top panel) and planet mass (bottom panel) are affected by the changes in the viscosity prescription for Run 3. The solid curves show Run 3 in which the $\alpha_{\rm v}$ parameter is time-dependent as in the method of Cullen & Dehnen (2010), and is allowed to vary between $0\leq\alpha_{\rm v}\leq 1$. The different colours in the solid curves indicate different values of the coefficient $\beta_{\rm v}$, which we varied in a broad range, from $\beta_{\rm v}=0.2$ to $\beta_{\rm v}=10$. The dashed curves in fig. 8 show simulations with the same range in $\beta_{\rm v}$ but which now use a fixed value for $\alpha_{\rm v}=1$.

First, without reference to the different artificial viscosity values in the figure, we note that the larger the planet mass, the more rapidly the planet migrates, at least until it opens a gap and switches to type II migration. Such a trend simply reflects the fact that more massive planets migrate more rapidly in the Type I regime (eq. 7).

Another trend obvious through all of the curves is that the higher artificial viscosity simulations tend to yield smaller gas accretion rates onto the planet. The least viscous run (red solid curve) shows the the largest gas accretion rate onto the sink and the most rapid migration. The most viscous run (green dashed curve) shows the slowest migration and the smallest gas accretion rate. The rest of the runs show a continuous transition between these two extremes.

This gas accretion trend with artificial viscosity is most likely due to the artificial viscosity heating of the gas inside the Hill radius. The larger the gas viscosity, the larger the dissipation rate within the Hill sphere, making the gas hotter. Such sensitivity of gas accretion rate onto the planet to heating within the Hill sphere was seen in the previous literature although for different reasons. Nayakshin & Cha (2013) and Stamatellos (2015) included planet radiative feedback on the surrounding gas, and found that when the feedback is present, it keeps the gas hotter in the planet’s Hill sphere, stifling gas accretion onto it. Nayakshin (2017b) found that slower radiative cooling rates within the Hill sphere, which also makes the gas hotter in that region, likewise leads to a reduction in the gas accretion rate.

In greater detail, we see that the runs with $\beta_{\rm v}=0.2$ and $\beta_{\rm v}=2$ are virtually indistinguishable, implying that the quadratic term in the artificial viscosity prescription is negligible for these small values of $\beta_{\rm v}$ for the given problem. Higher values of $\beta_{\rm v}$ however definitely affect the results. We also see that the fixed $\alpha_{\rm v}$ simulations lead to less massive and less rapidly migrating planets that tend to open a gap sooner.

The range of migration rates and planet masses in fig. 8 is large enough to conclude that although the artificial viscosity is not the only reason for differences in the results from the four runs explored in this paper, it is one of the major reasons for these differences. For example, GADGET’s planet separation versus time track for Run 3 is similar to the green dashed curve in fig. 8 for PHANTOM obtained with a fixed $\alpha_{\rm v}=1$, as used by GADGET. However, by default GADGET uses $\beta_{\rm v}=2\alpha_{\rm v}$, which is much smaller than $\beta_{\rm v}=10$ for the green dashed curve. Clearly, other code differences, both in viscosity implementation (GADGET uses the Balsara, 1995, switch; PHANTOM does not), and in how artificial softening and gas accretion onto the sink is implemented must be at play. A recent study by Stamatellos & Inutsuka (2018) found that the artificial viscosity coefficient $\alpha_{\rm v}$ can also drive differences in planet accretion/migration.

On the other hand, while PHANTOM simulations suggest a higher artificial viscosity might suppress accretion via spurious heating of the gas surrounding a sink particle, the trend shown by the GIZMO-MFM results in this paper suggest the role of numerical viscosity might be more complex. Indeed, as shown in Deng et al. (2017), the MFM method, which does not employ any artificial viscosity, at variance with all SPH methods, minimizes spurious transport of angular momentum inside self-gravitating disks and results in a lower accretion onto sink particles (see Appendix B in Deng et al., 2017). Indeed MFM solves the fluid equations via Riemann solver as in Godunov-type finite volume methods, which removes the need of an artificial viscosity term in the hydro equations (Hopkins, 2015).

Artificial viscosity implementations in SPH can induce enhanced angular momentum transport, and thus accretion, in non-shocking rotating flows inside fluid disks, owing to the contribution of the linear in $\Delta v$ term (even with correction terms such as the Balsara switch, e.g., Kaufmann et al., 2007). Spurious heating and artificial angular momentum transport are thus two different unwanted effects of artificial viscosity which affect accretion in opposite ways. Quantifying the interplay of these two effects warrants further investigation. Nevertheless, it is noteworthy that, in the GIZMO-MFM runs, the reduced accretion limits asymptotically the mass growth of the protoplanet to less than $10M_{J}$, namely within the gas giant planet regime.

4 Comparison to population synthesis

At the time of writing, there are three detailed population synthesis models that address the evolution of clumps formed by gravitational instability at distances of tens to 100 AU. Such population synthesis is a necessary step to correctly interpret the results of large observational surveys (e.g., Vigan et al., 2017) with respect to how often disc fragmentation might result in the formation of massive planets and/or brown dwarfs.

The population synthesis models differ in assumptions about the initial state of the disc and the clumps, disc dissipation, clump radiative cooling, dust dynamics and core formation, clump migration and accretion. It is of course not possible for us to examine these different approaches here. However, we can investigate a more limited but better defined question: how well would these models reproduce the evolution of the clumps that we see in our numerical models given the same disc and clump properties as our simulations?

To facilitate the population synthesis comparison to the simulations presented in this paper, we shall utilize the fact that the disc surface density profiles evolve relatively weakly in Run 3 as the planet remains in the Type I migration regime for most of the codes until it stalls not very far from the disc inner edge. We can therefore use the initial disc surface density profile for this comparison. For Run 4, there is a stronger surface density evolution, but we shall use the same approach (since two of the three population synthesis codes make such an approximation too), hoping that it will capture the essentials of the problem.

We first overview the clump migration approaches. Forgan & Rice (2013b) use the simplified migration scheme from Nayakshin (2010), in which the Type I migration timescale is

[TABLE]

This is derived from the Tanaka et al. (2002) formula (eq. 7) by requiring additionally a marginally unstable self-gravitating disc for which the Toomre parameter $Q\approx 1$ everywhere. For type II, the migration time scale is given by the disc viscous time,

[TABLE]

The switch between Type I and Type II migration occurs when $M_{\rm p}>M_{\rm t}$, where $M_{\rm t}$ is the transition mass given by,

[TABLE]

as used by Bate et al. (2003). We note that Forgan et al. (2018) have recently presented an updated population synthesis model. We do not include this study in our code comparison here because its migration module is similar to the Müller et al. (2018) treatment, which is discussed below. Furthermore, Forgan et al. (2018) also consider multiple gas clumps and model their N-body interactions. These effects can be very important in modifying the outcome of disc fragmentation (Hall et al., 2017) but is beyond the scope of our one-clump study.

Nayakshin & Fletcher (2015) use the Tanaka et al. (2002) expression for type I migration written as

[TABLE]

where $M_{\rm d}=2\pi\Sigma(R)R^{2}$ is a measure of the local disc mass, and $f_{\rm mig}$ is a dimensionless factor, set between 0.5 and 2 for different models. The factor is introduced to mimic the stochastic kicks from spiral density waves or other clumps. The Type II migration time is also set to the viscous time but with a correction multiplicative factor,

[TABLE]

The factor $(1+M_{\rm p}/M_{\rm d})$ takes into account planet inertia when the disc is less massive than the planet (Syer & Clarke, 1995). The correction is not very important for outer massive discs but may become large in the inner disc ($R\lesssim 10$ AU). The Crida parameter ($C_{\rm p}$, eq. 12) is used to model the transition between Type I and Type II migration. To prevent a sharp transition when $C_{\rm p}=1$, an exponential function of the form, $f={\rm min}(1,\exp[-(C_{\rm p}-1)])$ is used to smooth the transition out. Note that $C_{\rm p}$ is a function of the viscosity parameter $\alpha$, which is poorly known for protoplanetary discs. Nayakshin & Fletcher (2015) assumed that $\log\alpha$ is a random uniform variable in the limits between $\log(0.005)$ and $\log(0.05)$. We shall evaluate the results for these minimum and maximum values of $\alpha$. Finally, Nayakshin & Fletcher (2015) use a time-dependent 1D viscous disc model to evolve the disc surface density and other disc properties, and to conserve the angular momentum in the interactions between the disc and the planet, but for comparison below we shall assume the initial disc properties to be consistent with the two other models.

Müller et al. (2018) use a third set of equations to control planet migration, based on Baruteau et al. (2011), see eq. 8. For type II migration, eq. 25 is used but without the $(1+M_{\rm p}/M_{\rm d})$ correction, which however is unimportant for this paper as it is close to unity. Müller et al. (2018) also use the Crida parameter to determine when the planet switches to the type II migration, but consider two additional requirements for gap opening based on the work of Malik et al. (2015). They define three timescales, $\tau_{\rm visc}=R^{2}/\nu$, $\tau_{\rm cross}=2.5R_{\rm H}v_{r}^{-1}$, where $v_{r}$ is the radial velocity of the planet, and $\tau_{\rm gap}=q^{2}(H/R)^{5}\Omega^{-1}$. The additional requirements demand that $\eta\tau_{\rm gap}<\tau_{\rm cross}$ and $\tau_{\rm visc}<\tau_{\rm cross}$, where $\eta$ is a dimensionless factor varied from 10 to 1000, with $\eta=100$ used as a baseline model. Here we test only the first of these two additional criteria since it was the one used for most of the models in Müller et al. (2018).

Finally, population synthesis models also differ in how they treat the gas accretion onto clumps. Two of the population synthesis models (Forgan & Rice, 2013b; Nayakshin & Fletcher, 2015) neglected gas accretion onto the clumps, assuming a fixed gas mass unless the clumps are tidally disrupted. Müller et al. (2018) prescribed a gas accretion rate onto the clumps based on earlier simulations of Galvagni & Mayer (2014). Since our gas clumps accrete a significant amount of gas as they migrate, for a proper comparison with the population synthesis prescriptions we need all of them to take accretion into account. We previously found that the Bate et al. (2003) expressions for gas accretion rates, when reduced down to account for a smaller accretion efficiency of our slowly cooling discs, yields a reasonable match to the ccretion time scales of our simulation (fig. 7). We therefore use eq. 16 with $b=1/3$ here to let the planets gain mass when investigating Forgan & Rice (2013b); Nayakshin & Fletcher (2015) models.

We also need to take into account the decrease in the accretion rate when the planet switches from type I to the type II regime, which is clearly seen in fig. 7. To this end we write

[TABLE]

where $\dot{M}_{\rm acc}$ is the accretion rate estimate given by eq. 16 where $\Delta C=0.2$. We shall see below that this yields a decent fit to the planet mass evolution for both Run 3 and Run 4. With this approach, the comparison of population synthesis models to hydrodynamical simulations isolates just the planet migration and gap opening aspects.

Figs. 9 and 10 show such comparisons for Run 3 and Run 4, respectively. The shaded region represents approximately the range of numerical results obtained for these runs with the different numerical codes. In particular, PHATOM and GIZMO-MFM are selected to show the fastest and the slowest migrating planets for Run 3 in the left panel; the SEREN and GIZMO-MFM curves to show the range of models in the middle and right panels. For Run 4, PHANTOM and SEREN are selected as the extremes for the both planet accretion and migration tracks.

We see that there is a significant difference in how the three population synthesis models compare to the numerical results. The Müller et al. (2018) study appears to over-estimate somewhat how quickly and how far the planets migrate before they switch into the Type II regime. This seems to be because Müller et al. (2018) formulae are based on Baruteau et al. (2011) and yield too rapid migration by a factor of a few in the type I regime, as was seen in fig. 4. Also, the planet opens a gap a little closer to the star than it does in the simulations. This depends on the parameter $\eta$ which is set to 100 and 1000 for the solid and the dashed curves, respectively.

The Forgan & Rice (2013b) approach appears to yield too slow a migration rate. This is because the planet switches into the Type II migration rate immediately, as the transition mass in this approach is set to $2(H/R)^{3}M_{*}\approx 2{\,{\rm M}_{\rm J}}$, and the planet is already this massive in the beginning of the Run 3. This under-estimates the planet transition mass, which is found to be in the range of $\sim 7{\,{\rm M}_{\rm J}}$ to $\sim 30{\,{\rm M}_{\rm J}}$ (cf. fig. 6). We found that a much better fit to Run 3 is obtained with the Forgan & Rice (2013b) formulae if the transition mass is increased by a factor of $\sim 5$.

The Nayakshin & Fletcher (2015) formulae used Tanaka et al. (2002) expression for the migration rate with a dimensionless factor $f_{\rm migr}$ in front. The factor was a logarithmically uniform random variable in the limits $0.5<f_{\rm migr}<2$ and was meant to mimic possible stochastic kicks that the clumps obtain when interacting with the spiral density waves of the disc (see Baruteau et al., 2011). In the interest of figure clarity we use $f_{\rm migr}=1$ in fig. 9 for this model. Further, Nayakshin & Fletcher (2015) use the Crida et al. (2006) switch for gap opening, with the $\alpha$ parameter being a sum of two parts, a constant $\alpha$ and a part driven by self-gravity. We neglect the latter contribution to $\alpha$ here, and show two cases with $\alpha=0.005$ and $0.05$ in fig. 9. It is apparent that the smaller $\alpha$ curve (red solid) opens a gap in the disc far more easily than expected. The $\alpha=0.05$ curve (red dashed) seems more reasonable. However, we must remember that Nayakshin & Fletcher (2015) neglected gas accretion onto the planet. The agreement of their prescriptions with Run 3 would have been worse if we kept the planet mass fixed at $2{\,{\rm M}_{\rm J}}$.

Fig. 10 shows that for a more massive gas clump none of the population synthesis prescriptions fare particularly well. The Forgan & Rice (2013b) model and the low viscosity model of Nayakshin & Fletcher (2015) open a gap in the disc too early, as for Run 3. The higher viscosity model of Nayakshin & Fletcher (2015) does relatively well in terms of gap opening mass but over-estimates the speed with which the planet migrates in initially. The Müller et al. (2018) equations also yield clumps migrating in too rapidly, and the gap is opened too close in compared with numerical simulations.

We therefore conclude that matching numerical results with analytic expressions remain a problem. What is particularly alarming is that seemingly benign changes in the parameters of the population synthesis prescriptions (such as a factor of a few change in the planet transition mass) can yield planet migration rates different by $\sim$ two orders of magnitude as the planet switches into the Type II migration regime prematurely.

5 Discussion and Conclusions

5.1 Numerics

In this paper we set up four different simulations of a gas planet starting at an initial separation of 120 AU in a massive gaseous disc. These 4 Runs differed in treatment of gas accretion onto the planet and the initial planet mass. We then performed these simulations with seven different numerical codes in order to compare their results.

We find differences by a factor of $\sim 2$, and sometimes as large as 3, between different codes in the accretion and migration rates. A more detailed analysis using PHANTOM indicates that these differences are to a large degree due to variations in the artificial viscosity prescriptions between the codes, although other factors such as gravitational softening and sink particle treatment probably also contribute.

We also compared our results with the planet migration and accretion prescriptions from three previous population synthesis studies (§4 and figs. 9 & 10). The Forgan & Rice (2013b) approach is found to open deep gaps in the disc prematurely. Since planets migrate very slowly in the type II regime, this implies that this study may over-estimate the population of gas giants remaining at wide separation after gas discs are dispersed. The Müller et al. (2018) study, on the other hand, over-estimates the rate of inward migration of planetary mass clumps. The Nayakshin & Fletcher (2015) study fits the Run 3 results relatively well in the high viscosity case but not for the low viscosity case. In the latter case, clumps open deep gaps in the disc and tend to stall on wide orbits when they should migrate to smaller radii via Type I migration. However, all three population synthesis prescriptions fare poorly for Run 4 in which a more massive planet is considered. Additionally, Forgan & Rice (2013b) and Nayakshin & Fletcher (2015) neglect gas accretion onto clumps.

5.2 Observational implications

Recent observational surveys of solar type stars show that only a few % of such stars are orbited by massive planets or brown dwarfs on orbits larger than $\sim 10$ AU (e.g., Biller et al., 2013; Chauvin et al., 2015; Vigan et al., 2012, 2017). Let us call this fraction $N_{\rm present}$. This is a key constraint on the theory of planet and brown dwarf formation via gravitational instabilities of large massive gas discs. However, it is even more important to consider the frequency of such objects in a time-integrated sense, that is, the number of gas clumps formed by disc fragmentation per star. Let this fraction be $N_{\rm birth}$. The two fractions are clearly connected via

[TABLE]

where $P_{\rm surv}<1$ is the probability for a gas clump to survive to the present day at a wide separation.

Detailed calculations and population synthesis approaches are necessary to calculate $P_{\rm surv}$ accurately. Forgan & Rice (2013b) obtained $P_{\rm surv}\sim 1$, Nayakshin & Fletcher (2015) had $P_{\rm surv}\la 0.1$ (Nayakshin, 2016, found a yet smaller value, $P_{\rm surv}\sim 0.05$, when feedback effects of the luminous core onto the clump are included), and Müller et al. (2018) found $P_{\rm surv}\ll 1$ but noted that this depends strongly on model assumptions. Rice et al. (2015) in addition showed that N-body interactions with secondary stars may remove a number of wide separation planets, lowering the fraction of $P_{\rm surv}$ further in the post-disc dispersal phase.

Our simulations and population synthesis comparison (figs. 9 and 10) demonstrate that just varying the assumptions about the underlying physics of the disc or clumps by a factor of a few may influence the results very strongly. One has to also add to this that the exact birth mass of the fragments and the mass of the disc at which it fragments are not known to better than a factor of a few (e.g., Kratter & Lodato, 2016), and the evolution of the clump strongly depends on uncertain disc cooling and dust physics (Nayakshin, 2017b), radiative feedback from the clump (Nayakshin & Cha, 2013; Stamatellos, 2015; Mercer & Stamatellos, 2017), etc. Therefore, the uncertainty in $P_{\rm surv}$ at present is uncomfortably large. At this time we cannot rule out a survival probability that would imply $N_{\rm birth}>1$.

What is the best way forward in resolving these uncertainties? Clearly, theoretical and simulation efforts to constrain $P_{\rm surv}$ from first principles should continue. However, other indirect approaches can also help. If the migration processes allow GI planets to populate the whole range of separations between the stellar radius and their birth place, what would be the differences between the objects left behind from this migration and those made by Core Accretion? If we are able to understand these differences more robustly, then discovering (or not) such unusual objects at separations less than 10 AU may yield independent constraints on $P_{\rm surv}$ and $N_{\rm birth}$.

Acknowledgements

Theoretical astrophysics research at the University of Leicester is supported by a STFC grant. The work performed at the University of Leicester used the ALICE High Performance Computing Facility, and the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure.

FM acknowledges support from The Leverhulme Trust, the Isaac Newton Trust and the Royal Society Dorothy Hodgkin Fellowship. This work was undertaken on the COSMOS Shared Memory system at DAMTP, University of Cambridge operated on behalf of the STFC DiRAC HPC Facility. This equipment is funded by BIS National E-infrastructure capital grant ST/J005673/1 and STFC grants ST/H008586/1, ST/K00333X/1. This work also used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by a BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/K00087X/1, DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.