Short-lived gravitational instability in isolated irradiated discs

ABSTRACT

Irradiation from the central star controls the temperature structure in protoplanetary discs. Yet simulations of gravitational instability typically use models of stellar irradiation with varying complexity, or ignore it altogether, assuming heat generated by spiral shocks is balanced by cooling, leading to a self-regulated state. In this paper, we perform simulations of irradiated, gravitationally unstable protoplanetary discs using 3D hydrodynamics coupled with live Monte-Carlo radiative transfer. We find that the resulting temperature profile is approximately constant in time, since the thermal effects of the star dominate. Hence, the disc cannot regulate gravitational instabilities by adjusting the temperatures in the disc. In a |$0.1M_\odot$| disc, the disc instead adjusts by angular momentum transport induced by the spiral arms, leading to steadily decreasing surface density, and hence quenching of the instability. Thus, strong spiral arms caused by self-gravity would not persist for longer than ten thousand years in the absence of fresh infall, although weak spiral structures remain present over longer time-scales. Using synthetic images at 1.3 mm, we find that spirals formed in irradiated discs are challenging to detect. In higher mass discs, we find that fragmentation is likely because the dominant stellar irradiation overwhelms the stabilizing influence of |$P\mathrm{d}V$| work and shock heating in the spiral arms.

1 INTRODUCTION

A thin gas disc is unstable to axisymmetric perturbations when

where |$c_\text{s}$| is the sound speed, |$\kappa \approx \Omega$| is the epicyclic frequency assumed to be equal to the orbital frequency, |$\Sigma$| is the surface density (Toomre 1964), and G is the Gravitational constant. The disc is gravitationally stable at warmer temperatures (higher |$c_\text{s}$|⁠) or when the disc is less massive or more spread out (lower |$\Sigma$|⁠). For |$Q \lesssim 1.7$|⁠, the disc can become unstable to non-axisymmetric perturbations to form spiral arms (Durisen et al. 2007).

Complicating matters, neither |$c_{\rm s}$| nor |$\Sigma$| remain constant in response to the perturbation. Gammie (2001) used shearing box simulations with cooling to show that the non-linear outcome of gravitational instability is either fragmentation or self-regulation, depending on |$\beta$|⁠, the ratio of the local cooling time-scale to the orbital time, such that

This assumes that discs cool from their warmer midplane to their upper layers. Self-regulation occurs when heating from |$P\mathrm{d}V$| work and shocks is balanced by cooling, leading to a thermostat where gravitational instability controls the temperature and the disc hovers around |$Q\approx 1$|⁠, resulting in gravitoturbulence (Rice et al. 2003; Lodato & Rice 2004). Fragmentation occurs when cooling is faster than heating (⁠|$\beta \approx 3$|⁠; Gammie 2001). The exact criterion depends on a number of factors, such as the equation of state (Rice, Lodato & Armitage 2005), and numerical resolution (Meru & Bate 2010; Lodato & Clarke 2011; Meru & Bate 2012). A constant value of |$\beta$| implies that cooling is scale-free. This prescription has been used to investigate the angular momentum transport associated with a particular value of |$\beta$| (Gammie 2001; Cossins, Lodato & Clarke 2009), how dust collects in the spiral arms (Rice et al. 2004, 2006; Gibbons, Rice & Mamatsashvili 2012; Booth & Clarke 2016; Baehr & Zhu 2021; Baehr, Zhu & Yang 2022; Rowther et al. 2024), and how planets interact with gravitationally unstable discs (Baruteau, Meru & Paardekooper 2011; Malik et al. 2015; Rowther et al. 2020; Rowther & Meru 2020; Rowther, Nealon & Meru 2023).

But can such a self-regulated state exist in nature? Recent observations with the Atacama Large Millimetre/Submillimetre Array (ALMA) have revealed a subset of apparently massive protoplanetary discs with spiral arms, including Elias 2–27 (Pérez et al. 2016; Meru et al. 2017; Forgan, Ilee & Meru 2018; Hall et al. 2018; Huang et al. 2018; Veronesi et al. 2021), IM Lupi, and GM Aur (Lodato et al. 2023).

The problem with applying the |$\beta$|-cooling prescription to observed protoplanetary discs is that it ignores the radiation from the central star (Kratter & Lodato 2016). Observations reveal thermally stratified discs passively irradiated by their central star (Adams, Lada & Shu 1987; Kenyon & Hartmann 1987; Pinte et al. 2018; Law et al. 2021), with surface layers hotter than the midplane. It is unclear whether self-regulation is possible in this case, since one would infer ‘|$\beta$|-heating’, rather than ‘|$\beta$|-cooling’.

Previous authors have investigated the evolution of gravitationally unstable discs using improved thermodynamics (Boss 2002; Mejia 2004; Boley et al. 2007; Stamatellos et al. 2007; Forgan et al. 2009; Steiman-Cameron et al. 2023). A general conclusion of these studies is that the evolution is sensitive to the thermodynamic model. Some included irradiation: Cai et al. (2008) modelled irradiation from an envelope heated by the central star, by irradiating the disc with a blackbody flux of a characteristic temperature. Kratter & Murray-Clay (2011) derived a critical temperature above which a disc would be ‘irradiation-dominated’ and gravitoturbulence is impossible because the viscous heating is negligible compared to the irradiation. Meru & Bate (2010) used 3D radiation hydrodynamics with flux-limited diffusion, and prescribed the stellar irradiation by fixing the temperature at the disc surface. Rice et al. (2011) modified the cooling prescription in Gammie (2001) to include irradiation as an additional heating source. Forgan & Rice (2013), Haworth et al. (2020), and Cadman et al. (2020a) modelled stellar irradiation as a power-law temperature profile to which the disc cools. Hall et al. (2016) included irradiation as an additional term in the cooling rate. A general result of these works is that irradiation weakens gravitational instabilities, and none found thermal stratification. However, using shearing box simulations Hirose & Shi (2017) found that irradiation did not have a major impact on gravitoturbulence.

In this paper, we use 3D hydrodynamics coupled with live Monte-Carlo radiative transfer to simulate irradiated gravitationally unstable discs. This method allows us to model the disc thermodynamics realistically. Section 2 describes the simulations, Section 3 presents our results and comparison with |$\beta$|-cooling. We discuss in Section 4.

2 METHODS

We use phantom, a smoothed particle hydrodynamics (SPH) code developed by Price et al. (2018) coupled with mcfost, a Monte-Carlo radiative transfer code developed by Pinte et al. (2006, 2009).

2.1 Fiducial disc setup

We modelled the disc using |$10^6$| SPH particles between |$R_{\mathrm{in}} = 4$| au and |$R_\mathrm{out} = 100$| au in an |$0.1{\rm M}_{\odot }$| disc around a |$1{\rm M}_{\odot }$| star, modelled as a sink particle with an accretion radius of 4 au. We set up an initial surface density profile |$\Sigma = \Sigma _{0} \left({R}/{R_{\mathrm{in}}} \right)^{-0.6} f_{s}$|⁠, where |${\Sigma _{0} = 1.97 \times 10^2\, \mathrm{g\ cm}^{-2}}$|⁠, and |${f_{s} = 1-\sqrt{R_\mathrm{in}/R}}$| smooths the surface density at the inner boundary of the disc. The temperature in our simulations was set from mcfost, but we set the initial vertical extent of the disc according to |$H/R = 0.05 \left({R}/{R_{\rm in}} \right)^{0.25}$|⁠. The Cullen & Dehnen (2010) switch detects any shocks that form and generates the correct shock dissipation depending on the proximity to the shock. Near the shock, the linear shock viscosity parameter |$\alpha _{\text{AV}}\rightarrow 1$|⁠. Further away, |$\alpha _{\text{AV}}$| decreases to a minimum value of 0.1. We set the quadratic artificial viscosity parameter |$\beta _{\text{AV}}=2$| (Nealon, Price & Nixon 2015; Price et al. 2018). The disc is resolved with initial approximate mean smoothing length over disc scale height of |${\lt }h{\gt }/H = 0.2$|⁠.

2.2 Live radiative transfer

Temperatures in our simulation were calculated on-the-fly using mcfost. We included |$P\mathrm{d}V$| work and shock heating as source terms for the radiative transfer, thus allowing self-regulation of the disc in principle since spiral shocks lead to heating.

We set the stellar luminosity from the sink particle mass and a 3 Myr isochrone from Siess, Dufour & Forestini (2000). This corresponds to a stellar luminosity of |$L=0.88L_\odot$|⁠, temperature of |$T_\star = 4264$|K, and a stellar radius of |$R_\star = 1.72R_\odot$|⁠. Since we simulated only gas, for radiative transfer we assumed perfectly coupled dust with a constant dust-to-gas mass ratio of 0.01. Dust in mcfost was distributed between 0.03 and 1000μm using 100 grain sizes and a power-law exponent of |$\mathrm{d}n(s) \propto s^{-3.5} \mathrm{d}s$|⁠. In our disc these dust sizes gave Stokes numbers below 1, so the assumption that the dust is coupled to the gas is valid. We assumed spherical, homogenous, astrosilicate dust grains (Weingartner & Draine 2001). Optical properties were calculated from Mie theory.

We illuminated the disc using |$5 \times 10^9$| photons. Each mcfost Voronoi cell corresponded to one SPH particle. Our number of photons was larger than in previous studies with phantom + mcfost (Nealon, Price & Pinte 2020; Borchert et al. 2022a, b). A ratio of 100 photons for each SPH particle is usually enough to obtain the temperatures. We found more photons necessary to ensure at least some photons reach the densest part of the spiral arms in the disc. Otherwise, the particles hit the floor temperature of 2.73 K, triggering gravitational collapse, becoming optically thicker and harder for photons to reach in subsequent time-steps. The resulting negative feedback loop led to artificial fragmentation. The number of photons chosen was adjusted until every cell received enough photons for the temperatures to converge such that the floor temperature was unnecessary.

We assumed radiative equilibrium; heating and cooling equilibrate on a shorter time-scale than we simulate (see Discussion). We updated temperatures every |${1/(100\sqrt{2}) \approx 0.007071}$| of an orbit at |$R_\mathrm{out}$|⁠, i.e. every 7.071 yr. The high frequency ensured that the dynamical time was always longer than the time between temperature updates. In these simulations, the dynamical time in the inner-most regions (⁠|$R \sim 20$| au) is roughly 12 times longer than the time between temperature updates. Hence, despite the energy of the particles not evolving between temperature updates, the particles always have the correct temperatures based on their current location since they have not had the time to undergo significant dynamical evolution.

2.3 Investigating fragmentation

We modelled a |$0.25M_\odot$| disc to investigate whether fragmentation occurs when stellar irradiation is included. The disc was set up between |$R_{\mathrm{in}} = 8$| au and |$R_\mathrm{out} = 200$| au with an initial surface density profile of |$\Sigma = \Sigma _{0} \left({R}/{R_{\mathrm{in}}} \right)^{-1} f_{s}$|⁠, where |$\Sigma _0 = 3.45 \times 10^2\, \text{g cm}^{-2}$|⁠. Additionally, for the live radiative transfer we used a 1 Myr isochrone, and included accretion luminosity. All other parameters were the same as the fiducial simulation.

The changes from the fiducial simulation were numerically motivated by the need to avoid artificial fragmentation in transient features caused in the disc relaxation phase. These transients were common to both the live radiative transfer and |$\beta$|-cooling simulations. Physically, all the above changes to the setup favour stability by decreasing the amount of mass in the outer disc, and by using a more luminous star.

2.4 |$\beta$|  - Cooling

For comparison to the fiducial simulation, we also simulated an identical disc cooled by |$\beta$| - cooling. The internal energy equation in this case is given by

where we assume an adiabatic equation of state, P is the pressure, |$\rho$| is the density, and u is the specific internal energy. The first term is the |$P\mathrm{d}V$| work, |$\Pi _{\mathrm{shock}}$| is the shock viscosity heating term, and we model |$t_{\rm cool}$| using equation (2) with |$\beta =15$|⁠.

2.5 Synthetic images

We used mcfost to calculate synthetic continuum images at 1.3 mm assuming the same stellar parameters as in Section 2.2, and assuming a distance of 140 pc. The images were convolved using a 2D Gaussian with a beam size of 0.05 arcsec × 0.05 arcsec (or 7 au × 7 au at 140 pc).

3 RESULTS

3.1 Disc evolution

The top half of Fig. 1 shows the surface density of the disc evolved with radiative transfer every two orbits from 4 to 18 outer orbits. Within four outer orbits, the disc develops spiral arms due to gravitational instabilities (top left panel). Spirals are strongest in the outer regions of the disc. This behaviour is consistent with expectations of realistic discs where gravitational instabilities are weaker in the inner warmer regions of the disc (Rafikov 2005; Clarke 2009; Rice & Armitage 2009; Stamatellos & Whitworth 2009). For the next few orbits, there are no major changes to the morphology of the disc (middle and top right panels). After 12 orbits, the spiral arms are weaker (bottom left panel). As more time passes, the spiral arms continue to become weaker (middle and bottom right panels).

The weakening of gravitational instabilities is in contrast to the evolution of the disc with |$\beta$|-cooling. The bottom half of Fig. 1 shows that the |$\beta$|-cooled disc has reached a steady state with gravitational instabilities present throughout the disc. The density contrast of the spirals is larger with |$\beta$|-cooling compared to radiative transfer.

Fig. 2 compares the azimuthally averaged Toomre Q parameter (bottom panels), sound speed (middle panels), and surface density (top panels) for the two simulations (left, radiative transfer; right, |$\beta$|-cooling). We plot azimuthally averaged quantities for every outer orbit between 2 and 20 orbits. Darker lines represent later times. Fig. 3 shows the time evolution of surface density, sound speed, and Q (top, middle, and bottom panels, respectively) at |$R=77$| au (the location of minimum Q at 10 outer orbits for the disc modelled using radiative transfer). The left and right panels represent the disc evolved with radiative transfer and |$\beta$|-cooling, respectively.

Figs 2 and 3 demonstrate why the evolution of the disc is different with radiative transfer compared to |$\beta$|-cooling. With radiative transfer, the irradiation of the star dominates the temperature evolution. Thus, the temperature is fairly constant since the timeframe of these simulations is too short for the star to evolve. Whereas the surface density steadily decreases as seen in Fig. 3. The decrease in surface density is mirrored in the Q profile by its steady increase. Hence, gravitational instabilities in the disc steadily weaken. In contrast, with |$\beta$|-cooling the evolution of the disc temperature is more dynamic. From Fig. 3, the disc temperature continuously decreases before the onset of spiral structures since there is no source of heating. Once spirals form, |$P\mathrm{d}V$| work and shock heating is able to balance the cooling resulting in a disc with gravitational instabilities.

Fig. 4 shows the rate of change of the surface density between 5 and 20 outer orbits at all radial locations for the radiative transfer simulation. The regions that lose the most mass are also the most gravitationally unstable regions (⁠|${R = 80\pm 10}$|au). Hence, the disc becomes more stable as seen by Q increasing in the outer parts of the disc in Figs 2 and 3.

3.2 Disc temperatures

Fig. 5 shows the temperature profiles of the discs modelled with phantom + mcfost (left) and |$\beta$|-cooling (right). Top and bottom panels show face-on and edge-on views, respectively. Self-consistently calculating the temperatures results in two important differences compared to |$\beta$|-cooled simulations. First, the temperature is controlled by the irradiation of the star rather than the balance between cooling, and |$P\mathrm{d}V$| and shock heating. This is seen in the top left panel by the smooth temperature profile and the lack of spiral arms when the disc is modelled with radiative transfer. Whereas with |$\beta$|-cooling, the spirals are the only source of heating as seen by the spiral structure in the top right panel. Secondly, the temperature of the disc varies both radially and vertically with radiative transfer as seen in the bottom left panel. However, with |$\beta$|-cooling, the midplane temperature is radially constant (right panels). The unusual structures in the bottom right panel are due to poorly resolved regions in |$\beta$|-cooled simulations. The bright ‘x’ shape in the centre is due to a larger cavity in the inner disc. Additionally, unlike with the radiative transfer simulations, |$\beta$|-cooling also spreads the disc vertically resulting in the upper-most regions of the disc being more poorly resolved. Thus, heating due to artificial viscosity also becomes relevant in the upper parts of the disc with |$\beta$|-cooling.

3.3 Importance of the star

Fig. 6 compares the evolution of the disc modelled with radiative transfer using the three different isochrones to set the stellar luminosity. The reference simulation using a 3 Myr isochrone with |$L = 0.88L_\odot$|⁠, |$T_\star = 4264$|K, and |$R_\star = 1.72R_\odot$| is shown in the left panel. The middle and right panels show the simulations using a 2 Myr isochrone with |$L = 1.19L_\odot$|⁠, |$T_\star = 4277$|K, and |$R_\star = 1.99R_\odot$| and using a 1 Myr isochrone with |$L = 2.06L_\odot$|⁠, |$T_\star = 4280$|K, and |$R_\star = 2.62R_\odot$|⁠, respectively. The stellar parameters impact the gravitational instabilities in the disc. The more luminous the star, the hotter the disc. Hence, gravitational instability becomes weaker with increasing stellar luminosity.

3.4 Detectability of the spirals

Fig. 7 shows the synthetic images at 1.3 mm of the disc simulated with radiative transfer (left) and with |$\beta$|-cooling (right) at eight outer orbits. With radiative transfer the spirals are too weak to be detectable, even with high resolution, despite the spiral structures seen in the simulations (see Fig. 1). Whereas, with |$\beta$|-cooling, the spiral structures are apparent. The increased challenge in observing spirals when irradiation is included is consistent with Hall et al. (2016). However, if dust were included in the simulations, concentration of solids in the spirals could make it easier to detect spiral structures (Dipierro et al. 2015; Cadman et al. 2020b).

3.5 Fragmentation in higher mass discs

Fig. 8 shows the final state of a |$0.25M_\odot$| disc modelled with radiative transfer. Unlike the fiducial simulation, the disc is now massive enough that the disc cannot avoid fragmentation. In a |$\beta$|-cooled disc, fragmentation could be prevented as the disc could regulate its temperatures through |$P\mathrm{d}V$| work and shock heating in the spiral arms. However, this is negligible in an irradiated disc where the temperatures are set by the star since the thermal effects of stellar irradiation dominates. Hence, massive discs cannot stabilize through |$P\mathrm{d}V$| work and shock heating, confirming the predictions in Kratter & Murray-Clay (2011). This is in contrast to previous simulations where stellar irradiation was modelled with a power-law which found fragmentation can be completely suppressed in discs irradiated by a |$1M_\odot$| star. (Cadman et al. 2020a).

4 DISCUSSION

4.1 A different fate

With |$\beta$|-cooling, the disc cools until gravitational instability produces spiral arms. The spiral shocks heat the disc. As the disc heats up and becomes stable, |$P\mathrm{d}V$| work and shock heating lessen. The constant |$\beta$|-cooling then cools and destabilizes the disc. Eventually, an equilibrium state occurs where heating balances cooling represented by a consistent spiral structure.

Stellar irradiation changes the evolution and fate of a gravitationally unstable disc (Kratter & Murray-Clay 2011; Kratter & Lodato 2016). For a lower mass gravitationally unstable disc, the equilibrium state is now determined by the irradiation of the star and by changes to the surface density of the disc. The star supplies a constant source of heat which dominates over other forms of heating, making the irradiated disc behave more like an isothermal simulation (Kratter & Lodato 2016) where disc evolution is instead driven by changes in surface density. Angular momentum transport by the spiral arms results in viscous spreading which reduces the surface density and stabilizes the disc over time. The strongly gravitationally unstable phase is relatively short-lived (around 10 outer orbits, or 10 000 yr). Weak gravitational instabilities can persist for longer. This makes observing the short gravitationally unstable phase difficult, at least in the absence of infall, consistent with the rarity of such discs in observations.

Cadman et al. (2020a) found stellar irradiation suppressed fragmentation even for discs around solar mass stars even with disc-to-star mass ratios of 1. In contrast, we find fragmentation remains a possible fate for more massive discs (⁠|${\gt }\, 0.25M_\odot$|⁠). The dominance of stellar irradiation prevents the disc from stabilizing due to |$P\mathrm{d}V$| work and shock heating in the spiral arms. Hence, even in more massive discs, spiral structures due to gravitational instabilities may be short-lived as clump (planet)–disc interactions will determine the morphology of the disc (Rowther et al. 2020; Rowther et al.2023).

Although the above time-scale for the spiral structures is similar to Hall et al. (2019), a caveat to these results is that we do not model disc formation. In reality, discs will be more massive in the beginning and lose mass over time. Whether the initial larger disc masses lead to longer lived spiral structures or rapid fragmentation is beyond the scope of this work.

4.2 Self-regulation by infall?

As outlined by Kratter et al. (2010), the relevant thermal parameter is |$\xi \equiv \dot{M}_{\rm in} G/ c_{\rm s}^3$|⁠, where |$\dot{M}_{\rm in}$| is the mass infall rate from the environment. Thus while an irradiated disc cannot self-regulate by changing the sound speed, it may be possible to reach a steady state by fresh supply of mass (Kratter & Lodato 2016). Kratter & Murray-Clay (2011) suggested that high mass, thick discs could transport sufficient angular momentum to avoid fragmentation so long as the mass infall remains below |$\dot{M}_{\rm in} \sim c_{\rm s}^3/G$|⁠, potentially leading to steady mass accretion from excitation of global spiral modes (see Laughlin, Korchagin & Adams 1997). Examining if self-regulation could occur with mass infall in irradiation-dominated discs would make for worthwhile follow-up.

4.3 Consequences for planet formation

The stabilizing nature of stellar irradiation has consequences for planet formation in gravitationally unstable discs. Planets forming via gravitational instability will be limited to only the most massive discs. Although gas fragmentation may be limited in lower mass discs by the weaker spiral structures, the dust could instead play a bigger role in forming planets in gravitationally unstable discs. Spiral arms are pressure maxima where dust can be trapped and grow to form planetesimals (Rice et al. 2004, 2006; Gibbons et al. 2012; Booth & Clarke 2016; Elbakyan et al. 2020; Baehr & Zhu 2021; Baehr et al. 2022; Rowther et al. 2024). The amplitude and number of spiral arms impact the dynamics of dust in the spiral arms. Dust becomes more excited in stronger spiral arms, or if they are more numerous (Longarini et al. 2023a, b). Hence, weaker spiral arms with irradiation could be favourable to planetesimal formation as the dust could be less excited and collapse into clumps more easily. If planet formation occurs in these discs through gravitational collapse of gas or dust, then it must happen very early during the Class 0/1 phase when discs are massive enough to be gravitationally unstable (Xu 2022).

The weaker spiral arms could also impact the fate of any planet that does form. Rowther & Meru (2020) and Rowther et al. (2020, 2023) showed that thermodynamics strongly influences planet–disc interactions in gravitationally unstable discs. Weaker spiral arms driven by stellar irradiation could make it easier for a planet’s spiral wakes to control the evolution of the disc by more easily suppressing the spirals due to gravitational instabilities. Thus, giant planets may have better survival odds than in |$\beta$|-cooled simulations where they rapidly migrate inwards (Baruteau et al. 2011; Malik et al. 2015).

Protoplanetary discs do not evolve in isolation. The star-forming environment can shape the evolution of a disc through flybys or infall (Bate 2018). Warping can also alter disc evolution (Rowther et al. 2022). In |$\beta$|-cooled discs, the disc returns to being unstable as it cools back down until it reaches a steady state with spiral shocks. With stellar irradiation, changes in surface density determine the fate of the disc. Flybys and infall change how mass is distributed. Thus, with stellar irradiation, any interactions with the star-forming environment that change the surface density profile of the disc could permanently alter the disc’s evolution.

4.4 Importance of the star

Since stellar irradiation dominates the disc thermodynamics, the choice of stellar parameters is important for the evolution of the disc. The luminosities in this work (⁠|$0.88-2.06L_\odot$|⁠) are similar to observed stellar luminosities in discs thought to be gravitationally unstable such as Elias 2–27 and IM Lup (Huang et al. 2018) where the luminosity of the star is |$0.91L_\odot$| and |$2.57L_\odot$|⁠, respectively (Andrews et al. 2018).

4.5 Validity of assuming radiative equilibrium

The main limitation of the radiative transfer simulations is the assumption of radiative equilibrium. In reality, discs do not radiate away their energy instantly. This assumption is valid if the time-scale to reach temperature equilibration between the midplane and disc surface is sufficiently short. The time-scale for temperature equilibration is given by |$t_\text{rad} = \kappa \rho H^2/c$| (Laughlin & Bodenheimer 1994), giving

where we assume a scale height of 3 au (at R = 77 au) to match our simulations.

Compared to the dynamical time (676 yr at 77 au) and the time between temperature updates (7.071 yr), |$t_\text{rad}$| (1.1 d) is insignificant. Therefore, we can approximate radiative transfer by assuming radiative equilibrium and keeping each particle’s temperature fixed between temperature updates.

5 CONCLUSIONS

We use 3D hydrodynamics coupled with live Monte-Carlo radiative transfer to investigate the evolution of massive protoplanetary discs with irradiation from the central star. Our main findings are

• An irradiated gravitationally unstable disc does not regulate gravitational instabilities by adjusting the temperatures in the disc. The radiation of the star is the dominant source of energy, and is nearly constant over the time-scale explored here (⁠|${\lesssim }\, 100\, 000$| yr).

• In a lower mass gravitationally unstable disc (⁠|$0.1M_\odot$|⁠) the spiral structures are weaker compared to the spiral arms from |$\beta \text{-cooling}$|⁠, consistent with previous works which included stellar irradiation, despite the differing methodologies (Cai et al. 2008; Meru & Bate 2010; Forgan & Rice 2013). The weaker gravitational instabilities are due to the minimum temperature of the disc set by the irradiation of the star. Unlike with |$\beta$|-cooling, the disc cannot continuously cool until strong spiral structures form.

• Fragmentation can occur in a higher mass gravitationally unstable disc (⁠|$0.25M_\odot$|⁠). Any additional heating from |$P\mathrm{d}V$| work and shock heating in the spiral arms is negligible compared to the thermal effects of stellar irradiation. Hence, the disc cannot heat up and avoid fragmentation.

• In both scenarios, the spiral structures due to gravitational instabilities are short-lived. In lower mass discs, the evolution of the surface density is more important. The steady decrease in density due to angular momentum transport by the spiral arms results in the disc becoming more stable over time with weakening spiral structures. Whereas with higher mass discs, the disc quickly fragments resulting in a regime where clump–disc interactions will determine the morphology of the disc.

• Development of spiral arms remains possible in isolated irradiated discs, but is short-lived (⁠|${\sim }10\, 000$| yr) and thus an unlikely explanation for observed spiral structure in protoplanetary discs. The weaker spirals are also more challenging to detect.

ACKNOWLEDGEMENTS

DP thanks Kaitlin Kratter, Giuseppe Lodato, and Alison Young for useful discussions. SR thanks Monash for hosting his visit. SR and RA acknowledge support funding from the Science & Technology Facilities Council (STFC) through Consolidated Grant ST/W000857/1. FM acknowledges a Royal Society Dorothy Hodgkin Fellowship. RN acknowledges UKRI/EPSRC support via a Stephen Hawking Fellowship (EP/T017287/1). We acknowledge EU Horizon 2020 funding under the Marie Skłodowska-Curie grants 210021 and 823823 (DUSTBUSTERS). DP and CP acknowledge Australian Research Council funding via DP220103767 and DP240103290. We used the ALICE High Performance Computing Facility and the DiRAC Data Intensive service (DIaL), managed by the University of Leicester Research Computing Service on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC service at Leicester was funded by BEIS, UKRI, and STFC capital funding and STFC operations grants. DiRAC is part of the UKRI Digital Research Infrastructure.

DATA AVAILABILITY

Select snapshots of the simulations along with the data to recreate any of the simulations or the figures in the paper are available at https://doi.org/10.5281/zenodo.13760732. We utilized the following public software:

References