On the fragmentation boundary in magnetised self-gravitating discs

Accretion discs are prone to the gravitational instability if (Toomre, 1964)

where $c_{s}$ is the sound speed of the gas, $\Sigma$ is the gas surface density. The quantity $\kappa$ is the epicyclic frequency, which is equal to the angular frequency $\Omega$ if the disc is Keplerian. The ‘Toomre $Q$’ criterion (Eq. 1) refers to axisymmetric perturbations — non-axisymmetric perturbations are gravitationally unstable for $Q<1.5$–$1.7$ (Durisen et al., 2007). Unstable non-axisymmetric perturbations develop into spiral density waves thanks to the differential rotation of the disc (Goldreich & Lynden-Bell, 1965). These waves are trans-sonic, with Mach numbers of order unity (Cossins et al., 2009).

Depending on the available cooling mechanisms, the disc can become marginally stable (Paczynski, 1978). This is a self-regulated quasi-steady state where $Q$ is maintained near the stability limit, and the local weak shock heating through spiral density waves (increasing $Q$) approximately balances the local radiative cooling (reducing $Q$).

This quasi-steady state can be expressed in terms of the dimensionless cooling parameter $\beta_{c}$, defined according to

where $t_{cool}$ is the cooling time of the gas, and if the angular momentum transport is locally determined, by defining the dimensionless stress parameter $\alpha$ in terms of a pseudo-viscosity (Shakura & Sunyaev, 1973; Gammie, 2001):

where $H$ is the disc scale height. This approximates the heating of the disc as due to a turbulent pseudo-viscosity. This is an acceptable approximation if the spiral density wave modes couple and decay into gravito-turbulence and $Q\sim 1$ (Balbus & Papaloizou, 1999; Gammie, 2001). When the disc is marginally stable and in thermal equilibrium, there is a direct relationship between this viscous stress heating and the radiative cooling, and thus

where $\gamma$ is the ratio of specific heats. If the radiative cooling is sufficiently vigorous (or the accretion of material into the disc is sufficiently vigorous), then this marginally stable state cannot be maintained. This has been interpreted as a maximum stress that the gravitational instability can support, with $\alpha_{\rm max}\approx 0.06$ (Rice et al., 2005). The gravitational stress saturates at this maximum value, and if density perturbations can continue to grow (either via strong cooling or rapid mass loading), then these perturbations become unstable to gravitational collapse, resulting in disc fragmentation.

There remains debate as to the precise value of this maximum (Meru & Bate, 2011, 2012; Michael et al., 2012; Rice et al., 2012, 2014), and whether fragmentation can occur under low stresses, due to stochastic density fluctuations well above the RMS value that can occur given sufficient time (Paardekooper, 2012; Young & Clarke, 2015, 2016). It remains the case that fragmentation is favoured in discs that can no longer sustain thermal equilibrium through heating via spiral density waves and gravito-turbulence. Equivalently, disc fragmentation is favoured when the local Jeans mass inside a spiral density wave is decreasing rapidly (Forgan & Rice, 2011), either through decreases in the local sound speed, or increases in the local gas density in the case of irradiated discs (Kratter & Murray-Clay, 2011).

Magnetic fields have been absent from most studies of disc fragmentation to date. The usual argument for protostellar discs is that the low ionisation fractions mean that magnetic fields can be ignored. This is in stark contrast to studies of protostellar disc formation, which show that magnetic fields play a crucial role in extracting angular momentum from the system via magnetic braking. In ideal MHD simulations, this braking is so efficient that it can suppress disc formation entirely (Allen et al., 2003; Price & Bate, 2007; Hennebelle & Ciardi, 2009; Commerçon et al., 2010). Observations indicate that discs of order 100 au can exist around Class 0 protostars (e.g. Tobin et al., 2015), and the first direct observations of disc fragmentation proceed in this regime (Tobin et al., 2016). This apparent discrepancy between efficient magnetic braking and observed extended discs is resolved by either invoking turbulence, non-ideal MHD effects, or both (Seifried et al., 2013; Tomida et al., 2015; Wurster et al., 2016). Other surveys have shown that massive extended discs are perhaps less common (Maury et al., 2010), suggesting that some limited magnetic braking is occurring, and consequently limiting the frequency of disc fragmentation.

Discs around AGN are also likely to possess their own fields, through intense irradiation from the central engine and the surrounding environs. The presence of fields in both types of disc is evidenced by observations of jets and outflows, which are thought to be both launched and collimated by tightly wound field configurations (Pudritz et al., 2012; Marti-Vidal et al., 2015; Rodríguez-Kamenetzky et al., 2016).

There are few studies to date which consider the influence of magnetic fields on self-gravitating discs. During the same period that Gammie (2001) established equilibrium relationships for unmagnetised gravito-turbulence, Kim & Ostriker (2001) investigated magnetised self-gravitating gas in 2D shearing box simulations, determining that the magnetic Toomre $Q$

for non-axisymmetric perturbations to grow in the presence of a magnetic field (see also Kim et al. 2003). Linear stability analysis of magnetised self-gravitating discs (Lin, 2014) shows that the magneto-rotational instability (MRI) can be either stabilised or enhanced by gravitational instability (GI). If the resistivity is uniform, MRI modes with small radial lengths are stabilised. If the resistivity is a function of altitude from the midplane, MRI modes with radial length of order the disc thickness can be enhanced, which can significantly enhance axisymmetric density perturbations if the field is strong and toroidal. Interestingly, unstable modes can transition between MRI and GI, when the perturbations in gravitational potential and magnetic energy are comparable.

Fromang and collaborators (Fromang et al., 2004a; Fromang et al., 2004b; Fromang, 2005) performed 3D global simulations (assuming an isothermal disc), which suggested that weak magnetic fields (with a plasma parameter $\beta_{p}\sim 10^{3}$) would not prevent fragmentation. The growth of MRI induces turbulence, creating spiral modes that are out of phase with those produced by GI, modulating the gravitational stress and subsequently the star’s accretion rate. We will address here whether these features are due to the isothermal equation of state, or a generic feature of MRI-GI interaction in global discs.

More recently, local shearing box simulations were used to investigate the nature of 3D gravito-turbulence under the effect of weak ($\beta_{p}>>1$) and strong fields ($\beta_{p}=1$, Riols & Latter 2016). They identify three regimes, delineated by the energy balance between kinetic, gravitational and magnetic components. In the regime where magnetic energy is sub-dominant, a gravito-turbulent state can still be maintained, even with a relatively high Toomre-$Q$.

In the second regime, the magnetic energy exceeds the local gravitational energy (but not the kinetic energy). A gravito-turbulent state is still maintained, but the stresses in the system are now produced by magnetic fields, which use self-gravity as a dynamo. The weak turbulence maintained by gravitational instability is amplified by the formation of elongated current sheets, being dissipated as heat through magnetic reconnection. Interestingly, the deformation of the field results in “plasmoids”, magnetic islands which are remarkably similar to disc fragments, but created by fundamentally different physics. In the final regime, gravitational instability is completely suppressed by magnetic tension. Riols & Latter (2016) also studied the fragmentation boundary in the presence of magnetic fields. They find for a range of magnetic field strengths, the disc can still fragment if $\beta_{c}<5$.

To date, there is no corresponding study of magnetic fields on disc fragmentation using 3D global simulations at sufficient resolution. To this end, we conduct a series of numerical experiments on magnetised self-gravitating accretion discs using smoothed particle hydrodynamics (SPH). We adopt a simple parametrisation to initialise both the local cooling rate and the initial magnetic field configuration, which we describe along with the code and our analysis methods in §2. In §3 we describe how magnetic fields alter the structure of self-gravitating discs and the properties of fragments that they produce. In §4 we discuss the implications of our results, and in §5 we summarise the work.

Our principal result — that magnetic fields can suppress fragmentation for cooling times above five times the dynamical time— is in broad agreement with the local simulations conducted by Riols & Latter (2016). They also show fragmentation is suppressed for $\beta_{p}=8$ if $\beta_{c}>5$. We should however note that global and local simulations must be compared carefully. Our global simulations show fragmentation suppression due not only to direct action of the magnetic field, but also to redistribution of the disc material diluting the gravitational instability. Local simulations typically strictly conserve mass within the domain of interest, in effect preventing global redistribution of mass and energy as we observe.

What values of $\beta_{p}$ might we expect in real discs? Wurster et al. (2016) report values between 10-1000, from the inner and outer disc respectively. This is true both for their ideal and non-ideal MHD simulations of protostellar disc formation, and appears to be confirmed by other work (e.g. Tomida et al., 2015). This would suggest that disc fragmentation is not prevented by magnetic fields, provided the disc cooling rate is sufficiently strong locally. If the system retains enough angular momentum to form an extended disc, this seems quite possible (Forgan & Rice, 2012).

This suppression of fragmentation at what is quite low $\beta_{c}$, by typical values of $\beta_{p}$, is important to considerations of convergence in unmagnetised simulations. If the boundary for prompt fragmentation is $\beta_{c}\sim 7-8$ for zero field, and potentially higher for stochastic fragmentation, then the addition of what is quite a weak field reduces the critical $\beta_{c}$ significantly, and requires quite strong cooling for fragmentation to occur at all. This suggests that discussions of whether the prompt fragmentation boundary truly is $\beta_{c}=7,8$ or even 9 might not be particularly relevant to real self-gravitating discs, even if their magnetic fields are weak.

Our simulations indicate that there is both a minimum and maximum radius for fragmentation to occur when magnetic fields are active, but we should note that this could be an artifact of our idealised cooling. As $\beta_{c}$ declines with radius in realistic protostellar discs, we should therefore expect fragmentation to be further restricted to regions where $\beta_{c}$ is very low and $Q$ is maintained at the marginal stability limit, i.e. radii beyond 30 au. Further work is required to show whether there is a maximum radius for disc fragmentation, which will depend on the equilibrium disc structure produced by magnetised self-gravitating discs with realistic thermodynamics.

Our increased fragment masses are expected from analytic calculations of the magnetic Jeans mass (Strittmatter, 1966). For a medium with a uniform field along a single axis, a density perturbation will only be resisted along the two axes perpendicular to the field vector. The Jeans length $\lambda_{J}$ along these axes is modified thus (Chandrasekhar & Fermi, 1953):

A more accurate estimate of the fragment mass should consider the properties of the disc material inside a spiral density wave (Forgan & Rice, 2011). Future studies of expected initial fragment mass should consider the impact of magnetic fields in this manner.

We should be careful about extrapolating fragment behaviour from numerical experiments like this - the fragments are interacting with a disc structure unlikely to exist in real discs. In practice, we expect $\beta_{c}$ to increase with increasing proximity to the star, as the cooling time of the gas increases with increasing $\Sigma$. This being said, we do see rapid inward migration, which is a common feature of fragments in self-gravitating discs, not only due to the increased torques exerted on the fragment from a more massive disc (Baruteau et al., 2011), but also the fragment’s inability to open a gap and enter the slower Type II migration mode (Malik et al., 2015). We recommend CLUMPFIND analysis techniques be applied to simulations of fragmenting discs with realistic thermodynamics to investigate the dynamical properties of fragments further (see e.g. Hall et al, in prep).

Of particular importance to fragment survival is the evolution of spin angular momentum. Kim et al. (2003), in their study of giant molecular cloud formation in magnetised galactic discs, noted that the GMCs formed in the disc were subject to strong magnetic braking, losing a significant portion of their initial spin angular momentum. In fragments, this will have important consequences for the discs they may host (Forgan, 2016). It will also affect the fragment’s ability to collapse into bound objects, preventing Roche lobe overflow and tidal disruption (cf Boley & Durisen, 2010; Nayakshin, 2010; Zhu et al., 2012; Forgan & Rice, 2013).

We must note that none of our simulations resolve the magneto-rotational instability (MRI). In protostellar discs, this instability primarily acts in the inner disc, most especially in regions where artificial viscosity dominates our simulations. In the case of a global SPH disc simulation, we find that resolving the MRI is equivalent to another vertical resolution criterion, which in practice is computationally prohibitive. For the fastest growing MRI mode to be resolved, we must satisfy:

We do not satisfy this condition anywhere in our discs, even for relatively large $\beta_{p}$. As a result, we only resolve a fraction of the MRI modes present, and our description of the magnetic turbulence is incomplete. We therefore cannot comment on the possible interaction between MRI and GI, in particular whether GI modes can encourage the activation of MRI in the inner disc, which semi-analytic models predict sustains a limit cycle of episodic accretion (Armitage et al., 2001; Martin & Lubow, 2014). Conversely, simulations have already indicated that MRI modes can weaken the gravitational instability and prevent collapse (Fromang, 2005; Riols & Latter, 2016), in much the same way that fragmentation has been suppressed in our simulations. The fact that we find qualitative agreement with studies that fully resolve the MRI suggests that our results regarding disc fragmentation are robust.

Finally, we should acknowledge the critical role that non-ideal MHD effects may play in self-gravitating discs, especially in the coupling of gas and dust grains. Spiral density waves in unmagnetised discs have been shown to collect centimetre-metre size grains, in some cases showing significant overdensities (Rice et al., 2004; Dipierro et al., 2015). Whether such overdensities should be expected in observations of self-gravitating discs depends on the ability of grains to grow to sizes such that the Stokes number is of order unity. This in turn requires the velocity dispersion of the grains to remain sufficiently low during growth, which tends to restrict this effect to large radii (Booth & Clarke, 2016). This is where magnetic fields are most efficient at reducing gravitational stresses and spiral perturbation amplitudes, so magnetic fields may also suppress grain growth in self-gravitating discs.

Depending on the grain charge, non-ideal effects are likely to play an important role in mediating the separation of the gas and dust phases. In particular, Ohmic heating is likely to alter the structure of the spiral density wave and the weak shocks it induces.

We select sufficiently small values of $\beta_{c}$ that fragmentation occurs for all our discs in the absence of magnetic fields. We find that fragmentation in magnetised protostellar discs can still occur provided $\beta_{c}<5$, for both weak and strong fields (i.e. $\beta_{p}=1$). For $\beta_{c}=5$, fragmentation can be suppressed even by relatively weak fields ($\beta_{p}=10$).

Discs that fragment in the presence of magnetic fields produce fragments of increased mass compared to unmagnetised discs. Their orbital eccentricity and inclination are very similar to unmagnetised fragments, although our simulations show that there may be a minimum and maximum radius at which fragmentation occurs, as opposed to merely a minimum for unmagnetised discs.

We recommend that future studies of disc fragmentation must consider even weak magnetic fields, as they are likely to push the fragmentation zone even further from the central object, with obvious implications for planet/brown dwarf formation in protostellar discs, and star formation in discs around active galactic nuclei.

DHF and IAB gratefully acknowledge support from the ECOGAL project, grant agreement 291227, funded by the European Research Council under ERC-2011-ADG. DJP gratefully acknowledges funding via grants DP130102078 and FT130100034 and via Future Fellowship FT130100034 from the Australian Research Council. This work relied on the compute resources of the St Andrews MHD Cluster. Surface density plots were created using SPLASH (Price, 2007).