Chemical footprints of giant planet formation

Haochang Jiang; Yu Wang; Chris W. Ormel; Sebastiaan Krijt; Ruobing Dong

doi:10.1051/0004-6361/202346637

Home

All issues

Volume 678 (October 2023)

A&A, 678 (2023) A33

Full HTML

Open Access

Issue		A&A Volume 678, October 2023


Article Number		A33
Number of page(s)		18
Section		Interstellar and circumstellar matter
DOI		https://doi.org/10.1051/0004-6361/202346637
Published online		03 October 2023

A&A 678, A33 (2023)

Role of planet accretion in shaping the C/O ratio of protoplanetary disks^★

Haochang Jiang (蒋昊昌)¹^,2, Yu Wang (王雨)², Chris W. Ormel², Sebastiaan Krijt³ and Ruobing Dong (董若冰)⁴

¹ European Southern Observatory, Karl-Schwarzschild-Str 2, 85748 Garching, Germany
e-mail: hjiang@eso.org
² Department of Astronomy, Tsinghua University, 30 Shuangqing Rd, Haidian DS 100084, Beijing, PR China
e-mail: wang-y21@mails.tsinghua.edu.cn
³ School of Physics and Astronomy, University of Exeter, Stocker Road, Exeter EX4 4QL, UK
⁴ Department of Physics and Astronomy, University of Victoria, Victoria, BC V8P 5C2, Canada

Received: 9 April 2023
Accepted: 16 July 2023

Abstract

Context. Protoplanetary disks, the birthplaces of planets, commonly feature bright rings and dark gaps in both continuum and line emission maps. Accreting planets interact with the disk, not only through gravity, but also by changing the local irradiation and elemental abundances, which are essential ingredients for disk chemistry.

Aims. We propose that giant planet accretion can leave chemical footprints in the gas local to the planet, which potentially leads to the spatial coincidence of molecular emissions with the planet in the ALMA observations.

Methods. Through 2D multi-fluid hydrodynamical simulations in Athena++ with built-in sublimation, we simulated the process of an accreting planet locally heating up its vicinity, opening a gas gap in the disk, and creating the conditions for C-photochemistry.

Results. An accreting planet located outside the methane snowline can render the surrounding gas hot enough to sublimate the C-rich organics off pebbles before they are accreted by the planet. This locally elevates the disk gas-phase C/O ratio, providing a potential explanation for the C₂H line-emission rings observed with ALMA. In particular, our findings provide an explanation for the MWC 480 disk, where previous work identified a statistically significant spatial coincidence of line-emission rings inside a continuum gap.

Conclusions. Our findings present a novel view of linking the gas accretion of giant planets and their natal disks through the chemistry signals. This model demonstrates that giant planets can actively shape their forming chemical environment, moving beyond the traditional understanding of the direct mapping of primordial disk chemistry onto planets.

Key words: astrochemistry / protoplanetary disks / planet-disk interactions / accretion / accretion disks / hydrodynamics / stars: variables: T Tauri / Herbig Ae/Be

Note to the reader: the co-author's name "Chris W. Orme" was incorrect in the HTML version, it was corrected to "Chris W. Ormel" on 23 October 2023.

^★

Movie associated to Fig. 4 is available at https://www.aanda.org

© The Authors 2023

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Protoplanetary disks are the “cradles” of planets. Imaging observations at high resolution and high contrast using instruments such as ALMA and Gemini/GPI have found that most large disks (r > 50 au), if not all, exhibit substructures such as bright rings and dark gaps (e.g., Andrews et al. 2018; Long et al. 2018; van der Marel et al. 2019). It is widely speculated that planet-disk interactions are responsible for many of these substructures (e.g., Zhu et al. 2011; Dong et al. 2015; Rosotti et al. 2016; Zhang et al. 2018; Paardekooper et al. 2022). Nevertheless, detecting planets in these young systems is challenging (Xie et al. 2020; Huélamo et al. 2022; Follette et al. 2023). Except for PDS 70 (Keppler et al. 2018; Haffert et al. 2019) and AB Aur (Currie et al. 2022; Zhou et al. 2022), most claims of planet detections are indirect (e.g., Pinte et al. 2019, 2020; Speedie et al. 2022; Izquierdo et al. 2022) and might be debatable (Speedie & Dong 2022; Fedele et al. 2023).

Among the indirect planet detection method, the role of astrochemistry is gaining increasing prominence with the development of facilities such as ALMA and JWST (e.g., Öberg & Bergin 2021). The classical core accretion planet formation dictates that a Neptune-mass planet quickly grows into a Jovian-mass gas giant during a runaway growth phase, in which gas freely falls into its gravitational potential. During the runaway phase the planet's mass and accretion rate Ṁ_p exponentially grow, making the protoplanet an irradiating source in the disk (e.g., Rafikov 2006; Mordasini et al. 2012).

Luminous young planets locally heat up their vicinity, resulting in efficient thermal sublimation of molecular species that are originally frozen out. Therefore, localized submillimeter molecular emission might appear around the planet, as proposed in Cleeves et al. (2015). This scenario is particularly compelling for the recently identified local emission enhancements of SO in HD 100546 (Booth et al. 2023), and ¹³CO in AS 209 (Bae et al. 2022). However, the calculation of Cleeves et al. (2015) used a static gas density profile, where the differential rotation of material around the planet is ignored. For a gap-opening planet at distance r₀, as the planet opens a wider gap with width Δw as it is accreting, the synodic timescale around the planet gap r₀/ΔwΩ_K, where Ω_K is the Keplerian frequency, gets shorter. If this timescale is shorter than the typical chemical timescale, material can be sheared around the co-orbital region of the accreting planet, and the asymmetric planet-heated chemical emission might be elongated and become arc-like or even ring-like.

Thanks to the high spatial resolution and high sensitivity of ALMA, ring-like substructures are found to be ubiquitous in line-emission maps as well (e.g., Öberg et al. 2021). Among the species with line-emission detection, the radial C₂H, one of the most well-studied hydrocarbons, is believed to indicate an active photochemistry, driven by high UV fluxes and gas with a C/O ratio exceeding unity (Bergin et al. 2016; Kama et al. 2016; Miotello et al. 2019; Bergner et al. 2019; Alarcón et al. 2021; van der Marel et al. 2021; Bosman et al. 2021a). C₂H can serve as a precursor for larger molecule formation through chemical reaction networks (e.g., Bast et al. 2013; Öberg & Bergin 2021). Alternatively, C₂H may also be a product of top-down photodissociation hydrocarbon chemistry (Guzmán et al. 2015).

Due to the specific formation conditions required for C₂H (a high C/O ratio and high UV penetration) and its bridging role in the chemical reaction network, we explore whether a young giant planet can be responsible for achieving these conditions, and discuss how the chemical footprint of the accretion of the planet manifests in protoplanetary disks.

In this study we conduct 2D multi-fluid hydrodynamical simulations that include built-in sublimation to simulate an accreting planet in Athena++ (Huang & Bai 2022; Wang et al. 2023). The planet is modeled as a sink particle in the simulation with self-consistent accretion prescription. As the planet accretes, the potential energy of the accreted gas is released as the planet’s luminosity, locally heating the gas surrounding the planet and releasing C-rich volatiles, whose abundance we follow in Athena++. We then post-process our simulation outputs with an empirical formula derived from chemical network simulations on C₂H formation. In addition, we compare our simulation results with line-emission observations from ALMA.

The plan of the paper is as follows. The theoretical background and numerical setup are described in Sect. 2. We present and explain the numerical results in Sect. 3. The results are further discussed in Sect. 4. We summarize the main results and conclusions in Sect. 5.

2 Model

We simulate the planet and disk interaction with multi-fluid 2D hydrodynamics simulations in polar coordinates (radius r and azimuth ϕ). The disk consists of three mass components: noncondensable gas, ice-rich pebbles, and sublimated vapor. To model them, we introduce gas fluid (g), particle fluid (p) and tracer fluid (tr) respectively. The gas fluid is the mixture of non-condensable gaseous species (mostly H and He) and sublimated vapor. There is only one gas fluid. In addition, the density of each vapor species is followed by a tracer fluid, which inherits the hydrodynamical properties (i.e., velocity field) of the gas fluid. Pebbles are modeled as compounds of multiple particle fluids, each of which represents a solid component (i.e., ice and refractory). With the recently developed multi-dust fluid module of Athena++ (Huang & Bai 2022), the gas and ice can co-evolve in a self-consistent manner. The governing equations read (see also Sect. 2.2 of Wang et al. 2023)

$\frac{\partial ρ_{g}}{\partial t} + \nabla \cdot (ρ_{g} υ_{g}) = 0,$ ${{\partial {\rho _{\rm{g}}}} \over {\partial t}} + \nabla \cdot \left( {{\rho _{\rm{g}}}{{\bf{\upsilon }}_{\rm{g}}}} \right) = 0,$ (1)

$\frac{\partial (ρ_{g} υ_{g})}{\partial t} + \nabla \cdot (ρ_{g} υ_{g} υ_{g} + P_{g} I + Π_{v}) = ρ_{g} f_{g,src},$ ${{\partial \left( {{\rho _{\rm{g}}}{{\bf{\upsilon }}_{\rm{g}}}} \right)} \over {\partial t}} + \nabla \cdot \left( {{\rho _{\rm{g}}}{{\bf{\upsilon }}_{\rm{g}}}{{\bf{\upsilon }}_{\rm{g}}} + {P_{\rm{g}}}{\bf{I}} + {{\bf{\Pi }}_v}} \right) = {\rho _{\rm{g}}}{{\bf{f}}_{{\rm{g,src}}}},$ (2)

$\frac{\partial ρ_{p}}{\partial t} + \nabla \cdot (ρ_{p} υ_{p} + ℱ_{p,dif}) = 0,$ ${{\partial {\rho _{\rm{p}}}} \over {\partial t}} + \nabla \cdot \left( {{\rho _{\rm{p}}}{{\bf{\upsilon }}_{\rm{p}}} + {{\rm{F}}_{{\rm{p,dif}}}}} \right) = 0,$ (3)

$\begin{array}{l} \frac{\partial ρ_{p} (υ_{p} + υ_{p,dif})}{\partial t} + \nabla \cdot (ρ_{p} υ_{p} υ_{p} + Π_{dif}) = \\ ρ_{p} f_{p,src} + ρ_{p} \frac{υ_{g} - υ_{p}}{t_{s}}, \end{array}$ $\matrix{ {{{\partial {\rho _{\rm{p}}}\left( {{{\bf{\upsilon }}_{\rm{p}}} + {{\bf{\upsilon }}_{{\rm{p,dif}}}}} \right)} \over {\partial t}} + \nabla \cdot \left( {{\rho _{\rm{p}}}{{\bf{\upsilon }}_{\rm{p}}}{{\bf{\upsilon }}_{\rm{p}}} + {{\bf{\Pi }}_{{\rm{dif}}}}} \right) = } \hfill \cr {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{\rho _{\rm{p}}}{{\bf{f}}_{{\rm{p,src}}}} + {\rho _{\rm{p}}}{{{{\bf{\upsilon }}_{\rm{g}}} - {{\bf{\upsilon }}_{\rm{p}}}} \over {{t_{\rm{s}}}}},} \hfill \cr }$ (4)

$\frac{\partial ρ_{tr}}{\partial t} + \nabla \cdot (ρ_{tr} υ_{g} + ℱ_{tr,dif}) = 0.$ ${{\partial {\rho _{{\rm{tr}}}}} \over {\partial t}} + \nabla \cdot \left( {{\rho _{{\rm{tr}}}}{{\bf{\upsilon }}_{\rm{g}}} + {{\rm{F}}_{{\rm{tr,dif}}}}} \right) = 0.$ (5)

In the above, p is the density, v is the velocity, P_g is the gas pressure, and I is the identity tensor. The gas viscosity is described by the viscous stress tensor

$Π_{v} \equiv ρ_{g} v_{g} [\nabla υ_{g} + {(\nabla υ_{g})}^{T} - \frac{2}{3} (\nabla \cdot υ_{g}) I],$ ${{\bf{\Pi }}_v} \equiv {\rho _{\rm{g}}}{v_{\rm{g}}}\left[ {\nabla {{\bf{\upsilon }}_{\rm{g}}} + {{\left( {\nabla {{\bf{\upsilon }}_{\rm{g}}}} \right)}^T} - {2 \over 3}\left( {\nabla \cdot {{\bf{\upsilon }}_{\rm{g}}}} \right){\bf{I}}} \right],$ (6)

where the kinematic viscosity v_g is specifically expressed as a function of the dimensionless constant α (Shakura & Sunyaev 1973)

$v_{g} = α c_{s} H_{g},$ ${v_{\rm{g}}} = \alpha {c_{\rm{s}}}{H_{\rm{g}}},$ (7)

where H_g is the gas scale height and c_s is the local sound speed. With this, the particle concentration diffusion flux is given by

$ℱ_{p,dif} = - ρ_{g} D_{p} \nabla (\frac{ρ_{p}}{ρ_{g}}) \equiv ρ_{p} υ_{p,dif},$ ${{\rm{F}}_{{\rm{p,dif}}}} = - {\rho _{\rm{g}}}{D_{\rm{p}}}\nabla \,\left( {{{{\rho _{\rm{p}}}} \over {{\rho _{\rm{g}}}}}} \right) \equiv {\rho _{\rm{p}}}{{\bf{\upsilon }}_{{\rm{p,dif}}}},$ (8)

which meanwhile defines the effective diffusive velocity υ_p,dif. The diffusivity for dust can be prescribed as (Youdin & Lithwick 2007)

$D_{p} \equiv \frac{v_{g}}{1 + {St}^{2}}$ ${D_{\rm{p}}} \equiv {{{v_{\rm{g}}}} \over {1 + {\rm{S}}{{\rm{t}}^2}}}$ (9)

with the Stokes number St = t_sΩ_K. In addition, Π_dif is the momentum diffusion flux tensor, which, combined with υ_p,dif, ensures the Galilean invariance (Huang & Bai 2022) in Eq. (4). Similar to Eq. (8), the tracer concentration diffusion flux is defined as

$ℱ_{tr,dif} = - ρ_{g} D_{tr} \nabla (\frac{ρ_{p}}{ρ_{g}})$ ${{\rm{F}}_{{\rm{tr,dif}}}} = - {\rho _{\rm{g}}}{D_{{\rm{tr}}}}\nabla \,\left( {{{{\rho _{\rm{p}}}} \over {{\rho _{\rm{g}}}}}} \right)$ (10)

with the tracer diffusivity D_tr = v_g, where we implicitly assume that the Schmidt number is 1. Without loss of generality, we assume the gas viscosity coefficient α to be 5 × 10⁻⁴, and initialize the pebbles with constant Stokes number St = 0.01 throughout the domain. By assuming spherical pebbles, the particles size

$s_{p} = \frac{2}{π} \frac{{StΣ}_{g}}{ρ_{•}}$ ${s_{\rm{p}}} = {2 \over \pi }{{{\rm{St}}{{\rm{\Sigma }}_{\rm{g}}}} \over {{\rho _ \bullet }}}$ (11)

is ~0.5 mm at the planet orbit, with internal density ρ_•. = 1.67g cm⁻³ (Birnstiel et al. 2018). These values match the latest constraint approaches based on the ALMA observational results (see the recent review by Rosotti 2023, and references therein). In our simulations setup, the particle size s_p is kept constant throughout the simulation. The final term in Eq. (4) represents the aerodynamic drag experienced by particles, which follow the Epstein regime throughout our simulation, with the backreaction on gas omitted in our formulation. External source terms are denoted by f_src, including stellar and planetary gravity (see Eq. (18)).

2.1 Disk model

We assume that the initial gas and pebble surface densities follow power-law distributions:

$\begin{array}{l} Σg = Σ_{g, 0} \times {(\frac{r}{r_{0}})}^{- γ}, \\ Σp = Σ_{p, 0} \times {(\frac{r}{r_{0}})}^{- γ} . \end{array}$ $\matrix{ {{{\rm{\Sigma }}_{\rm{g}}} = {{\rm{\Sigma }}_{{\rm{g}},0}} \times {{\left( {{r \over {{r_0}}}} \right)}^{ - \gamma }},} \hfill \cr {{{\rm{\Sigma }}_{\rm{p}}} = {{\rm{\Sigma }}_{{\rm{p}},0}} \times {{\left( {{r \over {{r_0}}}} \right)}^{ - \gamma }}.} \hfill \cr }$ (12)

Here r is the distance from the star and Σ_g,0 is the gas surface density at the reference location r₀. The gas disk is parameterized based on the observational constraints of MWC 480, a typical Herbig Ae/Be star with a clear annular gap in its continuum at ~73 au (Long et al. 2018; Liu et al. 2019). From thermochemical modeling based on line emission of CO isotopes, Zhang et al. (2021) suggest a gas surface density profile of Σ_g,0 = 11 g cm⁻², γ = 1 at r₀ = 73 au. The densities of the ice-rich pebble at r₀ is Σ_p,0 = ƒ_mΣ_g,0, where the mass fraction f_m depends on the volatile species of interest in each simulation (see Table 1 for the values of each species). The stellar masses have been dynamically determined to be ~2.1 M_⊙ as described in Teague et al. (2021).

We assume that the disk is isothermal in the vertical direction and employ a locally isothermal equation of state. The ambient temperature is determined by the disk background, which is influenced by the irradiation of the host star, and the heating from the accreting planet. Specifically, the disk background temperature is described as

$T_{d} = T_{0} \times {(\frac{r}{r_{0}})}^{- q},$ ${T_{\rm{d}}} = {T_0} \times {\left( {{r \over {{r_0}}}} \right)^{ - q}},$ (13)

where T₀ = 20 K, and q = 0.45 according to the power-law fit of the dust temperature at the midplane (Sierra et al. 2021), based on the 2D (R-Z) dust density and temperature profile in the thermochemical models in Zhang et al. (2021). The aspect ratio of the gas disk therefore follows

$h_{g} \equiv \frac{H_{g}}{r} \equiv \frac{c_{s}}{Ω_{K} r} = h_{0} \times {(\frac{r}{r_{0}})}^{0.5 - q / 2}$ ${h_{\rm{g}}} \equiv {{{H_{\rm{g}}}} \over r} \equiv {{{c_{\rm{s}}}} \over {{{\rm{\Omega }}_K}r}} = {h_0} \times {\left( {{r \over {{r_0}}}} \right)^{{{0.5 - q} \mathord{\left/ {\vphantom {{0.5 - q} 2}} \right. \kern-\nulldelimiterspace} 2}}}$ (14)

with h₀ = 0.05 based on the temperature profile. We note that this value is inconsistent with the best-fit aspect ratio based on the SED fitting. Zhang et al. (2021) fit the vertical temperature profile at each radius directly by matching the emission surface of CO isotope lines. Since our model focuses on the ice sublimation at the midplane, we opt for the midplane temperature in our simulations.

As the planet is accreting, the lost gravitational energy of the accreting materials can be translated into accretion luminosity

$L_{acc} = \frac{G M_{p} \dot{M}}{R_{p}},$ ${L_{{\rm{acc}}}} = {{G{M_{\rm{p}}}\dot M} \over {{R_{\rm{p}}}}},$ (15)

where M_p is the planet mass, R_p is the radius of the planet accretion shock front, and Ṁ is the mass accretion rate of the planet. We fix R_p = 2 R_Jupiter in all of our simulations, which is the typical shock front radius of an accreting planet in the disk (1.5– 4 R_Jupiter·; see, e.g., Marleau et al. 2022, and references therein). A detailed explanation for the prescription of planet mass M_p and accretion rate Ṁ is given in Sect. 2.3.

We utilize the accretion luminosity to prescribe the temperature surrounding the planet. We consider the radiation near the planet to be optically thin at infrared wavelength, which is typically the case in the planet-opened gaps at the outer regions of the protoplanetary disk. The temperature is thus prescribed as (e.g., Rafikov 2006)

$T = {(\frac{L_{acc}}{16 π σ d^{2}} + T_{d}^{4})}^{0.25},$ $T = {\left( {{{{L_{{\rm{acc}}}}} \over {16\pi \sigma {d^2}}} + T_{\rm{d}}^4} \right)^{0.25}},$ (16)

where d is the distance to the planet.

Table 1

Abundances n(X)/n(H) for different molecules X of interest.

2.2 Volatile sublimation

At each grid point, volatile ice can be sublimated into vapor via a recently developed phase change module implemented in Athena++ (Wang et al. 2023). This module solves mass transfer due to ice sublimation after each hydrodynamic step, thus enabling the self-consistent simulation of volatile transport in both ice and vapor states. The volatile follows the saturated (or equilibrium) vapor pressure, which is given by the Clausius-Clapeyron equation,

$P_{eq} = P_{eq, 0} e^{- T_{a} / T},$ ${P_{{\rm{eq}}}} = {P_{{\rm{eq}},0}}{e^{{{ - {T_{\rm{a}}}} \mathord{\left/ {\vphantom {{ - {T_{\rm{a}}}} T}} \right. \kern-\nulldelimiterspace} T}}},$ (17)

where T_a and P_eq,0 are constants depending on the molecular species. We take T_a = 1110 K and P_eq,0 = 3.67 × 10⁹ g cm⁻¹ s⁻² for methane, which is the major C-carrier icy volatile (where we neglect higher order correction terms; Fray & Schmitt 2009). When the partial pressure of the volatile vapor is lower than the saturated vapor pressure, ice is assumed to be instantly converted into vapor. In reality, the sublimation timescale of the volatile ∝ ds/υ_th, where ds is the depth of the volatile-rich surface layer of the pebble and υ_th is the thermal velocity of vapor particles (e.g., Schoonenberg & Ormel 2017). Thus, the instant sublimation of ice-rich pebbles is validated for the volatiles of interest to us, which are assumed to cover the surface of the pebbles. As the released vapor will mainly stay inside the gap but outside the co-orbital region where there are few pebbles for vapor to condense on, for simplicity we ignore condensation in the simulations (see Sect. 4.1 for further discussion).

In each simulation, pebbles are assumed to be composed of refractories, which dominate by mass, and a single ice species, which can sublimate from ice to vapor phase. Therefore, we assume that the size of the pebbles is not affected by sublimation. The abundances of all species included in our simulation are listed in Table 1. The values are marginally the same as those in Bosman et al. (2021a), and are consistent with cometary ices in forming protoplanetary disk midplanes (Drozdovskaya et al. 2016). CH₄ are ice-rich pebbles that are initially in the ice phase; in other words, the planet locates beyond their snowlines. Other materials, including CO, are assumed to be non-condensable gases throughout the simulation, which means that the planet is located inside their snowlines. For instance, as suggested in Zhang et al. (2021), the midplane CO snowline for MWC480 is around 100 au.

Fig. 1

Sketch of 2D mesh grid used in our simulations. The circle represents the location of the planet (at r₀ = 73 au) and the star is the central star. The radial domain ranges from 30 au to 150 au. The static mesh refinement level is 2, which means the smallest mesh is four times smaller than the unrefined mesh at the same location.

2.3 Growth of the planet

In our 2D simulation, we modeled the planet as a sink particle with a fixed circular orbit at a radius of r₀. The planet mass is initialized as 10 M_⊕, a typical critical core mass for a planet starting runaway gas accretion (e.g., Perri & Cameron 1974; Rafikov 2006; but see also recent works considering atmosphere pollution that suggest even lower values Hori & Ikoma 2011; Ormel et al. 2021). The static mesh refinement was employed throughout the simulation to better capture the gas flow and ice sublimation in the vicinity of the planet, as depicted in Fig. 1. The simulation domain extended radially from 30 au to 150 au, and the smallest mesh size was dx_min = 0.1R_Hill, where R_Hill = (M_p/3M_★)^1/3r₀ is the Hill radius of the planet. We used the value of R_Hill for a planet mass of M_p = M₀ = 2.3M_J, as suggested by previous studies reproducing the continuum profile of MWC 480 (Liu et al. 2019). We tried a smaller dx_min for a convergent test, but there were no changes to the simulation results. A softening parameter ϵ = 0.1R_Hill was chosen for the gravitational potential

$Φ = - \frac{G M_{p}}{\sqrt{d^{2} + ϵ^{2}}} .$ ${\rm{\Phi }} = - {{G{M_{\rm{p}}}} \over {\sqrt {{d^2} + {^2}} }}.$ (18)

The choice of ϵ = 0.1R_Hill ensures that the accretion rate on the planet is high; this feature is essential to capture with our model as it heats the vicinity of the planet, allowing pebbles to sublimate. In the literature, however, a larger softening parameter is favored for 2D simulations in order to correctly capture the wave excitation and gap opening by the planet (e.g., Dong et al. 2011; Müller et al. 2012). Our choice for ϵ in 2D simulations could therefore expedite the gap-opening process somewhat. However, conducting a run with ϵ = 0.6H_g, we found little difference in the gap profile, and therefore do not expect that the choice of ϵ would much affect our conclusions.

Starting with the initial planet mass of 10 M_⊕, we ran the simulation with the fixed-initial-mass planet for t_inin = 300 orbits, which is about 10⁵ yr at the distance of the planet. The time is sufficiently long compared with a total disk lifetime, and enough for the disk to adapt itself to the planet embedded in the disk. We tried different t_inin, which did not affect our main conclusion. After t_inin, we assigned a mass fraction of the ice-rich pebbles (Table 1), then switched on mass accretion in the simulations. We removed all pebbles entering 0.9 R_Hill from the simulation. The gas accretion prescription follows Crida & Bitsch (2017) and Bergez-Casalou et al. (2020), and the accretion rate is based on the 3D isothermal shearing box simulations results from Machida et al. (2010), which is determined by the local gas density within 0.9 R_Hill of the planet

${\dot{M}}_{p} = f_{A} {〈 Σg 〉}_{0.9 R_{Hill}} H_{g}^{2} Ω_{K} \times \min [0.14, 0.83 \times {(\frac{R_{Hill}}{H_{g}})}^{\frac{9}{2}}],$ ${\dot M_{\rm{p}}} = {f_{\rm{A}}}{\left\langle {{{\rm{\Sigma }}_{\rm{g}}}} \right\rangle _{0.9{R_{{\rm{Hill}}}}}}H_{\rm{g}}^2{{\rm{\Omega }}_K} \times \min \,\left[ {0.14,\,0.83 \times {{\left( {{{{R_{{\rm{Hill}}}}} \over {{H_{\rm{g}}}}}} \right)}^{{9 \over 2}}}} \right],$ (19)

where H_g is the gas scale height. The prefactor f_A equals unity in the default setup; we also tried different values (see Sect. 3.2). At each timestep, the accreted mass of gas was removed from the domain and added to the planet accordingly. Following Bergez-Casalou et al. (2020), the distribution of the mass removing rate ${\dot{Σ}}_{g,acc}$ ${{{\rm{\dot \Sigma }}}_{{\rm{g,acc}}\,}}$ (d, t) reads

$\begin{array}{l} {\dot{M}}_{p} = \iint_{d < 0.9 R_{Hill}} {\dot{Σ}}_{g,acc} (d, t) d A, \\ = \iint_{d < 0.9 R_{Hill}} f_{red} (d) {\dot{Σ}}_{g} (t) d A, \end{array}$ $\matrix{ {{{\dot M}_{\rm{p}}} = \int\!\!\!\int_{d lt; 0.9{R_{{\rm{Hill}}}}} {{{{\rm{\dot \Sigma }}}_{{\rm{g,acc}}}}} \left( {d,t} \right)\,{\rm{d}}A,} \hfill \cr {\,\,\,\,\,\,\,\, = \int\!\!\!\int_{d lt; 0.9{R_{{\rm{Hill}}}}} {{f_{{\rm{red}}}}\left( d \right){{{\rm{\dot \Sigma }}}_{\rm{g}}}\left( t \right)\,{\rm{d}}A} ,} \hfill \cr }$ (20)

where the smooth distribution function f_red reads

$f_{red} = {\begin{array}{l} 1 & d < 0.45 R_{Hill} \\ \cos^{4} π (\frac{d}{R_{Hill}} - 0.45) & 0.45 R_{Hill} < d < 0.9 R_{Hill} \end{array}$ ${f_{{\rm{red}}}} = \left\{ {\matrix{ 1 \hfill & {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,d lt; 0.45{R_{{\rm{Hill}}}}} \hfill \cr {{{\cos }^4}\pi \,\left( {{d \over {{R_{{\rm{Hill}}}}}} - 0.45} \right)} \hfill & {0.45{R_{{\rm{Hill}}}} lt; d lt; 0.9{R_{{\rm{Hill}}}}} \hfill \cr } } \right.$ (21)

following Robert et al. (2018). The fraction of accreted mass increases as it gets closer to the planet, but turns flat within the innermost 0.45R_Hill. The gas removal rate within 0.45R_Hill, ${\dot{Σ}}_{g} (t)$ ${{{\rm{\dot \Sigma }}}_{\rm{g}}}\left( t \right)$ , is obtained by equating Eqs. (19) and (20). In other words, the mass accretion rate of the planet is calculated through Eq. (19), and the accreted gas is removed from the surrounding protoplanetary disk according to Eq. (20).

2.4 Sketch of the chemical evolution

We propose a chemical model to interpret a potential spatial correlation between the chemical ring and dust continuum profile caused by a planet exterior to the methane snowline. This model hinges on the local increase in the C/O ratio, which is attributed to localized volatile sublimation caused by the luminosity of the accreting planet. A sketch of the model can be found in Fig. 2.

High concentrations of organics, for example C₂H and HCN, indicate an active photochemistry, driven by high UV fluxes and gas with a C/O ratio exceeding unity (Bergin et al. 2016; Bosman et al. 2021a). In our model a young giant planet is responsible for achieving these conditions. First, pebbles in the planet’s orbital region are accreted by the growing protoplanet.

Since the accreting planet is hot, carbon-rich volatile ices (e.g., methane) with sublimation temperatures around 25 K are readily liberated from the ice, elevating the C/O ratio of the gas in the vicinity of the planet’s orbit above unity. On the other hand, the sublimation of CO₂ (H₂O), which is the major oxygen carrier, requires a temperature in excess of ~55 K (120 K). Such a high temperature at 73 au can only be reached in the immediate vicinity of the planet (<R_Hill), where the water and carbon dioxide vapor can no longer escape the planet’s gravity (e.g., Johansen et al. 2021; Barnett & Ciesla 2022). In other words, the sublimation temperature difference between C-rich ices (e.g., CH₄, C₂H₂) and O-rich ices (e.g., CO₂, H₂O) renders the gas-phase C/O>1. Then, as the gap region becomes optically thinner due to the removal of material, stellar UV photons penetrate the gap more deeply and trigger photochemistry reactions. These processes may produce ring-like structures in C₂H (Bergin et al. 2016), one of the typical carbon-bearing radicals, and the more complex gas-phase formation of organic molecules compounds (Bast et al. 2013; Calahan et al. 2023).

Fig. 2

Schematic showing how a gap-opening planet triggers the C₂H ring associated with the continuum gap. (a) Before runaway growth, a small planet is embedded inside the gaseous disk and carbon-rich ice is frozen out on pebbles. (b) During runaway gas accretion, the planet grows rapidly, consumes the pebbles in its vicinity, and opens a gas gap in the disk. Since the accreting planet is hot, carbon-rich volatile ice (e.g., methane) with lower sublimation temperature is readily liberated from the ice. On the other hand, the sublimation of H₂O or CO₂, which are the major oxygen carriers, can only be reached in the immediate vicinity of the planet (<R_Hill), at which point the material is lost to the planet. The discrepancy in sublimation temperature renders the gas C/O>1. In addition, stellar UV photons penetrate the gas gap. The combination of C/O>1 and high UV fluxes trigger photochemistry reactions to form the C₂H ring. (c) After the gap opens, the supply of sublimated pebbles terminates due to the gap opening. The C/O ratio is smoothed out across the disk through radial mixing by diffusion.

2.5 C₂H chemistry

We focus on C₂H rings in this work. As a radical, C₂H can act as a precursor for the formation of larger molecules through photochemical reactions (e.g., Bast et al. 2013; Öberg & Bergin 2021). On the other hand, it has also been proposed that C₂H can be a product of top-down photodissociation hydrocarbon chemistry (Guzmán et al. 2015). From the observational side, Bergner et al. (2019) has suggested that C₂H and HCN seem to share a common driver. More research is needed to fully understand the connection between C₂H and the formation of complex organics in these environments, where a close link between a high C/O ratio and molecular formation appears.

To estimate the column density of C₂H, we adopt an empirical approach instead of relying on computationally expensive full chemical network modeling methods (e.g., Bruderer et al. 2012; Du & Bergin 2014). Specifically, we utilize the DALI models presented in Bosman et al. (2021a), which were calibrated to reproduce the observed C2H emission levels of AS209, HD 163296, and MWC 480 using source-specific gas-grain thermochemical models (see their Figs. 2 and 5). Based on the model results of n(C₂H) in the three sources and the input parameters in Bosman et al. (2021a), we summarize an empirical relation (see Appendix B)

$\log_{10} n (C_{2} H) = F_{C} - \log_{10} n (H) + 2 \log_{2} (C / O),$ ${\log _{10}}n\left( {{{\rm{C}}_2}{\rm{H}}} \right) = {F_C} - {\log _{10}}n\left( {\rm{H}} \right) + 2{\log _2}\left( {{{\rm{C}} \mathord{\left/ {\vphantom {{\rm{C}} {\rm{O}}}} \right. \kern-\nulldelimiterspace} {\rm{O}}}} \right),$ (22)

where F_C = 37 ± 1 is a fitting constant. Equation (22) allows us to estimate the C₂H column density based on our simulation outputs. The dependence of n(H) is the result of two effects: (a) lower CO column density, and therefore less self-shielding, allowing more dissociation of CO; and (b) the depletion of the small dust, which scales with the total gas density, increasing the UV penetration. The third term on the RHS of the formula directly reflects the strong dependence of C₂H on the carbon-to-oxygen ratio C/O, as shown in Fig. 2 of Bosman et al. (2021a). Nevertheless, more comprehensive chemical modeling may rely on a series of additional assumptions and prescriptions (see discussion in Sect. 4.1).

Finally, we remark that although we chose the fiducial values of our model parameters based on the observational constraint on MWC 480, the model we propose here is in principle applicable to any planet located exterior to its natal disk methane snowlines (see discussion in Sect. 4.2).

Fig. 3

Time evolution of different parameters in the fiducial model. (a) Planet mass. The horizontal gray line indicates the planet mass M₀ = 2.3M_J inferred by Liu et al. (2019), which we use as a reference mass. (b) Planet’s gas accretion. (c) Planet’s accretion luminosity (Eq. (15)). (d) Disk temperature at the Hill radius. The sublimation temperature of CH₄ (green) and CO₂ (orange) are indicated by horizontal lines for reference (see discussion in Sect. 4.2). In each panel, the vertical dashed lines indicate the time when the mass of the planet reaches M₀.

3 Results

In this section, the results of the 2D hydrodynamical simulations that implement the model introduced in Fig. 2 are presented. The fiducial model is described in Sect. 3.1, and a parameter variation of the gas accretion rate of the planet is conducted in Sect. 3.2.

3.1 Fiducial model

In our fiducial model, we include methane as the only ice-rich pebble. The time-dependent changes in both the mass and the accretion rate of the planet following the initiation of gas accretion are presented in panels a and b of Fig. 3.

For the first 20 000 yr (≈50 orbits) after the onset of runaway gas accretion, the accretion rate shows an exponential increase. During this period, the gas surface density at the location of the planet does not significantly decrease. The reason is that the timescale for the opening of the gas gap by the planet, which is the viscous timescale over the gap width (Kanagawa et al. 2017)

$t_{gap} = 2. \times 10^{3} \frac{2 π}{Ω_{K}} (\frac{M_{p} / M_{⋆}}{10^{- 4}}) {(\frac{h}{0.05})}^{- 3.5} {(\frac{α}{5 \times 10^{- 4}})}^{- 1.5}$ ${t_{{\rm{gap}}}} = 2. \times {10^3}{{2\pi } \over {{{\rm{\Omega }}_K}}}\left( {{{{{{M_{\rm{p}}}} \mathord{\left/ {\vphantom {{{M_{\rm{p}}}} {{M_ \star }}}} \right. \kern-\nulldelimiterspace} {{M_ \star }}}} \over {{{10}^{ - 4}}}}} \right)\,{\left( {{h \over {0.05}}} \right)^{ - 3.5}}{\left( {{\alpha \over {5 \times {{10}^{ - 4}}}}} \right)^{ - 1.5}}$ (23)

is longer than the orbital timescale. The slope of the accretion rate as a function of time slightly changes around t = 20 000 yr, because of the transition from Hill radius scaling toward constant prefactor 0.14 in Eq. (19), causing the accretion rate to flatten. Since the planet has reached a mass above 1 M_J at t ≃ 30 000 yr, the planet starts to sculpt the gas gap more strongly. Therefore, as the gas available around the planet decreases with the deeper and wider gap, the accretion rate of the planet goes down rapidly with a scaling of Ṁ_p ∝ t⁻². A dot-dashed line is presented in Fig. 3b for reference to show the trend.

As depicted in panels c and d in Fig. 3, by using the planet’s mass and gas accretion rate, we calculate the corresponding luminosity of the planet (Eq. (15)) and the temperature at its Hill radius

$\begin{array}{l} T_{Hill} = {(\frac{L}{16 π σ R_{Hill}^{2}} + T_{d}^{4})}^{\frac{1}{4}} = {(\frac{G {(3 M_{⋆})}^{2 / 3} M_{p}^{1 / 3} {\dot{M}}_{p}}{16 π σ R_{in} r^{2}} + T_{d}^{4})}^{\frac{1}{4}}, \\ = 43 K {(\frac{M_{⋆}}{M_{⊙}})}^{\frac{1}{6}} {(\frac{M_{p}}{M_{J}})}^{\frac{1}{12}} {(\frac{{\dot{M}}_{p}}{10^{- 4} M_{J} {yr}^{- 1}})}^{\frac{1}{4}} {(\frac{R_{in}}{2 R_{J}})}^{\frac{1}{4}} {(\frac{r}{73 au})}^{- \frac{1}{2}}, \end{array}$ $\matrix{ {{T_{{\rm{Hill}}}} = {{\left( {{L \over {16\pi \sigma R_{{\rm{Hill}}}^2}} + T_{\rm{d}}^4} \right)}^{{1 \over 4}}} = {{\left( {{{G{{\left( {3{M_ \star }} \right)}^{{2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-\nulldelimiterspace} 3}}}M_{\rm{p}}^{{1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-\nulldelimiterspace} 3}}{{\dot M}_{\rm{p}}}} \over {16\pi \sigma {R_{{\rm{in}}}}{r^2}}} + T_{\rm{d}}^4} \right)}^{{1 \over 4}}},} \hfill \cr {\,\,\,\,\,\,\,\,\,\,\, = 43\,{\rm{K}}\,{{\left( {{{{M_ \star }} \over {{M_ \odot }}}} \right)}^{{1 \over 6}}}{{\left( {{{{M_{\rm{p}}}} \over {{M_{\rm{J}}}}}} \right)}^{{1 \over {12}}}}{{\left( {{{{{\dot M}_{\rm{p}}}} \over {{{10}^{ - 4}}{M_{\rm{J}}}{\rm{y}}{{\rm{r}}^{ - 1}}}}} \right)}^{{1 \over 4}}}{{\left( {{{{R_{{\rm{in}}}}} \over {2{R_{\rm{J}}}}}} \right)}^{{1 \over 4}}}{{\left( {{r \over {73\,{\rm{au}}}}} \right)}^{ - {1 \over 2}}},} \hfill \cr }$ (24)

where R_Hill is the Hill radius of the planet. The second equation assumes that T_d can be neglected. Due to the high accretion rate, the temperature inside the Hill radius is notably elevated compared to the background disk temperature, rising from 20 K to approximately 45 K at the Hill sphere (see also Fig. 4e). Between t = 20000 yr (≈50 orbits) and t = 150000 yr (≈350 orbits), balancing the increased planet mass and the reduced amount of available gas, the accretion rate of the planet remains in the range of 10⁻⁵ to 10⁻⁴ M_J yr⁻¹ during its growth, which roughly corresponds to the temperature range where T_Hill above the sublimation temperature of CH₄ (the horizontal green line in Fig. 3d, see Sect. 4.2 and Fig. 5), while it is still below the sublimation temperature of CO₂ (the horizontal orange line in Fig. 3d, see also Sect. 4.2 and Appendix A). In addition, we show the temperature purely contributed from the planet by taking T_d = 0 in Eq. (24), which is indicated by the red line in Fig. 3d. As the figure shows, in the region where T > 25 K, the temperature at Hill radius T_hill is primarily determined by the planet’s luminosity, rendering the results insensitive to the background temperature once it is below the CH₄ sublimation temperature.

The simulation snapshot of the surface densities of all the components, the temperature map, and the C/O ratio map in the gas phase at t = 64 560 yr (176 orbits) are displayed in Fig. 4, along with the corresponding azimuthally averaged radial density profiles shown in Fig. 6. At this time, the planet has reached a mass of 2.3 M_J, consistent with the planet mass suggested by Liu et al. (2019), which we utilized as an educated estimate and point of reference. The Hill radius of the planet is indicated by red circles in each panel. The top panel of Fig. 4 shows the azimuthally averaged radial profile of different components. A moderate gas gap has been opened by the planet (upper left panel of Fig. 4), and pebbles have followed the pressure gradient, resulting in dust being trapped at the outer edge of the gas gap (upper right panel of Fig. 4). The dust gap is quickly emptied, while a small fraction of pebbles remain inside the co-orbital region of the planet, where a faint arc appears. This amount of dust is also shown in the radial profile Fig. 6. The opened gap allows for increased UV penetration, and therefore possibly higher temperature (Broome et al. 2023), which may help with the sublimation of the ice trapped inside the co-orbital region, neither of which is included in our simulation (but see Sect. 4.1 for more discussion).

The CH₄ ice sublimates and transfers to vapor phase close to the planet location, where the temperature is higher compared with the disk background temperature. The majority of CH₄ vapor is transported away from the planet through the gas spiral arm, as shown in panels a and c of Fig. 4. However, a significant fraction of vapor is transported from the vicinity of the planet into the gap as well, though not exactly into the co-orbital region. When considering the CH₄ vapor mass fraction over the total gas density, the vapor fraction reaches a relatively high value of ~10⁻⁵, as shown in Fig. 4c. Although the amount of CH₄ vapor that ends up in the planet gap region is rather low, it is nevertheless sufficient to raise the C/O ratio above unity. Using the assumed constant CO abundance of 8.9 × 10⁻⁶ and the CH₄ vapor density, we calculate the C/O ratio in the gas phase in Fig. 4f. The azimuthally averaged C/O ratio is indicated by the black solid line in the top panel of Fig. 6. Initially, the disk is cold and there is no methane released by the planet. Thus, both carbon and oxygen in the gas phase are determined by the non-condensable CO in our model. Therefore, the gas-phase C/O ratio is unity throughout the entire simulation box at the beginning. As CH₄, the major C-carrier volatile in our simulation, is released from ice-rich pebbles into the gas, the C/O ratio increases accordingly.

Significant vapor release happens over a relatively short timescale at the beginning where T_Hill is above the sublimation temperature of CH₄, as shown in Fig. 3d, and the gap remains shallow. In Fig. 7 we show the time evolution of the radial profiles of the C/O ratio in the gas phase. As the planet rapidly grows and gradually opens the gap, ice-rich pebbles can no longer reach the planet, and the sublimated vapor will be smeared out by radial mixing as shown in Fig. 7. In reality, the CH₄ vapor may eventually be destroyed and converted into other species (e.g., Eistrup et al. 2016) or condense back into ice. Neither process is included in the simulations (see Sect. 4.1).

Both high UV radiation and high C/O ratio are beneficial for C₂H formation. With the simulation outputs, we estimate C₂H column densities via Eq. (22) in the bottom panel of Fig. 6. A clear ring profile peaks inside the dust gap. We labeled the peak C₂H column densities of AS 209, HD 163296, and MWC480 estimated from fitting the C₂H lines hyperfine structure for reference (Guzmán et al. 2021), where the observational results are from the Molecules with the ALMA at Planet-forming Scales (MAPS) large program (Öberg et al. 2021). For all three systems, the column density of C₂H typically peaks ~10¹⁵ cm⁻² at the local maximum. Due to the simplification of the empirical formula Eq. (22), we also show the 1 mag uncertainty of the fitting constant F_C by the gray dashed region. Our model result predicts a peak C₂H column density value that is by-and-large consistent with the MAPS observation results. Moreover, the scenario described by the model is testable in the high spatial resolution result of ALMA, as a significant spatial correlation between the C₂H rings and the continuum gap is expected. We discuss this in Sects. 4.4 and 4.5.

3.2 Influence of the gas accretion rate

The precise gas accretion rates onto planets remain an open question. Studies differ on the accretion rates attained during the runaway phase, with values ranging from ~5 × 10⁻⁹ to ~3 × 10⁻⁴ M_J yr⁻¹ (Machida et al. 2010; Szulágyi et al. 2014; Tanigawa et al. 2014; Tanigawa & Tanaka 2016; Homma et al. 2020; Szulágyi et al. 2022). To investigate how the accretion rate changes our results, we varied the gas accretion formula (Eq. (19)) by adopting different prefactors (ƒ_A = 0.05,0.1,0.2,0.5,1,2,5) in the equation.

The planet masses and accretion rates are shown in Fig. 8a and b separately, with different colors indicating different pref-actors ƒ_A. We also show the gap depth (the gas surface density at the planet-opened gap Σ_gap over the unperturbed surface density Σ_g,0) as functions of time in Fig. 8c. At early times, when the gap remains shallow (Σ_gap/Σ_g,0 ~ 1), all simulations show an exponentially growing accretion rate. The accretion rate increases with larger ƒ_A, resulting in a more massive planet within the same period after the onset of gas accretion. At the same time, the simulations with higher accretion rates open the gas gap faster through a combination of increased accretion and stronger torques from the more massive planet pushing the gas away. As a result, it reaches the maximum accretion rate at an earlier turning point. In Figs. 8d and e, the corresponding planet luminosity and temperature at the Hill radius T_Hill at each run are calculated. The increased mass and accretion rate leads to a higher temperature around the planet, which reaches 50 K in simulation with f_A = 2 and 60 K in simulation with f_A = 5, causing more sublimation of ice-rich pebbles in the planet’s vicinity within the same time. Since the accretion rates decay in a similar fashion, Ṁ_p ∝ t⁻², after passing the peak, a higher accretion rate peak means that the temperature around the Hill radius of the planet remains higher than the CH₄ sublimation temperature for a longer time.

We present the maximum gas-phase C/O ratio versus time in Fig. 8f. Initially, when the planet is accreting slowly, carbon in the gas phase is dominated by CO, the C/O ratio is unity, and UV photons cannot penetrate the shallow gap. These conditions are unfavorable for the formation of C₂H rings. In simulations with faster accretion rates (more massive planets), the C/O ratio in the gas phase rises faster and earlier. However, the maximum C/O ratio in the gas phase is limited by the (fixed) amount of carbon ice in the pebble reservoir. By virtue of the assumption that the total CO and CH₄ abundances (gas+dust) are comparable, the maximum C/O ratio becomes ~2 when methane ice fully sublimates into vapor. Additionally, a more massive planet opens a deeper gap, allowing more UV radiation penetration and promoting chemistry, which both facilitate C₂H formation. Hence, a higher gas accretion rate is preferred to generate C₂H rings.

Conversely, in simulations with f_A < 0.2, the required growth timescale of the planet is much longer. The planet’s mass grows slowly, resulting in a relatively low accreting luminosity, even at the accretion rate peak (Fig. 8e). Consequently, only a fraction of methane ice sublimates to vapor, leading to a peak C/O ratio of only 1.6 in the run with f_A = 0.2, and even lower, barely exceeding unity, in the runs with f_A = 0.1 and f_A = 0.05. As the gap becomes deeper, pebbles are pushed away from the planet and are trapped at the outer edge of the gap. In all simulations, the gas-phase C/O ratio around the continuum gap is then gradually smoothed out due to the viscous evolution of the gas, as shown in Fig. 7, and the peak value decreases.

We then compared simulations with different prefactors f_A > 0.5 at various time steps: t = 40, 88, 176, and 402. These timesteps correspond to the point when the mass of the planet has just reached 2.3 M_J in each run. Figure 9 displays the radial profiles of gas, ice, and vapor separately, as well as the correspondingly estimated C₂H column density profiles. The gas profile and C/O ratio exhibit significant differences among the simulations, but no apparent difference in C₂H production is observed in simulations with f_A > 0.5. Planets with higher accretion rates reach M₀ more quickly, and therefore have shallower gaps. Shallower gas gaps mean less UV penetration, which disfavors C₂H formation. However, the higher luminosity and shallower gap allow for more methane to be released from ice-rich pebbles via sublimation, resulting in a higher C/O ratio. These two effects balance each other, leading to similar C₂H column densities across different accretion rate setups. As a result, inferring a planet’s mass or luminosity from the radial profile of C₂H can be challenging.

For the other low accretion rate runs with f_A < 0.2, the planet mass does not reach 2.3 M_J at the end of the simulation time of 1000 orbits; the results at the end of the simulation are shown in Fig. 10, where the planet masses are 1.66 M_J, 0.98 M_J, and 0.47 M_J in the simulation with f_A = 0.2, 0.1, 0.05, respectively. Once the planet opens a sufficient gap, the accretion rate passes its peak turning point, as shown in Fig. 8b, and the temperature at the Hill radius follows suit, as seen in Fig. 8e. As the temperature around the planet cools down and the gas gap prevents further pebble supply, running the simulation further does not aid in methane release after the accretion rate peak is passed, as shown in Fig. 8f. All of our simulations have passed the peak accretion rate and temperature point.

As discussed above and shown in Fig. 10, if the accretion rate drops below f_A ≤ 0.2, the results will be qualitatively different. Specifically, if the peak luminosity of the planet is too small, the disk will remain too cold to release any methane throughout the gas accretion period (Fig. 8f). Although the increased UV penetration might be helpful with increasing the C₂H density inside the gap, the absolute column density of C₂H will still be insufficient to account for the observed values in the ALMA observations, unless a very deep gap is opened (gap depth Σ_gap/Σ_g,0 ≪ 0.01), as seen in Fig. 10 for the simulation with f_A = 0.2. We acknowledge that this result is sensitive to the empirical formula used in Eq. (22), and high C/O ratios are essential for reproducing the observed value in the literature with more detailed chemical networks (Bosman et al. 2021b,a). The possibility that C₂H formation is mainly due to high UV penetration has been proposed by recent work (e.g., Calahan et al. 2023), which should be particularly interesting in transitional disks and some deep gas gaps. Since the majority of gas gaps found for example in MAPS samples are relatively shallow (gap depth Σ_gap/Σ_g,0 ≳ 0.1; see Table 4 in Zhang et al. 2021), this routine is not the primary focus of this work.

Thus, a dichotomy exists between “hot” and “cold” planets, determined by the maximum accretion rate the planet can reach.

For the disk setup used in this study, a planet accretion rate above ~10⁻⁵M_J yr⁻¹ is necessary to create a prominent C₂H ring whose peak intensity is comparable to those detected in MAPS (Guzmán et al. 2021). This threshold is only weakly dependent on the planet mass M_p, as Eq. (24) shows that $T_{Hill} \propto M_{p}^{1 / 12} {\dot{M}}_{p}^{1 / 4}$ ${T_{{\rm{Hill}}}}\, \propto \,M_{\rm{p}}^{{1 \mathord{\left/ {\vphantom {1 {12}}} \right. \kern-\nulldelimiterspace} {12}}}\dot M_{\rm{p}}^{{1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-\nulldelimiterspace} 4}}$ . In addition, planet mass is always around M_J at the accretion rate peak. Therefore, the maximum accretion rate that a planet can achieve will primarily determine the evolution of volatiles in its vicinity. Conversely, the observed spatial correlation between a C₂H ring and a continuum ring can, in principle, constrain the accretion history of the putative planet.

Fig. 4

Snapshot of fiducial model at t = 75 751 yr (176 orbits), showing different components from the simulation. (a) Surface density of the gas (including vapor). The arrows indicate the velocity of the gas. (b) CH₄ ice surface density. (c) Surface density ratio of CH₄ vapor to total gas. (d) CH₄ sublimated vapor surface density. (e) Disk midplane temperature. (f) C/O ratio in the gas phase. In each panel, the red dashed circle indicates the Hill radius of the planet. A video showing the time evolution can be found in online.

Fig. 5

Location of snowlines at midplane for the molecules of interest in this paper. Based on the midplane gas density at each radius, equating the partial pressure with Eq. (17), the saturation abundance of molecules (solid lines) can be derived. The CO and CH₄ abundance used in this work (8.9 × 10⁻⁴, gray horizontal line) is given as a reference.

Fig. 6

Azimuthally averaged radial profiles of the fiducial model at t = 75 751 yr (176 orbits). (a) Radial profile of total gas density (blue), CH₄ ice (green dashed), and CH₄ vapor (green solid) at 176 orbits. The gas phase C/O ratio (black) is calculated based on the abundance assumed in Table 1. (b) The expected C₂H column density based on Eq. (22). The shaded region indicates the one-magnitude uncertainty of the fitting constant F_C in Eq. (22). We mark the typical C₂H peak density (orange) from the MAPS data constraint as a reference value for illustration.

Fig. 7

Radial profile of gas-phase C/O ratio as a function of time in the fiducial model. The horizontal white dashed line indicates t = 176 orbits. The orange lines represent r₀ ± R_Hill(t), where R_Hill(t) is the Hill radius of the planet at each time.

Fig. 8

Time evolution of different parameters in models with different gas accretion rates (ƒ_A factor in Eq. (19)). (a) Planet mass. The horizontal gray line indicates the planet mass M₀ = 2.3M_J inferred by Liu et al. (2019), which we use as a reference mass. (b) Planet’s gas accretion. (c) Gas surface density at the gap over the unperturbed surface density Σ_g,0. (d) Planet’s accretion luminosity (Eq. (15)). (e) Disk temperature at the Hill radius. (f) Peak gas-phase C/O ratio. The prefactors ƒ_A used in Eq. (19) are listed in the legend. The vertical dashed line indicates 176 orbits, when the planet mass has reached 0.04,0.05,0.12,1.2,2.3,3.8, 6.8 M_J in simulation with ƒ_A = 0.05,0.1,0.2,0.5,1, 2,5, respectively.

Fig. 9

Azimuthally averaged radial profiles of simulations with different accretion rates f_A > 0.5 at M_p = 2.3 M_J. (a) Total gas density (solid), CH₄ vapor density (dashed) at the indicated planet mass over the unperturbed density gas and pebble profile Eq. (12). (b) C/O ratio in the gas (midplane) in our simulations based on the abundance assumed in Table 1. (c) Expected C₂H column density.

4 Discussion

Our findings suggest that the formation of a giant planet could have a substantial impact on the chemistry of protoplanetary disks by locally altering the elemental abundance in the gas and ice phases. In the following we list the main caveats of this work and we broadly discuss some applications and implications of our model.

Fig. 10

Azimuthally averaged radial profiles of simulations with different accretion rates f_A ≤ 0.2 at t = 1000 orbits. The panels are the same as in Fig. 9. The amount of the released carbon is relatively low compared with the high accretion rate runs. High C₂H column density might still be reached when the gap is very deep.

4.1 Limitations

One major limitation of our study is that the hydrodynamical simulation is constrained within a 2D framework. To fully capture the complexity of planet accretion, necessitates the use of high-resolution 3D hydrodynamical simulations. Previous studies have demonstrated that the accretion rate and flow pattern may vary significantly in different model setups for both gas (e.g., Machida et al. 2010; Tanigawa et al. 2012; Fung et al. 2019) and solid material (e.g., Tanigawa et al. 2014; Homma et al. 2020; Szulágyi et al. 2022). In particular, some major differences can be caused by the lack of meridional flow, which can only be probed in 3D simulation (e.g., Szulágyi et al. 2014; Fung & Chiang 2016; Teague et al. 2019). The meridional flow would trigger strong vertical and radial mixing, further complicating the elemental mixing and chemistry. Nevertheless, the computational cost of 3D high-resolution hydrodynamical simulations of gas accretion is substantial, hindering our ability to study the evolution of these systems using current computing resources.

In addition, the lack of a vertical dimension in the simulations limits our ability to account for the stellar UV radiation, and renders the related disk irradiation temperature and UV photochemistry uncertain. As recent works have shown, the planet-opened gap may also become warmer due to higher irradiation from the star (Booth et al. 2021; Pirovano et al. 2022; Broome et al. 2023). The parameterized temperature used in our simulation might only be a reference lower limit. Moreover, the growth and settling of pebbles can affect the vertical transport of ice reservoirs and the depletion of small dust in the upper layer of the disk, which can also potentially elevate the C/O ratio of the gas and increase the UV penetration, as demonstrated by van Clepper et al. (2022). Due to the above reasons, we opted for an empirical C₂H column density estimation in this work (see Sect. 2.4 and Appendix B). The absence of chemical reactions does not impact our conclusion regarding C/O ratio elevation as long as the reaction product remains in the gas phase. To better reproduce the ALMA observational results and understand the astrochemical evolution of protoplanetary disks, future follow-up tests of our model featuring simultaneous hydrodynamic simulations and chemical networks are worthwhile (for some pioneer efforts, see, e.g., Gong et al. 2017, 2023; Ilee et al. 2017; Booth & Ilee 2019; Krijt et al. 2020; Bergner & Ciesla 2021; Hu et al. 2023).

Another caveat is that we do not include condensation of vapor in our simulation. Based on the simulation output in our fiducial model, we calculate the condensation timescales for methane vapor back to particles with size of ~0.5 mm, the pebble size used in the hydrodynamical simulations (see Sect. 2). The condensation timescale reads (e.g., Schoonenberg & Ormel 2017)

$τ_{con} \equiv \frac{Σ_{vapor}}{{\dot{Σ}}_{vapor \to ice}} = \frac{1}{8} \sqrt{\frac{μ_{Z}}{k_{B} T}} \frac{m_{p}}{s_{p}^{2}} \frac{H_{g}}{Σ_{pebble}},$ ${\tau _{{\rm{con}}}} \equiv {{{{\rm{\Sigma }}_{{\rm{vapor}}}}} \over {{{{\rm{\dot \Sigma }}}_{{\rm{vapor}} \to {\rm{ice}}}}}} = {1 \over 8}\sqrt {{{{\mu _Z}} \over {{k_B}T}}} {{{m_{\rm{p}}}} \over {s_{\rm{p}}^2}}{{{H_{\rm{g}}}} \over {{{\rm{\Sigma }}_{{\rm{pebble}}}}}},$ (25)

where µ_Z = 16m_H is the mean molecular weight of methane vapor, m_H is the proton mass, Σ_pebble is the surface density of the pebbles, H_g is the gas scale height, and the typical particle mass follows $m_{p} = (4 / 3) π s_{p}^{3} ρ •$ ${m_{\rm{p}}} = \left( {{4 \mathord{\left/ {\vphantom {4 3}} \right. \kern-\nulldelimiterspace} 3}} \right)\pi s_{\rm{p}}^3\rho \bullet$ . A map of the condensation timescale based on our simulation outputs is shown in Fig. 11a. In regions of high pebble density, outside the gap region, the condensation timescale is shorter than the orbital timescale. In addition, the co-orbital region of the planet, where a moderate amount of pebbles linger, is characterized by short condensation timescales.

Overall, however, due to the depletion of solids, the condensation timescale in the gap region is much longer than the orbital period, with maximum values of ~10⁵ yr. In particular, the C/O ratio is higher in the region with long condensation timescales, as illustrated by the green contour in Fig. 11, which indicates C/O = 1.5. Regions with condensation timescales exceeding 250 orbits (white) are associated with C/O>1.5. Therefore, we conclude that accounting for the condensation of methane in the simulations will not prevent high C/O ratios in the gap, and therefore our major conclusions remain justified.

Similarly, the sublimation timescale reads (see a similar formula for the desorption timescale in Piso et al. 2015)

$τ_{sub} \equiv \frac{Σ_{ice}}{{\dot{Σ}}_{ice \to vapor}} = \frac{1}{8 \sqrt{2 π}} \sqrt{\frac{k_{B} T}{μ_{Z}}} \frac{m_{p}}{s_{p}^{2}} \frac{1}{P_{eq}},$ ${\tau _{{\rm{sub}}}} \equiv {{{{\rm{\Sigma }}_{{\rm{ice}}}}} \over {{{{\rm{\dot \Sigma }}}_{{\rm{ice}} \to {\rm{vapor}}}}}} = {1 \over {8\sqrt {2\pi } }}\sqrt {{{{k_B}T} \over {{\mu _Z}}}} {{{m_{\rm{p}}}} \over {s_{\rm{p}}^2}}{1 \over {{{\rm{P}}_{{\rm{eq}}}}}},$ (26)

which strongly depends on the equilibrium pressure P_eq (Eq. (17)). In Fig. 11b, the sublimation timescale displays a sharp transition around the planet Hill radius, within which the sublimation timescale becomes ≪1 orbit, much shorter than the local spreading time, validating our assumption of instantaneous sublimation. The quick sublimation around the planet and the long condensation timescale in the gap indicate that the methane vapor will first spread through the co-orbital region (which happens on a synodical timescale of several orbital periods) before the vapor has a chance to condense back onto the cold pebble at the midplane.

In addition, the reduced dust optical depth in the gap enables radiation to penetrate deeper into the disk, more effectively warming the gap (e.g., Calahan et al. 2021; Broome et al. 2023). Finally, through strong vertical transport such as meridional flows (Tanigawa et al. 2012; Morbidelli et al. 2014) the released vapor may reach the upper layer of the disk, where temperatures are known to be higher (e.g., D’Alessio et al. 1998). Therefore, the recondensation of CH₄ vapor onto settled pebbles is likely to be ineffective in the gap region. It is worth noting that our results are not sensitive to vapor condensation onto grains of size ≪1 mm. Assuming the condensation timescale is short and the grains are perfectly coupled with the gas¹, the ice re-condensed on the surface of these small grains would nevertheless be exposed to sufficiently high doses of UV photons at the surface of the disk, releasing the carbon back into the gas.

On the other hand, condensation is bound to happen in the colder regions exterior to the continuum gap, where densities are high and temperatures low. As we do not include condensation processes in our simulation, the methane vapor seen beyond 100 au and the corresponding plateau in the C/O ratio are artificial (see, e.g., Fig. 6). Hence, the peak in CH₄ vapor at 100 au seen in Fig. 6a and the corresponding shoulder seen in Fig. 6b will, in reality, be absent. However, we have verified that the conclusions of this study, in particular the peak C₂H, are unaffected by its presence.

Finally, the model neglects any potential migration of the planet (see reviews of Kley & Nelson 2012; Paardekooper et al. 2022). Given the short chemical (Bosman et al. 2021b) and sublimation (desorption) timescale (Piso et al. 2015) of 10⁴ yr in the disk region of our interests, and the short gas accretion windows shown in our simulation, our results should be unaffected by the planet migration happening over a longer timescale. Once the planet opens the gas gap, the migration speed of the planet is slower (e.g., Kanagawa et al. 2018), which further validates our assumption of a slow migration.

Fig. 11

Maps of (a) condensation timescale and (b) sublimation timescale of methane ice in the default simulation at t = 75 751 yr (176 orbits), expressed in units of the orbital period. With the exception of the co-orbital region, the condensation timescale of methane in the gap region is much longer than the synodical timescale ~10 orbits. This allows the gas to spread azimuthally. Since moderate amounts of pebbles may remain inside the co-orbital region, the condensation timescale could still be low there. The C/O ratio is higher in the region with long condensation timescales, as shown in Fig. 4f (see also the green contour that indicates the C/O = 1.5 boundary). Regions with condensation timescales exceeding 250 orbits (white) are associated with C/O>1.5. The red dashed circles indicate the Hill radius of the planet.

4.2 Location of the planet relative to different snowlines

In our scenario, the ability for an accreting planet to elevate the C/O ratio rests on two conditions. First, the planet should be located outside the disk CH₄ snowline, and second, the planet’s accretion luminosity should be able to sublimate CH₄, but not the O-carriers, such as CO₂ or H₂O. To understand how general these conditions are, we plot in Fig. 5 the saturation profiles of said molecules for the disk model adopted in this work. In Fig. 5, the solid curves show the gas-phase abundance (n_i/n(H)) obtained by assuming complete saturation (i.e., when the partial pressure equals the equilibrium pressure of the vapor, P_Z = P_eq), where the equilibrium pressure for different species are taken from Fray & Schmitt (2009) as described in Appendix C. In other words, the curves give the maximum vapor abundance for each species at a given radius or corresponding temperature. It shows that at the initial background temperature of 20 K (73 au), all CO is in the gas phase, whereas CO₂ sublimation remains insignificant even at the point where the highest T_Hill is reached in the fiducial model.

In particular, a local C/O>1 gas enhancement will always appear whenever the planet is located beyond the methane snowline (≈40 au). As was shown in Eq. (24) and Fig. 3d, the peak temperature is determined by the accretion luminosity of the planet, independent of the background disk temperature. Furthermore, the presence of the CO snowline does not influence the outcome of a C/O>1 gas.

On the other hand, accreting planets located between snowlines of species with higher sublimation temperatures may leave a different chemical footprint. For instance, an accreting planet located between the methane and carbon dioxide snowlines may, in turn, lower the C/O ratio in its vicinity by unleashing more oxygen than carbon through the local sublimation of carbon dioxide The basic tenet of this paper, that an accreting planet can locally alter the gas chemistry by selective sublimation of ices, applies generally.

4.3 Exoplanet atmosphere and disk chemistry link

The relationship between the astrochemistry in protoplanetary disks and planetary atmospheres was first highlighted by Öberg et al. (2011). The authors proposed that examining the C/O ratio in planetary atmospheres may allow us to trace the location where the majority of the planetary gas accretion happened as the C/O ratio in both the gas and solid phases vary in step with the increasing radius in the protoplanetary disk. This picture, in which giant planets inherited their atmosphere directly from the gaseous protoplanetary disk, was actively developed in the past decade from both theoretical and observational aspects (e.g., Eistrup et al. 2018; Öberg & Wordsworth 2019; Booth & Ilee 2019; Cridland et al. 2019, 2020; GRAVITY Collaboration 2020; Mollière et al. 2022).

Conversely, our results show that giant planet formation may have significant local impacts on the astrochemistry of the disks. The gap opening and accretion of the planet can locally alter the chemical composition of the disks. Specifically, with a planet located outside the methane snowline, the water and carbon dioxide ice component of the accreted pebbles end up in the atmosphere of the planet, which would therefore become O-rich. In contrast, CH₄ sublimates from pebbles to vapor, which renders the gap region C-rich. Our results indicate that giant planets actively shape their environment not only by opening gaps, but also by changing the chemical inventory of disks, which extends the conventional understanding that the initial chemistry of natal disks is directly inherited by gas giants, as previously suggested by Öberg et al. (2011).

4.4 Chemical rings in MWC 480

Among the five targets of MAPS, MWC 480 stands out as a particularly interesting object for our research. In the continuum, MWC 480 features an annular gap, D76, at 76 au, potentially indicative of a planet (Liu et al. 2019). What makes the D76 gap special is that it is spatially coincident with almost all detected emission lines (14 out of the 18 species investigated in Law et al. 2021, see Fig. 12). Jiang et al. (2022) have found no statistically significant spatial correlations between line emission and continuum in other sources, except MWC 480. Intriguingly, a similar concentration of line-emission rings in a dust gap (cavity) configuration is also found in PDS 70 (Facchini et al. 2021).

For these reasons, we tailored the simulation parameters based on the observational constraints of MWC 480. Our findings can therefore be directly compared to the observations. In Fig. 13, we show the radial profiles of the ALMA continuum emission (Sierra et al. 2021) and the line emission map of C₂H 3−2 (Law et al. 2021). By applying our model to MWC 480, the normalized C₂H column density profile at t = 176 orbits in the fiducial model is shown in orange for comparison. We convolve the profile with a 0.15 arcsec beam, which is used in the MAPS data reduction process. Our results are therefore consistent with the observational data. It indicates that an accreting planet can be responsible for the formation of the observed C₂H ring, in the D76 continuum gap of MWC 480 and possibly other molecular rings. Further investigation to solidify these findings (e.g., by including a chemical network in our simulations) is worthwhile to confirm these preliminary results.

Fig. 12

ALMA Observations of the MWC 480 disk, (a) 257 GHz mm dust continuum emission (Sierra et al. 2021). (b) Line emission moment 0 maps (Law et al. 2021) of six molecular transitions (selected for illustrative clarity; in total over ten lines feature similar ring-like structures). Thf lines are labeled above each panel. The dashed ellipse indicates the D76 continuum gap, and the solid ellipse the B98 continuum ring on the outside in all panels.

Fig. 13

Azimuthally averaged continuum emission profile (gray, normalized to unity at B98) and intensity profile of C₂H 3−2 (blue, normalized to unity at the peak). The normalized pebble surface densities and C₂H column densities at t = 176 orbits in the default run are shown in orange and green for comparison; they are respectively convolved by 0.1 arcsec and 0.15 arcsec beams used in the MAPS observation. The dashed and solid vertical lines coincide with the D76 gap and the B98 ring in the continuum, respectively.

4.5 C₂H rings in other systems

Apart from MWC48O, ring-like C₂H line-emission substructures are identified and studied in a number of other sources (Bergin et al. 2016; Facchini et al. 2021; Miotello et al. 2019; Bergner et al. 2021; Law et al. 2021). In these sources, dust rings and gaps are spatially resolved in the ALMA continuum as well. This raises the interesting question of whether the model proposed in this paper can also explain other systems.

Within the MAPS samples, HD 163296 is the only other system where both C₂H 3−2 and 2−1 lines present a sharp local maximum inside the D49 continuum gap. In addition, a C₂H 3−2 ring peaks inside the D117 continuum gap of IM Lup, where Pinte et al. (2020) suggest a planet candidate depending on its localized kinematics perturbation (kink). And C₂H 1−0 emission radially peaks inside the D100 continuum gap of AS 209, where a planet candidate is hinted at by its kinematic signals (Fedele et al. 2023). Therefore, our model is possibly applicable to these systems as well.

By statistically studying the existence of C₂H emission in disks and the radial profile of their continuum counterpart in 26 disks, van der Marel et al. (2021) suggests a possible correlation between the presence of dust rings outside the CO snowline and the detection of C₂H. It should be noted that MWC 480 is an outlier to the correlation as there is no pebble ring outside its CO snowline, which is at approximately 100 au.

Under the scenario proposed in this work, the accreting planet may sculpt a pebble ring outside its orbit at the outer edge of the gap it opened. This ring is therefore located exterior to the methane snowline. Moreover, it is likely to be situated near or beyond the CO snowline as the scenario proposed in this study works at any location beyond the CH₄ snowline. Therefore, our model, which examines the effect of an accreting planet outside the methane snowline on the chemistry of the disk, could provide a possible explanation for and further complete the tentative correlation identified by van der Marel et al. (2021) between the presence of dust rings outside the CO snowline and the detection of C₂H.

The sample size of disks with detected C₂H remains small, and the available spatial resolution for testing the correlation is limited, which is only possible in the MAPS samples. Future observations, such as the ongoing ALMA Disk-Exoplanet C/Onnection (DECO) large program, will shed more light on this issue by assisting in sample selection. For instance, in the low spatial resolution data of Miotello et al. (2019) two sources, Sz71 (GWLup) and Sz 129, exhibit tentative features in their C₂H 3−2 emission; the emission plausibly radially peaks at the D74 and D64 continuum gaps, respectively. Interestingly, localized kinematic perturbations consistent with the presence of planet candidates have been detected within these gaps (Pinte et al. 2020). Deeper observations with a higher spatial resolution are required to further confirm the link and will be useful for further study.

4.6 Detectablity of the hypothesized planet

Observed rings and gaps in protoplanetary disks are often interpreted as evidence of planet formation (e.g., Pérez et al. 2019; Toci et al. 2020). Correspondingly, a growing number of embedded planets are being claimed within these disks as well (e.g., Hammond et al. 2023; Pinte et al. 2023). Therefore, it is important to determine where the planet is located. A compelling feature of the planet proposed in our model is that, in order to heat up the surrounding area, it is currently or was very recently in the process of accreting.

For instance, MWC 480 provides an apt illustration of this challenge. Although the Strategic Exploration of Exoplanets and Disks with Subaru (SEEDS) survey has succeeded in directly imaging this object (Uyama et al. 2017), its relatively bright Herbig Ae/Be star status, with an H-band magnitude of 6.26, means that detection is still an arduous task. Moreover, identifying any planets in close proximity to the dust ring around MWC 480 is made all the more difficult due to the possibility of confusion with disk features that may be influenced by advanced image analysis techniques (Rameau et al. 2017; Haffert et al. 2019).

Detecting planetary accretion signals provides other direct evidence of planets. If a planet accretes gas from the disk, the Hα line emission originates from the accretion front at the planetary surface. According to the ESO archive, VLT/MUSE (Bacon et al. 2010) observations have been carried out for more than 50 disks in the past 5 yr to search for Hα emission (e.g., Xie et al. 2020); however, so far PDS 70 is the only one in which accreting planets have been identified in MUSE (Haffert et al. 2019). A simple explanation, perhaps, is that many planets in disks are accreting too slowly. Therefore, they do not produce detectable Hα line emission (Marleau et al. 2022). We might have been fortunate with PDS 70 b,c, where the two planets are just above the detection limit. Alternatively, the formation of rings and gaps structure can be independent of the existence of planets (e.g., Bai & Stone 2014; Flock et al. 2015; Okuzumi et al. 2016; Owen & Kollmeier 2019; Ohashi et al. 2021; Jiang & Ormel 2021; Kuznetsova et al. 2022; Tominaga et al. 2022), making the problems more complex.

As discussed in Sect. 3.2, our rough calculation may face challenges in distinguishing different accretion rate histories above 10⁻⁵ M_J yr⁻¹. Additionally, a series of caveats should be acknowledged, as outlined in Sect. 4.1. Despite these limitations, making an educated guess for further investigation is worthwhile. By selecting a mass of 2.3 M_J and the corresponding accretion rate of 3 × 10⁻⁵ M_J yr⁻¹ in our fiducial model, we estimate an Hα flux from the planet of f_Hα = 1.2 × 10⁻¹⁵ erg s⁻¹ cm⁻². This estimate is based on the relationship between Hα emission and accretion rates presented in Fig. 10 of Marleau et al. (2022) and is reasonably above the MUSE detection limit. Our models, therefore, propose a novel testable approach for searching for planets in disks through their chemistry signals. A searching effort for the potential accreting planets in the future might be also been intriguing in other instruments, such as Mag AO-X (Males et al. 2018) SPHERE/ZIMPOL (Cugno et al. 2019), SCExAO/VAMPIRES (Uyama et al. 2020), and HST (Zhou et al. 2021).

5 Conclusions

In this work with the newly developed phase change module (Wang et al. 2023) of the Athena++ code, we describe global multi-fluid hydrodynamic simulations accounting for the sublimation process of the volatile components of (pebble-)accreting planets in a typical protoplanetary disk. The gas accretion of the planet was parameterized with the Machida et al. (2010) formula. Our simulations demonstrate how an accreting planet locally heats its surroundings, creates a gap in the disk, and establishes conditions conducive to C-photochemistry. Starting with a 10 M_⊕ core, the planet grows into a Jovian-mass planet with accretion rate ~10⁻⁵ M_J yr⁻¹ in our fiducial setup. A sufficiently deep gas gap with a gap depth of ~30 % of the unperturbed disk is opened at the moment when the planet reaches 2.3 M_J (Fig. 6). Such a bright planet may locally heat its vicinity, elevating the temperature at its Hill radius to ~45 K. A sketch of our model is shown in Fig. 2. Our main findings are the following:

The enhanced temperature sublimates the major C-carrier ice, methane, surrounding the planet, but is not high enough for water and carbon dioxide ice to sublimate. The released methane locally raises the carbon-to-oxygen ratio inside the gap from unity (CO-dominant gas) to ~2. The opened gap allows more UV penetration toward the planet’s neighborhood. A high C/O ratio and high UV radiation both assist with C-bearing radical (C₂H) formation and the subsequent assembly of organic compounds;
As a result, a C₂H ring is expected to emerge at the dust-gas gap location. The results offer an explanation for the emission of organic molecular lines observed by ALMA in the MAPS large program, where the majority of detected line emission (including C₂H; see Fig. 12) peaks inside the D76 continuum gap of MWC 480. This model may also apply to the C₂H ring peak at the D49 continuum gap of HD 163296;
The outlined model operates when an accreting planet is located beyond the disk CH₄ snowline. It predicts the formation of a C₂H ring at the planet’s location and a pebble ring exterior to it. As the planet may be located beyond the CO snowline as well, the model can potentially account for the correlation between the presence of dust rings outside the CO snowline and the detection of C₂H, as observed by van der Marel et al. (2021) (see Sect. 4.5);
This work demonstrates that the accretion of a giant planet can significantly impact its surrounding astrochemistry environment by locally altering the chemical composition of disks. Due to the different sublimation temperatures, the hot accreting planets can release some volatiles (e.g., methane) while incorporating others (e.g., water), leading to different chemical impacts depending on where the planet is located. This study extends the conventional understanding that the initial chemistry of disks is simply inherited by gas giants, by postulating that planets actively shape their chemical environment of formation;
The finding that a high C/O ratio is the consequence of a gap-opening planet can be inverted to give a rough constraint on the accreting luminosity of the planet. We found that an accretion rate above 10⁻⁵ M_J yr⁻¹ is necessary to ensure that the planet is hot enough to release methane, elevate the C/O ratio, and create a prominent C₂H ring;
For reference, in the fiducial model, the planet mass reaches 2.3 M_J mass when the accreting rate is 3 × 10⁻⁵ M_J yr⁻¹. The Hα flux from the planet is around f_Hα = 1.2 × 10⁻¹⁵ erg s⁻¹ cm⁻², which is potentially detectable via instruments such as VLT/MUSE, SCExAO/VAMPIRES, SPHERE/ZIMPOL, and MagAO-X.

The link between accreting planets and disk astrochemistry proposed in this work is testable with current optical (see Sect. 4.6), IR (JWST), and (sub)millimeter (ALMA) wavelength facilities. Conducting planet-hunting campaigns in disks with pronounced chemical rings and continuum gap correlation might therefore be a new avenue to hunt for and confirm (proto)planets. At the same time, further investigation into the spatially resolved disk chemistry and the characterization of the atmosphere of the protoplanet in confirmed planet-forming (e.g., for PDS 70, Facchini et al. 2021; Cridland et al. 2023) are both of importance to understand how the composition of the planet is built upon and how it influences the disk chemistry in reality because it can possibly reveal the formation history of the planet.

Simultaneously, we plan to extend the model to 3D and apply it to more sources. We also aim to couple our hydro- and ther-modynamical code with a (reduced) on-the-fly chemical network to obtain more quantitative constraints on chemical interactions between the planet and disk. Further investigation will focus on how the accretion of the planet shapes their atmospheric chemistry, how the effect of planet migration influences the chemical signatures, and the impact of multiple planets. Overall, our proposed model opens up exciting new avenues for studying the co-evolution of the disks and protoplanets in planet-forming disks.

Movie

Movie associated with Fig. 4 (Fig4) Access here

Acknowledgements

We thank the anonymous referee for their thoughtful and constructive report that improved the initial manuscript. H.J. and Y.W. would like to thank Pinghui Huang for fruitful discussions about Athena++ code. H.J. thanks Bin Ren and Chen Xie for insightful discussions. H.J. highly appreciates valuable comments from Ilse Cleeves, Alex Cridland, Stefano Facchini, Enrique Macías, Anna Miotello, Giovanni Rosotti, and Claudia Toci. C.W.O. acknowledges support by the National Natural Science Foundation of China (grants no. 12250610189 and 12233004). The authors acknowledge the Tsinghua Astrophysics High-Performance Computing platform at Tsinghua University for providing computational and data storage resources that have contributed to the research results reported within this paper. This work has used Astropy (Astropy Collaboration 2013), Matplotlib (Hunter 2007), Numpy (Harris et al. 2020), Scipy (Virtanen et al. 2020) software packages.

Appendix A Impact of O-carrier volatiles

The most important part of our work relies on the local elevation of the C/O ratio triggered by the accretion luminosity of the planet. In our default setup, we only consider methane in the volatile component. However, oxygen is present in pebbles and in the form of CO₂ and H₂O ice. Consequentially, we expect that oxygen-bearing carriers may also be released from the volatile to the gas phase, which would in contrast lower the local C/O ratio.

Although it is unlikely to happen due to the relatively high sublimation temperature (see Sect. 4.2), we conducted a comparison run by assuming that all icy volatiles in our model are CO₂, while keeping the other setups the same as the fiducial model. The equilibrium pressure (Eq. (17)) of CO₂ follows T_a = 2571 K and P_eq,0 = 2.57 × 10¹² g cm⁻¹ s⁻² (Fray & Schmitt 2009) in the simulation. The outcome is straightforward. The typical sublimation temperature of the CO₂ ice is ≈55 K (orange line in Fig. 3(d), and see Fig. 5), significantly higher than methane. These temperatures can only be reached in the inner regions of the Hill sphere, but never at the Hill radius. Therefore, as T_Hill < 50 K, CO₂ is almost entirely accreted by the protoplanet and cannot be released into the disk.

Since water ice has even higher sublimation temperatures (Fig. 5) H₂O will also stay on the pebbles in ice form, and be lost to the planet before it sublimates. As a result, we can confidently assert that our key findings remain unchanged regardless of the water content of the pebble, or indeed any other O-carrier with a significantly higher sublimation temperature than methane.

Appendix B The C₂H column density formula

In the study by Bosman et al. (2021a) a gas-grain chemical network was employed to calculate the 2D (r − z) chemical makeup of three MAPS disks. However, due to the lack of a vertical dimension in our simulation, it is challenging to directly compare our results with those of Bosman et al. (2021a), where the C₂H abundance is dominated by photochemistry in the upper layer of the disk (see, e.g., their Appendix C). Nevertheless, we outline the key features of the Bosman et al. (2021a) model outputs, and motivate an empirical fit formula that links the C/O ratio to the C₂H abundance, Eq. (22).

Bosman et al. (2021a) concluded that the resolved observations of high but localized C₂H column densities suggest the need for locally higher C/O ratios at C₂H ring locations. Additionally, a correlation is observed between the C/O ratio and the CO abundance, indicating that a higher CO abundance corresponds to a lower C/O ratio, and vice versa.

Here, we provide a description of how the empirical formula Eq. (22) is determined:

$\log_{10} n (C_{2} H) = F_{C} - \log_{10} n (H) + 2 \log_{2} (C / O) .$ ${\log _{10}}n\left( {{{\rm{C}}_2}{\rm{H}}} \right) = {F_C} - {\log _{10}}n\left( {\rm{H}} \right) + 2{\log _2}\left( {{{\rm{C}} \mathord{\left/ {\vphantom {{\rm{C}} {\rm{O}}}} \right. \kern-\nulldelimiterspace} {\rm{O}}}} \right).$

The dependency of n(C₂H) on n(H), and therefore log₁₀ n(CO) = log₁₀ n(H) + 5, is motivated by the anti-correlation observed in Figure 1, Figure 2, and Figure 5 of Bosman et al. (2021a). These figures demonstrate an anti-correlation between CO column densities and C₂H column densities. In particular, Figure 5 shows that the CO column density at the gap location in AS 209 (~ 10¹⁸ cm⁻²) is lower than in HD 163296 and MWC 480 (~ 10²⁰ cm⁻²). However, in Figure 2, the C₂H column density in AS 209 (~ 10¹⁴ cm⁻²) is higher than in HD 163296 and MWC 480 (~ 10¹² cm⁻²) when C/O = 1. Thus, a value of F_C = 37 is obtained as the sum of these contributions.

The dependency of n(C₂H) on log₂(C/O) is estimated based on Figure 2 of Bosman et al. (2021a), which shows that as the C/O ratio increases (or decreases) from 1.0 to 2.0 (or to 0.47), the C₂H column densities rise (or drop) by approximately two orders of magnitude.

Considering the moderate variation in C₂H column density observed in Figure 2 when using different small dust depletion factors, and given that the fitting values were visually identified, we leave an uncertainty of approximately one order of magnitude in F_C.

Appendix C Thermodynamic properties of different molecules

The Clausius–Clapeyron equation is derived from the equality of the chemical potentials between a pure gas-phase compound and its solid phase, which occurs when they are in equilibrium. Taking the perfect gas approximation and neglecting the molar volume of solids, the Clausius–Clapeyron equation can be written as (Fray & Schmitt 2009)

$In (P_{eq}) = In (P_{0}) + \int_{T_{0}}^{T} \frac{Δ H_{eq} (T)}{R T^{2}} d T,$ ${\rm{In}}\left( {{P_{{\rm{eq}}}}} \right) = {\rm{In}}\left( {{P_0}} \right) + \int_{{T_0}}^T {{{{\rm{\Delta }}{H_{{\rm{eq}}}}\left( T \right)} \over {R{T^2}}}} {\rm{d}}T,$ (C.1)

where P₀ is the equilibrium (sublimation) pressure at temperature T₀ and

$\frac{dΔ H_{eq} (T)}{d T} = Δ C_{p} (T) = Δ C_{p, gas} (T) - Δ C_{p, ice} (T)$ ${{{\rm{d\Delta }}{H_{{\rm{eq}}}}\left( T \right)} \over {{\rm{d}}T}} = {\rm{\Delta }}{C_p}\left( T \right) = {\rm{\Delta }}{C_{p,{\rm{gas}}}}\left( T \right) - {\rm{\Delta }}{C_{p,{\rm{ice}}}}\left( T \right)$ (C.2)

is the equilibrium (sublimation) enthalpy which depends on the heat capacities of gas C_p,gas and ice C_p,ice. If the equilibrium enthalpy is supposed to be independent of temperature, then the equilibrium pressure follows the simple relation

$In (P_{eq}) = A_{0} + A_{1} / T,$ ${\rm{In}}\left( {{P_{{\rm{eq}}}}} \right) = {{{A_0} + {A_1}} \mathord{\left/ {\vphantom {{{A_0} + {A_1}} T}} \right. \kern-\nulldelimiterspace} T},$ (C.3)

which can be linked to Eq. (17) by $P_{eq,0} = e^{A_{0}}$ ${P_{{\rm{eq,0}}}} = {e^{{A_0}}}$ and T_a = −A₁. We therefore take the value of A0 and A₁ from Table 5 of Fray & Schmitt (2009) to calculate the values of P_eq,0 and T_a used in our simulation. The original empirical extrapolation formula in Fray & Schmitt (2009) follows ln(P_eq) $= \sum_{i = 0}^{n} A_{i} T^{- i}$ $= \sum\nolimits_{i = 0}^n {{A_i}{T^{ - i}}}$ , and yet we only take the first two terms to keep the formula physically meaningful.

References

Alarcón, F., Bosman, A. D., Bergin, E. A., et al. 2021, ApJS, 257, 8 [CrossRef] [Google Scholar]
Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41 [NASA ADS] [CrossRef] [Google Scholar]
Astropy Collaboration (Robitaille, T. P., et al.) 2013, A & A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bacon, R., Accardo, M., Adjali, L., et al. 2010, SPIE Conf. Ser., 7735, 773508 [Google Scholar]
Bai, X.-N., & Stone, J. M. 2014, ApJ, 796, 31 [Google Scholar]
Bae, J., Teague, R., Andrews, S. M., et al. 2022, ApJ, 934, L20 [NASA ADS] [CrossRef] [Google Scholar]
Barnett, M. N., & Ciesla, F. J. 2022, ApJ, 925, 141 [NASA ADS] [CrossRef] [Google Scholar]
Bast, J. E., Lahuis, F., van Dishoeck, E. F., & Tielens, A. G. G. M. 2013, A & A, 551, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bergez-Casalou, C., Bitsch, B., Pierens, A., Crida, A., & Raymond, S. N. 2020, A & A, 643, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101 [Google Scholar]
Bergner, J. B., & Ciesla, F. 2021, ApJ, 919, 45 [NASA ADS] [CrossRef] [Google Scholar]
Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2019, ApJ, 876, 25 [Google Scholar]
Bergner, J. B., Öberg, K. I., Guzmán, V. V., et al. 2021, ApJS, 257, 11 [NASA ADS] [CrossRef] [Google Scholar]
Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45 [CrossRef] [Google Scholar]
Booth, R. A., & Ilee, J. D. 2019, MNRAS, 487, 3998 [CrossRef] [Google Scholar]
Booth, A. S., Walsh, C., Terwisscha van Scheltinga, J., et al. 2021, Nat. Astron., 5, 684 [NASA ADS] [CrossRef] [Google Scholar]
Booth, A. S., Ilee, J. D., Walsh, C., et al. 2023, A & A, 669, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bosman, A. D., Alarcón, F., Bergin, E. A., et al. 2021a, ApJS, 257, 7 [NASA ADS] [CrossRef] [Google Scholar]
Bosman, A. D., Alarcón, F., Zhang, K., & Bergin, E. A. 2021b, ApJ, 910, 3 [NASA ADS] [CrossRef] [Google Scholar]
Broome, M., Kama, M., Booth, R., & Shorttle, O. 2023, MNRAS, 522, 3378 [NASA ADS] [CrossRef] [Google Scholar]
Bruderer, S., van Dishoeck, E. F., Doty, S. D., & Herczeg, G. J. 2012, A & A, 541, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Calahan, J. K., Bergin, E., Zhang, K., et al. 2021, ApJ, 908, 8 [NASA ADS] [CrossRef] [Google Scholar]
Calahan, J. K., Bergin, E. A., Bosman, A. D., et al. 2023, Nat. Astron., 7, 49 [Google Scholar]
Cleeves, L. I., Bergin, E. A., & Harries, T. J. 2015, ApJ, 807, 2 [NASA ADS] [CrossRef] [Google Scholar]
Crida, A., & Bitsch, B. 2017, Icarus, 285, 145 [Google Scholar]
Cridland, A. J., van Dishoeck, E. F., Alessi, M., & Pudritz, R. E. 2019, A & A, 632, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cridland, A. J., van Dishoeck, E. F., Alessi, M., & Pudritz, R. E. 2020, A & A, 642, A229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cridland, A. J., Facchini, S., van Dishoeck, E. F., & Benisty, M. 2023, A & A, 674, A211 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cugno, G., Quanz, S. P., Hunziker, S., et al. 2019, A & A, 622, A156 [CrossRef] [EDP Sciences] [Google Scholar]
Currie, T., Lawson, K., Schneider, G., et al. 2022, Nat. Astron., 6, 751 [NASA ADS] [CrossRef] [Google Scholar]
D'Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411 [Google Scholar]
Dong, R., Rafikov, R. R., Stone, J. M., & Petrovich, C. 2011, ApJ, 741, 56 [NASA ADS] [CrossRef] [Google Scholar]
Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015, ApJ, 809, L5 [Google Scholar]
Drozdovskaya, M. N., Walsh, C., van Dishoeck, E. F., et al. 2016, MNRAS, 462, 977 [NASA ADS] [CrossRef] [Google Scholar]
Du, F., & Bergin, E. A. 2014, ApJ, 792, 2 [NASA ADS] [CrossRef] [Google Scholar]
Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2016, A & A, 595, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2018, A & A, 613, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Facchini, S., Teague, R., Bae, J., et al. 2021, AJ, 162, 99 [NASA ADS] [CrossRef] [Google Scholar]
Fedele, D., Bollati, F., & Lodato, G. 2023, A & A, 672, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A & A, 574, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Follette, K. B., Close, L. M., Males, J. R., et al. 2023, AJ, 165, 225 [NASA ADS] [CrossRef] [Google Scholar]
Fray, N., & Schmitt, B. 2009, Planet. Space Sci., 57, 2053 [Google Scholar]
Fung, J., & Chiang, E. 2016, ApJ, 832, 105 [NASA ADS] [CrossRef] [Google Scholar]
Fung, J., Zhu, Z., & Chiang, E. 2019, ApJ, 887, 152 [NASA ADS] [CrossRef] [Google Scholar]
Gong, M., Ostriker, E. C., & Wolfire, M. G. 2017, ApJ, 843, 38 [NASA ADS] [CrossRef] [Google Scholar]
Gong, M., Ho, K.-W., Stone, J. M., et al. 2023, arXiv e-prints [arXiv:2385.84965] [Google Scholar]
GRAVITY Collaboration (Nowak, M., et al.) 2020, A & A, 633, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Guzmán, V. V., Pety, J., Goicoechea, J. R., et al. 2015, ApJ, 800, L33 [CrossRef] [Google Scholar]
Guzmán, V. V., Bergner, J. B., Law, C. J., et al. 2021, ApJS, 257, 6 [CrossRef] [Google Scholar]
Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749 [Google Scholar]
Hammond, I., Christiaens, V., Price, D. J., et al. 2023, MNRAS, 522, L51 [NASA ADS] [CrossRef] [Google Scholar]
Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
Homma, T., Ohtsuki, K., Maeda, N., et al. 2020, ApJ, 903, 98 [NASA ADS] [CrossRef] [Google Scholar]
Hori, Y., & Ikoma, M. 2011, MNRAS, 416, 1419 [NASA ADS] [CrossRef] [Google Scholar]
Hu, X., Li, Z.-Y., Wang, L., Zhu, Z., & Bae, J. 2023, MNRAS, 523, 4883 [NASA ADS] [CrossRef] [Google Scholar]
Huang, P., & Bai, X.-N. 2022, ApJS, 262, 11 [NASA ADS] [CrossRef] [Google Scholar]
Huélamo, N., Chauvin, G., Mendigutía, I., et al. 2022, A & A, 668, A138 [CrossRef] [EDP Sciences] [Google Scholar]
Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
Ilee, J. D., Forgan, D. H., Evans, M. G., et al. 2017, MNRAS, 472, 189 [NASA ADS] [CrossRef] [Google Scholar]
Izquierdo, A. F., Facchini, S., Rosotti, G. P., van Dishoeck, E. F., & Testi, L. 2022, ApJ, 928, 2 [NASA ADS] [CrossRef] [Google Scholar]
Jiang, H., & Ormel, C. W. 2021, MNRAS, 505, 1162 [NASA ADS] [CrossRef] [Google Scholar]
Jiang, H., Zhu, W., & Ormel, C. W. 2022, ApJ, 924, L31 [NASA ADS] [CrossRef] [Google Scholar]
Johansen, A., Ronnet, T., Bizzarro, M., et al. 2021, Sci. Adv., 7, eabc0444 [NASA ADS] [CrossRef] [Google Scholar]
Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A & A, 592, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Kanagawa, K. D., Tanaka, H., Muto, T., & Tanigawa, T. 2017, PASJ, 69, 97 [NASA ADS] [CrossRef] [Google Scholar]
Kanagawa, K. D., Tanaka, H., & Szuszkiewicz, E. 2018, ApJ, 861, 140 [Google Scholar]
Keppler, M., Benisty, M., Müller, A., et al. 2018, A & A, 617, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Kley, W., & Nelson, R. P. 2012, ARA & A, 50, 211 [NASA ADS] [CrossRef] [Google Scholar]
Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134 [NASA ADS] [CrossRef] [Google Scholar]
Kuznetsova, A., Bae, J., Hartmann, L., & Low, M.-M. M. 2022, ApJ, 928, 92 [NASA ADS] [CrossRef] [Google Scholar]
Law, C. J., Loomis, R. A., Teague, R., et al. 2021, ApJS, 257, 3 [NASA ADS] [CrossRef] [Google Scholar]
Liu, Y., Dipierro, G., Ragusa, E., et al. 2019, A & A, 622, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17 [Google Scholar]
Machida, M. N., Kokubo, E., Inutsuka, S.-I., & Matsumoto, T. 2010, MNRAS, 405, 1227 [Google Scholar]
Males, J. R., Close, L. M., Miller, K., et al. 2018, SPIE Conf. Ser., 10703, 1070309 [NASA ADS] [Google Scholar]
Marleau, G. D., Aoyama, Y., Kuiper, R., et al. 2022, A & A, 657, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Miotello, A., Facchini, S., van Dishoeck, E. F., et al. 2019, A & A, 631, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Mollière, P., Molyarova, T., Bitsch, B., et al. 2022, ApJ, 934, 74 [CrossRef] [Google Scholar]
Morbidelli, A., Szulágyi, J., Crida, A., et al. 2014, Icarus, 232, 266 [Google Scholar]
Mordasini, C., Alibert, Y., Georgy, C., et al. 2012, A & A, 547, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Müller, T. W. A., Kley, W., & Meru, F. 2012, A & A, 541, A123 [CrossRef] [EDP Sciences] [Google Scholar]
Öberg, K. I., & Bergin, E. A. 2021, Phys. Rep., 893, 1 [Google Scholar]
Öberg, K. I., & Wordsworth, R. 2019, AJ, 158, 194 [Google Scholar]
Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16 [Google Scholar]
Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1 [CrossRef] [Google Scholar]
Ohashi, S., Kobayashi, H., Nakatani, R., et al. 2021, ApJ, 907, 80 [NASA ADS] [CrossRef] [Google Scholar]
Okuzumi, S., Momose, M., Sirono, S.-I., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82 [NASA ADS] [CrossRef] [Google Scholar]
Ormel, C. W., Vazan, A., & Brouwers, M. G. 2021, A & A, 647, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Owen, J. E., & Kollmeier, J. A. 2019, MNRAS, 487, 3702 [Google Scholar]
Paardekooper, S.-J., Dong, R., Duffell, P., et al. 2022, arXiv e-prints [arXiv:2203.09595] [Google Scholar]
Pérez, S., Casassus, S., Baruteau, C., et al. 2019, AJ, 158, 15 [Google Scholar]
Perri, F., & Cameron, A. G. W. 1974, Icarus, 22, 416 [Google Scholar]
Pinte, C., van der Plas, G., Ménard, F., et al. 2019, Nat. Astron., 3, 1109 [NASA ADS] [CrossRef] [Google Scholar]
Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJ, 890, L9 [CrossRef] [Google Scholar]
Pinte, C., Hammond, I., Price, D. J., et al. 2023, MNRAS, 526, L41 [NASA ADS] [CrossRef] [Google Scholar]
Pirovano, L. M., Fedele, D., van Dishoeck, E. F., et al. 2022, A & A, 665, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Piso, A.-M. A., Öberg, K. I., Birnstiel, T., & Murray-Clay, R. A. 2015, ApJ, 815, 109 [Google Scholar]
Rafikov, R. R. 2006, ApJ, 648, 666 [Google Scholar]
Rameau, J., Follette, K. B., Pueyo, L., et al. 2017, AJ, 153, 244 [Google Scholar]
Robert, C. M. T., Crida, A., Lega, E., Méheut, H., & Morbidelli, A. 2018, A & A, 617, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Rosotti, G. P. 2023, New Astron. Rev., 96, 101674 [NASA ADS] [CrossRef] [Google Scholar]
Rosotti, G. P., Juhasz, A., Booth, R. A., & Clarke, C. J. 2016, MNRAS, 459, 2790 [Google Scholar]
Schneider, A. D., & Bitsch, B. 2021, A & A, 654, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schoonenberg, D., & Ormel, C. W. 2017, A & A, 602, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Shakura, N. I., & Sunyaev, R. A. 1973, A & A, 24, 337 [NASA ADS] [Google Scholar]
Sierra, A., Pérez, L. M., Zhang, K., et al. 2021, ApJS, 257, 14 [NASA ADS] [CrossRef] [Google Scholar]
Speedie, J., & Dong, R. 2022, ApJ, 940, L43 [NASA ADS] [CrossRef] [Google Scholar]
Speedie, J., Booth, R. A., & Dong, R. 2022, ApJ, 930, 40 [NASA ADS] [CrossRef] [Google Scholar]
Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65 [CrossRef] [Google Scholar]
Szulágyi, J., Binkert, F., & Surville, C. 2022, ApJ, 924, 1 [CrossRef] [Google Scholar]
Tanigawa, T., & Tanaka, H. 2016, ApJ, 823, 48 [NASA ADS] [CrossRef] [Google Scholar]
Tanigawa, T., Ohtsuki, K., & Machida, M. N. 2012, ApJ, 747, 47 [NASA ADS] [CrossRef] [Google Scholar]
Tanigawa, T., Maruta, A., & Machida, M. N. 2014, ApJ, 784, 109 [NASA ADS] [CrossRef] [Google Scholar]
Teague, R., Bae, J., & Bergin, E. A. 2019, Nature, 574, 378 [NASA ADS] [CrossRef] [Google Scholar]
Teague, R., Bae, J., Aikawa, Y., et al. 2021, ApJS, 257, 18 [NASA ADS] [CrossRef] [Google Scholar]
Toci, C., Lodato, G., Fedele, D., Testi, L., & Pinte, C. 2020, ApJ, 888, L4 [Google Scholar]
Tominaga, R. T., Tanaka, H., Kobayashi, H., & Inutsuka, S.-I. 2022, ApJ, 940, 152 [NASA ADS] [CrossRef] [Google Scholar]
Uyama, T., Hashimoto, J., Kuzuhara, M., et al. 2017, AJ, 153, 106 [Google Scholar]
Uyama, T., Norris, B., Jovanovic, N., et al. 2020, J. Astron. Telescopes Instrum. Syst., 6, 045004 [NASA ADS] [Google Scholar]
van Clepper, E., Bergner, J. B., Bosman, A. D., Bergin, E., & Ciesla, F. J. 2022, ApJ, 927, 206 [NASA ADS] [CrossRef] [Google Scholar]
van der Marel, N., Dong, R., di Francesco, J., Williams, J. P., & Tobin, J. 2019, ApJ, 872, 112 [Google Scholar]
van der Marel, N., Bosman, A. D., Krijt, S., Mulders, G. D., & Bergner, J. B. 2021, A & A, 653, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
Wang, Y., Ormel, C. W., Huang, P., & Kuiper, R. 2023, MNRAS, 523, 6186 [NASA ADS] [CrossRef] [Google Scholar]
Xie, C., Haffert, S. Y., de Boer, J., et al. 2020, A & A, 644, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]
Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47 [Google Scholar]
Zhang, K., Booth, A. S., Law, C. J., et al. 2021, ApJS, 257, 5 [NASA ADS] [CrossRef] [Google Scholar]
Zhou, Y., Bowler, B. P., Wagner, K. R., et al. 2021, AJ, 161, 244 [NASA ADS] [CrossRef] [Google Scholar]
Zhou, Y., Sanghi, A., Bowler, B. P., et al. 2022, ApJ, 934, L13 [NASA ADS] [CrossRef] [Google Scholar]
Zhu, Z., Nelson, R. P., Hartmann, L., Espaillat, C., & Calvet, N. 2011, ApJ, 729, 47 [NASA ADS] [CrossRef] [Google Scholar]

¹

Dust coagulation may be ineffective inside the gap due to the low density. However, small particles can still be aerodynamically large, preventing vertical diffusion and allowing for sweep-up by pebbles. A detailed model of the balance between vertical transport and the thermal process is beyond the scope of this work (but see, e.g., van Clepper et al. 2022).

All Tables

Table 1

Abundances n(X)/n(H) for different molecules X of interest.

In the text

All Figures

	Fig. 1 Sketch of 2D mesh grid used in our simulations. The circle represents the location of the planet (at r₀ = 73 au) and the star is the central star. The radial domain ranges from 30 au to 150 au. The static mesh refinement level is 2, which means the smallest mesh is four times smaller than the unrefined mesh at the same location.
In the text

Fig. 2

Schematic showing how a gap-opening planet triggers the C₂H ring associated with the continuum gap. (a) Before runaway growth, a small planet is embedded inside the gaseous disk and carbon-rich ice is frozen out on pebbles. (b) During runaway gas accretion, the planet grows rapidly, consumes the pebbles in its vicinity, and opens a gas gap in the disk. Since the accreting planet is hot, carbon-rich volatile ice (e.g., methane) with lower sublimation temperature is readily liberated from the ice. On the other hand, the sublimation of H₂O or CO₂, which are the major oxygen carriers, can only be reached in the immediate vicinity of the planet (<R_Hill), at which point the material is lost to the planet. The discrepancy in sublimation temperature renders the gas C/O>1. In addition, stellar UV photons penetrate the gas gap. The combination of C/O>1 and high UV fluxes trigger photochemistry reactions to form the C₂H ring. (c) After the gap opens, the supply of sublimated pebbles terminates due to the gap opening. The C/O ratio is smoothed out across the disk through radial mixing by diffusion.

In the text

Fig. 3

Time evolution of different parameters in the fiducial model. (a) Planet mass. The horizontal gray line indicates the planet mass M₀ = 2.3M_J inferred by Liu et al. (2019), which we use as a reference mass. (b) Planet’s gas accretion. (c) Planet’s accretion luminosity (Eq. (15)). (d) Disk temperature at the Hill radius. The sublimation temperature of CH₄ (green) and CO₂ (orange) are indicated by horizontal lines for reference (see discussion in Sect. 4.2). In each panel, the vertical dashed lines indicate the time when the mass of the planet reaches M₀.

In the text

Fig. 4

Snapshot of fiducial model at t = 75 751 yr (176 orbits), showing different components from the simulation. (a) Surface density of the gas (including vapor). The arrows indicate the velocity of the gas. (b) CH₄ ice surface density. (c) Surface density ratio of CH₄ vapor to total gas. (d) CH₄ sublimated vapor surface density. (e) Disk midplane temperature. (f) C/O ratio in the gas phase. In each panel, the red dashed circle indicates the Hill radius of the planet. A video showing the time evolution can be found in online.

In the text

	Fig. 5 Location of snowlines at midplane for the molecules of interest in this paper. Based on the midplane gas density at each radius, equating the partial pressure with Eq. (17), the saturation abundance of molecules (solid lines) can be derived. The CO and CH₄ abundance used in this work (8.9 × 10⁻⁴, gray horizontal line) is given as a reference.
In the text

Fig. 6

Azimuthally averaged radial profiles of the fiducial model at t = 75 751 yr (176 orbits). (a) Radial profile of total gas density (blue), CH₄ ice (green dashed), and CH₄ vapor (green solid) at 176 orbits. The gas phase C/O ratio (black) is calculated based on the abundance assumed in Table 1. (b) The expected C₂H column density based on Eq. (22). The shaded region indicates the one-magnitude uncertainty of the fitting constant F_C in Eq. (22). We mark the typical C₂H peak density (orange) from the MAPS data constraint as a reference value for illustration.

In the text

	Fig. 7 Radial profile of gas-phase C/O ratio as a function of time in the fiducial model. The horizontal white dashed line indicates t = 176 orbits. The orange lines represent r₀ ± R_Hill(t), where R_Hill(t) is the Hill radius of the planet at each time.
In the text

Fig. 8

Time evolution of different parameters in models with different gas accretion rates (ƒ_A factor in Eq. (19)). (a) Planet mass. The horizontal gray line indicates the planet mass M₀ = 2.3M_J inferred by Liu et al. (2019), which we use as a reference mass. (b) Planet’s gas accretion. (c) Gas surface density at the gap over the unperturbed surface density Σ_g,0. (d) Planet’s accretion luminosity (Eq. (15)). (e) Disk temperature at the Hill radius. (f) Peak gas-phase C/O ratio. The prefactors ƒ_A used in Eq. (19) are listed in the legend. The vertical dashed line indicates 176 orbits, when the planet mass has reached 0.04,0.05,0.12,1.2,2.3,3.8, 6.8 M_J in simulation with ƒ_A = 0.05,0.1,0.2,0.5,1, 2,5, respectively.

In the text

Fig. 9

Azimuthally averaged radial profiles of simulations with different accretion rates f_A > 0.5 at M_p = 2.3 M_J. (a) Total gas density (solid), CH₄ vapor density (dashed) at the indicated planet mass over the unperturbed density gas and pebble profile Eq. (12). (b) C/O ratio in the gas (midplane) in our simulations based on the abundance assumed in Table 1. (c) Expected C₂H column density.

In the text

	Fig. 10 Azimuthally averaged radial profiles of simulations with different accretion rates f_A ≤ 0.2 at t = 1000 orbits. The panels are the same as in Fig. 9. The amount of the released carbon is relatively low compared with the high accretion rate runs. High C₂H column density might still be reached when the gap is very deep.
In the text

Fig. 11

Maps of (a) condensation timescale and (b) sublimation timescale of methane ice in the default simulation at t = 75 751 yr (176 orbits), expressed in units of the orbital period. With the exception of the co-orbital region, the condensation timescale of methane in the gap region is much longer than the synodical timescale ~10 orbits. This allows the gas to spread azimuthally. Since moderate amounts of pebbles may remain inside the co-orbital region, the condensation timescale could still be low there. The C/O ratio is higher in the region with long condensation timescales, as shown in Fig. 4f (see also the green contour that indicates the C/O = 1.5 boundary). Regions with condensation timescales exceeding 250 orbits (white) are associated with C/O>1.5. The red dashed circles indicate the Hill radius of the planet.

In the text

Fig. 12

ALMA Observations of the MWC 480 disk, (a) 257 GHz mm dust continuum emission (Sierra et al. 2021). (b) Line emission moment 0 maps (Law et al. 2021) of six molecular transitions (selected for illustrative clarity; in total over ten lines feature similar ring-like structures). Thf lines are labeled above each panel. The dashed ellipse indicates the D76 continuum gap, and the solid ellipse the B98 continuum ring on the outside in all panels.

In the text

Fig. 13

Azimuthally averaged continuum emission profile (gray, normalized to unity at B98) and intensity profile of C₂H 3−2 (blue, normalized to unity at the peak). The normalized pebble surface densities and C₂H column densities at t = 176 orbits in the default run are shown in orange and green for comparison; they are respectively convolved by 0.1 arcsec and 0.15 arcsec beams used in the MAPS observation. The dashed and solid vertical lines coincide with the D76 gap and the B98 ring in the continuum, respectively.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.

[1] Alarcón, F., Bosman, A. D., Bergin, E. A., et al. 2021, ApJS, 257, 8 [CrossRef] [Google Scholar]

[2] Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41 [NASA ADS] [CrossRef] [Google Scholar]

[3] Astropy Collaboration (Robitaille, T. P., et al.) 2013, A & A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[4] Bacon, R., Accardo, M., Adjali, L., et al. 2010, SPIE Conf. Ser., 7735, 773508 [Google Scholar]

[5] Bai, X.-N., & Stone, J. M. 2014, ApJ, 796, 31 [Google Scholar]

[6] Bae, J., Teague, R., Andrews, S. M., et al. 2022, ApJ, 934, L20 [NASA ADS] [CrossRef] [Google Scholar]

[7] Barnett, M. N., & Ciesla, F. J. 2022, ApJ, 925, 141 [NASA ADS] [CrossRef] [Google Scholar]

[8] Bast, J. E., Lahuis, F., van Dishoeck, E. F., & Tielens, A. G. G. M. 2013, A & A, 551, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[9] Bergez-Casalou, C., Bitsch, B., Pierens, A., Crida, A., & Raymond, S. N. 2020, A & A, 643, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[10] Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101 [Google Scholar]

[11] Bergner, J. B., & Ciesla, F. 2021, ApJ, 919, 45 [NASA ADS] [CrossRef] [Google Scholar]

[12] Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2019, ApJ, 876, 25 [Google Scholar]

[13] Bergner, J. B., Öberg, K. I., Guzmán, V. V., et al. 2021, ApJS, 257, 11 [NASA ADS] [CrossRef] [Google Scholar]

[14] Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45 [CrossRef] [Google Scholar]

[15] Booth, R. A., & Ilee, J. D. 2019, MNRAS, 487, 3998 [CrossRef] [Google Scholar]

[16] Booth, A. S., Walsh, C., Terwisscha van Scheltinga, J., et al. 2021, Nat. Astron., 5, 684 [NASA ADS] [CrossRef] [Google Scholar]

[17] Booth, A. S., Ilee, J. D., Walsh, C., et al. 2023, A & A, 669, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[18] Bosman, A. D., Alarcón, F., Bergin, E. A., et al. 2021a, ApJS, 257, 7 [NASA ADS] [CrossRef] [Google Scholar]

[19] Bosman, A. D., Alarcón, F., Zhang, K., & Bergin, E. A. 2021b, ApJ, 910, 3 [NASA ADS] [CrossRef] [Google Scholar]

[20] Broome, M., Kama, M., Booth, R., & Shorttle, O. 2023, MNRAS, 522, 3378 [NASA ADS] [CrossRef] [Google Scholar]

[21] Bruderer, S., van Dishoeck, E. F., Doty, S. D., & Herczeg, G. J. 2012, A & A, 541, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[22] Calahan, J. K., Bergin, E., Zhang, K., et al. 2021, ApJ, 908, 8 [NASA ADS] [CrossRef] [Google Scholar]

[23] Calahan, J. K., Bergin, E. A., Bosman, A. D., et al. 2023, Nat. Astron., 7, 49 [Google Scholar]

[24] Cleeves, L. I., Bergin, E. A., & Harries, T. J. 2015, ApJ, 807, 2 [NASA ADS] [CrossRef] [Google Scholar]

[25] Crida, A., & Bitsch, B. 2017, Icarus, 285, 145 [Google Scholar]

[26] Cridland, A. J., van Dishoeck, E. F., Alessi, M., & Pudritz, R. E. 2019, A & A, 632, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[27] Cridland, A. J., van Dishoeck, E. F., Alessi, M., & Pudritz, R. E. 2020, A & A, 642, A229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[28] Cridland, A. J., Facchini, S., van Dishoeck, E. F., & Benisty, M. 2023, A & A, 674, A211 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[29] Cugno, G., Quanz, S. P., Hunziker, S., et al. 2019, A & A, 622, A156 [CrossRef] [EDP Sciences] [Google Scholar]

[30] Currie, T., Lawson, K., Schneider, G., et al. 2022, Nat. Astron., 6, 751 [NASA ADS] [CrossRef] [Google Scholar]

[31] D'Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411 [Google Scholar]

[32] Dong, R., Rafikov, R. R., Stone, J. M., & Petrovich, C. 2011, ApJ, 741, 56 [NASA ADS] [CrossRef] [Google Scholar]

[33] Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015, ApJ, 809, L5 [Google Scholar]

[34] Drozdovskaya, M. N., Walsh, C., van Dishoeck, E. F., et al. 2016, MNRAS, 462, 977 [NASA ADS] [CrossRef] [Google Scholar]

[35] Du, F., & Bergin, E. A. 2014, ApJ, 792, 2 [NASA ADS] [CrossRef] [Google Scholar]

[36] Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2016, A & A, 595, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[37] Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2018, A & A, 613, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[38] Facchini, S., Teague, R., Bae, J., et al. 2021, AJ, 162, 99 [NASA ADS] [CrossRef] [Google Scholar]

[39] Fedele, D., Bollati, F., & Lodato, G. 2023, A & A, 672, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[40] Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A & A, 574, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[41] Follette, K. B., Close, L. M., Males, J. R., et al. 2023, AJ, 165, 225 [NASA ADS] [CrossRef] [Google Scholar]

[42] Fray, N., & Schmitt, B. 2009, Planet. Space Sci., 57, 2053 [Google Scholar]

[43] Fung, J., & Chiang, E. 2016, ApJ, 832, 105 [NASA ADS] [CrossRef] [Google Scholar]

[44] Fung, J., Zhu, Z., & Chiang, E. 2019, ApJ, 887, 152 [NASA ADS] [CrossRef] [Google Scholar]

[45] Gong, M., Ostriker, E. C., & Wolfire, M. G. 2017, ApJ, 843, 38 [NASA ADS] [CrossRef] [Google Scholar]

[46] Gong, M., Ho, K.-W., Stone, J. M., et al. 2023, arXiv e-prints [arXiv:2385.84965] [Google Scholar]

[47] GRAVITY Collaboration (Nowak, M., et al.) 2020, A & A, 633, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[48] Guzmán, V. V., Pety, J., Goicoechea, J. R., et al. 2015, ApJ, 800, L33 [CrossRef] [Google Scholar]

[49] Guzmán, V. V., Bergner, J. B., Law, C. J., et al. 2021, ApJS, 257, 6 [CrossRef] [Google Scholar]

[50] Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749 [Google Scholar]

[51] Hammond, I., Christiaens, V., Price, D. J., et al. 2023, MNRAS, 522, L51 [NASA ADS] [CrossRef] [Google Scholar]

[52] Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]

[53] Homma, T., Ohtsuki, K., Maeda, N., et al. 2020, ApJ, 903, 98 [NASA ADS] [CrossRef] [Google Scholar]

[54] Hori, Y., & Ikoma, M. 2011, MNRAS, 416, 1419 [NASA ADS] [CrossRef] [Google Scholar]

[55] Hu, X., Li, Z.-Y., Wang, L., Zhu, Z., & Bae, J. 2023, MNRAS, 523, 4883 [NASA ADS] [CrossRef] [Google Scholar]

[56] Huang, P., & Bai, X.-N. 2022, ApJS, 262, 11 [NASA ADS] [CrossRef] [Google Scholar]

[57] Huélamo, N., Chauvin, G., Mendigutía, I., et al. 2022, A & A, 668, A138 [CrossRef] [EDP Sciences] [Google Scholar]

[58] Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]

[59] Ilee, J. D., Forgan, D. H., Evans, M. G., et al. 2017, MNRAS, 472, 189 [NASA ADS] [CrossRef] [Google Scholar]

[60] Izquierdo, A. F., Facchini, S., Rosotti, G. P., van Dishoeck, E. F., & Testi, L. 2022, ApJ, 928, 2 [NASA ADS] [CrossRef] [Google Scholar]

[61] Jiang, H., & Ormel, C. W. 2021, MNRAS, 505, 1162 [NASA ADS] [CrossRef] [Google Scholar]

[62] Jiang, H., Zhu, W., & Ormel, C. W. 2022, ApJ, 924, L31 [NASA ADS] [CrossRef] [Google Scholar]

[63] Johansen, A., Ronnet, T., Bizzarro, M., et al. 2021, Sci. Adv., 7, eabc0444 [NASA ADS] [CrossRef] [Google Scholar]

[64] Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A & A, 592, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[65] Kanagawa, K. D., Tanaka, H., Muto, T., & Tanigawa, T. 2017, PASJ, 69, 97 [NASA ADS] [CrossRef] [Google Scholar]

[66] Kanagawa, K. D., Tanaka, H., & Szuszkiewicz, E. 2018, ApJ, 861, 140 [Google Scholar]

[67] Keppler, M., Benisty, M., Müller, A., et al. 2018, A & A, 617, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[68] Kley, W., & Nelson, R. P. 2012, ARA & A, 50, 211 [NASA ADS] [CrossRef] [Google Scholar]

[69] Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134 [NASA ADS] [CrossRef] [Google Scholar]

[70] Kuznetsova, A., Bae, J., Hartmann, L., & Low, M.-M. M. 2022, ApJ, 928, 92 [NASA ADS] [CrossRef] [Google Scholar]

[71] Law, C. J., Loomis, R. A., Teague, R., et al. 2021, ApJS, 257, 3 [NASA ADS] [CrossRef] [Google Scholar]

[72] Liu, Y., Dipierro, G., Ragusa, E., et al. 2019, A & A, 622, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[73] Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17 [Google Scholar]

[74] Machida, M. N., Kokubo, E., Inutsuka, S.-I., & Matsumoto, T. 2010, MNRAS, 405, 1227 [Google Scholar]

[75] Males, J. R., Close, L. M., Miller, K., et al. 2018, SPIE Conf. Ser., 10703, 1070309 [NASA ADS] [Google Scholar]

[76] Marleau, G. D., Aoyama, Y., Kuiper, R., et al. 2022, A & A, 657, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[77] Miotello, A., Facchini, S., van Dishoeck, E. F., et al. 2019, A & A, 631, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[78] Mollière, P., Molyarova, T., Bitsch, B., et al. 2022, ApJ, 934, 74 [CrossRef] [Google Scholar]

[79] Morbidelli, A., Szulágyi, J., Crida, A., et al. 2014, Icarus, 232, 266 [Google Scholar]

[80] Mordasini, C., Alibert, Y., Georgy, C., et al. 2012, A & A, 547, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[81] Müller, T. W. A., Kley, W., & Meru, F. 2012, A & A, 541, A123 [CrossRef] [EDP Sciences] [Google Scholar]

[82] Öberg, K. I., & Bergin, E. A. 2021, Phys. Rep., 893, 1 [Google Scholar]

[83] Öberg, K. I., & Wordsworth, R. 2019, AJ, 158, 194 [Google Scholar]

[84] Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16 [Google Scholar]

[85] Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1 [CrossRef] [Google Scholar]

[86] Ohashi, S., Kobayashi, H., Nakatani, R., et al. 2021, ApJ, 907, 80 [NASA ADS] [CrossRef] [Google Scholar]

[87] Okuzumi, S., Momose, M., Sirono, S.-I., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82 [NASA ADS] [CrossRef] [Google Scholar]

[88] Ormel, C. W., Vazan, A., & Brouwers, M. G. 2021, A & A, 647, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[89] Owen, J. E., & Kollmeier, J. A. 2019, MNRAS, 487, 3702 [Google Scholar]

[90] Paardekooper, S.-J., Dong, R., Duffell, P., et al. 2022, arXiv e-prints [arXiv:2203.09595] [Google Scholar]

[91] Pérez, S., Casassus, S., Baruteau, C., et al. 2019, AJ, 158, 15 [Google Scholar]

[92] Perri, F., & Cameron, A. G. W. 1974, Icarus, 22, 416 [Google Scholar]

[93] Pinte, C., van der Plas, G., Ménard, F., et al. 2019, Nat. Astron., 3, 1109 [NASA ADS] [CrossRef] [Google Scholar]

[94] Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJ, 890, L9 [CrossRef] [Google Scholar]

[95] Pinte, C., Hammond, I., Price, D. J., et al. 2023, MNRAS, 526, L41 [NASA ADS] [CrossRef] [Google Scholar]

[96] Pirovano, L. M., Fedele, D., van Dishoeck, E. F., et al. 2022, A & A, 665, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[97] Piso, A.-M. A., Öberg, K. I., Birnstiel, T., & Murray-Clay, R. A. 2015, ApJ, 815, 109 [Google Scholar]

[98] Rafikov, R. R. 2006, ApJ, 648, 666 [Google Scholar]

[99] Rameau, J., Follette, K. B., Pueyo, L., et al. 2017, AJ, 153, 244 [Google Scholar]

[100] Robert, C. M. T., Crida, A., Lega, E., Méheut, H., & Morbidelli, A. 2018, A & A, 617, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[101] Rosotti, G. P. 2023, New Astron. Rev., 96, 101674 [NASA ADS] [CrossRef] [Google Scholar]

[102] Rosotti, G. P., Juhasz, A., Booth, R. A., & Clarke, C. J. 2016, MNRAS, 459, 2790 [Google Scholar]

[103] Schneider, A. D., & Bitsch, B. 2021, A & A, 654, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[104] Schoonenberg, D., & Ormel, C. W. 2017, A & A, 602, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[105] Shakura, N. I., & Sunyaev, R. A. 1973, A & A, 24, 337 [NASA ADS] [Google Scholar]

[106] Sierra, A., Pérez, L. M., Zhang, K., et al. 2021, ApJS, 257, 14 [NASA ADS] [CrossRef] [Google Scholar]

[107] Speedie, J., & Dong, R. 2022, ApJ, 940, L43 [NASA ADS] [CrossRef] [Google Scholar]

[108] Speedie, J., Booth, R. A., & Dong, R. 2022, ApJ, 930, 40 [NASA ADS] [CrossRef] [Google Scholar]

[109] Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65 [CrossRef] [Google Scholar]

[110] Szulágyi, J., Binkert, F., & Surville, C. 2022, ApJ, 924, 1 [CrossRef] [Google Scholar]

[111] Tanigawa, T., & Tanaka, H. 2016, ApJ, 823, 48 [NASA ADS] [CrossRef] [Google Scholar]

[112] Tanigawa, T., Ohtsuki, K., & Machida, M. N. 2012, ApJ, 747, 47 [NASA ADS] [CrossRef] [Google Scholar]

[113] Tanigawa, T., Maruta, A., & Machida, M. N. 2014, ApJ, 784, 109 [NASA ADS] [CrossRef] [Google Scholar]

[114] Teague, R., Bae, J., & Bergin, E. A. 2019, Nature, 574, 378 [NASA ADS] [CrossRef] [Google Scholar]

[115] Teague, R., Bae, J., Aikawa, Y., et al. 2021, ApJS, 257, 18 [NASA ADS] [CrossRef] [Google Scholar]

[116] Toci, C., Lodato, G., Fedele, D., Testi, L., & Pinte, C. 2020, ApJ, 888, L4 [Google Scholar]

[117] Tominaga, R. T., Tanaka, H., Kobayashi, H., & Inutsuka, S.-I. 2022, ApJ, 940, 152 [NASA ADS] [CrossRef] [Google Scholar]

[118] Uyama, T., Hashimoto, J., Kuzuhara, M., et al. 2017, AJ, 153, 106 [Google Scholar]

[119] Uyama, T., Norris, B., Jovanovic, N., et al. 2020, J. Astron. Telescopes Instrum. Syst., 6, 045004 [NASA ADS] [Google Scholar]

[120] van Clepper, E., Bergner, J. B., Bosman, A. D., Bergin, E., & Ciesla, F. J. 2022, ApJ, 927, 206 [NASA ADS] [CrossRef] [Google Scholar]

[121] van der Marel, N., Dong, R., di Francesco, J., Williams, J. P., & Tobin, J. 2019, ApJ, 872, 112 [Google Scholar]

[122] van der Marel, N., Bosman, A. D., Krijt, S., Mulders, G. D., & Bergner, J. B. 2021, A & A, 653, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[123] Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]

[124] Wang, Y., Ormel, C. W., Huang, P., & Kuiper, R. 2023, MNRAS, 523, 6186 [NASA ADS] [CrossRef] [Google Scholar]

[125] Xie, C., Haffert, S. Y., de Boer, J., et al. 2020, A & A, 644, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[126] Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]

[127] Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47 [Google Scholar]

[128] Zhang, K., Booth, A. S., Law, C. J., et al. 2021, ApJS, 257, 5 [NASA ADS] [CrossRef] [Google Scholar]