# Molecules with ALMA at Planet-forming Scales (MAPS) VIII: CO Gap in AS
209–Gas Depletion or Chemical Processing?

Felipe Alarcón  Department of Astronomy, University of Michigan, 323 West Hall, 1085 S University, Ave., Ann Arbor, MI 48109, USA  Arthur D. Bosman Department of Astronomy, University of Michigan, 323 West Hall, 1085 S University, Ave., Ann Arbor, MI 48109, USA  Edwin A. Bergin Department of Astronomy, University of Michigan, 323 West Hall, 1085 S University, Ave., Ann Arbor, MI 48109, USA  Ke Zhang NASA Hubble Fellow  Department of Astronomy, University of Wisconsin-Madison, 475 N Charter St, Madison, WI 53706  Department of Astronomy, University of Michigan, 323 West Hall, 1085 S University, Ave., Ann Arbor, MI 48109, USA  Richard Teague Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA  Jaehan Bae NASA Hubble Fellowship Program Sagan Fellow  Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA  Department of Astronomy, University of Florida, Gainesville, FL 32611, USA  Yuri Aikawa Department of Astronomy, Graduate School of Science, The University of Tokyo, Tokyo 113-0033, Japan  Sean M. Andrews Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA  Alice S. Booth Leiden Observatory, Leiden University, 2300 RA Leiden, the Netherlands  School of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK  Jenny K. Calahan Department of Astronomy, University of Michigan, 323 West Hall, 1085 S University, Ave., Ann Arbor, MI 48109, USA  Gianni Cataldi National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Department of Astronomy, Graduate School of Science, The University of Tokyo, Tokyo 113-0033, Japan  Ian Czekala NASA Hubble Fellowship Program Sagan Fellow  Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA  Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA  Center for Astrostatistics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA  Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA  Department of Astronomy, 501 Campbell Hall, University of California, Berkeley, CA 94720-3411, USA  Jane Huang NASA Hubble Fellowship Program Sagan Fellow  Department of Astronomy, University of Michigan, 323 West Hall, 1085 S University, Ave., Ann Arbor, MI 48109, USA  Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA  John D. Ilee School of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK  Charles J. Law Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA  Romane Le Gal Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA  IRAP, Université de Toulouse, CNRS, CNES, UT3, 31400 Toulouse, France  Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France  IRAM, 300 rue de la piscine, F-38406 Saint-Martin d’Hères, France  Yao Liu Purple Mountain Observatory & Key Laboratory for Radio Astronomy, Chinese Academy of Sciences, Nanjing 210023, China  Feng Long Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA  Ryan A. Loomis National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA  François Ménard Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France  Karin I. Öberg Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA  Kamber R. Schwarz NASA Hubble Fellowship Program Sagan Fellow  Lunar and Planetary Laboratory, University of Arizona, 1629 E University Blvd, Tucson, AZ 85721, USA  Merel L. R. van ’t Hoff Department of Astronomy, University of Michigan, 323 West Hall, 1085 S University, Ave., Ann Arbor, MI 48109, USA  Catherine Walsh School of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK  David J. Wilner Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA

(Received XXX; Revised XXX; Accepted XXX)

###### Abstract

Emission substructures in gas and dust are common in protoplanetary disks.
Such substructures can be linked to planet formation or planets themselves. We
explore the observed gas substructures in AS 209 using thermochemical modeling
with RAC2D and high-spatial resolution data from the Molecules with ALMA at
Planet-forming Scales(MAPS) program. The observations of C18O $J=2-1$ emission
exhibit a strong depression at 88 au overlapping with the positions of
multiple gaps in millimeter dust continuum emission. We find that the observed
CO column density is consistent with either gas surface-density perturbations
or chemical processing, while C2H column density traces changes in the C/O
ratio rather than the H2 gas surface density. However, the presence of a
massive planet ($>$ 0.2 MJup) would be required to account for this level of
gas depression, which conflicts with constraints set by the dust emission and
the pressure profile measured by gas kinematics. Based on our models, we infer
that a local decrease of CO abundance is required to explain the observed
structure in CO, dominating over a possible gap-carving planet present and its
effect on the H2 surface density. This paper is part of the MAPS special issue
of the Astrophysical Journal Supplement Series.

astrochemistry-protoplanetary disks-planet–disk interactions

††journal: AAS Journals††software: RADMC3D (Dullemond et al., 2012), RAC2D (Du
& Bergin, 2014), astropy (Astropy Collaboration et al., 2013, 2018)

##  1 Introduction

Surveys of protoplanetary disks have shown that substructures in the dust
continuum emission are ubiquitous (e.g. Andrews et al., 2018; Long et al.,
2018; Cieza et al., 2019). The presence of substructures such as gaps and
rings is seen as a signpost of incipient and ongoing planet formation that
alters the local disk’s physical conditions (Pinilla et al., 2012; Dong et
al., 2015). A natural question is whether the emission of common observable
species traces the local disks chemistry or physical changes caused by
planet–disk interactions.

The AS 209 system is an ideal laboratory to test the difference in gas
emission substructure caused by planet–disk interactions and local chemistry
as it presents multiple substructures in dust continuum and line emission
(Huang et al., 2016, 2018; Guzmán et al., 2018; Law et al., 2021a). These
substructures, in particular the millimeter dust continuum emission, have been
previously analyzed for a possible association with hidden planets (Dong et
al., 2018; Fedele et al., 2018; Favre et al., 2019). We analyze the gas
chemical and physical state using high-resolution ($\sim$0$\farcs$15) ALMA
observations of C18O and C2H from Molecules with ALMA at Planet-forming Scales
(Czekala et al., 2021; Oberg et al., 2021) 111www.alma-maps.info. The C18O
emission profile shows a wide emission depression centered at 88 au with a
width of 47 au, while C2H shows an emission ring centered at 70 au with a
width of 68 au (Law et al., 2021a). The gap in C18O emission translates into a
decrease in the CO column density, being depleted by at least 47% when
compared with a smooth surface–density profile (Zhang et al., 2021). Because
of its relatively simple chemistry and strong lines (Kamp et al., 2011;
Williams & Best, 2014; Molyarova et al., 2017), CO is frequently used as a H2
gas tracer assuming a uniform abundance across the disk. Therefore, such CO
depletion can be explained by the presence of a sub-Jovian planet carving a H2
gap (Favre et al., 2019).

In general, local CO abundance variations at substructures are considered to
be a minor problem compared to the question of the global CO abundance
relative to H2 and its possible deviations from the interstellar medium(ISM)
value (e.g. $\sim$10-4), necessary for the use of CO emission as a calibrated
tracer of the disk gas mass (Bergin & Williams, 2018). However, detailed
models suggest that standard scaling between CO and H2 surface density does
not hold locally inside the dust substructures due to CO chemical processing
(van der Marel et al., 2018; Alarcón et al., 2020).

Given the limitations of CO emission as a tool in probing disk structure in
planet-forming regions, additional chemical tracers observed in the MAPS ALMA
Large Program provide new constraints between planet–disk interactions and
local disk chemistry cases (Oberg et al., 2021). In this regard, the emission
from C2H stands out as its emission is often as bright as 13CO in some disk
systems (Kastner et al., 2014) and C2H emission is widely detected in gas-rich
systems often exhibiting emission rings (Kastner et al., 2015; Bergin et al.,
2016; Bergner et al., 2019). Numerical models suggest that the explanation for
such high C2H column densities is a localized high C/O ratio (Cleeves et al.,
2018; Miotello et al., 2019) and a photon-dominated chemistry necessary for
C2H production (Bergin et al., 2016).

In this work, we analyze two main scenarios explaining the CO column density
structure, CO chemical processing or H2 gas depletion caused by a giant
planet. Our aim is to disentangle the degeneracy between CO abundance and H2
surface density via a combined observational and numerical study of AS 209,
understanding the possible footprint of a newborn planet. Thus, we explore a
solution where chemical processing locally alters the CO abundance across the
dust gaps and another scenario in which planet–disk interactions carve a gap
in the H2 surface density assuming a constant CO abundance. We also explore
whether the conditions producing structure in the CO surface density also lead
to concurrent production of C2H.

This paper is organized as follows: We describe the setup of our
thermochemical models in Section 2. Section 3 presents the degeneracy between
CO abundance and gas surface density, and Section 4 shows the role of C2H as a
tracer of active chemistry. Then, we discuss the main conclusions of our work
in Section 5. A brief summary of our work is presented in Section 6.

##  2 Chemical and Physical Modeling

In our goal of reproducing the CO and C2H column-density profiles in AS 209,
we use a chemical and a physical approach to understand the origin of the
line-emission substructure. Our chemistry-dominated approach assumes that a
depleted CO abundance is required to reproduce structure in C18O emission in
the context of a smooth H2 gas surface-density profile. The physical approach
assumes a constant CO abundance and that the structure is the result of a
local decrease in the H2 surface density. We built the models so they match
the inferred CO column densities, and we compared the expected column density
of C2H and other tracers of the chemical and physical structure of the disk,
such as emission heights and kinematic deviation, to disentangle the
degeneracy between both solutions to the CO radial profile.

Additionally, we introduced a variable C/O ratio to both models keeping the CO
abundance constant as well. We also account for multiple dust substructures in
the disk for both large and small dust grains. For completeness, we provide a
detailed review of the AS 209 system as a planet formation laboratory in
Appendix A.

###  2.1 Chemical Network and Thermochemical Code

We used the 2D thermochemical code RAC2D, described in Du & Bergin (2014), to
model the chemical evolution in AS 209 for 1 Myr and to reproduce the inferred
CO and C2H column densities in the disk. The chemical network includes 524
species with 6425 chemical reactions. It uses the reaction rates from the
UMIST 2012 database (McElroy et al., 2013) for the gas-phase chemistry. RAC2D
also includes reactions that take place on the surface of dust grains using
the formalisms of Hasegawa et al. (1992). For the photodesorption of H2O and
OH by Ly$\alpha$ photons, the code uses yield values from Öberg et al. (2009a)
and the yield values from Öberg et al. (2009b) for CO2 and CO.

###  2.2 AS 209 Model

The AS 209 disk structure includes multiple gaps and rings beyond 20 au. We
adopted the disk parameters for AS 209 from Zhang et al. (2021) with small
changes due to the inclusion of an inner gap at 7 au with a 90% depletion of
dust and H2 gas, and a width of 3.8 au. We list the parameter of our models in
Table 1. The location and depletion of this inner gap is uncertain. However,
the analysis of Bosman et al. (2021b) using the MAPS data of AS 209 requires a
strong H2 gas depletion in the inner 10 au of at least one order of magnitude,
which is consistent with the presence of a gap or a gas cavity.

Table 1: Parameters for the AS 209 disk model.  Parameter | Value | Parameter | Value
---|---|---|---
Disk mass | 0.0045 $M_{\odot}$ | $r_{c}$ | 80 au
Large dust mass | 4 x 10-4 $M_{\odot}$ | $\gamma$ | 1.0
Small dust mass | 5 x 10-5 $M_{\odot}$ | $h_{c}$ | 6 au
$\zeta$ ionization rate | 1.36 x 10-18 s-1 | $\psi$ | 1.25

The grid in our models is defined on the $(r,z)$ plane. We have a radial
logarithmic spacing with 250 radial cells ranging from 1 to 250 au covering
our radial range of interest, which is between 20 to 120 au. The vertical
spacing in the grid is variable and depends on the local physical conditions
of the disk, with a spatial refinement in the layers where CO and C2H have
abundances with respect to H higher than 10-7 and 10-13, respectively. The
cells are never larger than 3 au.

####  2.2.1 Gas Structure

We modeled the disk density structure using two different approaches, calling
them Model A and Model B. For Model A we used a smooth self-similar
distribution (Lynden-Bell & Pringle, 1974) for the gas surface-density
profile, i.e,

$\Sigma_{\rm
smooth}(r)=\Sigma_{c}\left(\frac{r}{r_{c}}\right)^{-\gamma}\exp\left[-\left(\frac{r}{r_{c}}\right)^{2-\gamma}\right],$
(1)

where $r_{c}$ is the characteristic radius for the disk and $\Sigma_{c}$ is
the characteristic surface-density of the disk. Figure 1 shows the gas density
structure in Models A and B as well as the dust surface-density structure,
which is explained further in Section 2.2.2. The main difference between the
surface-density profiles in Models A and B is the presence of a H2 surface-
density drop scaled from the inferred CO column density at the location of the
C18O emission gap. We show the thermal structure and the differences between
Models in Appendix B.

We set the vertical distribution with a Gaussian centered at the midplane
(Armitage, 2020):

$\rho(r,z)=\rho_{0}(r)\exp\Big{(}-\frac{z^{2}}{2h^{2}}\Big{)},$ (2)

with $\rho_{0}=\frac{\Sigma(r)}{\sqrt{2\pi}h}$ and $h$ being the scale height,
which is different for the gas and the dust. The parameters and gas
distribution of Model A are the same as the one used in Zhang et al. (2021).
The scale height also has a radial dependence through a power-law with a
flaring index $\psi$:

$h(r)=h_{c}\Big{(}\frac{r}{r_{c}}\Big{)}^{\psi},$ (3)

with $\psi$ the flaring index and $h_{c}$ the characteristic scale height at
the characteristic radius, $r_{c}$.

Figure 1: Two sets of models with varying surface density for the small
grains. The small-grain surface density was changed by varying the depletion
in the 59 au gap with three different depletion factors, 1x (nominal), 10x and
100x. Top: In Model A, the surface density has no gas substructures, i.e., it
follows a smooth H2 surface-density profile, where CO is not a H2 mass tracer.
Bottom: Model B, showing the surface density using the CO structure as a H2
gas mass tracer, scaling with a constant CO abundance. Model B shows a wide
gas gap at 59 au. The shaded areas in yellow and blue are the dust gaps and
dust ring, respectively.

In Model A we assume a smooth surface-density profile for the H2 gas. To
account for the C18O $J=2-1$ emission depression at $\sim$80 au, we require a
localized change in the CO abundance. Here, we adopt the CO depletion profile
in Figure 2, i.e., the CO abundance is locally weighted according to that
profile. This depletion profile implies that the CO chemical processing in the
disk is taken into account effectively as an initial condition.

In Model B, we scaled the observed CO column-density profile to a H2 column-
density profile. We assumed a constant CO abundance of $2.8\times 10^{-4}$ and
then calculated the disk chemical evolution with that gas mass distribution.
In this scenario, the structure is caused principally by dynamical evolution
of the disk rather than chemically evolution of the CO reservoir. Model B
implicitly assumes that the H2 gas surface density has a wide gap ($w_{\rm
gap}$= 13 au) at 59 au. As a benchmark, we also compare whether or not Model B
matches the observed column densities for C2H and the effects the gas
depletion may have on the C2H chemical equilibrium.

The assumed initial abundances for Model B with respect to the total number of
hydrogen atoms are provided in Table 2. In order to match the CO column
density and the constraints given by the disk mass, our simulations started
with a global depletion of C and O that is reduced by one order of magnitude
compared to ISM values (Nieva & Przybilla, 2012). We did not include any
atomic carbon or water ice initially as they will be added to change the C/O
ratio locally for further exploration of the C2H chemistry.

Table 2: Initial Abundances of the simulations with respect to the number of H
Atoms.

Species | Abundance | Species | Abundance
---|---|---|---
H2 | 5(-1) | He | 9(-2)
HD | 2(-5) | S | 8(-9)
CO | 2.8(-4) | N | 7.5(-5)
Si+ | 8(-9) | Na+ | 2(-9)
Mg+ | 7(-9) | P | 3(-9)
Fe+ | 3(-9) | F | 2(-9)
Cl | 4(-9) |  |

**footnotemark:

In our simulations we have C/O=1. We only added atomic carbon or water ice to
change the C/O ratio locally, either to increase it or decrease it
respectively.

Note. — All the abundances are written in the standard form A(B) = A$\times
10^{B}$

####  2.2.2 Dust Structure

We have two dust populations in our models. The opacity of each dust
population is set by its composition and respective sizes. The large dust
population consists of dust grains with sizes $a$ ranging from 0.005 $\mu$m to
1 mm with a dust size distribution following the standard size distribution
from Mathis et al. (1977) and a total mass of 4.5x10-4 $M_{\odot}$. The large
grain population is composed of a mixture of 40% refractory organics (Henning
& Stognienko, 1996), 33% silicates (Draine, 2003), 20% water (Warren & Brandt,
2008) and 7% troilite (Henning & Stognienko, 1996) by mass. The X-ray
opacities follow the prescription from Bethell & Bergin (2011), which depends
on the size (cross section) of each dust grain. The surface-density
distribution for the large dust is different from that for the gas and the
small dust. In the large dust-surface-density distribution, the substructures,
i.e., gaps and rings, were added into the underlying dust surface-density
profile with the self-similar solution profile (Lynden-Bell & Pringle, 1974).
Both large dust rings have Gaussian shapes, as in Alarcón et al. (2020), and
their amplitudes were chosen so that the mass in the rings matches the
inferred mass from Dullemond et al. (2018); see Table 3 for the substructure
parameters in our models. The large dust surface density has a lower $\gamma$,
0.2 instead of 1. This $\gamma$ is different from our gas and small-dust
distribution as it provides a better fit for the dust continuum emission. A
lower $\gamma$ in the large dust grains has a steeper cutoff, implying that
the pebble or millimeter disk is more concentrated in the inner regions as
they suffer from more efficient radial drift (Birnstiel et al., 2010, 2012).
We used a different approach for the pebble disk than Zhang et al. (2021). We
used a prescription that analytically characterizes the gaps and rings using
previous models from Fedele et al. (2018) and Dullemond et al. (2018), instead
of the iterative approach of Zhang et al. (2021).

Table 3: Parameters for the substructures present in the dust surface density
of Millimeter-sized particles in our models.

Substructure | Location (au) | Width (au) | Depletion / Enhancement
---|---|---|---
Gap 1 | 7 | 3.8 | 0.9
Gap 2 | 62 | 8 | 0.991
Ring 1 | 74 | 3.5 | 7.3
Gap 3 | 100 | 16 | 0.975
Ring 2 | 120 | 4.11 | 26

*We use values from DSHARP (Andrews et al., 2018; Dullemond et al., 2018). The prescription used for the substructure is the same as in Alarcón et al. (2020), where the width and depletion/enhancement represent the width of a gaussian and the amplitude of the modulation of the surface density.

The small grains have sizes ranging from 0.005 $\mu$m to 1 $\mu$m with a total
mass of 5 x 10-5 $M_{\odot}$. We assume that the composition of the small
grains is different from that of the large grains with a composition of 50%
silicates and 50% refractory organics. Instead of following the large dust
distribution, the small grains follow the H2 gas surface density as they are
more dynamically coupled with the molecular hydrogen gas. Since the small-dust
surface density is less constrained in terms of its mass distribution (in
comparison to the large millimeter-sized dust), we used different localized
small-dust depletion factors where the CO column-density gap is located to
assess their effect on different CO and C2H, in particular their vertical
distribution and column density.

For both models, A and B, we ran three different simulations changing the
small-dust surface density. The only parameter that changes between each
simulation is the depth of the gap at $r=59$ au, while the rest of the
features of the gap, i.e. width and radius, remain the same. Changes in the
radiation field are expected as the small dust is the main source of UV
opacity in the upper layers, playing a significant role in the chemistry of
photochemical tracers such as C2H or HCN (Fogel et al., 2011). Moreover, small
dust also changes the temperature of the disk in intermediate layers as the
large dust is more settled and concentrated in the midplane.

Figure 2: CO depletion used in Model A to fit the CO intensity profile. This
depletion is used to compare the case with CO chemical processing (Model A)
and the case with CO used as a tracer of H2 surface density (Model B).

###  2.3 Variable C/O Ratio

In order to reproduce the high column densities of C2H (Bosman et al., 2021c;
Guzmán et al., 2021), we require deep UV penetration and an increase in the
C/O elemental ratio (Bergin et al., 2016; Cleeves et al., 2018; Miotello et
al., 2019). For each model, we tried two different approaches for the C/O
ratios to find if, under depleted conditions, high column densities at the
level observed by Guzmán et al. (2021) are recoverable with just a low C/O
ratio. In our first approach, we started with C/O=1 instead of 0.4, meaning
that the C/O ratio was already higher than the ISM value (Wilson & Rood,
1994); this is the expected value for the gas outside the CO ice line, which
is found between 15 and 20 au in our model. In the second approach, we added
extra carbon to the abundances listed in Table 2, which are constant in the
vertical direction. We therefore kept the CO abundance constant and added
extra neutral atomic carbon in the CO emission gap to increase the C/O ratio
to 2. Beyond 100 au, where the C2H emission drops significantly, we added
extra water ice to decrease the C/O ratio to 0.4, as water ice is the main
carrier of oxygen in disks (van Dishoeck et al., 2021). Figure 3 illustrates
the C/O ratio for the models and its radial dependence. The transitions in the
radial profile of the C/O ratio were empirically created to fit the C2H
inferred column density as best as possible. Irrespective of the location of
the C/O ratio transitions, a high C/O ratio is required to match the peak in
the C2H column-density profile. In addition to changing elemental abundances,
we also reduced the abundance of small grains; this locally increases the UV
field promoting the C2H production (Bosman et al., 2021a). We then compared
our results for the two sets of models, with a particular focus on C2H since
CO column density should not have a significant increase as the starting point
is C/O=1.

Figure 3: C/O ratio for our models with radial changes in the amount of C. The
C/O starts at 1, then it reaches a value of 2 in the gas gap where extra
atomic carbon is added. In this region most of the carbon grain destruction
would occur. It also coincides with the region where the C2H column density
reaches its highest values. Beyond 100 au we added extra water ice to decrease
the C/O ratio toward an ISM value C/O=0.4.

###  2.4 Emission Heights

We compared our models to the derived column-density structure for a given
species. As an additional constraint, we also compare the emission heights
provided in Law et al. (2021b) with the emission heights of CO isotopologue
lines from our models. Given the vertical structure in the CO abundance and
gas temperature, we used the region where the optical depth $\tau$ of each
line ranged between 0.67 and 2 as a proxy of each emitting layer. To produce
the different optical depth layers for each molecular transition we used the
temperature structure and CO distribution from RAC2D as input for the
radiative transfer code RADMC-3D (Dullemond et al., 2012).

We generated the emitting heights in the $J=2-1$ rotational lines of the two
most abundant CO isotopologues: CO and 13CO. As we did not include
isotopologue-selective photodissociation, we assumed a constant isotopologue
ratio for 13CO, setting CO 13CO=69 (Wilson, 1999). The possible effects of
isotopologue-selective photodissociation are further discussed in Section 5.2.
The C18O $J=2-1$ line is mostly optically thin in AS 209, so the emission
height is uncertain. Thus, we did not include it in the analysis of AS 209 in
this work.

##  3 CO Abundance and H2 Surface-Density Degeneracy

###  3.1 CO column densities

We show the CO column densities from the thermochemical simulation in Figure
4. Both sets of models, A and B, are able to recover the CO column density
from Zhang et al. (2021) with subtle differences, but still within the
uncertainties. Therefore, CO column densities alone cannot discriminate
between CO processing or the presence of a gas gap carved by a massive planet.
As we do not have a full understanding of the CO depletion in the disk, the
fact that both models match the CO column-density profiles illustrates the
degeneracy between the two solutions to explain the structure in CO emission
(Calahan et al., 2020). This is not unexpected, but we confirm that the gas
thermal physics does not change substantially between these solutions to favor
one model over the other as temperature is more dependent on the dust
structure in the midplane and intermediate layers of the disk. We show the
differences between the temperature structures of both in Appendix B.2; the
differences are usually less than 10%, so they may present some small radial
or vertical variation for the abundance of given species, but it does not
change the results of our models significantly. Therefore, we need to look at
other features in the emission to break the degeneracy. We provide a deeper
discussion of the 2D structure of the models and possible implication for line
emission in Appendix B.

Figure 4: CO column-density profiles for both sets of simulations. Top: smooth
gas surface density applying the CO depletion profile in Figure 2. Bottom:
Simulation set with CO as a gas mass tracer. The vertical shaded areas are the
location of dust substructures, while the radial shaded regions represent the
CO column densities derived in Zhang et al. (2021) using 13CO and C18O
emission. In both cases, Models A and B, we observe that whether including CO
chemical processing or the gas depletion match the recovered CO column
densities. The shaded regions are the locations of the dust gaps and dust
rings, in yellow and blue, respectively.

By changing the amount of small dust in the gap at 59 au, we do not observe
noticeable differences in the CO column densities. However, changes in the CO
emission are still expected as the small dust is a key player in setting the
disk temperature, particularly in the warm molecular layer on top of the large
dust.

###  3.2 CO Emission Heights: Model versus Data

The observed emission heights present an additional constraint, in particular
for the CO abundance structure. The comparison between the emitting heights in
our models and the ones inferred from the data are shown in Figure 5 (Law et
al., 2021b). Here, it is shown that for the optically thick CO, the emitting
heights for both models, A and B, match, within the uncertainties, the data.

For 13CO, the observed emitting height is constrained in a narrower radial
range and closer to the midplane. Nevertheless, Models A and B still agree
with the constraints given by the observations. Because both models are within
the uncertainties and match emission heights, the degeneracy between these
models remains.

Figure 5: Emitting heights measured in Law et al. (2021b) compared to the
layers in the disk where the optical depth ranges from 0.67 to 2 for each
line. This layer serves as a probe of the emitting height as most of the
contribution of the line should come from these heights. The comparison seems
to coincide well with the surfaces of CO and 13CO. Models A and B are within
the range given by the observations.

###  3.3 Pressure Profile and Planets, Breaking the Degeneracy

CO column densities and emission heights do not distinguish between chemical
processing and planet–disk interactions. Therefore, we explore the thermal H2
gas thermal pressure in our models along the CO-emitting layers and compare
them to that obtained by Teague et al. (2018) from gas kinematics in AS 209.
We expect the pressure profiles to be different around the CO emission gap as
the H2 gas surface density is different in Models A and B.

####  3.3.1 Pressure Profiles and Kinematic Deviations

We made two radial cuts enclosing the emission heights of 12CO and measured
the pressure profile along those layers. One cut is assumed to follow the
emitting height:

$h=23\Big{(}\frac{r}{100}\Big{)}^{1.05}\mathrm{au},$ (4)

consistent with the layer explored by Teague et al. (2018). The second cut was

$h=26.5\left(\frac{r}{121}\right)^{1.292}\exp\left(-\left(\frac{r}{216}\right)^{4.9}\right)\
\mathrm{au},$ (5)

which is the parametric CO $J=2-1$ emission height found by (Law et al.,
2021b). A comparison between each layer is illustrated in Figure 16.

We show a comparison between the pressure profiles in our models along the
layers given by Eq. 4 and 5, and the observed pressure profile from Teague et
al. (2018) in Figure 6. As we do not include hydrodynamics in our models, we
look at the trend of each model rather than matching the empirically derived
pressure profile, which is more directly constrained by its slope and
structure than its absolute magnitude. The thermal variation between models,
in particular, along the CO emission layer is $\sim$10% (see Appendix B.2).
Therefore, the variations in the pressure profile are mostly induced by radial
density variations, as thermal variations at a 10% level do not cause
significant changes to the pressure profiles. The pressure profile from Model
A is smooth. Conversely, Model B exhibits a considerable pressure dip in the
region where CO is being depleted in Model A, tracing a one-order-of-magnitude
density drop in gas surface density. When we compare our pressure profiles
with the one inferred from gas kinematics (Teague et al., 2018), the pressure
drops seen in Model B are not reproduced. Instead the pressure analysis
suggests a smoother H2 surface-density profile (which could be interpreted as
smooth CO changes) as assumed in Model A. We note that the actual pressure
values do not exactly match Teague et al. (2018). This is particularly the
case when we use the emitting surfaces from Law et al. (2021b). Regardless,
radial variation in the pressure profile probes changes in the H2 gas surface
density.

The MAPS program has obtained higher-resolution and higher-S/N data than used
by Teague et al. (2018) and we can revisit this question with independent
data. The deviations from the Keplerian field as derived from MAPS CO J=2–1
data are shown in Appendix C (Fig. 17) and are of order 2% in a relative sense
($\delta v/v$). In Fig. 7 we illustrate the expected kinematic deviations in
our models on this same scale using the relations for the gas structure from
Rosotti et al. (2020) applied to our models. The Keplerian deviations in Model
B range over 30% (peak to peak) within the framework of the gap and are
clearly inconsistent with the data that show only small $\sim 2$% deviations.
Model A on the other hand does show some structure but is limited to
deviations of order 8 % at the gap edges. Fitting the kinematic deviations in
the AS 209 disk goes beyond the scope of this paper, but a large drop in H2
surface density of an order of magnitude would induce almost a factor of 10
larger than observed in the data.

Kinematical analysis within gas surface-density gap scenarios have also been
explored by Rab et al. (2020); in this case for the HD 163296 disk. This is a
different disk; however, their findings about kinematic deviations in gas-
depleted gaps are consistent with our results for AS 209. Keplerian deviations
in gas-depleted gaps are significant ,$\sim$ 10%, compared to dust gaps
without gas depletion. Thus, analyzing the kinematic data in two independent
data sets supports smoother transitions at the location of the CO column-
density gap and the C18O emission gap. A significant H2 gas depletion causes
more abrupt variations in both profiles that are not observed in the data.

Associating the gap in CO emission with gas depletion also represents a
disconnect with the planetary mass required to carve the gaps in the dust and
gas mass distribution. To reach a $\sim$ 50% CO depletion in AS 209 through a
decrease in the H2 surface density (Zhang et al., 2021), a 0.2
$M_{\mathrm{Jup}}$ planet would be required. Such a planet carving the dust
gaps would also need to perturb the gas in a much wider region encompassing
the two dust gaps without a common center. Thus, the disconnect between the
dust and the gas, in terms of the mechanisms for gap opening, invalidates a
general H2 depletion in the disk, supporting the idea that the CO distribution
in AS 209 is the result of a chemical effect on CO.

We are not able to rule out gaps in the H2 gas entirely, as there are small
fluctuations in the observed pressure profile that hint at gas pressure
differences corresponding to the gaps seen in the dust distribution. However,
any perturbation to the gas surface density must be rather shallow and less
broad that what is observed for CO emission.

Figure 6: Thermal pressure profile cuts in our two models compared to the
profile from Teague et al. (2018) inferred from AS 209 disk kinematics using
CO line emission. The other cut was taken following the parametric emission
surface of CO inferred by Law et al. (2021b). When the gas is depleted by
roughly an order of magnitude, as in Model B, the pressure profile shows
significant variations with the radius, these are inconsistent with those from
Teague et al. (2018). The lack of strong fluctuations in the observed pressure
profile disfavors the presence of a deep gas gap, at least at $\sim$ 59 au,
breaking the degeneracy between chemical processing and gas depletion. Thus,
chemical processing is the dominant mechanism for CO column variations in AS
209. Figure 7: Expected velocity deviations in our models inferred from the
pressure profiles. We inferred the velocity deviations using the relationships
described by Rosotti et al. (2020). The shaded area illustrates the range of
kinematic deviations from the data (see Appendix C). In Model A the deviations
are significantly smaller than in Model B. Model A does not present strong
kinematic deviations from its pressure profile, while Model B with a gas gap
shows deviations of the order of 10% in both surfaces.

####  3.3.2 Testing Kinematic Deviations in Hydrostatic Equilibrium Models

We also test whether or not our result, i.e., differences in the kinematic
deviation in models with and without a H2 density structure, stands in models
when the hydrostatic equilibrium is considered. We look at the third iteration
of a thermochemical model with hydrostatic equilibrium to see its effects on
the kinematic structure. Doing a detailed analysis of the effect of
hydrostatic equilibrium is not the goal of this work; a more complex and
detailed analysis is required. Nevertheless, our test allows us to understand
if accounting for the hydrostatic equilibrium changes our findings. We test
our results by comparing the 12CO $J=2-1$ radial emission profiles and
estimated temperature of the gas. This represents a change from our previous
practice of comparison to the CO column densities. For this test we match the
12CO $J=2-1$ emission profile. We start our vertical hydrostatic equilibrium
runs with the models matching the CO $J=2-1$ emission as they get us closer to
a smoother solution without impactful changes in fewer iterations, while
matching the CO emission.

We show the result of our test in Figure 8. In order to reproduce the
spatially resolved CO $J=2-1$ emission profile we had to increase the
thickness of the disk with respect to Zhang et al. (2021) to $h/r$=0.8 and the
flaring index to $\psi$=1.35. Using those values we are able to match the CO
emission in our region of interest. We refer to these models with hydrostatic
equilibrium iterations as Model Ah and Model Bh. The temperature along the CO-
emitting height in these models is within 10 K from the values inferred from
the data (Law et al., 2021a) in Models Ah and Bh. By slightly changing the
thickness and flaring of the disk in the model to reproduce 12CO emission,
these models find the warm molecular layer, the vertical zone where CO is not
frozen onto grains or photodissociated, to be vertically higher in the disk.
This also changes the overall CO abundance structure (see Appendix C) which
may not fully reproduce the overall CO column density. Nonetheless, we are
able match the 12CO emission profiles even in models with hydrostatic
equilibrium turned on.

To achieve consistency and uniformity between models we use the emitting
layers for the hydrostatic test as the ones that fit the observed CO $J=2-1$
brightness temperature in the disk (Law et al., 2021a) as shown in Fig. 8
(middle panel). Even with 10 K deviations in temperature, the major player in
the pressure profile is the density along the emitting layer, which dominates
any changes in the pressure profile. Thus in models with a smooth density
profile (Models A and Ah) we see small kinematical deviations, which is not
the case for Models B and Bh where the H2 surface density is assumed to
contain a large gap. We observe that the kinematic deviations are smoother due
to the hydrostatic balance (e.g. compare kinematic residuals in Fig. 7 and
Fig. 8); however, we still observe that the deviations in the model with H2
gas depletion, while having values $\sim$ 4 %, cover a range of the order of
7%-8% peak to peak, which is stronger than that seen in the CO kinematical
analysis of $\sim$ 2% at most ($>$3.5 $\sigma$). Thus, we can conclude that
our results stand in our test with hydrostatic equilibrium models, i.e., CO
depletion is the main cause of the C18O $J=2-1$ emission gap in AS 209.

Figure 8: Comparison of data and expected values of thermochemical models with
three iterations of hydrostatic equilibrium. Top: CO $J=2-1$ emission profile.
Middle: thermal cuts for both models matching the brightness temperature of
the CO $J=2-1$ line (Law et al., 2021a). Bottom: expected kinematic deviations
in Models Ah and Bh with the range from the observations represented by the
shaded area. The models show that the CO $J=2-1$ is, to a small uncertainty,
well traced in both cases. The kinematic deviations in Models Ah and Bh show
that even when a few iterations of hydrostatic equilibrium are considered, the
results stand, i.e., the kinematic deviations in Model Bh with a H2 gas
depletion are stronger and cover a larger range than the ones observed in the
data.

##  4 C2H as a Tracer of Active Chemistry

###  4.1 C2H column densities

If the abundance of CO is being reduced in planetary gaps due to changes in
the physical conditions of the disk (Alarcón et al., 2020), we expect that
other species might also be affected as CO is the main carrier of carbon and
oxygen in the disk gas. So, we look at possible fingerprints of CO processing
using C2H, as it has bright emission lines that are commonly observed in
protoplanetary disks (Kastner et al., 2014; Bergner et al., 2020; Miotello et
al., 2019). In Figure 9 we compare the C2H column densities in both Models A
and B with the retrieved values in Guzmán et al. (2021) for AS 209. In Model
A, we show solutions with C/O = 1 and C/O = 2 between 20 and 110 au to
demonstrate that an elevated C/O ratio is a requirement to increase the C2H
production matching the observations, with a subsequent reset at $r\sim$ 110
au to C/O = 0.4. This result is reproduced in Model B as well. Thus, the C/O
ratio is more important in determining the C2H column density than the actual
surface-density profile, agreeing with the results in Bosman et al. (2021c)
for the case of a smooth disk. Model A has a larger C2H column density
compared to Model B for the elevated C/O = 2 ratio profile, while that trend
is not clear when C/O = 1. C2H production depends on the local density and UV
field (Bosman et al., 2021a). Thus, the C2H gas is localized in a narrow layer
in the disk (see Figure 15 tracing the local condition in that region.
Therefore, C2H is not a reliable radial surface-density tracer as it does not
show a consistent correlation between Models A and B when the C/O ratio or the
small-grain abundance is changed. Our models also assumed that the C/O ratio
does not have vertical variations at a given radius. It is possible that the
same C2H can be reproduced by vertically localized increments in the C/O ratio
without changes in the H2 surface density. Such increases are consistent with
results from numerical models (Krijt et al., 2020). An elevated C/O ratio
dependence for C2H insensitive to surface-density changes is also consistent
with the findings of Miotello et al. (2019) that C2H flux does not strongly
correlate with the disk or dust mass. We note that even with a high C/O ratio
and large depletion in the abundance of small grains, we still underestimate
the C2H column density in Model B.

The C2H column-density peak seems to require a radial trend in the C/O ratio
(Figure 3), with a higher C/O ratio where there is a higher C2H column
density. We note that these increases in the C/O ratio occur in the same disk
regions where the CO abundance is being lowered. Even though both effects may
be related, we are not able to state a clear link between them as there are
other possible chemical processes producing the radial variation in the C/O
ratio.

Figure 9: C2H column-density profiles compared with the inferred values from
observations. Top: simulation set for Model A. Bottom: simulation Set for
Model B. The shaded areas in yellow and blue are the dust gaps and dust ring,
respectively. Variable C/O models are shown in continuous lines, while the
dashed lines are the models with a constant C/O ratio equal to 1. Our
simulations show that depleting the number of small dust grains increases the
C2H column density by a factor of a few. However the main player in the C2H
column density is the C/O ratio.

##  5 Discussion

###  5.1 planet–disk Interaction vs Chemical Processing

The CO column-density gap morphology, i.e., width, depth, and location, allows
us to put an upper bound on the mass of the planet by using the relationships
described by Kanagawa et al. (2016, 2017). Any planet more massive than this
upper bound would carve a gas gap deeper than the CO drop in column density.
Using the values from Zhang et al. (2021) for the gap in CO column density;
the gap width, $w_{\rm gap}=13$ au, the gap location $R_{\rm gap}=59$ au and
the depletion $\delta_{\rm gap}=0.47$ with an $\alpha$-viscosity = 10-3, the
mass of the gap-opening planet would be that of a sub-Jovian planet with a
mass of 0.2 $M_{\rm Jup}$, thus this is still consistent with the values
reported by Fedele et al. (2018), Zhang et al. (2018) and Favre et al. (2019).
Therefore, CO being chemically processed does not disagree with the presence
of a sub-Jovian planet carving the dust gaps in AS 209, similar to the
scenario proposed by Dong et al. (2017).

###  5.2 Carbon and CO Chemistry in Planet-forming Disks

There are some significant CO chemical processes that support the CO depletion
scenario discussed here. One such process that is not explicitly included in
our modeling is CO isotopologue-selective photodissociation. Miotello et al.
(2014) showed that, particularly for the case of C18O, isotopologue-selective
photodissociation would decrease the inferred CO column density and the line
emission as well. However, the CO column-density profiles inferred by Zhang et
al. (2021) using 13CO and C18O line emission independently are in agreement,
showing that isotopologue-selective photodissociation is not significant in AS
209. Even if its effect were to be significant, a strong isotopologue-
selective photodissociation effect supports the scenario in which the emission
gap in C18O is explained by changes in the abundance of C18O rather than a
strong gas depletion caused by a massive planet.

Based on our mass-independent analysis of the C18O emission structure using
the pressure profile, we find that the abundance of CO is lowered locally in
the AS 209 disk. The localized abundance depletion of CO in AS 209 reduces the
overall oxygen in the system, likely enabling C2H production (Bosman et al.,
2021a). Water ice is likely frozen and locked in the midplane (Hogerheijde et
al., 2011; Du et al., 2017) and it is CO that is the main carrier of oxygen in
the outer disk beyond the CO2 ice line at $\sim$ 5 au. However, our results
also require extra available carbon. This appears to be a common result as
Bosman et al. (2021c) show that at least three of the MAPS disks require
C/O$>$ 1 and a depletion of the small-grains abundance. The origin of this
excess carbon is a matter of debate. Carbon grain destruction might provide
the source term (Anderson et al., 2017; Gail & Trieloff, 2017; Klarmann et
al., 2018); this would be consistent with a reduction in small-grains optical
depth and extra gaseous carbon. Alternately, it might come from a two-step
process where CO is destroyed via ubiquitous He+ ions followed by water-ice
formation removing oxygen from the gas phase (Bergin et al., 2014; Reboussin
et al., 2015; Yu et al., 2017; Schwarz et al., 2019). Our model is not a self-
consistent exploration of these processes and we cannot distinguish between
these two scenarios. Whether the chemical processing of CO leads to abundance
depletion or it is anticorrelated with the production of C2H is uncertain and
requires detailed models that couple dust growth to chemistry, such as those
shown by Krijt et al. (2020).

From an observational perspective, a previous survey of a small sample of
circumstellar disks did not show any anticorrelation between the integrated
13CO flux and C2H luminosity (Miotello et al., 2019). Bergner et al. (2020)
also showed that there is no global anticorrelation between C18O and C2H
column densities, which is consistent with the behavior of the other MAPS
sources (Law et al., 2021a; Guzmán et al., 2021). We attribute this lack of
correlation to the radial structure of the CO abundance, hiding this chemical
processing in unresolved observations. C2H becomes abundant in narrower
regions in the disk, while CO is present throughout the disk.

Figure 10: Normalized line emission of select species. Top: radial profiles of
hydrocarbons and CO. We observe a multiple peak structure within the gas gap.
Bottom: radial profiles of nitriles in AS 209. They peak close to 25 au, one
of the inner continuum-emission rings. The black shaded area represents the CO
ice line in the midplane.

We also show the normalized radial intensity profiles of other carbon tracers
in the upper panel of Figure 10 (Law et al., 2021a). While C18O exhibits a
broad emission gap, other carbon tracers such as c-C3H2 and H2CO have emission
rings that coincide with that of C2H. Thus the abundance reduction of CO
powers the production of these organic compounds. Moreover, C2H and c-C3H2
show a broad and coherent ring and have similar formation paths (Henning et
al., 2010), while H2CO has multiple structures. H2CO peaks just outside the
inferred midplane CO ice line, while H2CO also has additional structure just
beyond the C2H ring. The outermost gap for H2CO at $\sim$ 150 au might be
associated with the second CO ice line located near the edge of the pebble
disk (Cleeves, 2016).

Another important chemical tracer of the local gas-to-dust ratio in the disk
is DCO+. Smirnov-Pinchukov et al. (2020) show that a dust-poor but gas-rich
gap produces more DCO+ and HCO+, increasing their column density by several
orders of magnitude. Observing these tracers at high spatial resolution, i.e.,
resolving the dust gaps, would provide extra support to disentangle gas- and
CO-depleted scenarios required to explain the CO emission in dust gaps. From
our current analysis, we expect strong HCO+ and DCO+ emission in the CO
column-density gap, which is observed in Favre et al. (2019) and Aikawa et al.
(2021).

The structure in the inner disk ($<$0$\farcs$4) is intriguing as this would be
chemistry associated within the inner tens of astronomical units of this
system. The bottom panel in Figure 10 shows the distribution of some nitriles
in AS 209. It is clear the emission from nitriles and CO falls off at $\sim$
50 au while other species such as C2H behave differently, starting to rise up
after the inner 50 au. This distinct behavior in AS 209 points at a different
chemistry between nitriles and hydrocarbons. The origin of this chemistry is
uncertain. This bright nitrile emission hints at an active chemistry
associated with the inner disk that favors the production of these species in
the disk’s primary planet-forming region, with implications for the C/N ratio
of compounds in the gas and within solids.

##  6 Summary

We investigated the gas structure of AS 209 by analyzing high-resolution data
from the MAPS ALMA Large Program (Oberg et al., 2021). These data show a broad
depression in C18O emission and a concurrent rise in C2H emission around
$\sim$ 60 au. We compared two possible scenarios to explain the observed C18O
emission gap: a local depletion in the gaseous H2 surface density or a local
CO abundance depletion. Both scenarios provided good fits of the observed CO
column density and reproduced the 12CO and 13CO $J=2-1$ emission heights.
However, the pressure profile as estimated by local velocity deviations from
the Keplerian flow by Teague et al. (2018) and our own reanalysis of the
kinematical deviations from the Keplerian flow in independent MAPS CO data,
shows that the presence of a broad (13 au) order-of-magnitude deep gas gap is
not detected. Thus we conclude that the CO abundance in AS 209 is locally
reduced within a nearly smooth H2 density profile and that chemical processing
is active in this system. Based on the lack of significant structure in the
overall gas density distribution, the presence of a Jovian-mass planet($>$0.2
$M_{Jup}$) in AS 209 is disfavored for $\alpha\leq 10^{-3}$. We cannot rule
out the presence of a less-massive planet, which would be consistent with
previous modeling of the CO emission substructure on the source (Favre et al.,
2019), although at a lower spatial resolution.

We also show that the chemistry that forms C2H is insensitive to H2 gas
depletion, but in the AS 209 disk it appears to be correlated with the local
reduction of the CO abundance. We speculate that this is due to the impact of
the oxygen carried by CO which, if present, would hinder C2H production by
maintaining a C/O ratio lower than 1. Moreover, the processing associated with
CO depletion can potentially produce C/O ratios that exceed unity. Beyond this
association, there appears to be a rich and active chemistry associated with
the early stages of planet formation in this system that remains to be
understood.

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01055.L.
ALMA is a partnership of ESO (representing its member states), NSF (USA) and
NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI
(Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA
Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio
Astronomy Observatory is a facility of the National Science Foundation
operated under cooperative agreement by Associated Universities, Inc. We thank
the anonymous referee for the feedback and constructive comment on this paper.
F.A, A.D.B. and E.A.B. acknowledge support from NSF AAG Grant #1907653. K.Z.
acknowledges the support of the Office of the Vice Chancellor for Research and
Graduate Education at the University of Wisconsin–Madison with funding from
the Wisconsin Alumni Research Foundation, and support of NASA through Hubble
Fellowship grant HST-HF2-51401.001. awarded by the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., for NASA, under contract NAS5-26555. J.B. acknowledges
support by NASA through the NASA Hubble Fellowship grant #HST-HF2-51427.001-A
awarded by the Space Telescope Science Institute, which is operated by the
Association of Universities for Research in Astronomy, Incorporated, under
NASA contract NAS5-26555. R.T. and F.L acknowledge support from the
Smithsonian Institution as a Submillimeter Array (SMA) Fellow. K.R.S.
acknowledges the support of NASA through Hubble Fellowship Program grant HST-
HF2-51419.001, awarded by the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
for NASA, under contract NAS5-26555. J.K.C. acknowledges support from the
National Science Foundation Graduate Research Fellowship under grant No. DGE
1256260 and the National Aeronautics and Space Administration FINESST grant,
under grant no. 80NSSC19K1534. Y.A. acknowledges support by NAOJ ALMA
Scientific Research Grant CODE 2019-13B and Grant-in-Aid for Scientific
Research (S) 18H05222, AND gRANT-IN-aID FOR transformative Reseacrh Areas (A)
20H05844 and 20H05847. J.D.I. acknowledges support from the Science and
Technology Facilities Council of the United Kingdom (STFC) under ST/T000287/1.
K.I.Ö. acknowledges support from the Simons Foundation (SCOL #321183) and an
NSF AAG Grant (#1907653). I.C. was supported by NASA through the NASA Hubble
Fellowship grant HST-HF2-51405.001-A awarded by the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., for NASA, under contract NAS5-26555. C.J.L. acknowledges
funding from the National Science Foundation Graduate Research Fellowship
under Grant DGE1745303. C.W. acknowledges financial support from the
University of Leeds, STFC and UKRI (grant Nos. ST/R000549/1, ST/T000287/1,
MR/T040726/1). M.L.R.H. acknowledges support from the Michigan Society of
Fellows. R.L.G. acknowledges support from a CNES fellowship grant. F.M.
acknowledges support from ANR of France under contract ANR-16-CE31-0013
(Planet-Forming-Disks) and ANR-15-IDEX-02 (through CDP “Origins of Life”).
G.C. is supported by NAOJ ALMA Scientific Research grant code 2019-13B. S.M.A.
and J.H. acknowledge funding support from the National Aeronautics and Space
Administration under grant No. 17-XRP17 2-0012 issued through the Exoplanets
Research Program. J.H. acknowledges support for this work provided by NASA
through the NASA Hubble Fellowship grant #HST-HF2-51460.001-A awarded by the
Space Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., for NASA, under contract
NAS5-26555. Y.L acknowledges the financial support by the Natural Science
Foundation of China (Grant No. 11973090).

##  Appendix A AS 209 as a planet formation laboratory

In this work we focused on the C18O emission gap at $\sim$80 au and the C2H
emission ring around the same location, while the dust continuum presents one
ring at this location with two gaps at each side of it (see Fig 11). It is
noteworthy that understanding the origin of substructures in line emission is
more complex than in dust continuum. The dust substructures are linked to the
dynamical evolution of the disk, while substructures in line emission can also
be linked to the disk’s chemical evolution.

Several studies have proposed the presence of one or more planets at different
locations within AS 209 as the origin of the dust rings and gaps. Zhang et al.
(2018) proposed that a single sub-Jovian planet with an $\alpha\sim$ 10-4 (0.2
$M_{\rm Jup}$) is able to carve the multiple gaps in AS 209 (Bae et al., 2017;
Dong et al., 2017). With this planet located at 100 au, models suggest that
the resulting interactions can explain the 100 au gap and all the interior
dust substructure. Fedele et al. (2018) explored the possibility of a second
planet in the inner dust gap at 62 au. They found that the presence of a 0.05
$M_{\rm Jup}$ planet also fits the dust substructure, although such a planet
is not necessary to explain the observations. Dong et al. (2018) suggested
that a 0.09 $M_{\rm Jup}$ planet embedded in the 80 au continuum ring could
explain the location and depth of the dust gaps at each side of the ring.

Figure 11: C2H and C18O radial profile for the velocity-integrated emission
lines in ALMA Band 6, while the gray line displays the continuum-emission
profile from DSHARP (Andrews et al., 2018; Guzmán et al., 2018) on the left
and from MAPS (Oberg et al., 2021) on the right. We observe that C18O shows a
wide gap before showing a second peak after the outer ring at 120 au. C2H
shows an inverse correlation to C18O, peaking at the location of the C18O gap,
making AS 209 a unique source for the study of CO chemical processing among
the sources in the MAPS survey.

Favre et al. (2019) were the first to constrain on a possible planet scenario
in AS 209 using gas emission from CO at a spatial resolution of $\sim$
0$\farcs$25\. They claimed the need of gas depletion to match the C18O gap in
the emission profile, which can be explained by a 0.2 MJup planet at the outer
dust gap at $\sim$100 au. Independently, Teague et al. (2018) find velocity
perturbations around three locations in the disk associated with radial
changes in the pressure profile from CO observations. One of these
perturbations was associated with a wide gas gap at $\sim$ 50 au with 80% gas
depletion.

###  A.1 Comparison with DSHARP

Figure 11 shows the radial-velocity-integrated intensity profile in C18O
$J=2-1$ and C2H $N=3-2$ from the MAPS survey in ALMA Band 6, along with the
dust continuum emission from MAPS and DSHARP at a higher spatial resolution
(Andrews et al., 2018). The C18O emission profile shows a wide emission
depression centered at 88 au with a width of 47 au, while C2H shows an
emission ring centered at 70 au with a width of 68 au (Law et al., 2021a). The
gap in C18O emission translates into a decrease in the CO column density,
being depleted by at least 30% when compared with a smooth surface-density
profile (Zhang et al., 2021). Even at the MAPS spatial resolution
(0$\farcs$15), we still observe that the gap in the C18O $J=2-1$ emission is
broader than the dust continuum-emission profile, showing that the structure
is resolved.

We also ran a model using the gas surface-density profile matching the dust
continuum emission from DSHARP using hydrodynamical simulations (Zhang et al.,
2018). The gas surface density was scaled with the disk mass used in our
simulations and the same abundances that were used in Model B (see Table 2).
We show the comparison in Figure 12 illustrating that CO depletion is needed
to obtain the inferred CO column densities. Despite the hydrodynamical model
in DSHARP fitting the dust continuum emission at a higher resolution, the dips
and peaks in the CO column density do not match the ones in the profile from
DSHARP, which supports the solution of dominant CO chemical processing rather
than H2 gas depletion.

Figure 12: Comparison between the inferred CO column densities from Zhang et
al. (2021) and a model using the gas surface-density from the hydrodynamic
simulations in DSHARP (Zhang et al., 2018). The gas surface density profile
was scaled to the disk mass used in our simulations and with the same
abundances listed in Table 2 without any extra CO depletion. The profile using
the the DSHARP gas surface density overestimates the CO column density by an
order of magnitude, and the location of the dips and peaks do not match the
ones observed in CO, supporting the CO-processing solution.

##  Appendix B Models Setup and Structure

###  B.1 Substructures

The gaps and rings were included following the same prescription as Alarcón et
al. (2020), i.e., gaps and rings are a Gaussian modulation of the local
surface density for a smooth parametric disk following the self-similar
solution on an $\alpha$-viscosity disk (Lynden-Bell & Pringle, 1974).

Gaps have the following functional form:

$\Sigma(r)=\Sigma_{\rm basic}(1-(1-\delta_{\rm gap})\exp(-(r-r_{\rm
gap})^{2}/2w_{\rm gap}^{2})),$ (B1)

where $\delta_{\rm gap}$ is the depletion factor, $r_{\rm gap}$ the gap’s
location and $w_{\rm gap}$ the width of the gap. Rings are enhancements of
dust grains so we parameterize them with this equation:

$\Sigma(r)=\Sigma_{\rm dust}(1+\delta_{\rm ring}\exp(-(r-r_{\rm
ring})^{2}/2w_{\rm ring}^{2})),$ (B2)

with $\delta_{\rm ring}$ the enhancement factor, $r_{\rm ring}$ the ring’s
location, and $w_{\rm ring}$ the width of the ring.

###  B.2 2D Thermal Structure Comparison

Figure 13: Thermal structure for Models A and B. Top row: dust temperature.
Bottom row: gas temperature. The right column shows that the normalized
residuals between models with $\Delta T=T_{A}-T_{B}$ is the difference between
the temperature of Model A and Model B. The residuals are on the order of 10%,
which corresponds to 3–4 K degrees at the emission layer.

We show the 2D temperature field for each model in Figure 13 and their
respective difference. The differences in both cases are coming from the
differences in the CO abundance between each model and small variations in the
gas structure which are more significant in the gap. Nevertheless, the
differences are usually less than 10% between each model, so they may present
some small radial or vertical variation for the abundance of given species,
but it does not change the results of our models significantly. However, such
differences between gas depletion and CO processing could potentially be
traced in a deep analysis of the line-emission profiles of molecular tracers.

###  B.3 2D Abundance Structure of CO and C2H

We present the abundance structure in Figure 14 showing that in Model A, when
CO has been depleted, the absolute CO abundance is lower than in Model B,
which is compensated by a higher gas density. When we deplete CO, it has a
more uniform abundance with height when compared to the second scenario.
Therefore, even if in both scenarios the CO column-density radial profile is
similar, the changes in the vertical distribution of CO could lead to changes
in the emission radial profiles of CO isotopologues. By having different
vertical structure, a detailed comparison of CO isotopologue emission profiles
could provide differences between the emission heights of both models.
Nevertheless, such approach goes beyond the scope of this paper.

If we compare the abundance structure of C2H for each case (see Fig. 15), we
observe that even though the vertically integrated column densities show
changes of less than an order of magnitude, depleting the small grains causes
the C2H to be produced closer to the midplane, consistent in both Models A and
B. Nevertheless, even though Models A and B produce C2H at lower heights, the
layers at which C2H is being produced differ in the models. In Model B, there
is more production of C2H, i.e., higher C2H abundances at lower heights. We
link the higher C2H abundances to a deeper UV penetration and a longer C2H
prevalence due to a lower gas density, i.e., fewer gas-phase reactions.

Figure 14: CO abundance structure for both model sets, A and B, with different
localized depletion factors of small grains from left to right.Top row: CO
abundance structure for Model A, where CO has been depleted and the gas
surface density has a smooth profile.Bottom row: CO abundance structure for
Model B, where CO is considered to trace the gas surface density by a constant
scaling or CO/H2 ratio. In both models, there are not significant changes with
different degrees of small-grains depletion. However, we observe that in Model
A, CO has a more uniform vertical distribution than in Model B. Figure 15: C2H
2D abundance structure in each set of simulations.Top row: C2H abundance
structure for Model A.Bottom row: C2H abundance structure for Model B. In the
regions where the small dust is depleted, C2H can be produced closer to the
midplane. This causes the slight increases in C2H column density observed in
Fig. 9.

###  B.4 CO Emission Surfaces in AS 209

We show the difference between the CO emission layer for the $J=2-1$
transition in the AS 209 disk in Figure 16 from Teague et al. (2018) and Law
et al. (2021b). There are subtle differences between each surface. The surface
from Teague et al. (2018) is slightly above the Law et al. (2021b) one, but
they have a reasonable agreement considering uncertainties.

Figure 16: Emission heights found by Teague et al. (2018) and Law et al.
(2021b) for the CO $J=2-1$ emission coming from the AS 209 disk. These
surfaces were used to extract the Keplerian deviation in our models and
compare them with the observational ones.

##  Appendix C Keplerian deviations in the AS 209 gas kinematics

###  C.1 Keplerian deviations in the MAPS data

We show the kinematics deviations in the AS 209 disk using MAPS data at
0$\farcs$15 resolution in Figure 17 using methods described in Teague et al.
(2018). The data shows that the inferred deviation in the disk is at the
$\sim$ 1% level in small spatial scales in the inner 100 au. Therefore, there
is not a strong H2 depletion in the disk at 59 au; otherwise, there should be
a strong footprint in the kinematic deviations in the disk, as shown in Figure
7.

Figure 17: Kinematic deviations in the AS 209 disk using the CO data from MAPS
(Oberg et al., 2021). The deviations show variations $<1$% at the location of
the CO column-density gap, which is illustrated by the vertical dotted line.
The continuous profile is the expected profile through a convolution with a
0$\farcs$15 resolution ($\sim$ 18 au). The recovered deviations seem
consistent with the profile obtained by Teague et al. (2018) with lower-
resolution data, reassuring that the kinematic deviation in AS 209 is not
strong enough to cause a deep gas gap at 59 au.

###  C.2 Hydrostatic Equilibrium

We show the 2D fields of CO and the emitting heights in the hydrostatic
equilibrium runs in Figure 18. The emitting heights were taken at a different
layer considering the thermal structure of the disk and the changes associated
with the different geometry of the disk to match the CO $J=2-1$ line emission.
Moreover, the CO 2D abundance structure shows changes behind the dusty ring at
120 au. Those changes are probably associated with the changes in the vertical
distribution of the dust, changing the self-shielding in the disk beyond 120
au. Understanding the effect of hydrostatic equilibirum requires a deeper and
more sophisticated analysis that goes beyond the scope of this paper, but our
test proves that our results remain in those conditions.

Figure 18: Left: 2D CO abundances of the model with hydrostatic equilibrium.
Right: Layers at which we do the test for kinematic deviations. The thicker
and more flared disk geometry changes the CO abundance structure in the models
with less freeze-out in the midplane beyond 120 au.

## References

  * Aikawa et al. (2021) Aikawa, Y., Cataldi, G., Yamato, Y., et al. 2021, arXiv e-prints, arXiv:2109.06419.  https://arxiv.org/abs/2109.06419
  * Alarcón et al. (2020) Alarcón, F., Teague, R., Zhang, K., Bergin, E. A., & Barraza-Alfaro, M. 2020, ApJ, 905, 68, doi: 10.3847/1538-4357/abc1d6
  * Anderson et al. (2017) Anderson, D. E., Bergin, E. A., Blake, G. A., et al. 2017, ApJ, 845, 13, doi: 10.3847/1538-4357/aa7da1
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Armitage (2020) Armitage, P. J. 2020, Astrophysics of Planet Formation, 2nd edn. (Cambridge University Press), doi: 10.1017/9781108344227
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bae et al. (2017) Bae, J., Zhu, Z., & Hartmann, L. 2017, ApJ, 850, 201, doi: 10.3847/1538-4357/aa9705
  * Bergin et al. (2014) Bergin, E. A., Cleeves, L. I., Crockett, N., & Blake, G. A. 2014, Faraday Discussions, 168, 61, doi: 10.1039/C4FD00003J
  * Bergin et al. (2016) Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101, doi: 10.3847/0004-637X/831/1/101
  * Bergin & Williams (2018) Bergin, E. A., & Williams, J. P. 2018, arXiv e-prints, arXiv:1807.09631.  https://arxiv.org/abs/1807.09631
  * Bergner et al. (2019) Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2019, ApJ, 876, 25, doi: 10.3847/1538-4357/ab141e
  * Bergner et al. (2020) —. 2020, ApJ, 898, 97, doi: 10.3847/1538-4357/ab9e71
  * Bethell & Bergin (2011) Bethell, T. J., & Bergin, E. A. 2011, ApJ, 740, 7, doi: 10.1088/0004-637X/740/1/7
  * Birnstiel et al. (2010) Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79, doi: 10.1051/0004-6361/200913731
  * Birnstiel et al. (2012) Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148, doi: 10.1051/0004-6361/201118136
  * Bosman et al. (2021a) Bosman, A. D., Alarcon, F., Zhang, K., & Bergin, E. A. 2021a, arXiv e-prints, arXiv:2101.12502.  https://arxiv.org/abs/2101.12502
  * Bosman et al. (2021b) Bosman, A. D., Bergin, E. A., Loomis, R. A., et al. 2021b, arXiv e-prints, arXiv:2109.06223.  https://arxiv.org/abs/2109.06223
  * Bosman et al. (2021c) Bosman, A. D., Alarcón, F., Bergin, E. A., et al. 2021c, arXiv e-prints, arXiv:2109.06221.  https://arxiv.org/abs/2109.06221
  * Calahan et al. (2020) Calahan, J., Bergin, E., Zhang, K., et al. 2020, arXiv e-prints, arXiv:2012.05927.  https://arxiv.org/abs/2012.05927
  * Cieza et al. (2019) Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698, doi: 10.1093/mnras/sty2653
  * Cleeves (2016) Cleeves, L. I. 2016, ApJ, 816, L21, doi: 10.3847/2041-8205/816/2/L21
  * Cleeves et al. (2018) Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2018, ApJ, 865, 155, doi: 10.3847/1538-4357/aade96
  * Czekala et al. (2021) Czekala, I., Loomis, R. A., Teague, R., et al. 2021, arXiv e-prints, arXiv:2109.06188.  https://arxiv.org/abs/2109.06188
  * Dong et al. (2017) Dong, R., Li, S., Chiang, E., & Li, H. 2017, ApJ, 843, 127, doi: 10.3847/1538-4357/aa72f2
  * Dong et al. (2018) —. 2018, ApJ, 866, 110, doi: 10.3847/1538-4357/aadadd
  * Dong et al. (2015) Dong, R., Zhu, Z., & Whitney, B. 2015, ApJ, 809, 93, doi: 10.1088/0004-637X/809/1/93
  * Draine (2003) Draine, B. T. 2003, ARA&A, 41, 241, doi: 10.1146/annurev.astro.41.011802.094840
  * Du & Bergin (2014) Du, F., & Bergin, E. A. 2014, ApJ, 792, 2, doi: 10.1088/0004-637X/792/1/2
  * Du et al. (2017) Du, F., Bergin, E. A., Hogerheijde, M., et al. 2017, ApJ, 842, 98, doi: 10.3847/1538-4357/aa70ee
  * Dullemond et al. (2012) Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool.  http://ascl.net/1202.015
  * Dullemond et al. (2018) Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46, doi: 10.3847/2041-8213/aaf742
  * Favre et al. (2019) Favre, C., Fedele, D., Maud, L., et al. 2019, ApJ, 871, 107, doi: 10.3847/1538-4357/aaf80c
  * Fedele et al. (2018) Fedele, D., Tazzari, M., Booth, R., et al. 2018, A&A, 610, A24, doi: 10.1051/0004-6361/201731978
  * Fogel et al. (2011) Fogel, J. K. J., Bethell, T. J., Bergin, E. A., Calvet, N., & Semenov, D. 2011, ApJ, 726, 29, doi: 10.1088/0004-637X/726/1/29
  * Gail & Trieloff (2017) Gail, H.-P., & Trieloff, M. 2017, A&A, 606, A16, doi: 10.1051/0004-6361/201730480
  * Guzmán et al. (2018) Guzmán, V. V., Huang, J., Andrews, S. M., et al. 2018, ApJ, 869, L48, doi: 10.3847/2041-8213/aaedae
  * Guzmán et al. (2021) Guzmán, V. V., Bergner, J. B., Law, C. J., et al. 2021, arXiv e-prints, arXiv:2109.06391.  https://arxiv.org/abs/2109.06391
  * Hasegawa et al. (1992) Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167, doi: 10.1086/191713
  * Henning & Stognienko (1996) Henning, T., & Stognienko, R. 1996, A&A, 311, 291
  * Henning et al. (2010) Henning, T., Semenov, D., Guilloteau, S., et al. 2010, ApJ, 714, 1511, doi: 10.1088/0004-637X/714/2/1511
  * Hogerheijde et al. (2011) Hogerheijde, M. R., Bergin, E. A., Brinch, C., et al. 2011, Science, 334, 338, doi: 10.1126/science.1208931
  * Huang et al. (2016) Huang, J., Öberg, K. I., & Andrews, S. M. 2016, ApJ, 823, L18, doi: 10.3847/2041-8205/823/1/L18
  * Huang et al. (2018) Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42, doi: 10.3847/2041-8213/aaf740
  * Kamp et al. (2011) Kamp, I., Woitke, P., Pinte, C., et al. 2011, A&A, 532, A85, doi: 10.1051/0004-6361/201016399
  * Kanagawa et al. (2016) Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2016, PASJ, 68, 43, doi: 10.1093/pasj/psw037
  * Kanagawa et al. (2017) Kanagawa, K. D., Tanaka, H., Muto, T., & Tanigawa, T. 2017, PASJ, 69, 97, doi: 10.1093/pasj/psx114
  * Kastner et al. (2014) Kastner, J. H., Hily-Blant, P., Rodriguez, D. R., Punzi, K., & Forveille, T. 2014, ApJ, 793, 55, doi: 10.1088/0004-637X/793/1/55
  * Kastner et al. (2015) Kastner, J. H., Qi, C., Gorti, U., et al. 2015, ApJ, 806, 75, doi: 10.1088/0004-637X/806/1/75
  * Klarmann et al. (2018) Klarmann, L., Ormel, C. W., & Dominik, C. 2018, A&A, 618, L1, doi: 10.1051/0004-6361/201833719
  * Krijt et al. (2020) Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134, doi: 10.3847/1538-4357/aba75d
  * Law et al. (2021a) Law, C. J., Loomis, R. A., Teague, R., et al. 2021a, arXiv e-prints, arXiv:2109.06210.  https://arxiv.org/abs/2109.06210
  * Law et al. (2021b) Law, C. J., Teague, R., Loomis, R. A., et al. 2021b, arXiv e-prints, arXiv:2109.06217.  https://arxiv.org/abs/2109.06217
  * Long et al. (2018) Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17, doi: 10.3847/1538-4357/aae8e1
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603, doi: 10.1093/mnras/168.3.603
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425, doi: 10.1086/155591
  * McElroy et al. (2013) McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36, doi: 10.1051/0004-6361/201220465
  * Miotello et al. (2014) Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96, doi: 10.1051/0004-6361/201424712
  * Miotello et al. (2019) Miotello, A., Facchini, S., van Dishoeck, E. F., et al. 2019, A&A, 631, A69, doi: 10.1051/0004-6361/201935441
  * Molyarova et al. (2017) Molyarova, T., Akimkin, V., Semenov, D., et al. 2017, ApJ, 849, 130, doi: 10.3847/1538-4357/aa9227
  * Nieva & Przybilla (2012) Nieva, M. F., & Przybilla, N. 2012, A&A, 539, A143, doi: 10.1051/0004-6361/201118158
  * Öberg et al. (2009a) Öberg, K. I., Linnartz, H., Visser, R., & van Dishoeck, E. F. 2009a, ApJ, 693, 1209, doi: 10.1088/0004-637X/693/2/1209
  * Öberg et al. (2009b) Öberg, K. I., van Dishoeck, E. F., & Linnartz, H. 2009b, A&A, 496, 281, doi: 10.1051/0004-6361/200810207
  * Oberg et al. (2021) Oberg, K. I., Guzman, V. V., Walsh, C., et al. 2021, arXiv e-prints, arXiv:2109.06268.  https://arxiv.org/abs/2109.06268
  * Pinilla et al. (2012) Pinilla, P., Benisty, M., & Birnstiel, T. 2012, A&A, 545, A81, doi: 10.1051/0004-6361/201219315
  * Rab et al. (2020) Rab, C., Kamp, I., Dominik, C., et al. 2020, A&A, 642, A165, doi: 10.1051/0004-6361/202038712
  * Reboussin et al. (2015) Reboussin, L., Wakelam, V., Guilloteau, S., Hersant, F., & Dutrey, A. 2015, A&A, 579, A82, doi: 10.1051/0004-6361/201525885
  * Rosotti et al. (2020) Rosotti, G. P., Teague, R., Dullemond, C., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 495, 173, doi: 10.1093/mnras/staa1170
  * Schwarz et al. (2019) Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2019, ApJ, 877, 131, doi: 10.3847/1538-4357/ab1c5e
  * Smirnov-Pinchukov et al. (2020) Smirnov-Pinchukov, G. V., Semenov, D. A., Akimkin, V. V., & Henning, T. 2020, arXiv e-prints, arXiv:2009.09962.  https://arxiv.org/abs/2009.09962
  * Teague et al. (2018) Teague, R., Bae, J., Birnstiel, T., & Bergin, E. A. 2018, ApJ, 868, 113, doi: 10.3847/1538-4357/aae836
  * van der Marel et al. (2018) van der Marel, N., Williams, J. P., & Bruderer, S. 2018, ApJ, 867, L14, doi: 10.3847/2041-8213/aae88e
  * van Dishoeck et al. (2021) van Dishoeck, E. F., Kristensen, L. E., Mottram, J. C., et al. 2021, A&A, 648, A24, doi: 10.1051/0004-6361/202039084
  * Warren & Brandt (2008) Warren, S. G., & Brandt, R. E. 2008, Journal of Geophysical Research: Atmospheres, 113, doi: 10.1029/2007JD009744
  * Williams & Best (2014) Williams, J. P., & Best, W. M. J. 2014, ApJ, 788, 59, doi: 10.1088/0004-637X/788/1/59
  * Wilson (1999) Wilson, T. L. 1999, Reports on Progress in Physics, 62, 143, doi: 10.1088/0034-4885/62/2/002
  * Wilson & Rood (1994) Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191, doi: 10.1146/annurev.aa.32.090194.001203
  * Yu et al. (2017) Yu, M., Evans, Neal J., I., Dodson-Robinson, S. E., Willacy, K., & Turner, N. J. 2017, ApJ, 841, 39, doi: 10.3847/1538-4357/aa6e4c
  * Zhang et al. (2021) Zhang, K., Booth, A. S., Law, C. J., et al. 2021, arXiv e-prints, arXiv:2109.06233.  https://arxiv.org/abs/2109.06233
  * Zhang et al. (2018) Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47, doi: 10.3847/2041-8213/aaf744