The Dependence of the Structure of Planet-opened Gaps in Protoplanetary Disks on Radiative Cooling

Planets can excite density waves and open annular gas gaps in protoplanetary disks. The depth of gaps is influenced by the evolving angular momentum carried by density waves. While the impact of radiative cooling on the evolution of density waves has been studied, a quantitative correlation to connect gap depth with the cooling timescale is lacking. To address this knowledge gap, we employ the grid-based code Athena++ to simulate disk-planet interactions, treating cooling as a thermal relaxation process. We establish quantitative dependencies of steady-state gap depth (Equation 36) and width (Equation 41) on planetary mass, Shakura–Sunyaev viscosity, disk scale height, and thermal relaxation timescale (β). We confirm previous results that gap opening is the weakest when the thermal relaxation timescale is comparable to the local dynamical timescale. Significant variations in gap depth, up to an order of magnitude, are found with different β. In terms of width, a gap is at its narrowest around β = 1, approximately 10%–20% narrower compared to the isothermal case. When β ∼ 100, it can be ∼20% wider, and higher viscosity enhances this effect. We derive possible masses of the gas gap-opening planets in AS 209, HD 163296, MWC 480, and HL Tau, accounting for the uncertainties in the local thermal relaxation timescale.


INTRODUCTION
The gravitational perturbation induced by an orbiting planet within a protoplanetary disk (PPD) excites density waves and triggers a redistribution of material in the disk (Paardekooper et al. 2022).In the presence of a sufficiently massive planet, the surrounding material in its orbit can be effectively cleared, resulting in the formation of a gap.
The accurate characterization of gaps in PPDs holds significant importance due to its impact on various crucial physical quantities.For instance, the migration rate of a planet within a disk is governed by the planetary gravitational torque, sensitive to the surface density in the gap (Σ gap ) (Kanagawa et al. 2018;Paardekooper et al. 2022).Likewise, the accretion rate onto a planet is also influenced by Σ gap (Chen et al. 2021).Furthermore, a quantitative correlation between the planet mass, and the gap depth and width, can be employed to infer the former when the latter is measured (Kanagawa et al. 2017;Zhang et al. 2018) Previous studies have made significant progresses in exploring the structure of gaps opened by planets through 2D and 3D numerical simulations as well as semi-analytical models.These investigations have explored the evolution of angular momentum and the dependence of gap depth and width on disk and planetary characteristics, including the planet-star mass ratio (q ≡ M p /M ⋆ ), viscosity (α), and disk aspect ratio (h p ). Fung et al. (2014) derived the analytical expression for gap depth by balancing viscous torque and one-sided planetary torque, providing Σ gap /Σ 0 ∼ q 2 α −1 h −5 p .This relation was validated by numerical simulations and initiated parameter explorations in the gap-opening process.Subsequent studies have expanded upon these findings.Kanagawa et al. (2015b,a) and Duffell (2015) developed 1D semi-analytical models for steady-state gap structures, presenting more accurate formulas for gap depth based on angular momentum conservation.When considering massive planets, nonlinear effects become important in calculating excitation torques.Kanagawa et al. (2016) proposed that the mass of a planet in a PPD can be constrained based on gap width.Kanagawa et al. (2017) investigated the gap-opening process of deep gaps and angular momentum deposition using numerical simulations and semi-analytical models, incorporating nonlinear density waves.Furthermore, Kanagawa et al. (2017) defined the cumulative deposited torque (Γ d ) to study the ability of density waves to exchange angular momentum with local substances and proposed an iterative method for generating gap profiles in the nonlinear regime, reducing the deviation between numerical gap depth and previous analytic formulas.
In addition to these advancements, other explorations have yielded significant results.For example, Duffell & Dong (2015) compared gap depth and width between single and multiple planets, finding that multiple planets of the same mass created wider but shallower gaps.Fung & Chiang (2016) conducted 3D numerical simulations, demonstrating that 3D and 2D simulations produced similar gap structures.Combining hydrodynamics simulations with radiative transfer, Dong & Fung (2017) shows that the observed gap depth and width can constrain the quantity M 2 p /α. Duffell (2020) derived an empirical formula for rapidly generating synthetic disk profiles, which shows excellent agreement with simulations across a wide range of parameters.Recent research by Sánchez-Salcedo et al. (2023) examined the relationship between gap depth and planetary orbit eccentricity.
However, most previous works on gaps opened by planets in PPDs adopted an isothermal equation of state (EOS).In recent years, the importance of thermodynamics has been emphasized.The pioneering work of Miranda & Rafikov (2020a) first conducted linear theory and numerical simulations to reveal the importance of thermodynamics in understanding the evolution of density waves in accretion disk systems.Subsequently, Zhang & Zhu (2020) conducted numerical investigations to elucidate the behavior of density waves under different thermal relaxation timescales.Miranda & Rafikov (2020b) investigated the impact of in-plane radiative diffusion by introducing a combined cooling timescale from the effects of both surface and in-plane cooling.Furthermore, Huang & Yu (2022) explored the effect of cooling on the Rossby wave instability (RWI) in PPDs through linear analysis.Wang et al. (2023) initiated studies on how thermodynamics affects the accretion of binary systems using numerical simulations.In addition to the parameterized β−cooling method mentioned above, there are also studies that adopt realistic terms for heating and cooling in the energy equation (Ziampras et al. 2020(Ziampras et al. , 2023)).Non-negligible discrepancies have been identified between the former approach and the latter (Ziampras et al. 2023).
In this paper, we study the influence of radiative cooling on the depth, width, and internal structures of gaps opened by giant planets in PPDs.When searching for planets in PPDs, our models can provide more precise information about the properties of embedded planets.The lack of such a model (or even an empirical formula) hampers the determination of planet-induced gaps.We introduce a thermal-relaxation treatment with a constant timescale to investigate its impact on the gapopening process.While one planet may open multiple gaps (Dong et al. 2017(Dong et al. , 2018;;Bae et al. 2017;Bae & Zhu 2018), we focus on the primary gap around the planet's orbit, mainly because the viscosity utilized in our investigation generally exceeds the typical viscosity associated with multiple gap opening.Gaps form when density waves transfer angular momentum to the local disk, and the depth of the gap depends on the evolution of angular momentum carried by these waves.Considering a non-isothermal EOS leads to distinct dissipation processes for density waves and different gap structures compared to the locally isothermal condition.
The paper is structured as follows.In Section 2, we present the setups and outcomes of a numerical parameter exploration focusing on the gap-opening process.Section 3 discusses the propagation of density waves and the evolution of torques and angular momentum.Furthermore, we derive empirical formulas for gap depth and width as a function of planetary mass q, viscosity coefficient α, disk aspect ratio h p , and cooling time scale β in the quasi-steady state.Section 4 contains discussions and applications to observed gas gaps related to the obtained results.

Method
To investigate the structure of steady-state gaps opened by planets in disks with radiative cooling, we employ the grid-based code Athena++ (Stone et al. 2020) to solve the conservative equations for mass, momentum, and energy.These equations governing the behavior of the system are expressed in a polar coordinate system (R, ϕ): The continuity equation: The momentum equation: The energy equation: (3) In these equations, Σ represents the surface density of the disk, v is the velocity, P denotes the isotropic pressure, and E is the total energy.The term I represents the identity tensor.The viscous stress tensor is given by: Here, ν = αc s H is the kinematic viscosity, where α is the Shakura-Sunyaev viscosity parameter and c s is the sound speed.The gravitational potential is given by Φ(R, ϕ, t) = Φ ⋆ (R)+Φ p (R, ϕ, t), where Φ ⋆ (R) represents the stellar gravitational potential and Φ p (R, ϕ, t) represents the planetary gravitational potential.We adopt the adiabatic equation of state, where the total energy E is defined as E = 1 2 Σv 2 + P/(γ − 1), with the adiabatic index γ = 1.4.Lastly, the term ∂ϵ ∂t cool represents the thermal relaxation term, which will be discussed in more detail later in the paper.

Initial and boundary conditions
Our 2D simulation domain ranges from R min = 0.2 to R max = 2.4 radially (in units where the planet-star distance R p = 1), and from 0 to 2π in azimuth.The aspect ratio of the disk is where H is the disk scale height.The initial surface density and temperature profiles of the disk are where Σ 0 and T 0 are the unperturbed values at R p , set to 1 and h 2 p .The initial velocity profiles are the sub-Keplerian velocity modified by a radial pressure gradient along the azimuthal direction: and the velocity of viscous inflow along the radial direction: The boundary conditions are set to match the initial conditions at the respective inner and outer ghost cells, which are positioned outside the computational mesh.Additionally, we incorporate a wave damping zone near the inner and outer boundaries to mitigate unphysical wave reflections at these boundaries (de Val-Borro et al. 2006).

Gravitational terms
To mitigate the initial strong shock from the massive planet, we gradually grow the planetary mass to its full extent M p (Bae et al. 2017): Here, t grow represents the growth time: where t p = 2πΩ −1 p is the planetary orbital period, with Ω p = Ω(R p ).The gravitational potential of the planet is denoted by Φ p (R, ϕ, t), where R p and Ω p represent the radial position and orbital frequency of the planet, respectively, and R s = 0.6R H is the softening length.Here, is the Hill radius of the planet.The gravitational potential of the star is denoted as Φ ⋆ (R, ϕ, t) = −GM ⋆ /R.If the planet has a significant mass, the center of mass is shifted away from the star.To maintain the star as the origin of the co-rotating non-inertial reference frame, we introduce the indirect potential as an extra correction of the effective gravitational potential (Fung et al. 2014;Fung & Chiang 2016).

Thermal relaxation term
A common approach to treat cooling and heating is to introduce thermal relaxation with a characteristic timescale t cool to achieve a background temperature profile the same as the initial condition (Zhu et al. 2015;Zhang & Zhu 2020;Paardekooper et al. 2022).We follow this practice and introduce thermal relaxation as: In the above equation, the cooling timescale t cool can be conveniently expressed in a non-dimensional form as t cool = βΩ −1 K , where Ω K represents the Keplerian angular velocity and β is the dimensionless cooling timescale.We treat β as a constant in each simulation that does not vary with time or space.

Simulation Runs
Athena++ provides two reconstruction methods, piecewise linear method (PLM, xorder = 2), and piecewise parabolic method (PPM, xorder = 3).We found no significant difference in the steady state gap profiles using these two methods (see Figure 1).As PLM is much less computation time cost, we adopt it as the default reconstruction method in this work.In all simulations, we selected the Riemann solver HLLC and the Van Leer 2 (VL2) as the time integrator.

Parameter Space
The main simulation parameters are summarized in Table .1.Specifically, we vary • the planet-star mass ratio q from 3 × 10 −4 to 3 × 10 −3 , • the viscosity parameter α from 3×10 −4 to 1×10 −2 , • the disk aspect ratio at the planet location h p from 0.05 to 0.1, and • the cooling time scale β from 0.01 to 100.

Convergence Tests
The resolution convergence tests are shown in Figure 1.We find that simulations with both β = 1 and β = 100 converge at the grid resolution (N R , N ϕ ) = (384, 1152), and adopt this as the default.With this resolution, approximately 9 grid cells per scale height are resolved at the planet location for h p = 0.05, and around 17 grid cells per scale height for h p = 0.1.
Each simulation runs until the surface density inside the gap reaches a quasi-steady state.We employ the relation, as proposed in Kanagawa et al. (2018)  To characterize the radial gap structure and quantify the gap depth, we introduce three measurements.First, we define the radial surface density profile, denoted as Σ(R), which is obtained by azimuthally averaging the surface density while excluding the azimuth around the planet: We set δϕ = π/12.The azimuthal position of the planet is ϕ p = 0. Next, we determine the gap depth, denoted as Σ gap , which corresponds to the average surface density within a narrow ring centered at the planet's location.Specifically, we compute Σ gap using the following expression: Here, ∆ represents the half-width of the narrow ring, set to be max(R H , h p ), the maximum value between the Hill radius R H and the scale height at the planet's location h p .This is half the value set in Fung et al. (2014).We do so to obtain the gap depth closer to the center of the gap.
In addition, we determine the radial width of the gap by locating its inner (R in ) and outer edges (R out ) corresponding to the positions where the surface density reaches the geometric mean of the minimum and undepleted values (Dong & Fung 2017), The gap width is then
Figure 4 shows the numerically-derived gap depth (Σ gap ) and width (w gap ) as a function of cooling timescale (β).The simulated gap depth and width with K ∼ 10 1 to K ∼ are listed in Table 4.

EMPIRICAL MODELS FOR GAP DEPTH AND WIDTH
In this section, we propose an empirical model that incorporates thermodynamic considerations to describe the gap depth.Drawing upon previous studies, we aim to generalize the gap depth formula to include the impact of the cooling time scale, β.To verify and provide insights into this model, we numerically investigate the time evolution of torques and angular momentum flux (AMF) with their dependence on β.

Torque Balance and Angular Momentum Transfer
The evolution of angular momentum carried by density waves plays a crucial role in shaping the disk.Therefore, understanding the dynamics of angular momentum flux (AMF) and torque balance is essential for comprehending the origins of different gap structures.
Consider a region of the disk spanning from the planetary orbit at R p to an arbitrary radius R near the gap's outer edge (see Figure 5 for details).The total AMF, denoted as F J , can be divided into three distinct components (Kanagawa et al. 2017): (1) AMF resulting from accretion flow:  Azimuthally averaged surface density at the planetary orbit, Σ(Rp), as a function of time for 4 simulations with (q, α, hp) = (3 × 10 −3 , 10 −3 , 0.08) and different β.The viscous time, as given by Eq. ( 15), is marked with a red vertical line.
where 2πR⟨v r Σ(R)⟩ ϕ represents the accretion rate at radius R.
(2) AMF originating from viscous diffusion: where ν denotes the kinematic viscosity, Σ represents the surface density, and Ω and v R denote the local angular velocity and radial velocity, respectively.The angled brackets ⟨⟩ denote azimuthal averaging.
(3) AMF arising from non-axisymmetric advection due to density waves (Zhang & Zhu 2020;Miranda & Rafikov 2020a,b): where δv ϕ and δv R represent the deviations in azimuthal and radial velocities from initial condition, respectively, caused by the presence of the planet.
The process of gap formation can be understood as a competition between the planetary torque, which aims to open gaps, and viscous diffusion, which tends to close gaps (Armitage 2020).Considering the region from R p to R (the red region in Figure 5), the conservation of angular momentum can be expressed as: where L represents the angular momentum of the gas in this region.The cumulative torque within this radius range is denoted as Γ p (R): The total angular momentum outflow from this region is given by F J (R) − F J (R p ).In a steady state, angular momentum conservation dictates that dL dt = Γ p (R) − F J (R) + F J (R p ) = 0. Kanagawa et al. (2017) pointed out that it is the cumulative planetary deposition torque, defined as: rather than Γ p (R), that governs the gap-opening process.The former represents the absorbed portion of the latter within the annular region (R p , R) of the disk.
To examine the dependence of gap depth on planetary mass and disk properties, we employ the conservation of angular momentum in the annulus extending from R p , the radius of the planet, to R + , the outer edge of the gap (Kanagawa et al. 2015a).Once a steady state is achieved, the deposition torque in this annulus should balance the net angular momentum flux entering the region from the inner and outer boundaries, which includes the contributions from accretion flow (F acc J ) and For each value of K, the gap depth and width are normalized with respect to the values obtained at β = 0.01, approximating the isothermal case.The viscous parameter α varies from 10 −4 to 3 × 10 −3 .In the range of K from approximately 10 2 to 10 4 , we observed that the gap reaches its maximum depth and minimum width at around β = 1.As K increases, the change in gap depth becomes more pronounced, whereas the change in width becomes less significant.The data used in this figure is tabulated in Table 4.

Planet
Figure 5.A sketch map depicting the components of the total AMF at an arbitrary radius R. The conservation of angular momentum is investigated within the red ring region (taking the outer disk for example).The total AMF flow out the region is viscous diffusion (F vis J ).This can be expressed as: A classic assumption, as proposed in Armitage (2020), is that the gap is wide enough to absorb all of the onesided planetary torque (Γ d (R + ) ≈ Γ p (R + )), and it is narrow enough that R + ≈ R p and Ω K (R + ) ≈ Ω p .Under this assumption, the accretion-driven AMF entering and leaving the region cancels out (F acc J (R + ) − F acc J (R p ) ≈ 0).The cumulative deposition torque balances the viscous diffusion at the inner (R p ) and outer (R + ) boundaries of the annular gap region.This can be written as: Considering a disk rotating at Keplerian velocity , we can approximate the viscous diffusion term F vis J (R) as: By evaluating Equation ( 28) at R + and R p , we obtain F vis J (R + ) and F vis J (R p ) in Equation ( 27): Note that the unperturbed surface density at , and they cancel out with R 4 + .Then, the viscous AMF at R + becomes We introduce the one-sided Lindblad torque as an approximation for the cumulative deposition torque at R + (Armitage 2020): where f 0 = 0.12π is a constant derived from the integration of the WKB torque formula (Lin & Papaloizou 1979).Substituting Equations ( 29) and ( 31) into Equation ( 27), we obtain: which provides a more precise formula for the gap depth: where  33as an accurate predictor of gap depth for given values of (q, α, h p ).Because in our sample most planetary masses are relatively large, the excitation of nonlinear density waves by these planets poses a challenge.As a result, the gaps predicted by Eq. ( 33) tend to underestimate the actual gaps observed in simulations.To address this issue, modifications to Eq. ( 33) are necessary.

Nonlinear effects
When K > 10 3 , due to the nonlinearity of density waves excited by the planets, the gap will be much wider than in linear cases, so that Eq. (31) will no longer be a good approximation for Γ d (R + ).As shown in Figure 15, the numerically simulated Γ d at the edge of the gap will be greater than the estimated value from Eq. ( 31), and this effect becomes more pronounced as K increases.As a result, if nonlinear effects are not taken into account, the predicted gap will appear shallower.
To tackle this problem, Kanagawa et al. (2017) proposed a model that incorporates a nonlinear factor, obtained by multiplying the torque density, and employs Eq. ( 33) Eq. ( 35) Eq. ( 37) , Linear, with cooling Eq. ( 36) , Nonlinear, with cooling an iterative algorithm to generate the gap profile.However, for the sake of simplicity in our fitting formula, we adopt an exponential modification e λK for the term representing the planetary deposition torque (Γ d (R + )), which is recognized as the driving force for gap-opening at the outer edge of the gap.Eq. ( 33) thus becomes This modification allows us to incorporate the nonlinear effects when K > 10 3 .

Formula for gap depth considering the cooling time scale
In order to incorporate thermal effects, we introduce a dependence of the cooling timescale β on λ 0 and f 0 .The modified formula for the gap depth, considering the cooling time scale, is given by: Here, the constant coefficent f 0 = 0.12π in Eq. ( 31) corresponding to linear, isothermal cases is replaced by f (β)e λ(β)K .This coefficient represents the ratio between the cumulative deposition torque Γ d (R + ) acting on the region from the planet's position to the periphery of the gap, and We utilize Equation (36) to fit the gap depth Σ gap , obtained from approximately 100 numerical simulations considering the cooling effect.
In the pursuit of a fitting formula, we initially employ Equation 36 to fit the gap depth for a specific β and derive the corresponding f and λ values.
In Figure 6, we present two fitted curves for gap depth: Equation ( 36), which incorporates a correction for the nonlinear effect, and  which shares the same form as Equation ( 33) but replaces f 0 with f (β), which depends on β. Figure 6 illustrates that even when considering the cooling effect, Equation (36) continues to provide a satisfactory fit but with a new coefficient f (β)e λ(β)K instead of f 0 .In addition, Equation (37) gives good performance when K < 10 3 .The comparison of different deposition coefficient are shown in Fig. 15.
The fitted values of f and λ are marked by circles in Figure 7. Subsequently, we adopt the functional form f (β) = log 10 β −1 +c1 c2β+β −1 ϵ +f 0 to fit the coefficients f and λ as functions of the cooling timescale β.The results are presented in Equations ( 38) and (39).
The deposition factor peaks at β ≈ 1, indicating that a cooling timescale comparable to the dynamical timescale enhances the absorption of excitation torque by the disk substance compared to the adiabatic and locally isothermal cases, as previously shown by Miranda & Rafikov (2020a,b); Zhang & Zhu (2020); Ziampras et al. (2023).
Figure 8 illustrates the relationship between the gap depth and the cooling time scale β for planets with varying gap-opening abilities.As the gap-opening capacity K increases from ∼ 10 1 to ∼ 10 4 , the cooling time scale corresponding to the minimum value of Σ gap shifts from 0.3 to 3. We also compared the performance of the fitted gap depth curves between low-viscosity and highviscosity disks.For slow cooling cases (β ≥ 10) with a viscous parameter α ≥ 3 × 10 −3 , the surface density in the gap may be underestimated due to viscous heating.
All fitting formulas and their applicable ranges are summarized in Table 2.The torque balance (i.e., Eq. ( 27)) revealed in the gap depth formula Eq. ( 36) are numerically verified.Γ p , Γ d , F wave J , and F vis J are performed according to Equations ( 24), ( 25), ( 22

Gap Width
An empirical formula for the dependence of gap width on β is provided.Using isothermal simulations, Dong & Fung (2017) found where η is approximately a constant, 5.8.To take into account the cooling effect, we retain the format of the fitting function, but treat η as a function of β, i.e., Figure 9 shows the fitted η values at a variety of β.
Finally, η(β) is fitted as and shown in Figure 10.Note that it has the same functional form as in f (Eq.38) and λ (Eq.39).The comparison between the gap width predicted by the fitting formula and numerical simulations is presented in Figure 11 for three different h p .

Applications to real PPDs
By incorporating a thermal relaxation term during our numerical simulations, we have generalized the results regarding the gap depth and width for massive planets under the locally isothermal condition (Fung et al. 2014;Kanagawa et al. 2015aKanagawa et al. , 2016Kanagawa et al. , 2017)).We provide approximate analytical expressions for the gap depth and width in terms of α, q, h p , and β, which accurately describe the gap depth and width within the range 10 −4 ≤ α ≤ 3 × 10 −3 , 3 × 10 −4 ≤ q ≤ 3 × 10 −3 , 0.05 ≤ h p ≤ 0.1, and 10 −2 ≤ β ≤ 10 2 .
Our empirical formulas for gap depth and width carries implications for observations.In recent years, numerous gap and cavity structures have been discov-  18)) versus the disk aspect ratio hp.The cooling time scales are 0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 100, from left to right.Simulations are marked with circles, while the dashed lines are the fitted linear relation between the two quantities (Eq.( 41)), whose slope (η) is labeled at the bottom.
ered by ALMA, not only in dust but also in gas (van der Marel et al. 2016;Yen et al. 2016;Isella et al. 2016).These observational measurements, when combined with theoretical models, allow for the constraints of various properties, such as planetary mass, for the planets residing within the disk.
For example, utilizing spatially resolved CO isotopologue emission observed by ALMA, van der Marel et al. (2015,2016); Dong et al. (2017) derived gas surface density distributions and identified significant drops up to several orders of magnitude in gas surface density inside central cavities in a few disks.This discovery implies the presence of close-in Jovian planets.While recent observations have achieved order of magnitude accuracy in characterizing the gap or cavity depth in gas, our study reveals that the gap depth can differ by roughly an order of magnitude when considering different thermal relaxation time scales.This finding highlights the significant 10 3 10 2 10 1 10 0 10 1 10 2 10 3 0.2 0.3 0.4 0.5 0.6 0.7 0.8
impact of the radiative cooling in theoretical models on constraining planetary mass.
The depth and width of CO gaps in five PPDs were reported by Zhang et al. (2021) from the ALMA large program Molecules with ALMA at Planet-forming Scales (MAPS), among which three (AS 209, HD 163296, and MWC 480) were reported to exhibit significant CO gaps.Gaussian functions were used in Zhang et al. (2021) to fit the CO gap profiles and measure the gap depth and width in these three PPDs.In the remaining two PPDs, GM Aur has a cavity, while IM Lup has a very wide inner gap (see Figure 18 in Zhang et al. (2021) for details.)In addition, an HCO+ gas gap in HL Tau was found in Yen et al. (2019).
Under the assumption that each gap is opened by one planet, we applied our fitting formula to the gas gaps mentioned above (see Table 3, Column (1)).We derive the expected planetary mass for 5×5 = 25 combinations of α and β, shown in Figure 12, based on Eq. 36 with observed gap depth and disk scale height in Table 3.Our α ranges from 10 −2 to 10 −4 , and β ranges from 0.01 to 100.
Furthermore, the fitting formula for gap width (i.e., Eq. 41) can be used to constrain β.Among the observed gas gaps listed in Table 3, we found it possible to estimate β in 3 gaps in AS 209, HD 163296, and MWC 480 by finding the intersection between the observed w gap (dashed line in Fig. 13) and Eq.41 (solid curve in Fig. 13).For the remaining gaps, our assumption may not be valid due to insufficient data to accurately measure the gap width.
It is crucial to emphasize that the observed gap properties pertain to CO and HCO + gaps, and whether they differ from those of H 2 gaps remains uncertain.Besides, our estimates should be regarded as indicative values, considering that real disks may deviate from the assumptions in our models, e.g., uniform viscosity and cooling timescale across the gap region.
Accounting for cooling in observations, a basic consequence is that the isothermal models (e.g., Fung et al. (2014); Duffell (2015); Kanagawa et al. (2017), etc.) lead to overestimation of a planet's mass when cooling timescale is around the dynamical timescale (t cool ≈ t dyn ) and underestimation when it exceeds 10t dyn .For instance, in a disk with (α, h p ) = (10 −3 , 0.05), the isothermal model predicts a 1M J mass planet responsible for a deep gap with a 1.9 × 10 −3 drop of gaseous surface density.However, our thermal relaxation model suggests that a gap with the same depth will be created by a 1.3M J planet with β = 100, or a 0.9M J planet with β = 1.
Additionally, when considering cooling, the temperature distribution of the disk deviates from the locally isothermal scenario due to radiative cooling, viscous heating, and other factors.Muley et al. (2021) demonstrates that temperature structures in spiral arms can serve as an additional theoretical reference for observations.The temperature profile across a gap (e.g., Fig. 3) may also encode information on the thermodynamics in the disk.

Limitations and outlook
Firstly, we have simplified our simulations by focusing solely on the evolution of gas distribution and neglecting the influence of dust feedback on the gas.Previous studies have shown differences in gaseous structures due to dust feedback, especially for wider and deeper gaps Observed gas gap width defined by the ratio between the FWHM and the location of the gap.Colume(5): disk aspect ratio at the gap location.
Note-The errors of the observed quantities are not considered in our applications.(Kanagawa 2019).Disk self-gravity is another factor that could affect the evolution of the angular momentum flux (AMF) (Zhang & Zhu 2020).Therefore, an important future direction based on our research is to explore the effects of disk self-gravity and radiative cooling on gas-dust interactions.
Besides, we adopted a simple assumption of constantβ thermal relaxation for the sake of parametric simplicity.In reality, β may be a function of location in disks (Zhang & Zhu 2020;Muley et al. 2023).
On top of that, our β-cooling models, although convenient for parameter exploration, has limitations when applied to ALMA disks.Recent studies have highlighted the strong influence of in-plane cooling on spiral wave AMF and gap opening induced by planets (Miranda & Rafikov 2020b).After that, Ziampras et al. (2023) found discrepancies between fully radiative models with direct cooling and models considering only surface cooling.By setting β as an effective cooling timescale that incorporates both surface and in-plane cooling (e.g., Miranda & Rafikov 2020b), our model is also able to take in-plane cooling into consideration, which could give a rough correction and make the results of AMFs, torques, and gap depth similar to those of fully radiative models with direct cooling term.However, even with this modification, the output of our adjusted model will still differ from those with realistic cooling terms.As direct constraints on cooling timescales in observations are still limited (e.g., Teague et al. 2022), we defer detailed studies in  3. The planetary mass is determined using Eq. ( 36) and inferred Σgap/Σ0 and hp from observations, with viscosity (α) ranging from 10 −4 to 10 −2 and cooling time scale (β) ranging from 10 −2 to 10 2 .
the effect of anisotropic cooling on gap properties to the future.

SUMMARY
Through extensive 2D hydrodynamical simulations, we have explored the effect of radiative cooling on gas structures within protoplanetary disks (see Figure 3 for a direct visualization).We have quantified the dependence of gap depth and width on the cooling timescale β (Eq.36, Fig. 4, and Table 4), a dimensionless quantity normalized by the dynamical timescale (t dyn = Ω −1 ), for planets with masses ranging from Saturn to 3 Jupiter masses, roughly equivalent to 0.6 to 6 disk thermal masses in a disk with a local aspect ratio h p = 0.08.The gap is narrowest and deepest at β ≈ 1, as found in Miranda & Rafikov (2020a,b) and Zhang & Zhu (2020), and widest and shallowest at β = 100.Significant variations in the gap depth, up to an order of magnitude, were observed for different β values.Additionally, a deviation of 10% ∼ 20% in width is observed compared to the isothermal case.A line chart comparing the depth and width of gaps in disks with β cooling and those in nearly locally isothermal disks are shown in Figure 4.A table listing the simulated gap depth and width for a wide range of gap-opening capacity can be found in Table 4.
Our numerical calculations have allowed us to evaluate the radial distribution of angular momentum fluxes and torques, confirming the torque balance achieved within the system as it reaches a quasi-steady state (see Figure 14).We have also discovered that the discrepancy in steady-state gap depth corresponding to different cooling time scales can be attributed to the varying ability of density waves in propagation (Figure 16).Notably, around the critical value of β = 1, the density wave experiences the strongest damping closest to the planet's position, leading to a more rapid decay of its carried angular momentum flux with increasing radius.This suggests that when the cooling timescale is comparable to the dynamical timescale, the disk material shows the highest absorption efficiency for the planetary excitation torque, as opposed to cases of immediate cooling (isothermal limit) and slow cooling (adiabatic limit).As a result, this leads to the formation of the deepest and narrowest gap.
In addition, by fitting the surface density near the planet's orbit obtained from our numerical simulations, we have generalized the formula for gap depth, Eq. ( 33), derived from torque balance for a steady-state gap in the linear and isothermal case.By introducing a correction factor e λK for nonlinear effects and another correction factor f (β) to account for the effect of cooling (Table 2), we have formulated an analytical expression (Eq.( 36)) that effectively captures the relationship between gap depth and (q, α, h, β) in both linear and nonlinear regimes.Gap depth as a function of β is shown in Figs. 8.This is an extension of previous works concerning the dependence of gap depth on (q, α, h) alone (Fung et al. 2014;Kanagawa et al. 2017).
An empirical formula for gap width with cooling effect (Eq. 41), as an extension to the isothermal case reported in Dong et al. (2017), is also provided.The gap edges are defined as the locations where the surface density reaches the geometric mean of the minimum inside the gap and the undisturbed value.Gap width as a function of β is shown in Fig. 11.
We applied our empirical model (Eq. 36 & 41) to the observations of 6 CO gaps reported in the ALMA large program MAPS by Zhang et al. (2021) and 1 HCO + gap observed in Yen et al. (2019).We inferred the masses of the possible planets responsible for opening these gaps (Figure 12) and estimated the cooling time scale for 3 of them (Figure 13).Over time, the deposition torque Γ d decreases to the same magnitude as the viscous torque F vis J , and the region of overlap between deposition torque and viscous torque gradually expands outward from the planetary orbit.This expansion continues until approximately 5000 orbits when the deposition torque and viscous torque of the planet balance within a radial region of approximately 0.5R p on either side of the planetary orbit.
The right column of Figure 16 shows the AMF and torques in a quasi-steady state after evolving for ∼ 5000 orbits.The angular momentum stored in the ring region from R p to any arbitrary radius R is precisely balanced by viscous diffusion at the boundaries R and R p .This indicates that the system approaches a steady state and torque balance prevails throughout the gap.Moreover, by comparing the simulated gap depth with the theoretical predictions, we observe that the gap depth remains constant thereafter.
As illustrated in Figure 14, when the system reaches a steady state, the profile of viscous AMF (F vis J (R)) closely resembles the cumulative deposited torque (Γ d (R)).This suggests that in a protoplanetary disk with a stable gaseous gap formed by a planet, the density wave excited by the planet's gravity locally deposits a portion of its angular momentum into the disk material.The angular momentum gained by the ring region from R p to R is promptly dissipated at the ring's boundary through viscosity.The remaining angular momentum continues to propagate with the density wave until it is fully damped.For a relatively deep gap (Σ(R p )/Σ 0 << 1), the viscous diffusion near the gap (F vis J (R p ) ∼ 3πνΣ(R p )Ω p R p ) can be neglected compared to F vis J (R + ) ∼ 3πνΣ 0 Ω p R p .
B.2.The dependence of quasi-steady-state AMFs and torques on β We investigate the steady-state radial profiles of the planetary excitation torque (Γ p ), deposition torque (Γ d ), and AMF due to density wave (F wave J ) to better understand the empirical formula in Equation 36.As mentioned in section 3.3, the deposition coefficient represents the ratio between the theoretical value of the deposition torque Γ d (R + ) and the normalized flux F J0 , i.e., Γ d (R + )/F J0 = f (β)e λ(β)K . (B1) To verify this physical significance of the correction terms in Eq. 36, we numerically calculate the deposition torque (Γ d /F J0 ) near the gap edge (R + ) and compare it with three cases: (1) the linear, isothermal coefficient f 0 , (2) the corrections f (β)e λ(β)K .Results are shown in Figure 15.We observed that the non-linear effect becomes pronounced at K ≳ 10 4 , where f 0 e λ0K can be ∼ 10 × f 0 .In addition, the coefficient f (β)e λ(β)K , incorporating both non-linear and cooling adjustments, displays an initial increase followed by a subsequent decrease as β varies.This behavior mirrors the increasing-decreasing pattern observed in the simulated Γ d (R + ).Because the gap width varies with β, the position of R + will also become closer to the planetary orbit when β ≈ 1 and farther when β < 1 or β > 1, which is consistent with Figure 15.(2)f 0 e 0 K (1)f 0 R/R p =1.40,1.45,1.50,1.55,1.60 simulations (q, , h p ) = (10 3 , 3 × 10 4 , 0.08), K = 1017 By conducting simulations with different cooling timescales, we observe that the density wave corresponding to β ∼ 1 peaks at the position closest to the planet's orbit and fails to propagate further (see Figure 16).Similar phenomena regarding the dependence of AMF damping on the cooling time scale have been explained in Miranda & Rafikov (2020b); Zhang & Zhu (2020); Speedie et al. (2022); Ziampras et al. (2023), among others.However, the excitation torque Γ p in the β = 1 case is slightly larger than in other cases, resulting in a more significant increase in deposition torque Γ d with increasing R. As Γ d (R + ) peaks at β = 1, the gap depth Σ gap reaches its minimum value according to Equation 27.

Figure 1 .
Figure 1.Radial surface density profiles in simulations with different resolutions and spatial difference orders xorder.The disk and planet parameters are listed on the top.The left two panels have xorder = 2 and right panel shows comparison between xorder = 2 and xorder = 3.
is a parameter that quantifies the planet's ability to open a gap.Numerical simulations conducted by Kanagawa et al. (2015a), Kanagawa et al. (2017), Fung et al. (2014), and others have verified Equation

Figure 6 .
Figure 6.Averaged value of the azimuthally averaged surface density in the gap area Σgap versus the parameter K.The cooling time scales are 0.01,0.1,1,10,100,from left to right.Simulations are marked with circles, where different colors represent varying cooling time scales.The red and gray dotted lines are fitted curves according to Equations (36) & (37) without and with the correction factor for nonlinearity, respectively.

Figure 8 .
Figure 8.The dependence of gap depth (left) and gap profile (right) on the thermal relaxation timescale β for different gap-opening parameters K. Left panel: gap depth (Σgap/Σ0) as a function of cooling time scale (β).Circles are simulation results.The black dotted curve is Eq.(36).Right panel: the azimuthally averaged gap profile with β ranging from 0.01 to 100.

Figure 13 .
Figure 13.Predicted values of β for three PPDs according to Eq. (41) and observed gap width (wgap) of (a) the 94 au gap of AS 209, (b) the 98 au gap of HD 163296, and (c) the 63 au gap of MWC 480, provided in Table3.

Figure 15 .
Figure15.The cumulative deposition torque evaluated near the gap edge with respect to β. Circles are Γ d at various R near the gap edge yielded from numerical simulations.The gray, black, and blue dotted lines symbolizes the deposition factors related to the (1) linear, isothermal (Eq.33), (2) nonlinear, isothermal (Eq.35), and (3) nonlinear, with cooling (Eq.36) cases, respectively, symbolizing the theoretical deposition torque at the periphery of the gap (Γ d (R+)).
(Dong & Fung 2017)scous gap-opening timescale t ν , assuming the gap width w gap /R p ≈ 5.8h p(Dong & Fung 2017).Considering substance in the gap (R ≈ R p ), t ν is approximately 1.3t p /α.For each model, the total running time t run , greater than t ν in all cases, is listed in Table1.Figure2shows how Σ(R p ) varies with time in a run with (q, α, h p ) = (3 × 10 −3 , 10 −3 , 0.08).In this case, the surface density at the planet's orbit, Σ(R p ), converges at ∼ 10 3 t p .Note that the corresponding viscous time for gap-opening according to Equation 15 is ∼ 1.3 × 10 3 t p .

Table 2 .
Fitting formulas for gap depth

Table 3 .
Applying Our Empirical Formulas to Real PPDs

Table 4 .
Simulated gap depth and width p ) experiences a relatively rapid decrease.Consequently, both the planet's excitation torque Γ p and the AMF due to density wave F wave J decrease as well.Meanwhile, the radial profile of viscous torque F vis J undergoes only minor adjustments as a result of changes in surface density.