2. Analytical framework to investigate gravitational instabilities

In this section, we outline the theoretical framework we have used to evaluate the conditions under which a disc may become gravitationally unstable and fragment. We describe the Toomre stability criterion, cooling efficiency conditions, and various models for the disc surface density and thermal structure. These provide the basis for evaluating whether a given disc structure is susceptible to planet formation by fragmentation.

Gravitational instabilities occur in discs if the self-gravity of the disc is able to overcome the effects of thermal pressure and rotation. This is quantified by the Toomre parameter Q (Toomre Reference Toomre1964):

(1) \begin{equation}Q(r)=\frac{c_\mathrm{s}(r)\Omega(r)}{\pi G \Sigma(r)} \leq 1\end{equation}

where $c_\mathrm{s}$ is the sound speed, $\Omega$ is the Keplerian angular frequency of the disc, $\Sigma$ is the disc’s surface density, G is the gravitational constant, and r is the radial distance from the centre of the disc. While the classical Toomre criterion requires $Q \leq 1$ for instability against axisymmetric perturbations, simulations indicate that non-axisymmetric instabilities and spiral structures can emerge for $Q \leq 1.4$ –1.7 (Helled et al. Reference Helled2014). The sound speed is related to the temperature, T(r), by the following expression:

(2) \begin{equation}c_\mathrm{s}(r)=\sqrt{\eta \frac{k_\mathrm{B}T(r)}{\mu m_\mathrm{p}}},\end{equation}

where $k_\mathrm{B}$ and $m_\mathrm{p}$ are the Boltzmann constant and proton mass, respectively, and we have assumed the adiabatic index $\eta=1$ .Footnote ^a We adopt a mean molecular weight of $\mu=2$ corresponding to molecular hydrogen gas.

For the optically thick regions of the disc:

(3) \begin{equation}T(r) = \left(\frac{\alpha}{4}\right)^{1/4}\left(\frac{R_\mathrm{star} }{r}\right)^{1/2}T_\mathrm{star}\end{equation}

(Chiang & Goldreich Reference Chiang and Goldreich1997), where $\alpha\approx 0.4R_\mathrm{star} /r$ is the grazing angle of the flux of stellar radiation incident on the disc and is typically much smaller than unity, $R_\mathrm{star}$ is the photospheric radius of the star and $T_\mathrm{star}$ is the surface temperature of the star. This gives $T\propto r^{-0.75}$ . This temperature profile is valid for optically thick discs in radiative equilibrium.

Using this, we propose two temperature profiles for two regions inside and outside of the sublimation radius defined as:

(4) \begin{equation} R_\mathrm{subl} =\frac{1}{2}R_\mathrm{star} \left(\frac{T_\mathrm{star}}{T_\mathrm{subl}} \right)^2,\end{equation}

where we assume $T_\mathrm{subl} =1\,200$ K is the dust sublimation temperature in post-AGB discs (Kluska et al. Reference Kluska2019). Assuming radiative equilibrium, the temperature of the disc in the midplane will be

(5) \begin{equation}\begin{aligned}T(r) = T_\mathrm{star} \left(\frac{R_\mathrm{star} }{2r}\right)^{0.5}\, \text{for} \, r < R_\mathrm{subl} ,\\[8pt]T(r) = \frac{A}{r^{0.75}}\, \text{for} \, r > R_\mathrm{subl} .\end{aligned}\end{equation}

where A is a constant that can be calculated by making use of the fact that the two equations should predict the same temperature, $T_\mathrm{subl}$ , at $r=R_\mathrm{subl}$ :

(6) \begin{equation} A= T_\mathrm{subl}R_\mathrm{subl}^{0.75}.\end{equation}

Using this we have:

(7) \begin{equation}\begin{aligned}T(r) = T_\mathrm{star} \left(\frac{R_\mathrm{star} }{2r}\right)^{0.5}\, \text{for} \, r < R_\mathrm{subl} ,\\[8pt]T(r) = T_\mathrm{subl} \left(\frac{R_\mathrm{subl}}{r}\right)^{0.75}\, \text{for} \, r > R_\mathrm{subl} .\end{aligned}\end{equation}

2.1 Cooling efficiency

In order for gravitational instabilities to lead to disc fragmentation, the fragments must cool sufficiently rapidly such that thermal pressure does not halt the collapse. This criterion is expressed by the following inequality (Gammie Reference Gammie2001):

(8) \begin{equation} t_\mathrm{cool}<\frac{\beta}{\Omega},\end{equation}

where $t_\mathrm{cool}$ is the cooling time and $\beta\approx3$ is a constant. It has been shown in simulations by Mayer et al. (Reference Mayer, Quinn, Wadsley and Stadel2002) that protoplanetary discs with minimum Toomre parameters of $Q_\mathrm{min}\sim 1.4$ start forming spiral structures which can fragment and form several planets within $\sim 1000$ yr, if the cooling is efficient enough. The ratio of the thermal energy of the disc to its emissivity can be used to estimate cooling times (Rafikov Reference Rafikov2005; Kratter, Murray-Clay, & Youdin Reference Kratter, Murray-Clay and Youdin2010):

(9) \begin{equation} t_\mathrm{cool} = \frac{ c_s^2 \Sigma}{(\eta -1)}\frac{f(\tau)}{\sigma T^4},\end{equation}

where $\sigma$ is the Stefan–Boltzmann constant, and $f(\tau)=\tau+1/\tau$ is a function to quantify the efficiency of disc cooling, which takes the vertical optical depth to be $\tau=\kappa \Sigma/2$ , where $\kappa$ is the opacity.

2.2 The surface density profile

The surface density profile, $\Sigma(r)$ , plays a key role in determining the gravitational stability of the disc. A commonly adopted empirical form is a truncated power-law (Schleicher & Dreizler Reference Schleicher and Dreizler2014):

(10) \begin{equation} \Sigma(r) = \Sigma_0\left(\frac{r_\mathrm{out}}{r}\right)^{n},\end{equation}

where r is the orbital distance, $r_\mathrm{out}$ is the outer disc’s radius, n is the power law index that determines the steepness of the surface density profile and $\Sigma_0$ is a normalisation constant set by the disc’s mass, $M_\mathrm{disc}$ , n, and the outer radius:

(11) \begin{equation} \Sigma_0 =\frac{(2-n)M_\mathrm{disc} }{2\pi r_\mathrm{out} ^{2}},\end{equation}

which is only valid for $n<2$ . This type of function is used to fit the surface densities of both YSO discs (Speedie et al. Reference Speedie2023), and post-AGB discs (Corporaal et al. Reference Corporaal, Kluska, Van Winckel, Kamath and Min2023b). For a fixed disc mass, decreasing $r_\mathrm{out}$ increases the overall surface density, effectively making the disc more compact (Figure 1).

Figure 1. Power-law surface density functions with $n=1$ as a function of orbital distance for $M_\mathrm{ disc}=0.15$ M $_\odot$ for different values of $r_\mathrm{out}$ , to demonstrate why the surface density is highly sensitive to $r_\mathrm{out}$ .

In Section 2.2.2, we introduce an alternative surface density formulation that addresses some of the limitations of this parametrisation, including the assumption of a sharply defined outer edge.

2.2.1 Discs with inner cavities

Post-AGB discs have inner cavities, evacuated regions, or holes centred on the central binary. These are not considered by Equation (10). Moreover, Equation (11) gives negative surface density values for power law indices $n\ge2$ . To account for this if our disc has a cavity, we use Equation (11) but set $\Sigma(r)=0$ for $r < r_\mathrm{in} $ , where $r_\mathrm{in} $ is the inner radius of the disc. The normalisation constant $\Sigma_0$ is consequently adjusted to account for a finite inner radius:

(12) \begin{equation}\begin{aligned}\Sigma_0 &=\frac{(2-n)M_\mathrm{disc} }{2\pi(r_\mathrm{out} ^{2}-r_\mathrm{in} ^{2}(\frac{r_\mathrm{out} }{r_\mathrm{in} })^n)},\: \mathrm{for} \: n\neq 2,\\[10pt]\Sigma_0 &=\frac{M_\mathrm{disc} }{2\pi \ln({r_\mathrm{out} /r_\mathrm{in} }) r_\mathrm{out} ^2},\: \mathrm{for} \: n=2.\end{aligned}\end{equation}

which can be used for all values of n. For $n<2$ and $r_\mathrm{in} =0$ we recover Equation (11). Physically motivated values for $r_\mathrm{in}$ include the dust sublimation radius or $\sim$ 3 times the binary separation, corresponding to the innermost stable orbit in circumbinary systems (Artymowicz & Lubow Reference Artymowicz and Lubow1994).

2.2.2 Surface density profiles of viscous discs

While power-law models are commonly used to describe disc surface density profiles due to their simplicity, they are not physically motivated and rely on poorly defined parameters such as the outer radius, $r_\mathrm{out}$ . Observationally, different diagnostics yield inconsistent estimates of $r_\mathrm{out}$ – for instance, optically thin dust continuum observations often indicate smaller disc sizes than optically thick molecular line tracers (e.g. Andrews et al. Reference Andrews, Wilner, Hughes, Qi and Dullemond2009).

A more physically motivated model, widely adopted for YSO discs, is based on viscous disc evolution theory. This approach originates from the seminal works of Lynden-Bell & Pringle (Reference Lynden-Bell and Pringle1974) and Pringle (Reference Pringle1981) and was formulated in a practical surface density profile by Hartmann et al. (Reference Hartmann, Calvet, Gullbring and D’Alessio1998). In this model, the disc’s kinematic viscosity is assumed to follow a radial power law with index $\gamma$ : $\nu \propto r^\gamma$ . Under this assumption, the surface density profile takes the form:

(13) \begin{equation} \Sigma(r) = \Sigma_0 \left( \frac{r}{r_\mathrm{ch}} \right)^{-\gamma} \exp \left[ -\left( \frac{r}{r_\mathrm{ch}} \right)^{2 - \gamma} \right],\end{equation}

where $r_\mathrm{ch}$ is the characteristic radius and $\Sigma_0$ is the normalisation constant. Unlike the outer radius in power-law models, $r_\mathrm{ch}$ is well defined within the viscous evolution framework and avoids the sensitivity associated with defining a sharp truncation.

Following Booth & Ilee (Reference Booth and Ilee2020), the normalization constant $\Sigma_0$ for this profile can be expressed in terms of the total disc mass $M_\mathrm{disc}$ and $r_\mathrm{ch}$ as:

(14) \begin{equation}\Sigma_0 = (2 - \gamma) \frac{M_\mathrm{disc}}{2\pi r_\mathrm{ch}^2} \exp \left[ \left( \frac{r_\mathrm{in}}{r_\mathrm{ ch}} \right)^{2 - \gamma} \right],\end{equation}

where $r_\mathrm{in}$ is the inner radius of the disc, if the disc has one.

This exponentially tapered surface density profile has been successfully applied to fit the data from the observation of many YSO discs; however, it has not yet been widely applied to characterise post-AGB discs.

3. Gravitational instability in YSO discs

Before applying our gravitational instability framework to the less explored environments of post-AGB discs, we first test its applicability using a set of well studied YSO systems. Many of these systems show direct or indirect features – such as spirals, clumps, or large gaps – that have been associated with gravitational instabilities (Cadman, Rice, & Hall Reference Cadman, Rice and Hall2021). By applying our analytical model to these systems, we aim to establish a benchmark for interpreting results in post-main sequence discs.

The YSO systems we consider are:

1. 1. AB Auriga: a system with a gravitationally unstable disc and a detected planet (Currie et al. Reference Currie2022).

2. 2. HL Tauri: a disc showing concentric rings with gaps, possibly indicating embedded planets (ALMA Partnership et al. Reference Partnership2015; Booth & Ilee Reference Booth and Ilee2020).

3. 3. Elias 2-27: a system with spiral structures and density waves, suggesting possible gravitational instability (Paneque-Carreño et al. Reference Paneque-Carreño2021).

4. 4. L1448 IRS3B: a triple protostar system embedded in a massive, gravitationally unstable disc with prominent spiral arms and a massive clump, as revealed by ALMA observations (Reynolds et al. Reference Reynolds2021).

5. 5. GQ Lupi: a T Tauri star with a 10–40 M $_\mathrm{Jup}$ companion at $\sim110$ au. While the current disc mass is relatively low ( $\sim$ 3–4 M $_\mathrm{Jup}$ ), the wide separation of the companion suggests it may have formed via gravitational instability during an earlier, more massive phase of the disc (Wu et al. Reference Wu2017). We note that the companion may lie in the brown dwarf regime, though its wide separation and possible disc origin make it relevant to gravitational instability-driven formation scenarios.

For comparison, we also include two YSO systems with active planet formation, but without prominent indicators of disc instability:

1. 1. PDS 70: a low-mass disc system with two directly imaged planets (Keppler et al. Reference Keppler2019).

2. 2. HD 169142: a system with a disc that shows gaps hinting at two planets, one of which has recently been confirmed through direct imaging and ALMA molecular line observations (Fedele et al. Reference Fedele2017; Hammond et al. Reference Hammond2023).

In Table 1 we summarise the relevant properties of the YSOs used in our analysis, based on studies that have modelled the disc using observational data. For these systems, we have adopted the density profiles based on what was used in the respective studies, where some used the surface density profiles of viscous discs (Equation 13), while others used a power law surface density profile (Equation 12; see Table 1 and Appendix A for more details). For HL Tauri and Elias 2-27, where stellar temperatures and radii are not available, we adopt plausible values consistent with their estimated stellar masses. We did not find an estimation for the inner radius of the disc of PDS 70 in the literature, hence we set its inner radius to $r_\mathrm{in}=0$ .

Table 1. Stellar and disc parameters for the benchmark YSO systems used in our gravitational instability analysis (see Appendix A for sources).

^a No planets confirmed via direct imaging or RV, though disc structures suggestive of embedded planets exist (ALMA Partnership et al. Reference Partnership2015).

^b The mass of the clump observed in L1448 IRS3B is high enough to be a protostar.

^c Two values exist in the literature; see Appendix A; in this analysis, we have plotted the Toomre curve of this disc using both disc mass values in Figure 2.

^d No value available – we set the inner radius to 0.

^e For all systems except AB Aurigae and L1448 IRS3B, the tabulated value in this row is the characteristic radius. For AB Aurigae and L1448 IRS3B, the tabulated values are their outer radii.

^f Indirectly calculated from the luminosity of the protostar, assuming a blackbody spectrum.

In Figure 2, we plot the Toomre parameter, Q (see Section 2.1), as a function of disc radius for each system. While our focus here is on the Toomre criterion as a first-order diagnostic, we note that gravitational fragmentation also requires rapid cooling (see Section 2.1), which we do not assess in this section. We find that only two of the systems are Toomre-unstable at present (AB Aurigae with the disc mass reported in Speedie et al. (Reference Speedie2023), and L1448 IRS3B) indicating that gravitational instability is not generally active under current disc conditions for most of the systems, although some of them such as HL Tauri, GQ Lup, and Elias 2-27 are close to the unstable regime. We note that Toomre analyses have been carried out in the past for HL Tauri (Booth & Ilee Reference Booth and Ilee2020) and Elias 2-27 (Pérez et al. Reference Pérez2016), and the resulting Toomre profiles in these studies match the profiles we found for them. As a further benchmark of our method, we compute the Toomre Q profiles of three protoplanetary discs with published Q profiles, as detailed in Appendix B.

Figure 2. Toomre parameter Q versus disc radius for the YSOs tabulated in Table 1. Only AB Aurigae with the disc mass reported in Speedie et al. (Reference Speedie2023), and L1448 IRS3B fall below $Q=1$ .

Our findings suggest that planet formation via gravitational instabilities may have occurred in these YSOs earlier in the disc’s evolution. We revisit this possibility in Section 5, where we explore whether these systems could have been gravitationally unstable at earlier stages.

4. Gravitational instability in typical post-AGB discs

Having used our gravitational instability framework on well-studied YSO systems, we now apply the same methodology to the post-AGB binary population.

Figure 3. The Toomre parameter Q as a function of radius for typical post-AGB discs, assuming $M_1+M_2=1.5$ M $_\odot$ , $R_\mathrm{star} =100$ R $_\odot$ , $T_\mathrm{star} =5\,000$ K, $r_\mathrm{in}=3$ au, and $r_\mathrm{out}=100$ au as a lower limit of post-AGB disc truncation radii. The index n is the power-law index of the surface density profile for each of the curves. Even in this relatively compact scenario, we see that the discs are generally stable.

Estimates of post-AGB disc gas masses, based on ALMA observations of molecular line emission (e.g. CO), typically fall in the range of $10^{-4}$ – $10^{-2}$ M $_{\odot}$ (e.g. Bujarrabal et al. Reference Bujarrabal, Alcolea, Van Winckel, Santander-Garca and Castro-Carrizo2013; Gallardo Cava et al. Reference Gallardo Cava2021). In contrast, disc dust masses inferred from radiative transfer modelling of near and mid-infrared interferometric data and SEDs imply significantly higher total disc masses – reaching up to $0.1$ – $0.2$ M $_\odot$ in a few individual cases (e.g. Kluska et al. Reference Kluska2018; Corporaal et al. Reference Corporaal, Kluska, Van Winckel, Kamath and Min2023b). These estimates are sensitive to assumptions about grain properties, vertical structure, and disc inclination and are often subject to model degeneracies. Interpreting such high dust masses with a canonical gas-to-dust ratio of $\sim$ 100Footnote ^b would imply total disc masses of order $10^{-1}$ M $_{\odot}$ , in tension with dynamical constraints and gas-based mass estimates (Hillen et al. Reference Hillen2015; Kluska et al. Reference Kluska2018). To account for the uncertainties in total disc mass estimates – arising from the discrepancy between gas and dust measurements and the unknown gas-to-dust ratio – we adopt a conservative range of $0.01$ – $0.1$ M $_\odot$ for the total (gas + dust) disc mass. We assess gravitational stability across this range using a power-law surface density profile, consistent with standard characterisations of post-AGB discs from observations (Section 2.2). Typical truncation radii lie between 100 and 500 au; for our fiducial case, we adopt $r_\mathrm{out} = 100$ au to represent the more compact end of the distribution, where surface densities are highest and instability most likely.

We assume a central binary mass of $M_1+M_2=1.5$ M $_\odot$ , consistent with typical post-AGB binaries (e.g. Oomen et al. Reference Oomen2018), and adopt $R_\mathrm{star} = 100$ R $_\odot$ and $T_\mathrm{star} = 5\,000$ K for the post-AGB star. Assuming a binary separation of 1 au, we set the disc inner radius to the expected dynamical truncation radius of $r_\mathrm{in}=3$ au following (Artymowicz & Lubow Reference Artymowicz and Lubow1994).

In Figure 3, we show the Toomre parameter Q as a function of radius for discs with varying surface density power-law indices. Post-AGB discs appear generally stable across this range. For power-law indices $n\geq 2$ , the surface density declines steeply with radius, and stability increases accordingly due to the rapidly decreasing mass at large radii. However, at smaller radii, discs with steeper profiles ( $n>2$ ) can approach marginal instability.

To reflect the observed structure and ensure physical consistency, we adopt a power-law (Equation 10) to parametrise the surface densities of post-AGB discs following how they are typically modelled in literature, and we also include an inner cavity in our analysis (Section 2.2, Equation 12) to avoid this divergence and to remain consistent with the typical observed structure of post-AGB discs. This is in contrast with how we have parametrised YSO discs (Equation 13) in the preceding section.

Even under compact, high-density configurations, the resulting Q values remain above the fragmentation threshold in most regions, indicating that typical post-AGB discs are gravitationally stable under present-day conditions.

4.1 Testing the stability of the massive disc in IRAS 08544-4431

IRAS 08544-4431 is a post-AGB system with an infrared excess, indicating the presence of an optically thick circumbinary disc. Its inner rim lies at the dust sublimation radius (Corporaal et al. Reference Corporaal, Kluska, Van Winckel, Kamath and Min2023b), consistent with a disc in an early evolutionary stage.Footnote ^c

For stellar parameters, we adopt $T_\mathrm{star}=7\,250$ K, $R_\mathrm{star} = 65$ R $_\odot$ , $M_1=0.75$ M $_\odot$ , and $M_2=1.65$ M $_\odot$ . For the disc, we adopt an outer radius of $r_\mathrm{out} =175$ au and the inner radius of three times the binary orbital separation ( $a=1.74$ au), giving $r_\mathrm{in} =5.22$ au (Kluska et al. Reference Kluska2018).

We explore a plausible range of disc masses: at the high end we adopt $M_\mathrm{disc} = 0.1$ M $_\odot$ which we adopt from modelling carried out by Corporaal et al. (Reference Corporaal, Kluska, Van Winckel, Kamath and Min2023b) based on their observations of the dust component, who found a dust mass of $M_\mathrm{disc} = 0.001$ M $_\odot$ under the assumption of a gas-to-dust ratio of 100.Footnote ^d At the lower end, we adopt the gas mass estimated by Bujarrabal et al. (Reference Bujarrabal2018) from ALMA CO Observations: $M_\mathrm{disc}=0.02$ M $_\odot$ .

Corporaal et al. (Reference Corporaal, Kluska, Van Winckel, Kamath and Min2023b) modelled the surface density using a piecewise power law, with a break at $r_\mathrm{ mid}=3~R_\mathrm{in}$ , and power-law indices of $n_\mathrm{in}=1.5$ and $n_\mathrm{out}=-1$ in the inner and outer regions, respectively:

(15) \begin{equation}\begin{aligned}\Sigma(r) = \Sigma_0\left(\frac{r_\mathrm{mid}}{r}\right)^{n_\mathrm{in}}\, \text{for} \, r_\mathrm{in} < r < r_\mathrm{mid} ,\\\Sigma(r) = \Sigma_0\left(\frac{r_\mathrm{mid}}{r}\right)^{n_\mathrm{out}}\, \text{for} \, r_\mathrm{out} > r > r_\mathrm{mid} .\end{aligned}\end{equation}

The surface density normalisation, $\Sigma_0$ , is then obtained from the total disc mass using (not valid for $n=2$ ):

(16) \begin{equation} M_\mathrm{disc}=\Sigma_0 \left( \frac{2\pi r_\mathrm{mid} ^{n_\mathrm{in}}(r_\mathrm{mid} ^{2-n_\mathrm{in}}-r_\mathrm{in} ^{2-n_\mathrm{in}})}{(2-n_\mathrm{in}) }+\frac{2\pi r_\mathrm{mid} ^{n_\mathrm{out}}(r_\mathrm{out} ^{2-n_\mathrm{ out}}-r_\mathrm{mid} ^{2-n_\mathrm{out}})}{(2-n_\mathrm{out}) }\right)\end{equation}

Figure 4 shows the resulting Toomre profiles for both disc mass estimates. If the higher mass inferred by Corporaal et al. (Reference Corporaal, Kluska, Van Winckel, Kamath and Min2023b) is accurate, the disc reaches $Q < 4$ in its outer regions – approaching but not exceeding the fragmentation threshold. This places IRAS 08544 $-$ 4431 among the most gravitationally unstable post-AGB discs in our sample, though it remains above the threshold for fragmentation. By contrast, the lower mass estimate from Bujarrabal et al. (Reference Bujarrabal2018) yields a Toomre profile that stays well within the stable regime across all radii. These results highlight the sensitivity of the Toomre parameter to disc mass assumptions and underscore the need for more robust constraints on gas and dust masses in evolved systems.

Figure 4. The Toomre parameter Q as a function of radius for the post-AGB system IRAS 08544-4431. The shaded region indicates the uncertainty in the total disc mass, based on estimates from Bujarrabal et al. (Reference Bujarrabal2018) and Corporaal et al. (Reference Corporaal, Kluska, Van Winckel, Kamath and Min2023b).

4.2 An exploration of the parameter space for gravitational instability across different disc masses and sizes

To complete our investigation, we explore the gravitational stability of discs across a range of mass and size configurations. Specifically, we plot the Toomre parameter Q as a function of disc mass and outer radius to assess whether a given system – such as a post-AGB binary with a more massive disc in the past – could have supported planet formation via gravitational instability.

We consider three central star cases: post-AGB binaries (Figure 5), main-sequence stars (Figure 6, top panel), and white dwarfs (Figure 6, bottom panel). These plots provide a reference framework for gauging the gravitational stability of discs around stars at different evolutionary stages.

Figure 5. Contours showing variation of the Toomre parameter Q as a function of disc radius and mass, evaluated at the outermost radius $r_\mathrm{out}$ for post-AGB discs with surface density exponent $n=0.5$ and $n=1$ (left and centre panels), and at the innermost radius $r_\mathrm{in}$ for $n=3$ (right panel). All calculations adopt $R_\mathrm{star}=100$ R $_\odot$ , $T\mathrm{star}=5\,000$ K, and $M_\mathrm{cnt}=1.5$ M $_\odot$ .

In all cases, we adopt a power-law surface density profile (Equation 10) with $n=1$ , and additionally show $n=0.5$ for post-AGB systems. For steep profiles ( $n=3$ ), where Q decreases with radius, we evaluate Q at the inner radius ( $r_\mathrm{in}$ ) instead of $r_\mathrm{out}$ , as shown in the right panel of Figure 5. We adopt fixed parameters of $r_\mathrm{in} =3$ au (corresponding to the dynamical truncation radius for a post-AGB binary system with an orbital separation of $a=1$ au). The contours represent the value of Q at the disc’s outer edge ( $r = r_\mathrm{out}$ ), which corresponds to the minimum Q for surface density slopes $n < 2$ . For $n > 2$ , where Q decreases outwards, the minimum occurs at the inner edge, and we report that value instead.

We find that only sufficiently massive and compact discs are prone to enter the Toomre-unstable regime. Discs with representative post-AGB disc masses ( $\sim0.01$ M $_\odot$ ) are unlikely to becomes unstable under present-day conditions and support planet formation via gravitational instability. This is in contrast with discs around main-sequence stars and white dwarfs (Figure 6) which can become gravitationally unstable with much lower disc masses.

We next examine whether the likely conditions of post-AGB discs at the time of their formation would have made them more susceptible to instability and whether planet formation could have occurred during that earlier phase.

5. Were discs less stable earlier in their lifetimes?

More massive and radially compact discs are more prone to gravitational instability. If the discs we observe today were more massive and compact in the past, they may have been unstable at earlier epochs – even if they appear stable today. Here we use viscous evolution theory of discs to determine their properties in the past and to determine whether, at that time, they may have been gravitationally unstable.

Viscous spreading naturally causes discs to expand and redistribute mass over time (Lynden-Bell & Pringle Reference Lynden-Bell and Pringle1974; Pringle Reference Pringle1981; Hartmann et al. Reference Hartmann, Calvet, Gullbring and D’Alessio1998). The evolution of the disc’s surface density under viscous processes is governed by the diffusion equation:

(17) \begin{equation} \frac{\partial \Sigma}{\partial t} = \frac{3}{r} \frac{\partial}{\partial r} \left[ r^{1/2} \frac{\partial}{\partial r} \left( r^{1/2} \nu \Sigma \right) \right]\end{equation}

where $\nu$ is the kinematic viscosity. Assuming $\nu$ varies with radius as a power law, $\nu \propto r^\gamma$ (index is $\gamma$ ), the radial evolution of the surface density can be solved analytically (Hartmann et al. Reference Hartmann, Calvet, Gullbring and D’Alessio1998):

(18) \begin{align} \Sigma(r,T)=\frac{C}{3\pi \nu_0} \left( \frac{r}{r_\mathrm{ch,0}}\right)^{-\gamma}T ^{-\frac{2.5-\gamma}{2-\gamma}} \exp \left (-\left( \frac{r}{r_\mathrm{ch,0}}\right)^{2-\gamma}\frac{1}{T}\right)\end{align}

where C is a scaling constant, $\nu_0=\nu(r=r_\mathrm{ch,0})$ , $r_\mathrm{ch, 0}$ is the characteristic radius at $t=0$ , and $T=1+t/t_\mathrm{visc}$ is a dimensionless time. The viscous timescale of the disc, $t_\mathrm{visc}$ , is given by:

Figure 6. Toomre parameter Q as a function of disc mass and outer radius, evaluated at $r = r_\mathrm{out}$ , for discs around a Sun-like star (top) and a white dwarf (bottom), assuming a surface density exponent $n = 1$ . Fixed parameters are $R_\mathrm{star}=1\,\textit{R}\odot$ , $T\mathrm{star}=6\,000$ K, $M_\mathrm{cnt}=1\,\textit{M}_\odot$ (top figure for the sun-like star), and $R_\mathrm{star} =0.02\,\textit{R}_\odot$ , $T_\mathrm{star} =30\,000$ K, $M_\mathrm{cnt}=0.5\,\textit{M}_\odot$ (bottom figure for the white dwarf).

(19) \begin{equation} t_\mathrm{visc}=\frac{r_\mathrm{ch, 0}^2}{3(2-\gamma)^2\nu_0}.\end{equation}

Equation (18) is the time-dependent analogue of Equation (13) from Section 2.2.2, showing that surface density profiles spread out from initially more centrally concentrated configurations. The characteristic radius evolves as (Andrews et al. Reference Andrews, Wilner, Hughes, Qi and Dullemond2009)

(20) \begin{equation} r_\mathrm{ch}(T)=r_\mathrm{ch, 0} T ^{\frac{1}{2-\gamma}}\end{equation}

Under the same assumptions, the total disc mass decreases over time due to accretion. Andrews et al. (Reference Andrews, Wilner, Hughes, Qi and Dullemond2009) give the mass evolution as:

(21) \begin{equation} M_\mathrm{disc}(T)=M_\mathrm{disc, 0} T ^{\frac{-1}{2(2-\gamma)}},\end{equation}

where $M_\mathrm{disc, 0}$ is the initial disc mass. In principle, Equations (20) and (21) can be used to infer past disc properties if the age is known. However, disc ages are poorly constrained for both YSOs and post-AGB systems. Stellar ages cannot be used as reliable proxies for disc ages, particularly in YSOs (Andrews et al. Reference Andrews, Wilner, Hughes, Qi and Dullemond2009).

Here, we assume that post-AGB disc evolution is also dominated by viscous processes, as in YSO discs, although the influence of binary companions may modify this picture. The gravitational torques from companions could affect the disc’s evolution, making these equations less valid (Hartmann et al. Reference Hartmann, Calvet, Gullbring and D’Alessio1998). However, since fragmentation is typically expected in outer disc regions, we assume binary effects are second-order for this analysis.

It is unclear how long discs – even in the case of YSOs – remain in the viscous evolution regime, that is, how far back in time the analytical relations of Equations (20) and (21) apply. However, the uncertainty in the viscosity parameter and disc age can be mitigated by expressing the evolution of the disc’s mass and radius in terms of the number of viscous timescales that have elapsed, where each viscous timescale depends on the initial radius via Equation (19).

In Figure 7, we plot the Toomre parameter Q as a function of radius for several YSO discs and a representative post-AGB disc ( $M_\mathrm{disc, now}=0.02$ M $_\odot$ , $r_\mathrm{ch, now}=100$ au, and assuming $\gamma=1$ ), evolved backwards in time for different numbers of elapsed viscous timescales ( $T=t/t_\mathrm{visc}+1$ ). We assume constant inner radii during viscous evolution. For post-AGB discs, the inner radius is set by the dynamical truncation radius, while for YSOs, it is set by the sublimation radius (taken here as the radius where $T_\mathrm{subl}=1\,200$ K). These boundaries depend on stellar properties rather than viscous evolution. Figure 7 shows that even at earlier epochs, some of the discs were likely gravitationally stable – particularly in the cases of PDS 70, HD169142, and the post-AGB example. However, in the case of GQ Lup, HL Tau, and Elias 2-27, they could have been gravitationally unstable in their pasts.

Our analysis includes all YSO systems discussed in Section 3, except AB Aurigae and L1448 IRS3B, whose previously estimated density profiles are not compatible with the form in Equation (13). For post-AGB discs, we assume they follow the viscously evolving form of Equation (13), thus assigning them a characteristic radius. This differs from the simpler power-law model used in earlier sections.

To convert viscous timescales into ages, we adopt the Shakura–Sunyaev (Shakura & Sunyaev Reference Shakura and Sunyaev1973) $\alpha$ -disc prescription for viscosity ( $\nu$ ):

(22) \begin{equation} \nu=\alpha c_s H\end{equation}

where $\alpha$ is the viscosity efficiency parameter, and $c_s$ and H are the sound speed and scale height, respectively. The scale height H of the disc can be found with the following equation:

(23) \begin{equation} H(r)=\frac{c_\mathrm{s}(r)}{\Omega(r)}.\end{equation}

Typical values of $\alpha$ in planetary discs range from 0.0005 to 0.08, based on observed accretion rates (Andrews et al. Reference Andrews, Wilner, Hughes, Qi and Dullemond2009). Combining this with the expressions for $c_s$ and H (Equations 2 and 23) and setting $\eta=1$ , the radial dependence of viscosity becomes:

(24) \begin{equation} \nu(r) =\alpha \left (\frac{k_\mathrm{B}T_\mathrm{sub}}{2m_\mathrm{p}}\right ) \left (\frac{R_\mathrm{ sub}^{0.75}}{(GM_\mathrm{cnt})^{0.5}}\right )r^{0.75}\end{equation}

Outside the sublimation radius, we evaluate $\nu$ at $r = r_\mathrm{ch, 0}$ to determine $\nu_0$ . Here, we assume that the stellar parameters have not evolved significantly over the relevant timescales, so the sublimation radius remains approximately constant. For convenience, we group all constants into a single term K, allowing us to write:

(25) \begin{equation} \nu_0 =Kr_\mathrm{ch, 0}^{0.75},\end{equation}

yielding a viscous timescale of (for $\gamma=1$ ):

(26) \begin{equation} t_\mathrm{visc}=\frac{r_\mathrm{ch, 0}^{1.25}}{3K}\end{equation}

Thus, the disc age can be expressed as $t_\mathrm{visc}\times(T-1)$ .

To calculate these parameters, we adopt representative post-AGB values similar to those used in Section 4 and Figure 3. Disc ages for various assumed viscosity parameters ( $\alpha$ ) are presented in Table 2. Importantly, varying $T=t/t_\mathrm{visc}+1$ alters both the inferred disc age and the initial characteristic radius $r_\mathrm{ch, 0}$ , but we find that the total system age remains relatively stable for each specific value of $\alpha$ .

We find that the disc age does not vary significantly with changes in the scaled time T for a given $\alpha$ , since the viscous timescale itself adjusts accordingly. However, for fixed T, the disc age is inversely proportional to $\alpha$ for $\gamma=1$ . The stellar ages compiled from the literature are included for comparison and should be considered upper limits to the corresponding disc ages. In general, lower $\alpha$ values yield disc ages that exceed the stellar ages – consistent with the findings of Hartmann et al. (Reference Hartmann, Calvet, Gullbring and D’Alessio1998), who argued that $\alpha \approx 0.01$ is typical for YSO discs. For Elias 2-27 specifically, Andrews et al. (Reference Andrews, Wilner, Hughes, Qi and Dullemond2009) derived $\alpha = 0.03$ from accretion rate measurements, which supports our finding that an assumed $\alpha = 0.01$ overestimates the disc age relative to the stellar age.

For post-AGB systems, the inferred ages are consistent with their expected lifetimes of $\sim$ $10^5$ yr (Bujarrabal et al. Reference Bujarrabal2016) for the higher end of viscosity parameters we adopted; $\alpha=0.01$ .

The estimated ages for the central stars of these discs are 1 Myr for HL Tauri (Booth & Ilee Reference Booth and Ilee2020), 1 Myr for Elias 2-27 (Andrews et al. Reference Andrews, Wilner, Hughes, Qi and Dullemond2009), 5.4 Myr for PDS 70 (Keppler et al. Reference Keppler2019), 6 Myr for HD 169142(Hammond et al. Reference Hammond2023), and 2-5 Myr for GQ Lupi (Wu et al. Reference Wu2017).

6. Planet formation in post-common envelope discs

In this section we apply our gravitational instability framework to the post common envelope binary NN Ser, where the presence of two proposed circumbinary planets has been interpreted as evidence for second-generation planet formation in a remnant disc. These planets were first inferred from eclipse timing variations (Beuermann et al. Reference Beuermann2010). Given their proximity to the central binary, it has been argued that they could not have formed during the system’s main sequence evolution and subsequently survived both the giant phase and the common envelope expansion. Instead, it has been proposed that they formed within a circumbinary disc created during the common envelope interaction.

Figure 7. Toomre parameter Q as a function of orbital radius for our YSO sample (excluding AB Aurigae and L1448 IRS3B, for which a characteristic radius is not available) and a representative post-AGB disc (central object parameters similar to those used in Figure 3; a typical value of $\gamma=1$ was assumed due to lack of empirical data). We show Q evaluated at the initial state of each disc for the different assumed viscous timescales.

Table 2. Disc age parameters estimated for each system by assuming different values of the viscosity parameter $\alpha$ . For each case, we list the corresponding viscous timescale ( $t_\mathrm{visc}$ ), number of viscous timescales elapsed ( $t/t_\mathrm{visc}$ ), inferred disc age (t), and the initial characteristic radius ( $r_\mathrm{ ch,0}$ ). For post-AGB discs, we assume $\gamma = 1$ due to a lack of empirical values; for YSOs, $\gamma$ is taken from Table 1. See Section 5 for more details.

The existence of these circumbinary planets, however, has been the subject of debate. While the planets were first inferred from eclipse timing variations (Beuermann et al. Reference Beuermann2010; Beuermann, Dreizler, & Hessman Reference Beuermann, Dreizler and Hessman2013), and later studies concluded the proposed orbits to be dynamically stable (Horner et al. Reference Horner, Wittenmyer, Hinse and Tinney2012; Marsh et al. Reference Marsh2013), the interpretation of the observations as being due to planets is nonetheless regarded with caution. In particular, Marsh et al. (Reference Marsh2013) noted that the planets are only indirectly detected, and that alternative mechanisms – most notably magnetic activity cycles (Applegate’s mechanism Applegate Reference Applegate1992) – can also reproduce the observed eclipse timing variations.

On the assumption of the planets’ existence, Schleicher & Dreizler (Reference Schleicher and Dreizler2014) developed an analytical framework to estimate the properties of a (non observed) circumbinary disc formed by a common envelope ejection. Then applying gravitational instability theory, they demonstrated that the planetary masses and orbital separations inferred for NN Ser are consistent with planet formation in such an environment. Here we connect our gravitational instability framework applied to (observed) post-AGB discs to their analysis of NN Ser, and in so doing we re-examine several aspects of their model.

Figure 8. The Toomre parameter Q as a function of orbital radius for a representative post-AGB disc for varying temperature profile index, p, and the corresponding $\gamma$ , assuming $\gamma=1.5+p$ . The central binary properties are similar to Figure 3, while assuming the surface density profile of a viscous disc (Equation 13). We adopt $r_\mathrm{ch}=100$ au and an inner dynamic truncation radius of $r_\mathrm{in}=3$ au.

By considering common envelope physics, Schleicher & Dreizler (Reference Schleicher and Dreizler2014) derived a disc mass of $M_\mathrm{ disk}=0.146$ M $_\odot$ .Footnote ^e With this disc mass, they went on to assume that the entire disc around NN Ser is marginally gravitationally unstable ( $Q(r)\sim 1$ ) in order to reproduce the observed planet masses.

Here we use our framework to determine whether their disc would be unstable to planet formation. We use their disc mass, their surface density (a power-law index of $n=1$ for the surface density profile – Equation 10), and their central binary mass, $M_\mathrm{core}+M_2=0.646$ M $_\odot$ . However, we use a more realistic temperature profile than used by them, in line with our analysis – Equation (7), where the disc is heated by irradiation from the central source, which in the case of NN Ser is a white dwarf. Finally, we use a disc outer radius of $r_\mathrm{out} =7.50$ au, calculated by conservation of angular momentum before and after disc formation – see Appendix C for details.

Figure 9 shows the radial variation of the Toomre parameter, Q, for NN Ser. The figure shows that even with a more realistic thermal structure, the disc becomes gravitationally unstable beyond $\sim3$ au – precisely where the proposed planets are located. This analysis therefore confirms that the gravitational instability mechanism remains viable under more physically motivated conditions than those selected by Schleicher & Dreizler (Reference Schleicher and Dreizler2014), but only if the disc mass is as high as they propose.

Figure 9. The Toomre parameter, Q, as a function of orbital radius for the progenitor disc of NN Ser, based on the values of Schleicher & Dreizler (Reference Schleicher and Dreizler2014), which are summarised in Table C1. The white dwarf parameters were $T_\mathrm{eff} = 57\,000$ K and $R_\mathrm{star} =0.0189$ R $_\odot$ , and the mass of the central binary was set to be $M_\mathrm{c}+M_2=0.646$ M $_\odot$ . We note that the outer radius of the disc used to reproduce this plot was $r_\mathrm{out} =7.50$ au, which we recomputed using the value we find for angular momentum in Appendix C.

7. Conclusion

In this study, we have carried out a comprehensive analytical investigation into whether post-AGB circumbinary discs can support planet formation via gravitational instability. Our key findings are:

1. 1. Under realistic assumptions of disc mass $\sim0.02$ M $_\odot$ , disc structure, and density and temperature profiles, post-AGB discs are generally gravitationally stable at present. Even systems with the most optimistic parameters – such as IRAS 08544-4431 with a disc mass of $0.1\,\textit{M}_\odot$ – remain above the fragmentation threshold, with Toomre $Q>1$ throughout the disc, largely due to thermal support from the luminous central star. In contrast, discs around lower-luminosity central objects (e.g. white dwarfs or Sun-like stars) can reach the gravitationally unstable regime at lower masses.

2. 2. Using self-similar viscous disc evolution models, in which the surface density profile retains its shape while the overall scale (mass and characteristic radius) evolves due to viscous spreading, we explored whether post-AGB discs may have been unstable earlier in their lifetimes. We found that, even after evolving them backward in time across several viscous timescales, most remain stable, unlike some protoplanetary discs (e.g. HL Tauri, Elias 2-27, GQ Lupi), which likely passed through an unstable phase during earlier epochs. Our analysis also shows that higher viscosity parameters ( $\alpha \sim 10^{-2}$ ) are more consistent with expected post-AGB disc lifetimes.

3. 3. Revisiting the post-common envelope binary NN Ser, we confirm that if the disc parameters proposed by Schleicher & Dreizler (Reference Schleicher and Dreizler2014) are assumed, the disc would become gravitationally unstable beyond $\sim\,3$ au. However recent developments in common envelope theory suggest that fallback discs may be far less massive than previously assumed.

It is important to note that using the Toomre parameter as a diagnostic for gravitational instability, while useful, does not take factors such as vertical structure, turbulence, and binary-driven dynamics into account. Post-AGB discs can differ from YSO discs in these regards as well, hence the results of this study should be treated as a first order analysis.

Overall, our analysis demonstrates that gravitational instability is highly unlikely to account for second-generation planet formation in typical post-AGB discs, as fragmentation would require unrealistically high disc masses ( $\gtrsim 0.4$ M $_\odot$ ), even under compact and cool conditions. The observed signatures of planets or planet-like companions in several post-AGB systems therefore point towards alternative formation channels – most likely a modified version of core accretion operating on shorter timescales in evolved environments. This study provides a first theoretical baseline for interpreting such systems. While future observations will better constrain disc masses, gas-to-dust ratios, temperatures, sizes, and lifetimes – and may uncover additional evidence for planet formation – we will, in parallel, explore additional theoretical models to investigate planet formation pathways in post-main-sequence environments.