Issue 
A&A
Volume 579, July 2015



Article Number  A32  
Number of page(s)  12  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201424528  
Published online  23 June 2015 
Modeling gravitational instabilities in selfgravitating protoplanetary disks with adaptive mesh refinement techniques
^{1}
Institute for Astronomy,
ETH Zürich, WolfgangPauliStrasse
27
8093
Zürich
Switzerland
email:
tim.lichtenberg@phys.ethz.ch
^{2}
Institute of Geophysics, ETH Zürich, Sonneggstrasse 5, 8092
Zürich,
Switzerland
^{3}
Departamento de Astronomía, Universidad de Concepción, Av. Esteban
Iturra s/n Barrio Universitario, 160C Casilla,
Chile
^{4}
Institut für Astrophysik, GeorgAugustUniversität
Göttingen, FriedrichHundPlatz
1, 37077
Göttingen,
Germany
Received: 4 July 2014
Accepted: 27 April 2015
The astonishing diversity in the observed planetary population requires theoretical efforts and advances in planet formation theories. The use of numerical approaches provides a method to tackle the weaknesses of current models and is an important tool to close gaps in poorly constrained areas such as the rapid formation of giant planets in highly evolved systems. So far, most numerical approaches make use of Lagrangianbased smoothedparticle hydrodynamics techniques or gridbased 2D axisymmetric simulations. We present a new global disk setup to model the first stages of giant planet formation via gravitational instabilities (GI) in 3D with the blockstructured adaptive mesh refinement (AMR) hydrodynamics code enzo. With this setup, we explore the potential impact of AMR techniques on the fragmentation and clumping due to largescale instabilities using different AMR configurations. Additionally, we seek to derive general resolution criteria for global simulations of selfgravitating disks of variable extent. We run a grid of simulations with varying AMR settings, including runs with a static grid for comparison. Additionally, we study the effects of varying the disk radius. The physical settings involve disks with R_{disk} = 10,100 and 300 AU, with a mass of M_{disk} ≈ 0.05 M_{⊙} and a central object of subsolar mass (M_{⋆} = 0.646 M_{⊙}). To validate our thermodynamical approach we include a set of simulations with a dynamically stable profile (Q_{init} = 3) and similar grid parameters. The development of fragmentation and the buildup of distinct clumps in the disk is strongly dependent on the chosen AMR grid settings. By combining our findings from the resolution and parameter studies we find a general lower limit criterion to be able to resolve GI induced fragmentation features and distinct clumps, which induce turbulence in the disk and seed giant planet formation. Irrespective of the physical extension of the disk, topologically disconnected clump features are only resolved if the fragmentationactive zone of the disk is resolved with at least 100 cells. The latter corresponds to a minimum requirement for all global disk setups. Our simulations illustrate the capabilities of AMRbased modeling techniques for planet formation simulations and underline the importance of balanced refinement settings to reproduce fragmenting structures. The clumps in our models are migrating inward and are eventually destroyed because of tidal disruptions, reflecting the absence of radiative feedback from the central star, which may stabilize the clumps on larger scales. We expect that the inclusion of additional physics, like a radiation transport mechanism and the formation of sink particles, will provide a detailed framework to study the formation of planets via gravitational instabilities in a global disk view.
Key words: protoplanetary disks / hydrodynamics / instabilities / planets and satellites: formation / methods: numerical
© ESO, 2015
1. Introduction
Since the dawn of observational exoplanet detection in recent decades the evergrowing 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 bottomup 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 metalpoor 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 precataclysmic 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 secondgeneration 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 selfgravitating disks. In the NN Serpentis scenario, the age of the white dwarf of 10^{6} 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 smoothedparticle 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 selfgravitating 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 builtin 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 gridbased 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 builtin singularity (i.e., the center of origin) and handle all computed regions equally. The latter means there is no builtin refinement in inner regions, which might avoid symmetrydependent 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 blockstructured AMR hybrid code enzo^{1} (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 nonideal magnetohydrodynamics, and Nbody dynamics, including selfgravity for fluids and gas chemistry. The fluid flow in enzo is evolved using a finitevolume discretization, the ideal gas dynamics calculations are treated with the piecewise parabolic method (PPM), which is a higherorder Godunov scheme (Bryan et al. 1995) with an accurate piecewise parabolic interpolation and a nonlinear Riemann solver for shock conditions. The solver is thirdorder accurate in space and second order accurate in time for fixed time stepping and formally secondorder 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 in which x and v refer to the onedimensional 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 onedimensional 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, onedimensional 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 (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 (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 selfgravity of the gas. To achieve this, Poisson’s equation (6)with the gravitational potential φ, the gravity constant G, and the gas density ρ, is solved in a multistep process for all cells. The cloudincell (CIC) interpolation method from Hockney & Eastwood (1988) is used to approximate the gas distribution as a spatiallydiscretized 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 (7)with the thermal energy density e governed by the equation of state (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)^{3}. 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., (9)with an arbitrary field quantity q, and with the velocity perpendicular to the boundary direction reversed (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 socalled 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 (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 superLagrangian (ϵ_{l}< 0) or subLagrangian (ϵ_{l}> 0).
The second approach refines the Jeans length (12)by a fixed number of cells (based on Truelove et al. 1998) if (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 powerlaw profile (Mestel 1963) (14)with r_{out} denoting the outer radius of the disk and Σ_{0} 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 (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 (16)Here, z denotes the distance to the midplane of the disk and H the scale height, which is evaluated as (Lodato 2008) (17)with the scale height for selfgravitating disks (18)and the scale height for Keplerian disks (19)In the equation above, the sound speed is given via the NewtonLaplace equation (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 (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 selfgravity, the azimuthal velocity satisfies (Lodato 2008) (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) (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 (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.
Overview of the settings of all performed simulation runs.
Fig. 1 Radial profiles for initial column density Σ and rotational velocity v_{rot} for simulations with R_{disk} = 100. The parameter g_{c} indicates the resolution of the coarsest grid, l_{max} the deepest level of the resolution by dynamic refinement. The continous lines show simulations with Q_{init} = 1, the dashed lines simulations with Q_{init} = 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 Q_{init} = 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. 

Open with DEXTER 
Fig. 2 Radial profiles for initial temperature T and effective Toomre Q parameter for simulations with R_{disk} = 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. 

Open with DEXTER 
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^{3}, while analogous simulations with the same effective resolution can be identified for the runs with g_{c} = 128^{3}. In the following we discuss and explain the differences in detail.
Fig. 3 Initial vertical density profile Σ, comparing disks with R_{disk} = 100 AU and Q_{init} = 3 (top), Q_{init} = 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 Q_{init} = 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 = l_{max} in the most inner box. 

Open with DEXTER 
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 R_{disk} = 100 and up to t = 1.5 T_{disk}. The maximum deviations from the initial state reach up to 2%. 

Open with DEXTER 
Fig. 5 Faceon column density profiles for simulations with g_{c} = 64^{3} and t = 1.5 T_{disk}. Top row: simulations with Q_{init} = 3, bottom row: simulations with Q_{init} = 1. The maximum refinement level l_{max} is increasing from left to right from 2−4. 

Open with DEXTER 
Fig. 6 Edgeon volume density profiles for simulations with g_{c} = 64^{3} and t = 1.5 T_{disk}. Top row: simulations with Q_{init} = 3, bottom row: simulations with Q_{init} = 1. The maximum refinement level l_{max} is increasing from left to right from 2−4. 

Open with DEXTER 
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 v_{rot} for all simulations with R_{disk} = 100 and t = 1.5 T_{disk}. The runs with Q_{init} = 3 approximately follow the initially chosen density distribution. The Q_{init} = 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. 

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

Open with DEXTER 
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 faceon 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 largescale 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 edgeon 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 edgeon slices of the disks with the faceon projections demonstrates the capabilities of threedimensional 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 selfgravitating 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 .
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 lowestlevel clumps for all simulations and t = 1.5 T_{disk}. The trend goes for just a few clumps for low maximum resolutions to many clumps for high maximum resolution. Different combinations of g_{c} and l_{max} but comparable resolution on the lowest level (i.e., green/cyan and blue/magenta) show approximately equal numbers of lowestlevel clumps. 

Open with DEXTER 
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 farend 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 T_{disk} for column density Σ and rotational velocity v_{rot} for all simulations with the highest resolution on the finest grid, normalized in radial direction. The simulation with Q_{init} = 3 is added for comparison. The density distributions for all radii show similarities. However, the simulation with r = 100 AU and Q_{init} = 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 Q_{init} = 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). 

Open with DEXTER 
Fig. 11 Radial profiles at t = 1.5 T_{disk} 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 Q_{init} = 3 is added for comparison. Again, we see a different behavior of the simulation with R_{disk} = 100 AU and Q_{init} = 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. 

Open with DEXTER 
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 faceon views provide detailed insight into the fragmented disk structures. The edgeon 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 selfgravitating, 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 Q_{init} = 1 and l_{max} = 4, t = 1.5 T_{disk}. Bottom row: vertical density slices of all disks with Q_{init} = 1, l_{max} = 4 and t = 1.5 T_{disk}. 

Open with DEXTER 
The Q_{init} = 1 disks display the onset of largescale 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 selfgravitating 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 (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 (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 gridbased 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 gravitydriven 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 selfgravitating 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 longlasting 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 clumpdisk 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 longlasting 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^{6} 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 selfgravitating protoplanetary disks with AMR simulations. We expect that this technique can provide additional insight into the formation of massive selfgravitating 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 gridbased 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, ytproject.org), which is based on the matplotlib Python library (Hunter 2007, matplotlib.org).
References
 Agertz, O., Moore, B., Stadel, J., et al. 2007, MNRAS, 380, 963 [NASA ADS] [CrossRef] [Google Scholar]
 Avenhaus, H., Quanz, S. P., Schmid, H. M., et al. 2014, ApJ, 781, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., Meru, F., & Paardekooper, S.J. 2011, MNRAS, 416, 1971 [NASA ADS] [CrossRef] [Google Scholar]
 Benz, W., Ida, S., Alibert, Y., Lin, D., & Mordasini, C. 2014, Protostars and Planets VI, 691 [Google Scholar]
 Berger, M. J., & Colella, P. 1989, J. Comput. Phys., 82, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Beuermann, K., Dreizler, S., & Hessman, F. V. 2013, A&A, 555, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bodenheimer, P., Hubickyj, O., & Lissauer, J. J. 2000, Icarus, 143, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Boley, A. C. 2009, ApJ, 695, L53 [NASA ADS] [CrossRef] [Google Scholar]
 Boley, A. C., Durisen, R. H., Nordlund, Å., & Lord, J. 2007, ApJ, 665, 1254 [NASA ADS] [CrossRef] [Google Scholar]
 Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Boss, A. P. 1997, Science, 276, 1836 [NASA ADS] [CrossRef] [Google Scholar]
 Boss, A. P. 2003, ApJ, 599, 577 [NASA ADS] [CrossRef] [Google Scholar]
 Bryan, G. L., Norman, M. L., Stone, J. M., Cen, R., & Ostriker, J. P. 1995, Comput. Phys. Commun., 89, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R., & Klessen, R. S. 2011, ApJ, 731, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Garufi, A., Quanz, S. P., Schmid, H. M., et al. 2014, A&A, 568, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gawryszczak, A. J., Mayer, L., Boley, A., & Tasker, E. 2010, in EAS Pub. Ser. 42, eds. K. Gożdziewski, A. Niedzielski, & J. Schneider, 267 [Google Scholar]
 Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Haefner, R., Fiedler, A., Butler, K., & Barwig, H. 2004, A&A, 428, 181 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Helled, R., Bodenheimer, P., Podolak, M., et al. 2014, Protostars and Planets VI, 643 [Google Scholar]
 Hessman, F. V., Beuermann, K., Dreizler, S., et al. 2011, in AIP Conf. Ser. 1331, eds. S. Schuh, H. Drechsel, & U. Heber, 281 [Google Scholar]
 Hockney, R. W., & Eastwood, J. W. 1988, Computer simulation using particles (Bristol: Hilger) [Google Scholar]
 Hopkins, P. F. 2013, MNRAS, 428, 2840 [NASA ADS] [CrossRef] [Google Scholar]
 Horner, J., Wittenmyer, R. A., Hinse, T. C., & Tinney, C. G. 2012, MNRAS, 425, 749 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Inayoshi, K., & Haiman, Z. 2014, MNRAS, 445, 1549 [NASA ADS] [CrossRef] [Google Scholar]
 Johnson, T. V., Mousis, O., Lunine, J. I., & Madhusudhan, N. 2012, in EGU General Assembly Conference Abstracts, 14, eds. A. Abbasi, & N. Giesen, 11532 [Google Scholar]
 Latif, M. A., & Schleicher, D. R. G. 2015, MNRAS, 449, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013a, MNRAS, 430, 588 [NASA ADS] [CrossRef] [Google Scholar]
 Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013b, MNRAS, 432, 668 [NASA ADS] [CrossRef] [Google Scholar]
 Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. C. 2013c, MNRAS, 436, 2989 [NASA ADS] [CrossRef] [Google Scholar]
 Lemaster, M. N., & Stone, J. M. 2009, ApJ, 691, 1092 [NASA ADS] [CrossRef] [Google Scholar]
 Lodato, G. 2008, New A Rev., 52, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Marsh, T. R., Parsons, S. G., Bours, M. C. P., et al. 2014, MNRAS, 437, 475 [NASA ADS] [CrossRef] [Google Scholar]
 Mayer, L., & Gawryszczak, A. J. 2008, in Extreme Solar Systems, eds. D. Fischer, F. A. Rasio, S. E. Thorsett, & A. Wolszczan, ASP Conf. Ser., 398, 243 [Google Scholar]
 Mayer, L., Quinn, T., Wadsley, J., & Stadel, J. 2003, in Scientific Frontiers in Research on Extrasolar Planets, eds. D. Deming, & S. Seager, ASP Conf. Ser., 294, 281 [Google Scholar]
 Mayer, L., Lufkin, G., Quinn, T., & Wadsley, J. 2007, ApJ, 661, L77 [NASA ADS] [CrossRef] [Google Scholar]
 Meru, F., & Bate, M. R. 2011, MNRAS, 410, 559 [NASA ADS] [CrossRef] [Google Scholar]
 Meru, F., Quanz, S. P., Reggiani, M., Baruteau, C., & Pineda, J. E. 2014, ApJ, submitted [arXiv:1411.5366] [Google Scholar]
 Mestel, L. 1963, MNRAS, 126, 553 [NASA ADS] [CrossRef] [Google Scholar]
 Michael, S., Durisen, R. H., & Boley, A. C. 2011, ApJ, 737, L42 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., Lunine, J. I., O’Brien, D. P., Raymond, S. N., & Walsh, K. J. 2012, Ann. Rev. Earth Planet. Sci., 40, 251 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., Szulagyi, J., Crida, A., et al. 2013, in AAS/Division for Planetary Sciences Meeting Abstracts, 45, 510.05 [Google Scholar]
 Nelson, A. F. 2006, MNRAS, 373, 1039 [NASA ADS] [CrossRef] [Google Scholar]
 O’Shea, B. W., Bryan, G., Bordner, J., et al. 2004, ArXiv eprints [arXiv:astroph/0403044] [Google Scholar]
 Paardekooper, S.J., & Mellema, G. 2006, A&A, 450, 1203 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Parsons, S. G., Marsh, T. R., Bours, M. C. P., et al. 2014, MNRAS, 438, L91 [NASA ADS] [CrossRef] [Google Scholar]
 Price, D. J. 2008, J. Comput. Phys., 227, 10040 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Price, D. J. 2012, J. Comput. Phys., 231, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Price, D. J., & Federrath, C. 2010, MNRAS, 406, 1659 [NASA ADS] [Google Scholar]
 Raymond, S. N., Barnes, R., & Kaib, N. A. 2006, ApJ, 644, 1223 [NASA ADS] [CrossRef] [Google Scholar]
 Raymond, S. N., Kokubo, E., Morbidelli, A., Morishima, R., & Walsh, K. J. 2014, Protostars and Planets VI, 914, 595 [NASA ADS] [Google Scholar]
 Read, J. I., & Hayfield, T. 2012, MNRAS, 422, 3037 [NASA ADS] [CrossRef] [Google Scholar]
 Schleicher, D. R. G., & Dreizler, S. 2014, A&A, 563, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schleicher, D. R. G., Bovino, S., Latif, M. A., Ferrara, A., & Grassi, T. 2015, A&A, submitted [arXiv:1504.06296] [Google Scholar]
 Setiawan, J., Klement, R. J., Henning, T., et al. 2010, Science, 330, 1642 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, B. D., Silvia, D. W., Turk, M. J., et al. 2010, in AIP Conf. Ser. 1294, eds. D. J. Whalen, V. Bromm, & N. Yoshida, 110 [Google Scholar]
 Toomre, A. 1964, ApJ, 139, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1998, ApJ, 495, 821 [NASA ADS] [CrossRef] [Google Scholar]
 Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Völschow, M., Banerjee, R., & Hessman, F. V. 2014, A&A, 562, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wise, J. H., & Abel, T. 2011, MNRAS, 414, 3458 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., Hartmann, L., Nelson, R. P., & Gammie, C. F. 2012, ApJ, 746, 110 [NASA ADS] [CrossRef] [Google Scholar]
 Zorotovic, M., & Schreiber, M. R. 2013, A&A, 549, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
Fig. 1 Radial profiles for initial column density Σ and rotational velocity v_{rot} for simulations with R_{disk} = 100. The parameter g_{c} indicates the resolution of the coarsest grid, l_{max} the deepest level of the resolution by dynamic refinement. The continous lines show simulations with Q_{init} = 1, the dashed lines simulations with Q_{init} = 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 Q_{init} = 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. 

Open with DEXTER  
In the text 
Fig. 2 Radial profiles for initial temperature T and effective Toomre Q parameter for simulations with R_{disk} = 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. 

Open with DEXTER  
In the text 
Fig. 3 Initial vertical density profile Σ, comparing disks with R_{disk} = 100 AU and Q_{init} = 3 (top), Q_{init} = 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 Q_{init} = 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 = l_{max} in the most inner box. 

Open with DEXTER  
In the text 
Fig. 4 Evolution of angular momentum for all simulations with R_{disk} = 100 and up to t = 1.5 T_{disk}. The maximum deviations from the initial state reach up to 2%. 

Open with DEXTER  
In the text 
Fig. 5 Faceon column density profiles for simulations with g_{c} = 64^{3} and t = 1.5 T_{disk}. Top row: simulations with Q_{init} = 3, bottom row: simulations with Q_{init} = 1. The maximum refinement level l_{max} is increasing from left to right from 2−4. 

Open with DEXTER  
In the text 
Fig. 6 Edgeon volume density profiles for simulations with g_{c} = 64^{3} and t = 1.5 T_{disk}. Top row: simulations with Q_{init} = 3, bottom row: simulations with Q_{init} = 1. The maximum refinement level l_{max} is increasing from left to right from 2−4. 

Open with DEXTER  
In the text 
Fig. 7 Radial profiles of column density Σ and rotational velocity v_{rot} for all simulations with R_{disk} = 100 and t = 1.5 T_{disk}. The runs with Q_{init} = 3 approximately follow the initially chosen density distribution. The Q_{init} = 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. 

Open with DEXTER  
In the text 
Fig. 8 Radial profiles of the temperature distribution and effective Toomre Q parameter for all simulations with R_{disk} = 100 and t = 1.5 T_{disk}. As the density distribution the temperature shows a spreadout for the inner and outer parts. This is much more pronounced for the simulations with Q_{init} = 1, whereas the runs with Q_{init} = 3 roughly follow the initial distribution. For Q_{init} = 3, the effective Q parameter stays well above the threshold for marginal stability. In comparison, the effective Q in the Q_{init} = 1 runs dips below the threshold at the positions with increased density (i.e., fragmenting regions, compare Fig. 7). 

Open with DEXTER  
In the text 
Fig. 9 Number of lowestlevel clumps for all simulations and t = 1.5 T_{disk}. The trend goes for just a few clumps for low maximum resolutions to many clumps for high maximum resolution. Different combinations of g_{c} and l_{max} but comparable resolution on the lowest level (i.e., green/cyan and blue/magenta) show approximately equal numbers of lowestlevel clumps. 

Open with DEXTER  
In the text 
Fig. 10 Radial profiles at t = 1.5 T_{disk} for column density Σ and rotational velocity v_{rot} for all simulations with the highest resolution on the finest grid, normalized in radial direction. The simulation with Q_{init} = 3 is added for comparison. The density distributions for all radii show similarities. However, the simulation with r = 100 AU and Q_{init} = 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 Q_{init} = 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). 

Open with DEXTER  
In the text 
Fig. 11 Radial profiles at t = 1.5 T_{disk} 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 Q_{init} = 3 is added for comparison. Again, we see a different behavior of the simulation with R_{disk} = 100 AU and Q_{init} = 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. 

Open with DEXTER  
In the text 
Fig. 12 Top row: density projections of all disks with Q_{init} = 1 and l_{max} = 4, t = 1.5 T_{disk}. Bottom row: vertical density slices of all disks with Q_{init} = 1, l_{max} = 4 and t = 1.5 T_{disk}. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.