Introducing two improved methods for approximating radiative cooling in hydrodynamical simulations of accretion discs

The evolution of many astrophysical systems depends strongly on the balance between heating and cooling, in particular star formation in giant molecular clouds and the evolution of young protostellar systems. Protostellar discs are susceptible to the gravitational instability, which can play a key role in their evolution and in planet formation. The strength of the instability depends on the rate at which the system loses thermal energy. To study the evolution of these systems, we require radiative cooling approximations because full radiative transfer is generally too expensive to be coupled to hydrodynamical models. Here we present two new approximate methods for computing radiative cooling that make use of the polytropic cooling approximation. This approach invokes the assumption that each parcel of gas is located within a spherical pseudo-cloud which can then be used to approximate the optical depth. The first method combines the methods introduced by Stamatellos et al. and Lombardi et al. to overcome the limitations of each method at low and high optical depths respectively. The second, the"Modified Lombardi"method, is specifically tailored for self-gravitating discs. This modifies the scale height estimate from the method of Lombardi et al. using the analytical scale height for a self-gravitating disc. We show that the Modified Lombardi method provides an excellent approximation for the column density in a fragmenting disc, a regime in which the existing methods fail to recover the clumps and spiral structures. We therefore recommend this improved radiative cooling method for more realistic simulations of self-gravitating discs.


INTRODUCTION
Radiative transfer is an important aspect of many hydrodynamical simulations.Accurate treatment of heating and cooling is particularly important for systems undergoing gravitational instability or fragmentation, such as star-forming clouds and young, self-gravitating protoplanetary discs, in which the strength of the gravitational instability is a function of the thermal balance between shock heating and radiative cooling (Gammie 2001).This determines whether the disc develops spiral arms, fragments or maintains a smooth, axisymmetric structure.
Early work using isothermal simulations suggested that a selfgravitating phase should be very short-lived, with the disc either fragmenting quickly or spiral instabilities acting to stabilize the disc's density profile (Laughlin & Bodenheimer 1994).However, simulations that implemented more detailed thermodynamics demonstrated that this was mostly a consequence of the isothermal assumption and that a self-gravitating phase may persist for many orbits (Lodato & ★ E-mail: alison.young@ed.ac.uk (AKY) Rice 2004).Now that we expect planet formation to begin very early when discs are likely to be self-gravitating (Nixon et al. 2018), we need robust models of self-gravitating discs to address key questions including how planetesimals grow in young discs.However, directly implementing radiative transfer in 3-D simulations is by no means trivial.A full treatment would require wavelength-dependent and potentially anisotropic absorption, scattering, and emission of photons; tracing of the photons' paths through the system; and variable opacities and emissivities for both dust and gas components.Additionally, explicit radiative transfer timesteps are prohibitively short relative to hydrodynamic timescales, particularly in optically thick locations, greatly increasing computation time (Castor 2004;Brandenburg & Das 2020).These difficulties have prompted the development of a number of approximate radiative transfer methods, each depending on some simplifying assumptions and approximations.Such methods are naturally limited in accuracy and scope compared to a full radiative transfer scheme, but allow for the simultaneous simulation of both hydrodynamic and radiative processes.
One such method which has been used to study the effects of cooling in protoplanetary discs is the so-called -cooling model, de-veloped by Gammie (2001).This approach parameterises the local cooling time   in terms of a  parameter and the local Keplerian angular velocity  K , such that   = / K .The simplicity of this parameterisation and ease of implementation in hydrodynamic codes has made it a popular choice for studying the relationship between cooling and fragmentation in discs (e.g.Rice et al. 2003;Cossins et al. 2009;Meru & Bate 2011;Boss 2017), but it does not paint a complete picture.There is no compelling argument for the  parameter remaining constant throughout the extent of the disc (Mercer et al. 2018), and more realistic cooling methods yield significantly different results (Johnson & Gammie 2003;Vorobyov et al. 2020).
There are two main classes of more physically-motivated cooling methods for protoplanetary disc simulations: those that impose an external radiation field and calculate radiation transport through the disc; and those that make use of local parameters to estimate the thermal balance at a given location in the disc.Examples of the first approach include codes that make use of the flux-limited diffusion approximation (Levermore & Pomraning 1981), such as Ayliffe & Bate (2011) and Meru & Bate (2010) with fixed prescribed temperatures at the edge of the model; or Monte-Carlo methods to transport photon "packets" throughout the disc, such as the coupled mcfost + phantom model (Pinte et al. 2006;Price et al. 2018;Nealon et al. 2020).For the coupled radiative transfer and hydrodynamics approach, a radiative equilibrium calculation is performed every few hydrodynamical timesteps (this is chosen to be as infrequent as possible to reduce the overheads), after which the temperatures are updated.Although the radiative transfer is more accurate than under the fluxlimited diffusion approximation and can, for example, reproduce the effects of shadowing, the coupled method is not time-dependent due to the radiative equilibrium assumption.The outcome of this is that temperatures are likely to be under-or over-estimated in different regions.The second approach -estimating the cooling rate based on local properties -is the focus of this paper.Existing methods are sufficient for modelling symmetrical systems, i.e. spherical cloud collapse, but are less accurate for modelling discs since they overestimate the optical depth in the disc.This then changes the cooling time, which in turn affects the evolution of self-gravitating discs.The use of Monte Carlo radiative transfer would allow for irregular structures but is very costly and currently imposes radiative equilibrium.As discussed above, it is crucial to have a robust cooling model for studying gravito-turbulent discs, therefore a more accurate approximation for the optical depth within discs and therefore the cooling rate is needed, without adding a large computational expense.
In section 2 we first describe the methods that are currently used to approximate radiative cooling by assuming that each parcel of gas, or SPH particle, is embedded within a spherical gas cloud ("polytropic cooling") and using this to estimate the optical depth.The limitations of these methods are detailed in 2.5 and we further motivate the development of our new techniques.In section 3.1 we introduce the first of our improved methods for approximating radiative cooling in discs, which combines the two.We present our second method in section 3.2 which provides an alternative means of regularising the Lombardi et al. (2015) method in the mid-plane by making use of disc-specific geometry.

EXISTING POLYTROPIC COOLING METHODS
Much of what we present here will be based on simulations using Smoothed Particle Hydrodyamics (SPH) (e.g., Monaghan 1992), but it is also applicable to other hydrodynamic formalisms.In SPH, each particle represents a parcel of gas and to estimate its radiative cool-ing rate we need to know its optical depth.Physically, this involves integrating the opacity and density along the path from the particle to the edge of the particle distribution (e.g., in the case of a disc, to the disc surface).The cooling approximations described below instead estimate a mean column density from the particle to the surface and combine this with mean opacity estimates.Two different approaches for estimating this mean column density are described, firstly using the gravitational potential (Stamatellos et al. 2007) and secondly using the pressure gradient (Lombardi et al. 2015).

The Stamatellos Method:
The gravitational potential approach Stamatellos et al. (2007) introduced a method for approximating the optical depth of an SPH particle from its gravitational potential.The gravitational potential energy is already calculated in SPH simulations that include the gas self-gravity, hence minimising additional computational expense.This approach treats each parcel of gas as if it were embedded within a spherically symmetric 'pseudo-cloud' with polytropic density and temperature profiles.The scale-length of the polytrope is determined by the gravitational potential at the location of the SPH particle.Because we don't know where this gas parcel lies within the pseudo-cloud, we must calculate a mass-weighted average of the pseudo-column density over all possible positions within it.
The mean column density for particle  is then given by where   is a constant that depends on the polytropic index (which is related to ; see Stamatellos et al. 2007 for the full derivation), and   and   are the gravitational potential and density at the position of the particle respectively.Following Stamatellos et al. (2007) and Forgan et al. (2009) we use  2 = 0.368.We now require a pseudo-mean mass opacity κ (, ), to estimate the mass weighted optical depth.The pseudo-mean opacity for particle  is calculated by averaging over locations in the pseudo-cloud in the same manner as for the mean column density.This allows us to estimate the optical depth using τ = Σ κ . (2)

The Lombardi Method: the pressure gradient approach
A slightly different approach was introduced by Lombardi et al. (2015) in which the pressure scale-height is used to estimate the optical depth rather than the gravitational potential.This method still assumes that each particle is embedded within a polytropic pseudocloud but by using the pressure gradient at the location of the particle, it performs better for non-spherical gas distributions.
The pressure scale-height,  p, , of particle  is given by: where   is the pressure at the location of particle , and ∇  is the pressure gradient.The hydrodynamic acceleration excluding the contributions from artificial (and real) viscosity and gravity is In smoothed particle hydrodynamics, the hydrodynamic acceleration  h, is already computed during the force calculation.For particle , this is calculated via the following sum over  neighbour particles: in which   is the particle mass,   (  ) is the gas density evaluated at particle  ( ), and the gradient of the smoothing kernel is ∇.Ω is related to the gradient of the smoothing length ℎ We can then obtain the pressure scale-height via Calculating the mass-weighted pseudo-mean column density in the same manner as Stamatellos et al. (2007), Lombardi et al. (2015) obtain with  ′ = 1.014.The optical depth can then be estimated using Equation 2.

Equation of state
We employ an identical approach to Stamatellos et al. (2007) in which values of the pseudo-mean opacity κ , Planck mean opacity   , mean molecular weight   , and specific internal energy   are pre-computed for the required density and temperature ranges and stored in a look-up table to reduce computational expense.These were calculated following Bell & Lin (1994).The Rosseland-mean and Planck-mean opacities are assumed to be the same and to depend only on density and temperature and are given by: where the constants  0 ,  and  are chosen for each main physical process contributing to the total opacity in that temperature and density range (see Bell & Lin 1994;Stamatellos et al. 2007 for further details).The tabulated values are then interpolated during the live hydrodynamical computation and the ideal gas approximation is used to obtain the gas pressure.

Updating the energy
Both methods give us a mass-weighted column density and tabulated values of the pseudo-mean opacity.These quantities are now used to calculate the radiative heating rate of gas parcel represented by particle : where  0 is a minimum temperature determined by background radiation and  −1  (  ,   ) is the Planck mean opacity, which is tabulated along with the mass-weighted opacity κ (  ,   ).External heating from nearby stars can be included by calculating  0 for each gas particle to prevent the gas cooling radiatively below this minimum temperature. 1The equilibrium temperature of particle  is found by assuming that the radiative cooling and the compressive and viscous heating are in balance: Substituting the equilibrium temperature  eq, for   in Equation 10gives: In the case that  4 eq,i ≤  4 0 (  ), we set  4 eq,i =  4 0 (  ).The equilibrium energy is then  eq, = ( eq, ,   ).The energy is updated by considering the approach towards the equilibrium temperature with thermal timescale For implementation within a typical SPH code, the cooling rate is then given by The energy equation is then updated via

Limitations of these methods
In both cases, given the estimate is based only on local parameters already used in hydrodynamic calculations in SPH, there is no need for communication of information with other fluid elements, and there is therefore a minor computational cost when compared with methods that utilise a radiation field or photon packets.However, any local approximation of optical depth will necessarily require assumptions about the system's global geometry.In the case where one wishes to study systems with multiple geometries co-occurring -for example, discs with spherical clumps and spiral waves, or with low-density envelopes or streamer tails -accurately estimating the cooling rate in all regions of the disc with just one method can be problematic.
The Stamatellos et al. (2007) method significantly over-estimates the optical depth in lower density regions of the disc but performs better towards the disc mid-plane and in spherical clumps.The Lombardi et al. (2015) method is more accurate overall in disc simulations (Mercer et al. 2018), but breaks down in regions where the pressure gradient approaches zero, such as in the disc mid-plane or the center of clumps.By combining the two methods, one can take advantage of the strengths of each approach to achieve more accurate cooling The column density from  to  out for an  = 3/2 polytrope (blue line, "Exact") compared with the estimate from Lombardi et al. (2015).Also shown are the estimate generated by the combined method (see subsection 3.1) and a modification of the Lombardi method for discs (see subsection 3.2).
throughout the whole of the disc.It is also possible to improve on the Lombardi et al. (2015) method in the regions where it fails by using another means of estimating the optical depth that makes use of a disc-specific density distribution assumption, as there is a simple relationship between mid-plane density, scale-height and column density in discs.

METHOD
In this section we detail two new methods for estimating the optical depth of a particle.First, in 3.1 we combine the estimates from the methods of Stamatellos et al. (2007) and Lombardi et al. (2015), and second in 3.2 we modify the method of Lombardi et al. (2015) specifically tailored for modelling self-gravitating discs.

The Combined Method: combining estimates of optical depth
While Lombardi et al. (2015) and Mercer et al. (2018) showed that the pressure scale-height and density can be used to provide a good estimate of the column density (Equation 8), this estimate breaks down in regions where the pressure gradient is close to zero, such as at the centre of a fragment or close to a disc mid-plane.Since a large part of the mass can be found under such conditions, the average cooling rate can be affected and we are motivated to find a way to improve the estimate under these conditions.
Here we explore replacing the Lombardi method with a surface density estimated by the Stamatellos method in regions where ∇ → 0 as the Stamatellos method remains well defined in such regions.
First, we define a scale-height based on the Stamatellos method via which we combine with the Lombardi estimate  P in inverse quadrature, i.e.

Surface density, Σ(Z)
Modified Lombardi: This selects the smaller of the two estimates, chosen because the Stamatellos method is known to overestimate the surface density in disc geometries and therefore generally we will only want to use this estimate when the Lombardi method is also failing, producing overestimates itself.
As a first test of this Combined method we consider a polytrope as an approximation to the structure of a self-gravitating fragment, and show Σ() for an  = 3/2 polytrope and also the different column density estimates in Figure 1.As previously stated, the Lombardi method works well away from the centre of the polytrope but overestimates the column density once the density and pressure gradients tend to zero close to the centre of the polytrope.The Combined method offers an improvement over the Lombardi method alone near the core, but it does lead to a moderate underestimate of the column.This follows because the Stamatellos method uses an estimate of the column density that is averaged over the entire polytrope and therefore must underestimate the column at the centre, where the column is highest.
Next, we consider how the Combined method works in a disc geometry.To this end, we consider a thin, self-gravitating disc.Highresolution simulations show that self-gravitating discs have a temperature structure that is typically close to vertically isothermal (Shi & Chiang 2014;Booth & Clarke 2019), so we choose a verticallyisothermal structure.Then vertical hydrostatic equilibrium implies where we have assumed that the disc's self-gravity can be calculated using a slab approximation, as in Wilkins & Clarke (2012) for example.Here  * =  s / K , where  s is the isothermal sound speed,  K is the Keplerian angular frequency, and Σ = ∫  0 ()/(0)d.For the last equality we have introduced  3D = 2 K /4 (0), which is the 3D definition of the Toomre Q parameter (e.g.Mamatsashvili & Rice 2010;Lin & Kratter 2016).
To compute the potential we need to extend this 1D model to a global disc because the potential of an infinitely extended selfgravitating slab is undefined.We choose a disc that has a constant  3D with radius so the vertical structure is given by the same solution to Equation 18 everywhere.The definition of  3D then enforces the mid-plane density scaling as  0 ∝  −3 .Specifying the scale-height,  * , fixes the surface density and the disc model is complete.Here we choose  * / = 0.1.We then compute the potential by direct integration of Green's function for Poisson's equation in 3D.We assume that the disc extends from 5 au to 100 au for consistency with the discs used in the SPH simulations.We note that while the final potential is not perfectly consistent with the disc structure, the difference is small when the disc is thin (e.g.Bertin & Lodato 1999;Wilkins & Clarke 2012) In Figure 2 we compare the surface density for the analytic disc model to the Lombardi estimate using the exact  P from Equation 18.The estimate works well at large  but, as expected, it overestimates the surface density at  = 0. To regularize the Lombardi method we evaluate the potential assuming a point 50 au from the star.Here, we see that the Combined method improves on the Lombardi estimate, as expected.In this case, the Combined method moderately overestimates the column at the mid-plane, which is expected since the Stamatellos method has been shown previously to overestimate the column in disc geometries (Wilkins & Clarke 2012;Lombardi et al. 2015;Mercer et al. 2018).
These two simple tests demonstrate that combining the Lombardi and Stamatellos methods should offer a significant improvement over using either method alone.In section 4 we explore this for more realistic scenarios that arise in SPH simulations.

The Modified Lombardi method: An improved Lombardi method for disc simulations
The Stamatellos method is not the only viable option for regularizing the Lombardi method.Since the Lombardi method only fails in regions where ∇ ≈ 0, such as the mid-plane of discs or the centre of polytropes, we only need to replace the estimate at such locations.
Here we explore an alternative that makes use of the fact that there is a simple relationship between the mid-plane density, scale-height and column density in a disc geometry.
For a non-self-gravitating disc, the (half) column density and mid-plane density are related via Σ/2 =  0 √︁ /2 * , and so  0 = √︁ /2 * could be used in the place of  S .Away from the mid-plane  P < √︁ /2 * , so the Lombardi method would take over.However, using √︁ /2 * overestimates the column density for self-gravitating discs or in fragments because self-gravity further compresses the gas, reducing the scale-height.Fortunately, it is possible to estimate  0 for a disc with any  3D .For  3D ≫ 1,  0 = √︁ /2 * suffices.In the regime  3D ≪ 1 the second term in Equation 18dominates and the equation can be made dimensionless by substituting  with  ′  where  =  * √︁  3D .We therefore expect  0 / * to scale with √︁  3D when  3D is small, and indeed find the fitting formula, describes Σ(0)/(0) well, where  = √︁ /2 (Figure 3) 2 .An advantage of Equation 19 is that it can be evaluated using quantities local to the disc mid-plane only (i.e., ,  s , and  K ).Further, while Equation 19 is only strictly meaningful at the mid-plane it can be evaluated anywhere as it only depends on local quantities.Since  0 increases away from the mid-plane (because  decreases) and  P decreases, replacing  S with Equation 19 in Equation 17 results in a method that produces similar estimates to the Lombardi method away from the mid-plane, which is generally quite accurate.Indeed, Figure 2 shows that this 'Modified Lombardi' method outperforms the normal Lombardi and Combined Lombardi + Stamatellos methods in discs.For the method to be robust, however, we need it to work well when the geometry deviates significantly from a disc.
Next, we apply this Modified Lombardi method to a polytrope as a model for a fragment in a self-gravitating disc.Under such conditions, we expect  3D ≪ 1, such that  0 ≈ 2 s / √︁ 4  0 .We apply this estimate of  0 to the polytrope in Figure 1, which shows that Equation 19 provides an estimate of Σ() that is accurate to within a factor of 2, in a spherical geometry.Since the size of a polytrope scales with  s √︁ ( + 1)/4 (0) (e.g.Clarke & Carswell 2014), the accuracy of this estimate is independent of  3D as long as  3D ≪ 1.Since the modified Lombardi method works well for disc geometries in general and for polytropic clumps when  3D ≪ 1, we expect it to be widely applicable to the structures present in self-gravitating discs and it also works well for spherical geometries.Finally, it should be noted that in deriving Equation 19we have assumed hydrostatic equilibrium.As a result, the Modified Lombardi method may not offer a significant improvement away from hydrostatic equilibrium.

TESTS
We now present tests using the SPH code phantom (Price et al. 2018).For comparisons with previous work, we first simulate the collapse of a uniform density sphere and then consider both low-and high-mass discs.For the collapse of the uniform density sphere, we only consider the Combined method, but for the low-and high-mass protoplanetary discs we consider both the Combined method and the Modified Lombardi method.The key factor that determines how accurately the cooling rate is estimated and the difference between the methods is how well the column density at each position is estimated.The column density estimates are therefore the focus of this section.

Collapse of a uniform density sphere
The gravitational collapse of a uniform density sphere represents the collapse of a protostellar cloud core and forms a useful benchmark (see Masunaga & Inutsuka 2000 for the first radiation hydrodynamic simulation of such a collapse).For this test, we use 2 × 10 6 SPH particles to simulate a cloud core of mass 1.5 M ⊙ , with an initial radius of 10 4 au, and an initially uniform temperature of 5 K.
Fig. 4 shows the evolution of the maximum density and temperature of the cloud core using the Stamatellos method, the Lombardi method, and the Combined method.For the Stamatellos and Combined methods, this is simply the density and temperature of the SPH particle with the highest density.For the Lombardi method there is considerable scatter so we average the 200 SPH particles with the highest density.The evolution of models with the Stamatellos, Lombardi and Combined methods are very similar.Moreover, Figure 4 shows that these methods produce results that are close to what is expected (black dash-dot line from Masunaga & Inutsuka 2000).As mentioned in Stamatellos et al. (2007), the differences with respect to Masunaga & Inutsuka (2000) are likely due to the different opacities implemented.
The small differences in the evolution between the models are due to the methods used to estimate the column density.Therefore, we now examine the column density estimates in each case.For this analysis we used a single simulation (using the Combined cooling method) and determined the column density estimate for all three cooling methods at different times in this simulation.
Figure 5 shows the resulting column density estimates as a function of radius for two snapshots during the collapse.We compare these estimates with the "true" column density (black dashed lines), calculated by integrating the density outwards from the specified cloud radius to the surface.Near the cloud core centre, the Lombardi method initially over-estimates the column density and Fig. 5 displays a lot of scatter.This is because the pressure gradient will initially be close to zero near the core centre and so the Lombardi method will not produce reliable estimates of the column density.It is because of this scatter that it is necessary to average the 200 particles with the highest density when producing the results for the Lombardi shown in Fig. 4.
On the other hand, near the cloud centre, the Stamatellos method, the Combined method (which tends to the Stamatellos method in these regions) and the Modified Lombardi method initially underestimate the column density, but avoid the large scatter seen in the Lombardi method.However, at larger radii, the Stamatellos method clearly over-estimates the column density, while the Lombardi and Modified Lombardi methods produce values much closer to the real column density.In these regions, the Combined method tends to the Lombardi method and produces a reliable estimate of the column density.
The right-hand panel of Figure 5 shows that as the cloud collapse progresses, the Combined method tends to the Lombardi method throughout most of the cloud core and both return reasonable estimates of the column density.The Modified Lombardi method also provides a reasonable estimate of the column density, despite being designed for a disc geometry.

Low-mass disc
We first consider a low-mass disc, similar to that presented in Mercer et al. (2018).The disc has an initial surface density profile of Σ ∝  −1 , extending from  = 5 au to  = 100 au, and an initial sound speed profile of  s ∝  −0.25 .The total disc mass is  disc = 0.01 M ⊙ , and the aspect ratio at  = 5 au is initially / =  s /Ω = 0.05.The disc gas is represented by 2 × 10 6 SPH particles and the central star has a mass of  * = 0.8 M ⊙ and is represented by a sink particle.The artificial viscosity parameters are  SPH = 1.0 and  SPH = 2.0.
We run the simulation for  = 1200 years to allow the disc to settle and specify a minimum disc temperature of  0 = 25 K.The column density snapshot of the disc at  = 1200 years (Fig. 6) shows  the disc has a smooth structure.We then consider a thin azimuthal ring within the disc at  = 75 au.Fig. 7 shows the vertical density structure at  = 75 au with a by-eye fit to the particle plot.Fitting with the expected profile of  =   exp(− 2 /2 2 ) demonstrates that the scale-height at  = 75 au is about  = 6.7 au.
We now compare this value obtained directly from the simulation to estimates obtained using each of the Stamatellos, Lombardi, Combined, and Modified Lombardi methods.Estimates of the scaleheight, , for a randomly selected sample of SPH particles are presented in Fig. 8, along with the value obtained directly from the simulation ( = 6.7 AU, dotted line, c.f., Fig. 7).Near the midplane, the Stamatellos method (open yellow circles) produces a reasonable estimate for the scale-height, but at higher altitudes is clearly far too large.The Lombardi method is, again, unreliable near the midplane, since this is a region where the pressure gradient is tending towards zero, but appears to produce more reasonable estimates at higher altitudes.Combining the two methods (Fig. 8, blue stars) then produces estimates that tend to the Stamatellos method near the midplane, and to the Lombardi method at higher altitudes.The Modified Lombardi method is an almost exact match to the expected scale-height in the midplane, and then follows the Lombardi method at higher altitudes.
Fig. 9 shows the column density estimates found using equations 1 and 8 at  = 75 AU as a function of the height above the mid- The column density estimates from the Stamatellos, Lombardi, and Combined methods are determined by multiplying the scaleheight, , in Figure 8 by the local gas density.It should be noted that the column density estimates from the four methods are one-sided column densities (i.e., from the  to the nearest disc surface) and, hence, at the midplane are about a factor of two smaller than the full column/surface density (dotted line in Fig. 9).
Comparing the column densities from the SPH simulation itself with the various estimates, shows that the Lombardi method is again unreliable near the midplane ( ∼ 0) but produces reliable estimates at higher , while the Stamatellos method produces reasonable estimates near the midplane, but produces column density estimates that are much too high at higher .Both the Combined method and the Modified Lombardi method produce reasonable estimates of the column density throughout the disc column.
Finally, Fig. 10 shows the midplane column density estimates Σ() and half the disc surface density determined directly from the SPH simulation, which we take as a reasonable estimate of the column density from the disc midplane to the surface.For the Lombardi, Stamatellos, Combined and Modified Lombardi method estimates, we plot a random sample of column density estimates that are within 0.05 AU of the midplane ( = 0).Again, the Lombardi method shows a lot of scatter, while the combined method removes most of this scatter and produces reasonable estimates for the column density at almost all radii.However, it is clear that the Stamatellos, Lombardi, and Combined methods all slightly over-estimate the midplane column densities, with the discrepancy becoming large at small disc radii.The Modified Lombardi method matches the expected column density at all radii, illustrating how this is a particularly reliable method for disc-like geometries.

High-mass disc
We now consider two different high-mass discs, both evolved using the Combined method.Both discs initially extend from  = 5 AU to  = 100 AU, have surface density profiles of Σ ∝  −1 , sound speed profiles of  s ∝  −1/4 , initial scale-heights at  = 5 AU of  = 0.275 AU (i.e., / = 0.055), and have central star masses of  * = 0.8 M ⊙ .The disc masses are  disc = 0.2 M ⊙ and  disc = 0.25 M ⊙ , and each disc is represented by 2 × 10 6 SPH particles.The artificial viscosity parameters are  SPH = 1.0 and  SPH = 2.0 and a floor temperature of  0 = 10 K was imposed.Fig. 11 shows that the  disc = 0.2 M ⊙ settles into a quasi-steady state with spirals, while the  disc = 0.25 M ⊙ disc becomes violently unstable and undergoes fragmentation.
We now consider how well the different methods for estimating the  Since the midplane is a region where the pressure gradient will tend to zero, the Lombardi method again shows a lot of scatter.The Combined method (blue stars) removes most of this scatter and seems to do reasonably well at capturing the spiral structure in the disc, for example it captures the trough at  ∼ 45 AU and the peak at  ∼ 55 AU.However, as in Fig. 10, the Stamatellos, Lombardi and Combined methods tend to over-estimate the midplane column density, generally by a factor of about 2. The Modified Lombardi method not only captures the midplane density structures, but also mostly matches the expected column density.
Turning now to fragmenting discs, we examine how well the methods can recover the structure of clumps within a disc.The clump considered here formed in the  disc = 0.25 M ⊙ disc (Fig. 11, right   panel) and is located at  ∼ 45 AU,  ∼ 45 AU in the snapshot.The column density as a function of the radius from the centre of this clump is presented in Fig. 13.Since we don't know the exact boundary of the clump, we calculate the "true " column density (black dashed line) by integrating the gas density out to a clump radius of 10 AU, which could slightly over-estimate the column density as it may include some of the surrounding gas disc.The Lombardi method over-estimates the column density near the clump centre, while at larger radii, the Stamatellos methods produces column density estimates that are too high.Combining the two methods seems to produce estimates that are reasonably close to that expected.In this case, the Modified Lombardi method produces estimates very similar to those from the Combined method.
Finally, Figure 14 shows a radial cut through the  disc = 0.25 M ⊙ disc (Fig. 11 right panel) that includes the clump at  ∼ 65 AU, the column density of which was shown in Fig. 13.All the methods capture the clump but, again, the Lombardi method over-estimates the column density near the clump centre, while the Stamatellos method slightly over-estimates the column density in the wings of the clump.The Combined and Modified Lombardi methods both closely recover the peak of the clump and the clump profile.

SUMMARY
We have described and tested two new methods for approximating radiative cooling in hydrodynamical models in discs.The Stamatellos method, estimating the optical depth via the gravitational potential, has been shown to overestimate the column density, and therefore underestimate the cooling rate in discs (Mercer et al. 2018).The Lombardi method, in which the column density derives from the local quantity of the pressure gradient, provides a more accurate estimate of the radiative cooling rate (Lombardi et al. 2015) but we show that the column density estimates become unreliable close to the disc midplane.These issues are resolved by combining the two methods, to achieve improved estimates in regions of both high and low optical depth.This method is particularly suited to spherical geometries.For disc geometries, we present the Modified Lombardi method, in which the Lombardi estimate is adjusted close to the disc midplane using an estimate of the scale height in a self-gravitating disc.Whereas Mercer et al. (2018) show that the Lombardi and Stamatellos methods both fail to provide accurate estimates of the optical depth in a disc with clumps, the Modified Lombardi method closely recovers the column density throughout a clumpy disc and traces the peaks and troughs associated with clumps and spirals well.Of course, since the optical depth is estimated from local disc properties, more complex geometries will not be modelled well.Full radiative transfer is required to include the effects of shadowing, for example.
It is straightforward to include the effects of a central star or stars, external stellar companions and/or a variable interstellar radiation field since the background temperature can incorporate arbitrary external heating.Radiative transfer can be incorporated via the flux limited diffusion approximation as formulated by Forgan et al. (2009), since the implementation of the new methods is identical to that of Stamatellos et al. (2007).A future development will be to consider the heating and gravitational influence of multiple sink particles.
The Modified Lombardi method is ideal for modelling selfgravitating discs and presents exciting prospects for studying young protoplanetary discs with a physical treatment of cooling.

Figure 2 .
Figure2.Top: The column density from  to infinity for the analytic selfgravitating disc model.Also shown are different estimates based on the Lombardi method and Combined method, as well as a Modified Lombardi method applicable to discs (see subsection 3.2).Bottom: The ratio of the estimated column density to the exact column density.

Figure 3 .
Figure 3.A comparison between the half thickness of a self-gravitating disc for different  3D and the fit formula given in Equation 19.

Figure 4 .
Figure 4.The evolution of the maximum density and temperature during the collapse of a uniform density sphere with three radiative cooling methods: Stamatellos method (yellow line), Lombardi method (red dotted line), and the Combined method (blue dashed line).The results from the simulation of Masunaga & Inutsuka (2000) are plotted (black dashed-dot line) for comparison.The initial mass and radius are 1.5 M ⊙ and 10 4 au respectively, and the temperature is initially 5 K throughout.The evolution under the Combined method lies almost exactly on top of the results obtained using the Stamatellos method.

Figure 5 .
Figure 5. Estimates of the column density for a collapsing sphere with an initially uniform density at  = 80000 years (left-hand panel) and at  = 145000 years (right-hand panel).Estimates are calculated using the Stamatellos method (yellow open circles), the Lombardi method (red open squares), the Combined method (blue stars), and the Modified Lombardi method (green triangles).The dashed line is the "true" column density profile, determined by integrating the density in the SPH simulation from the centre of the cloud to the surface.

Figure 6 .
Figure 6.Snapshot of the disc column density at  = 1200 years for the  disc = 0.01 M ⊙ (low-mass) disc discussed in Section 4.2.1.

Figure 7 .
Figure 7. By-eye fit to the vertical profile, at  = 75 au, of the SPH simulation of a 0.01 M ⊙ (low-mass) disc around a 0.8 M ⊙ central star.

Figure 8 .
Figure 8. Estimates for the disc scale-height at  = 75 AU of the low-mass disc calculated with the Stamatellos method (yellow circles), the Lombardi method (red open squares), the Combined method (blue stars) and the Modified Lombardi method (green triangles).The dotted line is  = 6.7 AU, the "true" scale-height obtained from the fit shown in Figure 7.

Figure 9 .
Figure 9. Estimates for the column density against  at  = 75 AU of the 0.01 M ⊙ (low-mass) disc using the Stamatellos method (yellow open circles), the Lombardi method (red open squares), the Combined method (blue stars) and the Modified Lombardi method (open green triangles).The dashed line shows the column density determined directly from the simulation by integrating the density from  to the nearest disc surface.The total column density at  = 75 AU (equivalent to the surface density) is indicated by the dotted line and will be twice the column density from the midplane to the nearest surface.We are looking to obtain estimates as close to the dashed line as possible.

Figure 10 .
Figure 10.Column density estimates against disc radius for the 0.01 M ⊙ (low-mass) disc.For the Stamatellos (yellow open circles), Lombardi (red open squares), Combined (blue stars) and Modified Lombardi (open green triangles) methods, we plot a random sample of estimates that are within 0.05 AU of the midplane ( = 0).

Figure 11 .
Figure 11.Column density renderings of the SPH simulations of high-mass discs with masses of  disc = 0.2 M ⊙ (left-hand panel) and  disc = 0.25 M ⊙ (right-hand panel).Both were evolved using the Combined method.

Figure 12 .
Figure 12.Column density against radius along a thin strip on the -axis of the  disc = 0.2 M ⊙ disc shown in the left-hand panel of Figure 11, close to the midplane.The black crosses show the column density determined by integrating the gas density in the SPH simulation, while the other symbols are column densities from the Stamatellos method (yellow open circles), the Lombardi method (red open squares), Combined method (blue stars), and Modified Lombardi method (open green triangles).

Figure 13 .
Figure 13.Column density of one of the clumps in the fragmenting high-mass disc shown in the right-hand panel of Figure 11.The dashed line is the SPH density integrated from the clump centre to  clump = 10 AU.The symbols are, again, a randomly selected sample of the column density estimates using the Stamatellos method (yellow open circles), Lombardi method (red open squares), Combined method (blue stars), and the Modified Lombardi method (open green triangles).The Combined and Modified Lombardi methods give very similar values near the centre of the clump and they agree with the Stamatellos method.

Figure 14 .
Figure 14.Cut through the disc along a radial line that includes the clump shown in Figure 13 (the  disc = 0.25 M ⊙ ) showing the column densities from the SPH simulation (black crosses) and the estimate from the Stamatellos method (yellow open circles), the Lombardi method (red open squares), the Combined method (blue stars), and the Modified Lombardi method (open green triangles).For  > 50 au, the points for the Modified Lombardi method mostly lie on top of those of the Combined method.