# Precision measurement of transverse velocity distribution of a strontium
atomic beam

F. Gao1,2 H. Liu1,2 P. Xu1 X. Tian1,2 Y. Wang1 J. Ren1 Haibin Wu3
<EMAIL_ADDRESS>Hong Chang1,3<EMAIL_ADDRESS>1 CAS Key Laboratory
of Time and Frequency Primary Standards, National Time Service Center, Xi’an
710600, China  2 University of Chinese Academy of Sciences, Beijing 100049,
China  3 State Key Laboratory of Precision Spectroscopy, Department of
Physics, East China Normal University, Shanghai 200062, China

###### Abstract

We measure the transverse velocity distribution in a thermal Sr atomic beam
precisely by velocity-selective saturated fluorescence spectroscopy. The use
of an ultrastable laser system and the narrow intercombination transition line
of Sr atoms mean that the resolution of the measured velocity can reach 0.13
m/s, corresponding to 90$\mu K$ in energy units. The experimental results are
in very good agreement with the results of theoretical calculations. Based on
the spectroscopic techniques used here, the absolute frequency of the
intercombination transition of 88Sr is measured using an optical-frequency
comb generator referenced to the SI second through an H maser, and is given as
434 829 121 318(10) kHz.

Atomic (or molecular) beams, which have now become standard laboratory tools,
play very important roles in the field of atomic, molecular and optical
physics. These beams have been widely used in the determination of atomic
structures atomstructure , measurement of physical constants physicsconstant ,
studies of chemical reactions reaction  and atomic frequency standards
frequencystandard . In all these applications, measurement of the velocity
distribution of the atomic beams is both necessary and highly important. Many
spectroscopic techniques have been developed to perform these measurements
beam1 ; beam2 ; beam3 ; beam4 ; beam5 ; beam6 ; beam7 ; beam8 . However, the
precise determination of the velocity distribution of an atomic beam, and
especially that of a well collimated atomic beam, remains challenging because
of the low signal-to-noise ratio and relatively low spectroscopy resolution of
these beams.

In this paper, we report the precise measurement of the transverse velocity
distribution of a well-collimated thermal strontium beam by velocity-selective
saturated fluorescence spectroscopy. The measurement accuracy of the
transverse velocity of the atomic beam is greatly increased by the combination
of an ultrastable laser system and high-resolution spectroscopy of the
intercombination transition of 88Sr. The measurement resolution for the
velocity can reach 0.13 m/s, corresponding to 90 $\mu K$ in energy units.
Because of the confinement of the cylinder nozzle of the atomic beam, the
transverse velocity distribution is no longer the well-known Maxwell-Boltzmann
distribution. The measured results show good agreement with the theoretically
predicted results theory . Using the spectroscopic techniques developed here,
the absolute frequencies of the isotopes of strontium atoms are measured using
an optical frequency comb generator referenced to the SI second through an H
maser. Particular attention is paid to the intercombination transition of
88Sr. The frequency is measured as 434 829 121 318(10) kHz.

Figure 1: Experimental setup. (a) The narrowed linewidth of the laser system
(689 nm). PBS: polarized beam splitter, ULE: ultra-low expansion optical
cavity, PDH: Pound-Drever-Hall frequency locking system. (b) Experimental
setup used for velocity-selective saturated fluorescence spectroscopy. The
frequencies of the two counterpropagating laser beams are scanned using a saw-
tooth waveform at 0.1 Hz. PMF: single-mode polarization-maintaining fiber,
AOM: acousto-optic modulator, VCO: voltage-controlled oscillator, PSC: power
stabilization servo controller, M-S: Earth magnetic shield. (c) Optical
frequency measurement using a fiber frequency comb (Menlo FC1500), with
respect to an H maser at a radio frequency signal of 10 MHz.

Velocity-selective saturation (fluorescence) spectroscopy is a powerful
technique that is widely used for high resolution measurements. We use it here
to measure the velocity distribution of the atomic beam precisely. The
measurement principle is as follows. Consider two-level atoms interacting with
the counterpropagating probe beam and pump beam; the probability of the atoms
being in excited states is proportional to the velocity of the atoms, the
light intensities and the scattering rate of the excited state $\gamma$. When
the probe and pump beams have the same frequency, the group of atoms with zero
velocity see the light beams at resonance and show a Lamb dip at the atomic
resonant frequency with a large Doppler background. The peak amplitude can be
used to determine the number of atoms. When the frequencies of the pump and
probe lasers are different, only the atoms with velocity $v=\Delta\omega/(2k)$
see the light beams at resonance, where $\Delta\omega$ is the frequency
difference between the beams and $k$ is the wave vector of the laser light.
The spectroscopic resolution of this method is limited by the natural
linewidth of the excited state, and therefore the velocity resolution is
$v=(\gamma/f)c$, where $f$ and $c$ are the resonant frequency and the light
speed in a vacuum, respectively. The atom used here is strontium, which has
four stable natural isotopes, comprising bosonic 88Sr (82.58%), 86Sr (9.86%),
84Sr (0.56%), with nuclear spin I=0, and fermionic 87Sr (7.0%), with I=9/2.
The intercombination line $5^{1}S_{0}-5^{3}P_{1}$, because of its high
frequency fraction ($f/\gamma$), has been widely studied in optical frequency
standards frequency1 ; frequency3 ; frequency4 . In this paper, we use this
transition to precisely measure the velocity distribution of the atomic beam.

The experimental setup is shown in Fig. 1. A commercial extended cavity diode
laser (ECDL) (Topical-110) is used as the spectroscopy laser, and can
typically deliver 16 mW at 689 nm. The laser linewidth is reduced by locking
the laser to an ultrastable optical cavity via the Pound-Drever-Hall scheme;
the phase modulation is produced by a resonant electro-optic modulator (EOM),
which is driven at 5 MHz. The cavity is made of an ultralow expansion glass
with a finesse of 12000. To avoid a residual standing wave in the EOM, which
induces spurious amplitude modulation (AM) on the locking signal, a 60 dB
optical isolator is placed between the EOM and the cavity. The power of the
locking beam is maintained at as low a level as possible to minimize the
frequency shift caused by light heating and optical feedback. The linewidth
and the fractional frequency drift are about 200 Hz and $2.8\times 10^{-13}$
at 1 s, respectively.

The strontium atomic beam is obtained from strontium metal heated to 823 K
(550 oC) in an oven. A bundle of stainless steel capillaries (Monel-400) are
used to collimate the beam. The length and the internal diameter of the
capillary are 8 mm and 200 $\mu m$, respectively. The residual atomic beam
divergence is 25 mrad, and the typical atomic flux is estimated to be
$10^{12}/s$. The vacuum for the atomic beam is maintained at $1\times
10^{-8}$Torr with a 40 l/s ion pump.

Velocity-selective saturation (fluorescence) spectroscopy is obtained by the
use of two counterpropagating laser beams (i.e. the probe and pump beams)
perpendicular to the atomic beam. These beams are generated by two acousto-
optical modulators (AOMs) with the same double-pass configurations. The AOMs
are controlled by the same oscillator but with different voltage-controlled
oscillators (VCOs). A homemade circuit is used to prevent the oscillator
frequencies from disturbing each other. The frequency resolution is
approximately 0.1 Hz. The power is stabilized by a servo, which is better than
$10^{-4}$. Both the pump and probe beams have been expanded to have 1.2 cm
waists (1/$e^{2}$ diameter). The fluorescence signal is collected by a large
diameter lens in a direction that is orthogonal to the atomic beam and the
laser light beams on a high-sensitivity detector. A magnetic shield was placed
in the interaction region to prevent the Zeeman effect being caused by the
Earth’s magnetic field.

Figure 2: Velocity-selective saturated fluorescence spectra: (a)
$\Delta=0MHz$, corresponding to the condition where a group of atoms with
velocity $v=0$ is measured; (b) $\Delta=5MHz$, corresponding to the condition
where a group of atoms with velocity $v=3.45m/s$ is measured.

By scanning the pump and probe light beams at different frequencies, we
observe the velocity-selective saturated fluorescence spectra. Fig. 2 shows
two typical spectra for $\Delta=0$ and $\Delta=5\,MHz$, where $\Delta$ is the
frequency difference between the probe and pump light beams. The linewidth of
the central peak in Fig. 2(a) is approximately 180 kHz, which is larger than
its natural linewidth of 7.5 kHz. This effect mainly stems from power
broadening. The beam intensity of 800 $\mu W/cm^{2}$ (the saturation intensity
is 3 $\mu W/cm^{2}$) is chosen to produce a large signal-to-noise ratio for
the spectroscopy signal. The contributions of second-order Doppler broadening
(0.55 kHz) and transit time broadening (1.06 kHz) are small and can be
neglected. With this linewidth for the spectroscopy, the transverse velocity
of the magnitude at 0.13 $m/s$ (corresponding to 90 $\mu K$ in energy units)
can be measured.

Figure 3:  The transverse velocity distribution of the atomic beam. The black
dots represent the data for the experimental results. The error bars denote
statistical fluctuations that arise from measurement uncertainty for the
amplitudes of the signals. The red solid curve shows the results of the
theoretical calculations. All parameters are taken from the experimental
measurements. There are no free parameters.

For different frequency detunings of the pump and probe light beams, groups of
atoms with different velocities see the light beams at resonance. The number
of atoms can be measured from the amplitude of the saturated fluorescence
spectroscopy signal, and therefore the velocity distribution can be obtained,
and is shown in Fig. 3. The number of atoms has been normalized with respect
to the number of atoms when $\Delta=0$. The maximum frequency scanning range
for our AOMs is approximately $20\,MHz$, corresponding to $13.78\,m/s$ for the
maximum measurable transverse velocity of the atoms. Fig. 3 clearly shows that
the transverse velocity distribution of the thermal atomic beam is no longer a
Maxwell-Boltzmann-like distribution. The confinement of the capillaries leads
to this new distribution, which was predicted in Ref. theory  as

$P(v,a)=\frac{|v|\exp(-v^{2}/v_{0}^{2})\Gamma(-1/2,v^{2}L^{2}/v_{0}^{2}a^{2})}{\sqrt{8\pi
v_{0}^{2}(1+a^{2}/L^{2})^{1/2}}},$ (1)

where $v_{0}\equiv\sqrt{2k_{B}T/m}$ is the most probable velocity for the
atoms, $k_{B}$ is the Boltzmann constant, T is the oven temperature, $m$ is
the mass of the atoms, $a$ is the internal diameter of the collimator and $L$
is the length of the collimator. $\Gamma(-1/2,v^{2}L^{2}/v_{0}^{2}a^{2})$ is
the incomplete Gamma function. Based on the experimental parameters, the
results of the theoretical calculation of the velocity distribution are
plotted as the solid curve (red curve) in Fig. 3. There are no free
parameters. The theoretical curve clearly shows very good agreement with our
measurement results. The influence of the three other isotopes, 86Sr, 87Sr and
84Sr, on the measurement results is small because of large frequency
differences and the low natural abundance of these isotopes.

Figure 4: (a) Frequency measurement of the ${}^{88}Sr$ $5^{1}S_{0}\to
5^{3}P_{1}$ transition. The error bars correspond to the standard deviation
for each data set. (b) Comparison of the optical frequency measurement of the
intercombination line of ${}^{88}Sr$ with the results of previous measurements
Refs mea1 ; mea2 . The black triangle represents the data from Ref mea1 , the
red square represents the data from Ref mea2 , and the blue dot represents the
data from our measurement.  Table 1:  Optical frequency measurement results
and uncertainties for all four natural isotopes of strontium atoms.

Isotopes | $5^{1}S_{0}\to 5^{3}P_{1}$ | Frequency (kHz)
---|---|---
${}^{88}Sr$ |  $J=0$$\to J^{\prime}=1$ | 434 829 121 318 (10)
${}^{87}Sr$ |  $F=9/2$$\to F^{\prime}=7/2$ | 434 830 473 227 (45)
$F=9/2\to F^{\prime}=9/2$ | 434 829 342 995 (55)
$F=9/2\to F^{\prime}=11/2$ | 434 827 879 835 (50)
${}^{86}Sr$ |  $J=0$$\to J^{\prime}=1$ | 434 828 957 500 (15)
${}^{84}Sr$ |  $J=0$$\to J^{\prime}=1$ | 434 828 769 730 (100)

We reduce the power of the pump and probe light beams and use our narrowed
linewidth laser (200 Hz), and the saturation fluorescence spectroscopy
linewidth is approximately 55 kHz. The optical frequency of the Sr atoms is
measured using a commercial optical-frequency comb (Menlo FC1500). The
repetition rate and the carrier offset envelope frequency are locked to an H
maser. Fig. 4(a) shows the results of measurements of the 88Sr transition
frequency taken over a period of several days. The error bars correspond to
the standard deviation. The absolute frequency of the intercombination
transition of 88Sr is 434 829 121 318(10) kHz. A comparison with the results
of previous measurements in Ref mea1 ; mea2  is presented in Fig. 4(b), and
shows reasonable agreement. The optical frequencies of 86Sr, 87Sr and 84Sr are
measured using high laser intensities (2.1 $mW/cm^{2}$, 4.3 $mW/cm^{2}$, and
12 $mW/cm^{2}$, respectively); because of the low natural abundance of these
isotopes, the high power is required to obtain a good signal-to-noise ratio.
The measurement results for the optical frequencies of all four natural
isotopes of atomic strontium are summarized in Table I.

In conclusion, we report on the precise measurement of the transverse velocity
distribution of a well-collimated thermal Sr atomic beam using the velocity-
selective saturation fluorescence spectroscopy technique. The combination of
the ultrastable laser system and the narrow linewidth of the intercombination
transition means that the velocity measurement resolution is greatly improved.
The detectable minimum of the velocity is approximately $0.13\,m/s$
(corresponding to $90\,\mu K$ in energy units). The velocity distribution is
no longer likely to be a normal Maxwell-Boltzmann distribution. Instead, it
shows a counter-intuitive umbrella shape. The measured data are in very good
agreement with the results of the theoretical calculations. We also measured
the optical frequency of the $5^{1}S_{0}\to 5^{3}P_{1}$ transition through an
optical-frequency comb generator referenced to the SI second via an H maser.
The measured absolute frequency of the intercombination transition of 88Sr is
comparable to the results of previous measurements. Accurate optical frequency
values for all other isotopes are also presented, and these results may
provide a benchmark for subsequent measurements.

This research is supported by the National Natural Science Foundation of China
(Grant Nos. 11074252 and 61127901), and the Key Projects of the Chinese
Academy of Sciences (Grant No. KJZD-EW-W02).

## References

  * (1) X. Li, P. Mooney, S. Zheng, C.R. Booth, M.Braunfeld, S. Gubbens, D. A. Agard, and Y. Cheng, Nature Methods 10, 584 (2013).
  * (2) M. C. George, L. D. Lombardi, and E. A. Hessels, Phys. Rev. Lett. 87, 173002 (2001).
  * (3) P. Casavecchia, Rep. Prog. Phys. 63, 355 (2000)
  * (4) J. J. McFerran, J. G. Hartnett and A. N. Luiten, Appl. Phys. Lett. 95, 031103 (2009)
  * (5) R. C. Miller, and P. Kusch, Physical Review, 99,1314 (1955).
  * (6) N. F. Ramsey, Molecular Beams (Oxford University Press, London, 1956).
  * (7) H. Hellwig, S. Jarvis Jr, D. Halford, and H. E. Bell, Metrologia 9, 107 (1973).
  * (8) U. Brinkmann, J. Kluge and K. Pippert, J. Appl. Phys. 51, 4612 (1980).
  * (9) K.Bergman W.Demtroder P.Hering, Appl.Phys. 8, 65 (1975).
  * (10) G. Di Domenico, G. Mileti, and P. Thomann, Phys. Rev. A 64, 043408 (2001).
  * (11) W. F. Holmgren, I. Hromada, C. E. Klauss, and A. D. Cronin, New J. Phys. 13, 115007 (2011).
  * (12) G. Bannasch, J. Castro, P. McQuillen, T. Pohl, and T. C. Killian, Phys. Rev. Lett. 109, 185008 (2012).
  * (13) P. T. Greenland, M. A. Lauder, and D. J. H. Wort, J. Phys. D: Appl. Phys. 18, 1223, (1985).
  * (14) N. Poli, F. Y. Wang, M. G. Tarallo, A. Alberti, M. Prevedelli, and G. M. Tino, Phys. Rev. Lett. 1069, 038501 (2011).
  * (15) T. Ido, and H. Katori, Phys. Rev. Lett. 91, 053001 (2003).
  * (16) X. Xu, T. Loftus, J. Hall, A. Gallagher, and J. Ye, J. Opt. Soc. Am. B, 20, 968 (2003).
  * (17) G. Ferrari, P. Cancio, R. Drullinger, G. Giusfredi, N. Poli, M. Prevedelli, C. Toninelli, and G. M. Tino, Phys. Rev. Lett. 91, 243002 (2003).
  * (18) I. Courtillot, A. Quessada-Vial, A. Brusch, D. Kolker, G. D. Rovera, and P. Lemonde, Eur. Phys. J. D., 33, 161 (2005).


# The Carbon Isotopic Ratio and Planet Formation

Edwin A. Bergin  Arthur Bosman  Department of Astronomy, University of
Michigan, 1085 S. University Ave, Ann Arbor, MI 48109, USA  Richard Teague
Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts
Institute of Technology, Cambridge, MA 02139  Jenny Calahan  Department of
Astronomy, University of Michigan, 1085 S. University Ave, Ann Arbor, MI
48109, USA  Karen Willacy  Jet Propulsion Laboratory, California Institute of
Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA  L. Ilsedore Cleeves
Department of Astronomy, University of Virginia, Charlottesville, VA 22904
Kamber Schwarz  Max Planck Institute for Astronomy, Königstuhl 17, Heidelberg,
Germany  Ke Zhang  Department of Astronomy, University of Wisconsin-Madison,
475 N. Charter St., Madison, WI 53706  Simon Bruderer  Max-Planck-Institut für
Extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, Germany

###### Abstract

We present the first detection of 13CCH in a protoplanetary disk (TW Hya).
Using observations of C2H we measure CCH/13CCH = 65$\pm$20 gas with a CO
isotopic ratio of 12CO/13CO = 21$\pm$5 (Yoshida et al., 2022a). The TW Hya
disk exhibits a gas phase C/O that exceeds unity and C2H is the tracer of this
excess carbon. We confirm that the TW Hya gaseous disk exhibits two separate
carbon isotopic reservoirs as noted previously (Yoshida et al., 2022a). We
explore two theoretical solutions for the development of this dichotomy. One
model represents TW Hya today with a protoplanetary disk exposed to a cosmic
ray ionization rate that is below interstellar as consistent with current
estimates. We find that this model does not have sufficient ionization in cold
(T $<$ 40 K) layers to activate carbon isotopic fractionation. The second
model investigates a younger TW Hya protostellar disk exposed to an
interstellar cosmic ray ionization rate. We find that the younger model has
sources of ionization deeper in a colder disk that generates two independent
isotopic reservoirs. One reservoir is 12C-enriched carried by
methane/hydrocarbon ices and the other is 13C-enriched carried by gaseous CO.
The former potentially provides a source of methane/hydrocarbon ices to power
the chemistry that generates the anomalously strong C2H emission in this (and
other) disk systems in later stages. The latter provides a source of gaseous
13C rich material to generate isotopic enrichments in forming giant planets as
recently detected in the super-Jupiter TYC 8998-760-1 b by Zhang et al.
(2021).

editorials, notices — miscellaneous — catalogs — surveys

††journal: ApJ

##  1 Introduction

The exoplanet revolution is now heading in a phase of greater characterization
where the composition of planetary atmospheres can be explored in greater
detail via high resolution spectroscopy (e.g., Line et al., 2021; Zhang et
al., 2021; Brogi & Birkby, 2021) or observations with new more sensitive
instruments (e.g., JWST; JWST Transiting Exoplanet Community Early Release
Science Team et al., 2023). A greater understanding of composition allows for
an exploration of how the formation conditions of the planet might be
reflected in its composition. Over the past decade the strongest link has been
to explore the bulk elemental abundance ratio of carbon to oxygen or C/O.
Öberg et al. (2011b) first presented this idea utilizing basic chemistry
(e.g., deposition/sublimation of primary molecular carriers of C and O) to
show that the gas phase and solid-state C/O ratios are predicted to change
with distance from the star and it might be possible to tie composition to the
formation location. This may be more complicated as planets migrate, can
accrete icy planetesimals/pebbles, and might have core-atmosphere mixing along
with gravitational settling (Cridland et al., 2016; Helled & Guillot, 2017;
Guillot et al., 2023). However, it captures a central element and is now
widely compared to exoplanet atmospheric composition retrievals (Barman et
al., 2015; Lavie et al., 2017; Oreshenko et al., 2017, to list a few).

Beyond the C/O ratio, the carbon isotopic ratio presents another interesting
avenue for exploring the links between disk and exoplanet composition. This
has been enabled by advances in high resolution spectroscopy of exoplanets
that have isolated the 12C/13C ratio in two systems. In one case for a young,
accreting super-Jupiter TYC 8998-760-1 b at 160 au, Zhang et al. (2021),
measure a ratio of 12C/13C = 31${}^{+17}_{-10}$ (90% confidence). This is a
significant 13C enrichment when compared to the local interstellar value of 68
(Langer & Penzias, 1993; Milam et al., 2005). A similar level of enrichment is
found towards the Hot Jupiter WASP-77Ab, 12C/13C = 10.2-42.6 at 68% confidence
(Line et al., 2021). This enrichment is posited to have originated in
molecular ices (Zhang et al., 2021). Curiously, the solar system betrays no
evidence of 13C enrichments at this level in any body, including comets that
are comprised of molecular ices (Nomura et al., 2023). Further, interstellar
ices towards quiescent lines of sight also show little evidence of enrichment
(McClure et al., 2023). This points towards later phases, i.e., disk chemical
evolution.

If molecular ices are enriched in 13C within the disk, we should see some
evidence of this in disk systems as the gas should exhibit the opposite (i.e.,
a depletion), provided the fractionation originates in the disk. The state-of-
the-art of these measurements can be found in TW Hya, the nearest young gas-
rich disk, and the best characterized system. Yoshida et al. (2022a) used the
line wing emission of optically thick CO isotopologues to derive the carbon
isotopic ratio in gaseous CO as a function of position. They find that 13C
isotopic enrichments are present in disk surface layers from $\sim$70-100 au.
The detection of 13C18O inside the CO snowline implies 12CO/13CO $\sim
40^{+9}_{-6}$ in sublimated CO ices (Zhang et al., 2017). In contrast, the
outer ($>$100 au) disk shows the opposite signature and appears to be depleted
in 13CO (Yoshida et al., 2022a).

However, the C/O ratio in TW Hya is measured to be $>$ 1 (Bergin et al., 2016;
Kama et al., 2016a). Thus, CO is not the only reservoir of carbon. Here we
report the detection of weak 13CCH emission within the TW Hya disk. The
chemistry of CCH is the tracer of the excess carbon not carried by CO
providing accurate measure of its isotopic ratio. To explore the origin of
isotopic enhancements considering the full extent of isotopic measurements we
use models of carbon fractionation of TW Hya today (a later stage few Myr
disk, i.e. class II) and in an early-stage (i.e. less evolved; class I) state.
In §2 we will present our observations and the technique used to obtain the
weak signal. §3 presents the derived isotopic ratio and the chemical models
that are used to explore fractionation in the disk environment. Finally, in §4
we discuss the implications of this result and our conclusions.

##  2 Observations

###  2.1 C2H

The C2H data we use was originally published in Bergin et al. (2016) using the
data from ALMA project 2013.1.00198.S. We refer the reader to the original
article for details on the data reduction and imaging. This included the
subtraction of the continuum by modeling the continuum contribution as a
linear component around the line emission and using the uvcontsub task to
remove this in the image plane. Note that as the continuum emission is
confined to inside of ${\sim}~{}1\arcsec$ (${\sim}~{}60$ au), the continuum
subtraction should have minimal effect on the line data presented here. The
resulting images have a beam of $0.50\arcsec\times 0.41\arcsec$ with a
position angle of $59.0\arcdeg$.

Figure 1 shows the integrated flux density for the $N=4-3$ $J=9/2-7/2$ $F=5-4$
transition made using a Keplerian mask assuming $i=6.8\arcdeg$, ${\rm
PA}=151\arcdeg$, $M_{\rm star}=0.6~{}M_{\rm sun}$ and $d=60.1~{}{\rm pc}$. The
mask was used rather than $\sigma$-clipping to remove any contribution from
the $F=4-3$ transition. The radial profile is well described by a Gaussian
centered at $1.1\arcsec$ with a width of $0.58\arcsec$ (Bergin et al., 2016).

Integrating over the whole disk yields a total integrated flux of $6.17\pm
1.16~{}{\rm Jy~{}km\,s^{-1}}$ where the uncertainties have been calculated
following Teague et al. (2019), while over the annulus with an inner radius of
$0.52\arcsec$ (31 au) and an out radius of $1.68\arcsec$ (101 au) we find
$4.20\pm 0.68~{}{\rm Jy~{}km\,s^{-1}}$. The peak intensity was $372\pm
6~{}{\rm mJy\,beam^{-1}}$.111There is a difference in this intensity and that
given in Table 1 of Bergin et al. (2016) where the peak intensity for this
transition is listed as $1.100\pm 0.007~{}{\rm Jy~{}beam^{-1}}$. In this work
we deblended the F = 5–4 and 4–3 which was not done by Bergin et al. (2016).

Figure 1: Zeroth moment of the $N=4-3$ $J=9/2-7/2$ $F=5-4$ transition of C2H
made using a Keplerian mask. The peak value is $148.7~{}{\rm
mJy~{}beam^{-1}~{}km\,s^{-1}}$. The beam size is shown in the bottom left of
the panel.

The non-Gaussian beam which results from the imaging process can introduce
subtle inaccuracies when calculated integrated fluxes. Following Czekala et
al. (2021), the non-Gaussianity of the beam can be quantified with $\epsilon$,
which is the ratio of the volumes of the CLEAN and dirty beams, where
$\epsilon=1$ indicates a purely Gaussian beam. Using the scripts released as
part of the MAPS Large Program (Öberg et al., 2021) to apply a correction to
this (the so-called ‘JvM correction’), the C2H data were found to have
$\epsilon=0.69$. Calculating the integrated flux density and peak intensity
with the JvM-corrected data found discrepancies of 5% and 2%, respectively.
These uncertainties were incorporated in the ratios described below by adding
this uncertainty in quadrature to the statistical uncertainty.

###  2.2 13CCH

The 13CCH line observations of TW Hya were carried out by ALMA on November 23,
2016 (2015.1.00308.S). The array included 42 antennas, with baselines lengths
spanning from 15 to 704 m. The correlators were set to have three-line windows
centered on 336.571, 351.060, and 349.346 GHz, and each window had a bandwidth
of 467 MHz. A fourth spectral window was dedicated to continuum observations,
centering on 338.769 GHz with a band width of 2 GHz. The Nearby Quasars
J1037-2934, J1058+0133, J1107+4449 were used for gain, bandpass, and flux
calibration, respectively. The total on-source time was 49 minutes. The raw
visibility data were calibrated by NRAO staff members in CASA version 4.7.0,
and we further applied phase-only self-calibration on the data using the
continuum spectral window and line-free channels in all three-line spectral
windows. The fully calibrated visibilities were then Fourier inverted to
generate images by using the CLEAN task with a Briggs weight of 0.5. The final
13CCH image has a noise level of $3.6~{}{\rm mJy\,beam^{-1}}$ per $210~{}{\rm
m\,s^{-1}}$ channel and a beam size of $0.32\arcsec\times 0.30\arcsec$ with a
position angle of $4.4\arcdeg$.

To extract a disk-averaged spectrum we use the Python package GoFish (Teague,
2019) which implements the method presented in Yen et al. (2016). First the
spectrum in each pixel was shifted to the systemic velocity by applying a
shift of $\delta v(r,\,\phi)=-v_{\rm kep}(r)\,\cos\phi\,\sin i$, where $v_{\rm
kep}(r)$ is the Keplerian rotation at radius $r$, $\phi$ is the polar angle of
the pixel, measured relative to the red-shifted major-axis, and $i$ is the
disk inclination. With all spectra centered on the systemic velocity, they are
stacked to increase the signal-to-noise of the spectrum. This method will
result in the average spectrum in a pixel within the averaged area. To obtain
an integrated flux density, this average spectrum is scaled by the total area
of the annulus over which the spectra were averaged. This approach of first
averaging and then scaling is mathematically the same as integrated over the
entire area, however the intermediate averaged spectrum allows for a more
accurate determination of the integration boundaries.

For the velocity shifting we assume the same parameters as for the Keplerian
mask in the previous section, $i=6.8\arcdeg$, ${\rm PA}=151\arcdeg$, $M_{\rm
star}=0.6~{}M_{\rm sun}$ and $d=60.1~{}{\rm pc}$. Assuming that the radial
emission profile of 13CCH follows that of C2H, we average from $0.52\arcsec$
(31 au) to $1.68\arcsec$ (101 au) to increase the signal-to-noise ratio of the
detection, weighting each individual annulus by their respective area. As each
annulus has a fixed width of 1/4 of the beam FWHM (${\sim}~{}80~{}{\rm mas}$),
the weights are simply proportional to $r^{2}$.

To verify that the image was centered, we vary the source center by small
pixel-scale offsets and measure the signal-to-noise ratio of spectrum. We find
a peak signal-to-noise with an offset of $(\delta x_{0},\,\delta
y_{0})=(-52~{}{\rm mas},-155~{}{\rm mas})$.

The resulting spectrum is shown in Fig. 2, showing the clear detection of the
F2 = 11/2 - 9/2 and F2 = 9/2 - 7/2 components of the N = 4 - 3, J = 9/2 - 7/2
F1 = 5 - 4 transition. Integrating this spectrum yields an integrated flux of
$56.9\pm 14.1~{}{\rm mJy~{}km\,s^{-1}}$. We have verified by changing the
inner and outer radii of the annuli we integrate over that this range gives
the highest signal-to-noise detection of the lines. Integrating over the whole
disk finds an integrated flux of $70.2\pm 23.2~{}{\rm mJy~{}km\,s^{-1}}$. From
Gaussian fits to the line, the peak intensity values for the bright annulus
region were $0.9\pm 0.3~{}{\rm mJy~{}beam^{-1}}$ and $1.2\pm 0.3~{}{\rm
mJy~{}beam^{-1}}$ for the F2 = 11/2 - 9/2 and F2 = 9/2 - 7/2 components,
respectively.

Figure 2: Integrated average flux density between 0.8″ and 1.8″, resampled to
a velocity resolution of 150 m s-1, clearly showing the F = 11/2 - 9/2 and F =
9/2 - 7/2 components of the N = 4 - 3, J = 9/2 - 7/2 transition. The black
lines show two Gaussian fits to the spectrum with a shared line width of 125 m
s-1.

To verify that the non-Gaussianity of the beam did not bias the results, the
JvM-correction (Czekala et al., 2021) was applied to the 13CCH data. These
data were found to have $\epsilon=0.70$, consistent with the C2H data, which
is unsurprising given they use the same antenna configuration. The integrated
flux density was not found to substantially vary compared to the native data,
finding only a 1% difference. The peak intensity varied more, up to 40%,
however this is because such a low-level signal is dominated by the noise
which is substantially suppressed with the JvM-correction. As such, this term
is ignored as the difference is equivalent to the difference in the noise
levels determined by the RMS of a line-free channel. The 1% uncertainty on the
integrated flux density is added in quadrature to the statistical
uncertainties.

##  3 Carbon Isotopic Ratio in TW Hya

###  3.1 12C12H/13C12H in TW Hya

Using the integrated intensities from above, we find a C2H / 13CCH ratio of
$74\pm 22$ in the bright annulus and $88\pm 3\mathbf{4}$ averaged over the
whole disk. Similarly, if we use the peak flux density for each line, after
correcting for the relative beam sizes and assuming that they’re optically
thin so that the two F2 components can be added, we find a C2H / 13CCH ratio
of $78\pm 15$, consistent with the integrated flux approach. The above ratios
are ratios in flux and not total column. Based on the fact that the flux ratio
is consistent with the interstellar 12C/13C isotopic ratio ($\sim$68 Langer &
Penzias, 1993; Milam et al., 2005) and hyperfine ratios (Bergin et al., 2016),
the C2H emission in TW Hya is optically thin. In this limit, for observations
that have been corrected for differences in the beam size (as above), the
integrated flux density ratio is essentially equivalent to the overall column
density ratio correcting for an 11% difference in the A-coefficents. We
therefore conclude that N(CCH)/N(13CCH) = 65$\pm$20 over spatial scales
commensurate with the TW Hya C2H ring.

Figure 3: Measurements of isotopic ratios in TW Hya placed into observational
context in terms of distance from the star. Measurements are from Yoshida et
al. (2022a, CO), Zhang et al. (2017, C18O), Hily-Blant et al. (2019, HCN), and
CCH (presented here). For illustrative purposes, we also include the
measurement of TYC 8998 b (Zhang et al., 2021) which is shown at the orbital
distance estimated from its host star.

###  3.2 Two Carbon Isotope Reservoirs

Fig. 3 provides the landscape of isotopic ratio measurements within TW Hya,
with explicit values given in Table 1. Inside of the CO snowline at 21 au
(Schwarz et al., 2016) an isotopic ratio of $\sim$40${}^{+9}_{-6}$ is measured
via detection of ^13C^18O by Zhang et al. (2017). Hily-Blant et al. (2019)
observe HCN and H^13CN and measure an average ratio of 86$\pm 4$ inside 60 au.
The differences between CO, HCN, and CN provide evidence of multiple isotopic
reservoirs in the system as noted by Hily-Blant et al. (2019) and Yoshida et
al. (2024). This is further traced when comparing the 12C/13C as estimated for
CO and CCH. In TW Hya CCH emits in a ring at radii of $\sim$40-110 au (Fig. 1;
Bergin et al., 2016) and our detection of ^13CCH must arise from within the
ring where we measure a ratio of 12C/13C in CCH of 65$\pm$15\. In contrast,
over a similar spatial extent, Yoshida et al. (2022a) estimate a 12C/13C ratio
of 21$\pm$5 in CO. If the isotopic ratio of 21 existed in CCH in our data then
we would have a 14$\sigma$ detection 13CCH in the integrated emission, which
is in contrast to the current 4$\sigma$ detection. This re-enforces the
evidence that there are two isotopic reservoirs for carbon present in TW Hya
as argued previously by Yoshida et al. (2022a) for carbon. We note two
isotopic reservoirs have also been isolated for nitrogen by Hily-Blant et al.
(2019).

Table 1:  12C/13C in TW Hya  Location (au) | Molecular Carrier |  12C/13C Ratio | Reference
---|---|---|---
0 – 20 |  ^12C^18O/^13C^18O | 40${}^{+9}_{-6}$ | Zhang et al. (2017)
$<$ 60 |  H^12CN/H^13CN | 86${}^{+4}_{-4}$ | Hily-Blant et al. (2019)
25 – 55 |  ^12CN/^13CN | 70${}^{+8}_{-6}$ | Yoshida et al. (2024)
70 – 110 |  ^12CO/^13CO | 21${}^{+5}_{-5}$ | Yoshida et al. (2022a)
40 – 110 |  ^12CCH/^13CCH | 65${}^{+20}_{-20}$ | This Paper
$>$ 130 |  ^12CO/^13CO |  $>$ 84 | Yoshida et al. (2022a)

##  4 Model of Disk Carbon Isotopic Chemistry

To examine the origins of two isotopic reservoirs for carbon in the TW Hya
disk we use a detailed model of the carbon isotopic chemistry for TW Hya. For
this purpose, we use the Dust And LInes thermochemical model (DALI) code
(Bruderer et al., 2012; Bruderer, 2013). Our aim is to explore a generic model
of the chemistry of carbon isotopic fractionation in two frameworks. One is an
extant model of the TW Hya disk that has been used in previous modeling
efforts (Trapman et al., 2017; Bosman & Banzatti, 2019). We call this model
the “late-stage disk” model as it is motivated by the TW Hya Class II disk as
it is today with numerous observational constraints. In addition, we present
an “early-stage” model that explores the origins of isotopic chemistry during
earlier phases more consistent with a Class I disk. In our efforts we will not
fit the overall emission levels; rather, we will explore whether two separate
carbon isotopic reservoirs are created. That is, we are not attempting to
fully match this system but rather attempting to explore the conditions that
may lead to carbon isotopic fractionation.

We use the extant TW Hya model published by Trapman et al. (2017) and Bosman &
Banzatti (2019). The basic parameters of the model are given in the appendix.
The chemical network including isotopic fractionation is created by Miotello
et al. (2014) and includes key fractionation reactions linking 12C and 13C
(Langer et al., 1984; Röllig & Ossenkopf, 2013). The critical reactions in
this case are the primary fractionation reaction:

${\rm{}^{13}C^{+}+CO\leftrightarrow C^{+}+^{13}CO+\Delta E=35~{}K},$ (1)

which favors 13CO at low ($\lesssim$40 K) temperature. The other key isotopic
fractionation pathway in the disk is isotopic selective photodissociation of
CO and we use the rate prescription of Visser et al. (2009). For a basic
introduction to carbon isotopic fractionation in the disk framework we refer
the reader to Woods & Willacy (2009) and to the review by Nomura et al.
(2023). In this model we include a grain surface chemistry that is limited to
simple hydrogenation reactions to make saturated species (e.g. water, ammonia,
methane). Further, the overall network does not extend to methanol.

###  4.1 Late-Stage Model

Here, we use the existing model of TW Hya by Trapman et al. (2017) and Bosman
& Banzatti (2019). This model has been designed to match several well
characterized aspects of this disk. First, the strong mm-continuum image and
scattered light surface requires grain growth and settling (van Boekel et al.,
2017). Further, the mm-continuum emission shows that the large grains are more
spatially confined compared to the gas (Andrews et al., 2016; Huang et al.,
2018). We therefore, assume the disk is several Myr-old (Sokal et al., 2018)
with 99% of the grain mass confined to the midplane to match the mm-wave dust
emission distribution.

The primary ionization sources towards this disk (ultraviolet, X-ray, cosmic
rays) have been constrained via observations. For thoroughness the UV field is
the stellar spectrum of TW Hya with UV excess from Cleeves et al. (2015). We
adopt LX = 1.4 $\times$ 1030 ergs s-1 with a plasma temperature of 3 MK
(Stelzer & Schmitt, 2004). The cosmic ray ionization rate ($\zeta_{cr}$) is
assumed to be 10-19 s-1 which is constrained via chemical analysis of
molecular ion emission (Cleeves et al., 2015). Additional information is found
in the original reference for this model (Trapman et al., 2017).

Analysis of numerous disk CO isotopologue observations, starting with TW Hya
(Favre et al., 2013), have found that the abundance of CO appears to be
reduced in disk gas (see discussion in Miotello et al., 2023). Essentially, CO
is believed to be provided to the disk with an “interstellar” abundance
relative to H2 $\sim$ 10-4 (Bergin & Williams, 2017). Some regions of the disk
have temperatures below the CO sublimation temperature ($\sim$20 K) and CO ice
forms. However, in regions where the dust temperature is $>$ 20 K, the
expectation is that the CO gaseous abundance will be interstellar. However, in
systems where strong constraints on the H2 mass exist, such as TW Hya, the CO
gas phase abundance is found to be depleted, i.e. [C/H]gas $<$ [C/H]⋆ (Favre
et al., 2013; Kama et al., 2016b; Zhang et al., 2017, 2019). This has now been
confirmed via a direct measurement of the H2 density through pressure
broadened line wings (Yoshida et al., 2022b). The depletion of CO (and also
water vapor; Du et al., 2017) is believed to occur via chemical processing of
CO into species with lower volatility, such as CO2 or CH3OH, and via
interactions with grain evolution (see Krijt et al., 2020, for comprehensive
model). In sum, the Trapman et al. (2017) model assumes that overall gas-phase
carbon is depleted relative to ISM levels, with the missing carbon present as
ices on the mm-grains that are settled to the midplane and more spatially
confined than the gas.

Table 2: Initial Chemical Abundancesa Species | Late Stage | Early Stage
---|---|---
H2 | 0.5 | 0.5
He | 7.6 $\times$ 10-2 | 7.6 $\times$ 10-2
Mg | 4.2 $\times$ 10-10 | 4.2 $\times$ 10-10
Si | 7.9 $\times$ 10-9 | 7.9 $\times$ 10-9
S | 1.9 $\times$ 10-9 | 1.9 $\times$ 10-9
Fe | 4.3 $\times$ 10-10 | 4.3 $\times$ 10-10
CH4 | 2.7 $\times$ 10-6 | $\cdots$
13CH4 | 3.5 $\times$ 10-8 | $\cdots$
CO | 2.7 $\times$ 10-6 | 2.7 $\times$ 10-4
13CO | 3.5 $\times$ 10-8 | 3.5 $\times$ 10-6
N2 | 1.1 $\times$ 10-5 | 1.1 $\times$ 10-5
N15N | 2.4 $\times$ 10-8 | 2.4 $\times$ 10-8
H2Oice | 2.6 $\times$ 10-6 | 2.6 $\times$ 10-6
aRelative to H

Finally, we assume C/O $\sim$2 based on the analysis from C2H (Bergin et al.,
2016; Kama et al., 2016a). This model assumes that 50% of the remaining
gaseous carbon is in CO and 50% originally in CH4. The presence of the CH4 has
been previously shown to generate strong C2H emission in disks (Bosman et al.,
2021a).222Specifically this paper states that the origin of C2H arises from
CH4 released from carbon in refractory organics. However, this is simulated in
the system by putting carbon in CH4. Our model is run for 1 Myr, which is much
less than the estimated age of TW Hya ($\sim$8 Myr; Sokal et al., 2018). This
difference is because we are not modeling the CO abundance evolution, which we
assume took place during earlier stages. Rather our goal is to explore whether
existing TW Hya chemical models can readily fractionate carbon carriers to
create two interdependent and stable isotopic reservoirs. Table 2 provides the
overall initial chemical abundances.

Figure 4:  Predicted vertical distribution of CO, CCH, and CH4 ice isotopic
ratios in the “late-stage disk” model at a radial location commensurate with
the center of the C2H emission ring (63 au). (a) Overall gas density structure
and modeled vertical gas thermal profile at this radius. (b,c,d) Predicted
vertical distribution of CO/13CO, C2H/13CCH, and CH4,ice/13CH4,ice from the
thermochemical model. In these panels the chemical abundances are referenced
to the left axis while the dashed line provides the predicted isotopic ratio
that is referenced to the axis on the right. For all species the initial
isotopic ratio is 77 which is shown as the dotted grey line.

With this calibrated model we explored whether this physical perspective could
produce two separate isotopic reservoirs. For this purpose, we include a
network that includes fractionation of carbon (Miotello et al., 2014). The
results are shown in Fig. 4. In this Figure, panel (a) provides the vertical
density and thermal profile. In panels (b-d) we show model predictions as a
vertical cut taken at 63 au. This location traces the C2H ring and is
commensurate with where C2H and CO appear to trace two separate isotopic
reservoirs (Fig. 3). The figure shows the C2H, CO, and CH4 vertical
distributions in abundance alongside the 13C isotopologues. In our models the
initial isotopic ratio in both CO and CH4 is 77 and we seek to determine if
the overall isotopic chemistry produces two separate reservoirs. The full
chemical evolution from this model with the major carriers shown is given in
Fig. 5.

Exploring panel (b) in Fig. 4 we find an interesting fractionation pattern for
CO. At the upper surface, near z = 25 au, the effects of CO isotopic selective
photodissociation are found with a deficit in 13CO. Deeper in the disk there
is a dip in the 13CO abundance near z = 17 au that produces a sharp increase
in 12CO/13CO over a short vertical distance. The gaseous CO isotopic ratio
estimated from column density ratios at 63 au is 79 which is well above the
observed value of 21.

Figure 5: Plot of (top) major carriers of 12C and (bottom) 13C. Abundances are
normalized to the total amount of isotopic carbon with $x$(12C) = $5.2\times
10^{-6}$ and $x$(13C) = $7.0\times 10^{-8}$ (abundances given relative to
total H). This model is sampled at r=63 au, the location of the C2H emission
ring. We note that HC13CH2 is not shown here but has the same distribution and
abundance as H13CCH2.

In layers above the midplane (4 au $<$ z $<$ 20 au) where CH4 ice sublimates,
much of the excess carbon provided by CH4 is processed into hydrocarbon ices
(C2H2 and HCCH2; Fig. 5). Above the edge of this vertical zone near 20 au
gaseous C2H production peaks. For C2H, the model predicts a complex
fractionation pattern. However, it is important to focus on vertical layers
where C2H is found in abundance between 14 au $<$ z $<$ 24 au. Here, some of
the excess carbon resulting from isotopic selective photodissociation of CO
finds its way into C2H raising the 13CCH abundance above the initial value,
with a sharp decrease in CCH/13CCH at z = 18 au as the 13C finds its way into
hydrocarbon ices. The gaseous C2H isotopic ratio estimated from column density
ratios is 78 which agrees with observations (Table 1). The predicted total C2H
column density at this radius is 6 $\times$ 1014 cm-2. This would give an
integrated flux density for the N=4–3 J=9/2–7/2 F-5-4 transition of $\sim$0.2
Jy km/s in a 0.45′′ beam at 60 pc (assuming a gas temperature of 40 K). This
is close to the observed peak integrated flux density within the emission ring
(see Fig. 1). We also confirm that the gas phase C/O ratio in the model is 2
in the C2H emitting layer, which is consistent with previous work (Bergin et
al., 2016).

In summary, there is a narrow range of physical parameters where 13C
fractionation is found, which limits the impact over the entire column. This
particular “late-stage disk” model is incapable of producing two separate
isotopic ratios for CO and C2H (and HCN). As a final check, to explore whether
the excess CH4 in the gas could hide the presence of isotopic fractionation,
we ran an identical model with no excess CH4 as an initial condition. This
model did not produce significant carbon isotopic fractionation of CO or C2H.
For this set of initial conditions, the C/O ratio is effectively solar at
these distances (i.e. C/O = 1) and, as a result, the C2H column was reduced by
3 orders of magnitude. This is consistent with earlier analyses (Bergin et
al., 2016; Cleeves et al., 2018; Bosman et al., 2021b).

###  4.2 Early-Stage Disk Model

####  4.2.1 Motivation

Our model of TW Hya as it exists today is incapable of producing two separate
carbon isotopic reservoirs. We therefore also explore a model of a younger
disk where conditions might be significantly different. Our hypothesis is as
follows. Numerous works have shown that, if cosmic rays are present, CO can be
chemically processed into a variety of less volatile species which may relate
to the reduced abundance of CO as seen in numerous disk systems (Bergin et
al., 2014; Furuya & Aikawa, 2014; Reboussin et al., 2015; Eistrup et al.,
2017; Schwarz et al., 2018; Bosman et al., 2018). Perhaps this process might
also lead to isotopic fractionation of CO and thereby simultaneously set the
stage where another carbon-bearing molecule would carry the opposite isotopic
signature. If this carbon-bearing molecule were a hydrocarbon then it might
provide the excess carbon needed to generate strong C2H emission at later
stages.

The chemical processing of CO appears to be time dependent as two independent
analyses of C18O and lesser abundant isotopologues find a decline in the
gaseous CO abundance from the initial interstellar level as a function of
evolutionary state corresponding to $\sim$1 Myr timescales (Bergner et al.,
2020; Zhang et al., 2020). Along these lines the gaseous C/O ratio might also
change as the disk evolves. Interstellar ices have some CO and CH4 (Öberg et
al., 2011a). These ices would sublimate at 30 K (Minissale et al., 2022) to
mix with CO gas that has C/O = 1. However, interstellar CH4 ice is only
$\sim$10% of the CO ice content (McClure et al., 2023). When both sublimate to
the gas at 30-40 K (with water and CO2 frozen), the resulting composition is
effectively C/O $\sim$ 1\. Given that observational analyses are suggesting
that TW Hya has C/O $\sim$ 1.5-2 there also must be some additional chemical
evolution in the volatile inventory and its C/O ratio.

In this light, cosmic rays are known to be present in well characterized pre-
stellar cores with measured ionization rates consistent with expectations for
the dense interstellar medium (Maret & Bergin, 2007; Redaelli et al., 2021).
In contrast, Class II systems appear to have reduced ionization that are
argued to be the result of cosmic ray exclusion by magnetized winds (analogous
to the solar system’s heliosphere) or tangled magnetic fields within the disk
environment (Cleeves et al., 2015; Aikawa et al., 2021; Seifert et al., 2021).
Class I systems still have residual envelopes that limit the reach of
magnetized stellar winds, restricting their effect to the extent of the cavity
carved by the protostar’s energetic flows (Arce et al., 2007; Hsieh et al.,
2023). Within the outflow cavity it is possible that cosmic rays are excluded.
However, energetic particles might still penetrate through envelope and into
the protostellar disk through equatorial zones. Overall, it is possible that
changes in the ionization rate are associated with the clearing of the
envelope from the protostellar to the protoplanetary disk stage. It has long
been recognized that cosmic ray ionization is central towards the isotopic
fractionation of carbon in the interstellar medium (Langer et al., 1984;
Langer & Graedel, 1989; Furuya et al., 2011; Roueff et al., 2015). Perhaps
ionization evolution might also provide some imprint on isotopic ratios
carried through subsequent evolutionary stages.

####  4.2.2 Model and Results

To model an earlier stage TW Hya protostellar  disk (e.g., Class 0/I), we
assume the grains have grown but only 90% of their mass is in large grains
(i.e., mm-sized). We also include the presence of cosmic rays ($\zeta_{\rm
H_{2}}=5\times 10^{-17}$ s-1). This value is at a level a factor of 50 above
that currently estimated towards TW Hya by Cleeves et al. (2015) which was
used in the late-stage disk simulation. The initial chemical abundances, given
in Table 2, place the initial CO abundance at 2.7 $\times$ 10-4 (i.e.
interstellar) with no additional carbon in any form. Thus, we assume that
there has been no chemical processing of CO as assumed in the earlier model.
The initial CO isotopic ratio is also set to an interstellar isotopic ratio of
77. At the start of the simulation the only reservoir of 12C and 13C is CO.
This model has some water ice included with an abundance 2.7 $\times$ 10-6 and
is run for 1 Myr. Here we assume that a significant amount of water ice has
accumulated in the midplane where most of the dust mass still resides (Krijt
et al., 2016). The baseline model structure parameters are given in the
appendix.

Figure 6: Predicted vertical distribution of CO and CCH isotopic ratios in the
“early-stage disk” model at a radial location commensurate with the center of
the C2H emission ring (63 au). (a) Overall gas density structure and modeled
vertical gas thermal profile at this radius. (b), (c), (d) Predicted vertical
distribution of 12CO/13CO, 12C2H/13C12CH, and 12CH4/13CH4 from the
thermochemical model. In these three panels the chemical abundances are
referenced to the left axis while the dashed line provides the predicted
isotopic ratio that is referenced to the axis on the right. In this model all
13C is placed in CO at the start with an initial isotopic ratio of 77, which
is shown as the dotted grey line. This model is sampled at r=63 au, the
location of the C2H emission ring.

Fig. 6 presents the results from this model at r = 63 au in the form of
vertical cuts of the abundance distribution of key isotope carriers. Fig. 7
provides the overall carbon inventory over this same region. Comparison of the
thermal structure between the late-stage (Fig. 4) and early-stage (Fig. 6)
disk model show that the less evolved system, with a larger mass in small
grains, traps radiation higher in the disk leading towards slightly cooler
deeper layers. Further the gas phase abundance of CO is distributed within the
warm molecular layer (i.e., above the CO freeze-out layer in the midplane and
below the photodissociation layer; Aikawa et al., 2002) which exists over a
wider range in vertical distance (Fig. 6b). As in the late-stage disk model
(Fig. 4), C2H forms at intermediate heights albeit in slightly higher layers
in the early-stage disk model shown in panel (c) in Fig. 6. One additional
difference for C2H is the peak abundance of C2H is reduced by about an order
of magnitude as the amount of free carbon is reduced in this model. The total
C2H column density from this model (at r = 63 au) is 5 $\times$ 1012 cm-2.
This is two orders of magnitude reduced from the late-stage model. In this
model the C/O ratio in the gas is near unity and a reduction in the C2H
column/line flux density is expected.

Figure 7: Plot of (top) major carriers of 12C and (bottom) 13C. Abundances are
normalized to the total amount of isotopic carbon with $x$(12C) = $2.7\times
10^{-4}$ and $x$(13C) = $3.5\times 10^{-6}$ (abundances given relative to
total H). This model is sampled at r=63 au, the location of the C2H emission
ring.

In the early-stage disk model at high altitude (z=35 au; Fig. 6b) CO is
destroyed via photodissociation. Within the layer where CO is photodestroyed
the effects of isotopic selective photodissociation of CO are observed and
12CO/13CO exceeds the assumed initial value shown as the dotted line. More
substantive changes are seen deep within the warm molecular layer where cosmic
rays chemically process about 40% of gas phase CO into methane ice (Fig. 6b
and c). The chemical processing of CO is initiated via reactions with He+
atoms which free 12C+ and 13C+ alongside O atoms (Aikawa et al., 1996). The
oxygen atoms collide with cold dust grains and hydrogenate to form water ice.
In this model about 40% of the O in CO is placed into water ice leaving excess
C+ in the gas. This causes the drop in the CO abundance near z = 14 au Fig.
6b).

This excess ionized carbon would normally find its way to CO, but the lack of
oxygen means it sits in a sea of H2 gas. As noted by Du et al. (2015) this
leads to the formation of unsaturated and saturated hydrocarbons that work
towards species that can freeze onto grains when their individual sublimation
temperature is above the temperature of the dust. In this instance the gas-
phase chemistry favors the formation of mostly CH4 ice but also C2H2 and HC2H2
ices.

In the layer where gas phase CO chemical processing occurs, cosmic rays are
also powering carbon isotopic fractionation. In this case, the He+ destruction
of CO and 13CO produces an excess 13C+ atoms which are favored to be implanted
into 13CO via reaction (1) due to the low gas temperature in these layers
(T${}_{gas}<$ 40 K). In all, 13CO gas becomes modestly enriched with 12C/13C
$\sim$ 50 compared to the initial value of 77. Because the 13C becomes locked
in gaseous CO, the hydrocarbon ices carriers have 12C/13C $>$ (12C/13C)ISM.

Figure 8: (top) predicted 12CO/13CO in CO and CH4 ice as a function of radius
from the early-stage disk model. Isotopic ratios are derived through the ratio
of the total column densities of the respective isotopologue. (bottom)
Predicted total column density profiles for gas phase 12CO and solid state CH4
as a function of radius for both molecules. Inside $\sim 12$ au the CH4
iceline is reached in the midplane and other more complex hydrocarbons with
assumed higher sublimation temperatures (e.g. HCCH2) are created. The drop in
the gas phase CO abundance near 4 au is due to an assumed inner dust cavity.

Similar behavior is seen throughout the disk, but at variable height about the
midplane. In Fig. 8 we provide the 12C/13C ratio in CO and methane ice as a
function of radius. For this plot we have computed the ratio from the
estimated vertical column densities. This demonstrates that this model
effectively produces two isotopic reservoirs for carbon across much of the
disk. In the bottom panel of the Figure, we also provide the column densities
of the main isotopologues. The comparable column density between gaseous CO
and CH4 ice shows that this effect is present both radially and vertically. We
note that the midplane (z $<$ 3 au) remains unchanged by this evolution with
CO ice as the dominant carbon carrier for both isotopologues and the ratio
fixed at 77.

One issue that we have not addressed is fact that the CO isotopic ratio
appears to be elevated above ISM values beyond 130 au (see Fig. 3). Earlier
work by Zhang et al. (2019) inferred extreme CO depletion values beyond 100 au
(i.e., abundances reduced by factors of 100 more compared to interstellar
levels). They suggested that this may be the result of gaseous dissipation of
the disk. If this is the case, then the elevated CO isotopic ratio at these
large distances could be the result of isotopic selective photodissociation.
If this is the case, then perhaps isotopic ratios might be used as a probe of
gas disk dissipation.

####  4.2.3 Critical Dependencies

We have made a number of assumptions in this model including the presence of
cosmic rays, increasing the amount of small grains, and reducing the amount of
water ice present in upper layers. Of these we find that the first two are
important for the generation of the isotopic enrichments. Cosmic ray
ionization has the highest penetration power, and if present, clearly will
dominate the ion-molecule chemistry leading to fractionation. However, high
levels of hard X-rays ($>$ 1 keV), with scattering, might also be important
(Igea & Glassgold, 1999). Increasing the ratio of small to large grains (i.e.
more mass in small grains) leads to a more flared surface and a higher mass
present below 35 K (e.g. compare Figs 4a and Fig. 6a). In this model this
effect is not large ($\sim$20% more mass below 35 K in the early-stage model),
but can be expected to be higher with less dust evolution and higher flaring
angles. Another aspect is the greater amount of small grains increases the
effective surface area ($<n\sigma>_{gr}$) of grains, enhancing the dust/gas
interaction and enabling greater freeze-out of methane. For the last factor,
the assumed water ice abundance, we find that if we assume an initial water
ice abundance in the early-stage model of 10-4 (relative to H) as opposed to
the original value of 2.6 $\times$ 10-6 then there is no major difference in
the overall chemistry compared to the results presented in Figs. 6 and 7.
While fractionation levels and abundances of key carriers (CO, CH4) are
unaffected, the C2H column is altered and is reduced by an order of magnitude.

##  5 Discussion

For the past decade observational analyses of resolved gas-rich disks have
revealed two puzzles. One regards the overall abundance of CO which appears to
be below the expected value of 10-4 relative to H2 (Favre et al., 2013;
Schwarz et al., 2016; Bergin & Williams, 2017; Miotello et al., 2017, 2023).
The other is the strong emission from C2H which rivals that of 13CO (Kastner
et al., 2014; Bergin et al., 2016; Miotello et al., 2019) and has been
attributed to a C/O ratio that exceeds unity. In the latter case, we cannot
rule out a source term for the excess CO in TW Hya from the photoablation of
refractory carbon grains as discussed by Bosman et al. (2021a). However, the
detection of enriched 13CO within an exoplanetary atmosphere (Zhang et al.,
2021) requires a source term for the enrichment. That source term is
apparently not carbon isotopic fractionation in the cold interstellar medium
based on initial results of ice abundances from JWST (McClure et al., 2023).

In this initial exploration, the evolution of an early-stage disk exposed to
cosmic ray ionization offers a potential avenue to generate a 13C enrichment
for CO. It simultaneously provides a source term for gas phase hydrocarbons in
evolved disk systems whose chemical formation pathways can potentially be
powered by photodesorption of hydrocarbon ices. Bosman et al. (2021a) shows
that if this carbon is released into the gas via photodesorption, it can
generate a long-lived photochemical equilibrium cycle.

###  5.1 Disk Ionization

The central need to generate the carbon isotopic fractionation points to the
presence of an ionization source deep in the disk. To reach significant mass
this source of ionization likely needs the penetrating power of interstellar
cosmic rays (Eistrup et al., 2017; Schwarz et al., 2018; Bosman et al., 2018).
Based on an analogy to the Sun, and detailed models of cosmic ray-stellar wind
interactions it is thought that the interstellar cosmic rays responsible for
ionization of H2 are not present within disk systems that have dissipated
their natal envelopes (e.g., class II systems; Cleeves et al., 2013; Fujii &
Kimura, 2022). Analysis of the observed molecular ion chemistry in class II
systems also appear to support reduced interstellar cosmic ray ionization
rates (Cleeves et al., 2015; Aikawa et al., 2021; Seifert et al., 2021).

However, the timescales of the CO abundance evolution in disk systems suggest
a shift between the protostellar phase, where the envelope is present (e.g.,
Class 0 and I), and the protoplanetary disk phase (Zhang et al., 2020; Bergner
et al., 2020). Thus, chemical processing might occur via ionization during
earlier evolutionary phases. Indeed, numerous combined observational efforts,
over much larger spatial scales, suggest that particle ionization rates might
be quite high ($>$10-15 s-1) during some portion of the protostellar phase
(Ceccarelli et al., 2014; Favre et al., 2018; Cabedo et al., 2023; Lattanzi et
al., 2023). This is supported by theoretical work that suggest cosmic rays
could be generated via shocks within the jets (Padovani et al., 2016; Gaches &
Offner, 2018). It is worth nothing that these observations are obtained on a
spatial scale much larger than the young disk. Evidence for enhanced
ionization rates is not universal as van’t Hoff et al. (2022) measures a
cosmic ray ionization rate below the ISM value in the L1527 protostellar disk.
However, if ionization rates are elevated, even for a short timescale, then
chemical processing of CO within infalling material or in the protostellar
disk might provide the needed ionization to generate multiple isotopic
reservoirs.

###  5.2 Caveats and the Link to Planets

From the perspective of planet formation and the detected 13C enrichment in an
exoplanet by Zhang et al. (2021), the combination of this work and Yoshida et
al. (2022a) characterizes the carbon isotopic ratios within the main carbon
reservoirs in gas with C/O $>$ 1\. In this gas the 13C carrier is gas-phase
CO. If the early-stage disk composition persists to later stages, as observed
in TW Hya today, then perhaps gas phase carbon provides the 13C enrichment and
not pebbles. One complication is that our early-stage disk model places the
13C deficit as ice coatings of small grains in the upper layers (z = 5-15 au).
The evolution of these grains somehow must be separated from CO via dust
evolution. Since CO carries the opposite isotopic signature if both reservoirs
were available, they would effectively cancel out the enrichment when supplied
to a planet. How this might occur is uncertain and more detailed models of
coupled gas/grain evolution are needed.

Our models also do not fully match observations. The level of 13CO enrichment
in our early-stage disk model (assuming no additional) evolution is only
12CO/13CO = 50 and not the observed level of 20 (Yoshida et al., 2022a).
Further, the level of 12C enrichment (12C/13C $\sim$ 250) in hydrocarbon ices,
if they are source terms for C2H in later stages, is a factor of 3 or more
above the carbon isotopic ratio inferred from C2H via the detection of 13CCH
presented here. Our models suggest that, in the context of an assumed low
cosmic ray ionization rate, little CO fractionation occurs in the dense disk.
Thus, over time the active chemistry in surface layers, perhaps with weak
turbulent mixing, would tend to lower the 12C/13C ratio in C2H towards the
value found for CO. Thus, time evolution could be an answer for C2H. For CO
the higher levels of 13C enrichment would require a different model that is
colder and/or with higher ionization rates than the early-stage model adopted
here. Regardless the detection of two isotopic reservoirs in carbon (Yoshida
et al., 2022a) and nitrogen (Hily-Blant et al., 2019) is providing information
on the overall evolution of disk gas.

###  5.3 Ionization Evolution: Model Predictions

This work highlights the importance of sources of ionization in the dense
portion of the disk. Our late-stage model, which matches other aspects of the
TW Hya disk (Trapman et al., 2017; Bosman & Banzatti, 2019), assumes that the
cosmic ray flux impinging on the disk is reduced based on detailed analyses of
molecular ion emission in the TW Hya disk (Cleeves et al., 2015). Under that
assumption this model cannot create two disparate carbon isotopic reservoirs.
In our work we highlight the potential for ionization evolution as a possible
solution to some aspects of this conundrum. This does lead to some testable
predictions. If the methane (or hydrocarbon-rich) ice created via CO chemical
processing is the source term for the excess carbon that powers C2H formation,
then that ice is predicted to have an extreme carbon isotopic ratio. Further,
protostellar disks would have higher levels of methane/hydrocarbon ice
(relative to CO ice) when compared to interstellar ices. If this 13C poor
material (with 12CH4,ice/13CH4,ice $\sim$150–300) is the source term for
today’s TW Hya isotopic reservoir then some mixing, as discussed above, is
required to lower the isotopic ratio.

Another possibility is that the destruction of refractory carbon powers the
formation of C2H as suggested by Bosman et al. (2021a). In the solar system,
refractory carbon carries little carbon isotopic enrichment and is consistent
with the interstellar medium (Nomura et al., 2023). This source term would be
consistent with our late-stage model for C2H, but this model cannot
fractionate CO. In this case detailed models will be needed to explore how the
isotopic signature of this carbon would mix with that of CO to match
observations.

##  6 Summary

We present the first detection of 13CCH in the TW Hya protoplanetary disk
system and a detailed analysis of carbon isotopic fractionation in disk
systems. Our primary results are as follows.

  1. 1.

We determine an isotopic ratio for C2H of CCH/13CCH $=65\pm 20$ in gas
sampling the same radii where Yoshida et al. (2022a) find CO/13CO $=21\pm 5$.
As C2H is posited as the primary tracer of excess carbon in disks with C/O
$\geq 1$, this confirms a carbon isotopic dichotomy in this system as
suggested by Yoshida et al. (2022a). This is comparable to a similar result
found for nitrogen by Hily-Blant et al. (2019).

  2. 2.

We explore the origins of the carbon isotopic dichotomy through the use an
extant thermochemical model of the disk (Trapman et al., 2017; Bosman &
Banzatti, 2019) with a full model of carbon isotope chemistry (Miotello et
al., 2014) in a disk with C/O $>$ 1 with a reduced carbon content as
consistent with observational constraints (Bergin et al., 2016; Kama et al.,
2016a; Bosman & Banzatti, 2019; Yoshida et al., 2022b). We find that this
model is incapable of reproducing two independent reservoirs primarily due to
the lack of cosmic ray ionization in deep disk layers where the temperature
would be low enough to power the primary gas phase carbon isotopic
fractionation reaction. Additional efforts are needed to model the system
adopting all available constraints on both the gas temperature (Calahan et
al., 2021) and ionization (Cleeves et al., 2015) simultaneously with the
isotopic fractionation.

  3. 3.

We also explore a model of the TW Hya disk in an earlier evolutionary state
where the initial C/O ratio is unity, interstellar cosmic rays are present,
and there is no carbon depletion. We find that this model can create two
independent isotopic reservoirs one 12C-rich carried by methane and other
hydrocarbon ices and the other 13C-rich carried by gas phase carbon monoxide.
This dichotomy is predicted to be present throughout the disk and may lay the
seeds for the future chemical evolution of the disk as CH4 is a potential
source term for the strong C2H emission in disk systems (Bosman et al.,
2021a).

  4. 4.

We discuss the implications of this result for exoplanetary systems where gas
phase CO appears to be the carrier of 13C enrichments. Additional
observational constraints regarding the carbon isotopic ratio within both
Class I and Class II disks are sorely needed. Further, an understanding of the
ionization evolution of disk systems is critical.

We are grateful to the thorough referee who provided comments that improved
this contribution. E.A.B. acknowledges support from NSF grant No. 1907653 and
NASA NASA’s Emerging Worlds Program, grant 80NSSC20K0333, and Exoplanets
Research Program, grant 80NSSC20K0259. KW’s research was conducted at the Jet
Propulsion Laboratory, California Institute of Technology under contract with
the National Aeronautics and Space Administration. KW was supported by a grant
from NASA/Emerging Worlds program (18-EW182-0083). L.I.C. acknowledges support
from NASA ATP 80NSSC20K0529, the David and Lucille Packard Foundation, and the
Research Corporation for Scientific Advancement Cottrell Scholar Award. This
paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00198.S,
ADS/JAO.ALMA#2015.1.00308.S. 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.

## References

  * Aikawa et al. (1996) Aikawa, Y., Miyama, S. M., Nakano, T., & Umebayashi, T. 1996, ApJ, 467, 684
  * Aikawa et al. (2002) Aikawa, Y., van Zadelhoff, G. J., van Dishoeck, E. F., & Herbst, E. 2002, A&A, 386, 622
  * Aikawa et al. (2021) Aikawa, Y., Cataldi, G., Yamato, Y., et al. 2021, ApJS, 257, 13
  * Andrews et al. (2016) Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, in press.
  * Arce et al. (2007) Arce, H. G., Shepherd, D., Gueth, F., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 245
  * Barman et al. (2015) Barman, T. S., Konopacky, Q. M., Macintosh, B., & Marois, C. 2015, ApJ, 804, 61
  * Bergin et al. (2014) Bergin, E. A., Cleeves, L. I., Crockett, N., & Blake, G. A. 2014, Faraday Discuss., 168, 61.  http://dx.doi.org/10.1039/C4FD00003J
  * Bergin et al. (2016) Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101
  * Bergin & Williams (2017) Bergin, E. A., & Williams, J. P. 2017, Formation, Evolution, and Dynamics of Young Solar Systems, 445, in press
  * Bergner et al. (2020) Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2020, ApJ, 898, 97
  * Bosman et al. (2021a) Bosman, A. D., Alarcón, F., Zhang, K., & Bergin, E. A. 2021a, ApJ, 910, 3
  * Bosman & Banzatti (2019) Bosman, A. D., & Banzatti, A. 2019, A&A, 632, L10
  * Bosman et al. (2018) Bosman, A. D., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 618, A182
  * Bosman et al. (2021b) Bosman, A. D., Bergin, E. A., Loomis, R. A., et al. 2021b, ApJS, 257, 15
  * Brogi & Birkby (2021) Brogi, M., & Birkby, J. 2021, in ExoFrontiers; Big Questions in Exoplanetary Science, ed. N. Madhusudhan, 8–1
  * Bruderer (2013) Bruderer, S. 2013, A&A, 559, A46
  * Bruderer et al. (2012) Bruderer, S., van Dishoeck, E. F., Doty, S. D., & Herczeg, G. J. 2012, A&A, 541, A91
  * Cabedo et al. (2023) Cabedo, V., Maury, A., Girart, J. M., et al. 2023, A&A, 669, A90
  * Calahan et al. (2021) Calahan, J. K., Bergin, E., Zhang, K., et al. 2021, ApJ, 908, 8
  * Ceccarelli et al. (2014) Ceccarelli, C., Caselli, P., Bockelée-Morvan, D., et al. 2014, Protostars and Planets VI, 859
  * Cleeves et al. (2013) Cleeves, L. I., Adams, F. C., & Bergin, E. A. 2013, ApJ, 772, 5
  * Cleeves et al. (2015) Cleeves, L. I., Bergin, E. A., Qi, C., Adams, F. C., & Öberg, K. I. 2015, ApJ, 799, 204
  * Cleeves et al. (2018) Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2018, ApJ, 865, 155
  * Cridland et al. (2016) Cridland, A. J., Pudritz, R. E., & Birnstiel, T. 2016, MNRAS, submitted
  * Czekala et al. (2021) Czekala, I., Loomis, R. A., Teague, R., et al. 2021, ApJS, 257, 2
  * Du et al. (2015) Du, F., Bergin, E. A., & Hogerheijde, M. R. 2015, ApJ, 807, L32
  * Du et al. (2017) Du, F., Bergin, E. A., Hogerheijde, M., et al. 2017, ApJ, 842, 98
  * Eistrup et al. (2017) Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2017, A&A, in press
  * Favre et al. (2013) Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJ, 776, L38
  * Favre et al. (2018) Favre, C., Ceccarelli, C., López-Sepulcre, A., et al. 2018, ApJ, 859, 136
  * Fujii & Kimura (2022) Fujii, Y. I., & Kimura, S. S. 2022, ApJ, 937, L37
  * Furuya & Aikawa (2014) Furuya, K., & Aikawa, Y. 2014, ApJ, 790, 97
  * Furuya et al. (2011) Furuya, K., Aikawa, Y., Sakai, N., & Yamamoto, S. 2011, ApJ, 731, 38
  * Gaches & Offner (2018) Gaches, B. A. L., & Offner, S. S. R. 2018, ApJ, 861, 87
  * Guillot et al. (2023) Guillot, T., Fletcher, L. N., Helled, R., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 947
  * Helled & Guillot (2017) Helled, R., & Guillot, T. 2017, in Handbook of Exoplanets
  * Hily-Blant et al. (2019) Hily-Blant, P., Magalhaes de Souza, V., Kastner, J., & Forveille, T. 2019, A&A, 632, L12
  * Hsieh et al. (2023) Hsieh, C.-H., Arce, H. G., Li, Z.-Y., et al. 2023, ApJ, 947, 25
  * Huang et al. (2018) Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42
  * Igea & Glassgold (1999) Igea, J., & Glassgold, A. E. 1999, ApJ, 518, 848
  * JWST Transiting Exoplanet Community Early Release Science Team et al. (2023) JWST Transiting Exoplanet Community Early Release Science Team, Ahrer, E.-M., Alderson, L., et al. 2023, Nature, 614, 649
  * Kama et al. (2016a) Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016a, A&A, in press.
  * Kama et al. (2016b) Kama, M., Bruderer, S., Carney, M., et al. 2016b, A&A, in press
  * Kastner et al. (2014) Kastner, J. H., Hily-Blant, P., Rodriguez, D. R., Punzi, K., & Forveille, T. 2014, ApJ, 793, 55
  * Krijt et al. (2020) Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134
  * Krijt et al. (2016) Krijt, S., Ciesla, F. J., & Bergin, E. A. 2016, ApJ, in press
  * Langer & Graedel (1989) Langer, W. D., & Graedel, T. E. 1989, ApJS, 69, 241
  * Langer et al. (1984) Langer, W. D., Graedel, T. E., Frerking, M. A., & Armentrout, P. B. 1984, ApJ, 277, 581
  * Langer & Penzias (1993) Langer, W. D., & Penzias, A. A. 1993, ApJ, 408, 539
  * Lattanzi et al. (2023) Lattanzi, V., Alves, F. O., Padovani, M., et al. 2023, A&A, 671, A35
  * Lavie et al. (2017) Lavie, B., Mendonça, J. M., Mordasini, C., et al. 2017, AJ, 154, 91
  * Line et al. (2021) Line, M. R., Brogi, M., Bean, J. L., et al. 2021, Nature, 598, 580
  * Maret & Bergin (2007) Maret, S., & Bergin, E. A. 2007, ApJ, 664, 956
  * McClure et al. (2023) McClure, M. K., Rocha, W. R. M., Pontoppidan, K. M., et al. 2023, Nature Astronomy, 7, 431
  * Milam et al. (2005) Milam, S. N., Savage, C., Brewster, M. A., Ziurys, L. M., & Wyckoff, S. 2005, ApJ, 634, 1126
  * Minissale et al. (2022) Minissale, M., Aikawa, Y., Bergin, E., et al. 2022, ACS Earth and Space Chemistry, 6, 597
  * Miotello et al. (2014) Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96
  * Miotello et al. (2023) Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. C., & Kataoka, A. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 501
  * Miotello et al. (2017) Miotello, A., van Dishoeck, E. F., Williams, J. P., et al. 2017, A&A, 599, A113
  * Miotello et al. (2019) Miotello, A., Facchini, S., van Dishoeck, E. F., et al. 2019, A&A, 631, A69
  * Nomura et al. (2023) Nomura, H., Furuya, K., Cordiner, M. A., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 1075
  * Öberg et al. (2011a) Öberg, K. I., Boogert, A. C. A., Pontoppidan, K. M., et al. 2011a, ApJ, 740, 109
  * Öberg et al. (2011b) Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011b, ApJ, 743, L16
  * Öberg et al. (2021) Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1
  * Oreshenko et al. (2017) Oreshenko, M., Lavie, B., Grimm, S. L., et al. 2017, ApJ, 847, L3
  * Padovani et al. (2016) Padovani, M., Marcowith, A., Hennebelle, P., & Ferrière, K. 2016, A&A, 590, A8
  * Reboussin et al. (2015) Reboussin, L., Wakelam, V., Guilloteau, S., Hersant, F., & Dutrey, A. 2015, A&A, 579, A82
  * Redaelli et al. (2021) Redaelli, E., Sipilä, O., Padovani, M., et al. 2021, A&A, 656, A109
  * Röllig & Ossenkopf (2013) Röllig, M., & Ossenkopf, V. 2013, A&A, 550, A56
  * Roueff et al. (2015) Roueff, E., Loison, J. C., & Hickson, K. M. 2015, A&A, 576, A99
  * Schwarz et al. (2016) Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 823, 91
  * Schwarz et al. (2018) —. 2018, ApJ, 856, 85
  * Seifert et al. (2021) Seifert, R. A., Cleeves, L. I., Adams, F. C., & Li, Z.-Y. 2021, ApJ, 912, 136
  * Sokal et al. (2018) Sokal, K. R., Deen, C. P., Mace, G. N., et al. 2018, ApJ, 853, 120
  * Stelzer & Schmitt (2004) Stelzer, B., & Schmitt, J. H. M. M. 2004, A&A, 418, 687
  * Teague (2019) Teague, R. 2019, The Journal of Open Source Software, 4, 1632
  * Teague et al. (2019) Teague, R., Bae, J., & Bergin, E. A. 2019, Nature, 574, 378–381
  * Trapman et al. (2017) Trapman, L., Miotello, A., Kama, M., van Dishoeck, E. F., & Bruderer, S. 2017, A&A, 605, A69
  * van Boekel et al. (2017) van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132
  * van’t Hoff et al. (2022) van’t Hoff, M. L. R., Leemker, M., Tobin, J. J., et al. 2022, ApJ, 932, 6
  * Visser et al. (2009) Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323
  * Woods & Willacy (2009) Woods, P. M., & Willacy, K. 2009, ApJ, 693, 1360
  * Yen et al. (2016) Yen, H.-W., Koch, P. M., Liu, H. B., et al. 2016, ApJ, 832, 204
  * Yoshida et al. (2022a) Yoshida, T. C., Nomura, H., Furuya, K., Tsukagoshi, T., & Lee, S. 2022a, ApJ, 932, 126
  * Yoshida et al. (2022b) Yoshida, T. C., Nomura, H., Tsukagoshi, T., Furuya, K., & Ueda, T. 2022b, ApJ, 937, L14
  * Yoshida et al. (2024) Yoshida, T. C., Nomura, H., Furuya, K., et al. 2024, arXiv e-prints, arXiv:2403.00626
  * Zhang et al. (2017) Zhang, K., Bergin, E. A., Blake, G. A., Cleeves, L. I., & Schwarz, K. R. 2017, Nature Astronomy, 1, 0130
  * Zhang et al. (2019) Zhang, K., Bergin, E. A., Schwarz, K., Krijt, S., & Ciesla, F. 2019, ApJ, 883, 98
  * Zhang et al. (2020) Zhang, K., Schwarz, K. R., & Bergin, E. A. 2020, ApJ, 891, L17
  * Zhang et al. (2021) Zhang, Y., Snellen, I. A. G., Bohn, A. J., et al. 2021, Nature, 595, 370

##  Appendix A DALI model for TW Hya

This model is based on the code and parameters as published by Trapman et al.
(2017) and modified by Bosman & Banzatti (2019). We adopt this model as a
baseline and change a few parameters to investigate the chemical evolution at
a slightly earlier evolutionary stage.

Table 3: TW Hya Model Parametersa Stellar parameters
---
Stellar mass | 0.74 $M_{\odot}$
Stellar luminosity | 1 $L_{\odot}$
X-ray luminosity |  $10^{30}$ erg s-1
UV Luminosity | 2.7 $\times 10^{31}$ erg s-1
Disk parameters
Disk mass | 0.025 $M_{\odot}$
Critical radius ($R_{c}$) | 35 AU
Surface density slope ($\gamma$) | 1
scale-height at $R_{c}$ ($h_{c}$) | 0.1 rad
Flaring angle ($\psi$) | 0.3
Inner radius | 0.05 AU
Gap inner radius | 0.3 AU
Gap outer radius | 2.4 AU
Inner disk gas-to-dust | 100
Inner disk gas depletion | $10^{-2}$
Small dust size | 0.005–1 $\mu$m
Small dust fraction | 0.01 (late-stage), 0.10 (early-stage)
Large dust size | 0.005–1000 $\mu$m
Large dust fraction | 0.99 (late-stage), 0.90 (early-stage)
Large dust settling factor | 0.2
Cosmic Ray Ionization Rate | 10-19 s-1 (late-stage), 5 $\times$ 10-17 s-1 (early-stage)
a Trapman et al. (2017); Bosman & Banzatti (2019)