© ESO, 2015

1. Introduction

Since the dawn of observational exoplanet detection in recent decades the ever-growing number and variety of planetary species gives rise to renewed theoretical efforts for understanding the formation of these bodies (e.g., Boss 1997; Bodenheimer et al. 2000). The ultimate goal of achieving a consistent explanation for the planetary population (e.g., Benz et al. 2014) can only be adressed by pushing observational and theoretical tools to their limits.

The two major theories of planet formation, core accretion (CA) and formation via gravitational instabilities (GI), have the potential to explain a wide variety of planetary species. Whereas it is widely accepted that terrestrial planets are formed via a bottom-up model like CA (e.g., Raymond et al. 2006, 2014; Morbidelli et al. 2012), the formation of more massive objects, like gas giants, can in principle be understood in terms of CA (Helled et al. 2014) or GI (Mayer et al. 2003). Recent efforts (e.g., Boley 2009) are dedicated to a synthesis of both models in a more general framework of planet formation. Scenarios based on GI, however, can be particularly relevant for the formation of planets around metal-poor stars, where dust grain coagulation becomes increasingly difficult (e.g., Setiawan et al. 2010; Johnson et al. 2012), and the buildup of supermassive and Population III stars (e.g., Inayoshi & Haiman 2014; Latif & Schleicher 2015; Schleicher et al. 2015).

In this work we aim for an investigation of the capabilities of an alternative technical treatment of GI to achieve a better understanding of the buildup of massive planets. This includes gaseous planets like Jupiter and Saturn in our own solar system, as well as similar or even more massive planets in extrasolar systems (e.g., Borucki et al. 2010) or planets around exotic configurations like in the pre-cataclysmic binary system NN Serpentis. The latter system is an example for the class of evolved binaries, for which planets were detected (Zorotovic & Schreiber 2013; Beuermann et al. 2013). NN Serpentis recently underwent detailed investigations (e.g., Haefner et al. 2004; Hessman et al. 2011; Horner et al. 2012; Marsh et al. 2014; Völschow et al. 2014; Parsons et al. 2014), indicating a second-generation scenario where the planets were formed after a common envelope (CE) event. Subsequently, Schleicher & Dreizler (2014) argued for these planets to be consistent with the formation via GI and introduced a semianalytical model, which further motivates studies in compact self-gravitating disks. In the NN Serpentis scenario, the age of the white dwarf of 10⁶ yr provides an upper limit on the timescale for planet formation, therefore favoring GI with respect to CA and indicating a formation scenario in a very compact environment (<10 AU).

So far, the hydrodynamical treatment of GI in protoplanetary disks mostly relies on a description via smoothed-particle hydrodynamics (SPH; e.g., Mayer et al. 2007) or axisymmetric (e.g., Boss 2003; Boley et al. 2007; Meru & Bate 2011) static grid approaches without adaptive mesh refinement (AMR). Only a minority of grid approaches employ a mechanism for adaptive and dynamic refinement of the initial mesh (e.g., Paardekooper & Mellema 2006). In the case of Cartesian grids, the only attempt we know about is the project of Mayer & Gawryszczak (2008) and Gawryszczak et al. (2010), who compared flash and gasoline simulations of a self-gravitating disk. In general, the use of more than one numerical approach yields the advantage of ruling out systematic errors, introduced by the chosen numerical treatment (e.g., Agertz et al. 2007). Advantages of grid codes employing finite volume methods in general are their built-in option to conserve mass and momentum. Regarding the treatment of astrophysical phenomena they offer a precise modeling of turbulence in hydrodynamical instabilities (Agertz et al. 2007) and a more robust treatment of shocks in comparison with SPH approaches (Price & Federrath 2010). Additionally, classical SPH approaches tend to overrate viscosity parameters (Price 2008, 2012), which is not the case for grid-based implementations. However, recent developments in the SPH community extenuate many of these problems (see, e.g., Read & Hayfield 2012; Hopkins 2013). In all numerical codes the chosen geometry dictates the implementation and discretization of the physical equations. In the case of disk simulations, one could argue that the natural geometry of choice would be axisymmetric. However, we see some principle advantages in the use of Cartesian geometry, which motivates comparative studies. For example such codes do not include a built-in singularity (i.e., the center of origin) and handle all computed regions equally. The latter means there is no built-in refinement in inner regions, which might avoid symmetry-dependent solutions. Including the use of AMR schemes enables us to specificially choose regions of greater interest and neglect others. This can be particularly relevant when fragmentation occurs because of turbulent motion in the outer disk parts and when the deviations from axisymmetry are more prominent.

We use the Eulerian block-structured AMR hybrid code enzo¹ (O’Shea et al. 2004; Bryan et al. 2014), which enables massively parallel simulations and supports a wide variety of astrophysical problems, i.e., hydrodynamics, ideal and non-ideal magnetohydrodynamics, and N-body dynamics, including self-gravity for fluids and gas chemistry. The fluid flow in enzo is evolved using a finite-volume discretization, the ideal gas dynamics calculations are treated with the piecewise parabolic method (PPM), which is a higher-order Godunov scheme (Bryan et al. 1995) with an accurate piecewise parabolic interpolation and a non-linear Riemann solver for shock conditions. The solver is third-order accurate in space and second order accurate in time for fixed time stepping and formally second-order accurate for variable time stepping. enzo has been tested on a variety of typical fluid flow test problems (Bryan et al. 2014).

The code employs the implementation of a (structured) AMR scheme, first developed by Berger & Colella (1989), which enables higher resolution than affordable via uniform grids by introducing additional finer mesh structures on a coarse uniform grid when necessary (which can be defined with a variety of conditions). This scheme utilizes an adaptive hierarchy of rectangular grid patches, which cover a snippet of space with a certain resolution, refined from the top to the bottom of the hierarchy. When the solution evolves and interesting regions develop, finer meshes are placed within the coarse grid, enabling a higher resolution of these structures and not draining too much computing power for less interesting regions. This facilitates a variety of options for modeling GI in protoplanetary disks, in which we need to deal with dynamical effects on a large range of spatial scales.

We perform 3D simulations of gravitational instabilities in protoplanetary disks with varying disk radius and resolution settings using AMR techniques. In Sect. 2 we describe the physical and numerical setup of the simulations. In Sect. 3 we compare the initial and evolved state of stable and unstable disks and investigate the effect of different refinement settings. In Sect. 4 we verify the setup for various disk radii and demonstrate its ability to model systems of a wide range of disk extensions. Finally, we summarize our findings in Sect. 5.

2. Numerical experiments

2.1. Numerical setup

For a complete description of enzo’s code structure we refer to Bryan et al. (2014). To provide a sufficient explanation of our numerical methods, however, we briefly outline the major numerical recipes used in the calculations.

Hydrodynamic equations

In our simulations, we use the piecewise parabolic method (PPM), which was originally developed by Colella & Woodward (1984) and was adapted for the enzo code by Bryan et al. (1995) in a direct Eulerian fashion. In this implementation, the governing equations are dimensionally split and rewritten (as example in one dimension) in conservative form as $\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial \rho v}}{\mathit{\partial x}}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial \rho v}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial \rho }{\mathit{v}}^{\mathrm{2}}}{\mathit{\partial x}}\mathrm{+}\frac{\mathit{\partial p}}{\mathit{\partial x}}\mathrm{=}\mathit{\rho g,}\\ & & \frac{\mathit{\partial \rho E}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial \rho vE}}{\mathit{\partial x}}\mathrm{=}\mathit{\rho vg,}\end{array}$in which x and v refer to the one-dimensional position and velocity of the gas, g to the acceleration at the cell center, ρ to the gas density, and E to the total fluid energy density. The one-dimensional solution to these equations can be found by computing the effective left and right states at each grid boundary. This is done by constructing a piecewise parabolic solution of density, velocity and total energy and then averaging over the distance that each wave can travel in the specific time step. Then a Riemann problem with these effective left and right states is solved and all quantities are updated. The standard Riemann solver for the PPM method in enzo approximates the appearing rarefaction waves, which are traveling left and right, as shocks. The solution is then found using an iterative approach. For the case of occurring numerical problems due to the Riemann solver we used, enzo employs a fallback mechanism. This allows the code to switch the Riemann solver to the more diffusive HLL solver at particular interfaces, where negative densities or energies occur (see, e.g., Lemaster & Stone 2009). Finally, the fluxes in between the grid patches are constructed. The full, one-dimensional mathematical description and the interpolation methods are given in Colella & Woodward (1984) and Bryan et al. (1995).

Dual energy formalism

As long as the structures in the simulations are well resolved and the Mach numbers are on a reasonable scale (i.e., Ma< 100), the aforementioned solver works well. However, if hypersonic bulk flows with E_therm/E_tot ~ 10^-3 appear in the simulation, the numerical situation becomes disastrous: the ratio of kinetic energy E_k to gas internal energy e can approach numbers that are much too high. Hence, the pressure (proportional to E − E_kin) is the difference between two extremely large numbers, which would be a rather disadvantageous numerical situation with major sources for errors, especially for the temperature distribution. To overcome this problem, enzo separately solves the internal energy equation $\begin{array}{ccc}\frac{\mathit{\partial e}}{\mathit{\partial t}}\mathrm{+}{v}\mathrm{·}\mathrm{\nabla }\mathit{e}\mathrm{=}\mathrm{-}\frac{\mathit{p}}{\mathit{\rho }}\mathrm{\nabla }\mathrm{·}{v}\mathit{,}& & \end{array}$(4)with thermal energy density e and pressure p, which is then used in hypersonic flows for calculating the pressure and thus the temperature distribution only. The gas dynamics (i.e., the hydrodynamic routines) should be unaffected to avoid introducing new errors. This is achieved by choosing a selection criterion for the pressure and thus separating both formalisms via $\begin{array}{ccc}\mathit{p}\mathrm{=}\left\{\begin{array}{c}\\ & \\ & \end{array}\right\}\mathit{.}& & \end{array}$(5)If the barrier parameter η is chosen to be small enough, the dynamical effect of the dual energy formalism is approximately zero. Bryan et al. (2014) list η = 10^-3 as the used standard value, consistent with usual truncation errors in the simulations.

Gravity

Of particular importance for our planet formation model is the inclusion of the self-gravity of the gas. To achieve this, Poisson’s equation $\begin{array}{ccc}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\phi }\mathrm{=}\mathrm{4}\mathit{\pi G\rho ,}& & \end{array}$(6)with the gravitational potential φ, the gravity constant G, and the gas density ρ, is solved in a multistep process for all cells. The cloud-in-cell (CIC) interpolation method from Hockney & Eastwood (1988) is used to approximate the gas distribution as a spatially-discretized density field, from which the gravitational potential is generated via a fast Fourier transform. To get to know the acceleration at each (sub)grid a finite difference scheme with Dirichlet boundary conditions is used and protected against oscillations through buffer zones surrounding the active grid zones. Because of approximations in this process, Bryan et al. (2014) estimates the resolution of the gravitational force to be approximately twice as coarse as the corresponding refinement level.

Equation of state

The total fluid energy density is given by $\begin{array}{ccc}\mathit{E}\mathrm{=}\mathit{e}\mathrm{+}\frac{\mathit{\rho }{\mathit{v}}^{\mathrm{2}}}{\mathrm{2}}\mathit{,}& & \end{array}$(7)with the thermal energy density e governed by the equation of state $\begin{array}{ccc}\mathit{e}\mathrm{=}\frac{\mathit{p}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathrm{·}& & \end{array}$(8)This is the equation of state for an ideal gas with the ratio of specific heats γ.

Computational domain

To prevent interactions of the disk material with the boundaries of the computational domain, we choose the physical size of the box to be (R_disk × 10)³. Thus, there is enough space for the disk to evolve and to be embedded in a diffuse medium. The boundaries of the domain behave as reflecting walls, i.e., $\begin{array}{ccc}\mathit{q}\mathrm{(}\mathrm{-}\mathit{x}\mathrm{)}\mathrm{=}\mathit{q}\mathrm{(}\mathrm{+}\mathit{x}\mathrm{)}\mathit{,}& & \end{array}$(9)with an arbitrary field quantity q, and with the velocity perpendicular to the boundary direction reversed $\begin{array}{ccc}{\mathit{v}}_{\mathit{x}}\mathrm{(}\mathrm{-}\mathit{x}\mathrm{)}\mathrm{=}\mathrm{-}{\mathit{v}}_{\mathit{x}}\mathrm{(}\mathrm{+}\mathit{x}\mathrm{)}\mathit{.}& & \end{array}$(10)

Refinement criteria

The refinement to the next (deeper) level of the grid patches in our simulations depends on the corresponding cell mass and a so-called Jeans length criterion. The cell mass refinement works such that it mimics a Lagrangian method in trying to keep a fixed mass resolution. Thus, if the mass of a specific cell is $\begin{array}{ccc}{\mathit{M}}_{\mathit{g}}\mathrm{=}\mathit{\rho }\mathrm{(}\mathrm{\Delta }\mathit{x}{\mathrm{)}}^{\mathit{d}}\mathit{>}{\mathit{\rho }}_{\mathrm{flag}}\mathrm{(}\mathrm{\Delta }{\mathit{x}}_{\mathrm{root}}{\mathrm{)}}^{\mathit{d}}{\mathit{r}}^{{\mathit{ϵ}}_{\mathit{l}}\mathit{l}}\mathrm{=}{\mathit{M}}_{\mathrm{Jupiter}}\mathit{/}\mathrm{1000}\mathit{,}& & \end{array}$(11)the cell is refined to a deeper level (if allowed by the maximum refinement level). Here, ρ_flag is the required density on the root grid, Δx_root is the root grid cell spacing, r the refinement factor (usually 2), l the level and ϵ_l is an aggressiveness parameter, to make the refinement super-Lagrangian (ϵ_l< 0) or sub-Lagrangian (ϵ_l> 0).

The second approach refines the Jeans length $\begin{array}{ccc}{\mathit{\lambda }}_{\mathrm{J}}\mathrm{=}{\left(\frac{\mathit{\pi }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{\mathit{G\rho }}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}& & \end{array}$(12)by a fixed number of cells (based on Truelove et al. 1998) if $\begin{array}{ccc}\mathrm{\Delta }\mathit{x}\mathit{<}{\left(\frac{\mathit{\pi }{\mathit{k}}_{\mathrm{B}}\mathit{T}}{{\mathit{N}}_{\mathrm{J}}^{\mathrm{2}}\mathit{G\rho }{\mathit{m}}_{\mathit{H}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(13)with N_J = 32 the required number of cells per Jeans length. This is especially valuable since we deal with fragmentation effects, where contraction plays a key role in the later building of clumps.

2.2. Initial conditions

3D density profile

To begin with, we set the column density structure of the disk according to a classical Mestel power-law profile (Mestel 1963) $\begin{array}{ccc}\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{0}}\frac{{\mathit{r}}_{\mathrm{out}}}{\mathit{r}}\mathit{,}& & \end{array}$(14)with r_out denoting the outer radius of the disk and Σ₀ the column density at this radius. Since the contribution to mass and angular momentum at small disk radii can be neglected we assume that the disk extends to r = 0. This leads to a normalization of the column density as $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\frac{{\mathit{M}}_{\mathrm{disk}}}{{\mathit{r}}_{\mathrm{out}}^{\mathrm{2}}}\mathit{,}& & \end{array}$(15)with the overall disk mass chosen to be M_disk ≈ 0.05 M_⊙. Using hydrostatic balance the analytical solution to resolve the Gaussian density distribution in the vertical direction is $\begin{array}{ccc}\mathit{\rho }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{H}\mathrm{\left(}\mathit{r}\mathrm{\right)}\sqrt{\mathrm{2}\mathit{\pi }}}\exp \left(\mathrm{-}\frac{{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}\mathit{H}\mathrm{(}\mathit{r}{\mathrm{)}}^{\mathrm{2}}}\right)\mathrm{·}& & \end{array}$(16)Here, z denotes the distance to the midplane of the disk and H the scale height, which is evaluated as (Lodato 2008) $\begin{array}{ccc}\mathit{H}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& \end{array}$(17)with the scale height for self-gravitating disks $\begin{array}{ccc}{\mathit{H}}_{\mathrm{sg}}\mathrm{=}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}\mathit{\pi G}\mathrm{\Sigma }\mathit{,}& & \end{array}$(18)and the scale height for Keplerian disks $\begin{array}{ccc}{\mathit{H}}_{\mathit{\star }}\mathrm{=}{\mathit{c}}_{\mathrm{s}}\mathit{/}\mathrm{\Omega }\mathit{.}& & \end{array}$(19)In the equation above, the sound speed is given via the Newton-Laplace equation $\begin{array}{ccc}{\mathit{c}}_{\mathrm{s}}\mathrm{=}\sqrt{\mathit{\gamma }\mathrm{·}\frac{\mathit{p}}{\mathit{\rho }}}\mathit{,}& & \end{array}$(20)and the angular velocity Ω (see next subitem). As mentioned by Lodato (2008), this approximation holds for a comparable gravitational influence of the disk and the central object, H_sg≈H_⋆.

Central object

We choose the host star of the system to be of subsolar mass $\begin{array}{ccc}{\mathit{M}}_{\mathit{\star }}\mathrm{=}\mathrm{0.646}{\mathit{M}}_{\mathrm{\odot }}\mathit{,}& & \end{array}$(21)which is represented in the simulations as a point mass particle without spatial extension (denoted in the later figures with a single black dot in the middle) and is included in the gravity calculations.

Orbital velocity profile

Including the disk’s self-gravity, the azimuthal velocity satisfies (Lodato 2008) $\begin{array}{ccc}\frac{{\mathit{v}}_{\mathit{\phi }}^{\mathrm{2}}}{\mathit{r}}& \mathrm{=}\frac{\mathit{\partial }{\mathrm{\Phi }}_{\mathrm{tot}}}{\mathit{\partial r}}\mathrm{+}\frac{\mathrm{1}}{\mathit{\rho }}\frac{\mathrm{d}\mathit{p}}{\mathrm{d}\mathit{r}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathit{\star }}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{+}\mathrm{2}\mathit{\pi G}\mathrm{\Sigma }\mathrm{+}\frac{\mathrm{1}}{\mathit{\rho }}\frac{\mathrm{d}\mathit{p}}{\mathrm{d}\mathit{r}}\mathit{,}& \end{array}$(22)with Φ_tot the gravitational potential of central object and disk material.

Thermal profile

Finally, we need to address the thermal profile of the disk. We characterize the disk through the Toomre parameter (Toomre 1964) $\begin{array}{ccc}\mathit{Q}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{s}}\mathit{\kappa }}{\mathit{\pi G}\mathrm{\Sigma }}\mathrm{\approx }\frac{{\mathit{c}}_{\mathrm{s}}\mathrm{\Omega }}{\mathit{\pi G}\mathrm{\Sigma }}\mathit{,}& & \end{array}$(23)with Boltzmann’s constant k_B. Fixing Q to a constant value then translates into a condition for the sound speed and therefore the temperature of the gas. The epicyclic frequency κ can be approximated with the angular velocity Ω = v_φ/r, when the velocity profile throughout the disk is mostly dominated by the Keplerian velocity (Lodato 2008). We assume a gas with composition similar to the solar system, i.e., the mean molecular weight μ is set to 2.4 (Mayer et al. 2007). From this assumption we infer the temperature profile using Eq. (20) to be $\begin{array}{ccc}\mathit{T}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{\mu }{\mathit{m}}_{\mathrm{p}}}{\mathit{\gamma }{\mathit{k}}_{\mathrm{B}}}\mathrm{·}& & \end{array}$(24)The ratio of specific heats is fixed to be γ = 1.0001 ≈ 1, according to an isothermal equation of state (EOS).

3. Resolution study

In this chapter, we investigate the imprint of the refinement level on the physical conditions within the disk. We benchmark the simulations for disks with R_disk = 100 AU for the initial conditions and at t = 1.5 T_disk, i.e., 1.5 orbital times of the most outer radius of the disk. An overview of all of our featured simulations is given in Table 1.

Fig. 1 Radial profiles for initial column density Σ and rotational velocity vrot for simulations with Rdisk = 100. The parameter gc indicates the resolution of the coarsest grid, lmax the deepest level of the resolution by dynamic refinement. The continous lines show simulations with Qinit = 1, the dashed lines simulations with Qinit = 3. The density distribution converges toward the analytic solution for higher refinement levels, i.e., the density gap in the inner part of the disk is smaller for higher refinement levels. Additionally, it is smaller for the simulation runs with Qinit = 1, which is due to the scale height dependency on the Q parameter (compare Eq. (23) and Fig. 3). The velocity, only nonzero within the disk matter, behaves analogous to the density.

Fig. 2 Radial profiles for initial temperature T and effective Toomre Q parameter for simulations with Rdisk = 100. The temperature increases outward because of its dependence on Q (see Eq. (24)). Q is approximated correctly in the inner parts of the disk. In outer parts, where the density decreases for numerical reasons, Q rises, and thus artificially stabilizes the disk matter.

3.1. Initial state

Figures 1 and 2 show comparisons of the radial profiles of the values of Σ,v_rot,T and Q of disks with R_disk = 100 AU and Q_init = 1, Q_init = 3, respectively. All of them converge for higher maximum refinement levels toward the analytical solution, as given in Sect. 2.2. Some differences between the simulations with Q_init = 1 and Q_init = 3 arise because of the change of scale height H by a change in Q, according to Eq. (19). It is important to say that the simulations with the same effective resolution behave very similarly, despite any differences in resolution for more coarse grids. Therefore, we restrict the discussion on the simulations with g_c = 64³, while analogous simulations with the same effective resolution can be identified for the runs with g_c = 128³. In the following we discuss and explain the differences in detail.

Fig. 3 Initial vertical density profile Σ, comparing disks with Rdisk = 100 AU and Qinit = 3 (top), Qinit = 1 (bottom), and increasing maximum refinement level from 2−4 (left to right), respectively. The scale height H, according to Eq. (19), changes with Q, which is why the Qinit = 1 disks show a flatter profile. The gray boxes indicate the deepness of the refinement level, from level l = 0 outside of the boxes to level l = lmax in the most inner box.

The profiles for the surface density Σ yield two conclusions. First, higher refinement means that the disk is much better resolved, especially in the inner part. Second, it displays the dependency of H on Q. This is illustrated in Fig. 3, which shows slices of the initial density profile in a cut through the vertical profile of the disk, comparing disks with R_disk = 100 AU and Q_init = 1, Q_init = 3, respectively. The different refinement levels are indicated with gray boxes. These show that the disk midplane is better resolved than the outer parts, thus enabling a higher resolution in the regimes important for fragmentation. The vertical disk extension is no longer resolved when the scale height drops below the cell size at small radii.

All disks represent the analytical solution for the rotational velocity very well, with only minor differences. They are better resolved (i.e., better converge toward the analytic solution in the inner regions) for a higher refinement level. This is a direct consequence from the density distribution (Figs. 1 and 3), since the rotational velocity is only defined in the disk material.

The temperature profiles show fundamental differences for the Q_init = 3 and Q_init = 1 cases as a result of the Q dependence of the thermal profile (Eq. (24)). Fixing the Toomre parameter to Q_init = 1 results in lower temperatures throughout the disk and thus an increased chance for developing fragmentation. As a result, the temperature distributions show a positive slope with increasing radius, which is adopted to explore disks in a marginally stable state. Outside of the disk the temperature rapidly increases to the temperature of the background gas T_background = 200 K.

Finally, the effective Toomre Q distributions show the numerical imprint of Q_init on the disk as a consequence of discretization effects. Thus, Q drastically increases in the inner and outer parts of the disk, where the density per cell decreases for numerical reasons. Despite these limitations Q stays constant in the dynamically relevant disk parts, where the density is high enough to support fragmentation effects.

3.2. Evolved state

A potential weakness of Cartesian grid codes is angular momentum dissipation, which is especially relevant for modeling systems with high rotational velocities. Therefore, to demonstrate the numerical robustness of our setup, we show the evolution of the total angular momentum L_total in the whole simulation box in Fig. 4 for simulation times up to t = 1.5 T_disk.

Fig. 4 Evolution of angular momentum for all simulations with Rdisk = 100 and up to t = 1.5 Tdisk. The maximum deviations from the initial state reach up to 2%.

Fig. 5 Face-on column density profiles for simulations with gc = 643 and t = 1.5 Tdisk. Top row: simulations with Qinit = 3, bottom row: simulations with Qinit = 1. The maximum refinement level lmax is increasing from left to right from 2−4.

Fig. 6 Edge-on volume density profiles for simulations with gc = 643 and t = 1.5 Tdisk. Top row: simulations with Qinit = 3, bottom row: simulations with Qinit = 1. The maximum refinement level lmax is increasing from left to right from 2−4.

Overall, the deviations for all simulations go at maximum up to +2% of the initial angular momentum. For the simulations with Q_init = 3, the angular momentum shows a very stable configuration and all simulations show the same general trend. However, the simulations with Q_init = 1 differ from each other, depending on the deepest refinement level in the simulation. Here, the simulations with the highest resolutions (i.e., with deepest refinement level l_max = 3 or l_max = 4) show a decrease of up to − 1.5% total angular momentum.

Fig. 7 Radial profiles of column density Σ and rotational velocity vrot for all simulations with Rdisk = 100 and t = 1.5 Tdisk. The runs with Qinit = 3 approximately follow the initially chosen density distribution. The Qinit = 1 disks show peaks at radii with distinct clumps or spiral structure (compare Figs. 5 and 6). The velocity profiles show only minor deviations from the general trend, which are due to the turbulent velocity in the most unstable regions.

Fig. 8 Radial profiles of the temperature distribution and effective Toomre Q parameter for all simulations with Rdisk = 100 and t = 1.5 Tdisk. As the density distribution the temperature shows a spread-out for the inner and outer parts. This is much more pronounced for the simulations with Qinit = 1, whereas the runs with Qinit = 3 roughly follow the initial distribution. For Qinit = 3, the effective Q parameter stays well above the threshold for marginal stability. In comparison, the effective Q in the Qinit = 1 runs dips below the threshold at the positions with increased density (i.e., fragmenting regions, compare Fig. 7).

Now we focus on the imprint of changes in resolution on the evolved state of the disks for t = 1.5 T_disk to spot deviations for different refinement depths. Figure 5 visualizes the density distribution for a face-on view of the disks for various refinements for Q_init = 3 and Q_init = 1. The disks with Q_init = 3 do not differ qualitatively with increasing refinement level, but show increased resolution and details in the gas streams.

However, the differences for the Q_init = 1 disks are dramatic. Whereas the disk with l_max = 2 does not show any signs of clumping or spiral arms, these characteristics start to appear for l_max = 3 and are sharply defined in the simulation with l_max = 4. The differences due to l_max do not only change the overall density distribution of the disk. In addition, the planetary (clump) wake, associated with each forming clump, differs dramatically from l_max = 3 to l_max = 4. For l_max = 3, the clump only marginally disturbs the large-scale structure of the spiral arms, forming regions of higher density at approximately 50 and 80−90 AU separation. In contrast, the simulation with l_max = 4 shows very distinct spots, where fragmentation occurs. These clumps massively disturb the spiral arm structure, caused by the turbulent motion in the fragmenting areas.

The edge-on views of the density distribution of the disks (Fig. 6) reveal structural differences in the vertical gas distribution with increasing refinement level. The Q_init = 3 disks again show relatively weak differences, however, the inner part of the disk flattens with increasing refinement and thus the flared disk profile is better resolved.

The Q_init = 1 disks show this flattening as well. Moreover, the density peaks in the midplane parts of the disks are more pronounced with higher refinement level. Comparing the edge-on slices of the disks with the face-on projections demonstrates the capabilities of three-dimensional simulations. In areas of strong fragmentation the disk’s vertical extension is lower compared to more stable regions. Thus, the flared profile of the disk is disturbed and shows strong deviations from vertical and axisymmetry.

The visual observations are supported by the quantitative analysis of the radial profiles in Figs. 7 and 8.

Figure 7 shows the density distributions of the evolved state. The simulations with Q_init = 3 roughly follow the initally chosen density profile. In the inner and outer parts the profile is marginally spread out. The increase of density in the inner parts is the result of our analytic velocity approximation. Therein we assume the disk to be extended all the way down to the center. However, since the disk’s inner extension is cut off at some radius (see Sect. 3.1) the pressure gradient at the inner radius of the density distribution is too high. Therefore, the effective velocities in the inner part of the numerical implementation are too small. Thus, some of the gas on lower orbits rapidly migrates toward the center until an equilibrium state is reached. A similar effect occurs in the outer parts, where the gas is spread out until a smooth transition from high to low density is reached.

In comparison, the density distributions of disks with Q_init = 1 show deviations from the initially chosen profile. The turbulent and self-gravitating disk material induces stochastic fluctuations in the velocity. This, in turn, initiates the development of spiral arms and fragmentation, which is reflected in the radial profiles as peaks. Overall, the most massive clump can be seen for the simulation with l_max = 4 at ~52 AU, reaching column densities up to ${\mathrm{10}}^{\mathrm{3}}\frac{\mathit{g}}{\mathrm{c}{\mathrm{m}}^{\mathrm{2}}}$.

As discussed above, the velocity profiles for all disks reach an equilibrium state, which follows the analytical solution in most parts of the disk. The small deviations in the distribution for the Q = 1 runs emerge from the spiral arms and fragmentation, which are more pronounced for higher maximum refinement.

The temperature distribution in Fig. 8 behaves analogously to the density in Fig. 7. Thus, it features a spread in the inner and outer parts, more remarkable for lower resolutions. Higher resolutions however are able to sustain the initial profile better, which is true for both Q_init = 3 and Q_init = 1 simulations. All of them reflect the spread of the gas and feature an approximately similar temperature distribution to the outer edge.

The inital Q state is well preserved for the Q_init = 3 runs. However, the strong instabilities and the resulting fragmentation in the simulations with Q_init = 1 lead to minima in the effective Q below the unity threshold where the density reaches its highest values.

Fig. 9 Number of lowest-level clumps for all simulations and t = 1.5 Tdisk. The trend goes for just a few clumps for low maximum resolutions to many clumps for high maximum resolution. Different combinations of gc and lmax but comparable resolution on the lowest level (i.e., green/cyan and blue/magenta) show approximately equal numbers of lowest-level clumps.

As an attempt to quantify the behavior of the fragmentation for different refinement settings, we use the clump finding method of Smith et al. (2010) to detect topologically disconnected structures within the dataset. Its principle mechanism is to create a single contour in between ρ_min and ρ_max over the whole computational domain. The lower value is then continously increased, until it reaches the maximum value. During that process disconnected structures are identified as separate contours and are treated as individual structures in which the routine continues recursively. To get an idea of the smallest structures in our simulations, we only print the far-end leaves (“clumps”) of this treelike structure, i.e., topologically disconnected regions featuring approximately the highest density in their specific contours.

Figure 9 illustrates our findings, showing a histogram with statistics for all simulations with Q_init = 1. The number of clumps increases with the maximum level of refinement, i.e., we find more clumps for higher resolution. This underlines the importance of a thorough resolution to achieve correct structures throughout the grid.

4. Parameter study

Fig. 10 Radial profiles at t = 1.5 Tdisk for column density Σ and rotational velocity vrot for all simulations with the highest resolution on the finest grid, normalized in radial direction. The simulation with Qinit = 3 is added for comparison. The density distributions for all radii show similarities. However, the simulation with r = 100 AU and Qinit = 1 shows its single peak at ~52 AU and the negative slope at the outer edge of the disk is much shallower than the slope of the other simulations with Qinit = 1. The velocity distributions for all radii show a similar trend. The simulation with r = 10 AU abruptly goes down in velocity. This behavior might be due to the very sharp transition from disk material (with approximately Keplerian velocity profile) to surrounding material (with approximately zero velocity).

Fig. 11 Radial profiles at t = 1.5 Tdisk for temperature T and effective Toomre Q parameter for all simulations with the highest resolution on the finest grid, normalized in radial direction. The simulation with Qinit = 3 is added for comparison. Again, we see a different behavior of the simulation with Rdisk = 100 AU and Qinit = 1, displaying a smaller increase in temperature because of its smoother transition zone. The dips where the Q parameter goes below 1 are consistent with the density peaks in Figs. 10 and 12.

In this section we study the influence of the disk extension on the resolution dependent outcome of our model. In general, observations of protoplanetary disks show a wide variety of disk morphologies (e.g., Avenhaus et al. 2014; Garufi et al. 2014). Therefore, disk characteristics are likely to vary in radial extent. Additionally, some exotic systems like NN Ser are expected to host very compact systems of only a few AU radius (Schleicher & Dreizler 2014). Thus, we test our setup on its ability to model disks of varying radius. Note that all our simulated disks have the same overall mass M_disk and because of the Q dependency of our model, the temperature is adjusted to be able to reproduce global Q unstable configurations.

Figures 10 and 11 show radial profiles for the parameters Σ, v_rot, T and effective Q for all simulations with the highest resolution, i.e., with maximum refinement level l_max = 4. In principle, the simulations with different initial disk radii reproduce similar results and the distributions of the various physical parameters show similar trends. However, there are some differences, which mainly arise because of the single peak in the density distribution of the simulation with R_disk = 100 AU and Q_init = 1 at around 52 AU separation from the central object. This peak, resulting from a vary massive clump (compare Fig. 12), influences the whole disk structure, and drains material from the other parts of the disk, which results in a different average distribution of matter in the disk, compared to the simulations with R_disk = 10 AU and R_disk = 300 AU. While we expect this peak essentially to represent the stochastic nature of fragmentation, its presence leads to differences for many other parameters as well. For example, the velocity decrease is shallower, especially in comparison with the R_disk = 10 AU run, which shows a sharp transition at its outer edge from disk material to surrounding material. Addtionally, the temperature distribution shows an approximately constant distribution up to 1.3 R_disk, because the (cold) disk material (in comparison with the surrounding gas) is extended over a larger radius and therefore the transition from disk gas to surrounding gas happens at larger radii.

The Toomre Q radial profiles can be matched very well with the profiles for density and temperature. Thus, the most pronounced dip in the Q radial profile is found at approximately 0.5 R_disk separation for the simulation with R_disk = 100 AU. The simulations with R_disk = 10 AU and R_disk = 300 AU show a very similar (radially) normalized behavior. Their most pronounced Q dips are at separations 0.5 R_disk and 0.8−0.9 R_disk. Matching this observation with Fig. 12, we see that the highest amount of fragmentation can be spotted in these regions.

Figure 12 shows a final impression of the evolved structures for all simulations. The turbulent structures of the disk can be seen in the column density, as well as in the volume density plots. The projected face-on views provide detailed insight into the fragmented disk structures. The edge-on views represent very well the flared disk structures and reveal the regions of the highest density, where the most massive clumps are located.

5. Conclusions

We performed simulations of self-gravitating, massive protoplanetary disks using the AMR code enzo. Our physical setup was motivated by the semianalytical approach of Schleicher & Dreizler (2014) for modeling the characteristics of a formation via GI in diverse planet formation environments. We calculated the time evolution of disk configurations close to the Q_init = 1 threshold with disk radii of R_disk = 10,100 and 300 AU, varying grid settings, and the time evolution of a stable disk (Q_init = 3) with R_disk = 100 AU to validate our approach.

Fig. 12 Top row: density projections of all disks with Qinit = 1 and lmax = 4, t = 1.5 Tdisk. Bottom row: vertical density slices of all disks with Qinit = 1, lmax = 4 and t = 1.5 Tdisk.

The Q_init = 1 disks display the onset of large-scale GI on the orbital timescale, yielding gas fragmentation and the formation of clumps throughout the disk material. The imprint of AMR effects plays a major role in correctly modeling the fragmentation in the disk. Only the highest resolution levels yield pronounced signals of fragmentation and single peaks in the developing spiral arms. Additionally, the structure of the spiral arms and the fragmentation process is qualitatively different for different refinement levels. Whereas for low resolution the spiral arms are less pronounced and the forming clumps differ only slightly from their surrounding medium, with higher resolution the clumps build sharp column density peaks.

To summarize these findings, we are able to define a lower limit for the resolution in global AMR simulations of GI induced fragmentation in self-gravitating protoplanetary disks. A comparison of the imprint of the maximum refinement level on the fragmentation structure and the building of clumps (Figs. 5 and 9) yields a minimum ratio of fragmenting disk radius R_disk to resolution level r_lvl (in physical units). Thus, to induce the development of spiral arms in the disk this ratio has to be $\begin{array}{ccc}{\mathit{R}}_{\mathrm{disk}}\mathit{/}{\mathit{r}}_{\mathrm{lvl}}⪆\mathrm{50}\mathit{,}& & \end{array}$(25)corresponding to our simulations with l_max = 3. However, to be able to reveal further fragmentation and the building of distinct clumps in the disk spiral structures and resolve the planetary wake around them, one has to satisfy the more rigid criterion of $\begin{array}{ccc}{\mathit{R}}_{\mathrm{disk}}\mathit{/}{\mathit{r}}_{\mathrm{lvl}}⪆\mathrm{100}\mathit{,}& & \end{array}$(26)which corresponds to our simulations with l_max = 4. These critera, giving the number of horizontal cells the disk radius should contain, should be understood as a lower limit to resolve fragmentation structures in a global disk view. Running higher resolutions will most likely yield even better and more resolved structures.

As another valuable feature of the 3D implementation of our setup, disks of all radii conserve their flared structure very well and resolve deviations in scale height during the fragmentation process. This might be especially important to investigate deviations from vertical and axisymmetry in the planet formation process. This is in agreement with the findings of Mayer & Gawryszczak (2008), who underline the importance of (initial) resolution in the vertical direction to resolve clumping effects. Even if it is in general hard to compare the outcome of SPH calculations with those from grid-based models, our resolution dependent results are similar to those of Nelson (2006) concerning the importance of vertical resolution in SPH simulations.

The total angular momentum in all simulations is approximately conserved (with total errors of up to 2%) during the time evolution of the simulations. The marginal loss in angular momentum, shown in Fig. 4, can be explained with discretization effects due to the grid structure of the simulations. From a numerical point of view the simulations presented here can be understood as a first step toward a more realistic coverage of turbulence effects in accretion disk scenarios with high Mach numbers. As found by Federrath et al. (2011) the energy injection scale of gravity-driven turbulence is close to the local Jeans scale. Therefore, in AMR simulations in which gravitational energy is converted into turbulent energy, it is crucial to resolve the Jeans length by (at least) 32 (Federrath et al. 2011) or even 64 (Latif et al. 2013a,b) cells. Therefore, the maximum refinement level chosen in our simulations is a constraint for even better angular momentum conservation.

Additional ideas to speed up simulations of self-gravitating accretion disks and to enable simulations on a wider range of spatial scales are necessary to eventually cover the range from the global disk view down to smaller scales as, e.g., the physics of circumplanetary disks (Morbidelli et al. 2013). A first step toward this goal might be the introduction of a more specific refinement criterion, which only covers areas that can be associated with long-lasting clumps, which might eventually end up on a planetary mass scale.

A particular interesting feature of the simulations from a physical point of view is the formation and rapid inward migration of clumps. They migrate inward within a few orbital times, which is roughly consistent with the timescale of Type I migration (Helled et al. 2014). Apart from the relatively small masses, the rapid inward migration is consistent with the findings of Michael et al. (2011) and Baruteau et al. (2011). The migration in unstable and turbulent disks is driven by several forces, such as clump-disk interactions (like Lindblad torques) or stochastic torques. Identifying the major driving forces in our scenario and to find a stabilization mechanism for the migrating clumps are subjects for further research.

To overcome these fast migration scenarios, Meru et al. (2014) proposed the sculpting of long-lasting dust rings through pressure maxima in the disk. However, a scenario to explain the migration stalling purely based on the thermal evolution of the gas might be the influence of radiative feedback from the central star. Following Chiang & Goldreich (1997), the photoheating would transform a significant part of the disk to be stabilized against fragmentation instabilities (Schleicher & Dreizler 2014). Clumps migrating toward such stabilized regions might still be able to grow and evacuate the disk in this region. This stalling of the inward migration and possible opening of gaps in the disk would prevent them from being destroyed by tidal interactions (Zhu et al. 2012).

The timescale of the disk’s evolution, comparable to the orbital timescale, and the fast fragmentation and formation of clumps is of major interest for the formation of planets in exotic systems like the NN Serpentis binary system. Since its evolutionary timescale is determined by the cooling age of its white dwarf (Hessman et al. 2011), the formation of the planets must happen on a timescale of 10⁶ yr. Thus, the rapid fragmentation and clumping opens up a possible way to create the proposed NN Ser planets in the appropriate time. However, further observational coverage and more realistic theoretical investigations of this system are needed.

We have presented here a systematic study exploring the GI in self-gravitating protoplanetary disks with AMR simulations. We expect that this technique can provide additional insight into the formation of massive self-gravitating clumps in future simulations, which will complement the existing numerical approaches well. This technique can be combined with additional physics modules like the chemistry package krome (Grassi et al. 2014), radiation transport techniques (Wise & Abel 2011) and a sink particle algorithm (Latif et al. 2013c) for an improved modeling of the formation of planets. The simulations presented here will be particularly valuable as a reference model to understand the evolution of the GI in the absence of additional processes that can complicate the situation, but already shed light on the important role and necessity of high resolution in grid-based calculations.

Acknowledgments

We thank the anonymous referee for valuable suggestions, which helped to improve the quality of the paper. Moreover, we thank Stefan Dreizler, Robi Banerjee, Rick Hessman, Farzana Meru, Marcel Völschow, Christiane Diehl and Muhammad Latif for stimulating discussions. The computer simulations featured in this work were run on the GWDG computing cluster (gwdg.de). The simulation results were analyzed using the visualization toolkit for astrophysical data yt (Turk et al. 2011, yt-project.org), which is based on the matplotlib Python library (Hunter 2007, matplotlib.org).

References

All Tables

All Figures

Fig. 1 Radial profiles for initial column density Σ and rotational velocity vrot for simulations with Rdisk = 100. The parameter gc indicates the resolution of the coarsest grid, lmax the deepest level of the resolution by dynamic refinement. The continous lines show simulations with Qinit = 1, the dashed lines simulations with Qinit = 3. The density distribution converges toward the analytic solution for higher refinement levels, i.e., the density gap in the inner part of the disk is smaller for higher refinement levels. Additionally, it is smaller for the simulation runs with Qinit = 1, which is due to the scale height dependency on the Q parameter (compare Eq. (23) and Fig. 3). The velocity, only nonzero within the disk matter, behaves analogous to the density.

In the text

	Fig. 2 Radial profiles for initial temperature T and effective Toomre Q parameter for simulations with Rdisk = 100. The temperature increases outward because of its dependence on Q (see Eq. (24)). Q is approximated correctly in the inner parts of the disk. In outer parts, where the density decreases for numerical reasons, Q rises, and thus artificially stabilizes the disk matter.
In the text

Fig. 3 Initial vertical density profile Σ, comparing disks with Rdisk = 100 AU and Qinit = 3 (top), Qinit = 1 (bottom), and increasing maximum refinement level from 2−4 (left to right), respectively. The scale height H, according to Eq. (19), changes with Q, which is why the Qinit = 1 disks show a flatter profile. The gray boxes indicate the deepness of the refinement level, from level l = 0 outside of the boxes to level l = lmax in the most inner box.

In the text

	Fig. 4 Evolution of angular momentum for all simulations with Rdisk = 100 and up to t = 1.5 Tdisk. The maximum deviations from the initial state reach up to 2%.
In the text

	Fig. 5 Face-on column density profiles for simulations with gc = 643 and t = 1.5 Tdisk. Top row: simulations with Qinit = 3, bottom row: simulations with Qinit = 1. The maximum refinement level lmax is increasing from left to right from 2−4.
In the text

	Fig. 6 Edge-on volume density profiles for simulations with gc = 643 and t = 1.5 Tdisk. Top row: simulations with Qinit = 3, bottom row: simulations with Qinit = 1. The maximum refinement level lmax is increasing from left to right from 2−4.
In the text

Fig. 7 Radial profiles of column density Σ and rotational velocity vrot for all simulations with Rdisk = 100 and t = 1.5 Tdisk. The runs with Qinit = 3 approximately follow the initially chosen density distribution. The Qinit = 1 disks show peaks at radii with distinct clumps or spiral structure (compare Figs. 5 and 6). The velocity profiles show only minor deviations from the general trend, which are due to the turbulent velocity in the most unstable regions.

In the text

Fig. 8 Radial profiles of the temperature distribution and effective Toomre Q parameter for all simulations with Rdisk = 100 and t = 1.5 Tdisk. As the density distribution the temperature shows a spread-out for the inner and outer parts. This is much more pronounced for the simulations with Qinit = 1, whereas the runs with Qinit = 3 roughly follow the initial distribution. For Qinit = 3, the effective Q parameter stays well above the threshold for marginal stability. In comparison, the effective Q in the Qinit = 1 runs dips below the threshold at the positions with increased density (i.e., fragmenting regions, compare Fig. 7).

In the text

	Fig. 9 Number of lowest-level clumps for all simulations and t = 1.5 Tdisk. The trend goes for just a few clumps for low maximum resolutions to many clumps for high maximum resolution. Different combinations of gc and lmax but comparable resolution on the lowest level (i.e., green/cyan and blue/magenta) show approximately equal numbers of lowest-level clumps.
In the text

Fig. 10 Radial profiles at t = 1.5 Tdisk for column density Σ and rotational velocity vrot for all simulations with the highest resolution on the finest grid, normalized in radial direction. The simulation with Qinit = 3 is added for comparison. The density distributions for all radii show similarities. However, the simulation with r = 100 AU and Qinit = 1 shows its single peak at ~52 AU and the negative slope at the outer edge of the disk is much shallower than the slope of the other simulations with Qinit = 1. The velocity distributions for all radii show a similar trend. The simulation with r = 10 AU abruptly goes down in velocity. This behavior might be due to the very sharp transition from disk material (with approximately Keplerian velocity profile) to surrounding material (with approximately zero velocity).

In the text

Fig. 11 Radial profiles at t = 1.5 Tdisk for temperature T and effective Toomre Q parameter for all simulations with the highest resolution on the finest grid, normalized in radial direction. The simulation with Qinit = 3 is added for comparison. Again, we see a different behavior of the simulation with Rdisk = 100 AU and Qinit = 1, displaying a smaller increase in temperature because of its smoother transition zone. The dips where the Q parameter goes below 1 are consistent with the density peaks in Figs. 10 and 12.

In the text

	Fig. 12 Top row: density projections of all disks with Qinit = 1 and lmax = 4, t = 1.5 Tdisk. Bottom row: vertical density slices of all disks with Qinit = 1, lmax = 4 and t = 1.5 Tdisk.
In the text