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 β=tcool/tdyn. 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 Mdisc/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.
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 tcool 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 haveBinney & 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)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 kuns (i.e. where the RHS of equation 2 is at a minimum) occurs at a radial wavenumber 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
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
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 TReynRφ is such thatLynden-Bell & Kalnajs (1972) as produced by this viscous torque is then simply given by (Pringle 1981; Frank, King & Raine 2002) Equation (11) therefore links the transport of angular momentum and the associated rate of work done by torques in the case of a local process, as historically modelled via the α-viscosity parameter.
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)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):1 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.
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:equation 11) and can therefore be represented using the α-formalism. The first term however, equal to the same angular momentum term times Ω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 thatSection 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 thatGammie (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 asand , 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=c2s/γ(γ− 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)., can therefore be estimated via
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: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 R0, and thus, with G= 1, the unit time is the dynamical time tdyn=Ω−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 Mdisc. Although the bulk of the simulations have been conducted with a disc to central object mass ratio q=M*/Mdisc 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 tvisc≫tdyn=Ω−1, this profile remains roughly unchanged throughout the simulations. The initial temperature profile is c2s∼R−1/2 and is such that the minimum value of the Toomre parameter Qmin= 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 Rin= 0.25 to Rout= 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 ttherm for each simulation as the cooling time evaluated at the initial outer edge of disc, taken to be at R= 25– thus ttherm=tcool(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=Mdisc/M*||No. of particles||Duration|
|4||0.10||500 000||4.0 ttherm|
|5||0.10||500 000||10.0 ttherm|
|6||0.10||500 000||10.0 ttherm|
|7||0.10||500 000||10.0 ttherm|
|8||0.10||500 000||10.0 ttherm|
|9||0.10||500 000||10.0 ttherm|
|10||0.10||500 000||10.0 ttherm|
|5||0.050||500 000||10.0 ttherm|
|5||0.075||500 000||10.0 ttherm|
|5||0.100||500 000||10.0 ttherm|
|5||0.125||500 000||10.0 ttherm|
|β||q=Mdisc/M*||No. of particles||Duration|
|4||0.10||500 000||4.0 ttherm|
|5||0.10||500 000||10.0 ttherm|
|6||0.10||500 000||10.0 ttherm|
|7||0.10||500 000||10.0 ttherm|
|8||0.10||500 000||10.0 ttherm|
|9||0.10||500 000||10.0 ttherm|
|10||0.10||500 000||10.0 ttherm|
|5||0.050||500 000||10.0 ttherm|
|5||0.075||500 000||10.0 ttherm|
|5||0.100||500 000||10.0 ttherm|
|5||0.125||500 000||10.0 ttherm|
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.Fig. 2 for different values of β (top panel) and 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 Nneigh particles, where Nneigh≈ 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 mav≈ 15 where q= 0.05 to mav≈ 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 kav 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 kav∼R−3/2, which remains constant with varying mass ratio. This is easily understood by noting that since the sound speed cs 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),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, 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
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
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Σ/c2s. 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.
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.
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 kav and mav, 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.
|Mdisc/M*||Annuli||Initial width||Final width||Sectors|
|Mdisc/M*||Annuli||Initial width||Final width||Sectors|
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:
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 Nsectors (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 〈Ak〉:
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 Am within each annulus are then computed via
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|
|Annuli||Initial width||Final width||Sectors|
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 mmax and kmax as a function of radius throughout our discs, such thatequations (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.