The response of self-graviting protostellar discs to slow reduction in cooling timescale: the fragmentation boundary revisited
aa r X i v : . [ a s t r o - ph ] A ug Mon. Not. R. Astron. Soc. , 000–000 (0000) Printed 31 October 2018 (MN L A TEX style file v2.2)
The response of self-graviting protostellar discs to slow reduction incooling timescale: the fragmentation boundary revisited
C.J. Clarke , E. Harper-Clark , G. Lodato Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA Department of Astronomy and Astrophysics, 50 St. George Street, Toronto, Ontario, Canada, M5S 3H4 Department of Physics and Astronomy, University of Leicester, University Road, Leicester, LE1 7RH
31 October 2018
ABSTRACT
A number of previous studies of the fragmentation of self-gravitating protostellar discshave involved suites of simulations in which radiative cooling is modeled in terms of a cool-ing timescale ( t cool ) which is parameterised as a simple multiple ( β cool ) of the local dynamicaltimescale. Such studies have delineated the ‘fragmentation boundary’ in terms of a criticalvalue of β cool ( β crit ) such that the disc fragments if β cool < β crit . Such an approach howeverbegs the question of how in reality a disc could ever be assembled in a state with β cool < β crit .Here we adopt the more realistic approach of e ff ecting a gradual reduction in β cool , as mightcorrespond to changes in thermal regime due to secular changes in the disc density profile. Wefind that the e ff ect of gradually reducing β cool (on a timescale longer than t cool ) is to stabilisethe disc against fragmentation, compared with models in which β cool is reduced rapidly (overless than t cool ). We therefore conclude that the ability of a disc to remain in a self-regulated,self-gravitating state (without fragmentation) is partly dependent on the disc’s thermal history,as well as its current cooling rate. Nevertheless, the e ff ect of a slow reduction in t cool appearsonly to lower the fragmentation boundary by about a factor two in t cool and thus only permitsmaximum ‘ α ’ values (which parameterise the e ffi ciency of angular momentum transfer in thedisc) that are about a factor two higher than determined hitherto. Our results therefore do notundermine the notion that there is a fundamental upper limit to the heating rate that can be de-livered by gravitational instabilities before the disc is subject to fragmentation. An importantimplication of this work, therefore, is that self-gravitating discs can enter into the regime offragmentation via secular evolution and it is not necessary to invoke rapid (impulsive) eventsto trigger fragmentation. Key words: accretion, accretion discs – star: formation – gravitation – instabilities – stars:formation
Following the seminal work of Gammie (2001), there has been con-siderable progress in recent years in understanding the behaviourof self-gravitating accretion dics (see Durisen et al. 2007 and refer-ences therein). A number of simulations (Gammie 2001; Rice et al.2003; Lodato & Rice 2004, 2005) have demonstrated that if thethermodynamic properties of the disc are evolved according to athermal equation (involving a cooling term parameterised in termsof a cooling timescale, t cool ), then the disc may be able to establisha self-gravitating, self-regulated state. In this state, the Toomre Q parameter: Q = c s κπ G Σ , (1)(where c s is the sound speed, κ is the epicyclic frequency (equal tothe angular velocity Ω in a Keplerian disc) and Σ is the disc surface density) hovers at a value somewhat greater than unity over an ex-tended region of the disc. Whereas the state Q =
1, corresponds toa situation of marginal stability against axisymmetric perturbations,in the self-regulated state the disc is instead subject to a variety ofnon-axisymmetric self-gravitating modes whose e ff ect, through theaction of weak shocks, is to dissipate mechanical energy (i.e. ki-netic and potential energy of the accretion flow) as heat. Thermalequilibrium is then attained through the balancing of such heatingby the prescribed radiative cooling: in essence, self-regulation re-sults when the amplitude of these modes is able to self-adjust so asto maintain thermodynamic equilibrium against the relevant energyloss processes.The above studies have all found, however, that such self-regulation is only possible in the case that the cooling timescaleis not too short: stability demands that β cool = t cool / Ω − exceeds a c (cid:13) C. Clarke, E. Harper-Clark & G. Lodato critical value which, for discs with adiabatic index of 5 /
3, is ∼ β cool . Discs with β cool < β crit then fragment on the localcooling timescale (i.e. a few times the local dynamical timescale),thus begging the question of how such unstable initial conditionscould ever have been set up in the first place.A more likely scenario for disc fragmentation is that the discis instead set up in self-regulated, self-gravitating state and thenconditions gradually change so that β cool is lowered. (For example,continued infall of material onto a disc or secular re-arrangement ofmaterial in the disc due to the action of gravitational torques couldalter the surface density profile of the disc and allow it to enter anew cooling regime with lower β cool ). It is not however clear thatthe fragmentation boundary would be the same in the case that β cool is gradually reduced.In this paper we conduct a suite of idealised simulations inwhich we explore whether the fragmentation boundary just de-pends on the instantaneous value of β cool (as has been assumedhitherto) of whether the system ‘remembers’ the history of howit evolved to a point of given β cool . Such a (‘toy model’) approach,is complementary to studies (Boley et al. 2006; Mayer et al. 2007;Stamatellos et al. 2007, see also the analytical estimates by Rafikov2005; Rafikov 2007) which attempt to achieve ever-increasingverisimilitude via the incorporation of more realistic treatments ofradiative transfer. Here, instead, we make no claims that the sim-plified cooling law (for example, the assumption that β cool is spa-tially uniform) actually corresponds to a situation encountered ina real disc, because our aim is to isolate a particular physical ef-fect (i.e. the timescale on which the fragmentation boundary isapproached). The computational expense of ‘realistic’ simulationshowever prevents their use to study secular e ff ects: even in the caseof the present ‘toy’ simulations, it is impracticable to run simula-tions over the long timescales on which the Σ profile changes due togravitational torques or infall. We can nevertheless assess the e ff ectof relatively slow changes in β cool on the fragmentation boundarythrough imposing an ad hoc reduction in the value of β cool and canapply this insight to the secular evolution of real discs.In particular, we want to examine the cause of the fragmenta-tion for β cool < β crit , that has been found in previous simulations. Isthis (i) due to the disc’s inability to maintain - under any circum-stances - a gravitational heating rate that can match the imposedhigh cooling rate? This is the hypothesis of Lodato & Rice (2005),who identify the minimum value of β cool with a maximum value ofthe gravitationally induced angular momentum transfer that can bedelivered by a disc without its fragmenting. They parameterise thisstate of maximal angular momentum transfer in terms of the ratioof the r , φ component of the stress tensor to the thermal pressure,i.e., by analogy with the equivalent expression for a viscous disc,in terms of a maximum in the well known viscous ‘ α ’ parameter(Shakura & Sunyaev 1973). A critical value of β cool of ∼ α of ∼ . on the shorttimescale ( t cool ) on which the disc is cooling? If this were the case,then with su ffi ciently gradual approach to the regime of low β cool ,the disc could in principle deliver a value of α that exceeded theabove limit by a generous margin. We can obviously distinguish between these alternatives byinvestigating the case in which β cool is reduced on a timescale τ that is longer than t cool , since in this case the disc temperature willfall via a sequence of thermal equilibrium states (on timescale τ ),rather than dropping on timescale t cool . The aim of this investigationis thus to see whether the disc is more resistant to fragmentation inthe regime that τ > t cool . If it is not , then the manner in which thedisc approaches the fragmentation boundary is unimportant. If, onthe other hand, it is found that rapid changes in cooling regime arerequired, then it may be necessary to invoke impulsive events (suchas an external dynamical interaction) to trigger fragmentation.In Section 2 we describe the numerical setup, discuss our re-sults in Section 3 and in Section 4 we present some conclusions. Our three-dimensional numerical simulations are carried outusing SPH, a Lagrangian hydrodynamic scheme (Benz 1990;Monaghan 1992). The general implementation is very similarto Lodato & Rice (2004), Lodato & Rice (2005) and Rice et al.(2005). The gas disc is modeled with 250,000 SPH particles(500,000 in a run used as a convergence test) and the local fluidproperties are computed by suitably averaging over the neighbour-ing particles. The disc is set in almost Keplerian rotation (allowingfrom slight departures from it to account for the e ff ect of pressureforces and of the disc gravitational force) around a central pointmass onto which gas particles can accrete if they get closer thanthe accretion radius, taken to be equal to 0.5 code units.The gas disc can heat up due to p d V work and artificial vis-cosity. The ratio of specific heats is γ = /
3. Cooling is here imple-mented in a simplified way, i.e. by parameterizing the cooling ratein terms of a cooling timescale: d u i d t ! cool = − u i t cool , (2)where u i is the internal energy of a particle and the coolingtimescale t cool is assumed to be proportional to the dynamicaltimescale, t cool = β cool Ω − , where β cool is varied according to a time-dependent prescription (see Section 2.3 below).Artificial viscosity is introduced using the standard SPH for-malism. The actual implementiation is very similar to the one usedin Rice et al. (2005), that is we set the two relevant numerical pa-rameters to α SPH = . β SPH = . The main physical properties of the disc at the beginning of the sim-ulation are again similar to those of Lodato & Rice (2004, 2005).The disc surface density Σ is initially proportional to R − (where R is the cylindrical radius), while the temperature is initially propor-tional to R − / . Given our simplified form of the cooling function,the computations described here are essentially scale free and canbe rescaled to di ff erent disc sizes and masses. For reference, we willassume that the unit mass (which is the mass of the central star) is1 M ⊙ and that the unit radius is 1 AU . In this units the disc extendsfrom R in = . AU to R out = AU . The normalization of the sur-face density is generally chosen such as to have a total disc mass of c (cid:13) , 000–000 he fragmentation boundary revisited Simulation x β hold N fragmentationF1 10.5 — 250K yesF2 10.5 3 250K yesV 105-10.5 3 250K yesS1 105 3 250K noSh 105 3 500K yesS2 105 2.75 250K noS3 105 2.62 250K yesVS1 314 3 250K noVS2 314 2.75 250K yes Table 1.
Details of the various simulations discussed in this paper. The dif-ferent columns indicate: the name of the run, the value of the parameter x determining the speed of the reduction of the cooling time, the value of β hold (if any) at which the cooling time was held fixed after reduction, thenumber of particles used in the run N and whether fragmentation did oc-cur or not. Simulation V was performed with an initially slow reduction of β (with x = x = . β = M disc = . M ⊙ , while the temperature normalization is chosen soas to have a minimum value of Q =
2, which is attained at the outeredge of the disc.Initially, the disc is evolved with constant β cool = .
5, thisvalue of β cool being in the regime where previous work Gammie(2001), Rice et al. (2005) has shown that the disc does not frag-ment. The general features of this initial evolution is described indetail in Lodato & Rice (2004). The disc starts cooling down untilthe vertical scale-length H is reduced such that H / R ≈ M disc / M ⋆ = .
1. At this point the disc becomes Toomre unstable and developsa spiral structure that heats up the disc and maintains it close tomarginal stability. We have evolved the disc with this value of β cool for 7 . Q = R = −
23 A.U. ). β cool After evolution of the disc with β cool = β cool (0) = . ff ect a linear reduction of β cool on a timescale T , i.e. β cool ( t ) = β cool (0) (cid:18) − tT (cid:19) (3)where we set T = x Ω − ( R out ).Such a prescription implies that the timescale τ ( t ) on whichthe local instantaneous value of t cool (i.e. t cool ( R , t )) drops to zero is τ ( t ) = β cool ( t ) | ˙ β cool | = x β cool (0) RR out ! − . t cool ( R , t ) . (4)We adopt three values of x : x = . x =
105 (slow) and x =
314 (very slow). In the fast case, τ ( R , t ) ∼ t cool ( R , t ) in the outerdisc so t cool is changing faster than the disc can come into thermalequilibrium at that value of t cool . In the slow case, τ ( r , t ) > t cool ( r , t )so that the disc is everywhere able to come into thermal equilibriumat that t cool . This situation is even more amply satisfied in the veryslow case.For each value of x , we run the simulation until a fragmentforms (at β cool = β frag ). In those cases where the timing of frag-mentation suggests that the slow reduction of β cool is acting so as tostabilise the disc at lower β cool , we test this hypothesis by turning o ff the reduction in β cool when it attains a value equal to β hold ( > β frag ).We then experiment with values of β hold in order to find the min-imum value of β hold at which the disc does not fragment over the Figure 1.
Time evolution of the parameter β cool in the various models. Thethree solid lines refer to cases where β cool is first decreased and then heldat β cool = β cool (V). Anasterisk at the end of the line indicates fragmentation at this time, whereasruns without an asterisk are unfragmented at the end of the simulation. duration of the numerical experiment. Table 1 and Fig. 1 summa-rize the main details of the various runs we have performed, wherethe ‘F’ simulations are the ‘fast’ ones, the ‘S’ are the ‘slow’ onesand the ‘VS’ are the ‘very slow’ ones (see discussion in Section 3below). The simulation named V was performed with an initiallyslow reduction of β (with x = x = . β = ff ect our results (see below). One important aspect that needs to be taken into account is whetherthe resolution of our simulation is high enough to reproduce frag-mentation, when it occurs. Resolution criteria for fragmentationwith SPH codes have been discussed by Bate & Burkert (1997).They obtained that SPH correctly reproduces fragmentation if therelevant Jeans mass contains at least 100 SPH particles, that istwice the typical number of neighbours ( N neigh =
50) within onesmoothing region. More recently, Nelson (2006) has revisited thisissue focussing on fragmentation in self-gravitating discs and hasfound a slightly more stringent criterion, requiring that the Jeansmass is resolved with three times as many particles as required byBate & Burkert (1997). In a gravitationally unstable disc, the mostunstable wavelength is given by λ = c / G Σ . The Jeans mass (or,as Nelson 2006 calls it, the “Toomre mass”) is then given by: M J = π Σ λ = π c s G Σ = π Q (cid:18) HR (cid:19) Σ R . (5) c (cid:13) , 000–000 C. Clarke, E. Harper-Clark & G. Lodato
Figure 2.
Images of Sh at fragmentation (left) and S1 (centre) at the same time. Although Sh has fragmented, only one fragment is seen at large radius andaside from the immediate area around the fragment the discs are very similar. This should be contrasted with the profusion of fragments in F2 at the point offragmentation
The cumulative disc mass at radius R is given by: M disc ( R ) = π Σ R = M disc RR out , (6)since in our setup Σ ∝ R − (see above). We can then rewrite theJeans mass using eq. (6), as: M J = π Q (cid:18) HR (cid:19) RR out ! m p N tot , (7)where we have also used M disc = m p N tot , where N tot is the to-tal number of particles used and m p is the mass of an individualSPH particle. In order to properly resolve fragmentation, we re-quire that M T > m p N reso , where N reso = N neigh = N reso = N neigh =
300 according toNelson’s more restrictive criterion, and recalling that in our simu-lations the mean number of neighbours per particle is 50. We thenobtain that we have enough resolution at radii R that satisfy: RR out & π Q q N reso N tot ! / , (8)where q = M disc / M ⋆ and where we have aso used the relationshipbetween disc thickness and the parameter Q : HR = Q M disc ( R ) M ⋆ , (9)that can be easilty derived from Eq. (1). Based on equation (8) wecan then conclude that for q = .
1, as used in the present paper,fragmentation is well resolved at radii R &
5. Note that, since N reso only enters eq. (8) to the power of one third, if we had used the morerestrictive condition of Nelson (2006), we would only increase ourminimum radius by a factor 1.4. As shown in Fig. 2, whenever weobserved fragmentation, this occurred outside R ≈
5, so that we canbe confident that we do resolve the relevant mass and length scalesfor fragmentation.A second aspect related to resolution is that we require arti-ficial viscosity to play a role only when modeling shocks. In or-der to ensure this, we then require that the velocity di ff erence ac-cross a smoothing kernel is subsonic, i.e. h Ω < c s , where h is thesmoothing length. This in turn requires that the smoothing length issmaller than the disc thickness H = c s / Ω . We have indeed checkedthat, even at the lower resolution of 250,000 particles, the averagesmoothing length is a fraction ≈ . We find that in the fast case ( x = . β cool = .
75, i.e. about 8 . β cool commenced. Since fragmentation alwaystakes about a dynamical timescale to get under way, it follows that,as expected, the ‘fast’ case behaves like the usual case where a fixed β cool is imposed.We however see di ff erent behaviour in the slow case: here wefind that when β cool is reduced to, and then held at, β hold = . β hold =
3, even after integra-tion for 63 outer dynamical timescales at this β cool value. On theother hand, when β cool was instead held at 2 .
62, it fragmented af-ter a further ∼
18 outer dynamical timescales, so it would appearthat the fragmentation boundary is at around 2 .
7. This is in strongcontrast with the value of ∼ β cool is imposed. We hesitate to say that we have proved thata disc will never fragment when brought to such a low value of β cool value at this slow rate, since our experience shows that whereone is close to the limit of marginally stable β cool , fragmentationmay ensue after long timescales, and that its timing may dependon numerical noise that can be a ff ected by resolution. Indeed, wefound that when we re-ran the x =
105 simulation at higher res-olution ( N = , β hold =
3, it eventually didform a fragment at large radius . We however show the disc struc-ture in this simulation at the point of fragmentation and contrast itwith the corresponding situation when the cooling time is rapidlyreduced and then held at constant β cool = ff erent inthe two cases, with the ‘rapid’ simulation containing a number of We have tested whether this resistance to fragmentation in the slow caseis simply because the disc takes longer to reach β cool = β cool = x = .
5) reduction in β cool between β cool = β cool = x which controls fragmentation. c (cid:13) , 000–000 he fragmentation boundary revisited regions that are on the point of fragmentation at the moment thatthe first fragment appears. Our interest here is not in defining pre-cise boundaries at which fragmentation will or will not occur (sincethe definition of such a boundary is always contingent on the dura-tion of the simulation) but in demonstrating that the structure of thedisc is indeed a ff ected not just by the instantaneous value of β cool (and hence on the heating rate that has to be delivered through theaction of the self-gravitating modes) but also on the history of howthe disc arrived at such a value of β cool .This then raises the possibility (which we discussed in Section1) that the lower limit on β cool for self-regulation might representthe di ffi culties that a disc might have in achieving a self-regulatedstate on an appropriately short timescale, rather than a fundamentalupper limit on the dissipation rate that can be provided by gravita-tional modes in the absence of fragmentation. In principle, then, wecould envisage a situation where the disc might be self-regulated atan arbitrarily low value of β cool (i.e. where the non-linear develop-ment of the spiral modes delivered an arbitrarily high heating ratewithout the disc fragmenting) provided that the disc approachedthis state su ffi ciently slowly. Such a conclusion would contradictthat of Lodato and Rice 2005, who interpreted the fragmentationboundary in terms of a (history independent) limit on the maximum α value delivered by such instabilities.In order to explore this further, we ran the very slow ( x = β cool than for the x = ff erence in the results for the x =
105 and x =
314 case, the lowest values at which β cool could beheld being respectively 2 .
75 and 3 .
00 for the two cases (for N = ,
000 in both cases).
We have found that the rate at which the cooling timescale ischanged indeed a ff ects the minimum value of β cool at which thedisc can exist in a stable, self-regulated state. As expected, this ef-fect is only manifest when the the cooling timescale is varied ona timescale ( τ ) that is longer than the cooling timescale, since for τ < t cool , the temperature always falls on a timescale t cool , irre-spective of τ . We find that when τ > t cool , the self-regulated stateis sustainable at cooling times that are about a factor two less thanthose that are possible when a fixed cooling timescale is imposed atthe outset of the simulation. This implies that (in the slow coolingcase) the gravitational instabilities are able to deliver about twicethe heating rate without the disc fragmenting. In terms of the ‘vis-cous alpha’ description of such instabilities (Shakura and Sunyaev1973, Gammie 2001, Lodato & Rice 2005), the maximum α de-liverable by such a disc is then increased from ∼ .
06 to ∼ . . . M ⋆ ), which is the case for our simulations.Mejia et al. (2005), using a constant cooling time, argue that globale ff ects might be present, but do not explicitly calculate such globaltorques. On the other hand, recent calculation by Boley et al. (2006)(see in particular their Fig. 13), which employ more realistic cool-ing properties, confirm that in the limit of small disc mass, thetransport induced by gravitational instabilities is essentially local.We have thus found that thermal history can a ff ect the abil- ity of the disc to exist in a self-regulated state without fragmen-tation but that this a ff ects the location of the stability boundary atonly the factor two level. The fact that there was negligible changein the fragmentation boundary when even slower changes in t cool were employed, demonstrates that thermal history is only part ofthe story. Our results suggest that, however slowly the disc is cooledthrough a sequence of thermal equilibria, there is still a fundamen-tal upper limit to the heating that can be provided by gravitationalinstabilities in a non-fragmenting disc. Thus it would appear thatthe initiation of fragmentation in a self-gravitating disc does not require that the disc enter the regime of rapid cooling on a short timescale. It is thus unnecessary to invoke sudden events (e.g. im-pulsive interactions with passing stars, Lodato et al. 2007) to tipa previously self-regulated disc into the fragmenting regime. In-stead our results suggest that fragmentation can in principle be ap-proached via the secular evolution of a self-gravitating disc. ACKNOWLEDGEMENTS
The simulations presented in this work have been performed atthe UK Astrophysical Fluid Facility (UKAFF). We thank RichardDurisen and Ken Rice for interesting discussions
REFERENCES
Balsara D. S., 1995, Journal of Computational Physics, 121, 357Bate M. R., Burkert A., 1997, MNRAS, 288, 1060Benz W., 1990, in Buchler J., ed., The Numerical Modeling of NonlinearStellar Pulsations Kluwer, DordrechtBoley A. C., Mej´ıa A. C., Durisen R. H., Cai K., Pickett M. K., D’AlessioP., 2006, apj, 651, 517Durisen R. H., Boss A. P., Mayer L., Nelson A. F., Quinn T., RiceW. K. M., 2007, in Reipurth B., Jewitt D., Keil K., eds, Protostars andPlanets V p. 607Gammie C. F., 2001, ApJ, 553, 174Lodato G., Meru F., Clarke C. J., Rice W. K. M., 2007, MNRAS, 374, 590Lodato G., Rice W. K. M., 2004, MNRAS, 351, 630Lodato G., Rice W. K. M., 2005, MNRAS, 358, 1489Mayer L., Lufkin G., Quinn T., Wadsley J., 2007, ApJ, 661, L77Mejia A. C., Durisen R. H., Pickett M. K., Cai K., 2005, ApJ, 619, 1098Monaghan J. J., 1992, ARA&A, 30, 543Nelson A. F., 2006, MNRAS, 373, 1039Rafikov R., 2005, ApJ, 621, 69Rafikov R. R., 2007, ApJ, 662, 642Rice W. K. M., Armitage P. J., Bate M. R., Bonnell I. A., 2003, MNRAS,338, 227Rice W. K. M., Lodato G., Armitage P. J., 2005, MNRAS, 364, L56Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337Stamatellos D., Whitworth A. P., Bisbas T., Goodwin S., 2007, ArXive-prints, 705c (cid:13)000