## Abstract

In this paper we perform a systematic analysis of the structure induced by the onset of gravitational instabilities in cooling gaseous accretion discs. It is well known that for low enough cooling rates the disc reaches a quasi-steady configuration, with the instability saturating at a finite amplitude such that the disc is kept close to marginal stability. We analyse the dependence of the saturation amplitude on the imposed cooling rate, and we find that it scales with the inverse square root of the cooling parameter β=*t*_{cool}/*t*_{dyn}. This indicates that the heating rate induced by the instability is proportional to the energy density of the density waves excited by the disc self-gravity. In particular, we find that at saturation the energy dissipated per dynamical time by weak shocks due to the gravitational perturbations is of the order of 20 per cent of the wave energy. We further perform a Fourier analysis of the disc structure, and subsequently determine the dominant radial and azimuthal wavenumbers of the density waves. While the number of spiral arms (corresponding to the azimuthal wavenumber) is fairly constant with radius, we find that the disc displays a range of radial wavenumbers whose mean increases with increasing radius. The dominant modes closely match the locally most unstable wavelength as predicted by linear perturbation analysis. As a consequence, we numerically demonstrate that the density waves excited in relatively low-mass discs *M*_{disc}/*M*_{*}∼ 0.1 are always close to corotation, deviating from it by approximately 10 per cent. This result can be understood in terms of the constancy of the Doppler-shifted phase Mach number of the flow – the pattern speed self-adjusts so that the flow into spiral arms is always sonic. This has profound effects on the degree to which the extraction of energy and angular momentum from the mean flow through density waves can be modelled as a viscous process. Our results thus provide (i) a detailed description of how the self-regulation mechanism is established for low cooling rates, (ii) a clarification of the conditions required for describing the transport induced by self-gravity through an effective viscosity, (iii) an estimate of the maximum amplitude of the density perturbation before fragmentation takes place and finally (iv) a simple recipe to estimate the density perturbation in different thermal regimes.

## 1 INTRODUCTION

In cold relatively massive accretion discs, the effects of the disc self-gravity can become dynamically important, and gravitational instabilities may play a major role in determining the long-term evolution of the disc. On the one hand, the instability can induce fragmentation of the disc, potentially resulting in the formation of massive stars in active galactic nuclei accretion discs (Nayakshin, Cuadra & Springel 2007) or, in the case of protoplanetary discs, low-mass companions (Stamatellos, Hubber & Whitworth 2007) or even giant planets (Boss 1997, 1998). On the other hand, the instability is very effective in transporting and redistributing angular momentum within the disc (Lodato & Rice 2004, 2005) and might therefore promote accretion in a variety of different contexts and scales, from supermassive black hole growth (Shlosman, Begelman & Frank 1990; Lodato & Natarajan 2006) to star formation (Vorobyov & Basu 2005).

The balance between the heating provided through such gravitational instabilities and the cooling as energy is radiated away is known to have major implications for the fate of the disc. Fragmentation into bound objects occurs when the cooling dominates, whereas the disc can persist in a quasi-equilibrium marginally stable condition if the cooling rate is sufficiently low (Gammie 2001; Lodato & Rice 2004; Rice, Lodato & Armitage 2005).

Such a quasi-steady state is characterized by the presence of spiral density waves, which transport energy and angular momentum through the disc. Where the discs are weakly ionized, and therefore where the magneto-rotational instability is likely to be ineffective, this angular momentum transport mechanism may be the primary driver of accretion. Additionally, these waves extract rotational energy from the disc, some of which is then returned as heat as the waves steepen into shocks. This heat then stabilizes the disc against further gravitational collapse into bound fragments. As the amplitude of the waves increases, so also does the wave energy density, increasing the reservoir of energy available to be returned to the disc as heat. However, numerical experiments (Gammie 2001; Rice et al. 2005) have shown that if the cooling time *t*_{cool} is less than a few times the local dynamical time Ω^{−1}, this feedback process breaks down and fragmentation ensues. The quasi-steady marginal stability state therefore represents a restricted regime of dynamic thermal equilibrium, where the heating through gravitational instabilities is balanced by the cooling rate.

In this paper, we seek to characterize the relationship between the strength of the cooling and the amplitude of the spiral density waves excited within the disc through self-gravity, while remaining within this dynamic thermal equilibrium state. To this end, we use a smoothed particle hydrodynamics (sph) code to run global numerical simulations of self-gravitating gaseous discs. From such controlled numerical experiments, we measure the amplitude of the density perturbations over a range of cooling times, down to the limit where the disc fragments, and investigate the overall structure formation. The Fourier analysis allows us to characterize the mode spectra and pattern speeds associated with this structure, and to associate the dynamics of the spiral density waves excited through self-gravity with the thermodynamics of the disc self-regulation process. We begin by discussing some relevant properties of density waves and transport processes in gaseous discs in Section 2, and link this to the disc thermodynamics in Section 3. Details of the numerical simulations and initial conditions are given in Section 4, and the results of these simulations are presented in Section 5. In Section 6, we discuss the implications of these results and the conclusions that may be drawn from them.

## 2 DYNAMICS OF SELF-GRAVITATING GASEOUS DISCS

In the following two sections, we analytically derive the basic relations that link the disc quantities in the presence of perturbations to the induced transport properties of self-gravitating discs. These standard results have been mostly derived in the context of stellar dynamics (Shu 1970; Bertin 2000; Binney & Tremaine 2008). In the context of gaseous (accretion) discs, qualitatively similar (but rather more involved) analyses have been presented in Balbus & Papaloizou (1999) and Balbus (2003).

### 2.1 Dispersion relations for self-gravitating discs

The linear development of the gravitational instability is usually described by the dispersion relation *D*(ω, *k*, *m*) between the wave (angular) frequency ω and the radial and azimuthal wavenumbers of excited modes, *k* and *m*, respectively. For an infinitesimally thin (i.e. two-dimensional) disc, this can be obtained using the standard WKB method of perturbation analysis, and in the tight-winding limit (where the radial wavelength is small in comparison with the azimuthal wavelength, |*kR*| ≫*m*, with *R* as the cylindrical radius), we have

*R*) is the angular speed, κ(

*R*) is the epicyclic frequency,

*c*

_{s}(

*R*) is the sound speed and Σ(

*R*) is the surface density of the unperturbed disc (Binney & Tremaine 2008). Introducing the pattern speed of the wave Ω

_{p}(

*R*), where ω=

*m*Ω

_{p}, we may express this as follows:

The well-known stability criterion for axisymmetric disturbances (Toomre 1964)

identifies the region of parameter space where the right-hand side (RHS) of equation (2) is positive definite, and therefore where the disc is stable at all wavelengths. For the case of an unstable disc, the most unstable wavelength*k*

_{uns}(i.e. where the RHS of equation 2 is at a minimum) occurs at a radial wavenumber For a marginally stable disc, where

*Q*≈ 1, we would therefore expect only modes with

*k*≈

*k*

_{uns}to be excited, and in this instance, equation (2) tells that Ω

_{p}≈Ω, i.e. all excited modes are expected to be close to corotation. Note that the most unstable wavenumber is exactly the inverse of the disc thickness for a self-gravitating disc where the last approximation is exact for a Keplerian rotation curve.

### 2.2 The stress tensor

In non-axisymmetric self-gravitating discs, torques arising due to the perturbed gravitational potential play an important role in transporting both energy and angular momentum.

Traditionally, viscous accretion discs are described by the α-formalism of Shakura & Sunyaev (1973), where (for an infinitesimally thin disc) the only non-vanishing component of the vertically integrated stress tensor * T* is the shear term, defined as

*c*

^{2}

_{s}, and indicates that the α-formalism is fundamentally a local relationship. Furthermore, since α is not necessarily constant, this represents a completely general description for any purely local process. Also note that for Keplerian rotation, d ln Ω/d ln

*R*=−3/2, implying that the stress is negative, i.e. it acts to oppose rotation, and therefore allows for inwards accretion flows.

The accretion disc theory identifies the origin of the stress with torques arising due to perturbations in a ‘turbulent’ disc. Crucially for the gas dynamics, these perturbations manifest themselves as fluctuations in the mean flow velocity, the gravitational potential and the magnetic field threading the disc. For non-magnetized self-gravitating discs such as we consider here, the stress tensor can therefore be broken down into a Reynolds stress term, associated with velocity fluctuations, and a gravitational stress term, associated with fluctuations in the gravitational potential. The Reynolds stress term *T*^{Reyn}_{Rφ} is such that

*v*

_{R}and δ

*v*

_{φ}are the velocity fluctuations about the mean flow velocity in the

*R*and φ directions, respectively, and the brackets indicate azimuthal averaging) and the gravitational stress term

*T*

^{grav}

_{Rφ}is given by Lynden-Bell & Kalnajs (1972) as where again

*g*

_{R}and

*g*

_{φ}are the accelerations due to the perturbed gravitational potential of the disc in the

*R*and φ directions, respectively.

The viscous torque per unit area is related to the vertically integrated stress tensor * T* through

### 2.3 Wave energy and angular momentum densities

In this section, we turn our attention to the transport of energy and angular momentum through the propagation of the spiral density waves expected to arise in a self-gravitating disc. To this end, it is convenient to introduce the wave action density for density waves. Within the WKB approximation used to derive the dispersion relation given in equation (2), the wave action per unit surface is given by (Toomre 1969; Shu 1970; Fan & Lou 1999)

where δΦ is the perturbed gravitational potential, itself given by where δΣ is the surface density perturbation. The wave energy per unit surface and the wave angular momentum per unit surface are obtained in a straightforward way from the wave action through the standard wave dynamics relations (Shu 1970; Bertin 2000):1Combining the first of these relations with equations (12) and (13), we obtain

where are the radial phase speed and Doppler-shifted radial phase speed of the wave, respectively. Note that equations (16) and (18) together explain why self-induced density waves are launched at corotation – since the energy density changes sign at corotation, waves that propagate away from corotation extract no net energy from the flow.Looking again at equations (14) and (15), we see that the relationship between the energy per unit surface and the angular momentum per unit surface in a density wave is given by

In the case of quasi-stationary waves (where Ω_{p} is constant) propagating in a disc in dynamic thermal equilibrium, the rate at which energy is lost per unit surface due to cooling must be balanced by the power per unit surface dissipated by the waves. In order to maintain the amplitude of the wave, the instability has to keep extracting energy and angular momentum from the background flow. The fluxes of energy (angular momentum) carried by the wave are simply times the local group velocity. Hence, when a wave dissipates, it adds energy and angular momentum to the flow in the ratio of to , i.e. in the ratio Ω_{p}. We therefore conclude that

Equation (20) is the wave analogue of equation (11). Comparing these two equations, we note a fundamental difference with respect to the viscous model – for a given torque , waves extract energy from the flow at a rate proportional to the wave pattern speed Ω_{p}, whereas the rotation speed Ω is the underlying rate in the local (viscous) case.

Balbus & Papaloizou (1999) have noted similarly that in general, energy transport through the gravitational instability cannot be described purely in viscous terms, and indeed this is only possible at corotation, when Ω_{p}=Ω. This can be readily understood if we consider that the wave energy per unit surface (equation 16) can be decomposed into two separate terms, as follows:

_{p}−Ω, is a non-local term. In fact the energy flux associated with this non-local transport term is precisely that identified by Balbus & Papaloizou (1999) as an ‘anomalous flux’, preventing self-gravitating discs from acting as pure α-discs.

We see that far from corotation, where Ω≠Ω_{p}, this term becomes significant, and thus non-local (global) transport becomes important within the disc. For trailing waves launched at corotation, the direction of energy (and angular momentum) transport is outwards throughout the disc (i.e. inwards travelling waves with negative energy density inside corotation and outwards travelling waves with positive energy density outside corotation), and thus the dissipation of such waves effects a net outwards transport of energy (and angular momentum). If such a wave dissipates at large radius (where Ω_{p}≫Ω), then the ratio in which energy and angular momentum are added to the disc (equation 20) is significantly greater than the equivalent ratio in the viscous case (equation 11). Consequently, under such conditions, the energy dissipated at large radii in a steady state disc with wave transport can significantly exceed that dissipated in an equivalent viscous disc – with this extra energy being extracted by the wave from deep in the potential and transported to large radii (Bertin & Lodato 2001; Lodato & Bertin 2001). However, if waves instead dissipate close to corotation, the wave transport is dominated by the local term in equation (21); since the energy and angular momentum transport are exchanged with the disc in roughly the same ratio as for a viscous process, then in this regime the α-formalism is a good approximation to the actual transport properties of the disc.

From equation (21) it is therefore possible to quantify a non-local transport fraction ξ, from the ratio of the two terms on the RHS, such that

Thus in order to assess the importance of non-local effects, we need to know the relationship between the angular frequency Ω and the pattern speed Ω_{p}. In Section 5.5, we use the dispersion relation along with the information extracted through the Fourier analysis in order to estimate the pattern speed, and hence to evaluate ξ directly.

## 3 SIMULATING THE DISC THERMODYNAMICS

Realistically simulating the thermodynamics of accretion discs is a complex undertaking and as such has received much attention, from the opacity-based treatment employed by Johnson & Gammie (2003) through the various convective and radiative transfer models of Boss (2004), Boley et al. (2007), Mayer et al. (2007), Stamatellos & Whitworth (2008a) and Stamatellos & Whitworth (2008b), the latter two of which also account for heating from the central star.

However, in this paper, one of our aims is to investigate the relationship between the properties of the density perturbations and the rate at which the disc cools. This purpose is served most readily by *imposing* a known cooling rate against which the heating rate and subsequent disc structure may be easily correlated. It is therefore not necessary to consider the exact physics of the cooling regimes found in astrophysical discs, and hence we can use a cooling law for the heat loss rate per unit mass *Q*^{−}, such that

*u*is the internal energy per unit mass and where the details of the cooling function (and possibly of additional

*external*heating) are absorbed into the simple parameter

*t*

_{cool}. As long as such a characteristic time-scale can be defined, it is therefore possible to use this formalism to represent a wide range of cooling mechanisms. Within this paper, we use a fixed ratio between the local dynamical and cooling time-scales, such that β=Ω

*t*

_{cool}is constant. This form of cooling has been used extensively in simulations of discs in various contexts, e.g. Gammie (2001), Lodato & Rice (2005) and Hobbs & Nayakshin (2008), and has proven useful in elucidating the properties of the gravitational instability in controlled numerical experiments.

In terms of the heating from the gravitational instability, we noted in the previous section that density waves extract energy from the disc. In addition to compression heating, in the case where the pattern speed differs from the rotation speed by more than the local sound speed, these waves will steepen into shocks, liberating further heat into the disc. We expect the rate at which energy is added to the disc to scale with the energy of the wave and the local dynamical time-scale, and we can then express the heating rate per unit mass due to the instability as

where we define the radial phase and Doppler-shifted phase Mach numbers to be and , respectively. In equation (24), we have also introduced a dimensionless proportionality factor ε, hereinafter referred to as the heating factor. If the relationship between the pattern speed and the angular speed is self-similar (i.e. it does not vary across the disc), we expect ε to be simply a constant, independent of radius.Once the gravitational instability has been instigated and has subsequently saturated, we may consider the disc to be in dynamic thermal equilibrium, such that the energy released through wave-driven shock and compression heating is balanced by the imposed cooling, i.e. *Q*^{−}+*Q*^{+}= 0. Recalling that *u*=*c*^{2}_{s}/γ(γ− 1), we can equate equations (23) and (24) and thereby determine the following relationship between the amplitude of the density perturbations and strength of the cooling, as measured by the β parameter:

In this paper, we will therefore use global numerical simulations to test the above energy balance and to investigate the relative magnitude of the local and non-local transport terms.

## 4 NUMERICAL SETUP

### 4.1 The sph code

All of the simulations presented hereafter were performed using a three-dimensional sph code and a Lagrangian hydrodynamics code capable of modelling self-gravity (see e.g. Benz 1990 and Monaghan 1992). All particles evolve according to individual time-steps governed by the Courant condition, a force condition (Monaghan 1992) and an integrator limit (Bate, Bonnell & Price 1995).

We note here that with sph, the integral of a physical quantity *A* over a given volume *V* is estimated by the sum over the individual particle values of this quantity, as below:

*m*

_{i}is the particle mass, and

*i*loops over all the particles within the volume

*V*. In a similar manner, we note that a volume-averaged value for

*A*, which we will call , can therefore be estimated via when, as in all our simulations, all particles have the same mass.

We have modelled our disc systems as a single point mass (on to which gas particles may accrete if they enter within a given sink radius, and satisfy certain boundness conditions – see Bate et al. 1995) orbited by 500 000 sph gas particles; a setup common to many other sph simulations of such systems, (e.g. Rice et al. 2003; Lodato & Rice 2004, 2005; Clarke, Harper-Clark & Lodato 2007) but with increased resolution. The central object is free to move under the gravitational influence of the disc. In order to ensure that the simulations were properly converged, resolution checks were undertaken with discs consisting of both 250 000 and 1000 000 particles – these are discussed briefly in Appendix A.

As described in Section 3 we use a simple cooling model, implemented in the following manner:

where the*u*

_{i}and

*t*

_{cool,i}are the internal energy per unit mass and the cooling time associated with each particle, respectively. Again as above, the functional form of the cooling time is kept simple, such that Ω

_{i}

*t*

_{cool,i}=β, where Ω

_{i}is the angular velocity of each particle, and where β is held constant throughout any particular simulation. All simulations have been run modelling the particles as a perfect gas, with the ratio of specific heats γ= 5/3, heat addition being allowed for via

*P*d

*V*work and shock heating and with the cooling implemented as specified above. Artificial viscosity has been included through the standard sph formalism, with α

_{SPH}= 0.1 and β

_{SPH}= 0.2. Note that these values are smaller than those commonly used in sph simulations; we use these values to limit the transport induced by artificial viscosity. As shown in Lodato & Rice (2004), with this choice of parameters, the transport of energy and angular momentum due to artificial viscosity is a factor of 10 smaller than that due to gravitational perturbations, while we are still able to resolve the weak shocks occurring in our simulations.

By using the cooling prescription outlined above, the rate at which the disc cools is governed by the dimensionless parameter β and the cooling is thus implemented scale-free. The governing equations of the entire simulation can likewise be recast in dimensionless form. In common with the previous sph simulations mentioned above, we define the unit mass to be that of the central object – the total disc mass and individual particle masses are therefore expressed as fractions of the central object mass. We can self-consistently define an arbitrary unit (cylindrical) radius *R*_{0}, and thus, with *G*= 1, the unit time is the dynamical time *t*_{dyn}=Ω^{−1} at radius *R*= 1.

### 4.2 Initial conditions

All our simulations model a central point object of unit mass *M*_{*}= 1, surrounded by a gaseous disc of mass *M*_{disc}. Although the bulk of the simulations have been conducted with a disc to central object mass ratio *q*=*M*_{*}/*M*_{disc} of 0.1, simulations have been run with *q*= 0.05, 0.075 and 0.125 to investigate the effects of the mass ratio on the non-local energy transport fraction ξ.

All simulations have used an initial mass surface density profile Σ∼*R*^{−3/2}, which implies that in the marginally stable state where *Q*≈ 1, the disc temperature profile should be approximately flat. Since the surface density evolves on the viscous time *t*_{visc}≫*t*_{dyn}=Ω^{−1}, this profile remains roughly unchanged throughout the simulations. The initial temperature profile is *c*^{2}_{s}∼*R*^{−1/2} and is such that the minimum value of the Toomre parameter *Q*_{min}= 2 occurs at the outer edge of the disc. In this manner, the disc is initially gravitationally stable throughout. Note that the disc is *not* initially in thermal equilibrium – heat is not input to the disc until gravitational instabilities are initiated.

Radially the disc extends from *R*_{in}= 0.25 to *R*_{out}= 25.0, as measured in the code units described above. The disc is initially in approximate hydrostatic equilibrium in a Gaussian distribution of particles with scaleheight *H*. The azimuthal velocities take into account both a pressure correction (Lodato 2007) and the enclosed disc mass. In both cases, any variation from dynamical equilibrium is washed out on the dynamical time-scale.

Given the dimensions above, one outer dynamical time-scale of the disc corresponds to 125 time units. To ensure that thermal equilibrium is reached and that the gravitational instability is saturated, all simulations are followed for at least 10 outer cooling times. To this end we will refer to the thermal time *t*_{therm} for each simulation as the cooling time evaluated at the initial outer edge of disc, taken to be at *R*= 25– thus *t*_{therm}=*t*_{cool}(25) =βΩ^{−1}(25).

### 4.3 Simulations run

In all a total of nine distinct simulations were run for various values of the cooling parameter β and the disc to central object mass ratio *q*, as detailed in Table 1. Although previous investigations with Σ∼*R*^{−1} have found that the fragmentation boundary is at β_{frag}≈ 6, (Rice et al. 2003, 2005), we find that in the case where Σ∼*R*^{−3/2}, the fragmentation boundary is slightly different, with 4 < β_{frag} < 5. The simulation where β= 4 therefore contains a fragment, and is included primarily for completeness. All results henceforth are given at the time quoted in the final column of Table 1 unless otherwise stated. The raw data are time-averaged over 500 unit times about these values to enhance the signal-to-noise ratio and to give the approximate steady-state values.

β | q=M_{disc}/M_{*} | No. of particles | Duration |

4 | 0.10 | 500 000 | 4.0 t_{therm} |

5 | 0.10 | 500 000 | 10.0 t_{therm} |

6 | 0.10 | 500 000 | 10.0 t_{therm} |

7 | 0.10 | 500 000 | 10.0 t_{therm} |

8 | 0.10 | 500 000 | 10.0 t_{therm} |

9 | 0.10 | 500 000 | 10.0 t_{therm} |

10 | 0.10 | 500 000 | 10.0 t_{therm} |

5 | 0.050 | 500 000 | 10.0 t_{therm} |

5 | 0.075 | 500 000 | 10.0 t_{therm} |

5 | 0.100 | 500 000 | 10.0 t_{therm} |

5 | 0.125 | 500 000 | 10.0 t_{therm} |

β | q=M_{disc}/M_{*} | No. of particles | Duration |

4 | 0.10 | 500 000 | 4.0 t_{therm} |

5 | 0.10 | 500 000 | 10.0 t_{therm} |

6 | 0.10 | 500 000 | 10.0 t_{therm} |

7 | 0.10 | 500 000 | 10.0 t_{therm} |

8 | 0.10 | 500 000 | 10.0 t_{therm} |

9 | 0.10 | 500 000 | 10.0 t_{therm} |

10 | 0.10 | 500 000 | 10.0 t_{therm} |

5 | 0.050 | 500 000 | 10.0 t_{therm} |

5 | 0.075 | 500 000 | 10.0 t_{therm} |

5 | 0.100 | 500 000 | 10.0 t_{therm} |

5 | 0.125 | 500 000 | 10.0 t_{therm} |

## 5 SIMULATION RESULTS

Common to all the simulations run is an initial phase in which the discs cool rapidly until the value of *Q* becomes approximately unity, at which point the gravitational instability is initiated and heat is liberated to balance the cooling. This stage is complete after approximately one thermal time, and from then on the discs settle into a quasi-steady state characterized by the presence of spiral arms throughout almost their entire radial range, with *Q*≈ 1. The quasi-static *Q* profiles to which the discs converge are shown in Fig. 1, with the cooling parameter β varying in the top panel, and the disc to central object mass ratio *q* varying in the bottom panel. Note that the data are plotted at the times given in Table 1. Throughout all the simulations, it can be seen that the discs self-regulate to the marginally stable *Q*≈ 1 condition over a large range of radii.

Once the disc has reached a quasi-steady state, the disc aspect ratio *H*/*R* also stabilizes to the value predicted by the self-regulation condition *Q*≈ 1:

*q*(bottom panel).

### 5.1 Saturation amplitude of the instability

Once the gravitational instability has been initiated, for the case where β≤β_{frag} (i.e. where β= 4) the amplitude of the perturbations required to balance the cooling rises to the point where non-linear effects dominate, leading to the fragmentation of the disc into bound objects. In the case where β > β_{frag}, however, the amplitude increases on the dynamical time-scale until the disc reaches dynamic thermal equilibrium, at which point the amplitude of the surface density fluctuations becomes constant, and the heating they provide balances the imposed cooling. This is observed in all our simulations where β≥ 5.

From the simulations we can now test the prediction for the saturation amplitude provided by equation (25) numerically. Fig. 3 shows images of the surface density of the disc for the two cases β= 5 and 10, respectively, where in both cases the mass ratio is *q*= 0.1. It can be seen that, while the overall disc structure remains essentially constant [as confirmed by a more detailed Fourier analysis (see below)], the spiral wave amplitude as characterized by the surface density contrasts appears to decrease with increasing β. Noting that the direction of rotation of the discs is anticlockwise, we find that throughout our simulations the waves excited are all trailing waves – they all point in opposition to the direction of rotation.

While the sph code allows us to conduct a global three-dimensional simulation of discs with relative ease, it does not permit the direct calculation of an intrinsically two-dimensional quantity, such as the surface density perturbation amplitude δΣ/Σ. Therefore, in order to calculate this quantity, we overlay a cylindrical grid on the disc, such that each cell contains approximately *N*_{neigh} particles, where *N*_{neigh}≈ 50 is the average number of neighbours within a smoothing kernel for our simulations. For each annulus of cells we can calculate an average surface density , and by comparing this to the value calculated for each cell within this annulus we can evaluate an annulus-averaged rms value for the perturbation amplitude . This is shown as a function of radius *R* and the cooling parameter β in Fig. 4.

From Fig. 4, it is clear that there is an increasing trend in with decreasing β and that furthermore, away from the disc boundaries, the saturation amplitude is approximately constant with radius. The low values for the perturbation amplitude at small radii (*R*≲ 5) are probably due to the increased number of particles per grid cell smoothing out the underlying variation. We can however characterize the strength of the perturbation by simply averaging over the self-regulated portion of the disc that we define as 5 ≤*R*≤ 25 (cf. Fig. 1). Fig. 5 shows the relation between the azimuthally and radially averaged amplitude, which we denote as , and the cooling parameter β. Each point represents a single simulation, while the curve shows our best fit to the data using the inverse square root dependence predicted by equation (25). From the simulations we therefore obtain the following empirical relationship for the case where *q*= 0.1:

In a similar manner we can also calculate the variation in with *q*, which is shown in Fig. 6. We see that the strength of the perturbation tends to increase with the mass ratio *q*, although it is clear that this dependence on *q* is rather less than linear.

### 5.2 Fourier analysis: azimuthal structure

From the simulations we have empirically found that the perturbation strength follows a β^{−1/2} relationship, as predicted by equation (25). However, this equation also shows a clear dependence on the wave modes excited within the disc through the action of the gravitational instability. To elucidate this relationship further, we have therefore conducted a Fourier analysis of the wave modes in the disc, a full description of which may be found in Appendix B. In this section, we therefore consider the effects of both the cooling (via β) and the disc to central object mass ratio *q* on the excitation of the azimuthal *m* wavenumbers – the next section will describe the excitation of the radial wavenumbers.

We find that in general, whatever the imposed cooling regime, for a given mass ratio *q*= 0.1, the distribution of the azimuthal wavenumbers determined by the gravitational instability remains approximately constant, with the dominant mode at around *m*≈ 5. The mode distributions at five radii throughout the disc for the cases where β= 4, 6, 8 and 10 are shown in Fig. 7. It is clear that the spectral distribution of the modes shows little variation with radius and with the imposed cooling, except that the amplitude of the modes decreases as β increases, as expected in view of the decreased perturbation amplitudes seen in Figs 3, 4 and 5. Additionally, Fig. 8 (top panel) shows the variation of the average wavenumber against radius for all the values of β that have been tested, and we note that, although some small variation is seen, it is uncorrelated with the imposed cooling.

The bottom panel of Fig. 8 shows the variation in the average azimuthal modes excited as the disc to central object mass ratio *q* varies, while the cooling is held constant with β= 5. It can be seen from the plot that variation of this parameter does have a marked effect on the power spectrum of the waves – the average mode number varies inversely with the mass ratio, from *m*_{av}≈ 15 where *q*= 0.05 to *m*_{av}≈ 10 where *q*= 0.125. This variation is also clearly seen in Fig. 9, where a large number of flocculent arms are present in the disc with *q*= 0.05, and fewer, rather more well-defined, spiral arms appear in the disc where *q*= 0.125 (cf. the similar result obtained in Lodato & Rice 2004). By comparison, the left-hand panel of Fig. 3 shows a disc with β= 5 and *q*= 0.1, and the pattern of spiral arms present is intermediate to those shown in Fig. 9.

### 5.3 Fourier analysis: radial structure

We now consider the radial wavenumbers *k* of the waves excited by the gravitational instability. In contrast to the azimuthal modes, it is clear from Figs 3 and 9 that there is significant variation in the radial wavenumber *k* with radius, and Fig. 9 suggests that there is an additional variation with the disc to the central object mass ratio.

In Fig. 10, we show the variation in the power spectrum of different radial wavenumbers for the cases where β= 4, 6, 8 and 10 and *q*= 0.1, at the same radii as the azimuthal wavenumbers shown in Fig. 7. As with the azimuthal modes, we see little overall change in the spectral distribution of the modes with varying β excepting that the amplitudes of the modes decrease as the cooling weakens. Conversely, however, these plots show a significant variation with radius, in that the peak wavenumber decreases with increasing radius, and thus the dominant wavelength similarly increases with radius.

Fig. 11, on the other hand, shows the power spectrum as a function of *kH*, where we normalize the wavenumber to the expected most unstable one, *H*^{−1} (see Section 2.1). We thus confirm the expectations from the linear WKB approach, as these plots show a clear peak at *kH*≈ 1. This, therefore, suggests that throughout the disc the excited waves are close to corotation. Fig. 12 (top panel) shows the average radial wavenumber as a function of radius for all the considered values of β, again confirming the trends already discussed and also further showing that, excepting the variation in amplitude discussed above, the structure excited by the simulation is essentially independent of the cooling imposed. It also shows that the simulation-to-simulation scatter is very small.

The bottom panel of Fig. 12 shows that, as with the azimuthal wavenumbers, there is clear variation in the average radial wavenumber *k*_{av} with the disc to central object mass ratio for a given β (in this case β= 5); increasing the mass ratio decreases the average wavenumber in approximate inverse proportion. The dashed line in the bottom panel of Fig. 12 plots a simple *R*^{−3/2} curve, indicating that the average wavenumber follows a power-law distribution with radius, such that *k*_{av}∼*R*^{−3/2}, which remains constant with varying mass ratio. This is easily understood by noting that since the sound speed *c*_{s} is approximately constant by construction, equation (4) indicates that *k*∼Σ∼*R*^{−3/2}.

### 5.4 Mach number of the spiral modes

Returning briefly to the dispersion relation given in equations (1) and (2), we note that this is only strictly valid for infinitesimally thin discs. As our simulations are fully three dimensional, we require a correction to the self-gravity term to account for this, and this is shown in equation (31) (Bertin 2000),

where we have also used the fact that our discs are approximately Keplerian, and thus κ≈Ω. The reduction factor of 1/(1 + |*k*|

*H*) arises from the vertical dilution of the gravitational potential due to the finite thickness

*H*of the disc (Vandervoort 1970; Bertin 2000; Binney & Tremaine 2008). Using this finite-thickness dispersion relation and our averaged values for

*k*and

*m*, it is possible to calculate a Doppler-shifted angular speed |Ω

_{p}−Ω|, noting that the sign of Ω

_{p}−Ω cannot be determined from equation (31). Since the average radial wavenumber is always very close to the most unstable one, and since the disc is almost exactly marginally stable, the resultant average pattern speed always turns out to be very close to corotation (as can be seen in Fig. 16). We quantify the deviation of the pattern speed from corotation later in Section 5.5.

We can further calculate the radial phase and Doppler-shifted phase Mach numbers, and these are shown in Fig. 13. The upper panel shows the wave phase Mach number (thick lines) and the Doppler-shifted phase Mach number (thin lines) as functions of radius for various values of β with *q*= 0.1. Similarly, the lower panel of Fig. 13 shows the variation of these Mach numbers with the mass ratio *q* for β= 5. We immediately see that both quantities are independent of the cooling rate as measured by β with very little scatter. Moreover, the Doppler-shifted phase Mach number is very close to unity. In a similar manner, this quantity remains unchanged with variations in the mass ratio, although the phase Mach number decreases with increasing *q*.

The above results essentially imply that the wave structure is determined by the requirement that the normal component of the flow into the shock is almost exactly sonic – a natural criterion for a quasi-steady system due to the dissipative nature of shocks. For waves with winding angle *i*, and radial Doppler-shifted phase speed , a sonic normal component of velocity into the shock implies that , leading to

*i*≈ 1, we should expect that , as indeed we find in Fig. 13. For completeness, Fig. 14 shows the winding angle

*i*as a function of radius for varying β (top panel) and mass ratio

*q*(bottom panel), using the definition tan

*i*=

*m*/

*kR*. In all cases,

*i*≲ 15°, so the waves are reasonably tightly wound throughout. Again there is no significant variation with cooling, but the structure becomes more open as the mass ratio increases, as expected from Figs 3 and 9.

We can also use equation (25) to estimate the amount of energy dissipated by the weak spiral shocks per dynamical time, as characterized by ε, which is shown in Fig. 15. We have seen through the constancy of the Doppler-shifted phase Mach number that the shock structure that forms in the disc is indeed self-similar, and thus the heating factor ε is also largely independent of the applied cooling, the mass ratio and the radial position. Note that the larger values for ε generated at low radii (*R*≲ 5) are probably due to the inaccuracies in calculating in this region rather than a breakdown in self-similarity.

### 5.5 On the locality of transport induced by self-gravity

In the previous section we noted that the (spectrally averaged) pattern speed of the waves Ω_{p} is always very close to the angular velocity of the flow Ω, thereby indicating that the waves are close to the corotation resonance as suggested earlier in the results of the radial mode decomposition. We can estimate more quantitatively how close to corotation the spiral waves lie by calculating the quantity ξ, as given in equation (22). This is shown in Fig. 16 as a function of radius for all the values of β and mass ratio simulated.

For the *q*= 0.1 case, we thus see that varying the cooling has no significant effect on the transport properties of the disc, with ξ≈ 0.1 throughout the radial range, albeit with some scatter. This means that in this configuration the disc is dominated by local transport processes, and is as such reasonably well described by the viscous α prescription of Shakura & Sunyaev (1973)– global effects, although not negligible, are smaller than local effects by an order of magnitude.

By varying the mass ratio, we see that the strength of non-local effects increases with *q*, rising to ξ≈ 15 per cent for the case where *q*= 0.125. This confirms the results of Lodato & Rice (2004), who similarly found that non-local effects (characterized by strong transient structures in the disc) become increasingly important as the disc mass ratio rises, although for the parameter range considered here the disc remains dominated by local effects. This non-local behaviour can be elucidated further by noting that from the definitions of ξ and (equations 22 and 18) we have

*kH*≈ 1, reduces to The non-locality of the transport is therefore directly linked to the openness of the structure that is induced in the disc through self-gravity. Since, as we have noted earlier, for larger disc masses the spiral structure tends to become more open and more dominated by low

*m*modes, we then see that more massive discs tend to become more subject to non-local effects.

## 6 DISCUSSION AND CONCLUSIONS

In this paper, we have undertaken three-dimensional global numerical simulations of gaseous, non-magnetized discs, evolving under the influence of a massive central object and their own self-gravity. We have modelled the gas as an ideal gas with γ= 5/3, together with a simple cooling prescription based on a local cooling time-scale. We have used these simulations to investigate the structure that forms once the discs have settled into a quasi-steady marginally stable state as a function of both the imposed cooling and the disc to central object mass ratio.

We have found that the amplitude of spiral arms induced in self-gravitating discs, as characterized by the rms surface density perturbations, can be described straightforwardly through the empirical relationship

(where β is the ratio between the local cooling and dynamical time-scales), with only a weak dependence on the disc to central object mass ratio. This is in fact closely linked to the result that the Doppler-shifted Mach number is very close to unity – by considering the entropy change Δ*S*across an adiabatic shock where the Mach number

*M*≈ 1, it can be shown that Δ

*S*∼ (

*M*

^{2}− 1)

^{2}. Thermal equilibrium in these discs is established between cooling, at a rate inversely proportional to β, and the irreversible conversion of mechanical energy into heat, at a rate proportional to the entropy jump Δ

*S*at the shock front – we therefore find that β∼ (

*M*

^{2}− 1)

^{−2}. Standard shock relations show that the density perturbation δρ/ρ∼ (

*M*

^{2}− 1), and hence simply from considering the properties of weak adiabatic shocks we can arrive at the relationship δρ/ρ∼β

^{−1/2}.

Additionally, we find that the heating factor ε– that fraction of the available wave energy that is liberated as heat back into the disc gas – remains essentially invariant at ≈20 per cent with both the imposed cooling regime and the mass ratio of the disc to the central object.

As expected, our simulations show that the dominant radial wavenumber is approximately equal to the reciprocal of the local scaleheight of the disc throughout the radial range, *k*≈π*G*Σ/*c*^{2}_{s}. We therefore find that the radial spacing of the arms is dependent only on the surface density and temperature profiles of the disc. Likewise, although further work is required to understand the relationship fully, the azimuthal disc structure is dependent on the disc to central object mass ratio, with more massive discs being characterized by more open structures than their lower mass counterparts for a given central object mass.

Our numerical results bear out the theoretical analysis of Balbus & Papaloizou (1999) and Gammie (2001), who suggest that discs in the *Q*≈ 1 marginally stable state may be modelled as predominantly local. Simulations of self-gravitating discs with radiative transfer by Boley et al. (2006) also found that close to corotation, angular momentum transport was well modelled by a local α-prescription even when global modes were present. Balbus & Papaloizou (1999) further predicted that non-local transport from an ‘anomalous flux’ proportional to Ω−Ω_{p} would become significant far from corotation, a result we have derived analytically using the WKB approximation for tightly wound waves. We have then used the WKB dispersion relation along with empirically determined information on the dominant wavenumbers to make an estimation of |Ω_{p}−Ω|. We find that, at least for low-mass discs, this is a small fraction of Ω (less than 15 per cent for discs with *q*≤ 0.125, regardless of the efficacy of the cooling). Our results on the magnitude of the non-local transport fraction ξ= |Ω_{p}−Ω|/Ω can furthermore be readily understood in terms of the empirical constancy of the Doppler-shifted radial phase Mach number, . We conclude that the importance of such non-local effects in gaseous self-gravitating discs is set by the self-adjustment of the pattern speed to ensure that the normal flow speed into the arms is sonic. We have then demonstrated that this condition implies that , where *i* is the opening angle of the spiral structure. Since the structure within the disc becomes more open as the disc to central object mass ratio increases, this implies that the importance of non-local transport also scales with *q*.

We also note that in collisionless systems such as stellar discs, this self-regulation process for the pattern speed breaks down as shocks cannot be formed. Hence, it is possible to excite global modes in such discs, and thus non-local transport of energy and angular momentum may be more significant dynamically. The results that we present here are therefore restricted to the case of predominantly collisional, gaseous discs. Our results provide a theoretical underpinning for the results of Lodato & Rice (2004, 2005 on how the importance of global transport depends on the disc to central object mass ratio in gaseous discs. In particular, we note that in cases (like those described here) where the disc mass is a small fraction of the central object mass (as could be the case for relatively evolved protostellar discs), the effects of self-gravity are expected to be well described as a pseudo-viscous process.

One of the most important applications of our study is that we can relate the amplitude of spiral modes in gaseous discs to the cooling regime. With the Atacama Large Millimetre Array (ALMA) coming online in the relatively near future, promising milli-arcsecond resolution in the millimetre/submillimetre range, it is possible that such observations of spiral structure in protoplanetary discs may become technically feasible.

*E*and angular momentum

*S*are found via

*E*=ℎω and

*S*=ℎ

*m*, where ω and

*m*are the angular frequency and spin quantum number, respectively.

We acknowledge the use of SPLASH (Price 2007) throughout this paper for the visualization of surface densities. We would also like to thank Jim Pringle for helpful discussions and a careful reading of the manuscript.

## REFERENCES

## Appendices

### APPENDIX A: RESOLUTION AND CONVERGENCE TESTS

In this appendix, we will briefly outline the tests that were undertaken to ensure the convergence of our results.

Three simulations were run, all with the cooling parameter β= 6 and mass ratio *q*= 0.1, using discs of 250 000, 500 000 and 1000 000 particles. These were otherwise identical to the simulations that were used for this paper, as described in detail in Section 4. The three values that are of most significance to our results are the rms surface density perturbation amplitude and the average radial and azimuthal wavenumbers *k*_{av} and *m*_{av}, respectively. Fig. A1 shows how varies with resolution, and although there is a considerable scatter, it is clear that there is no systematic variation with resolution. A similar result (with even less scatter) is also obtained when one conducts a Fourier analysis of the simulations – there is no systematic variation with resolution. We may therefore conclude that our simulations are converged, and that the resolution when using 500 000 particles is satisfactory for our purposes.

### APPENDIX B: FOURIER DECOMPOSITION METHODS

In this appendix, we detail how the Fourier mode analysis was conducted, using the sph particle positions as the input values. For simplicity, we begin by discussing how the radial *k*-mode amplitudes were computed, as this had practical implications for the azimuthal *m*-mode analysis.

#### B1 Radial mode analysis

To calculate the radial Fourier mode amplitudes within a disc presents certain problems, since even a cursory glance at Fig. 3 reveals that the radial wavenumber *k* varies significantly with radius. Also, unlike the azimuthal wavenumbers, the disc is neither uniform nor periodic in radius, and therefore the underlying Fourier distribution corresponding to the disc surface density profile also has to be taken into account. The following method addresses both of these problems while keeping the signal-to-noise ratio as high as possible.

The disc to be analysed is divided into numerous overlapping annuli of width Δ*R*, which varies with the central radius of the annulus, and into a number of sectors of fixed angular width Δφ. The Δ*R* values are chosen such that each annulus is of sufficient radial extent to resolve the greatest radial wavelength present at that radius, likewise each sector must be narrow enough to ensure that the wave crests are distinct and not smeared out across a wide range in *R*. In this manner, the smaller the winding angle θ= tan^{−1}|*m*/*kR*| of the waves, the wider the sectors can be for a given resolution.

Since the radial wavenumber profile depends on the disc to central object mass ratio, the radial extent Δ*R* of the annuli varied likewise, in order to capture all the relevant modes. The values used in our analyses are summarized in Table B1. Note that the widths of the annuli increase linearly across the disc, from the initial to the final widths quoted.

M_{disc}/M_{*} | Annuli | Initial width | Final width | Sectors |

0.050 | 25 | 2 | 8 | 60 |

0.075 | 25 | 2 | 8 | 60 |

0.100 | 25 | 2 | 10 | 60 |

0.125 | 25 | 2 | 10 | 60 |

M_{disc}/M_{*} | Annuli | Initial width | Final width | Sectors |

0.050 | 25 | 2 | 8 | 60 |

0.075 | 25 | 2 | 8 | 60 |

0.100 | 25 | 2 | 10 | 60 |

0.125 | 25 | 2 | 10 | 60 |

To calculate the underlying Fourier distribution due to the unperturbed surface density profile, the Fourier transform was taken over the whole of each annulus. This thereby smears out all the waves and takes the *average* distribution, and is evaluated according to the following relation:

*A*

_{k}is the

*k*-mode amplitude corresponding to the underlying disc distribution,

*N*

_{ann}is the number of particles per annulus,

*k*is the radial wavenumber and

*R*

_{i}are the radii of the individual particles.

The Fourier distribution of the waves overlaid on the disc is calculated by taking an equivalent transformation over each sector within the annulus, such that

where*A*

_{k,n}is the

*k*-mode amplitude of the waves and disc evaluated in the

*n*th sector and

*N*

_{sect}is the number of particles in that sector.

Finally, the Fourier distribution solely due to the waves in each sector is given by the difference between equations (B2) and (B1). Since each sector should be statistically similar to the others we may then average over all the sectors *N*_{sectors} (which in this case does not smear the wave component out, but reduces computational noise), to give the average radial Fourier mode amplitudes of the waves 〈*A*_{k}〉:

#### B2 Azimuthal mode analysis

For the azimuthal *m* wavenumbers, the analysis is more straightforward. The disc was initially divided into annuli of fixed width Δ*R* in such a manner that each of these annuli is narrow enough to ensure that the wave crests occupy only a small range in φ. In contrast to the radial modes, these annuli therefore need to become narrower with decreasing winding angle to maintain resolution. We found Δ*R*= 0.2 (in code units) to be sufficient for the purposes of this analysis. The azimuthal wavenumber amplitudes *A*_{m} within each annulus are then computed via

_{i}are the azimuthal angles of the individual particles,

*N*

_{ann}is the number of particles in each annulus and

*m*is the radial wavenumber of the wave, corresponding to the number of arms in the spiral.

However, to ensure that we have the azimuthal *m*-mode amplitudes specified at the same radii as the radial *k* modes, then an average value is taken of the *m*-mode amplitudes over all annuli where the central radius falls within that annulus in which the *k* modes are determined.

#### B3 Analysis checks

To ensure that the results of the Fourier analysis are accurate, we ran the following test case. A disc with an underlying surface density profile Σ∼*R*^{−3/2} and an analytically superimposed structure was created with five spiral arms, such that *m*= 5 throughout the entire disc and with the radial wavelength increasing linearly from λ_{min}= 2 to λ_{max}= 7.6. This gives a total of five full wavelengths across the face of the disc, which extends from *R*= 1 to 25. The surface density of the disc, clearly indicating the imposed structure, is shown in Fig. B1. The Fourier analysis was then conducted using the annulus widths quoted in Table B2, where, as described above, the width of the annuli increased linearly from the minimum to the maximum quoted value. The results from the azimuthal Fourier decomposition are shown in Fig. B2, and show that the azimuthal wavenumber is resolved extremely well. The fundamental frequency *m*= 5 is clearly dominant, with no other modes except higher harmonics present at any significant amplitude. Note that the results for the azimuthal modes show no sensitivity to the annuli used for the analysis, and are quoted for the first case in Table B2 where the annuli width correspond exactly to the radial wavelengths.

Annuli | Initial width | Final width | Sectors |

10 | 2.0 | 7.6 | 60 |

10 | 1.5 | 6.5 | 60 |

10 | 3.0 | 8.5 | 60 |

10 | 4.5 | 10.0 | 60 |

10 | 6.0 | 12.5 | 60 |

Annuli | Initial width | Final width | Sectors |

10 | 2.0 | 7.6 | 60 |

10 | 1.5 | 6.5 | 60 |

10 | 3.0 | 8.5 | 60 |

10 | 4.5 | 10.0 | 60 |

10 | 6.0 | 12.5 | 60 |

The results of the Fourier decomposition for the radial wavenumbers are shown in Fig. B3, which shows the actual distribution of wavenumbers (as calculated directly from the known distribution of wavelengths) and the distributions derived from analyses using the annuli given in Table B2. Note that the analysis using annuli that fit the wavelengths exactly correspondingly reproduces the exact result. Clearly, there is a scatter within the results for the radial wavenumbers, which arises from the fact that if the actual wavelength is not an integer divisor of the annulus over which the analysis is being conducted, more than one wavenumber appears to be excited. We note however that the scatter is never more than a factor of 1.5 above or below the true value, which we deem to be accurate enough for the purposes of this analysis.

#### B4 Resolution limits

Throughout this paper, we have used discs of 500 000 particles when undertaking the Fourier analysis. From the fundamental sph resolution limit of the local smoothing length *h* we can evaluate the maximum resolvable azimuthal and radial wavenumbers *m*_{max} and *k*_{max} as a function of radius throughout our discs, such that

*kH*≈ 1 and

*mH*/

*R*≈ 1, respectively, equations (B5) and (B6) show that the accuracy of the Fourier analysis is closely tied to the vertical resolution of the disc through

*H*/

*h*. Using the average smoothing length at each radius, these resolution limits are shown in Fig. B4. We have used data from the simulation where β= 10, as this gives the most conservative limits of all our experiments.

The vertical resolution of the disc as indicated by *H*/*h* is shown in Fig. B5 for simulations using 500 000 particles. We find that the disc thickness is covered by approximately two smoothing lengths throughout, and thus is adequately resolved. For the Fourier analysis we therefore see that the expected peak wavenumbers are resolved by a factor of approximately 4 π throughout the radial range. For the radial wavenumbers, we are primarily interested in *k* < 10, which is well resolved until at least *R*= 25, again adequate for the analyses we have undertaken. We conclude therefore that throughout the radial ranges of interest, the Fourier analyses we have presented are well resolved.