# Short-lived gravitational instability in isolated irradiated discs

Sahl Rowther, Daniel J. Price, Christophe Pinte, Rebecca Nealon, Farzana

Meru, Richard Alexander

► To cite this version:

Sahl Rowther, Daniel J. Price, Christophe Pinte, Rebecca Nealon, Farzana Meru, et al. Short-lived gravitational instability in isolated irradiated discs. Monthly Notices of the Royal Astronomical Society, 2024, 534, pp.2277-2285. 10.1093/mnras/stae2167. insu-04836781

HAL Id: insu-04836781

https://insu.hal.science/insu-04836781v1

Submitted on 13 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Short-lived gravitational instability in isolated irradiated discs

Sahl Rowther<sup>1,2,★</sup> Daniel J. Price<sup>1,2</sup> Christophe Pinte<sup>1,2,3</sup> Rebecca Nealon<sup>1,4,5</sup> Farzana Meru<sup>1,4,5</sup>
and Richard Alexander<sup>1</sup>

Accepted 2024 September 11. Received 2024 August 11; in original form 2024 May 30

## 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.1 M_\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 PdV work and shock heating in the spiral arms.

**Key words:** hydrodynamics – instabilities – planets and satellites: formation – protoplanetary discs.

## 1 INTRODUCTION

A thin gas disc is unstable to axisymmetric perturbations when

$$Q(R) = \frac{c_s \kappa}{\pi G \Sigma} < 1, \quad (1)$$

where $c_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_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_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

$$t_{\text{cool}} = \frac{\beta}{\Omega}. \quad (2)$$

This assumes that discs cool from their warmer midplane to their upper layers. Self-regulation occurs when heating from PdV 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.

<footnote n="1"><sup>1</sup>School of Physics and Astronomy, University of Leicester, Leicester LE1 7RH, UK</footnote>

<footnote n="2"><sup>2</sup>School of Physics and Astronomy, Monash University, VIC 3800, Australia</footnote>

<footnote n="3"><sup>3</sup>Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="4"><sup>4</sup>Centre for Exoplanets and Habitability, University of Warwick, Coventry CV4 7AL, UK</footnote>

<footnote n="5"><sup>5</sup>Department of Physics, University of Warwick, Coventry CV4 7AL, UK</footnote>

<footnote>* E-mail: sahlrowther@leicester.ac.uk</footnote>
2018; Veronesi et al. 2021), IM Lupi, and GM Aur (Lodato et al. 2023).

The problem with applying the β-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 ‘β-heating’, rather than ‘β-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 β-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_{\text{in}} = 4$ au and $R_{\text{out}} = 100$ au in an $0.1\text{M}_\odot$ disc around a $1\text{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(R/R_{\text{in}})^{-0.6} f_s$, where $\Sigma_0 = 1.97 \times 10^2 \text{ g cm}^{-2}$, and $f_s = 1 - \sqrt{R_{\text{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 (R/R_{\text{in}})^{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}} \to 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 $<h>/H = 0.2$.

### 2.2 Live radiative transfer

Temperatures in our simulation were calculated on-the-fly using MCFOST. We included $PdV$ 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.88 L_\odot$, temperature of $T_* = 4264$K, and a stellar radius of $R_* = 1.72 R_\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\mu\text{m}$ using 100 grain sizes and a power-law exponent of $dn(s) \propto s^{-3.5} \text{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_{\text{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.25 M_\odot$ disc to investigate whether fragmentation occurs when stellar irradiation is included. The disc was set up between $R_{\text{in}} = 8$ au and $R_{\text{out}} = 200$ au with an initial surface density profile of $\Sigma = \Sigma_0 (R/R_{\text{in}})^{-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 β-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 β - Cooling

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

$$ \frac{du}{dt} = -\frac{P}{\rho}(\nabla \cdot v) + \Pi_{\text{shock}} - \frac{u}{t_{\text{cool}}}, \quad (3) $$

where we assume an adiabatic equation of state, P is the pressure, ρ is the density, and u is the specific internal energy. The first term is the PdV work, Π$_{shock}$ is the shock viscosity heating term, and we model t$_{cool}$ using equation (2) with β = 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 β -cooling. The bottom half of Fig. 1 shows that the β -cooled disc has reached a steady state with gravitational instabilities present throughout the disc. The density contrast of the spirals is larger with β -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, β -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 β -cooling, respectively.

Figs 2 and 3 demonstrate why the evolution of the disc is different with radiative transfer compared to β -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 β -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, PdV 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 ± 10au). 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 β -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 β -cooled simulations. First, the temperature is controlled by the irradiation of the star rather than the balance between cooling, and PdV 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 β -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 β -cooling, the midplane temperature is radially constant (right panels). The unusual structures in the bottom right panel are due to poorly resolved regions in β -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, β -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 β -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$_{⊙}$, T$_{∗}$ = 4264K, and R$_{∗}$ = 1.72R$_{⊙}$ is shown in the left panel. The middle and right panels show the simulations using a 2 Myr isochrone with L = 1.19L$_{⊙}$, T$_{∗}$ = 4277K, and R$_{∗}$ = 1.99R$_{⊙}$ and using a 1 Myr isochrone with L = 2.06L$_{⊙}$, T$_{∗}$ = 4280K, and R$_{∗}$ = 2.62R$_{⊙}$, 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 β -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 β -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).
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Surface density of a 0.1M☉ disc around a 1M☉, 0.88 L☉ star, shown every two orbits from 4 to 18 outer orbits. Each orbit is 1000 yr. The top half shows the disc evolution using radiative transfer. Spiral arms develop due to gravitational instabilities within a couple of orbits. The disc remains strongly gravitationally unstable until around 10 outer orbits. However, this phase is short-lived. After ~14 orbits, the spiral structures weaken and are barely noticeable. Lower half shows the disc evolution using β-cooling with no stellar luminosity. The spiral arms are stronger, and the disc remains gravitationally unstable.</figcaption></figure>

## 3.5 Fragmentation in higher mass discs

Fig. 8 shows the final state of a 0.25M☉ disc modelled with radiative transfer. Unlike the fiducial simulation, the disc is now massive enough that the disc cannot avoid fragmentation. In a β-cooled disc, fragmentation could be prevented as the disc could regulate its temperatures through PdV 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 PdV 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☉ star. (Cadman et al. 2020a).

# 4 DISCUSSION

## 4.1 A different fate

With β-cooling, the disc cools until gravitational instability produces spiral arms. The spiral shocks heat the disc. As the disc heats up and
<figure><img src="image_3.png" /></figure>

**Figure 2.** Azimuthally averaged plots of the surface density (top row), sound speed (middle), and Toomre $Q$ (bottom). We show azimuthally averaged quantities for every outer orbit from 2 to 20 orbits. Darker lines represent later times. With $\beta$-cooling (right panels), the temperature of the disc evolves to compensate for changes in the surface density resulting in a steady gravitationally unstable state (see $Q$ profile). With radiative transfer (left panels) the temperature of the disc does not evolve. Since the dominant source of heating is due to stellar irradiation, which does not evolve in the timeframe of these simulations, the disc cannot cool down to compensate for the evolution of the disc surface density. Hence, as surface density decreases, $Q$ increases, resulting in a stable disc.

becomes stable, PdV 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 ($> 0.25 M_{\odot}$). The dominance of stellar irradiation prevents the disc from stabilizing due to PdV work and shock heating in the spiral arms. Hence, even in more

<figure><img src="image_4.png" /></figure>

**Figure 3.** Time evolution of the surface density, sound speed, and Toomre $Q$ at $R = 77$ au (top, middle, and bottom panels, respectively). With $\beta$-cooling (right panels), the disc cools until it forms spiral arms. Once spiral shocks are strong enough ($\sim 8$ orbits), the disc then heats, eventually reaching a steady state after 14 orbits in terms of spiral structure, see Fig. 1. With radiative transfer (left panels), the disc steadily becomes more stable, as seen by the increasing $Q$. While the disc does heat up during the strongly gravitationally unstable phase ($\sim 2–14$ orbits) due to PdV work and shock heating by the spiral arms, the temperature change is tiny compared to the decreasing surface density that stabilizes the disc.

<figure><img src="image_5.png" /></figure>

**Figure 4.** The rate of change of the surface density between 5 and 20 outer orbits at all parts of the disc for the radiative transfer simulation. The most gravitationally unstable parts of the disc lose the most amount of mass (at $R = 80 \pm 10$ au). This mass-loss results in the disc becoming more stable over time.

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
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Cross-section slices of the face-on and edge-on views of the temperature structure of a $0.1 M_{\odot}$ disc at nine outer orbits. In contrast to the density structure, spirals are barely visible in temperature with irradiation (top left panel). PdV work and shock heating is negligible compared to the radiation of the star. Whereas with $\beta$-cooling, the temperature structure is strongly correlated to the density structure (top right panel) since the spirals are the only source of heating. The irradiated disc is also thermally stratified and varies with radius, consistent with observations (bottom left panel). The vertical temperature structure with $\beta$-cooling (bottom right panel) is unusual due to the disc also spreading vertically. Thus, heating due to artificial viscosity is also relevant for the upper-most regions which is more poorly resolved. In the midplane, the disc temperature is uniform with radius unlike with radiative transfer. The bright ‘x’ shape in the inner-most regions is also due to a large inner cavity with very few particles.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Surface density of a $0.1 M_{\odot}$ disc after six orbits. From left to right, the luminosity of the star is 0.88, 1.19, and 2.06 $L_{\odot}$ using a 3, 2, and 1 Myr isochrone, respectively. Gravitational instabilities become weaker with increasing stellar luminosity.</figcaption></figure>

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 M_{\text{in}} G/c_s^3$, where $M_{\text{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}_{\text{in}} \sim c_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.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Synthetic images at 1.3 mm of the disc simulated with radiative transfer (left) and with β-cooling (right). The beam size (7 au × 7 au) is represented by the white circle. With radiative transfer, the spiral structures seen in the simulations would be challenging to detect in observations. Whereas with β-cooling, the spiral structures are much easier to detect.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Fragmentation of a $0.25M_{\odot}$ disc modelled with radiative transfer. Since the irradiation of the star controls the disc temperatures, the disc is unable to stabilize and avoid fragmentation through PdV work and shock heating in the spiral arms.</figcaption></figure>

## 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 β-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 β-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_☉$) 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_☉$ and $2.57L_☉$, 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

$$t_{\text{rad}} = 1.1 \, \text{d} \left( \frac{\kappa}{5 \, \text{cm}^2 \, \text{g}^{-1}} \right) \left( \frac{\rho}{2 \times 10^{-13} \, \text{g} \, \text{cm}^{-3}} \right) \left( \frac{H}{3.0 \, \text{au}} \right)^2, \quad (4)$$

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

(i) 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).

(ii) In a lower mass gravitationally unstable disc ($0.1M_☉$) the spiral structures are weaker compared to the spiral arms from $\beta$-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.

(iii) Fragmentation can occur in a higher mass gravitationally unstable disc ($0.25M_☉$). Any additional heating from $PdV$ 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.

(iv) 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.

(v) 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:

<table><tr><td>PHANTOM</td><td>https://github.com/danieljprice/phantom<br>(Price et al. 2018)</td></tr><tr><td>MCFOST</td><td>https://github.com/cpinte/mcfost<br>(Pinte et al. 2006, 2009)</td></tr><tr><td>SPLASH</td><td>https://github.com/danieljprice/splash<br>(Price 2007)</td></tr><tr><td>SARRACEN</td><td>https://github.com/ttricco/sarracen<br>(Harris & Tricco 2023)</td></tr></table>

## REFERENCES

Adams F. C., Lada C. J., Shu F. H., 1987, ApJ, 312, 788

Andrews S. M. et al., 2018, ApJ, 869, L41

Baehr H., Zhu Z., 2021, ApJ, 909, 135

Baehr H., Zhu Z., Yang C.-C., 2022, ApJ, 933, 100

Baruteau C., Meru F., Paardekooper S.-J., 2011, MNRAS, 416, 1971

Bate M. R., 2018, MNRAS, 475, 5618

Boley A. C., Durisen R. H., Nordlund Å., Lord J., 2007, ApJ, 665, 1254

Booth R. A., Clarke C. J., 2016, MNRAS, 458, 2676

Borchert E. M. A., Price D. J., Pinte C., Cuello N., 2022a, MNRAS, 510, L37

Borchert E. M. A. et al., 2022b, MNRAS, 517, 4436

Boss A. P., 2002, ApJ, 576, 462

Cadman J., Rice K., Hall C., Haworth T. J., Biller B., 2020a, MNRAS, 492, 5041

Cadman J., Hall C., Rice K., Harries T. J., Klaassen P. D., 2020b, MNRAS, 498, 4256

Cai K. et al., 2008, ApJ, 673, 1138

Clarke C. J., 2009, MNRAS, 396, 1066

Cossins P., Lodato G., Clarke C. J., 2009, MNRAS, 393, 1157

Cullen L., Dehnen W., 2010, MNRAS, 408, 669

Dipierro G., Pinilla P., Lodato G., Testi L., 2015, MNRAS, 451, 974
Durisen R. H. et al., 2007, Protostars and Planets V, eds, Reipurth B, Jewitt D, Keil K, p.607

Elbakyan V. G. et al., 2020, A&A, 637, A5

Forgan D., Rice K., 2013, MNRAS, 430, 2082

Forgan D., Rice K., Stamatellos D., Whitworth A., 2009, MNRAS, 394, 882

Forgan D. H., Ilee J. D., Meru F., 2018, ApJ, 860, L5

Gammie C. F., 2001, ApJ, 553, 174

Gibbons P. G., Rice W. K. M., Mamatsashvili G., 2012, MNRAS, 426, 1444

Hall C., Forgan D., Rice K., Harries T. J., Klaassen P. D., Biller B., 2016, MNRAS, 458, 306

Hall C. et al., 2018, MNRAS, 477, 1004

Hall C., Dong R., Rice K., Harries T. J., Najita J., Alexander R., Brittain S., 2019, ApJ, 871, 228

Harris A., Tricco T. S., 2023, J. Open Source Softw., 8, 5263

Haworth T. J., Cadman J., Meru F., Hall C., Albertini E., Forgan D., Rice K., Owen J. E., 2020, MNRAS, 494, 4130

Hirose S., Shi J.-M., 2017, MNRAS, 469, 561

Huang J. et al., 2018, ApJ, 869, L43

Kenyon S. J., Hartmann L., 1987, ApJ, 323, 714

Kratter K., Lodato G., 2016, ARA&A, 54, 271

Kratter K. M., Murray-Clay R. A., 2011, ApJ, 740, 1

Kratter K. M., Matzner C. D., Krumholz M. R., Klein R. I., 2010, ApJ, 708, 1585

Laughlin G., Bodenheimer P., 1994, ApJ, 436, 335

Laughlin G., Korchagin V., Adams F. C., 1997, ApJ, 477, 410

Law C. J. et al., 2021, ApJS, 257, 4

Lodato G., Clarke C. J., 2011, MNRAS, 413, 2735

Lodato G., Rice W. K. M., 2004, MNRAS, 351, 630

Lodato G. et al., 2023, MNRAS, 518, 4481

Longarini C. et al., 2023a, MNRAS, 519, 2017

Longarini C. et al., 2023b, MNRAS, 522, 6217

Malik M., Meru F., Mayer L., Meyer M., 2015, ApJ, 802, 56

Mejia A. C., 2004, PhD thesis, Indiana University, Bloomington

Meru F., Bate M. R., 2010, MNRAS, 406, 2279

Meru F., Bate M. R., 2012, MNRAS, 427, 2022

Meru F. et al., 2017, ApJ, 839, L24

Nealon R., Price D. J., Nixon C. J., 2015, MNRAS, 448, 1526

Nealon R., Price D. J., Pinte C., 2020, MNRAS, 493, L143

Pérez L. M. et al., 2016, Science, 353, 1519

Pinte C., Ménard F., Duchêne G., Bastien P., 2006, A&A, 459, 797

Pinte C. et al., 2009, A&A, 498, 967

Pinte C. et al., 2018, A&A, 609, A47

Price D. J., 2007, Publ. Astron. Soc. Austr., 24, 159

Price D. J. et al., 2018, Publ. Astron. Soc. Austr., 35, e031

Rafikov R. R., 2005, ApJ, 621, L69

Rice W. K. M., Armitage P. J., 2009, MNRAS, 396, 2228

Rice W. K. M. et al., 2003, MNRAS, 339, 1025

Rice W. K. M. et al., 2004, MNRAS, 355, 543

Rice W. K. M., Lodato G., Armitage P. J., 2005, MNRAS, 364, L56

Rice W. K. M. et al., 2006, MNRAS, 372, L9

Rice W. K. M., Armitage P. J., Mamatsashvili G. R., Lodato G., Clarke C. J., 2011, MNRAS, 418, 1356

Rowther S., Meru F., 2020, MNRAS, 496, 1598

Rowther S. et al., 2020, ApJ, 904, L18

Rowther S., Nealon R., Meru F., 2022, ApJ, 925, 163

Rowther S., Nealon R., Meru F., 2023, MNRAS, 518, 763

Rowther S. et al., 2024, MNRAS, 528, 2490

Siess L., Dufour E., Forestini M., 2000, A&A, 358, 593

Stamatellos D., Whitworth A. P., 2009, AIP Conf. Proc. Vol. 1094, 15th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun. Am. Inst. Phys., New York, p. 557

Stamatellos D. et al., 2007, A&A, 475, 37

Steiman-Cameron T. Y. et al., 2023, ApJ, 958, 139

Toomre A., 1964, ApJ, 139, 1217

Veronesi B. et al., 2021, ApJ, 914, L27

Weingartner J. C., Draine B. T., 2001, ApJ, 548, 296

Xu W., 2022, ApJ, 934, 156

This paper has been typeset from a TeX/LATEX file prepared by the author.