The Reaction of Self-Gravitating Protostellar Discs to Moderate Lessening In Cooling Time-Scale: the Fragmentation Boundary Revisited

Following the seminal work of Gammie (2001), there has been considerable progress in recent years in understanding the behaviour of self-gravitating accretion disc (see Durisen et al. 2007 and references therein). A number of simulations (Gammie 2001; Rice et al. 2003; Lodato & Rice 2004, 2005) have demonstrated that if the thermodynamic properties of the disc are evolved according to a thermal equation (involving a cooling term parametrized in terms of a cooling time- Toomre Q parameter: (where cs is the sound speed, κ is the epicyclic frequency (equal to the angular velocity _ in a Keplerian disc) and _ is the disc surface density) hovers at a value somewhat greater than unity over an extended region of the disc. Whereas the state Q = 1, corresponds to a situation of marginal stability against axisymmetric perturbations, in the self-regulated state the disc is instead subject to a variety of non-axisymmetric self-gravitating modes whose effect, through the action of weak shocks, is to dissipate mechanical energy (i.e. kinetic and potential energy of the accretion flow) as heat. Thermal equilibrium is then attained through the balancing of such heating by the prescribed radiative cooling: in essence, self-regulation results when the amplitude of these modes is able to self-adjust so as to maintain thermodynamic equilibrium against the relevant energyloss processes. Young stellar clusters are dynamic environments and encounters between cluster members can be quite common. Such encounters can have a significant effect on the structure and evolution of the gaseous discs that surround young stars (Bate, Bonnell & Bromm 2003). Some of these effects include a tidal truncation or disruption of the disc (Clarke & Pringle 1993; Hall, Clarke & Pringle 1996), a burst of accretion activity on to the central star (Bonnell & Bastien 1992), or the triggering of a gravitational instability in the disc, which might then lead to disc fragmentation (Boffin et al. 1998;Watkins et al. 1998a,b). The latter possibility has been often invoked as a possible formation mechanism for low-mass companions, such as brown dwarfs, or massive planets (for some recent reviews, see Goodwin et al. 2006; Whitworth et al. 2006). However, previous analyses of this process (Boffin et al. 1998; Watkins et al. 1998a,b) were limited by low resolution (which did not permit the proper resolution of the disc vertical structure) and by the replacement of the energy equation by a simple barotropic equation of state. All these analyses, in fact, assumed that the disc is isothermal, which might be a reasonable assumption at very large distances (the discs considered in these works are≈1000 au in size), but is certainly not adequate for smaller discs. Indeed, as clearly stated in Whitworth et al. (2006) ‘it is important that such simulations be repeated, with a proper treatment of the energy equation . . . to check whether low-mass companions can form at closer periastra’. The evolution of self-gravitating discs can therefore be neatly parametrised using the phenomenology of viscosity, much in the same way that the action of magneto-rotational instability (MRI) has been parametrised in less massive, more ionised discs (Blaes & Balbus 1994; Balbus & Papaloizou 1999). (Shakura & Sunyaev 1973), has been a useful tool in simplifying the computation of the structure of self-gravitating discs (Rafikov 2005; Clarke 2009; Rice & Armitage 2009), as well as computing the evolution of said discs (Lynden-Bell & Pringle 1974; Pringle 1981; Lin & Pringle 1990; Rice et al. 2010). Cold, dense molecular cloud cores are thought to be the sites of low-mass star formation (Terebey et al. 1984). Typically, these cores will contain an excess of angular momentum (compared with the rotational angular momentum of low-mass stars, cf Caselli et al. 2002). If these cores are sufficiently dense, they will collapse under their own gravity to form a protostar plus protostellar disc, with a substantial envelope surrounding both. In this pre-main sequence phase, we can expect that most of the star’s eventual mass will be processed by the protostellar disc (which in turn accretes this mass from a surrounding envelope). The star’s accretion is inextricably linked with the process of outward angular momentum transport - understanding the mechanisms by which angular momentum transport functions in these early phases is therefore essential to theories of star formation. While turbulence induced by the magneto-rotational instability (MRI) defines angular momentum transport at late times (Blaes & Balbus 1994; Balbus & Papaloizou 1999), the disc is unlikely to be even weakly ionised at early times. We must therefore seek another means by which to generate angular momentum transport. In this paper we conduct a suite of idealized simulations in which we explore whether the fragmentation boundary just depends on the instantaneous value of βcool (as has been assumed hitherto) of whether the system ‘remembers’ the history of how it 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 Rafikov 2005, 2007) which attempt to achieve ever-increasing verisimilitude via the incorporation of more realistic treatments of radiative transfer. Here, instead, we make no claims that the simplified cooling law (e.g. the assumption that βcool is spatially uniform) actually corresponds to a situation encountered in a real disc, because our aim is to isolate a particular physical effect (i.e. the time-scale on which the fragmentation boundary is approached). The computational expense of ‘realistic’ simulations however prevents their use to study secular effects: even in the case of the present ‘toy’ simulations, it is impracticable to run simulations over the long time-scales on which the profile changes due to gravitational torques or infall. We can nevertheless assess the effect of relatively slow changes in βcool on the fragmentation boundary through imposing an ad hoc reduction in the value of

Can the Jeans mass provide a timescale criterion for disc fragmentation (much as the cooling time has done in the past)? Let us now define three timescales, normalised to the orbital period: We use and to distinguish between variables of differing behaviour. is positive definite - this is to keep it in line with the conventional definition as it is currently used (and also to prevent taking imaginary values). measures the competition between disc accretion and stellar accretion. measures the timescale on which the local Jeans Mass changes - in orbital periods, we can expect the Jeans mass to either double or decrease to zero, depending on the sign of . and can be either positive or negative - for fragmentation to be favourable, the local Jeans mass must be decreasing. A small, negative therefore represents the most likely circumstances for disc fragmentation. The critical value of is less clear - any disc that can maintain a negative will proceed towards fragmentation. Rather than presenting us with two discrete regimes, fragmenting and non-fragmenting, we see a spectrum of potential outcomes, some fragmenting rapidly; some fragmenting on much longer timescales (possibly longer than the lifetime of the disc, and therefore effectively non-fragmenting); and others moving away from potential fragmentation (either slowly or rapidly). We suggest that −10 < < 0 gives highly favourable conditions for prompt fragmentation, but the lower limit is by no means fixed, and will require empirical confirmation. We wish to now derive : substituting for the self-gravitating scale height, we can rewrite the Jeans mass in terms of and : Instead of directly calculating , it is easier to describe the rate of change of specific internal energy u, which we can relate back to cs: We will assume that the sound speed varies due to two effects only - firstly, radiative cooling, given by and secondly, heating due to viscous dissipation, given by Using the chain rule on and rearranging gives: We now use equations, where we subsequently substitute back for MJ : Equivalently, This quickly gives the Jeans timescale as mass will decrease on a slower timescale, making fragmentation less likely. Also, small positive values of will rapidly increase the Jeans mass, suppressing fragmentation.

The SPH code: Our three-dimensional numerical simulations are carried out using smoothed particle hydrodynamics (SPH), a Lagrangian hydrodynamic scheme (Benz 1990; Monaghan 1992). The general implementation is very similar to Lodato & Rice (2004, 2005) and Rice et al. (2005). The gas disc is modelled with 250 000 SPH particles (500 000 in a run used as a convergence test) and the local fluid properties are computed by suitably averaging over the neighbouring particles. The disc is set in almost Keplerian rotation (allowing from slight departures from it to account for the effect of pressure forces and of the disc gravitational force) around a central point mass on to which gas particles can accrete if they get closer than the 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 viscosity. The ratio of specific heats is = 5/3. Cooling is here implemented in a simplified way, i.e. by parameterizing the cooling rate in terms of a cooling time-scale: where ui is the internal energy of a particle and the cooling timescale tcool is assumed to be proportional to the dynamical time-scale, , where βcool is varied according to a time-dependent prescription. Artificial viscosity is introduced using the standard SPH formalism. The actual implementation is very similar to the one used in Rice et al. (2005), that is we set the two relevant numerical parameters to and and we have not included here (consistent with Rice et al. 2005) the so-called Balsara switch (Balsara 1995) to reduce shear viscosity. Disc set-up : The main physical properties of the disc at the beginning of the simulation are again similar to those of Lodato&Rice (2004, 2005). The disc surface density is initially proportional to R−1 (where R is the cylindrical radius), while the temperature is initially proportional to R−1/2. Given our simplified form of the cooling function, the computations described here are essentially scale-free and can be rescaled to different disc sizes and masses. For reference, we will assume that the unit mass (which is the mass of the central star) is and that the unit radius is 1 au. In this generally chosen such as to have a total disc mass of Mdisc = , while the temperature normalization is chosen so as to have a minimum value of Q = 2, which is attained at the outer edge of the disc. Initially, the disc is evolved with constant βcool = 7.5, this value of βcool being in the regime where previous work Gammie (2001), Rice et al. (2005) has shown that the disc does not fragment. The general features of this initial evolution is described in detail in Lodato&Rice (2004). The disc starts cooling down until the vertical scalelength H is reduced such that . At this point the disc becomes Toomre unstable and develops a spiral structure that heats up the disc and maintains it close to marginal stability. We have evolved the disc with this value of βcool for 7.8 outer disc orbits. At this stage it is close to Q = 1 over most of the disc (i.e. over the radial range R = 3–23 au). Adding the perturber : Once the disc has reached the quasi-steady, self-regulated state of marginal stability, we add the perturber star. This is initially put at a distance of 100 au from the primary star, and is placed in a parabolic orbit. At this distance, the sudden addition of the perturber does not cause any significant perturbation to the disc structure. We have performed a large number of simulations, by varying the orbital parameters, considering coplanar encounters (both prograde and retrograde with respect to the disc rotation) and ‘head on’ encounters, where the orbital plane of the secondary is perpendicular to the disc plane.We have also performed a simulation of a hyperbolic (coplanar) encounter. We have varied the disc mass, the perturber mass Mpert and the pericentre of the perturber orbit rperi and we have also run simulations with different numerical resolutions.

We find that in the fast case (x = 10.5), a fragment forms when βcool = 0.75, i.e. about 8.9 outer disc dynamical times after the rapid reduction in βcool commenced. Since fragmentation always takes about a dynamical time-scale 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 different behaviour in the slow case: here we find that when βcool is reduced to, and then held at, βhold = 2.75, the disc does not fragment even when the disc is then integrated for a further 44 outer dynamical time-scales. Likewise, for the slowcase, the disc does not fragment when held at βhold = 3, even after integration for 63 outer dynamical time-scales at this βcool value.1 On the other hand, when βcool was instead held at 2.62, it fragmented after a furtherouter dynamical time-scales, so it would appear that the fragmentation boundary is at around

imposed. We hesitate to say that we have proved that a disc will never fragment when brought to such a low value of βcool at this slow rate, since our experience shows that where one is close to the limit of marginally stable βcool, fragmentation may ensue after long time-scales, and that its timing may depend on numerical noise that can be affected by resolution. Indeed, we found that when we reran the x = 105 simulation at higher resolution (N = 500 000), and held it at βhold = 3, it eventually did form a fragment at large radius.We however show the disc structure in this simulation at the point of fragmentation and contrast it with the corresponding situation when the cooling time is rapidly reduced and then held at constant βcool = 3 (i.e. model F2). Evidently, notwithstanding the fact that a fragment does eventually form in the former case also, the disc structure is quite different in the two cases, with the ‘rapid’ simulation containing a number of regions that are on the point of fragmentation at the moment that the first fragment appears. Our interest here is not in defining precise boundaries at which fragmentation will or will not occur (since the definition of such a boundary is always contingent on the duration of the simulation) but in demonstrating that the structure of the disc is indeed affected not just by the instantaneous value of βcool (and hence on the heating rate that has to be delivered through the action of the self-gravitating modes) but also on the history of how the disc arrived at such a value of βcool. As has been found many times previously, discs with outer radii less than around 40 au will typically not form fragments (Rafikov 2005; Matzner & Levin 2005; Whitworth & Stamatellos 2006; Boley et al. 2006; Stamatellos & Whitworth 2008; Forgan et al. 2009; Clarke 2009). While this is more a symptom of the imposed disc model (which in itself is not a particularly new implementation), we can use this result to satisfy ourselves that the model is indeed correctly distinguishing the fragmentation boundary. The dependence of the radial fragmentation boundary on stellar mass (i.e. lower mass stars have a lower boundary) is also an indication that the model is indeed performing as expected.

We have found that the rate at which the cooling time-scale is changed indeed affects the minimum value of βcool at which the disc can exist in a stable, self-regulated state. As expected, this effect is only manifest when the cooling time-scale is varied on a time-scale (τ ) that is longer than the cooling time-scale, since for τ < tcool, the temperature always falls on a time-scale tcool, irrespective of τ . less than those that are possible when a fixed cooling time-scale is imposed at the outset of the simulation. This implies that (in the slow-cooling case) the gravitational instabilities are able to deliver about twice the heating rate without the disc fragmenting. In terms of the ‘viscous alpha’ description of such instabilities (Shakura&Sunyaev 1973; Gammie 2001; Lodato & Rice 2005), the maximum deliverable by such a disc is then increased from . It should be noted that such ‘local’ description of the transport induced by gravitational instabilities is only possible in the limit in which global, wave-like transport does not play an important role. Lodato & Rice (2004, 2005), using a cooling prescription similar to ours, have shown that this is the case, as long as the total disc mass is small ( 0.2M*), which is the case for our simulations. Mejia et al. (2005), using a constant cooling time, argue that global effects might be present, but do not explicitly calculate such global torques. On the other hand, recent calculation by Boley et al. (2006), which employ more realistic cooling properties, confirm that in the limit of small disc mass, the transport induced by gravitational instabilities is essentially local. In this paper, we have analysed the effects of an encounter with a discless companion on the structure and evolution of a protostellar disc. We have improved on previous analyses by considering in more detail the energy balance for the gas, thus going beyond the approximation of isothermal evolution used in the past. This turns out to be very important, because the outcome of the encounter is strongly dependent on the tidal heating induced in the disc. Indeed, probably the main result of the present paper is the demonstration that an encounter with a companion inhibits rather than promotes fragmentation of a gravitationally unstable disc. This is in marked contrast with previous results (Boffin et al. 1998), who instead have shown effective fragmentation following an encounter. The main difference between our simulations and those of Boffin et al. (1998) is that while Boffin et al. (1998) consider much larger discs, for which the isothermal approximation is probably adequate, we instead follow the heating and cooling processes in the disc and therefore allow the growth of gravitational instabilities to feed back upon the thermodynamic state of the disc. We include heating from pdV work and shocks and cool the disc down with a cooling rate that is sufficiently small that the disc would not fragment in isolation. We have derived the local Jeans mass inside a spiral perturbation, and shown that it can be calculable using azimuthally averaged disc variables. This allows the expected mass of objects formed from disc fragmentation to be calculated simply, both in expression and derived a Jeans timescale , which describes the length of time required for the Jeans Mass to increase or decrease by its current value (normalised by the local angular velocity, in much the same way as the cooling time is often normalised). We show that the resulting expression for reduces to the standard cooling time criterion as was previously defined (Gammie 2001; Rice et al. 2005), and folds in related results on the disc’s thermal history (Clarke et al. 2007) and envelope accretion (Kratter et al. 2010). We confirm that subjecting the disc to extra shock heating promotes fragmentation, and that rapidly accretion encourages discs to fragment, provided that local angular momentum transport and cooling is efficient.

• Rice W. K. M., Armitage P. J., Bate M. R., Bonnell I. A., 2003, MNRAS, 338, 227

• Boley A. C., Mej´ıa A. C., Durisen R. H., Cai K., Pickett M. K., D’Alessio P., 2006, ApJ, 651, 517

• Stamatellos D., Whitworth A. P., Bisbas T., Goodwin S., 2007, ArXiv e-prints, 705

• Durisen 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 and Planets V.

• University of Arizona Press, Tuscon p. 607

• Zhu Z., Hartmann L., Gammie C., McKinney J. C., 2009, ApJ, 701, 620

• Goodwin S. P., Kroupa P., Goodman A., Burkert A., 2006, ArXiv Astrophysics e-prints

• Watkins S. J., Bhattal A. S., Boffin H. M. J., Francis N., Whitworth A. P.,

• 1998a, MNRAS, 300, 1205

• Mejia A. C., Durisen R. H., Pickett M. K., Cai K., 2005, ApJ, 619, 1098

• Cai K., Durisen R., Boley A. C., PickettM. K.,Mejia A. C., 2008, ApJ, 673, 1138

• Rice W. K. M., Lodato G., Pringle J. E., Armitage P. J., Bonnell I. A., 2004, MNRAS, 355, 543

• Walch S., Naab T., Whitworth A., Burkert A., Gritschneder M., 2010, MNRAS, 402, 2253

• Rice W. K. M., Mayo J. H., Armitage P. J., 2010, MNRAS, 402, 1740

• Kratter K. M., Matzner C. D., Krumholz M. R., 2008, ApJ, 681, 375