Convergence of the critical cooling rate for protoplanetary disk fragmentation achieved; the key role of numerical dissipation of angular momentum
DDraft version July 22, 2018
Typeset using L A TEX twocolumn style in AASTeX61
CONVERGENCE OF THE CRITICAL COOLING RATE FOR PROTOPLANETARY DISK FRAGMENTATIONACHIEVED; THE KEY ROLE OF NUMERICAL DISSIPATION OF ANGULAR MOMENTUM
Hongping Deng, Lucio Mayer, and Farzana Meru Center for Theoretical Astrophysics and Cosmology, Institute for Computational Science, University of Zurich, Winterthurerstrasse 190,8057 Zurich, Switzerland Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK (Received July 1, 2016; Revised September 27, 2016; Accepted July 22, 2018)
Submitted to ApJABSTRACTWe carry out simulations of gravitationally unstable disks using smoothed particle hydrodynamics(SPH) and thenovel Lagrangian meshless finite mass (MFM) scheme in the GIZMO code (Hopkins 2015). Our aim is to understandthe cause of the non-convergence of the cooling boundary for fragmentation reported in the literature. We runSPH simulations with two different artificial viscosity implementations, and compare them with MFM, which doesnot employ any artificial viscosity. With MFM we demonstrate convergence of the critical cooling time scale forfragmentation at β crit ≈ . . Non-convergence persists in SPH codes, although it is significantly mitigated withschemes having reduced artificial viscosity such as inviscid SPH (ISPH) (Cullen & Dehnen 2010). We show how thenon-convergence problem is caused by artificial fragmentation triggered by excessive dissipation of angular momentumin domains with large velocity derivatives. With increased resolution such domains become more prominent. Vorticitylags behind density due to numerical viscous dissipation in these regions, promoting collapse with longer cooling times.Such effect is shown to be dominant over the competing tendency of artificial viscosity to diminish with increasingresolution. When the initial conditions are first relaxed for several orbits, the flow is more regular, with lower shearand vorticity in non-axisymmetric regions, aiding convergence. Yet MFM is the only method that converges exactly.Our findings are of general interest as numerical dissipation via artificial viscosity or advection errors can also occurin grid-based codes. Indeed for the FARGO code values of β crit significantly higher than our converged estimate havebeen reported in the literature. Finally, we discuss implications for giant planet formation via disk instability. Keywords:
Gravitational instability —accretion, accretion disks — fragmentation —numerical viscos-ity
Corresponding author: Hongping [email protected] a r X i v : . [ a s t r o - ph . E P ] N ov Deng et al. INTRODUCTIONMassive rotating gas disks can become unstable dueto their own self-gravity. Using linear perturbation the-ory applied to an axisymmetric, razor-thin disk, Toomre(1964) showed that the exponential growth of local den-sity perturbations is controlled by the parameter Q ≡ c s κπG Σ (1)where c s is the sound speed in the disk, κ is theepicyclic frequency, which equals the angular frequencyΩ for a Keplerian disk, and Σ is the mass surface density( G is the gravitational constant). When Q < Q crit ≈ Q above thestability threshold. More specifically, whether or not Q will remain low enough for the instability to grow, andeventually lead to fragmentation, ultimately depends onthe ability of the gas to release the heat via radiativecooling. Gammie (2001) proposed to parameterize thelocal cooling time scale t cool as β = t cool Ω (2)where t cool = u ( dudt ) − (3) u is the specific internal energy. He showed that when β < β crit ≈ γ = 2.Physically this means cooling has to occur on a timescalecomparable to the local orbital time, which can be eas-ily understood as the spiral pattern that generates shockheating also evolves on the local orbital time. Using 3DSPH simulations, Rice et al. (2005) showed β crit ≈ − γ = 5 /
3. A similar value of β crit wasalso found by Mayer et al. (2005) using a different 3D SPH code while comparing the evolution of isolated andbinary protoplanetary disks. However, Meru & Bate(2011a) carried out similar SPH simulations to thosein Rice et al. (2005) and found non-convergence of thecritical cooling rate when the resolution of the simula-tions was increased. They found disks can fragment at ahigher β crit value in higher resolution simulations. Meru& Bate (2012) suggested that the non-convergent be-haviour was caused by extra heating due to artificial vis-cosity occurring at low resolution. At higher resolution,less heating is generated through artificial viscosity anda lower critical cooling rate (i.e. a larger β ) is capableof balancing the combined heating by gravitational in-stability and artificial viscosity. Earlier on Paardekooperet al. (2011) had run 2D simulations with the grid-basedcode FARGO, showing that the non-convergence prob-lem is not specific to SPH. By running a large set ofsimulations with both 3D SPH and the 2D grid-basedcode FARGO,Meru & Bate (2012) concluded that β crit should converge to an asymptotic value ≈ −
30. Thishas important implications since it suggests fragmenta-tion can occur even with relatively long cooling times,corresponding to almost an order of magnitude longerthan the local orbital time. If clumps produced by frag-mentation can later evolve into gas giant planets (see eg.Helled et al. (2014)), it follows that giant planet forma-tion via disk instability is more common than previouslythought.Paardekooper et al. (2011) suggested that the non-convergence is at least partly due to the emergence ofspecial locations in the disk, at the boundary betweenthe turbulent and the laminar region, which appearwhen the disk is still adjusting towards a well devel-oped gravito-turbulent state. In order to avoid this edgeeffect, Paardekooper (2012) carried out high resolution2D local shearing sheet FARGO simulations, but founda similar increase of β crit as resolution increased. Frag-mentation also exhibited a stochastic nature in his simu-lations. However, recently Klee et al. (2017) performeda numerical study arguing that stochastic fragmenta-tion is probably caused by oversteepening when per-forming slope limiting in grid-based codes. Moreover,Young & Clarke (2015) showed that the fragmentationboundary is strongly related to the gravitational soft-ening in 2D simulations. When the softening is compa-rable to the resolution scale β crit increases with resolu-tion, while when the softening is comparable to the diskscale height, β crit varies little. Young & Clarke (2015)also showed that clumps form from overdensities withlength-scales ∼ H . This is consistent with the earlierresult of Gammie (2001), who found no power at scalesmuch smaller than H , which were resolved in the simula- ragmentation of self-gravitating disks H .In this paper, we attempt to explore again the issueof non-convergence by comparing different Lagrangianhydro methods which have by design different numeri-cal dissipation. We are driven by the notion that non-convergence might be caused by numerical errors, per-haps associated with the boundary regions identified byPaardekooper et al. (2011). Our aim is to identify theexact source of non-convergence as well as a way to over-come the problem, in order to place self-gravitating disksimulations on a firm ground.We run simulations using exactly the same initialcondition of Meru & Bate (2012) with three differentLagrangian hydro methods. Specifically, we use thevanilla SPH implementation in the GIZMO code (Hop-kins 2015), which is equivalent to the original GAD-GET3 code ((Springel 2005) with standard Monaghanartificial viscosity and Balsara switch (Balsara 1995), amodified version of the same code where we implementedthe artificial viscosity scheme of Cullen & Dehnen (2010)(inviscid SPH, hereafter ISPH), and the new Lagrangianmeshless finite mass(MFM) method in GIZMO (Hop-kins 2015),which does not need any artificial viscosity.Finally, a smaller set of additional simulations exploringrecent variants of the SPH method are presented at theend of the paper.In section 2, we describe the setup of the simulations.In section 3, we present the main results. We discuss theinterpretation and implications of our results in section 4,and finally draw our conclusions in section 5. Results ofancillary numerical tests are also presented in AppendixA and Appendix B. SIMULATIONSIn order to compare with previous work, we use thesame disk model as in Meru & Bate (2012). A 0 . M (cid:12) disk spans a radial range, 0 . < R <
25 AU, surround-ing a solar mass star. The initial surface mass densityand temperature profiles are Σ ∝ R − and T ∝ R − / ,respectively. The temperature is normalized so that atthe outer edge of the disk the Toomre Q parameterequals 2. We use an adiabatic gas equation of state with γ = 5 /
3. We note that the problem is essentially scalefree and we 1 solar mass, 1 AU, and 1yr as the mass,length and time unit respectively in our calculation. Thesimulations are run with different cooling rates as wellas three different resolutions. The low resolution model(LR) comprises 250,000 particles, while 2 million and 16million particles are used, respectively, in the high reso-lution (HR) and ultra-high resolution model(UHR). The central star is modeled with a sink particle with sink ra-dius of 0.25 AU. Gas particles reaching the sink radiusare deleted and their mass and momentum are addedto it to ensure mass and momentum conservation. Weuse adaptive gravitational softening(Price & Monaghan2007).The three disk models are run with the three differ-ent hydro-methods, namely standard vanilla SPH as inMeru & Bate (2012) (hereafter TSPH), inviscid SPH(ISPH), which uses the implementation of artificial vis-cosity by Cullen & Dehnen (Cullen & Dehnen 2010), andthe MFM method. The first two methods differ only inthe way artificial viscosity is implemented. In vanillaSPH we solve the fluid equations using the standarddensity-entropy formulation(Springel 2005) and stan-dard Monaghan artificial viscosity with Balsara switch(Balsara 1995), designed to reduce viscosity in flows withhigh vorticity. In ISPH we use the same formulation ofthe hydro equations but we use the Cullen & Dehnenviscosity which uses the time derivative of the velocitydivergence as shock indicator in order to eliminate nu-merical viscous dissipation away from shocks. In thisformulation, particles have individual adjustable viscos-ity coefficients between an α max and an α min value,andare further assigned a coefficient l that sets the decaytime of viscosity away from shocks.Therefore the latterformulation requires to set three parameters instead ofthe conventional two coefficients α sph and β sph in thestandard Monaghan viscosity. While previous artificialviscosity schemes had appeared that used switches act-ing on individual particles Morris & Monaghan (1997),including with similar parameterization (Rosswog et al.2000), the ISPH method is a further improvement to-wards a truly inviscid method owing to its shock track-ing procedure. The MFM method is instead an entirelydifferent method to solve the hydrodynamical equations,but shares with SPH the Lagrangian approach and theuse of particles as tracers of the fluid. Most importantly,it does not need any artificial viscosity in order to gen-erate stable solutions, even in high Mach number flows,hence it is by nature inviscid. While we defer the readerto Hopkins (2015) for a detailed description of the codehere we recall that MFM uses the SPH-like kernel ofa particle distribution to construct a volume partitionof the domain. The fluid equations are then solved onthe unstructured mesh thus generated by the partitionusing a Riemann solver (HLLC in this paper, see Hop-kins (2015)). As the volume elements are constructedfrom the particles the smoothing length h of the assignedkernel function still defines the fundamental resolutionlength as in SPH. The method conserves mass, momen-tum and angular momentum by design as an SPH code, Deng et al. but this is even more strictly true for angular momen-tum relative to SPH since there is no added artificialviscosity, as pointed out and demonstrated with a seriesof tests in Hopkins (2015). Being an integral formu-lation, in this sense an hybrid with finite-volume grid-based methods, it should also not suffer from the er-rors associated with poor estimates of gradients properof SPH. Gravity is solved using the same tree code inall three methods, inherited from the progenitor GAD-GET3 code.For all codes we employ two different techniques tostart the simulations. In the first one we start directlyfrom the initial disk model without any relaxation phase,while in the second one we first relax the system by run-ning the disk for 4 outer rotation periods(ORPs) withrelatively long cooling time, using β = 12. By relax-ing the initial condition in the second case we allow thedevelopment of a gravoturbulent state, with sustainedspiral structure but no fragmentation, and then switchto the desired β value. The latter setup, which we in-dicate with TIC(turbulent initial condition) resemblesthat adopted by Paardekooper et al. (2011). It avoidsthe transition from a smooth disk to a turbulent disk asin figure 3. We will show this transition can lead to nu-merical fragmentation, as also found by Paardekooperet al. (2011). For both initial setups, once the desiredcooling timescale is turned on, simulations are run eitherfor at least 6 ORPs, or until the disk fragments.In order to facilitate comparisons we adopt standardvalues of artificial viscosity parameters in both TSPHand ISPH runs. For TSPH runs we use the Monaghanviscosity with α sph = 1 and β sph = 2, which is thedefault choice for flows in which high Mach numbers canarise (Hernquist & Katz 1989), and has been advocatedalso by Mayer et al. (2004) and Meru & Bate (2012). ForISPH we set α max = 2, α min = 0 and l = 0 . β . RESULTS3.1.
Overview
Table 1 summarizes our results on fragmentation. Ithighlights three main results borne out of our suite ofruns. First, the TSPH runs confirm the dependence
Table 1.
Simulations fragmented are marked with F, otherwiseare marked with NF. The boundary between fragmentation andnon-fragmentation is marked with green color.Simulation β =3 3.5 4 5 6 8 BoundaryTSPH-LR F F NF 6-8ISPH-LR F NF 5-6MFM-LR F F NF NF 3.5-4TSPH-HR F > of the critical β for fragmentation on resolution. Sec-ond, we can infer that non-convergence with resolutionis strongly mitigated by methods that reduce numeri-cal dissipation in shear flows. This includes improvedSPH artificial viscosity schemes designed to reduce it inshear flows, as in the case of ISPH, but is even more evi-dent for a new hydro schemes with no explicit numericalviscosity, such as MFM. Third, exact convergence is ob-tained only with MFM but demands the to start froma relaxed, gravoturbulent initial setup, confirming theclaim of Paardekooper et al. (2011) (see Table 1, “TIC”runs). Without relaxation a marginal increase of thecritical β with increasing resolution persists. In section3.2 we will provide an explanation of these main find-ings.At low resolution the fragmentation boundary inTSPH is found to lie between β = 6 and β = 8,which is slightly higher than the fragmentation bound-ary β = 5 . − . α sph = 0 . β sph = 0 .
2. If we compare instead with the subsetof their runs that adopted α sph = 1 and β = 2, consis-tent with ours, our results are in substantial agreement.Note that Meru & Bate (2012) show that using lowvalues of the coefficients α sph and β sph can result, coun-terintuitively, in high dissipation resulting from randomparticle motion. The noisier shear flow resulting in SPHsimulations with lowered coefficients in the standard ragmentation of self-gravitating disks α sph = 1 and β = 2 are the preferred values forhigh Mach number flows in a variety of astrophysicalregimes, and have been advocated for both galactic andprotoplanetary disk simulations (Mayer et al. (2004);Kaufmann et al. (2007)).Figure 1 shows runs using the three hydro implemen-tations for β values near critical. The markedly differ-ent behaviours at equivalent high resolution is evident.MFM needs β = 4 or lower to fragment,depending onthe initial conditions setup, while TSPH fragments fora β twice as large, and ISPH falls in between. We re-mind the reader that between TSPH and ISPH the onlydifference is the implementation of artificial viscosity.The Figure also shows how MFM simulations startingfrom turbulent initial conditions (TIC) are resilient tofragmentation even with very fast cooling at the highestresolution (UHR, see bottom right panel). The TICsimulations also exhibit a more flocculent spiral patterndominated by lower order modes relative to the non-relaxed simulations, and this is true for all codes.3.2. Numerical dissipation in MFM runs
When comparing TSPH and ISPH with MFM differ-ences in the hydro implementation are more subtle. In-deed MFM does not have explicit numerical viscosity,rather it solves the hydro equations analogously to a fi-nite volume method using a Riemann solver on volumeelements mapped from the particle distribution. Therecould be some residual numerical dissipation associatedwith the non-analytical Riemann problem solution, butat the same time MFM does not advect fluid elementsthrough a grid, a common source of numerical dissipa-tion in finite volume grid-based methods. In the Keple-rian ring test shown in Hopkins (2015) MFM conservesangular momentum better than all the other methodsit was compared to, including SPH with a modified hy-dro force (pressure-entropy,hereafter PSPH) that helpsto remove unwanted numerical effects, such as artificialtension, at fluid interfaces. This would suggest MFMis inherently less viscous than all variants of SPH , irre-spective of the artificial viscosity scheme adopted. How-ever in the latter test the disk is 2D, adiabatic and nonself-gravitating, which is very different from the con-figuration we are studying here. Therefore, in orderto reassess that MFM is indeed less viscous in a rele-vant albeit simple configuration we run two different testproblems.First we consider the rotating uniform isother-mal sphere collapse test in Boss & Bodenheimer (1979). The results are shown in appendix A. We find that an- gular momentum is transported slightly faster outwardsin TSPH runs as expected from stronger numerical vis-cous transport, albeit differences are small and overallthe angular momentum profiles are comparable betweenthe two codes at late times. Furthermore, in appendixB, we use an accreting sink particle in a shearing sheetconfiguration. Here the flow responds more abruptly tothe strong gravitational pull of the sink, and differencesbetween the three methods emerge more clearly, estab-lishing that MFM has better angular momentum con-servation in flows with large velocity gradients. Whilewe find that ISPH, as expected, has a lower numericalviscosity, the accretion rate, as measured by the growthof the central density, is comparable to TSPH. Instead,MFM leads to consistently smaller density increase, andthe differences remains over time.In addition to angular momentum dissipation, the evo-lution of the internal energy is another important proxyfor how strong the numerical dissipation by artificial vis-cosity is, especially in spiral shocks. Figure 2 shows theinternal energy measured at the disk mid-plane after1.07 ORPs in a non-fragmenting hi-res run (for whichwe used a large enough β = 8 to avoid fragmentation).A non-fragmenting run is a conservative choice for ourpurpose since spiral shocks are weaker. The Figure high-lights how the MFM run exhibits the cooler profile ev-erywhere because, at an equivalent cooling rate, it suf-fers less from artificial viscous heating, even in the regionnear the inner boundary where SPH is notoriously prob-lematic in handling the rapidly increasing shear (Mayeret al. 2004; Kaufmann et al. 2007). The TSPH and ISPHsimulations are almost indistinguishable. The similarinternal energy peak occurs slightly above 20 AU in allruns because that is where a strong ring-like overdensitydevelops (similar to those in the upper panel of Figure3), suggesting that where the flow is intrinsically morecompressive in all codes the dominant heating comesfrom PdV work rather than numerical viscosity.Having now assessed that MFM is the less viscousamong the methods considered, we can turn again tothe interpretation of the results in Figure 1. If, as sug-gested in Meru & Bate (2012), reduced numerical vis-cous heating as the resolution is increased is the sourceof non-convergence, one would expect that MFM runsshould be more prone rather than less prone to frag-mentation relative to TSPH and ISPH runs. This isbecause gas can cool more effectively at a given β ,as wehave just seen by comparing the internal energy evo-lution. Instead the opposite trend is seen – MFM isless prone to fragmentation and its behaviour is nearlyconvergent with resolution (exactly convergent in TICruns). Furthermore, ISPH and TSPH appear to be Deng et al.
Figure 1.
Surface mass density ([ g/cm ]) in logarithmic scale. The box size is 50 code units (axes in code units are shownin red at the bottom left corner).A subset of the simulations for the three different methods and with β near critical is shownto highlight that MFM and ISPH do not fragment with much lower values of β relative to TSPH at comparable resolution(compare bottom panel to top panel). A TIC run at even higher resolution (UHR) is also shown in the bottom right panel,which also does not fragment despite a very low β (see Table 1). From the first to the second row, left to right, snapshots aretaken at 3.7, 0.8, 5.8, 6, 6, 6 ORPs. Note the different times were chosen for the snapshots in the top panel according to whenfragmentation begins, which is much earlier than the time chosen in the bottom panel. The non-fragmenting runs have reacheda self-regulated state even before 6ORPs, which is thus a reliable reference time. Figure 2.
Azimuthally averaged radial internal energy pro-file in the disk mid-plane, in code units and logarithmic scaleat 1.07 ORPs. The HR simulations run with different hydromethods are shown. The MFM run has a much cooler diskdue to less numerical viscosity and thus reduced numericalheating (see Section 3). equivalent when we inspect the internal energy evolution in Figure 2, yet ISPH is less prone to fragmentation. Allthis suggests that, in order to understand the nature ofour results, among the effects associated with numericalviscous dissipation, enhanced angular momentum trans-port, as opposed to enhanced heating, is key. We discussthis crucial point, as well as dependence on resolution,in the next section.3.3.
Numerical dissipation and theresolution-dependent flow properties
In standard SPH methods artificial viscosity is neededfor the stability of the flow near discontinuities. Withoutthat particle-interpenetration occurs and shocks cannotbe properly resolved. The conventional form of artifi-cial viscosity employs two terms (Monaghan & Gingold1983; Springel 2005); one term linear in the divergenceof the flow, controlled by the α sph coefficient, and oneterm with quadratic dependence on the flow divergence,controlled by β sph . The quadratic term is dominant inshocks. When β sph = 0, the linear term can be re- ragmentation of self-gravitating disks η = (1 / α sph κhc s ρ and a bulkviscosity ζ = (5 / η , where h is the SPH smoothinglength, which controls the resolution, and κ is a constantthat should be chosen according to the kernel (Cullen &Dehnen 2010). The shear component of the viscosity willalways be present,namely also in non-shocking regions,and can generate artificial angular momentum transportin shearing flows such as those inherent to astrophysi-cal disks, even in presence of damping factors dependenton the local vorticity such as the Balsara switch (Kauf-mann et al. 2007). Cullen & Dehnen (2010), buildingon previous work (Morris & Monaghan 1997; Rosswoget al. 2000) proposed a new artificial viscosity switchwith a high order shock indicator which results, in amuch more effective damping of shear viscosity awayfrom shocks. It is thus expected that, at fixed resolution,ISPH runs, which employ the latter viscosity scheme,should be less affected by numerical dissipation of angu-lar momentum. At the same time, we caution that oneexpects mild shocks to arise as the spiral pattern growsin amplitude and the disk becomes marginally unstable,which could then result in significant viscous dissipationeven in ISPH.Having shown in more than one way, in the previ-ous section, that MFM is a less viscous method, andknowing that ISPH is by construction less viscous thanTSPH, we have inferred that the different susceptibilityto fragmentation of the three type of runs must be some-how caused by the different degree of angular momen-tum dissipation. The next step is to address how angularmomentum dissipation affects fragmentation. Moreover,there is a related puzzling fact emerging from Table 1.This is that one would expect the results of hydro meth-ods whose main difference is the viscous dissipation ofangular momentum should eventually converge as theresolution is increased because,owing to the dependenceof shear viscosity on the smoothing length h just high-lighted, viscous dissipation should decrease as resolutionincreases. Instead the opposite is seen; β crit keeps in-creasing for TSPH as the resolution is increased, hencedeparting more and more from the β crit value deter-mined with MFM. However, we caution already at thispoint that this conventional way of reasoning is correctonly if no other property of the flow changes with in-creasing resolution. If, for some reason, velocity gra-dients become intrinsically stronger as resolution is in-creased, as one could imagine that convergence mightbe hard to achieve irrespective of resolution.Here we argue that the nature of the flow in gravita-tionally unstable disks does change as resolution is in-creased . It develops regions of stronger shear and vortic-ity as sharper flow features become resolved, and such regions then become sites of higher numerical viscousdissipation. The resolution-dependent nature of the flowalso explains why the way the initial setup is preparedmatters for the fragmentation outcome at varying reso-lution , as we illustrate below. The changing nature ofthe flow with resolution is well illustrated in figure 3.The simulations using unrelaxed initial conditions offerthe best tool to understand what happens.In these thespiral pattern originates in the inner disk and propa-gates outwards(Meru & Bate 2011b). We refer to thestage before spirals reach the outer edge of the disk asthe transition phase . The duration of such transitionphase depends on the cooling rate, but in general lastsabout 2 ORPs. By visual inspection it is clear that athigher resolution non-axisymmetric modes are better re-solved. This is the result of an improved gravitationalforce resolution, a point already emphasized in Mayeret al. (2004),Mayer & Gawryszczak (2008) for both SPHand adaptive mesh refinement (AMR) codes.We refer to the boundary zones between the turbu-lent, highly non-axisymmetric flow component gener-ated by gravitational instability and the backgroundsmooth,axisymmetric disk as the interface region . Asseen in the upper panels of figure 3, these are regionsroughly a few code units in size, and occur during thetransition phase (see the rectangular region in Figure 3).Afterwards the entire disk becomes non-axisymmetricand gravito-turbulent. These interface regions are rel-evant because they are the location of the largest localvelocity gradients, associated with both strong shear andvorticity.We will use the analysis of vorticity as a proxy for highvelocity derivatives in the flow as vorticity also carriesinformation on the local angular momentum, which weposit plays a central role. The azimuthally averaged vor-ticity profile in the mid-plane is plotted in figure 4 forMFM runs at different resolution. In all of them it is ev-ident that vorticity can be an order of magnitude higherthan the Keplerian flow vorticity. These vorticity peaksoccur near interface regions and more pronounced andmore frequent as resolution increases. The latter is aneffect of increased force resolution as higher overdensi-ties appear that induce higher local vorticity. Indeed, ifnumerical dissipation of angular momentum does not af-fect the properties of the flow one would expect densityand vorticity to remain correlated. Let us consider anapproximately spherical cusp of the flow located at theinterface, namely an overdense region of radius R c andenclosed M c , so that its mean density is ρ = M c /R c .When the cusp density is large enough a fluid element atits boundary will move with a (circular) velocity V c suchthat √ Gρ ∝ V c R c . In other words, the cusp dominates Deng et al.
Figure 3.
Surface mass density ([ g/cm ]) in logarithmic scale for a subset of the MFM simulations run with different resolutionand β close to critical (axes in code length units are shown in red at the bottom left corner). The bottom panel shows the TICruns. Resolution increases by nearly 3 orders of magnitude from left to right. All snapshots are taken at 1.07 ORPs. A cusp ishighlighted in the red box in the top right panel. Cusps form in the interface regions between spirals and outer smooth flow.They cause strong numerical dissipation and are absent in TIC runs (compare bottom and top panels). Figure 4.
Azimuthally averaged radial vorticity profile inthe disk mid-plane (the modulus of vorticity is in code units)measured at 1.07 ORPs for MFM simulations with differentresolution. The vorticity of a pure Keplerian flow is alsoshown for comparison the local gravitational field, even if collapse has not en-sued, so that locally the vorticity of the flow should becorrelated with local density, since, dimensionally, thevorticity ∇ × v ∝ V c R c ∼ √ Gρ ( v is the velocity vector of the gas flow) We can now turn back to our simula-tions and check whether or not vorticity and density inthe cuspy regions grow at the relative rates implied bythe proportionality just highlighted. First of all, we no-tice that the peak density at interface regions in MFMruns in Figure 3 is almost an order of magnitude largerin the UHR simulation relative to the LR simulation.As expected, vorticity is correspondingly higher (Figure5). Note that TSPH runs share the same trend. ISPHruns yield a similar vorticity map as TSPH because atinterface regions the flow has high velocity derivativeshence artificial viscosity is not suppressed by the Cullen& Dehnen switch.We cannot compare the vorticity in the MFM andSPH simulations in figure 5 directly because the diskis much cooler in MFM(see figure 2) so that the spiralpattern evolves faster than in the SPH simulations. Butwe can still compare the evolution of maximum densityand maximum vorticity in a relative sense. We focuson the transition phase of a set of fast cooling simu-lations with β = 4. We plot the peak volume densityand peak vorticity in the r > t = 2 ORPs the spiral pattern reaches the outer edge ragmentation of self-gravitating disks ∇ × v ∝ √ Gρ .Indeed at 2 code units it has grown by about an orderof magnitude while density as grown by nearly two or-ders if magnitude.Conversely, a comparable increase invorticity occurs in TSPH and ISPH runs at this timebut there density increases much more, by about fourorders of magnitude. We argue that this remarkablemismatch reflects strong numerical dissipation of angu-lar momentum near cusps,which damps vorticity andpromotes collapse..We can attempt to understand the nature of theproblem by considering resolution requirements in self-gravitating disks. In order to properly follow the frag-mentation, the most unstable Toomre wavelength λ T must be resolved(Nelson 2006). This is λ T = 2 c s G Σ (4)where Σ ∼ ρH and in a marginally unstable disk, (Σand H are, respectively, surface density and disk scaleheight). Note that the most unstable wave number πλ T is exactly the inverse of the disk thickness, so that,dropping constant factors we can write λ T ∼ Q c s Ω ∼ H (Cossins et al. 2009). Since Q ∼ ρλ T in equation 4. We can then write λ T ∼ c s √ Gρ .We recall that the Toomre wavelength is borne out oflinear perturbation theory for axisymmetric local per-turbations. It is a conservative resolution marker for ourpurpose since it has been shown that the nonlinear stageof fragmentation in spiral arms occurs on a characteris-tic wavelength that is 5-6 times smaller than the Toomrewavelength (Boley 2009; Tamburello et al. 2015). In theinterface regions the spiral pattern dominates the flowhence the same considerations would apply. Recallingagain that, near cusps, we can write √ Gρ ∝ V c R c and V c R c ∼ ∇ × v , we obtain: λ T h ∼ c s h | ∇ × v | (5) We define q ≡ c s h | ∇ × v | (6)where h is the kernel smoothing length. Since theToomre wavelength should be resolved by at least onekernel(Nelson 2006), the resolution will not be adequateif q < h has to decrease, in order to maintain q large and thusfollow a self-gravitating fluid appropriately. The upperpanel of Figure 7 shows that, as we expected, q < q increases as resolu-tion is increased,as the dominant effect is now the de-crease of the smoothing length. Note that Figure 7 usesthe MFM runs. The lower panel of Figure 7 also showsthat q is always high at all resolutions when the ini-tial conditions are relaxed (TIC runs),which reflects thefact that, even for MFM, only in this case there is exactconvergence for the critical cooling time (see Table 1).Indeed in TIC runs the velocity field of the flow variesmore smoothly across the disk, with no transition phaseand thus no cusps (Figure 3). This confirms the claimof (Paardekooper et al. 2011) on the importance of theinitial setup for convergence studies. This is becauseartificial angular momentum transport during the tran-sition from a smooth disk to a fully developed spiral pat-tern is avoided. Yet table 1 also shows that only MFMshows exact convergence among TIC runs ( β crit = 3).The SPH runs still have excessive viscous dissipation ofangular momentum even in presence of more moderatevelocity gradients in shearing flows. While we did notrun ISPH with TIC conditions, we have already shownthat the viscosity switch alone becomes rather ineffec-tive at high resolution, resulting in a behaviour similarto TSPH (Figure 6). DISCUSSIONFragmentation in self-gravitating disks is known tobe a process highly sensitive to the numerical imple-mentation of the hydrodynamical equations and reso-lution. Earlier work carried out with both grid-basedand SPH codes highlighted the importance of numeri-cal approaches,including the effect of artificial viscosity,for the simple case of locally isothermal disks (Pickettet al. 2001; Mayer et al. 2004; Durisen et al. 2007; Mayer& Gawryszczak 2008). Strong dependence on numericshas also been seen in cloud fragmentation simulations.0
Deng et al.
Figure 5.
Vorticty maps. We show the modulus of vorticity in the disk mid-plane, in code units (axes in code length unitsare shown in red at the bottom left corner). In the innermost regions, vorticity is dominated by the Keplerian rotation, see thevorticity profile in figure 4. In the MFM simulations(top row), as the resolution increases the vorticity due to local overdenseregions increases, as such regions become better resolved and are characterized by strong departures from the background flowvelocity. The corresponding peak volume density in code units, in the interface region (located at 8-11 code length units), fromleft to right, upper to lower, is 0.0026, 0.0081, 0.019, 0.0013, 0.0048, 0.0031.
Figure 6.
Peak volume density(left) and peak vorticity(right) at r >
8, which includes the interface region, in the first 2.7ORPs. We show results for TSPH, ISPH ad MFM runs at different resolutions. The peak density of SPH simulations growsexponentially after the spirals reach the outer edge of the disk at about 2 ORPs. The peak vorticity is uncorrelated with thepeak density, see discussion in section 3
For example, Bate (2011) carried out SPH simulations ofmolecular cloud collapse to form protostars and disks,and found more fragmentation with increasing resolu- tion. Hu et al. (2014) ran galaxy simulations with theSPH code GADGET (Springel 2005) and showed differ- ragmentation of self-gravitating disks Figure 7.
The q factor at 1.07 ORPs in MFM runs at different resolution. The top panel shows the standard runs, while thebottom panel shows the TIC runs. As the resolution increases q increases mainly because of a smaller smooth length. However,in the interface regions, especially between spirals, q can become very small which indicates strong numerical dissipation (seedark spots, even more prominent at higher resolution in the top panel). These regions are not present in the TIC runs, wherecorrespondingly q is higher on average and increases with increasing resolution everywhere, in sharp contrast with the runs inthe top panel. The innermost regions have small q due to strong background shear but are Toomre stable and axisymmetric soour arguments in Section 3 do not apply. ent artificial viscosity switches can lead to very differentfragmentation results.In this paper we have shown that simply increasing theresolution does not yield a more correct answer for diskfragmentation, rather the nature of viscous dissipationand how it couples with the properties of the flow is thecrucial aspect. Since the properties of the flow changewith resolution, studying convergence as a function ofresolution alone is insufficient, Instead comparing dif-ferent hydro methods, and in particular different meth-ods with varying numerical dissipation, reveals the realnature of the problem. Angular momentum dissipationin strong shearing regions, which we study using vortic-ity, is the main reason behind non-convergence in SPHmethods. ISPH,which suppresses artificial viscosity inpure shearing flows, improves considerably but does notfully converge, likely because in regions of high vortic-ity the flow is also highly compressive so that viscosityturns on based on the shock tracking scheme.MFM is the only method for which we were ableto prove convergence of the critical cooling timescale among the methods explored so far, but even in thiscase exact convergence requires relaxation of the initialconditions in order to avoid transients that lead to sharpflow velocity gradients and high local vorticity at inter-face regions. The important role of the initial conditionsis supported also by a recent study of Backus & Quinn(2016) with locally isothermal disks using a modern for-mulation of the SPH hydro force, but still standard ar-tificial viscosity, using the ChaNGa code. It is also inline with previous findings by Paardekooper et al. (2011)with the FARGO code.Clearly we have only considered a small set of hydromethods in this paper.Even in the context of SPH manydifferent revised formulations of the hydro force have ap-peared in the literature over last few years which wouldbe worth exploring, such as SPHS (Hayfield et al. 2011),GDSPH (Keller et al. 2014; Tamburello et al. 2015)and the, pressure-entropy formulation(PSPH) in Hop-kins (2013). Thermal diffusion alone has been shown tobe effective at removing artificial surface tension (Price2008; Shen et al. 2010).It may play a role as it does in all2 Deng et al. problems where sharp fluid interfaces appears (Agertzet al. 2007). These modifications of other aspects of theSPH formulation could also affect β crit . The smoothedcooling approach of Rice et al. (2014) is also another ex-ample of how the fragmentation problem is sensitive tothe details of the numerical implementation within thesame category of methods. In order to at least partiallyaddress this we rerun the HR simulation from relaxedinitial conditions (TIC) using PSPH and thermal diffu-sion in the GIZMO code (Hopkins 2015), and with theCullen & Dehnen viscosity formulation. We find thatthe disk fragments for β = 4 −
5. Therefore using PSPHdoes not lead to an improvement relative to SPH, con-firming the central role of residual viscous dissipation asopposed to other aspects of the SPH scheme, for exam-ple the SPH “E0” zeroth-order errors(Read et al. 2010).Dehnen & Aly (2012) showed that the kernel functionsdoes affect numerical dissipation in shear flow simula-tions.The Wendland C4 kernel shows the best conserva-tion properties in their numerical tests. Although MFMsolves the fluid equations on the volume elements con-structed from the particle distribution, the volume ele-ments themselves are constructed using a specified ker-nel function hence it is relevant to ask what effect thechoice of the kernel function has on the flow solution.Therefore we rerun the MFM-HR-Beta3-TIC simulationwith Wendland C4, using 200 neighbours. We foundthat the simulation still fragments at β crit = 3. This isreassuring as it strengthens the conclusion that the re-sults have truly converged,namely we cannot push β crit even further by reducing numerical dissipation using amore conservative kernel function. However, a word ofcaution must be said. With such many neighbours theeffective gravitational force resolution is half the forceresolution in the standard HR simulation, being equalto that in the LR simulations, and the force resolutionat later stage is also not directly comparable due to itsfull adaptive nature.Meru & Bate (2012) could not prove exact conver-gence for a set of FARGO simulations, for which theyalso considered varying the coefficient of viscosity. In-deed, despite showing a trend with resolution that waspointing towards convergence, they did not find β crit to maintain the exact same value when they increasedtheir resolution by a factor of 2. The asymptotic valueof β crit suggested by their results was also higher thanthe corresponding value found in their SPH runs. Inlight of our findings that show how convergence occursto a very small vaue of β crit = 3,we argue that thepublished FARGO simulations might suffer from evenstronger numerical dissipation than SPH, resulting notonly from numerical viscosity but also from advection errors occurring in cuspy regions that are not at allaxisymmetric and hence cannot be captured properlyby the cylindrical geometry of the grid. For the samereason, significant improvements can be expected usingrelaxed initial conditions as deviations from axisymme-try are less severe due to the absence of the transitionphase. This is once again in line with the findings of(Paardekooper et al. 2011) Advection errors in the in-terface regions could be even more severe in Cartesiangrids,unless aggressive AMR is used to better capturethe curvature of the flow. We plan to reassess this ina future comparison employing more than one type ofgrid-based code.Finally,at fixed smoothing length the density increasein self-gravitating flow is strongly affected by gravita-tional softening. Adaptive softening is expected to playa key role as it increases the force resolution as par-ticles converge towards a cusp or sharp interface. Werun additional TSPH simulations with fixed softeningto quantify the importance of the choice of softening forthe problem under study.Figure 8 indicates that also fixed softening simulationsalso experience an exponential peak density growth. Inthe TSPH-LR run, when we choose a gravitational soft-ening of 0.1 code units, which equals half the smoothinglength in the disk middle plane at r = 9 ( see the curveof the TSPH-LR-Beta4-S0.1 run), the peak volume den-sity evolution is similar to that in the the adaptive soft-ening run. If we decrease the softening to 0.05(TSPH-LR-Beta-S0.05) , the peak volume density increase sig-nificantly. The TSPH-HR-Beta4-S0.05 has even higherpeak density than TSPH-LR-Beta4-S0.05, and also evenhigher than the peak density in the adaptive softeningrun, resulting from the denser cusps forming at higherresolution. As a result, runaway fragmentation in pres-ence of vorticity damping can be even more problematicthan in the reference adaptive softening run.The results found in this paper suggest that disksneed to cool really fast to fragment, on a timescalenearly identical to the local orbital time. This rein-forces the notion that it is in the outer disk regions,at R >
30 AU (Rafikov 2009; Clarke & Lodato 2009),that the conditions are favourable for fragmentation.However, our study focused on idealized isolated diskmodels,neglecting the effect of both external triggers,such as mass accretion onto the disk from the envelopein the early stages of disk evolution (Boley & Durisen2010; Vorobyov & Basu 2010) or internal triggers suchas perturbations by a previously-formed fragment (Meru2015). Mass loading has already been shown to bringdisks to fragmentation even when radiative cooling isslow (Boley 2009), or even formally absent as in the col- ragmentation of self-gravitating disks Figure 8.
The effect of fixed gravitational softening andresolution on peak density inside the interface region at r > lapse of polytropic turbulent cloud cores (Hayfield et al.2011). . The correct treatment of radiative cooling andheating originating from accretion ultimately requiresradiative transfer to be employed (eg (Meru & Bate2011b; Rogers & Wadsley 2011; Mayer et al. 2016) How-ever, attention should be paid to the initial conditionsetup and accretion boundary condition. Complex flowstructures and overdensities arising in turbulent collapsemay lead to strong numerical dissipation and cause spu-rious fragmentation for the same reasons described inthis paper. Therefore a weakly dissipative method suchas MFM holds promise to be most suitable for radiativedisk formation studies as well.The complexity inherent to model disk fragmentationemerging from our study suggests that it is too prema-ture to use existing simulations to compare with exo-planet detection via imaging surveys for wide orbit gasgiants (eg (Vigan et al. 2017)). Indeed, once the forma-tion stage of clumps is robustly modeled, their evolutionvia further collapse (Galvagni et al. 2012; Szulgyi et al.2017) and migration in the disk (Baruteau et al. 2011;Malik et al. 2015) is also not completely understood, andnot yet firm from the numerical side. Fast inward mi-gration of clumps is routinely found which alone wouldinvalidate any direct comparison with wide orbit gas gi-ants as a way to infer the role of disk instability, but themigration rate in self-gravitating disks might be affecteddirectly or indirectly by numerical dissipation effects. Acomparison of hydro methods is warranted in this caseas well.Finally, there is an interesting upshot of our study.This is that actual physical sources of viscosity in disks,such as the magnetorotational instability (MRI) (Balbus & Hawley 1991)might promote fragmentation if they candissipate angular momentum efficiently in overdense re-gions with strong shear. This will be explored in futurework. CONCLUSIONSWe performed 3D hydrodynamic simulations usingSPH with different artificial viscosity implementationsas well as with a new Lagrangian method, MFM, whichemploys a Riemann solver on a volume partition gener-ated by particles, thus not requiring any explicit numer-ical viscosity.Our TSPH simulations agree well with a previousstudy by Meru & Bate (2011a, 2012).They confirm the non-convergence of the critical cooling rate for disk frag-mentation in a self-gravitating disk. MFM instead at-tains convergence, and to a significantly lower value ofthe critical cooling time, β crit = 3. ISPH simulations,which adopt the more conservative artificial viscosity im-plementation by Cullen & Dehnen (Cullen & Dehnen2010), resulting formally in no shear viscosity, exhibitan intermediate behaviour. Finally,exact convergence,even in MFN, occurs only with relaxed initial conditions.We find a coherent explanation for our results,showingthat the driving effect is numerical dissipation of angularmomentum in strong shearing flows. This is exacerbatedas resolution is increased, despite shear viscosity in SPHis formally decreasing with resolution,because the flowdevelops stronger cusps with large velocity derivativeswhere vorticity is artificially damped. In those regionsthe characteristic wavelength associated with fragmen-tation, conservatively measured with the Toomre wave-length, remains poorly resolved even when the numberof particles is increased by nearly three orders of mag-nitude, hence the high sensitivity on the numerical dis-sipation. Relaxed initial conditions exhibit a more con-vergent behaviour because they avoid a transient phasein which overdense cusps and associated high velocityderivatives occur, thus suffering less the impact of arti-ficial angular momentum dissipation by numerical vis-cosity.We stress that the small value β crit = 3 in the con-verged MFM simulations is significantly smaller thanany value previously found in the literature when tryingto assess convergence with increasing resolution. Thisindicates that the trends suggestive of asymptotic con-vergence reported in the literature for both SPH codesand FARGO (Meru & Bate 2012; Rice et al. 2014) werejust reflecting saturation of numerical angular momen-tum dissipation rather than a decreasing numerical dis-sipation with increasing resolution.4 Deng et al.
We have further assessed the validity of our statementsby running additional simulations that explore other as-pects of SPH implementations,such as adaptive soften-ing, kernel function and form of the hydro force. Thesetests confirm our claim that only MFM attains exactconvergence.In summary, our results indicate that the use of hydromethods with minimal numerical viscosity is the firstnecessary step towards predictive simulations of diskinstability because the issues we highlighted would still play a role in complex simulations with radiative trans-fer and envelope accretion.We thank Jim Stone, Frederic Masset for useful dis-cussions. We are grateful to Philip Hopkins for helpingrun the code. Pynbody(Pontzen et al. 2013) is usedfor the visualization. We acknowledge support from theSwiss National Science Foundation via the NCCR Plan-etS. F.M. acknowledges support from The LeverhulmeTrust and the Isaac Newton Trust.REFERENCES
Agertz, O., Moore, B., Stadel, J., et al. 2007, MonthlyNotices of the Royal Astronomical Society, 380, 963Backus, I., & Quinn, T. 2016, MNRAS, 463, 2480Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214Balsara, D. S. 1995, Journal of Computational Physics, 121,357Baruteau, C., Meru, F., & Paardekooper, S.-J. 2011, inSF2A-2011: Proceedings of the Annual meeting of theFrench Society of Astronomy and Astrophysics, ed.G. Alecian, K. Belkacem, R. Samadi, & D. Valls-Gabaud,459–461Bate, M. R. 2011, MNRAS, 417, 2036Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS,277, 362Boley, A. C. 2009, ApJL, 695, L53Boley, A. C., & Durisen, R. H. 2010, ApJ, 724, 618Boss, A. P., & Bodenheimer, P. 1979, ApJ, 234, 289Clarke, C. J., & Lodato, G. 2009, Monthly Notices of theRoyal Astronomical Society: Letters, 398, L6Cossins, P., Lodato, G., & Clarke, C. J. 2009, MNRAS,393, 1157Cullen, L., & Dehnen, W. 2010, MNRAS, 408, 669Dehnen, W., & Aly, H. 2012, MNRAS, 425, 1068Durisen, R. H., Boss, A. P., Mayer, L., et al. 2007,Protostars and Planets V, 607Galvagni, M., Hayfield, T., Boley, A., et al. 2012, MNRAS,427, 1725Gammie, C. F. 2001, ApJ, 553, 174Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ,440, 742Hayfield, T., Mayer, L., Wadsley, J., & Boley, A. C. 2011,MNRAS, 417, 1839Helled, R., Bodenheimer, P., Podolak, M., et al. 2014,Protostars and Planets VI, 643Hernquist, L., & Katz, N. 1989, ApJS, 70, 419Hopkins, P. F. 2013, Mon. Not. R. Astron. Soc., 428, 2840—. 2015, MNRAS, 450, 53 Hu, C. Y., Naab, T., Walch, S., Moster, B. P., & Oser, L.2014, MNRAS, 443, 1173Kaufmann, T., Mayer, L., Wadsley, J., Stadel, J., & Moore,B. 2007, MNRAS, 375, 53Keller, B. W., Wadsley, J., Benincasa, S. M., & Couchman,H. M. P. 2014, MNRAS, 442, 3013Klee, J., Illenseer, T. F., Jung, M., & Duschl, W. J. 2017,ArXiv e-prints, arXiv:1704.02193Malik, M., Meru, F., Mayer, L., & Meyer, M. 2015, ApJ,802, 56Mayer, L., & Gawryszczak, A. J. 2008, in AstronomicalSociety of the Pacific Conference Series, Vol. 398,Extreme Solar Systems, ed. D. Fischer, F. A. Rasio, S. E.Thorsett, & A. Wolszczan, 243Mayer, L., Peters, T., Pineda, J. E., Wadsley, J., & Rogers,P. 2016, ApJL, 823, L36Mayer, L., Quinn, T. R., Wadsley, J., & Stadel, J. 2004, inBulletin of the American Astronomical Society, Vol. 36,AAS/Division of Dynamical Astronomy Meeting ragmentation of self-gravitating disks Paardekooper, S.-J., Baruteau, C., & Meru, F. 2011,MNRAS, 416, L65Pickett, B. K., Mejia, A. C., & Durisen, R. H. 2001, inBulletin of the American Astronomical Society, Vol. 33,American Astronomical Society Meeting Abstracts, 1397Pontzen, A., Roˇskar, R., Stinson, G. S., et al. 2013,pynbody: Astrophysics Simulation Analysis for Python, ,, astrophysics Source Code Library, ascl:1305.002Price, D. J. 2008, Journal of Computational Physics, 227,10040Price, D. J., & Monaghan, J. J. 2007, MNRAS, 374, 1347Rafikov, R. R. 2009, ApJ, 704, 281Read, J. I., Hayfield, T., & Agertz, O. 2010, MNRAS, 405,1513Rice, W. K. M., Lodato, G., & Armitage, P. J. 2005,MNRAS, 364, 56Rice, W. K. M., Paardekooper, S.-J., Forgan, D. H., &Armitage, P. J. 2014, MNRAS, 438, 1593 Rogers, P. D., & Wadsley, J. 2011, MNRAS, 414, 913Rosswog, S., Davies, M. B., Thielemann, F.-K., & Piran, T.2000, A&A, 360, 171Shen, S., Wadsley, J., Hayfield, T., & Ellens, N. 2010,Monthly Notices of the Royal Astronomical Society, 401,727Springel, V. 2005, MNRAS, 364, 1105Szulgyi, J., Mayer, L., & Quinn, T. 2017, Monthly Noticesof the Royal Astronomical Society, 464, 3158Tamburello, V., Mayer, L., Shen, S., & Wadsley, J. 2015,MNRAS, 453, 2490Toomre, A. 1964, ApJ, 139, 1217Vigan, A., Bonavita, M., Biller, B., et al. 2017, ArXive-prints, arXiv:1703.05322Vorobyov, E. I., & Basu, S. 2010, ApJ, 719, 1896Young, M. D., & Clarke, C. J. 2015, MNRAS, 451, 3987 Deng et al.
APPENDIX A. ISOTHERMAL SPHERE COLLAPSETo further assess the low numerical viscosity in MFM we run the standard isothermal test case for the collapse of arotating molecular cloud. This test, first described by Boss & Bodenheimer (1979), is a typical test case for numericalcodes studying fragmentation(Bate et al. 1995). An initially isothermal, spherical, self-gravitating hydrogen molecularcloud starts to collapse from rest. Physical viscosity and the magnetic field are not included, so the angular momentumof the system should be conserved. The gas is assumed to be ideal gas and initial temperature is 10 K . The mass ofthe cloud is 1 M sun , and radius of the cloud is 3 . × cm . The angular velocity is Ω = 1 . × − rads − . The ratiosof thermal energy to the magnitude of gravitational energy, and rotational energy to the magnitude of gravitationalenergy of α = 0 .
25 and β = 0 .
20, respectively. The mean density is ρ = 1 . × − gcm − . No density perturbationis imposed initially. The free-fall time for the initial cloud is t ff = 5 . × s .We run this setup with 10K, 100K, 1M particles, using both MFM and TSPH. All the simulations share the sameproperties discussed below. In figure 9, we plot the specific angular momentum profile of the sphere at different times. Figure 9.
The evolution of the specific angular momentum profile in the MFM collapsing rotating isothermal sphere run with1M particles
Initially the sphere is uniformly rotating, and the special angular momentum j ∝ r , where r is the radius. Theinitial cloud radius is about 0.6 code length units. When material starts to flow inwards, angular momentum must betransported outward in order to conserve the total angular momentum of the system. Particles at the outer part ofthe cloud drift further. We show the angular momentum profile out to only one code length unit because the outerpart is noisy due to low number of particles.The angular momentum profile evolves less at higher resolution in 11, which is a sign of less numerical dissipation.The is because of the nice continuous flow property and numerical viscosity scales with h . From figure 10, we seevery similar evolution of angular momentum profile by MFM and TSPH. However, at fixed time, the TSPH angularmomentum profile(dashed lines) has evolved a bit faster, ie, more angular momentum has been transported outwards,which is reflected by the fact that the MFM curve is above the TSPH curve inside 0.6 code length units and belowoutside.We caution that the changing nature of the flow with increasing resolution,which we have shown to be a centralaspect of unstable disk evolution, does not play a role here. This also means there are no sites of strong shear and ragmentation of self-gravitating disks Figure 10.
Specific angular momentum profile with both MFM and TSPH, using 100K particles. At fixed time, the angularmomentum profile evolves slightly faster in TSPH due to more significant numerical transport.
Figure 11.
Resolution dependence of the angular momentum profile in TSPH runs, using 10K (LR), 100K (HR) and 1 Mparticles (UHR) Deng et al. B. ACCRETING PLANET TESTWe are interested in assessing how numerical viscosity affects the angular moment transport, especially in cusps atinterface regions where high velocity derivatives occur. . To avoid the complexity of resolving a full disk but carry outa test that can address a configuration with velocity derivatives, we use a 2D shearing sheet (Hawley et al. 1995) withan active accreting planet in the center to mimic the collapse under self-gravity. The simulations are similar to those inOrmel et al. (2015a) while the gravity of the planet is added following Ormel et al. (2015b). We choose the local soundspeed and the disk scale height as the natural units of speed and length, respectively. We set the gravitational constantto 1. The planets dimensionless mass m ≡ GM p/ ( c s / Ω) equals its Bondi radius in dimensionless units R Bondi /H . Inthese units the dimensionless Hill sphere r Hill = R Hill /H = ( m/ / .The gravity of the planet is added F ( t ) = F { − exp [ − ( tt inj ) ] } . The force increases gradually to mimic thecollapse due to self-gravity. t inj is the injection time. The gravity from the planet is smoothed to avoid numericalinstability and increase the efficiency of calculation. F = − mr exp [ − A ( r in r ) p − B ( rr out ) q ]. We set A=10, p=8 and B=1,q=4. The force transits from a pure Newtonian force F = − m/r for r in < r < r out to zero continuously and quickly.We set r in = 0 . , r out = m to prevent boundary condition corruption by the gravity of the planet.We run simulations with m = 0 . , t inj = 2 at a resolution of 64 ×
64 particles (this means 64 particles per scale height,which is even higher resolution than the 16 million particles UHR 3D disk runs used in this paper). The artificialviscosity α coefficient is shown in figure12 to compare directly TSPH and ISPH. ISPH, as expected, has lower artificialviscosity. We use the time evolution of peak density to measure how much mass in accreted over time inside theplanet radius r in as shown in figure 13. The faster the growth of peak density the faster must be the outward transportof angular momentum resulting from any form of numerical dissipation, including, but not only, explicit artificialviscosity. Figure 13 highlights the lower dissipation of MFM. Indeed MFM always has the lowest peak density whileISPH and TSPH are almost identical despite the lower value of the α coefficient. We argue the faster accretion inSPH methods is due to higher numerical viscous dissipation. The Cullen & Dehnen switch is most effective away fromshocks and it is not able to effectively suppress numerical viscosity in a flow with a large velocity gradient as in thistest. Figure 12.