The effects of disc self-gravity and radiative cooling on the formation of gaps and spirals by young planets

ABSTRACT

We have carried out 2D hydrodynamical simulations to study the effects of disc self-gravity and radiative cooling on the formation of gaps and spirals. (1) With disc self-gravity included, we find stronger, more tightly wound spirals and deeper gaps in more massive discs. The deeper gaps are due to the larger Angular Momentum Flux (AMF) of the waves excited in more massive discs, as expected from the linear theory. The position of the secondary gap does not change, provided that the disc is not extremely massive (Q ≳ 2). (2) With radiative cooling included, the excited spirals become monotonically more open (less tightly wound) as the disc’s cooling time-scale increases. On the other hand, the amplitude and strength of the spirals decrease when the cooling time increases from a small value to ∼1/Ω, but then the amplitude starts to increase again when the cooling time continues to increase. This indicates that radiative dissipation becomes important for waves with T_cool ∼ 1. Consequently, the induced primary gap is narrower and the secondary gap becomes significantly shallower when the cooling time becomes ∼1/Ω. When the secondary gap is present, the position of it moves to the inner disc from the fast cooling cases to the slow cooling cases. The dependence of gap properties on the cooling time-scale (e.g. in AS 209) provides a new way to constrain the disc optical depth and thus disc surface density.

hydrodynamics, waves, planet–disc interactions, protoplanetary discs

1 INTRODUCTION

Protoplanetary discs from ALMA observations reveal many substructures (e.g. gaps, rings, and spiral arms) in dust continuum emission (Pérez et al. 2016; Andrews et al. 2018; Long et al. 2018). Among them, the most common features are concentric gaps and rings (Huang et al. 2018). There are many interpretations for these features, such as zonal flows (Flock et al. 2015), aggregate sintering (Okuzumi et al. 2016), secular gravitational instabilities (Takahashi & Inutsuka 2014), self-induced dust pile-ups (Gonzalez, Laibe & Maddison 2017), dust growth at snowlines (Zhang, Blake & Bergin 2015; Pinilla et al. 2017), planet–disc interactions (de Juan Ovelar et al. 2013; Dong, Zhu & Whitney 2015), and so on. Recent high-resolution observations for a relatively large sample of discs (Huang et al. 2018; Long et al. 2018; van der Marel et al. 2019) conclude that snowlines cannot consistently match the positions of these gap features, but it may still apply to individual objects.

Of all the exciting scenarios, this paper focuses on the interpretation of planet–disc interactions. Assuming these features are due to the planets, many authors infer the planet mass from the gap properties by comparing observations with hydrodynamical simulations. These properties include gap width and depth in dust emission (Zhang et al. 2018), and sub/super-Keplerian rotational velocity in CO channel maps (Teague et al. 2018; Zhang et al. 2018; Gyeol Yun et al. 2019). Most past works associate each gap with a planet. Recent simulations show more promising results that one single planet can explain up to five gaps in ALMA observations (Bae, Zhu & Hartmann 2017; Guzmán et al. 2018; Zhang et al. 2018). These multiple gaps that are induced by a single planet is due to the presence of multiple spirals, which results from constructive interference of the density waves at different m modes in the inner disc (Bae & Zhu 2018a,b; Miranda & Rafikov 2019a). The one-armed spiral opens the primary (major) gap at the position of the planet, and the secondary spiral opens the secondary gap in the inner disc and so on. The spacing of the gaps is also useful to constrain the gaseous disc scale height h/r since the ratio of the positions between the secondary gap and the primary gap mainly depends on the disc h/r (Dong et al. 2018; Zhang et al. 2018). If we can further constrain h/r using other observables like gap widths and depths, we can use the position of the secondary gap to even constrain the planet mass (Kanagawa et al. 2015, 2016; Bae et al. 2017; Zhang et al. 2018; Gyeol Yun et al. 2019).

However, most simulations that are used to infer the planet mass from the disc substructures neglect the (1) self-gravity from the disc and (2) radiative cooling. For the latter, they employ locally isothermal Equation of State (EoS) instead, which is equivalent to instant cooling. These two simplifications might not be valid to realistic discs. Consequently, the planet masses inferred from observations might be subject to systematic errors. Self-gravity and radiative cooling can play essential roles in some discs.

Several recent works suggest discs might be more massive than previously thought. Booth et al. (2019) measure the HD 163296 disc mass using ¹³C¹⁷O line, and find the gas mass is a factor of 2–6 higher than previously estimated using C¹⁸O. Using the DSHARP opacity (Birnstiel et al. 2018) and optical–thin assumption, Zhang et al. (2018) find that Toomre Q at the gap edges are ≲10 for most of the discs (see table 3 therein). Among the same series of papers, Dullemond et al. (2018) also find that most of the prominent rings are marginally gravitationally stable if the gas-to-dust ratio is 100 (see fig. 7 therein). Furthermore, Zhu et al. (2019) suggests that these discs might be very optically thick, even though dust scattering makes these discs look like optically thin. Without the assumption of the dust opacity, Powell et al. (2019) analyse seven discs by measuring the locations of ‘dust lines’ – the cutoff radii of continuum emission – at different wavelengths, and also find all of those discs have Q ≲ 10, with several approaching unity. For these reasons, disc self-gravity can become quite important for the interactions between planets and discs.

The disc should also have orders of magnitude difference in cooling time at different radii (Zhu et al. 2015; Miranda & Rafikov 2019b). Assuming the Minimum Mass Solar Nebula (MMSN), the cooling is faster at the outer disc and slower at the inner disc – 7 orders of magnitude of the cooling time between 1 au (10⁵/Ω) and 100 au (10⁻²/Ω) (Zhu et al. 2015). The fast cooling resembles the locally isothermal disc, whereas the slow cooling represents the adiabatic disc. However, it is unclear how the disc might look like beyond these two extremes. Radiative cooling with different cooling time-scales at different locations in discs might further change the disc substructures.

In this paper, we study spirals and gaps induced by planets in discs with non-negligible self-gravity and radiative cooling. In Section 2, we lay out the theoretical background of the linear theory and some basic quantities being used throughout the paper. In Section 3, we introduce the simulation setups. In Section 4, we present the results after adding these two physical processes separately, and also explore a situation that includes both self-gravity and radiative cooling. We also compare our simulation results with previous analytical results from Goldreich & Tremaine (1980) and Miranda & Rafikov (2019a). In Section 5, we quantify the change of AMF and discuss the observational implications in linear regime. Then, we take AS 209 disc as a test bed for these two processes, and address some limitations of our simulations. Finally, we conclude this paper in Section 6.

2 THEORETICAL BACKGROUND

Goldreich & Tremaine (1978,1979,1980) develop the linear theory for planet–disc interactions. This theory produces insightful analytical results which can be used to understand numerical simulations on planet migration and gap opening.

In order to study planet–disc interactions in the linear regime, the planet mass should be less than the thermal mass (Goodman & Rafikov 2001)

$$\begin{eqnarray} M_{\mathrm{ th}} = \Big (\frac{h}{r}\Big)_\mathrm{ p}^3 M_* = 1 M_\mathrm{ J} \Big (\frac{(h/r)_\mathrm{ p}}{0.1}\Big)^3 \frac{M_*}{M_\odot } \, , \end{eqnarray}$$

(1)

where (h/r)_p is the disc’s aspect ratio at the planet’s position. M_* is the mass of the central star. At the thermal mass, the gaseous disc scale height is comparable to both the Hill and Bondi radii of the planet, and the distance between the planet and spiral shock forming region is comparable to the planet’s Bondi radius so that spiral shocks begin to affect planet accretion (Béthune & Rafikov 2019). If M_p ≳ M_th, the density wave excited by the planet starts from the non-linear regime, so that the spiral waves immediately become spiral shocks after the wave excitation. Below the thermal mass, the spiral shock region can be well separated from the wave launching area. Most of the simulations in this paper have M_p ≲ M_th, but we also explore several cases in the high mass regime.

Angular momentum transport is one of the key aspects in the accretion theory. In the linear theory, the changing rate of the angular momentum in the disc enclosed within a radius r is equal to the planet’s torque Γ on that region minus the Angular Momentum Flux¹ (AMF) F_J that flows out of the region, that is

$$\begin{eqnarray} \frac{{\mathrm{ d}L}}{\mathrm{ d}t} = {\Gamma }(r) - F_J(r). \end{eqnarray}$$

(2)

2.1 Torque

The first term on the right-hand side of equation (2) is the planet’s torque to the disc region within r,

$$\begin{eqnarray} \Gamma = - \int _{\mathrm{ disc}} \Sigma (\vec{r} \times \vec{F}) \mathrm{ d}f = \int _r \frac{\mathrm{ d}T}{\mathrm{ d}r} \mathrm{ d}r, \end{eqnarray}$$

(3)

as in Kley & Nelson (2012), where

$$\begin{eqnarray} \frac{\mathrm{ d}T}{\mathrm{ d}r}(r) = r\int _0^{2\pi }\Sigma (r,\phi)\frac{\partial \Phi _p}{\partial \phi } \mathrm{ d}\phi \end{eqnarray}$$

(4)

is torque density and Φ_p is the planet’s potential in the disc coordinate,

$$\begin{eqnarray} \Phi _\mathrm{ p} = -\frac{GM_\mathrm{ p}}{\big \lbrace r_\mathrm{ p}^2+r^2-2r_\mathrm{ p}r\mathrm{cos}(\phi -\phi _\mathrm{ p})+s^2\big \rbrace ^{1/2}}, \end{eqnarray}$$

(5)

where r_p, ϕ_p are planet’s position in radial and azimuthal directions and s is the smoothing length used in simulations to avoid the singularity of the planetary potential. The planet feels the opposite torque, which leads to planet migration if the torque from the inner disc is not balanced by the torque from the outer disc.

Goldreich & Tremaine (1980) derive the one-sided torque as

$$\begin{eqnarray} F_{J0} = (M_\mathrm{ p}/M_*)^2h_\mathrm{ p}^{-3}\Sigma _\mathrm{ p} r_\mathrm{ p}^4\Omega _\mathrm{ p}^2 , \end{eqnarray}$$

(6)

where Ω_p is the angular frequency of the planet (all the ‘p’ as a subscript in the paper refers to that quantity evaluated at the position of the planet). This value is commonly used to normalize the torque and AMF.

2.2 AMF

The second term on the right-hand side of equation (2) is the AMF carried by the disturbance (e.g. spiral waves). The angular momentum carried by the spiral waves will be eventually deposited to the background disc to induce gaps in discs. In order to understand the gap properties, such as depth, width, and secondary gap positions, it is necessary to understand the AMF of the spiral waves. With disc self-gravity included in the analysis, AMF has two components, advective AMF, F_A, and gravitational AMF, F_G (Lynden-Bell & Kalnajs 1972; Goldreich & Tremaine 1979). The latter only occurs in self-gravitating discs. The total AMF is, F_J = F_A + F_G.

F_A is the advective transport of angular momentum due to the motion of the fluid and leads to a flow of angular momentum through circumference at radius r. It is also related to the Reynolds stress in turbulence studies. It can be calculated as,

$$\begin{eqnarray} F_\mathrm{ A}(r) = r^2\Sigma (r)\oint u_r(r,\phi)u_{\phi }(r,\phi)\mathrm{ d}\phi , \end{eqnarray}$$

(7)

where Σ is the gaseous disc surface density and u_r and u_ϕ are the velocity perturbations in the r and ϕ directions.

F_G is due to the non-axisymmetric part of the disc potential. The spiral structure produces a spiral gravitational field, which exerts torque and transfers angular momentum from one part of the disc to another. It is non-zero in self-gravitating discs and becomes important when the following Tommre Q parameter is small,

$$\begin{eqnarray} Q = \frac{c_s \Omega }{\pi G\Sigma }, \end{eqnarray}$$

(8)

where c_s is the sound speed and c_s = (h/r)v_ϕ. Given fixed c_s and Ω, a smaller Q means a higher disc surface density. When Q = 1, the disc becomes gravitationally unstable. The value of F_G is the torque exerted to the inner disc,

$$\begin{eqnarray} F_G(r) = \int _{0}^{r} \mathrm{ d}r^{\prime } r^{\prime } \int _{0}^{2\pi } \mathrm{ d}\phi \Sigma \frac{\partial \Phi _{\mathrm{ out}}}{\partial \phi }, \end{eqnarray}$$

(9)

where Φ_out is the disc potential at r (the formula is the same as equation (4), except that Φ_p is replaced with Φ_out). For trailing spiral arms, the angular momentum is transferred from the inner to the outer disc, since the inner disc exerts positive gravitational torque on the outer disc (Lynden-Bell & Kalnajs 1972; Binney & Tremaine 2008).

Note that, in globally isothermal discs, the total AMF is conserved (dF_J/dr = 0) based on the linear theory (Goldreich & Tremaine 1979). However, as Miranda & Rafikov (2019b) point out, this is not the case in locally isothermal discs due to the torque applied on to the wave by the background shear flow (this is also reported in Lin & Papaloizou 2011; Lin 2015). In locally isothermal discs, |$\mathrm{ d}(F_J/c_s^2)/\mathrm{ d}r$| = 0 is conserved.

In the linear theory, all quantities can be further Fourier decomposed into individual m harmonics,

$$\begin{eqnarray} \eta (r,\phi) = \sum _{m=-\infty }^{\infty } \eta _m(r)e^{im(\phi -\phi _p)} , \end{eqnarray}$$

(10)

where η can be u_r, u_ϕ, and Σ, and η_m is a complex number. Putting u_r(r) and u_ϕ(r) in equation (7),

$$\begin{eqnarray} F_{A}(r) &=& 2\pi r^2\Sigma (r) \sum _{m=-\infty }^{\infty }u_{r,m}(r)u_{\phi ,m}^{*}(r) \\ &=& 4\pi r^2\Sigma (r) \sum _{m=0}^{\infty } \Re [u_{r,m}(r)u_{\phi ,m}^{*}(r)]. \end{eqnarray}$$

(11)

The second equality holds because the Fourier transform of real signals is Hermitian. If F_A is expanded with positive integer harmonics (the m = 0 term is always zero),

$$\begin{eqnarray} F_{\mathrm{ A}}(r)= \sum _{m=1}^{\infty } F_{\mathrm{ A},m}(r), \end{eqnarray}$$

(12)

each m component of the advective AMF can be expressed as,

$$\begin{eqnarray} && F_{\mathrm{ A}\mathrm{ },m}(r) = 4\pi r^2\Sigma (r) \Re [u_{r,m}(r)u_{\phi ,m}^{*}(r)] \\ && = 4\pi r^2\Sigma (r) \Big (\Re [u_{r,m}(r)]\Re [u_{\phi ,m}(r)]+\Im [u_{r,m}(r)]\Im [u_{\phi ,m}(r)]\Big), \\ \end{eqnarray}$$

(13)

which is identical to the result using continuous Fourier transform in the shearing sheet geometry (Dong et al. 2011b; Rafikov & Petrovich 2012). Likewise, the total AMF can be expanded as,

$$\begin{eqnarray} F_{J}(r)= \sum _{m=1}^{\infty } F_{J,m}(r). \end{eqnarray}$$

(14)

The sum goes to infinity, but the normalized AMF (⁠|$F_{J,m}/F_{F_{J0}}$|⁠) reaches the maximum at m ∼ (1/2)/(h/r) and becomes ≪1 at several 1/(h/r). Thus, the sum can stop at some m (Goldreich & Tremaine 1980; Artymowicz 1993a,b; Ward 1997), which is known as the ‘torque cutoff’.

2.2.1 Relative importance between F_A and F_G

Here we discuss the relative importance between F_A and F_G. We borrow the WKB solutions in the tightly-wound limit (k ≫ 1/R) obtained by Goldreich & Tremaine (1979). Their analytical solution for the advective AMF is,

$$\begin{eqnarray} F_A = -\frac{\pi m r\Sigma k}{2\pi G\Sigma |k| - k^2c_s^2}\Big (1-\frac{c^2 |k|}{2\pi G\Sigma }\Big)^2\Phi (r)^2 , \end{eqnarray}$$

(15)

where Φ(r) is the amplitude of the perturbed disc gravitational potential in the form of WKB solution. The gravitational AMF is,

$$\begin{eqnarray} F_G = \mathrm{sgn}(k)\frac{mr\Phi ^2(r)}{4G}\, . \end{eqnarray}$$

(16)

Taking the ratio between the F_A and the F_G, and define the critical wavenumber k_crit (Binney & Tremaine 2008) as,

$$\begin{eqnarray} k_{\mathrm{ crit}} = \frac{\Omega ^2}{2\pi G \Sigma } = \frac{Q}{2}\frac{\Omega }{c_s} = \frac{Q}{2h}, \end{eqnarray}$$

(17)

we get,

$$\begin{eqnarray} \frac{F_\mathrm{ A}}{F_\mathrm{ G}} = 2\Big (\frac{Q^2}{4}\frac{|k|}{k_{\mathrm{ crit}}} - 1\Big) = 2\Big (\frac{|k|}{k_{\mathrm{ cross}}} - 1\Big), \end{eqnarray}$$

(18)

where k_cross = 4k_crit/Q² = 2/Qh. When |k| = (3/2)k_cross, advective AMF and gravitational AMF have equal contribution, and we denote that |k| as k_eq. k_eq = (3/2)k_cross = 3/Qh. When |k| < k_eq, F_G has larger contribution, whereas when |k| > k_eq, F_A contributes more. Hence, to compare the relative importance of F_A and F_G, we just need to compare the characteristic values of |k| and k_eq.

Plugging equations (8) and (17) into the dispersion relation,

$$\begin{eqnarray} {\tilde \omega }^2 = \Omega ^2 - 2\pi G\Sigma |k| + c_s^2k^2, \end{eqnarray}$$

(19)

where |${\tilde \omega }(r) = m[\Omega _\mathrm{ p} - \Omega (r)]$| and Ω(r) is the angular velocity of the disc, we get the absolute value of the wavenumber in unit of k_crit,

$$\begin{eqnarray} \frac{|k|}{k_{\mathrm{ crit}}} = \frac{2}{Q} \Bigg \lbrace m^2\Bigg (1-\frac{r^{3/2}}{r_\mathrm{ p}^{3/2}}\Bigg)^2 - \Big (1-\frac{1}{Q^2}\Big) \Bigg \rbrace ^{1/2} + \frac{2}{Q^2}. \end{eqnarray}$$

(20)

Here we choose the positive sign in front of the first term since it corresponds to the short wave that propagates beyond Lindblad resonances in discs (the negative sign solution corresponds to the long wave that is restrained inside Lindblad resonances and outside the forbidden region; see fig. 2 in Lovelace, Jore & Haynes 1997, Fig. 6.14 in Binney & Tremaine 2008 and equation 19 in Goldreich & Tremaine 1979 for details). Note that this equation reduces to the following (Ogilvie & Lubow 2002; Bae & Zhu 2018a) as Q → ∞,

$$\begin{eqnarray} \frac{|k|}{m} = \frac{\Omega }{c_s} \Bigg |\Big (1 - \frac{r^{3/2}}{r_p^{3/2}}\Big)^2 - \frac{1}{m^2}\Bigg |^{1/2}. \end{eqnarray}$$

(21)

With equation (20), we evaluate the ratio of the |k| and k_eq at the effective locations of Lindblad resonances. While both sound pressure and self-gravity shift the positions of the Lindblad resonances (Pierens & Huré 2005), we only account for the contribution from the pressure effect for simplicity. The locations of the resonances become,

$$\begin{eqnarray} r = r_\mathrm{ p} \Bigg \lbrace 1 \pm \frac{\big [1+m^2(h/r)^2\big ]^{1/2}}{m} \Bigg \rbrace ^{2/3}. \end{eqnarray}$$

(22)

At these radii, the ratio reads,

$$\begin{eqnarray} \frac{|k|}{k_{\mathrm{ eq}}} = \frac{Q}{3}\Bigg \lbrace m^2(h/r)^2+\frac{1}{Q^2}\Bigg \rbrace ^{1/2} + \frac{1}{3}. \end{eqnarray}$$

(23)

As Q → ∞, the ratio |k|/k_eq → ∞, which means that F_G has no contribution to the total AMF. However, when Q = 5 and we adopt m = (1/2)/(h/r), |k| ≈ k_eq. Hence, F_G contributes comparable amount of AMF as F_A. This simple calculation shows that as Q becomes smaller, the gravitational AMF becomes increasingly important and will exceed the contribution of the advective AMF at very low Q.

We caution that the detailed analysis on the relative importance of F_A and F_G requires numerically solving the linearized equations and is beyond the scope of this paper. In this paper, since we focus on discs with Q ≳ 5, we will only discuss F_A for simplicity. As shown in Section 5.4, we find that when Q ∼2, the gap edge quickly becomes unstable even though the disc itself is gravitationally stable.

2.2.2 Calculation of the total AMF

Here we summarize the calculation of the total AMF in the linear theory. In regions far from corotation (|r − r_p|/r_p → ∞), advective AMF dominates (F_A ≫ F_G; Goldreich & Tremaine 1980). Since F_J is conserved at different radii, F_J equals F_A far from corotation where WKB solutions (of equation 83 therein) are valid. The analytical solution can be obtained (Goldreich & Tremaine 1978) at the following condition (Goldreich & Tremaine 1980),

$$\begin{eqnarray} m (h/r) \ll \mathrm{min}\Big (\frac{Q^2-1}{3Q^2}, 1\Big). \end{eqnarray}$$

(24)

The resulting AMF is given by,

$$\begin{eqnarray} F_J^{WKB}(m) = \frac{4}{3}\frac{m^2\Sigma }{\Omega ^2}\Big (\frac{GM_p}{R}\Big)^2 \lbrace 2K_0(2/3)+K_1(2/3)\rbrace ^2, \end{eqnarray}$$

(25)

where the last term is around 6.35, and K₀ and K₁ are zeroth and first-order modified Bessel function of the second kind. The value of |$F_J^{WKB}(m)$| is approximately equal to 8.47 m²F_J0. However, the condition (24) is not always satisfied in the disc. The analytical solution deviates from the exact solution when m ≳ 1/(h/r), or disc self-gravity is important (⁠|$m \gtrsim \frac{Q^2-1}{3Q^2}$|⁠). At this regime, Goldreich & Tremaine (1980) numerically integrate the WKB trial solutions for each m harmonic at given Q and obtain F_J(m, Q), expressed in unit of |$F_J^{WKB}(m)$| (see figs 2 and 3 therein and Figs 3 and 4 in this paper). Summing up each m mode, the total AMF is given by,

$$\begin{eqnarray} F_{J} = f(Q) \frac{F_J^{WKB}(m)}{m^2} \approx 8.47 f(Q) F_{J0}, \end{eqnarray}$$

(26)

where F_J0 is expressed in equation (6), and f(Q) is a factor, which is a function of Q. From Goldreich & Tremaine (1980),

$$\begin{eqnarray} f(Q) = \frac{1}{3} \mu ^3_{\mathrm{ max}} = \int _{0}^{\infty }\mathrm{ d}\mu \mu ^2\frac{F_J(m, Q)}{F_J^{WKB}(m)}, \end{eqnarray}$$

(27)

where μ = mc_s/Ωr = m(h/r). The coefficient F_J/F_J0 ≈ 0.93 if Q = ∞ and F_J/F_J0 ≈ 8.6 if Q = 2, which means that more massive discs (with lower Toomre Q values) have higher AMF. We will compare this analytical solution with our simulation results in 4.1.1.

2.3 Pitch angle

The pitch angle β is defined as the angle between the direction of the line tangent to the spiral arm and the azimuthal direction in the disc. If the spiral arms have m-fold rotational symmetry, and k is the radial wavenumber under the WKB approximation,

$$\begin{eqnarray} \beta = \mathrm{ tan}^{-1}\left (\frac{m}{kr}\right). \end{eqnarray}$$

(28)

The dispersion relation for the spiral wave in the tightly wound limit (equation 19) can also be used to derive the pitch angle. In the absence of disc self-gravity, the middle term on the right-hand side of equation (19) becomes zero. Thus, by manipulating equation (21), the pitch angle for the m mode is,

$$\begin{eqnarray} \mathrm{ tan}(\beta) = \frac{m}{kr} = \frac{c_\mathrm{ s}}{\Big \lbrace |\Omega _\mathrm{ p} - \Omega |^2 - \frac{\Omega ^2}{m^2}\Big \rbrace ^{1/2}}. \end{eqnarray}$$

(29)

If m ≫ 1, different m modes have the same pitch angle and interfere with each other (Ogilvie & Lubow 2002), so that

$$\begin{eqnarray} \mathrm{ tan}(\beta) = \frac{m}{kr} = \frac{1}{r} \frac{c_\mathrm{ s}}{|\Omega _\mathrm{ p} - \Omega |}. \end{eqnarray}$$

(30)

When the disc self-gravity is non-negligible, we can rewrite equation (20) to derive the pitch angle

$$\begin{eqnarray} \mathrm{ tan}(\beta) = \frac{1}{r} \frac{c_\mathrm{ s}}{\Big \lbrace |\Omega _\mathrm{ p} - \Omega |^2 - \frac{\Omega ^2}{m^2}\big (1-\frac{1}{Q^2}\big)\Big \rbrace ^{1/2} + \frac{\Omega }{mQ}}\, . \end{eqnarray}$$

(31)

Compared to equation (29), two additional terms are (Ω/mQ)² within the square root and Ω/mQ outside the square root. When Q ≫ 1, equation (31) reduces to equation (29). Equation (31) also indicates that the pitch angle becomes smaller for higher disc masses (with smaller Q). However, this effect is not strong unless Q is very close to unity. For example, when m ≫ 1, equation (31) reduces to equation (30) even in massive discs. Thus, different m modes should still have the same pitch angle and interfere with each other unless Q is close to one. The effect should be stronger in the inner disc than the outer disc since Ω decreases with r. We will confirm this pitch angle calculation using our simulations in Section 4.1.

2.4 Orbital cooling

Miranda & Rafikov (2019b) carry out linear analysis on locally isothermal discs with temperature varying with the radius and run simulations to show that the spiral wave absorbs AMF from the background when it is propagating to hotter regions, resulting in higher AMF and higher density perturbation. The phenomenon is unique to locally isothermal discs. It would not happen even if the adiabatic disc is very close to locally isothermal (γ is very close to 1). They advocate that the inclusion of the energy equation would result in weaker density waves and shallower gaps. Thus, planet masses inferred from previous locally isothermal simulations might be systematically lower than their actual masses.

Besides running locally isothermal simulations, we further study this issue by carrying out adiabatic simulations with a simple orbital cooling prescription. The adiabatic EoS is

$$\begin{eqnarray} P = (\gamma - 1)E, \end{eqnarray}$$

(32)

where E is the internal energy per unit area and γ is the adiabatic index. We adopt γ = 1.4 in this study. Since the energy equation is solved in simulations with the adiabatic EoS, we can prescribe the effect of radiative cooling using the orbital cooling approach:

$$\begin{eqnarray} \frac{\mathrm{ d}E}{\mathrm{ d}t} = -\frac{E - c_v \Sigma T_{\mathrm{ irr}}}{t_{\mathrm{ cool}}}, \end{eqnarray}$$

(33)

where Σ is the disc surface density. T_irr is the disc initial temperature, which is determined by the stellar irradiation in realistic discs. c_v ≡ R/(μ(γ − 1)) is the heat capacity per unit mass, R ≡ k/m_H is the specific gas constant, k is the Boltzmann constant, μ is the mean molecular weight, m_H is the atomic unit mass, and t_cool is the cooling time. It is useful to define the dimensionless cooling time

$$\begin{eqnarray} T_{\mathrm{ cool}} = t_{\mathrm{ cool}}\Omega (r). \end{eqnarray}$$

(34)

Thus, T_cool is t_cool in the unit of orbital time over 2π. This prescription is identical to Zhu et al. (2015) and similar to the β-cooling in Gammie (2001). Small T_cool means fast cooling. Thus, the adiabatic disc with fast cooling should be closer to the isothermal disc, whereas large T_cool means slow cooling, and this disc is closer to the adiabatic disc without cooling. In the rest of the paper, we simply use adiabatic discs to refer to the adiabatic discs without cooling.

To estimate T_cool in realistic discs, we follow Zhu et al. (2015) which use the radiative cooling rate of

$$\begin{eqnarray} \frac{\mathrm{ d}E}{\mathrm{ d}t} = -\frac{16}{3}\sigma (T_{\mathrm{ mid}}^4 - T_{\mathrm{ irr}}^4)\frac{\tau }{1+\tau ^2}, \end{eqnarray}$$

(35)

where σ is the Stefan-Boltzmann constant, τ = (Σ/2)κ_R is the optical depth in the vertical direction, κ_R is the Rosseland mean opacity normalized to the gas surface density assuming dust-to-gas ratio is 1/100, κ_R = κ_R(a_max, T), and T_mid is the mid-plane temperature. Assuming E = c_vΣT_mid and using equations (33) and (35), we can derive

$$\begin{eqnarray} t_{\mathrm{ cool}} = \frac{3\Sigma c_v}{16\sigma (T_{\mathrm{ mid}}^2+T^2_{\mathrm{ irr}})(T_{\mathrm{ mid}}+T_{\mathrm{ irr}})}\frac{1+\tau ^2}{\tau }. \end{eqnarray}$$

(36)

Approximating T_mid = T_irr = T_d, where

$$\begin{eqnarray} T_\mathrm{ d}(r) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\Big (\frac{\phi L_*}{8\pi r^2 \sigma } \Big)^{1/4} & r \le r_{T_\mathrm{ f}} \\ T_{\mathrm{ f}} & r \gt r_{T_{\mathrm{ f}}} \\ \end{array}\right. , \end{eqnarray}$$

(37)

where ϕ = 0.02, representing flaring angle (Dullemond et al. 2018), and L_* is the stellar luminosity. We assume the disc temperature decreases as r and at a given radius it reaches a minimum floor temperature T_f, setting by the background heating processes (e.g. cosmic rays). This radius is

$$\begin{eqnarray} r_{T_\mathrm{ f}} = 39\ \mathrm{au}\ \Big (\frac{T_\mathrm{ f}}{20\ \mathrm{K}}\Big)^{-2} \Big (\frac{L_*}{L_\odot }\Big)^{1/2} \Big (\frac{\phi }{0.02}\Big)^{1/2}. \end{eqnarray}$$

(38)

The dimensionless cooling time becomes,

$$\begin{eqnarray} & T_{\mathrm{ cool}} = 0.015 \Big (\frac{f}{0.01}\Big)^{-1} \Big (\frac{\kappa _R(a_{\mathrm{ max}}, T_d)}{1\ \mathrm{cm^2 g^{-1}}}\Big)^{-1} \Big (\frac{L_*}{L_\odot }\Big)^{-3/4} \Big (\frac{\phi }{0.02}\Big)^{-3/4} \\ & \times \Big (\frac{M_*}{M_\odot }\Big)^{1/2} (1+\tau ^2) \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1 & r \le r_{T_\mathrm{ f}} \\ \big (\frac{r}{r_{T_\mathrm{ f}}}\big)^{-3/2} & r \gt r_{T_\mathrm{ f}} \end{array}\right., \end{eqnarray}$$

(39)

where f is the dust-to-gas mass ratio. For a given protoplanetary disc interior to |$r_{T_\mathrm{ f}}$|⁠, and assuming constant f, T_cool only depends on the Rosseland mean opacity in the optically thin limit, whereas T_cool ∝ τ²/κ_R in the optically thick limit. Within |$r_{T_\mathrm{ f}}$|⁠, it does not explicitly depend on r (the Rosseland mean opacity depends on the temperature which varies with the radius). Note that this optical depth is different from the optical depth in the dust continuum observation, since the Rosseland mean opacity is used for the cooling process. Aside from the optical depth effect, higher opacity, higher stellar radiation, larger flaring angle and smaller stellar mass would make T_cool smaller and vice versa. We also ignore energy diffusion in the radial direction (Goodman & Rafikov 2001), which might be important for very optically thick discs.

Here we provide several typical values of κ_R at different temperatures for references. Using the DSHARP opacity (Birnstiel et al. 2018), the Rosseland mean absorption opacity for the dust population with the maximum dust size a_max = 1 mm and the size distribution n(a) ∝ a^−3.5 is κ_R(20 K) = 0.21 cm²g⁻¹, κ_R(40 K) = 0.48 cm²g⁻¹, and κ_R(125 K) = 1.03 cm²g⁻¹.

For a solar-type star, suppose Σ = 100 g cm⁻², f = 0.01 at 10 au and T = 40 K, T_cool ∼ 20, which is approaching the adiabatic limit. Suppose Σ = 1 g cm⁻², f = 0.01 at 100 au, and T = 20 K, T_cool ∼ 0.02, which is on the locally isothermal limit. Thus, the T_cool estimated from the gap substructures can help constrain the optical depth thus the surface density of the disc. We will use this argument to constrain the surface density at ∼100 au of the AS 209 disc in Section 5.3.

3 METHOD

To test the effects of disc self-gravity and radiative cooling on planet–disc interactions, we run four sets of 2D hydrodynamical simulations. We use FARGO-ADSG (Baruteau & Masset 2008a,b; Baruteau & Zhu 2016) in the first three sets and Athena++ (Stone et al., in preparation) in the fourth set. For all the simulations, we fix the planet in circular orbits and do not consider planet migration. We apply the evanescent boundary condition to all the simulations. We choose not to include the indirect term for separating its effect from the effects studied in this paper. Our simulation setups are motivated and based on Bae & Zhu (2018a) and Bae et al. (2017). We keep all the parameters the same but add disc self-gravity or/and radiative cooling. We denote them as B18 (Set 1) and B17 (Set 2), respectively. We denote simulations with adiabatic EoS and radiative cooling as AD, and simulations with self-gravity as SG. For example, a simulation that is built on Bae et al. (2017) with both orbital cooling and self-gravity is written as B17ADSG. We also run a set of simulations AS209 (Set 3) including radiative cooling and self-gravity with massive planets by adopting the parameters in the AS 209 α varying model in Zhang et al. (2018) ({(h/r)_p, α, M_p} = {0.05, 3 × 10⁻⁴(r/r_p)², 1 M_th}). The Athena++ (Stone et al. 2020, in preparation) simulations (Set 4) have the same setup as B17 and B17AD to verify the results. All the simulations that appear in this paper are summarized in Table 1, and descriptions are detailed as follows.

1. B18 focuses on the linear regime of planet–disc interactions with very low mass planets. This setup is also adopted in Miranda & Rafikov (2019b). The planet mass M_p = 0.01 M_th. The aspect ratio at r_p is (h/r)_p = 0.1. The gas surface density is Σ ∝ r⁻¹. The disc temperature is T(r) ∝ r^−0.5 (h/r ∝ r^0.25) with the locally isothermal EoS. The disc viscosity α is 0 (inviscid). The inner boundary is 0.05 r_p, whereas the outer boundary is 5.0 r_p. The damping regions are inside 0.06 r_p and outside 4.6 r_p. We run five simulations with Q = ∞ (non-self-gravitating), 100, 10, 5, and 2, where Q is evaluated at r_p. The self-gravity smoothing length is 0.3 h. Given the density and temperature profiles, Q is larger at the inner disc and smaller at the outer disc. Q ∝ r^−4/3 and Q at the inner (outer) boundary is 9.5 times (0.30 times) the value at r_p. The resolution in (logarithmically spaced) r and ϕ directions is 4096 × 5580. There are ∼89 grid cells per scale height in the r direction. The planet potential is in the form of equation (5) and the smoothing length s is 0.6 h. The implementation of the Poisson equation solver for self-gravity can be found in Baruteau & Masset (2008b). To study the AMF in simulations with adiabatic EoS and temperature varying with the radius, we also run three simulations (B18AD) with lower resolution (2048 × 2790) for T_cool = 10⁻⁴, 10⁻³, and 0.01 cases.

2. B17 focuses on gap opening by more massive planets. These simulations start with globally isothermal discs. While one subset keeps the temperature constant, the other subset (AD) allows the disc to cool. The aspect ratio (h/r)_p = 0.07, with h/r ∝ r^0.5 so that T is a constant. The disc viscosity α = 5 × 10⁻⁵, and the gas surface density Σ ∝ r⁻¹. The inner boundary is 0.2 r_p and the outer boundary is 2.0 r_p. The damping regions are inside 0.24 r_p and outside 1.6 r_p. The planet smoothing length s is 0.1 h. If the self-gravity is included, the self-gravity smoothing length is 0.3 h. Q ∝ r^−1/2 and Q at the inner (outer) boundary is 2.2 times (0.71 times) the value at r_p. The resolution is 2048 (logarithmically spaced) × 5580 in the r and ϕ directions, so that dr:dϕ ≈ 1:1 for every grid in the domain. There are ∼ 62 grid cells per scale height in the r direction. The planet mass grows with time as M_p = M_{p, f}|$\mathrm{sin^2}(\frac{\pi }{2}\frac{t}{t_{\mathrm{ grow}}})$|⁠, and t_grow is 20 t_p for 0.1 M_th, 60 t_p for 0.3 M_th and 100 t_p for 1 and 3 M_th, where t_p = 2π/Ω is the orbital period of the planet. The planet grows to its full mass at t_grow and stays constant afterwards. We mainly use M_p = 0.3 M_th cases to study the gap opening, but some instability has developed for the Q = 2 disc before its planet grows to the full mass, complicating the analysis. Thus, we add four simulations (M_p = 0.3 M_th, and Q = 100, 10, 5, and 2) with t_grow = 10 t_p. This set of simulations with a shortened planet growing time is almost identical to the previous one at the time-step we choose to analyse, but now the planet grows to its full mass before the instability occurs. Overall, we have Q = ∞, 100, 10, 5, and 2 discs, with M_p = 0.1, 0.3, 1 and 3 M_th.

   This set of simulations also includes dust, represented by 200 000 dust super particles of different sizes. The Stokes numbers (St) of the particles at r_p ranges from 1.57 × 10⁻⁵ to 1.57. The setup for dust particles is identical to Zhang et al. (2018). The Stokes number St for particles (also called particles’ dimensionless stopping time) is

   $$\begin{eqnarray} St=t_{\mathrm{ stop}}\Omega =\frac{\pi s \rho _{\mathrm{ p}}}{2 \Sigma _{\mathrm{ gas}}}=1.57\times 10^{-3}\frac{\rho _{\mathrm{ p}}}{1\, \mathrm{g \, cm^{-3}}}\frac{s}{1 \mathrm{mm}}\frac{100 \mathrm{g\, cm^{-2}}}{\Sigma _{\mathrm{ g}}}\, , \\ \end{eqnarray}$$

   (40)

   where ρ_p is the density of the dust particle, s is the radius of the dust particle, and Σ_g is the gas surface density. We assume ρ_p = 1 g cm⁻³ in our simulations. The Stokes number mentioned above is the Stokes number at the beginning of the simulations, St₀. As the disc evolves, the dust sizes are fixed, but as they can drift in the disc, the Stokes number can change.

   For B17AD and B17ADSG discs, T_cool is constant everywhere in each simulation, but t_cool varies as radius since Ω_k ∝ r^−3/2. We explore T_cool = 0.01, 0.1, 1, 10, and 100, covering fast cooling to slow cooling. Note that for these simulations with cooling, the given h/r and the temperature profiles above are just the initial values. Their values are subject to change as the simulations evolve. We also run a set of B17ADSG discs, which includes both radiative cooling and self-gravity. These discs have Q = 5 at the planet’s position and T_cool = 0.01, 0.1, 1, 10, and 100.

3. AS209 is used to test the effects of radiative cooling in the massive planet regime (∼thermal mass). We run a set of varying α (α = 3 × 10⁻⁴(r/r_p)²) models (AS209AD) with T_cool = 0.01, 0.1, 1, 10, and 100. The (h/r)_p is 0.05 and M_p is 1 M_th or 0.1 M_J if M_* = M_⊙. The resolution in the r and ϕ directions is 1024 × 2790. With lower resolution, we are able to run the simulations longer. We also add a set of simulations that also include self-gravity (AS209ADSG). Q is chosen to be 5 at r_p. The planet grows to its full mass at 20 t_p. Other parameters are the same as (2).

4. Athena++ simulations are used to verify the results. They include an isothermal simulation and a suite of adiabatic simulations, with the same setup as M_p = 0.1 M_th, B17 and B17AD discs in (2), except for a shorter t_grow = 5 t_p. To avoid a numerical artefact due to the small smoothing length, we use s = 0.2 h for T_cool = 100 (twice as large as what is used in the FARGO simulation).

Table 1.

Models.

FARGO-ADSG
Run	Q	Tcool	Mp/Mth	α	(h/r)p	q	Domain (R)	Resolution (R × ϕ)	tgrow/tp
B18	∞	0	0.01	0	0.1	0.5	[0.05, 5]	4096 × 5580	1
B18SG	100	0	0.01	0	0.1	0.5	[0.05, 5]	4096 × 5580	1
B18SG	10	0	0.01	0	0.1	0.5	[0.05, 5]	4096 × 5580	1
B18SG	5	0	0.01	0	0.1	0.5	[0.05, 5]	4096 × 5580	1
B18SG	2	0	0.01	0	0.1	0.5	[0.05, 5]	4096 × 5580	1
B18AD	∞	10−4	0.01	0	0.1	0.5	[0.05, 5]	2048 × 2790	1
B18AD	∞	10−3	0.01	0	0.1	0.5	[0.05, 5]	2048 × 2790	1
B18AD	∞	0.01	0.01	0	0.1	0.5	[0.05, 5]	2048 × 2790	1
B17	∞	0	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
B17SG	100	0	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
B17SG	10	0	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
B17SG	5	0	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
B17SG	2	0	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
B17SG	100	0	0.3	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	10, 60
B17SG	10	0	0.3	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	10, 60
B17SG	5	0	0.3	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	10, 60
B17SG	2	0	0.3	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	10, 60
B17SG	100	0	1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	100
B17SG	10	0	1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	100
B17SG	5	0	1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	100
B17SG	2	0	1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	100
B17SG	100	0	3	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	100
B17SG	10	0	3	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	100
B17SG	5	0	3	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	100
B17SG	2	0	3	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	100
B17AD	∞	0.01	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
B17AD	∞	0.1	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
B17AD	∞	1	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
B17AD	∞	10	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
B17AD	∞	100	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
B17ADSG	5	0.01	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
B17ADSG	5	0.1	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
B17ADSG	5	1	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
B17ADSG	5	10	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
B17ADSG	5	100	0.1	5 × 10−5	0.07	0	[0.2, 2]	2048 × 5580	20
AS209AD	∞	0.01	1	3 × 10−4(r/rp)2	0.05	0	[0.2, 2]	1024 × 2790	20
AS209AD	∞	0.1	1	3 × 10−4(r/rp)2	0.05	0	[0.2, 2]	1024 × 2790	20
AS209AD	∞	1	1	3 × 10−4(r/rp)2	0.05	0	[0.2, 2]	1024 × 2790	20
AS209AD	∞	10	1	3 × 10−4(r/rp)2	0.05	0	[0.2, 2]	1024 × 2790	20
AS209AD	∞	100	1	3 × 10−4(r/rp)2	0.05	0	[0.2, 2]	1024 × 2790	20
AS209ADSG	5	0.01	1	3 × 10−4(r/rp)2	0.05	0	[0.2, 2]	1024 × 2790	20
AS209ADSG	5	0.1	1	3 × 10−4(r/rp)2	0.05	0	[0.2, 2]	1024 × 2790	20
AS209ADSG	5	1	1	3 × 10−4(r/rp)2	0.05	0	[0.2, 2]	1024 × 2790	20
AS209ADSG	5	10	1	3 × 10−4(r/rp)2	0.05	0	[0.2, 2]	1024 × 2790	20
AS209ADSG	5	100	1	3 × 10−4(r/rp)2	0.05	0	[0.2, 2]	1024 × 2790	20
Athena++
B17	∞	0	0.1	5 × 10−5	0.07	0	[0.2, 2]	1024 × 2796	5
B17AD	∞	0.01	0.1	5 × 10−5	0.07	0	[0.2, 2]	1024 × 2796	5
B17AD	∞	0.1	0.1	5 × 10−5	0.07	0	[0.2, 2]	1024 × 2796	5
B17AD	∞	1	0.1	5 × 10−5	0.07	0	[0.2, 2]	1024 × 2796	5
B17AD	∞	10	0.1	5 × 10−5	0.07	0	[0.2, 2]	1024 × 2796	5
B17AD	∞	100	0.1	5 × 10−5	0.07	0	[0.2, 2]	1024 × 2796	5

Note. The q is the exponent in the temperature profile, T(r) ∝ r^−q. In AD models, the (h/r)_p and the temperature profile are just initial values.

4 RESULTS

4.1 Self-gravitating discs

To study the effects of disc self-gravity in the linear regime, we first analyse B18SG discs with M_p = 0.01 M_th. Fig. 1 shows the density perturbations with different disc masses in the polar coordinate. The planet is at ϕ = 180° and r = 1 r_p. The density perturbation is

$$\begin{eqnarray} \delta \Sigma /\Sigma _0 = (\Sigma - \Sigma _0)/\Sigma _0, \end{eqnarray}$$

(41)

where Σ₀ is the gaseous surface density at the initial condition, and Σ is the density at the time-step (t = 10 t_p for this plot) that is analysed. From left to right, the disc mass increases. The Toomre Q values at the r = r_p are 100, 10, 5, and 2, respectively. As the disc mass increases, the spiral perturbations become stronger. This transition is the most evident from the Q = 5 to Q = 2 disc. The density perturbation changes slightly until Q = 5 and has a noticeable increase from the Q = 5 to Q = 2 disc. The underlying AMF that influences the density perturbation is shown in Fig. 2. The AMF increases as Q decreases. The AMF at Q = 100 (blue solid curve) is almost identical to the one in the non-self-gravitating disc (black curve) and identical to the result in fig. 1(a) in Miranda & Rafikov (2019b). As pointed out by Miranda & Rafikov (2019b), we can clearly see that the AMF increases towards the inner disc where the disc is hotter. However, it cannot be captured even with very small cooling time T_cool = 0.01 and 0.001, as shown by the short and long dashed curves. Different from the locally isothermal simulation, the spiral waves cannot pick up the AMF from the background disc when they propagate to hotter regions, even with a relatively short T_cool. Only extremely small cooling time (T_cool = 10⁻⁴, shown by the dotted-dashed curve) can capture this increase of the AMF towards inner disc (the slightly lower values are due to the lower resolution). The jump of AMF is the largest from the Q = 5 to Q = 2 disc, where the AMF value of the latter almost doubles. The increase of the AMF with a smaller Q is consistent with the fact that the density perturbation becomes stronger with a smaller Q.

Figure 1.

The gaseous density perturbations (δΣ/Σ₀) of the inviscid self-gravitating discs (B18SG) with a low mass planet M_p = 0.01 M_th at 10 orbits. The disc masses increase from left to right. At the position of the planet, the Toomre Q parameters are 100, 10, 5, and 2, respectively. Since the density profile goes as r⁻¹, Q is higher at the inner disc and lower at the outer disc. The spirals become stronger and tighter as the mass increases. This is the most evident in the Q = 2 disc.

Figure 2.

The Angular Momentum Flux (AMF) of the B18, B18SG, and B18AD discs calculated from equation (7) and normalized by equation (6). The black curves represent the AMF (solid curve) and torque (dotted curve) of the non-self-gravitating disc, which is almost identical to fig. 1(a) in Miranda & Rafikov (2019b). The AMF increases as the disc mass increases. The AMF of Q = 2 disc is ∼ 3 times that of Q = 100 disc. This results in a stronger density perturbation. The short dashed, long dashed, and dashed-dotted curves in black represent B18AD with T_cool = 0.01, 10⁻³, and 10⁻⁴, and with lower resolution (2048 × 2790). The T_cool = 10⁻⁴ curve is similar to the case without cooling at such resolution.

The spirals become more tightly wound as the disc mass increases. It is the most evident at the inner disc between the Q = 5 and Q = 2 disc. This effect is consistent with the prediction in Section 2 around equation (31). Only when Q becomes very small, the two additional terms become important, which makes the pitch angle smaller. Note that Q is larger than Q_p in the inner disc, given the surface density profile we choose. The changes in pitch angle would be larger if Q is constant across the disc. This trend is also reported in Pohl et al. (2015), where their simulations lie in the non-linear regime with higher planet masses. Hence, the decrease of pitch angle with the disc mass applies to a large range of planet masses.

Yu, Ho & Zhu (2019) measure the pitch angles from 10 protoplanetary discs in near-infrared scattered-light images and three discs in ALMA millimetre dust continuum images. They find a tight relation between the spiral arm pitch angle and the disc mass – more massive discs have smaller pitch angles. While we find that the disc self-gravity cannot lead to such significant changes in pitch angles, spirals in scattered-light images might show more dramatic effects with disc self-gravity which needs to be studied in future. Without the consideration of the self-gravity, the different pitch angles between scattered light and millimetre continuum images might be due to the vertical temperature gradient (Juhász & G. 2018; Rosotti et al. 2019). The near-infrared scattered-light observations probe higher atmosphere with higher temperature, thus the spirals have larger pitch angles than those in the mid-plane probed by ALMA.

4.1.1 Comparison to the linear theory

In this subsection, we provide a quantitative comparison with the linear theory. In Section 2, we briefly introduce the AMF in discs with and without self-gravity. The quantity we want to compare is |$F_J(m, Q)/F_J^{WKB}(m)$| (ratio of the numerically integrated m harmonics of the AMF and the m mode AMF solved analytically; see equation 27), which has been given in fig. 2 of Goldreich & Tremaine (1980). Similar comparisons have been performed before in non-self-gravitating shearing sheet simulations (Dong et al. 2011b). With our simulations, we measure the AMF at several h away from a low-mass planet (M_p ≪ M_th) in the inner disc. Note that the spirals will shock at a certain distance away from the planet. This distance is l_sh (Goodman & Rafikov 2001), given by

$$\begin{eqnarray} l_{\mathrm{ sh}} \approx 0.8 \Big (\frac{\gamma +1}{12/5}\frac{M_\mathrm{ p}}{M_{\mathrm{ th}}}\Big)^{-2/5} h. \end{eqnarray}$$

(42)

The linear theory of wave propagation fails for |r − r_p| > l_sh, so the position that we measure the AMF cannot be too far from the planet. The position should also be chosen closer to the planet as the planet mass increases, as l_sh becomes shorter with higher mass planets.

We start the comparison with a very low planet mass, M_p = 0.01 M_th in B18SG discs. Fig. 3 shows the AMF measured at 5h, where we find the values match the linear theory the best. The blue curve shows the ratio of the non-self-gravitating disc, and the orange curve shows the ratio of the Q = 2 disc. The ratios calculated from Goldreich & Tremaine (1980) are shown in dashed curves (the ratio of Q = 2, |$F_{J}(m, 2)/F_{J}^{WKB}(m)$| is shown above and that of Q = ∞, |$F_{J}(m, \infty)/F_{J}^{WKB}(m)$| is shown below). The Q = ∞ curve matches the linear calculation relatively well, especially at μ ≲ 1. As for the differences, one of the reasons is that the temperature is assumed to be constant throughout the disc in the linear calculation, whereas T ∝ r^−0.5 in the simulations. |$F_J(m, Q)/F_J^{WKB}(m)$| of the self-gravitating disc is higher than that of Q = ∞ disc, which is consistent with the linear theory.

Figure 3.

Comparison between |$F_J(m, Q)/F_J^{WKB}(m)$| measured from simulations (M_p = 0.01 M_th, B18SG) and calculated from the linear theory (Goldreich & Tremaine 1980). This ratio is measured at 5h (at 0.5 r_p) inside of the planet and is plotted on y-axis. The x-axis is |$\mu = \frac{m c_s}{\Omega R}$| (or m(h/r), where (h/r)_p = 0.1). The simulation results are shown in blue (Q = ∞) and orange (Q = 2) solid curves, whereas the analytical results are shown in dashed curves (Q = 2 curve is the higher).

To make a more proper comparison with the analytical theory, we also present B17SG discs that are globally isothermal with M_p = 0.1 M_th. |$F_{J}(m, Q)/F_{J}^{WKB}(m)$| at different distances from the planet are shown in Fig. 4 (Top panel: Q = 100 disc, Bottom panel: Q = 2 disc). Similar to the results of isothermal shearing sheet simulations in Dong et al. (2011b), we also find a close agreement with the linear theory when Q → ∞. However, the values in Q = 2 disc are still different from the linear theory expectation. They are lower by one order of magnitude when μ ≲ 1 and higher than the linear theory prediction at μ ∼ 6. The discrepancy in the high m (or high μ) is possibly due to the limited resolution in the azimuthal direction in simulations. We need very high resolutions to resolve the high m modes properly with a much smaller self-gravity smoothing kernel. As for the low m modes, it is actually not surprising to see discrepancy, since we only measure the advective AMF, and ignore the gravitational AMF (see Section 2.2). In Q = 2 discs, F_G should contribute a significant amount of AMF. It should be more important when Q is even smaller (close to 1). To match the analytical result of total AMF, one has to include F_G. On the other hand, we do see the increases of advective AMF as the Q becomes smaller. This indicates that as the disc becomes more massive, both F_A and FG increase to contribute for the increase of total AMF, F_J. Being part of the AMF, F_A alone is enough to explain the change of the density perturbations, i.e. the increases of F_A leads to the increases of the density perturbations. It is not necessary to add F_G when we want to explain the density perturbations qualitatively.

Figure 4.

Similar to Fig. 3, but for B17SG, M_p = 0.1 M_th. This ratio is measured at 2h (blue), 3h (orange), 4h (green), and 5h (red) inside of the planet. The upper panel shows the results in Q = 100 disc, whereas the lower panel shows the results in Q = 2 disc. The dashed curves are the same as in Fig. 3.

4.1.2 Gap opening

In the linear theory, if the planet mass M_p is large enough (close to or higher than M_th), the planet can open gaps in the disc. However, the gap opening mass can be much lower, since the non-linear effects such as shock can occur even if the system starts in a linear regime (Goodman & Rafikov 2001; Rafikov 2002). The study of gap opening is very important since gaps are observable in both gaseous and dusty discs. They have been found ubiquitously in recent high-resolution ALMA observations (Huang et al. 2018; Long et al. 2018). Assuming these gaps are due to the planet, one might infer the planet mass from the gap width and depth (e.g. Zhang et al. 2018). Now we use higher mass planets to demonstrate the effect of the self-gravity on gap opening.

We first show the results of B17SG discs with M_p = 0.3 M_th, where the evolution starts in linear regime, marginally opening gap for non-self-gravitating disc at longer time evolution. Figs 5 and 6 show the density perturbations in 2D and 1D at t = 20 t_p (we use the models whose t_grow = 10 t_p here). As the Q decreases, the gap becomes deeper. Pohl et al. (2015), Li & Li (2016) also find this effect for higher mass planets. The transition is the most evident from Q = 5 to Q = 2 discs, where the gap depth becomes a factor of 3 stronger at the major gap where the planet is located. The major gap is separated by the horseshoe region, thus resembles ‘w’ shape. The gap width changes only slightly comparing to the gap depth (see Section 5.2 for details). The secondary gap is around 0.5 to 0.6 r_p, the location of which does not change as the Q decreases.

Figure 5.

The gas density perturbations (δΣ/Σ₀) of the isothermal self-gravitating discs B17SG with a low mass planet M_p = 0.3 M_th at 20 orbits (with t_grow = 10 t_p). Compared to Fig. 1, higher planet mass can open deeper gap and even the secondary gap (see Fig. 6). The depth of the primary gap becomes deeper with the decreases of Q. This is the most evident in Q = 2 disc.

Figure 6.

The azimuthally averaged density perturbations of Fig. 5. The blue, orange, green, and red curves represent Q = 100, 10, 5, and 2 discs. As the disc mass increases, the perturbations become stronger and gaps become deeper, but the position of the secondary gap does not change significantly.

To understand the density perturbations, we plot the AMF of each simulation in Fig. 7. The left-hand panel shows M_p = 0.1 M_th discs and right-hand panel shows M_p = 0.3 M_th discs. Given the same Q, the normalized AMF decreases as the planet mass increases, which is due to the increasing non-linear shock dissipation. This is consistent with the result of non-self-gravitating discs in Miranda & Rafikov (2019b). The changes of the AMF at each Q has the same trend as the changes in density perturbations. From Q = 5 to Q = 2, the AMF doubles. We will quantify this trend in Section 5.1. Since the normalization factor F_J0 ∝ Mp², the absolute value of AMF is still higher in the discs with higher planet mass.

Figure 7.

The AMF of the isothermal discs at 20 orbits with M_p = 0.1 and 0.3 M_th. Given the same disc mass, as the planet mass increases, the normalized AMF becomes lower at any disc with different Q, as shock forms earlier and closer to the planet with the increases of the planet mass.

To study how disc self-gravity affects the dust continuum observations, we evolve simulations for a longer time-scale (500 t_p) and plot the density perturbations of both gas and dust with different Stokes numbers in Fig. 8. From top to bottom, it shows the δΣ/Σ₀ of gas, and dust with St = 0.001 (represented by St = [0.0005, 0.002]), 0.01 (represented by St = [0.005, 0.02]) and 0.1 (represented by St = [0.05, 0.2]). Notice that the St₀ is the Stokes number at the beginning of the simulations. As the simulations evolve, the particle sizes are fixed, so the dust grain might have different St after it drifts. From left to right, the planet masses are M_p = 0.1, 0.3, and 1.0 M_th. The blue, orange, and green curves represent Q = 100, 10, and 5 discs. As expected, the perturbations in the dust are stronger than those in gas, and dust with higher Stokes number has larger perturbations. This is because dust particles with higher St drift faster when St ≲ 1. For M_p = 0.1 M_th, St = 0.1 disc, the dust particles have almost drifted to the inner disc. Similar to the trend in the gas, the gap becomes deeper, and the ring becomes higher as the Q becomes smaller, which are consistent with Li & Li (2016). The positions of the secondary gaps do not change. The overall shapes of the profiles do not change, except the gaps outside the planet become weaker for M_p = 0.3 and 1.0 M_th discs at St₀ = 0.001. The positions of them also shift outside. The region beyond r = 1.6 r_p should not be trusted as it is in the wave damping zone. For small Q and large M_p, some curves for discs with Q = 5 are missing since the discs have already become unstable at 500 t_p. We will discuss this phenomenon in Section 5.4.

Figure 8.

The density perturbations of the gas and dusts with different Stokes number (St = 0.001, 0.01, and 0.1) at t= 500 t_p in globally isothermal disc with M_p = 0.1, 0.3, and 1 M_th (from left to right). From top to bottom panels, there are density perturbations of the gas, St = 0.001, St = 0.01 and St = 0.1 sized dusts. The blue, orange, and green curves represent Q = 100, 10, and 5 discs. Some curves for discs with Q = 5 are missing since the disc becomes unstable at this time-steps in the simulations. Overall, the position of the secondary gap does not change with the increase of disc masses. The gap becomes deeper for lower Q values.

4.2 Adiabatic discs with radiative cooling

In this subsection, we study the effects of radiative cooling on planet–disc interactions, focusing on B17AD. The simulations presented in this section are globally isothermal, with q = 0. The absorption of AMF from the background occurs in locally isothermal disc (in Figs 1 and 2), but not in globally isothermal discs here. From left to right, Fig. 9 shows the density perturbations for the M_p = 0.1 M_th isothermal disc, and adiabatic discs with cooling time T_cool = 0.01, 0.1, 1, 10, and 100. The top panels are results from FARGO-ADSG, whereas the bottom panels are results from Athena++. They are almost identical, except for the small difference for the corotation features in the T_cool = 100 cases. We analyse the simulations at t = 40 t_p for FARGO and 20 t_p forAthena++ due to a shorter t_grow adopted in Athena++ simulations. At these times, we find that all the discs are stable and have reached quasi-steady states. Discs with shorter dimensionless cooling times T_cool are more similar to isothermal discs, whereas discs with larger T_cool resemble adiabatic discs. The primary spirals are marked as ‘P’, the secondary spirals are marked as ‘S’, and the tertiary spirals are marked as ‘T’ on the figure.

Figure 9.

The density perturbations of the adiabatic simulations B17AD at 40 orbits (FARGO) and 20 orbits (Athena++) with M_p = 0.1 M_th. From left to right, the discs represent fast to slow cooling. The isothermal dick is shown in the leftmost panel, then the dimensionless disc cooling parameters T_cool are 0.01, 0.1, 1, 10, and 100, respectively. The upper panels show results from FARGO-ADSG and the lower panels show results from Athena++. The primary spirals are marked as ‘P’, the secondary spirals are marked as ‘S’, and the tertiary spirals are marked as ‘T’.

From the isothermal to the adiabatic limit (increasing T_cool from the left to the right-hand panels), the spirals become weaker as the cooling time-scale becomes longer until T_cool = 1 and then become stronger again as T_cool continues to increase. The secondary and tertiary spirals also become very weak at T_cool = 1 and stronger at T_cool = 10. Nevertheless, the secondary and tertiary spirals are still present at T_cool = 1 as indicated by the arrows.

The spirals become more open as the cooling time increases. This is due to the increases of the sound speed from c_s = (RT/μ)^1/2 in the isothermal limit to c_s = (γRT/μ)^1/2 in the adiabatic limit. In other words, the effective γ in equation (30) increases from 1 to 1.4 as T_cool becomes larger.

The change of the amplitude and position of the spiral arms with respect to the cooling time can be seen more clearly in 1D plots. Fig. 10 (left-hand panels) shows the density perturbations versus azimuthal angle (ϕ) cut at r = 0.45 r_p and r = 0.32 r_p. The density perturbations due to the primary spiral and the secondary spiral are marked as ‘P’ and ‘S’ on the figure. The one due to the tertiary spiral is also marked as ‘T’ in the r = 0.32 r_p plot. The positions of the primary, secondary, and tertiary spirals in T_cool = 0.01 discs are similar to those in the isothermal disc. At r = 0.45 r_p, the amplitude of the spiral in the T_cool = 0.01 disc is slightly lower than that in the isothermal disc

Figure 10.

The density perturbations along r = 0.45 r_p (upper panels) and r = 0.32 r_p (lower panels) for B17AD (left-hand panels) and B17ADSG (right-hand panels) discs at t = 40 t_p. The positions of the primary, secondary, and tertiary spirals for small T_cool are marked as ‘P’, ‘S’, and ‘T’ on the left-hand panels.

There is a significant change of the spiral’s amplitude when T_cool increases from ≲0.01 to 0.1. The spiral’s amplitude becomes much smaller at T_cool = 0.1 with a slight shift of the primary spiral position. The tertiary spiral almost disappears at r = 0.32 r_p, but we can still see its presence with a amplitude of |$\sim 1{{\ \rm per\ cent}}$|⁠. But the biggest change of the spiral’s position occurs when T_cool increases from 0.1 to 1. At r = 0.45 r_p, the primary spiral shifts ∼45° between T_cool = 0.1 to T_cool = 1 cases, although the amplitude does not change much. The tertiary spiral at r = 0.32 r_p is now the hump from 160° to 250° with a amplitude of |$\sim 0.5{{\ \rm per\ cent}}$|⁠. At T_cool = 10, the amplitude of the spiral becomes much higher. The positions are very similar to those in T_cool = 1. The tertiary spiral also becomes stronger. But the spirals are still weaker than those in the isothermal disc. From T_cool = 10 to T_cool = 100 discs, the perturbations seem to converge to the limit of an adiabatic disc. Overall, from isothermal discs (T_cool ≲ 0.01) to adiabatic discs (T_cool ≳ 10) the azimuthal angles of the spirals shift to smaller values by ∼50° at r = 0.45 r_p. This is consistent with the finding in 2D plots that the spirals become less tightly wound as T_cool increases.²

Since spirals are directly related to gap opening, we want to study how the properties of the induced gaps are related to radiative cooling. The top panel of Fig. 11 shows the azimuthally averaged density perturbations of discs with different cooling times. T_cool = 0.01 disc has a similar profile to the isothermal disc, but with smaller amplitudes. The secondary gap at T_cool = 0.1 and 1 disappears (the disappearance of the gap in gas does not necessarily lead to that in dust; see Section 5.3). It reappears at T_cool = 10, but the position changes from 0.5 r_p to 0.4 r_p and the amplitude becomes smaller. Since spirals carve gaps in discs, the inward shift of the secondary gap position is due to the more open secondary spiral for higher T_cool. On the other hand, the non-monotonically transition – the trend of high, low, and high amplitudes of the secondary gap depth from the fast cooling to the slow cooling cases – is due to weaker spirals at intermediate T_cool. Overall, for the intermediate massive planet (M_p ∼ 0.1 M_th), discs with T_cool ≲ 0.01 can be treated as isothermal discs and discs with T_cool ≳ 10 can be seen as purely adiabatic discs, whereas the discs undergo a sharp transition between these limits, characterized by low amplitude spirals and shallow secondary gaps. Note that this T_cool ≲ 0.01 criterion is more stringent with a low mass planet in the linear regime. As shown in Fig. 2, only when T_cool = 10⁻⁴ the AMF is close to the isothermal case with M_p = 0.01 M_th. Both the temperature gradient (T ∝ r^−0.5) and the linear damping contribute to the AMF profile, but the linear damping is the dominant effect here.

Figure 11.

Top panel: The azimuthally averaged density perturbations in Fig. 9, upper panels. Bottom panel: The normalized AMF (solid curves) and torque (dashed curves) for different cooling times, T_cool. The black curves represent isothermal disc, whereas blue, orange, green, red, and purple curves represent adiabatic discs with T_cool = 0.01, 0.1, 1, 10, and 100, respectively. The changes of the AMF distribution should be responsible for the change of the density perturbations.

Since such conditions with intermediate cooling times are very likely to occur in realistic protoplanetary discs (see Section 2), we are interested in whether these much shallower secondary gaps can still trap enough dust particles to be observable. Examining the previous models in Zhang et al. (2018), we find the current set of simulations has very similar setup to the h/r = 0.07, α = 10⁻⁴, M_p/M_* = 11M_⊕/M_⊙ model in Zhang et al. (2018) (the upper-left plots in 2D figures therein), despite the α is lower here. However, the secondary gap cannot be found in any dust configurations in that model, even though those simulations assume the locally isothermal EoS. Thus, we conclude that the models being analysed in this subsection with such low planet mass is not suitable for testing the observability of the secondary gap. To that end, we run simulations with higher planet mass, 1 M_th, close to the setup of the AS 209 simulation, α varying model (Zhang et al. 2018). We direct the discussion to Section 5.3.

We also plot the AMF (in solid curves) and the integrated torque (in dashed curves) of these discs in the lower panel of Fig. 11. The torque is integrated from the planet’s position towards either side of the disc. The AMF decreases monotonously close to the planet (|r − r_p| ≲ 0.1 r_p) due to the torque cutoff. The AMF far away from the planet (e.g. at r = 0.3 r_p and 1.5 r_p) decreases to the lowest values at T_cool = 0.1 and 1 and increases back to the isothermal values at T_cool = 10 and 100. On the other hand, from the fast cooling (and isothermal) to the slow cooling cases, the torque first increases and reaches the maximum at T_cool = 1 and then decreases and reaches the minimum at T_cool = 100. The AMF and torque at several h away from the planet (e.g. at r = 0.85 r_p or 1.2 r_p) have similar values in both the very small (e.g. T_cool = 0.01) and very large (e.g. T_cool = 100) T_cool cases. However, the torque there is much higher than the AMF at T_cool = 0.1, 1, and 10. This indicates that, for intermediate cooling cases, the waves fail to launch or are damped right after they are launched, indicating radiative cooling is important for these cases. Note that the normalization (characteristic value, F_J0) should have varied with different cooling times, since |$F_{J0} \propto h_p^{-3}$|⁠. In adiabatic EoS, the sound speed, |$c_s = \big (\gamma \frac{\mathrm{ d}p}{\mathrm{ d}\rho }\big)^{1/2}$| and the scale height h ∝ c_s ∝ γ^1/2. Thus, F_J ∝ γ^−3/2. The effective normalization is 0.6 times the F_J0 in equation (6) at the slow cooling limit when γ = 1.4. This value should be in between 1 and 0.6 in transition from small to large T_cool. Nevertheless, we just use F_J0 as equation (6) and compare the absolute values of the AMF. We will have more detailed discussions in Section 5.1.

4.3 Disks with both self-gravity and radiative cooling

We also explore a situation that both self-gravity and radiative cooling are included, even though we do not seek to carry out a parameter space study. The parameters we choose are Q = 5, B17ADSG discs with various T_cool. The density perturbations at two radii are shown in Fig. 10’s right-hand panels. The profiles of the spirals are very similar to the non-self-gravitating ones on the left-hand panels. All the trends found in Section 4.2 still apply. At the same cooling time, the perturbation amplitudes become stronger and the spirals become tighter if Q becomes smaller, which are consistent with the results in Section 4.1. This means that when both of these two physical processes (disk self-gravity and radiative cooling) are operating, they play similar roles as those that they play individually, at least for low mass planets. From the magnitudes of the changes due to these two processes, we can conclude that when the disc is not too massive (Q ≳ 5), the effects of the self-gravity on the discs are less important than those of the cooling.

5 DISCUSSION

5.1 AMF versus Q and T_cool

To quantify the relationships between AMF versus Q and AMF versus T_cool, we measure the AMF at given positions in Figs 7 and 11 for B17SG and B17AD discs with M_p = 0.1 M_th. We plot the AMF as a function of Q (the left-hand panel) and T_cool (the right-hand panel) in Fig. 12. The AMF is measured at 0.87 r_p for B17SG simulations, where it reaches the maximum at the inner disc. The AMF increases as the disc becomes more massive (with a smaller Q). This increase becomes quite fast when Q is approaching 2. The advective AMF in the Q = 2 simulation is almost five times the AMF in the Q = 100 simulation, which is roughly consistent with the analytical expectation (equation 26).

Figure 12.

AMF versus Q (left-hand panel) and AMF versus T_cool (right-hand panel). The AMF is measured at r = 0.87 r_p for the SG runs and at r = 0.6 r_p (marked as ‘X’s) and 0.9 r_p (marked as dots) for the AD runs.

As for the radiative cooling, we measure the AMF at 0.9 r_p and 0.6 r_p. At the position very close to the planet (e.g. 0.9 r_p), the AMF decreases as the cooling time increases. As mentioned in Section 4.2, this is most likely due to the transition of the effective γ (that enters the sound speed) from 1 to 1.4. The ratio between the AMF at T_cool = 100 and that at T_cool = 0.01 is ∼0.6, which is very close to the expected change due to our definition of F_J0. However, if the AMF is measured farther away from the planet (e.g. 0.6 r_p), it first decreases and then increases as T_cool increases, reaching the minima at 0.1 and 1 T_cool. The decrease of the AMF with the distance away from the planet is due to the wave damping as the wave propagates. When the wave propagates to 0.6 r_p, the damping is the strongest at T_cool ∼ 0.1 − 1, implying that radiative cooling plays an important role on the wave damping/dissipation when T_cool ≲ 1. We are expecting a stronger wave damping with a smaller T_cool. However, when T_cool becomes very small and approaches the isothermal limit, the temperature of the wave becomes the background disc temperature, and shock damping in this case behaves similar to the shock damping in the adiabatic limit (Goodman & Rafikov 2001; Dong, Rafikov & Stone 2011a). At the position that the secondary spiral should form, the wave becomes already too weak to perturb the gas for the T_cool = 0.1 − 1 cases. Thus, the secondary gap can even disappear at these intermediate cooling times.

5.2 Implications for observations

To quantify the properties of the primary gap, we measure the gap width and depth in the gas density perturbation maps for B17SG and B17AD, M_p = 0.1 M_th, using the definitions in Zhang et al. (2018). Different from Zhang et al. (2018), we only measure the gaps at ∼50 orbits and from the gas density perturbation profiles, because we do not run the simulations for very long time and the planet mass is still too small to open significant gaps in gas and dust. However, these measurements suffice to demonstrate the trend.

In the upper panels of Fig. 13, the gap width–Q relation is shown on the left-hand panel, whereas the depth–Q relation is shown on the right-hand panel. Both the x and y axes are in logarithmic scales. Since the exponents in the gap depth relations are ∼5 times larger than that of gap width relations (see Table 1 and 2 in Zhang et al. 2018 for details), the plotted y-axis range for the depth (δ − 1) is 5 times of the range for the width (Δ). In this way the relative magnitude read on both y-axes is roughly proportional to the planet mass ratio. As Q changes to lower values, the gap depth has more significant changes than the width. The planet mass inferred from the gap depth without considering disc self-gravity will be overestimated. Adopting common exponents 0.25 and 1.25 for the width–M_p and depth–M_p relations (more precisely, width–K′ and depth–K in Zhang et al. 2018), the ratio of the inferred masses between Q = 2 and Q = 100 discs is 1.5 using the width–M_p relation, and 6 using the depth–M_p relation. Thus, if one wants to infer the planet mass, the gap width is a better observable than the gap depth since it is less sensitive to the disc mass.

Figure 13.

Gap width/depth versus Q (top) and Gap width/depth versus T_cool (bottom). Δ and δ - 1 are proportional to M_p. The range for the gap (depth - 1) is five times than that in the width. Under this factor, the y-axes on the left-hand and right-hand panels have approximately the same scaling to M_p. Thus, the y-axes of the left-hand and right-hand panels can be directly taken as the planet mass in an arbitrary unit.

The lower panels of Fig. 13 show the width and depth as functions of the T_cool. The gap width and depth are not very well defined due to the irregular gap shape. We try to neglect the horseshoe region around r_p when we measure the gap width. It is better to use the gap depth as a proxy for the planet mass. If one uses the gap width to measure the planet mass in locally isothermal simulations, the planet mass will be underestimated if T_cool ≳ 0.01.

Note that these trends found in the low planet mass regime do not necessarily apply to the high planet mass regime, as we will discuss in the next subsection, where we find much less change in the gap shape for various cooling times.

5.3 AS 209: cooling in the high planet mass regime

We have two goals in this subsection. One is to explore the effects of radiative cooling on the gap properties when the planet mass is relatively large reaching the thermal mass. Meanwhile, we decide to choose a model that resembles a realistic disc so that we are able to constrain its surface density via estimating its T_cool. Of all the DSHARP discs, AS 209 features many intricate gaps and rings (Guzmán et al. 2018). Using simple (without radiative cooling and self-gravity) hydrodynamical simulations, Zhang et al. (2018) find an excellent match to the AS 209 disc with a model whose (h/r)_p = 0.05, α = 3 × 10⁻⁴(r/r_p)², and M_p/M_* = 10⁻⁴, or 1 M_th.³ The dust size distribution is a power law n(a) ∝ a^−3.5, a_max = 0.68 mm, and the gas surface density Σ_{g, 0} = 6.4 g cm⁻¹. This model can explain up to five gaps in the disc, not only their locations, but also the intensity. Since we have found that self-gravity and radiative cooling have the potential to change the gap shape and the secondary gap position in the low planet mass regime, we would like to explore how the planet mass for AS 209 will change considering these two physical processes.

Fig. 14 shows the azimuthally averaged density perturbation and brightness temperature of the AS209 models with different cooling times (T_cool = 0.01, 0.1, 1, and 100) and different disc masses (Q = ∞ and 5) at 400 t_p. We choose not to show T_cool = 10 case to avoid crowdedness of the figure, but its feature is similar to the T_cool = 100 case with shallower gaps in dust. The left-hand panels show the Q = ∞ cases, whereas the right-hand panels show the Q = 5 cases. The blue, orange, green, and red curves represent T_cool = 0.01, 1, 10, and 100 cases, respectively. From the top to bottom panels, it shows the density perturbations of the gas, the dust with St_max = 3 × 10⁻³, St_max = 3 × 10⁻², and the brightness temperature of the St_max = 3 × 10⁻² models. We assume that the temperature in the disc is proportional to r^−0.5 in the last row and the procedure to obtain the optical depth is detailed in Zhang et al. (2018). As the T_cool becomes larger, the perturbation at the major gap (where the planet locates) becomes shallower in gas. However, this is not the case in dust. The gap width and depth are almost the same across different cooling times. This result is different from that in the low planet mass regime, where the gap shapes are quite different across different cooling times. For dust profiles, the horseshoe region has the smallest density for T_cool = 1. The secondary gap disappears for T_cool = 1 in gas, but it is still noticeable in dust for St_max = 3 × 10⁻², with a lower amplitude than T_cool = 0.01, 0.1, and 100 cases. This indicates that, due to the gaseous bump at the inner edge of the primary gap, dust particles can still be trapped there, forming the secondary gap in dust even if the secondary gap disappears in the T_cool = 1 case. The location of the secondary gap shifts from ∼0.65 r_p to ∼0.6 r_p when T_cool increases. The tertiary gap is also present in T_cool = 0.01 and 100 cases. It moves inwards from ∼0.43 r_p to ∼0.37 r_p. Although the tertiary gap is not present in the T_cool = 1 case, a change of slope can been seen at ∼0.4 r_p (this is more obvious in Q = 5 case). The Q = 5 discs are similar to the Q = ∞ discs but with slightly stronger perturbations. The differences between different T_cool curves are slightly smaller in Q = 5 discs than those in Q = ∞ discs.

Figure 14.

Azimuthally averaged density perturbations and brightness temperature of AS209AD and AS209ADSG at 400 t_p. From top to bottom, they are the density perturbations of the gas, dusts with St_max = 3 × 10⁻³ and St_max = 3 × 10⁻². The bottom panels show the brightness temperature of the St_max = 3 × 10⁻² (a_max = 0.5 mm) dust models. The left-hand panels show the cases without self-gravity, whereas the right panels show the cases with Q = 5. The blue, orange, green, and red curves represent T_cool = 0.01, 0.1, 1, and 100, respectively. The T_cool = 0.01 curves on the left-hand panels are almost identical to the locally isothermal cases in Zhang et al. (2018) (the slight difference is mainly due to the different setups in resolution, temperature profile, indirect term, and smoothing length).

Note that aside from the EoS, the simulation setups in this subsection are different from those in Zhang et al. (2018) in several aspects (e.g. the resolution, domain size, temperature profile, indirect term, smoothing length, and planet growing time). Thus, the gap shape and position are expected to be slightly different from the simulations in Zhang et al. (2018). Also, as time evolves, the position of the secondary and tertiary gaps will move slightly inwards. Thus, we do not try to compare these profiles with AS 209 observations. Instead, we use this example to demonstrate the observable differences due to different cooling times.

We also want to emphasize that without the self-gravity, the location of the planet and the gaseous surface density from the simulation can be scaled to any value. However, when the disc self-gravity is included, the gas density and the location of the planet are related. The relation is |$\Sigma _{g,0} = (M_*/r_\mathrm{ p}^2)(h/r)_\mathrm{ p}/(\pi Q_\mathrm{ p})$|⁠. Taking M_* = 0.83 M_⊙, r_p = 99 au, Q_p = 5 and (h/r)_p = 0.05, we get Σ_{g, 0} = 2.4 g cm⁻². In Zhang et al. (2018), the best-fitting model has Q ≈ 2 at the planet’s position in the initial condition.

Assuming the gaps at 24, 35, and 61 au in AS 209 are related to the planet at 99 au, we can use these gaps to constrain the AS 209 disc mass. Using equation (39), and adopting the stellar luminosity as 1.4 L_⊙ and the mass as 0.87 M_⊙ from Andrews et al. (2018), ϕ = 0.02, T_f = 20 K, and κ_R = 0.21 cm²g⁻¹, we have

$$\begin{eqnarray} \begin{split} T_{\mathrm{ cool}} = 0.02 \Bigg [1+ 0.01\Big (\frac{\Sigma }{\mathrm{1\ g\ cm^{-2}}} \Big)^2 \Bigg ] \left\lbrace \begin{array}{@{}l@{\quad }l@{}}2.5 & r \le 46\ \mathrm{au} \\ \Big (\frac{r}{99\ \mathrm{au}}\Big)^{-3/2} & r \gt 46\ \mathrm{au} \\ \end{array}\right.\, . \end{split} \end{eqnarray}$$

(43)

To produce an observable tertiary gap, we need low or high cooling times (T_cool ≲ 0.1 or T_cool ≳ 10). However, T_cool ≳ 10 requires Σ ≳ 200 g cm⁻², which is impossible since Q would be ≪1. If T_cool ≲ 0.1, we can constrain Σ ≲ 20 g cm⁻². While it is an independent constraint, this does not give tighter constraint on the density for this disc, since from gravitational instability criterion (equation 8), the disc will become unstable if Σ_{g, 0} > 12 g cm⁻² at 99 au. Nevertheless, this constraint from the disc cooling is consistent with that from the disc instability.

The Q = 5 models with low T_cool might be plausible models for the AS 209 disc, since their density at 99 au is allowed given the constraints derived above. However, with the density Σ_{g, 0} = 2.4 g cm⁻² from the Q = 5 disc, we cannot restore the intensity measured from the DSHARP observation. It is much lower than that from the observation. Only Σ_{g, 0} ≳ 6 g cm⁻² can explain the intensity. This means that either the actual opacity is higher than what we adopt (e.g. the dust-to-gas mass ratio is larger) or the disc has even lower Q. If the latter is true, the T_cool can be larger (but still <1), since the low Q that is close to unity can produce much stronger perturbations, compensating for the weaker perturbations due to the larger value of the T_cool. Interestingly, Powell et al. (2019) estimate the gas surface density Σ ≈ 10 g cm⁻² at 100 au, using the ‘dust line’ method. This high surface density is actually compatible and consistent with the constraints from disc self-gravity and radiative cooling for the AS 209 disc. If this is the case, it indicates that the outer disc of the AS 209 disc is marginally gravitational stable.

5.4 Long time evolution for self-gravitating discs with planets

On the other hand, in our long time-scale simulations, we find that the self-gravitating discs with planets tend to become unstable at later times. This happens earlier if the Q is small or if the planet mass is high. Fig. 15 shows the density perturbations of B17SG discs when M_p = 0.3 and 1.0 M_th (top and bottom) at t = 100 t_p. Q = 2 discs have become unstable at this point. This is why low Q and high M_p curves are missing in Fig. 8. The disc would become more stable if the radiative cooling is included. For instance, the AS209ADSG simulations discussed in Section 5.3 can evolve for very long time without any sign of instability even though their Q = 5 and M_p = 1 M_th. Another interesting behaviour in the lower panels is that as the disc becomes more massive, the number of vortices at the gap edge becomes larger. This is consistent with the simulation results in Lin (2012), Zhu & Baruteau (2016), as self-gravity can suppress large-scale (small m) vortices produced by Rossby wave instability (Lovelace et al. 1999).

Figure 15.

The density perturbations of B17SG discs with M_p = 0.3 and 1.0 M_th (top and bottom) at t = 100 t_p. The disc becomes unstable at Q = 2. The lower panel shows that self-gravity can suppress small m vortices, so that the number of the vortices increases as Q becomes smaller.

6 CONCLUSIONS

We run 2D hydrodynamical simulations to explore the effects of self-gravity and radiative cooling on the gap and spiral formation in protoplanetary discs. We explore these effects with different Toomre Q parameters and the dimensionless cooling times T_cool for low mass planets. As the disc becomes more massive (smaller Q), we find:

• The spirals become slightly more tightly wound, especially when Q ∼ 2, which is consistent with the linear theory.

• The spirals become stronger and the primary and secondary gaps become deeper. This is due to the higher AMF in more massive discs, which is consistent with the linear theory. The advective AMF in the Q = 2 disc is almost five times the AMF in the Q = 100 disc.

• If Q ≳ 2, the secondary gap’s position does not change, despite gaps being deeper.

In discs with radiative cooling, we find:

• Even with a relatively short T_cool = 0.01, the spiral waves cannot pick up AMF from the background disc when they are propagating to hotter regions, different from locally isothermal discs.

• The spirals become less tightly wound (more open) as the cooling time increases, due to the increase of the disc’s effective sound speed. The spiral’s openness changes dramatically from discs with T_cool = 0.1 to T_cool = 1.

• The spirals become weaker as the cooling time increases from a small value to 1/Ω, and then start to become stronger as the cooling time increases. There is a significant change of the spiral’s amplitude when T_cool increases from ≲ 0.01 to 0.1.

• The AMF of the wave dissipates the fastest during its propagation when T_cool ∼ 0.1 to 1, implying that radiative damping may be important.

• With weaker spirals, the secondary gaps due to low mass planets disappear when T_cool = 0.1 to 1.

• If the secondary gap is present, its position moves from ∼ 0.5 r_p to ∼ 0.4 r_p in a disc with (h/r)_p = 0.07 in transition from the isothermal limit to the adiabatic limit.

• For Q ≳ 5 discs, the effect of radiative cooling is more critical than self-gravity.

Our simulations also have implications for the observations. We have run planet–disk simulations with massive planets (∼ thermal mass):

• One might overestimate the planet mass using the simulation that neglects self-gravity. The gap width is less sensitive to the disc self-gravity than the gap depth. Thus, it is better to constrain the planet mass using the gap width when disc self-gravity is important.

• One might underestimate the planet mass using the simulation that neglects the radiative cooling. This is especially the case when T_cool ∼ 1.

• For deep primary gaps in the high planet mass regime, their shapes are less affected by the cooling time, but the secondary gap’s position and depth for the gas depend on the cooling time. On the other hand, the secondary gap in dust or brightness temperature is less sensitive to the cooling time since dust drift fast to the primary gap edge forming secondary gaps. The tertiary gap becomes unnoticeable in both gas and dust when T_cool = 1.

• The dependence of the gap properties (e.g. gap depth, width, and secondary/tertiary gap’s position and depth) on the cooling time-scale provides a new way to constrain the gaseous disc surface density.

• Assuming the gaps at 24, 35, and 61 au in AS 209 are all associated with the planet at 99 au, the gas surface density at ∼100 au in AS 209 should be ≲20 g cm⁻² using the T_cool ≲ 0.1 cooling constraints, which is consistent with the surface density constraint from Q ≳ 1.

ACKNOWLEDGEMENTS

SZ and ZZ thank Chao-Chin Yang, Steve Lubow, and Jeffery Fung for helpful discussions. ZZ acknowledges support from the National Aeronautics and Space Administration through the Astrophysics Theory Program with Grant No. NNX17AK40G, the National Science Foundation under CAREER grant No. AST-1753168, and Sloan Research Fellowship. Simulations are carried out with the support from the Texas Advanced Computing Center (TACC) at The University of Texas at Austin through XSEDE grant TG-AST130002, and the NASA High-End Computing Program through the NASA Advanced Supercomputing Division at Ames Research Center.

The following softwares were used in this study: Dusty FARGO-ADSG (Baruteau & Masset 2008a,b; Baruteau & Zhu 2016), Athena++(Stone et al., in preparation), Matplotlib (Hunter 2007), Numpy (Van Der Walt, Colbert & Varoquaux 2011), Scipy (Jones et al. 01).

Footnotes

1

This is coined as angular momentum current in Binney & Tremaine (2008), denoted by capital ‘C’.

2

Aside from the interesting behaviour at T_cool ∼ 1, we want to note that at two extremes, the spiral at the adiabatic limit is still weaker than that at the locally isothermal limit as expected. This is related to the increase of the effective sound speed with a longer cooling time. When the sound speed increases, the thermal mass increases and ratio between the planet mass and the thermal mass decreases so that the planet is less capable of exciting spiral waves.

3

We also run simulations with constant α as the model (b) in Zhang et al. (2018). Same as previous models, they reproduce AS 209’s 1D profile, but are non-axisymmetric in 2D. We do not see more differences due to the viscosity profile.