FAUST XXXII. Chemically layered structures and ionisation in the disk and envelope of NGC 1333 IRAS 4C

Corresponding author: Shota Notsu

[email protected] 28 Department of Chemistry, University of Virginia, McCormick Road, PO Box 400319, Charlottesville, VA 22904, USA

29 Department of Applied Physics, Faculty of Engineering, Kanagawa University, 3-27-1 Rokukakubashi, Kanagawa-ku, Yokohama, Kanagawa 221-8686, Japan

30 Mizusawa VLBI Observatory, National Astronomical Observatory of Japan (NAOJ), 2-12 Hoshiga-oka, Mizusawa, Oshu-shi, Iwate 023-0861, Japan

31 The Graduate University for Advanced Studies, SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

32 Department of Physics and Astronomy, Rice University, 6100 Main Street, MS-108, Houston, TX 77005, USA

33 Centro de Astrobiología (CSIC-INTA), Ctra. de Torrejón a Ajalvir, km 4, 28850, Torrejón de Ardoz, Spain

34 NRC Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, BC, V9E 2E7, Canada

35 Department of Physics and Astronomy, University of Victoria, Victoria, BC, V8P 5C2, Canada

36 Laboratoire d’Astrophysique de Bordeaux, B18N allée Geoffroy St-Hilaire, 33615 Pessac, France

37 Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, 58341 Morelia, Mexico

38 Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA

39 Department of Physics, National Sun Yat-Sen University, No. 70, Lien-Hai Road, Kaohsiung City 80424, Taiwan, R.O.C.

40 Center of Astronomy and Gravitation, National Taiwan Normal University, Taipei 116, Taiwan, R.O.C.

41 Università degli Studi di Florence, Dipartimento di Fisica e Astronomia, via G. Sansone 1, 50019 Sesto Fiorentino, Italy

42 National Radio Astronomy Observatory, PO Box O, Socorro, NM 87801, USA

43 Institute for Advanced Study, Tsinghua University, Beijing 100084, People’s Republic of China

44 Astronomical Institute, Tohoku University, 6-3 Aramaki, Aoba-ku, Sendai 980-8578, Japan

45 Department of Physics, Tamkang University, S320, No.151, Yingzhuan Rd., Tamsui Dist., New Taipei City 251301, Taiwan, R.O.C.

46 Instituto de Astronomía, Universidad Nacional Autónoma de México, AP106, Ensenada CP 22830, B.C., México

47 Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria, 16, I-20133 Milano, Italy

48 NASA Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

49 National Astronomical Observatory of Japan (NAOJ), 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

50 Institute of Low Temperature Science, Hokkaido University, N19W8, Kita-ku, Sapporo, Hokkaido 060-0819, Japan

51 Department of Physics and Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

52 Yukawa Institute for Theoretical Physics, Kyoto University, Oiwake-cho, Kitashirakawa, Sakyo-ku, Kyoto-shi, Kyoto-fu 606-8502, Japan

53 Departament de Química, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain

54 Graduate School of Informatics and Engineering, The University of Electro-Communications, Chofu, Tokyo 182-8585, Japan

55 Department of Astronomy, University of Illinois, Urbana, IL, USA

56 Steward Observatory, 933 N Cherry Ave., Tucson, AZ 85721, USA

57 Materials Science and Engineering, College of Engineering, Shibaura Institute of Technology, 3-7-5 Toyosu, Koto-ku, Tokyo 135-8548, Japan

58 The Graduate University for Advanced Studies (SOKENDAI), Shonan Village, Hayama, Kanagawa 240-0193, Japan

59 Department of Physics, Graduate School of Science, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

60 Research Center for the Early Universe, Graduate School of Science, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

Introduction

The physical and chemical conditions within protoplanetary disks determine the properties of forming planets, including mass and chemical composition (e.g., Endres et al. [49]; Eistrup et al. [48]; Cridland et al. [36]; Öberg & Bergin [100]; Öberg et al. [101, 102]; Notsu et al. [95, 98]; Mollière et al. [92]). In order to understand the origin of the chemical composition in both the solar system and exoplanets, it is necessary to investigate disk chemical structure at high spatial resolution around forming stars, from Class 0/I to Class II sources.

Along the vertical direction in the disks, gases are in hydrostatic equilibrium and the density decreases towards the disk surface. Based on thermo-chemical models, protoplanetary disks can be simplified into three chemical layers: the hot disk surface, the inner warm molecular layer, and the central cold freeze-out layer (e.g., Aikawa et al. [3, 2]; Woitke et al. [185]; Dutrey et al. [46]; Walsh et al. [180, 182]; Notsu et al. [96, 97]; Öberg et al. [102], see also § 4.1). The most direct way to investigate the distribution and origin of molecules is to observe these vertical structure of disks so as to resolve the distinct chemical layers.

Recently, spatially resolved observations of some Class I-II edge-on disks (i.e. disks whose rotational axis is perpendicular to the line of sight) have been conducted with molecular lines of simple molecules such as CO, CS, SO, $\mathrm{H}_{2}\mathrm{CO}$, and c-$\mathrm{C}_{3}\mathrm{H}_{2}$ (e.g., Dutrey et al. [45, 44]; Louvet et al. [84]; Podio et al. [120]; Teague et al. [155]; van’t Hoff et al. 2020b, 2023; Lin et al. [78]; Foucher et al. [56]). Podio et al. [120] reported a resolved disk vertical structure with multiple molecular tracers, in the edge-on disk of the Class I protostar IRAS 04302+2247, using Atacama Large Millimeter/Submillimeter Array (ALMA) observations from the ALMA-DOT program (ALMA chemical survey of Disk-Outflow sources in Taurus, see e.g., Garufi et al. [60]). They showed that molecular emission is vertically stratified, with ${}^{12}\mathrm{CO}$ observed at greater disk heights compared to both CS and $\mathrm{H}_{2}\mathrm{CO}$, which are nearly co-spatial.

In the outer midplane, the line emission decreases due to molecular freeze-out onto dust grains (freeze-out layer). van’t Hoff et al. (2020b) also showed the vertically stratified structures in this disk from their ALMA observations that $\mathrm{C}^{17}\mathrm{O}$ line emission originates in the disk midplane and $\mathrm{H}_{2}\mathrm{CO}$ line emission traces more surface layers. They reported the presence of CO freeze-out in the outermost region of the disk midplane ($\gtrsim100\ \mathrm{au}$). Lin et al. [78] also confirmed these chemical layered structures in the same disk using ALMA eDisk observations, and noted that the dust has yet to settle significantly in this disk. They reported a dust scale height of $\sim6$ au at a disk radius of 100 au, which is comparable to the gas scale height at the same radius. Villanave et al. (2023) also reported a similar dust scale height based on the analyses of ALMA and VLA${}^{1}$ observations. Villanave et al. (2024) showed from their JWST${}^{2}$ imaging observations from 2 to 21 $\mu$m that the dark lane in dust continuum emission does not vary significantly with wavelength, which are interpreted as an indication for intermediate-sized ($\sim10\ \mu$m) grains in the upper layers of the disk.

van’t Hoff et al. (2023) conducted high-spatial resolution ($\sim24$ au) observations of molecular line emission toward the Class 0 protostar L1527 IRS (also known as IRAS 04368+2557; hereafter L1527) as part of the ALMA Large Program “Early Planet Formation in Embedded Disks (eDisk)” (e.g., Ohashi et al. [103]). Previous ALMA observations revealed that CCH and c-$\mathrm{C}_{3}\mathrm{H}_{2}$ emission originates primarily from the inner envelope and disappears inside the disk (Sakai et al. [134, 133]). Consistent with this picture, van’t Hoff et al. (2023) showed that $\mathrm{C}^{18}\mathrm{O}$, $\mathrm{H}_{2}\mathrm{CO}$, and SO trace the disk and inner envelope (see also Zhang et al. [193] for SO), with $\mathrm{H}_{2}\mathrm{CO}$ emission arising predominantly in the disk surface layers. They also confirmed that c-$\mathrm{C}_{3}\mathrm{H}_{2}$ is absent in the disk and that its emission originates in the envelope (see also van’t Hoff et al. 2026). It remains an open question when this kind of chemical layered structure forms during the evolution from Class 0 to Class I/II sources.

In addition to the vertical layered structure, the radial chemical gradients driven by the radial variations in disk physical structures are also crucial for understanding the disk chemical evolution processes and the chemical origin of the composition in both the solar system and exoplanets. The disk thermal structure profiles determine the snowline and sootline³ positions, which influence the distribution of elemental abundance ratios in the gas phase (such as C/O and N/O ratios) and complex organic molecules (e.g., Eistrup et al. 2016, 2018; van’t Hoff et al. 2020a; Ohno & Ueda 2021; Notsu et al. 2022; Aikawa et al. 2024; Houge et al. 2025; Sellek & van Dishoeck 2025). Moreover, the disk ionisation degree and ionisation rate $\zeta_{\mathrm{CR}}$, which affect gas-phase and grain-surface reactions, also have radial (and vertical) variations in disks and protostellar envelopes (e.g., Ceccarelli et al. 2014; Cleeves et al. 2014, 2015; Aikawa et al. 2021; Seifert et al. 2021; Fujii & Kimura 2022; van’t Hoff et al. 2022; Cabedo et al. 2023; Long et al. 2024; Schwarz et al. 2026). Recently, Cabedo et al. [17] reported the resolved map of $\zeta_{\mathrm{CR}}$ for the solar-type Class 0 protostar B335 using ALMA. They found an enhanced ionisation that increases at small envelope radii peaking at the central protostellar position (up to $\zeta_{\mathrm{CR}} \sim10^{-14}\ \mathrm{s}^{-1}$). They posit that the extreme $\zeta_{\mathrm{CR}}$ can be attributed to the presence of cosmic-rays accelerated close to the protostar (e.g., [53, 163, 121]). We note that Ceccarelli et al. [24] reported the enhanced ionisation rate of $>10^{-14}\ \mathrm{s}^{-1}$ for the whole envelope and $\sim10^{-12}\ \mathrm{s}^{-1}$ for the innermost envelope around the cluster of a few intermediate and low-mass protostars in the source OMC-2 FIR 4 with the total mass of $30M_{\odot}$ (e.g., [146, 83, 158]). Bergin et al. [7] suggested that such extreme cosmic-ray ionisation may help to generate the carbon isotope fractionation between molecules in younger disks, such as ${}^{12}\mathrm{C}$-enriched methane and hydrocarbon ices (see also [135]), and ${}^{13}\mathrm{C}$-enriched CO ([190]) and complex organic molecules (e.g., [67, 187, 16, 66, 69]). Therefore, it is important to resolve the radial distribution of the ionisation degree and ionisation rate from the envelope to the innermost disk scale ($<100\ \mathrm{au}$) in Class 0/I objects.

In this paper, we report the results of molecular line observations of the disk and envelope (Keplerian disk + pseudo disk, see § 3.2 around NGC 1333 IRAS 4C (hereafter IRAS 4C)⁴ in the Perseus Molecular Cloud, as part of the ALMA Large Program Fifty AU STudy of the chemistry in the disk/envelope system of solar-like protostars (FAUST, see e.g., [33])⁵ The FAUST program explores physical and chemical structures on a 50 au scale toward 13 nearby young protostellar sources with the same surface brightness sensitivity (e.g., Bianchi et al. 2020; Codella et al. 2021, 2022, 2024; Okoda et al. 2021, 2023; Imai et al. 2022; Ohashi et al. 2022; Vastel et al. 2022, 2024; Cacciapuoti et al. 2023, 2025; Evans et al. 2023; Mercimek et al. 2023, 2025; Chahine et al. 2024a,b; De Simone et al. 2024; Maureira et al. 2024, 2026; Podio et al. 2024; Sabatini et al. 2024, 2025; Zhang et al. 2024; Durán et al. 2025; Frediani et al. 2025; Hirota et al. 2025; Martínez-Henares et al. 2025; Oya et al. 2025; Cuello et al. 2026; Saury et al. 2026; Yang et al. 2026). Previous observations of IRAS 4C (ALMA/PEACHES, [191, 189, 5]) have shown that the disk rotation axis of this source is largely tilted to the line of sight, making it suitable for studying the molecular abundance distribution along the disk height (see also § 2.1). In this paper, we focus on the vertical and radial chemical structures in the disk and envelope. The molecular lines analyzed in this paper are listed in Table 1. We analyzed molecular line emission of CCH, c-C$_3$H$_2$, H$_2$CO, and C$^{18}$O for vertical chemical layers and H$^{13}$CO$^+$, DCO$^+$, and C$^{18}$O for ionisation.

Molecule	Transition	Frequency (GHz)	$E_l$ (K)	$E_u$ (K)	$A_{ul}$ (s$^{-1}$)	$n_{\rm cr}$ (cm$^{-3}$)
C$^{18}$O	$J = 2 - 1$	219.5603541	5.269	15.806	$6.01 \times10^{-7}$	$1 \times10^{4}$
CCH	$N=3$-2, $J=7/2$-5/2, $F=4$-3	262.0042600	12.575	25.149	$5.31 \times10^{-5}$	$8 \times10^{6}$
	$N=3$-2, $J=7/2$-5/2, $F=3$-2	262.0064820	12.574	25.148	$5.10 \times10^{-5}$	$8 \times10^{6}$
c-C$_3$H$_2$	$6_{0,6} - 5_{1,5}$	217.8221480	28.154	38.607	$5.40 \times10^{-4}$	$5 \times10^{7}$
	$6_{1,6} - 5_{0,5}$	217.8221480	28.154	38.607	$5.40 \times10^{-4}$	$5 \times10^{7}$
H$_2$CO	$3_{0,3} - 2_{0,2}$	218.2221920	10.483	20.956	$2.82 \times10^{-4}$	$6 \times10^{5}$
H$^{13}$CO$^+$	$J = 3 - 2$	260.2553390	12.491	24.981	$1.34 \times10^{-3}$	$2 \times10^{6}$
DCO$^+$	$J = 3 - 2$	216.1125822	10.372	20.744	$2.38 \times10^{-3}$	$2 \times10^{6}$

The paper is organized as follows. § 2 introduces the previous studies on our target, and outlines the observational setup and our data reduction process. The results and discussion of our observations are described in § 3 and § 4, respectively. In § 5, we provide our conclusions.

Observations

Target Object

Our target IRAS 4C is a low-mass protostar ($L_{\mathrm{bol}} = 1.1L_{\odot}$; Enoch et al. 2009; Sadavoy et al. 2014; Tobin et al. 2016)⁶ located in the NGC 1333 star-forming region in the Perseus Molecular Cloud at a distance of 293 pc ([108])⁷. It is classified as a Class 0 source based on its bolometric temperature of $T_{\mathrm{bol}} = 31$ K (e.g., [50, 126, 72, 156]). IRAS 4C lies 47″ northeast of IRAS 4A, and 54″ north-northeast of IRAS 4B. IRAS 4C does not show hot corino activity although both IRAS 4A and 4B are hot corino sources which are extremely rich in complex organic molecules (e.g., [116, 152, 82, 131, 41, 25 b] for IRAS 4A; [136, 12, 116, 41] for IRAS 4B). Koumpia et al. [72] suggested that IRAS 4C seems to be in a different evolutionary state than IRAS 4A and 4B based on its lower luminosity and colder nature than the other two sources, but they also discussed that neither the younger nor the older age scenario seems to be fully consistent with the body of observational data on IRAS 4C. The inclination $i$ between the disk plane and the line of sight is about $25^\circ$ ([157, 139]). In this paper, we define the inclination such that $i = 0^\circ$ for an edge-on disk and $i = 90^\circ$ for a face-on disk. Note that this inclination value is estimated by the image-plane 2D Gaussian fitting of the dust continuum emission ([139]) assuming geometrically thin disk, and thus it should be considered as lower limit. Since IRAS 4C is not fully edge-on disk ($i = 25^\circ$), in integrated intensity maps, even at the disk surface there are vertical contamination from the midplane. If the disk radius is around 200 au (based on continuum in Figure 1 and Figure 2), emission from the midplane side would be mixed in by up to around $200\ \mathrm{au} \times\sin(25^\circ) \sim84\ \mathrm{au}$, which would contaminate to the emission in the disk surface. However, such contamination will concentrate around the region close to the protostar only (small $R$, see details in § 3.2).

Zhang et al. [191] observed this object in the ALMA/PEACHES survey and resolved the outflow cavity structure in the CCH and CS line emission with clear signatures of rotation with respect to the outflow axis. The velocity pattern is highly symmetric with respect to the outflow axis, and the northern cavity wall is blue-shifted while the southern cavity wall is red-shifted. While the CCH emission appears to mostly trace the outflow cavity structure, the CS emission traces the envelope. The rotation is detected from about 150 au up to about 1750 au above the envelope/disk midplane. As the distance from the central source increases, the rotation velocity of the outflow decreases while the outflow width increases, which gives a flat specific angular momentum distribution along the outflow. From the CS and SO emission which trace the envelope and disk in this source, Zhang et al. [191] roughly estimated the centrifugal barrier radius of $0.2^{\prime\prime} - 0.6^{\prime\prime}$ ($\sim60 - 180$ au). We note that this centrifugal barrier radius is the boundary between the infalling envelope and the rotating disk (e.g., [134, 112]).

At the base of the outflow cavity, previous observations have found a decrement in the molecular emission (such as CCH, see [191]) along the disk midplane. However, due to the lack of spatial resolution and the limited number of observed molecules, it was unclear if this effect is due to a dark lane caused by optically thick dust emission or chemical structures (e.g., [120, 78] (Li et al. 2017)).

Observational Setup and Data Reduction

The observations of NGC 1333 IRAS 4C were conducted with ALMA (Cycle 6) between 2018 November and 2021 June as part of the large program FAUST (2018.1.01205.L, PI: Satoshi Yamamoto). The molecular lines analyzed in this paper (see Table 1) are included in two frequency settings in Band 6; 216–235 GHz (Setup 1: 1.3 mm) and 244–262 GHz (Setup 2: 1.2 mm). We used the 12 m array data from two different configurations (C43-6 and C43-3 for extended and compact configurations, respectively) and the 7 m array data of the Atacama Compact Array (ACA/Morita Array, Setup 2 only), combining these visibility data in the UV plane. The ACA was employed to ensure a maximum recoverable size of $8^{\prime\prime}.4$ ($= 2460$ au, at the source distance). We note that such maximum recoverable size is reached by the 12 m array for setup 1, therefore no ACA observation was performed for setup 1. The field center is $(\alpha_{2000}, \delta_{2000}) = (03^{\mathrm h}29^{\mathrm m}13^{\mathrm s}.550, +31^\circ13^{\prime}58^{\prime\prime}.110)$, which is close to the protostellar position. The spectral windows optimized for targeting specific molecular line observations were set to a bandwidth of 58.59 MHz and a resolution of 141 kHz. The channel width corresponds to a velocity resolution of $\sim0.19\ \mathrm{km\ s^{-1}}$ and $\sim0.17\ \mathrm{km\ s^{-1}}$ for Setup 1 and Setup 2, respectively. The spectral windows used for broad continuum coverage in the observations have a bandwidth of 1875 MHz and a channel width of 0.488 MHz (Setup 1) and 0.977 MHz (Setup 2). This channel width corresponds to a velocity resolution of $\sim1.25$ km s$^{-1}$ and $\sim1.38$ km s$^{-1}$ for Setup 1 and Setup 2, respectively.

The Common Astronomy Software Applications (CASA) package ([89, 20]) was used for data reduction and imaging. The data were calibrated using a modified version of the ALMA calibration pipeline in CASA v5.6.1-8, with an additional in-house calibration routine to correct for the system temperature and spectral line data normalization.⁸ Self-calibration was carried out using line-free continuum emission for each configuration. Two rounds of phase-only self-calibration were performed using solution intervals of 6 seconds (for the more compact configurations) or 18 seconds (for the extended configuration). This was followed by amplitude self-calibration using per-scan solutions and solnorm=True in the CASA task gaincal. Finally, once the source model was as complete as possible, offsets between spectral windows were corrected using a single amplitude solution per observation using solnorm=False. The complex gain corrections derived from the self-calibration were then applied to all channels in the data, and the continuum model derived from the self-calibration was subtracted from the data to produce continuum-subtracted line data. A similar self-calibration technique was also used to align both amplitudes and phases across the multiple configurations. The overall improvement in the dynamic range of the continuum images is more than a factor of two. The uncertainty in the absolute flux density scale is approximately 10%.

Maps of the spectral line emission were obtained by CLEANing the dirty images with a Briggs robustness parameter of 0.5. The synthesized beams of the molecular line data are summarized in Table 2.

Molecule	Transition	$\sigma$ (mJy beam$^{-1}$)	Synthesized Beam
C$^{18}$O	$J = 2 - 1$	1.5	$0.''26 \times0.''19$ (P.A. = $5.^\circ5$)
CCH	$N=3$-2, $J=7/2$-5/2, $F=4$-3	1.5	$0.''26 \times0.''15$ (P.A. = $4.^\circ6$)
	$N=3$-2, $J=7/2$-5/2, $F=3$-2	1.3	$0.''26 \times0.''15$ (P.A. = $4.^\circ6$)
c-C$_3$H$_2$	$6_{0,6} - 5_{1,5}$	0.9	$0.''27 \times0.''20$ (P.A. = $-1.^\circ8$)
	$6_{1,6} - 5_{0,5}$	0.9	$0.''27 \times0.''20$ (P.A. = $-1.^\circ8$)
H$_2$CO	$3_{0,3} - 2_{0,2}$	0.7	$0.''27 \times0.''19$ (P.A. = $6.^\circ3$)
H$^{13}$CO$^+$	$J = 3 - 2$	1.3	$0.''26 \times0.''15$ (P.A. = $4.^\circ5$)
DCO$^+$	$J = 3 - 2$	0.8	$0.''27 \times0.''21$ (P.A. = $-1.^\circ4$)

Results

Dust Continuum Emission

Figure 1 shows the 1.3 mm continuum map, whose beam size is 0\farcs247 \times0\farcs178 (P.A. = 5\farcs61). We adopt the 1.3 mm continuum map in the present study as it has a higher angular resolution than the 1.2 mm continuum map. The rms noise level ($\sigma$) of this continuum map is 20 $\mu$Jy/beam. The peak position determined by a 2D Gaussian fit algorithm performed in the image and provided by CASA in the image is ($\alpha_{\mathrm{ICRS}}$, $\delta_{\mathrm{ICRS}}$) = (03^\mathrm{h}29^\mathrm{m}13\fs551, +31^\circ13\arcmin58\farcs08). The integrated flux density and peak intensity are 63.75 $\pm$ 0.79 mJy and 32.24 $\pm$ 0.27 mJy/beam, respectively. The peak intensity corresponds to a brightness temperature of 17.6 $\pm$ 0.1 K. The source image size de-convolved from beam is 0\farcs332 \pm0\farcs006 \times0\farcs110 \pm0\farcs003 with a P.A. of 14\fdg93 \pm0\fdg58. This corresponds to a linear size of $(97.3 \pm1.8)$ au $\times$ $(32.2 \pm0.9)$ au, which is consistent with previous estimates ([157]). The continuum emission is elongated in the north-south direction, which might trace the outflow launching region ([191]). This elongation of the aspect ratio of the continuum emission is consistent with the fact that the disk rotation axis of this source is largely tilted to the line of sight ([157]). In addition, there is faint ($\sim5 - 10\sigma$) continuum emission extended toward the outflow direction.

The Vertical Layered Structures

Figure 2a and (Figure 2b) show the integrated intensity (moment 0, between $4 - 12$ km s$^{-1}$) maps of CCH ($N=3-2$, $J=7/2-5/2$, $F=4-3$) and c-C$_3$H$_2$ ($6_{0,6} - 5_{1,5}$), respectively. The black contours show the 1.3 mm continuum emission intensity (Figure 1). Here, we define $R$ and $Z$ axis in the disk major axis direction (the disk/envelope direction) and minor axis direction (the outflow direction), respectively, where $Z$ is positive in NW direction and $R$ is positive in NE direction (see Figure 2a). We note that in this paper “disk” does not necessarily refer to the Keplerian disk only, and it cannot be distinguished from a flattened envelope (also known as a pseudo disk). Figure 2a and (Figure 2b) show the dark lane of apparently negative (CCH) and almost zero (c-C$_3$H$_2$) intensity, respectively, around the continuum emission peak (Z \lesssim0\farcs3). Emission from the molecules seems to trace the upper surface structures and outflow cavity walls.

We note that the CCH and and c-C$_3$H$_2$ emission extend toward the outflow direction (perpendicular to the disk/envelope direction) from a position slightly offset from the continuum peak, and the offsets are larger along the northern sides of the cavity ($\sim1.0 - 1.5''$) than on the southern sides ($\sim0.7''$). In addition, the length of the disk-like structure along the major axis of the continuum estimated from the CCH and c-C$_3$H$_2$ line emission is also longer toward the north ($\sim1.0 - 1.5''$) Note—$^{a}$ Line parameters are taken the from Cologne Database for Molecular Spectroscopy (Endres et al. [49]). $E_l$ and $E_u$ are the energies of the lower and upper levels, respectively. $A_{ul}$ is the Einstein A coefficient. $n_{\rm cr}$ is the critical density for $T_{\rm ex} = 50$ K (Chandra & Kegel [27]; Schöier et al. [142]; van der Tak et al. [175]).

Note—$^{b}$ $\sigma$ is the rms noise level. The beam size for each line is obtained from the observations.

than toward the south ($\sim0.7''$), although the continuum emission ($> 5\sigma$) does not show such clear asymmetry (see Figure 1). The distribution of the CCH line emission is consistent with the results of the previous ALMA observations. Zhang et al. [191] analyzed the velocity structure of the CCH line emission toward the outflow direction carefully. They described from the moment 1 map that toward the outflow direction, CCH emission has constant velocity at least $Z \lesssim0.5''$. We note that the apparent starting radii of CCH extended emission in the $Z$ direction ($R \sim0-0.''1$) is not the same as the outflow launching radius ($R \lesssim0.''05$) derived from angular momentum measurements by Zhang et al. [191]. In this paper we mainly focus on the disk/envelope component.

Figure 2c and (Figure 2d) show the integrated intensity maps of C$^{18}$O ($J = 2 - 1$) and H$_2$CO ($3_{0,3} - 2_{0,2}$), respectively. Unlike the CCH and c-C$_3$H$_2$ line emission, the C$^{18}$O line emission is concentrated around the central source and it exists along the midplane of the disk/envelope structure. The spatial extent is similar to the continuum emission ($R \sim300$ au), although the line intensity decreases around the continuum peak position ($R \lesssim50$ au). The H$_2$CO line emission appears to trace mainly the warm molecular layer ($\sim50-150$ au from the central source), i.e. between the line emitting regions of C$^{18}$O and CCH/c-C$_3$H$_2$. It is also located along the outflow cavity walls. The length of the H$_2$CO line emission along the disk-envelope direction is also longer toward the north ($\sim1.0-1.5''$) than toward the south ($\sim0.7''$), similar to CCH and c-C$_3$H$_2$ line emission lengths. Evans et al. [51] and Martínez-Henares et al. [86] discussed the similar asymmetric structures in H$_2$CO line emission of the Class 0/I protobinary system [BHB2007] 11.

We note that the H$_2$CO emitting layer has a slightly parabolic X-shape (especially in the side with $Z > 0$), and the $Z$-coordinate of the emitting layer increases as $R$-coordinate increases (see Figure 2d). This result might suggest that contamination from emission originating in the disk midplane is less significant than emission from the disk surface layer, since the radius $R of the contaminating emission is expected to be largest in the region at $Z \sim0$ (see also Section Figure 2.1). We note that the parabolic X-shape structure of H$_2$CO line emission was also confirmed by previous ALMA observations in edge-on disks, such as the disk around the Class I protostar IRAS 04302+2247 (e.g., Podio et al. [120]).

Figures 3 (a), (b), and (c) show the integrated intensity maps of c-C$_3$H$_2$, H$_2$CO, and C$^{18}$O overlaid with the CCH emission, which focus on the disk/envelope component. The line emitting regions of CCH and c-C$_3$H$_2$ are similar to each other. The left panel of Figure 4 shows the profiles of the integrated intensities of C$^{18}$O, H$_2$CO, CCH, and c-C$_3$H$_2$ lines in the vertical (disk minor axis, $Z$) direction passing through the protostar position. The distances of the line intensity peaks from the central source are $\sim50$ au for C$^{18}$O, $\sim50 - 150$ au for H$_2$CO, and $\sim150 - 300$ au for CCH and c-C$_3$H$_2$.

In Figure 2a, (Figure 2b), (Figure 2c), and (Figure 2d), while the “hole” seen in the C$^{18}$O emission map could originate from high dust optical depth, the dark lane found in CCH and c-C$_3$H$_2$ emission map is unlikely to be solely due to the high optical depth. We focus on the individual channel maps after the continuum subtraction, i.e., (Line brightness temperature) - (Dust brightness temperature). For the C$^{18}$O and H$_2$CO, the fact that there is a dip in the emission coincident with the dust continuum but it does not go negative except at the continuum peak means that the C$^{18}$O and H$_2$CO emission has a brightness temperature similar to, or greater than, that of the dust. It is therefore more natural to consider that this dark lane seen in CCH and c-C$_3$H$_2$ maps is due to the scarcity of the chemical species (See § 4.1 for details). Indeed, the dark lane continues to the outer part of the disk/envelope system, $\sim1''$ ($\sim300$ au) away from the protostar in the disk/envelope direction, where the dust continuum emission is very weak. We note that if the emission from CCH was obscured by the dust, we would not expect to have C$^{18}$O at the same region.

We note that the reason why both C$^{18}$O and H$_2$CO have very weak emission in the central region where we have the peak continuum emission could be (1) because the molecules freeze onto the dust grains or (2) because the continuum emission is optically thick thus absorbing the line emission. The second hypothesis is more likely because we would see no emission if the C$^{18}$O molecules freeze onto dust grains. In addition, the first hypothesis can be denied because the line emission have negative components in Position-Velocity diagrams (Figure 6).

The Radial Profiles of C$^{18}$O, H$^{13}$CO$^+$, and DCO$^+$ Line Emission

Figure 2e and (Figure 2f) show the integrated intensity maps of H$^{13}$CO$^+$ ($J = 3 - 2$), and DCO$^+$ ($J = 3 - 2$) lines overlaid on the dust continuum emission contours. Figure 3 (e) and (Figure 3f) show the integrated intensity maps of H$^{13}$CO$^+$ ($J = 3 - 2$), and DCO$^+$ ($J = 3 - 2$) overlaid with the CCH line emission. Figure 5 shows the integrated intensity maps of H$^{13}$CO$^+$ ($J = 3 - 2$), and DCO$^+$ ($J = 3 - 2$) overlaid with the C$^{18}$O line emission. The H$^{13}$CO$^+$ line emission mainly exists along the midplane of the disk/envelope structure. The intensity peaks of the H$^{13}$CO$^+$ line ($\sim0.5^{\prime\prime}$ from the central source) are located outside the peaks of the continuum emission and C$^{18}$O line emission.

The DCO$^+$ line emission mainly exists along the midplane of the disk/envelope structure, although there is weak emission from DCO$^+$ which is located at the outflow cavity walls. The DCO$^+$ line emission from northern side ($> 0.01$ Jy beam$^{-1}$ km s$^{-1}$) is stronger than that from southern side ($\sim0.008$ Jy beam$^{-1}$ km s$^{-1}$ at most). In addition, the intensity peak of the DCO$^+$ line ($> 1.0^{\prime\prime}$ from the central source) is located outside the C$^{18}$O and H$^{13}$CO$^+$ line emission.

The right panel of Figure 4 shows the radial profiles of the integrated intensities of C$^{18}$O, H$^{13}$CO$^+$, and DCO$^+$ lines in the radial (disk major axis, $R$) direction. The distances of the line intensity peaks from the central source are $\sim50 - 100$ au for C$^{18}$O, $\sim100 - 200$ au for H$^{13}$CO$^+$, and $\gtrsim300$ au for DCO$^+$.

Figure 6 shows the position-velocity (PV) diagrams of C$^{18}$O, H$^{13}$CO$^+$, H$_2$CO, and DCO$^+$ line emission along the cuts perpendicular to the outflow axis and passing through the continuum peak, i.e., along the major axis of the continuum emission and therefore along the disk ($R$ direction). We set the averaging widths as $0.1^{\prime\prime}$, so that we extract the emission component around the midplane only. The white and gray curves correspond to the rotation profile assuming Keplerian rotation. In calculating the rotation curves, we assume the distance of the source as 293 pc, the central star mass as 0.18 M$_\odot$ (based on previous estimates in Ortiz-León et al. [108]; Zhang et al. [191], white curves) and 0.06 M$_\odot$ (gray curves), the disk inclination angle $i$ as 25° (Tobin et al. [157]; Segura-Cox et al. [139]), and the source systemic velocity as 7.7 km s$^{-1}$. Here, we adopted inclination between the disk plane and the line of sight, thus we assumed completely edge-on as 0° instead of 90° throughout this paper. We note that the rotation curves with lower central star mass (0.06 M$_\odot$) well fit the PV diagrams of C$^{18}$O and H$^{13}$CO$^+$ than those with 0.18 M$_\odot$, suggesting that the central star mass might be smaller than the previous estimates. We note that fitting infalling envelope models would also be useful to investigate the velocity profiles in detail (e.g., Oya et al. [110]). The detailed kinematic analyses of the disk/envelope system will be described in a separate

FAUST publication, combined with the outflow emission analyses.

In Figure 6, the emission peaks of the C$^{18}$O line are located at $\lesssim0.5$″ and $\gtrsim1$ km s$^{-1}$ away from the source systemic velocity position ($\sim7.7$ km s$^{-1}$; Higuchi et al. [62]; Zhang et al. [191]) in second and fourth quadrant, whereas those of H$^{13}$CO$^{+}$ and DCO$^{+}$ are located $\sim1.0$″ and $< 1$ km s$^{-1}$ away. We note that the differences between H$^{13}$CO$^{+}$ and DCO$^{+}$ are also clear in the PV diagrams. There are second emission peaks at $\sim1.5$ km s$^{-1}$ observed within the inner regions in H$^{13}$CO$^{+}$, which are not seen in the DCO$^{+}$ emission.

Derivation of Column Densities and Optical Depth

Figure 7 shows the distribution of column densities for C$^{18}$O, H$^{13}$CO$^{+}$, DCO$^{+}$, and H$_2$. We calculated the column densities $N_i$ except for H$_2$ based on the following equation which Persson et al. [116] and Cabedo et al. [17] adopted:

where $\lambda$ and $\nu$ are the wavelength and frequency of the transition, $A_{ul}$ is the Einstein A coefficient, $k_{\mathrm{B}}$ is the Boltzmann constant, $g_u$ is the upper state degeneracy, $Q_{\mathrm{rot}}$ is the partition function at the temperature $T = T_{\mathrm{ex}}$, $E_l$ is the energy of the lower level, $B_{\nu}(T)$ is the Planck function for the background ($T_{\mathrm{bg}} = 2.7$ K) and for the excitation temperature of the line ($T_{\mathrm{ex}}$), and $\int I_{\nu}\mathrm{d}\nu$ is the integrated intensity. This equation assumes optically thin emission. First of all, in Figure 7 and Figure 8 we assume $T_{\mathrm{ex}} = 50$ K throughout the disk/envelope system, since we mainly focus on the inner region where the SO line emission is detected on the scale of 1″ (300 au) radius from the central protostar (Zhang et al. [191]). The sublimation temperature of SO in disk condition is around 50 K (e.g., Sakai et al. [134]). We have confirmed that the column density values increase / decrease only by a factor of around two by changing $T_{\mathrm{ex}}$ to 20 K /100 K, respectively.

The bolometric luminosity $L_{\mathrm{bol}}$ is 1.1$L_{\odot}$ for IRAS 4C (See § 2.1). Thus, assuming that the dust temperature in the midplane $T_d(R)$ is determined by the balance of heating from the central star radiation and cooling by radiation of optically thick dust, we can estimate $T_d(R)$ by the following equation (e.g., Ceccarelli

et al. 2000),

where we assume that $T_d(R) = 50$ K at $L_{\mathrm{bol}= 1.1L_\odot}$ and $R = 300$ au (on the scale of $1^{\prime\prime}$) based on the spatial extent of the detected SO line emission (Zhang et al. [191]). In this case, the disk radii where $T(R) < 20$ K and $T(R) > 100$ K are around $R \gtrsim1.9 \times10^3$ au and around $R \lesssim76$ au. In the innermost region the $T_d(R)$ would be much larger because of the mass accretion heating (like the disk around Class I protostar R CrA IRS7B-a, Takakuwa et al. [151]). We note that Tobin et al. [159] described the empirical equation of the average dust temperature for the whole disk as $T_d = 43(L_{\mathrm{bol}}/L_\odot)^{0.25}$. Based on this empirical equation, the averaged dust temperature $T_d$ is estimated to be 44 K.

The peak column densities are $\sim10^{13}\ \mathrm{cm}^{-2}$ for H$^{13}$CO$^+$, $\sim2\times10^{12}\ \mathrm{cm}^{-2}$ for DCO$^+$, and $\sim(3-4)\times10^{16}\ \mathrm{cm}^{-2}$ for C$^{18}$O. The peaks of the H$^{13}$CO$^+$ column density ($\sim0.5^{\prime\prime}$ ($\sim150$ au) from the central source) are located outside the peaks of the continuum emission and C$^{18}$O column density (see also Figure 5). The DCO$^+$ column density in the northern side ($\sim(1-2)\times10^{12}\ \mathrm{cm}^{-2}$) are larger than those on the southern side ($\lesssim5\times10^{11}\ \mathrm{cm}^{-2}$). In addition, the peaks of the DCO$^+$ column density ($>1.0^{\prime\prime}$ from the central source) are located outside the C$^{18}$O and H$^{13}$CO$^+$ column density peaks.

Here we note that the derivation of column densities for the molecules is dependent on the assumption that the emission is optically thin (van’t Hoff et al. 2018b; Sharma et al. 2025). Figure A1 in Appendix § 6.1 shows the line optical depth $\tau$ for H$^{13}$CO$^+$, DCO$^+$, and C$^{18}$O. We calculated the line optical depth $\tau$ based on the following equation which van’t Hoff et al. (2018a) adopted:

Image not extracted

where $N_i$ is the column densities calculated in equation (1) (see also Figure 7) and we assume $T_{\mathrm{ex}} = 50$ K and line velocity width $\Delta V = 0.17\ \mathrm{km\ s}^{-1}$ throughout the disk/envelope system. According to Figure A1 in Appendix § 6.1, the peak positions of optical depth are coincident with those of column densities, and the peak optical depth are $\sim0.7$ for DCO$^+$, $\sim1.5$ for H$^{13}$CO$^+$, and $\sim3.5$ for C$^{18}$O. The H$^{13}$CO$^+$ line emission is marginally optically thick ($\tau\sim1-1.5$) at $r \sim100-200$ au, and C$^{18}$O is optically thick $\tau\gtrsim2$ in $r \lesssim150$ au. Thus, the column densities of C$^{18}$O might be underestimated compared to those of H$^{13}$CO$^+$. We have confirmed that the values of $\tau$ increase / decrease only by a factor of around two by changing $T_{\mathrm{ex}}$ to 20 K /100 K, respectively.

The column densities of H$_2$ are calculated from the dust continuum emission using the following equation:

where $\mu_{\mathrm{H}_2} = 2.8$, $m_{\mathrm{H}}$ is the atomic hydrogen mass, $\Omega_{\mathrm{beam}}$ is the beam solid angle, $B_\nu(T_d)$ is the Planck function at the dust temperature $T_d$ (here we assume $T_{\mathrm{ex}} = T_d = 50$ K), $S_\nu^{\mathrm{beam}}$ is the dust continuum emission flux density per beam (see Figure 1), and $\kappa_\nu$ is the dust mass opacity per gas mass at the frequency $\nu$. For the dust opacity we adopt a power-law fit given by the following equation,

with a dust mass absorption coefficient $\kappa_0 = 1.6\ \mathrm{cm}^2\ \mathrm{g}^{-1}$ at $\lambda_0 = 1.3$ mm, following Ossenkopf & Henning [109] and Cabedo et al. [17]. We assume a standard gas-to-dust ratio $\xi_d = 100$, and a dust emissivity exponent $\beta= 1.7$ (= the value of interstellar average). Cacciapuoti et al. [19] reported the values of $\beta$ in protostellar envelopes for the whole samples of FAUST observations, and they found $\beta> 1.4$ for IRAS 4C. We note that some previous (sub)millimeter observations of protoplanetary disks have reported much smaller values of $\beta$ (such as $< 0$), which may be explained by dust self-scattering (Li et al. 2017). The resulting H$_2$ column densities are $3\times10^{23}\ \mathrm{cm}^{-2}$ within $0.5^{\prime\prime}$ from the central protostar.

Discussion

The Vertical Chemical Layered Structures

As we reported in § 3.2, it was found that the line emission of C$^{18}$O, H$_2$CO, c-C$_3$H$_2$, and CCH show a layered distribution in this order, outward and starting closest to the disk midplane. Two possibilities can be considered as the origin of this layered structure. One is the excitation effect of molecular line emission (that is, molecular lines with higher critical densities trace the regions closer to the disk midplane), and the other is the varying abundance of molecules in each layer, also referred to as chemical stratification.

The $\mathrm{C}^{18}\mathrm{O}$ line, which is detected most strongly in the region closest to the midplane and has the lowest critical density ($1\times10^{4}\,\mathrm{cm}^{-3}$ for $T_{\mathrm{ex}}=50$ K, Chandra & Kegel [27]; Schöier et al. [142]) among the observed species (see e.g., Schöier et al. [142]; Shirley [147]), suggests that $\mathrm{C}^{18}\mathrm{O}$ is abundant near the midplane ($Z\sim50$ au). When we compare the $\mathrm{C}^{18}\mathrm{O}$ map with the dust continuum map, the $\mathrm{C}^{18}\mathrm{O}$ emission appears slightly weaker along the midplane ($Z\sim0$), which is considered to be in the highest density condition (It is easily recognized by the contour map in Figure 2 and vertical intensity profiles in Figure 4). Since we observed only one transition, we cannot rule out the possibility that the intensity of the spectral line changes due to temperature variations. In addition, since the emission components of $\mathrm{C}^{18}\mathrm{O}$ line at around 8 km s$^{-1}$ on the Position-Velocity diagram (see Figure 6) are negative, we cannot remove the possibility that highly optically thick dust in the foreground obscured the information of the exact disk midplane. However, it is also natural to consider that $\mathrm{C}^{18}\mathrm{O}$ line emission is depleted because the regions along the exact disk midplane ($Z\sim0$) have the highest density and lower temperature conditions. Such phenomena have indeed been reported in some other young disks, such as L1489 IRS and IRAS 04302+2247 (e.g., Pinte et al. [118]; Yamato et al. [186]). Note that the column density map in Figure 7 was created under the assumption of a constant temperature ($T_{\mathrm{ex}}=50$ K).

Intensity decreasing towards the disk height direction as approaching the disk midplane becomes much clearer in the $\mathrm{H}_{2}\mathrm{CO}$ line emission map. Recognizing the much higher critical density of the $\mathrm{H}_{2}\mathrm{CO}$ line ($6\times10^{5}\,\mathrm{cm}^{-3}$ for $T_{\mathrm{ex}}=50$ K, Shirley [147]) compared to the $\mathrm{C}^{18}\mathrm{O}$ line, with similar upper state energies (see Table 1), the difference seen between the $\mathrm{C}^{18}\mathrm{O}$ line and the $\mathrm{H}_{2}\mathrm{CO}$ line maps could originate from a different distribution for these molecules. Since we do not know the temperature distribution of the disk structure in this object, it is difficult to make a firm argument. However, considering that $\mathrm{H}_{2}\mathrm{CO}$ has a higher sublimation temperature ($\sim40-70$ K, see e.g., Ferrero et al. [54]) than CO ($\sim20$ K), it is possible that $\mathrm{H}_{2}\mathrm{CO}$ is more depleted than CO in the disk midplane where its temperature is lower than that in the warm molecular layer of the disk. The left panel in Figure 4 shows that in the region within 200 au from the disk midplane, emissions from both $\mathrm{C}^{18}\mathrm{O}$ and $\mathrm{H}_{2}\mathrm{CO}$ become significantly stronger than at $Z\sim0$. These components might reflect the distribution of the molecular gas in the disk.

The observed line of c-C$_3$H$_2$ has a higher critical density and excitation than that of CCH (the critical densities for $T_{\mathrm{ex}}=50$ K are $5\times10^{7}\,\mathrm{cm}^{-3}$ for c-C$_3$H$_2$ and $8\times10^{6}\,\mathrm{cm}^{-3}$ for CCH, Chandra & Kegel [27]; Schöier et al. [142]; Spielfield et al. 2012), but their distributions do not differ significantly (see Figure 3). Therefore, it is reasonable to consider that these molecules are inherently abundant in this layer, rather than being bright due to excitation effects. Considering that the outermost layer is exposed to irradiation from the protostar and even radiation from outside the system, it is reasonable to interpret that photochemistry occurs in this layer. These molecules are indeed representative of photochemistry, as also reported with spatially resolved imaging in photodissociation regions (e.g., Cuadrado et al. [38]; Guzmán et al. [61]).

Overall, the layered structure found in the $\mathrm{C}^{18}\mathrm{O}$, $\mathrm{H}_{2}\mathrm{CO}$, c-C$_3$H$_2$, and CCH lines appear to originate from the different abundance of molecules in each layer, which can be referred to as a chemical layered structure. This chemical layered structure is found previously in Class I-II edge-on disks where gas accretion from the outer envelope is almost completed (e.g., Podio et al. [120]; Öberg et al. [102]). On the other hand, in IRAS 4C, the $\mathrm{H}^{13}\mathrm{CO}^{+}$ and DCO$^{+}$ line emission seems to trace not layered structure but disk/envelope system, as we described in § 3.3 and Figure 6.

Based on thermo-chemical models, protoplanetary disks consist of three chemical layers in the vertical direction (e.g., Aikawa et al. [3, 2]; Woitke et al. [185]; Dutrey et al. [46]; Walsh et al. [180, 182]; Notsu et al. [97]; Oberg et al. [102]). The first layer is the hot surface layer, directly irradiated by the stellar radiation and thus warm, and molecules are photodissociated into atoms and radicals by (inter-)stellar UV and X-ray radiation. Thus, this disk surface is a photon-dominated region (PDR). CCH and c-C$_3$H$_2$ become abundant in this region, especially in conditions where C/O $\gtrsim1$ (e.g., Bergin et al. [9]; Cleeves et al. [32]; Bosman et al. [10]). The second layer is the warm ($\gtrsim20$ K) molecular layer, where UV radiation is moderately shielded and molecules are in the gas phase and gas-phase chemistry is active. Ion-molecule reactions are triggered by X-rays from the central stars and cosmic-rays. The third layer is the freeze-out layer, that is, the cold outer disk midplane where molecules freeze out onto dust grains (beyond the snowlines). Even in the outer cold midplane, reactions in the gas phase and within ices can be triggered by high-energy cosmic-rays, X-rays, and decay of radioactive nuclei (e.g., $^{26}$Al). Various complex organic molecules, such as CH$_3$OH, only form on the icy mantles of dust grains, whereas unsaturated molecules such as carbon chain molecules can form by gas-phase reactions in the warm molecular and hot surface layers.

Here we compare our observed images in IRAS 4C with the vertical chemical layered structures obtained in other objects. van’t Hoff et al. (2023) reported the results of high-spatial resolution ($\sim24$ au) observations of molecular line emission towards the Class 0 protostar L1527 in ALMA Large Program eDisk and also discussed the presence of similar chemical layered structures between C$^{18}$O and H$_2$CO. They showed that C$^{18}$O, H$_2$CO, and SO trace the disk and inner envelope (see also Zhang et al. [193] for SO in L1527), with H$_2$CO emission arising predominantly in the disk surface layers. In contrast, they discussed that c-C$_3$H$_2$ is absent in the disk and its emission originates solely in the envelope (see also van’t Hoff et al. 2026), which is consistent with previous ALMA observations of different transitions (Sakai et al. [134]) and the CCH distribution (Sakai et al. [132]). We consider the possible reasons why CCH and c-C$_3$H$_2$ are absent in the surface layer of the disk in L1527, unlike IRAS 4C, that in L1527 (1) The vertical density gradient is steep and the PDR layer is geometrically thin, (2) As exemplified by the presence of SO emission (Zhang et al. [193]), oxygen is supplied to the gas phase, creating unfavorable conditions for formation of carbon-chain molecules (Bergin et al. [9]; Anderson et al. [5]) (3) The temperature is high and CCH + H$_2$ (with energy barrier) reaction proceeds.

Recently, spatially resolved observations of some Class I-II edge-on (=small inclination angle) disks have been obtained with CO and simple molecules (e.g., Dutrey et al. [45, 44]; Louvet et al. [84]; Podio et al. [120]; Teague et al. [155]; van’t Hoff et al. 2020b). Podio et al. [120] reported from their observations for the edge-on disk around the Class I source IRAS 04302+2247 (the butterfly star) by ALMA-DOT program (ALMA chemical survey of Disk-Outflow sources in Taurus, see e.g., Garufi et al. [60]) that they could resolve, the vertical structure of a disk with multiple molecular tracers (see also § 1 and van’t Hoff et al. 2020b; Lin et al. [78]). They showed that molecular emission is vertically stratified, with $^{12}$CO observed at larger disk heights compared to both CS and H$_2$CO, which is nearly co-spatial. They found that the height of molecular emitting regions scaled by radius ($z/r$) are $\sim0.41-0.45$ for $^{12}$CO and $\sim0.21-0.28$ for H$_2$CO at $r \sim100-200$ au. van’t Hoff et al. (2020b) also showed the vertically stratified structures in this disk from their ALMA observations that C$^{17}$O line emission originates in the disk midplane and H$_2$CO line emission traces more surface layers. They reported the presence of CO freeze-out in the outermost region ($\gtrsim100$ au) of the disk midplane. Lin et al. [78] also confirmed these chemical layered structures in the same disk using ALMA eDisk observations.

At first glance, the chemical layer structures observed in Class I/II disks and those of the Class 0 IRAS 4C disk look similar, but also have some differences as well. First of all, the spatial extent of the IRAS 4C emission is $r \sim300$ au along the midplane of the disk and envelope structures. The H$^{13}$CO$^+$ and DCO$^+$ line emission traces the infalling envelope structure from the outer envelope to the disk-forming region of IRAS 4C, which is clearly different from the Class I/II disks. It is difficult to clarify with only single source observations how continuous infall motion affects the formation of the chemically layered structure. Furthermore, the radiation in the disk and envelope surface in the Class 0 sources are dominated by accretion luminosity and protostellar outflow shocks, and not a stellar photosphere like Class I/II sources. Because the disk+envelope (Keplerian disk + pseudo disk) in IRAS4C is less dense than Keplerian disk in Class I/II sources, UV radiation from outflow shocks is transmitted to more deeper region in the vertical direction, and the “surface layer” might be geometrically thicker. To clarify the difference of the temperature structures in the disk vertical direction between IRAS4C and L1527 ((van’t Hoff et al. 2023; Zhang et al. 2024)), where CCH and c-C$_3$H$_2$ are absent in the disk surface layer, it would be useful to analyze the distribution of multiple line transitions of SO and SO$_2$ gas, which are formed by SO and SO$_2$ ice sublimation where gas is heated by UV radiation, or by H$_2$O photodissociation ([176]).

The Radial Variations of the Disk Ionisation Degree and Ionisation Rates

Disk ionisation is thought to be a fundamental parameter in the physical and chemical evolution of protostellar envelopes and disks. The non-ideal MHD (magnetohydrodynamic) effects, which depend on ionisation degree, play a major role in angular momentum transfer and mass accretion (e.g., [150, 160, 162, 73, 94, 194]). Recently, Kuffmeier et al. [73] showed from their non-ideal MHD simulations of a collapsing protostar that a higher ionisation rate leads to stronger magnetic braking, and hence to the formation of smaller disks. In particular, they speculate that the protostellar clusters in the ONC and CrA cause a higher cosmic-ray ionisation rate and thereby smaller disks in these regions than in regions such as Taurus, Ophiuchus, Chamaeleon I, or Lupus (see also Tobin et al. [159]). In addition, ionisation also triggers ion-molecule reactions which, for example, convert CO to other species, such as CO$_2$, CH$_4$, and CH$_3$OH (e.g., [59, 11, 47, 143, 98, 141]). There are several possible ionisation sources in disks and envelopes: X-rays and gamma-ray emission from the central star, cosmic-rays and stellar energetic particles, and the decay of short-lived radioactive nuclei (SLRs) (e.g., [164, 165, 163, 28, 29, 30, 31, 121, 122, 114, 113, 99, 58, 179, 71, 80, 81, 183, 161, 15, 40, 144]). The typical cosmic-ray (CR) ionisation rate $\zeta_{\mathrm{CR}}$ in nearby star-forming regions is $\sim10^{-16}$–$10^{-17}\ \mathrm{s}^{-1}$ (e.g., [125, 124, 129, 117]), and it is attenuated in the denser locations such as the disk midplane ([164, 58]), with the attenuation length of cosmic-rays much larger than for X-rays.

Aikawa et al. [4] analyzed HCO$^+$, H$^{13}$CO$^+$, N$_2$H$^+$, and N$_2$D$^+$, line emission from five Class II disks obtained by the ALMA Large Program “Molecules with ALMA at Planet-forming Scales (MAPS)” (e.g., [101]). They reported that the ionisation degree is low in the inner dense gas, especially within the CO snowline, and the midplane ionisation rate $\zeta_{\mathrm{CR}}$ is $\gtrsim10^{-18}\ \mathrm{s}^{-1}$ for the disks around the Class II objects IM Lup, AS 209, and HD 163296. Seifert et al. [140] reported from independent ALMA observations of H$^{13}$CO$^+$ and N$_2$H$^+$ line emission for the IM Lup disk that the distribution of line emission is best reproduced using a radially increasing cosmic-ray gradient, with low cosmic-ray ionisation in the inner disk ($\zeta_{\mathrm{CR}} \lesssim10^{-18}\ \mathrm{s}^{-1}$), high cosmic-ray ionisation in the outer disk ($\zeta_{\mathrm{CR}} \gtrsim10^{-17}\ \mathrm{s}^{-1}$), and a transition at $\sim80$–$100$ au. In addition, Long et al. [80] reported from their ALMA observations of HCO$^+$, H$^{13}$CO$^+$, and N$_2$H$^+$ line emission for the Class II DM Tau disk that the cosmic-ray ionisation rate is modestly reduced ($\zeta_{\mathrm{CR}} \gtrsim10^{-18}\ \mathrm{s}^{-1}$) and there is a hard X-ray spectrum produced by flares. They also discussed that their best-fit model under-produced emission in the inner disk, suggesting that there may be an additional mechanism enhancing ionisation in DM Tau’s inner disk, such as stellar energetic particles (e.g., [163, 121, 15]). Trapman et al. [161] estimated the radially averaged disk midplane cosmic-ray ionisation rates for 10 disks in ALMA Large Program exoALMA samples (e.g., [154]), by assuming the kinematically derived disk gas mass and refitting the N$_2$H$^+$ and CO line emission spectra. They reported a wide range of midplane ionisation rates, from $\zeta_{\mathrm{CR}} \sim10^{-17}\ \mathrm{s}^{-1}$ for AS 209 and RXJ 1852.3-3700 to $\zeta_{\mathrm{CR}} \sim3 \times10^{-20}\ \mathrm{s}^{-1}$ for Elias 2-27. They also reported the median cosmic-ray ionisation rate of $1.5 \times10^{-18}\ \mathrm{s}^{-1}$, although the uncertainty of ionisation rates for individual sources is large, typically $\sim1.5$ dex, suggesting that N$_2$H$^+$ alone only provides a moderate constraint on the disk ionisation rate.

van’t Hoff et al. (2022) discussed from the HCO$^+$ and H$^{13}$CO$^+$ line emission obtained by ALMA that the cosmic-ray ionisation rate has a canonical value of $\sim10^{-17}\ \mathrm{s}^{-1}$ in the inner envelope and may be attenuated to $\sim10^{-18}\ \mathrm{s}^{-1}$ in the disk around L1527. Ceccarelli et al. [24] reported the enhanced ionisation rate of $>10^{-14}\ \mathrm{s}^{-1}$ for the whole envelope and $\sim10^{-12}\ \mathrm{s}^{-1}$ for the innermost envelope around the cluster of a few intermediate and low-mass protostars in the source OMC-2 FIR 4 with a total mass of $30M_{\odot}$ (e.g., Shi majiri et al. 2008; López-Sepulcre et al. [83]; Tobin et al. [158] from their Herschel/HIFI observations of HCO$^{+}$ and N$_{2}$H$^{+}$, suggesting the presence of stellar energetic particles whose energies are $> 10$ MeV. Recently, Cabedo et al. [17] reported the first resolved map of the cosmic-ray ionisation rate at $10^{2}-10^{3}$ au scales in a solar-type Class 0 protostar, and found an enhanced ionisation that increases at small envelope radii, reaching values up to the cosmic-ray ionisation rate of $\sim10^{-14}$ s$^{-1}$ at a few hundred astronomical units (au) from the central Class 0 protostellar object B335. They showed that the extremely high cosmic-ray ionisation rate can be attributed to the presence of cosmic-rays accelerated close to the protostar (e.g., Turner & Drake [163]; Rab et al. [121]; Brunn et al. [15]). It is therefore important to resolve the radial distribution of the ionisation degree and ionisation rate from the envelope to inner disk scales ($< 100$ au) in Class 0/I objects.

In this study, we estimate the distribution of ionisation degree and ionisation rates around IRAS 4C from the observed C$^{18}$O, H$^{13}$CO$^{+}$, and DCO$^{+}$ line emission. Figure 8a shows the distribution of the column density ratio $R_{\rm H}=N_{\rm HCO^{+}}/N_{\rm CO}$ (CO ionisation degree). To estimate $N_{\rm HCO^{+}}$ from $N_{\rm H^{13}CO^{+}}$ and $N_{\rm CO}$ from $N_{\rm C^{18}O}$, we assume the $^{12}$C/$^{13}$C ratio as 69 and the $^{16}$O/$^{18}$O ratio as 560 (Wilson [184]; Jørgensen & van Dishoeck [70]). We note that the carbon isotope fractionation in cold envelope with higher ionisation rate ($\gtrsim5 \times10^{-17}$ s$^{-1}$) would decrease $^{12}$C/$^{13}$C ratios in CO and HCO$^{+}$ by a factor of $2-3$ (e.g., Bergin et al. [7]). Thus, the values of $N_{\rm HCO^{+}}$ might be overestimated by a factor of $2-3$. In Figure 8a, $R_{\rm H}$ is $\sim10^{-5}-10^{-4}$ at $R \gtrsim0.2^{\prime\prime}-0.7^{\prime\prime}$ ($R \sim50$-$200$ au) and $(1-2) \times10^{-4}$ at $R \gtrsim0.7^{\prime\prime}$ ($R \gtrsim200$ au) from the central protostar, respectively (see also Figure 9). Aikawa et al. [4] estimated that $R_{\rm H}$ is $\sim10^{-5}-10^{-4}$ at $R \gtrsim100$ au, and $\lesssim10^{-5}$ at $R \lesssim100$ au in Class II disks. Thus, both Class II disks and the IRAS 4C disk/envelope system have the radial gradients of $R_{\rm H}$, although the values in the same radii are different.

We note that according to Figure A1 in Appendix § 6.1, the H$^{13}$CO$^{+}$ line emission is marginally optically thick ($\tau\sim1-1.5$) at $r \sim100-200$ au, and C$^{18}$O is optically thick $\tau\gtrsim2$ in $r \lesssim150$ au. Thus, $R_{\rm H}$ might be overestimated by a factor of around $2-3$ at $r \lesssim150$ au since the column densities of C$^{18}$O might be underestimated compared to those of H$^{13}$CO$^{+}$.

Based on van’t Hoff et al. (2018a), by assuming HCO$^{+}$ is formed by the reaction between CO and H$_{3}^{+}$ and $N_{e^{-}}=N_{\rm HCO^{+}}$, the ionisation rate $\zeta_{\rm CR}$ is estimated by the following equation,

where $L$ is the depth of the structure, $k_{1}$ is the reaction rate of HCO$^{+}$ and e$^{-}$ and $k_{1}=2.40 \times10^{-7} \times(T_{\rm ex}/300)^{-0.69}$ cm$^{3}$ s$^{-1}$ (= $8.3 \times10^{-7}$ cm$^{3}$ s$^{-1}$ at $T_{\rm ex}=50$ K, Brian & Mitchell [14]). We note that this assumption of the relation between HCO$^{+}$ and CO is basically reasonable outside the water snowline, whereas HCO$^{+}$ molecules are significantly destroyed by water gas within the water snowline (van’t Hoff et al. 2018a; Notsu et al. [99]). Thus, within the water snowline, $\zeta_{\rm CR}$ could be underestimated. We assumed that the size of the NGC 1333 cloud and the structure of IRAS 4C are the same, and assumed $L=6 \times10^{4}$ au based on the assumption by Pineda et al. [117] for the NGC 1333 cloud. Here we note that if we assume a much smaller value of $L$ (=size of the structure), the ionisation rates increase in inverse proportion to $L$ (see below). Figure 8c shows the distribution of $\zeta_{\rm CR}$ estimated from $N_{\rm H^{13}CO^{+}}$ and $N_{\rm H_{2}}$. In addition, Figure 8d shows the distribution of $\zeta_{\rm CR}$ estimated from $N_{\rm H^{13}CO^{+}}$ and $N_{\rm C^{18}O}$, by assuming the $^{16}$O/$^{18}$O ratio as 560 and $N_{\rm CO}/N_{\rm H_{2}}=10^{-4}$. In Figure 8c and (Figure 8d), we mask the innermost region with significantly higher dust optical depth by the black-filled ellipse, whose size and shape are the same as the beam of the dust continuum emission. We assumed that CO depletion has not efficiently happened in the Class 0 object IRAS 4C, based on previous observations (e.g., Zhang et al. [192]). Although the actual values of $\zeta_{\rm CR}$ are different by a factor of $2-3$ between these two maps, $\zeta_{\rm CR}$ is larger in the outer region and smaller in the inner region in both panels.

As described in § 3.1, the dust optical depth is larger in the inner region of the disk/envelope system, especially within 50 au of the protostar. We note that the value “50 au” has an uncertainty, since it is similar to the beam size of our observations (see Table 2). Considering the high dust optical depth at $R<50$ au and beam sizes of our observations ($\sim0.^{\prime\prime}2$, see Table 2), in this innermost region at $R<50$ au, it is not appropriate to discuss the ionisation rates by comparing the distribution of C$^{18}$O, H$^{13}$CO$^{+}$, and DCO$^{+}$ line emission (see Figure 4). Thus, in Figure 8 and Figure 9 we consider only $R>50$ au. Figure 9 shows the radial profile of $R_{\rm H}=N_{\rm HCO^{+}}/N_{\rm CO}$ and the ionisation rate estimated from the H$^{13}$CO$^{+}$ and H$_{2}$ with $L=6 \times10^{4}$ au along the radial (disk major axis) direction. We note that in deriving the CO ionisation degree, we assume that both C$^{18}$O and H$^{13}$CO$^+$ line emission trace the same region along the line of sight. Focusing on the region with $R \sim50 - 200$ au (see Figure 4 and Figure 9), we can see that the relationship between the line emission intensities of the C$^{18}$O line and ion molecules (H$^{13}$CO$^+$ and DCO$^+$) is different from that at $R > 200$ au, and thus we can also see that the $R_{\mathrm{H}}$ is changing around $R \sim200$ au ($R_{\mathrm{H}} \sim10^{-5} - 10^{-4}$ at $50 < R < 200$ au and $R_{\mathrm{H}} \sim(1 - 2) \times10^{-4}$ at $R > 200$ au, see Figure 8a and Figure 9). We consider that the gas densities traced by the CO emission line increase within $R \sim200$ au. It might suggest that the disk within $\sim200$ au is Keplerian disk whereas the disk outside $\sim200$ au is pseudo disk where infall motion is dominant. In Figure 6, the C$^{18}$O line emission on the PV diagram increases at $\lesssim0.7^{\prime\prime}$ ($\lesssim200$ au), suggesting a significant increases in gas density due to the presence of a Keplerian disk. However, further detailed analyses with infalling envelope models (e.g., Oya et al. [110]) are required to constrain the Keplerian disk radius more precisely.

In Figure 8 and Figure 9, the absolute values of ionisation rates $\zeta_{\mathrm{CR}}$ in the IRAS 4C disk/envelope system are estimated to be $\sim10^{-17} - 10^{-16}$ s$^{-1}$ at $R > 100$ au and $\lesssim10^{-18}$ s$^{-1}$ at $50 < R < 100$ au. Pineda et al. [117] reported the maps of ionisation rates in the larger scale of NGC 1333 star-forming region, and the ionisation rates of the clouds around IRAS 4C ($\sim$ a few $10^3$ au) is $\sim10^{-16}$ s$^{-1}$, which is consistent with the value for the outer envelope of IRAS 4C in Figure 8. We have confirmed that the absolute value of ionisation rates for IRAS 4C varies within a factor of $\sim2$ by changing $T_{\mathrm{ex}}$ by 20K and 100K (see Figure B1 and Figure B2 in Appendix § 6.2). We consider that the disk formation may have increased the local gas density and decreased the ionisation rate in the inner dense midplane regions of the disk/envelope system. In addition, apart from B335 (Cabedo et al. [17]), the additional cosmic-rays acceleration close to the protostar would not be effective for the disk ionisation at the regions with $\gtrsim50 - 100$ au of IRAS 4C.

Image not extracted

Image not extracted

However, for the estimated values of ionisation rates in IRAS 4C, it has an uncertainty up to orders of magnitude based on the assumptions used to estimate the ionisation fraction (Sabatini et al. [129]). The assumption of $N_{\mathrm{e}^{-}} = N_{\mathrm{HCO}^{+}}$ and the value of $L$ we adopted also brings uncertainty. If we assume a much smaller value of $L$ (=depth of the structure), the ionisation rates increase in inverse proportion to $L$. We do not know the exact value of the depth of IRAS 4C structure, but if we assume $L$ as around a few $\times10^3$ au, the ionisation rates increase by a factor of a few $\times10$. Figure B3 in Appendix § 6.2 shows the distributions of the ionisation rate by assuming $L = 1000$ au, which is the size of relatively smaller protostellar envelopes (such as in B335, Cabedo et al. [17]). In Figure B3, the values of ionisation rates are estimated to be $5 \times10^{-16} - 10^{-14}$ s$^{-1}$ at $R > 100$ au and $\lesssim10^{-16}$ s$^{-1}$ at $50 < R < 100$ au.

Image not extracted

Recent studies (e.g., [13, 128, 129, 123, 124, 148, 74]) proposed a new analytical approach to estimate the ionisation rates by using H$_2$D$^+$ lines as the main observational constraint to derive the amount of H$_3^+$. Future comparison of the line emission of multiple ionisation tracers (such as HCO$^+$, H$^{13}$CO$^+$, N$_2$H$^+$, H$_2$D$^+$, and N$_2$D$^+$) with more detailed disk models including chemical reaction networks and temperature structures will help to clarify the distribution of ionisation rates in the IRAS 4C disk-envelope system (see e.g., Luo et al. [85]).

The top right panel in Figure 8 shows the deuteration fraction $R_{\mathrm{D}}=N_{\mathrm{DCO}^{+}}/N_{\mathrm{HCO}^{+}}$. According to this Figure, $R_{\mathrm{D}}$ is $\sim0.4 - 0.5\%$ at $R \gtrsim0.5^{\prime\prime}$ ($R \gtrsim150$ au) and $\lesssim0.1\%$ at $R \lesssim0.5^{\prime\prime}$ ($R \lesssim150$ au). These values of $R_{\mathrm{D}}$ are around one order of magnitude smaller than those at the $\sim1000$ au scale envelope around B335 (Cabedo et al. [17]) and prestellar cores (Caselli et al. [22]). We try to estimate $\zeta_{\mathrm{CR}}$ from DCO$^+$, C$^{18}$O, and H$^{13}$CO$^+$ line emission based on the method of Caselli et al. [22] and Cabedo et al. [17], but the estimated value of $\zeta_{\mathrm{CR}}$ is unreasonably high ($\sim10^{-12} - 10^{-11}$ s$^{-1}$). Thus, apart from the lower density molecular clouds and prestellar cores (Caselli et al. [22]; Fontani et al. [55]), $R_{\mathrm{D}}=N_{\mathrm{DCO}^{+}}/N_{\mathrm{HCO}^{+}}$ would not be an appropriate probe of the ionisation degree in higher density disks and envelopes around protostars. The observations of other deuterated molecules (such as H$_2$D$^+$, D$_2$H$^+$, N$_2$D$^+$, and DCN) would be helpful to estimate the ionisation rates correctly (e.g., [21, 52, 13, 129, 148, 123, 124]), although the deuteration chemistry of H$_3^+$ (especially H$_2$D$^+$/D$_2$H$^+$ ratio) is not fully understood (Pagani et al. [115]). We note that Koumpia et al. [72] reported from their single-dish JCMT⁹ observations that the deuteration fraction estimated from DCO$^+$ and HCO$^+$ emission lines is $\sim12 \pm2\%$ toward IRAS 4A, $\sim3.5 \pm1\%$ toward IRAS 4B, and an upper limit of $\sim20\%$ toward IRAS 4C. They also pointed out that the outer envelope of IRAS 4C should show the higher degree of deuteration since it is the coldest among the sources and thus the enhancement of $\mathrm{H_2D^+}$ abundance contributes to the enhancement of other deuterated molecules.

Conclusion

We report the results of molecular line observations of the disk and envelope (Keplerian disk + pseudo disk) around the Class 0 low-mass protostar NGC 1333 IRAS 4C in the Perseus Molecular Cloud, as part of the ALMA Large Program FAUST. Previous observations (ALMA/PEACHES) have shown that the disk rotation axis of this source is largely tilted to the line of sight, making it possible for studying the molecular abundance distribution along the disk height. We found that the molecules generated by photochemical reactions such as CCH and c-$\mathrm{C_3H_2}$ are distributed in the low-density upper layers of the disk/envelope system, while the molecules such as $\mathrm{C^{18}O}$ and $\mathrm{H_2CO}$ are distributed in the high density regions around the protostar, indicating a chemically layered structure in the disk vertical direction. The chemical layered structure is similar to those known for Class I/II objects, but also had some differences. In addition, we compared the differences of chemical layered structure between IRAS 4C and another Class 0 object L1527, and the difference may suggest that compared with L1527, the surface layer of the disk in IRAS 4C are more active conditions for formation of carbon-chain molecules such as CCH and c-$\mathrm{C_3H_2}$ and/or the PDR layer is not geometrically thin. Moreover, the distribution of ionisation degree and ionisation rate in the radial direction of the disk/envelope system is estimated from the comparison of $\mathrm{H^{13}CO^+}$ and $\mathrm{C^{18}O}$ line emission. The $R_{\mathrm{H}}(=N_{\mathrm{HCO^+}}/N_{\mathrm{CO}})$ is $\sim10^{-5}-10^{-4}$ at $R\sim0.2''-0.7''$ ($R\sim50$–200 au) and $(1-2)\times10^{-4}$ at $r\gtrsim0.7''$ ($\gtrsim200$ au) from the central protostar, respectively. The estimated values of the ionisation rate $\zeta_{\mathrm{CR}}$ are $\sim10^{-17}-10^{-16}\ \mathrm{s^{-1}}$ at $R>100$ au, and $\lesssim10^{-18}\ \mathrm{s^{-1}}$ at $50<R<100$ au, although the absolute values of ionisation rates have uncertainty (up to orders of magnitude) based on the assumptions invoked to estimate the ionisation fraction and the depth of the structure. Disk formation may have increased the local gas density and decreased the ionisation degree and the ionisation rate in the inner dense midplane regions of the disk/envelope system.

This paper makes use of the following ALMA data set: ADS/JAO.ALMA# 2018.1.01205.L (FAUST, PI: Satoshi Yamamoto). ALMA is a partnership of the ESO (representing its member states), the NSF (USA) and NINS (Japan), together with the NRC (Canada) and the NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by the ESO, the AUI/NRAO, and the NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The authors thank the ALMA staffs for their excellent support. S.N. is grateful for support from RIKEN Special Postdoctoral Researcher Program (Fellowships), Start-up Research Grant for The University of Tokyo Excellent Young Researcher 2024, Grants-in-Aid for JSPS (Japan Society for the Promotion of Science) Fellows Grant Number JP23KJ0329, and MEXT/JSPS Grants-in-Aid for Scientific Research (KAKENHI) Grant Numbers JP20K22376, JP23K13155, JP24K00674, and JP26K00755. N.S. is supported by MEXT/JSPS Grants-in-Aid for Scientific Research (KAKENHI) Grant Numbers JP20H05844, JP20H05845, JP20H00182, and JP25H00676, and a pioneering project in RIKEN (Evolution of Matter in the Universe). Y.A. is supported by MEXT/JSPS Grants-in-Aid for Scientific Research (KAKENHI) Grant Number JP20H05847 and JP24K00674. Y.Z. acknowledges the support from the Yangyang Development Fund. H.B.L. is supported by the National Science and Technology Council (NSTC) of Taiwan (Grant Nos. 111-2112-M-110-022-MY3, 113-2112-M-110-022-MY3). C.C., L.P., and E.B acknowledge financial support under the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.1, Call for tender No. 104 published on 2.2.2022 by the Italian Ministry of University and Research (MUR), funded by the European Union - NextGenerationEU-Project Title 2022JC2Y93 Chemical Origins: linking the fossil composition of the Solar System with the chemistry of protoplanetary disks - CUP J53D23001600006 - Grant Assignment Decree No. 962 adopted on 30.06.2023 by the Italian Ministry of University and Research (MUR).

E.B. acknowledges financial support from the Italian Ministry of University and Research (MUR) under the Italian Science Fund (FIS 2 Call - Ministerial Decree No. 1236 of 1 August 2023).

Cl.Co., L.P., E.B., N.B., and G.S. also acknowledge the PRIN-MUR 2020 BEYOND-2p (Astrochemistry beyond the second period elements, Prot. 2020AFB3FX), the project ASI-Astrobiologia 2023 MIGLIORA (Modeling Chemical Complexity, F83C23000800005), the INAF-GO 2023 fundings PROTO-SKA (Exploiting ALMA data to study planet forming disks: preparing the advent of SKA, C13C23000770005), the INAF-MiniGrant 2022 “Chemical Origins” (PI: L. Podio), and the INAF-MiniGrant 2023 TRIESTE (“Tracing the chemical Heritage of our originS: from proTOstars to planEts”; PI: G. Sabatini). S.B.C. was supported by the NASA Planetary Science Division Internal Scientist Funding Program through the Fundamental Laboratory Research work package (FLaRe). ALS, CeCe, and CV acknowledge support by the Thematic Action “Physique et Chimie du Milieu Interstellaire” (PCMI) of INSU Programme National “Astro”, with contributions from CNRS Physique & CNRS Chimie, CEA, and CNES. N.C. was supported by the funding from the European Research Council (ERC) under the European Union Horizon Europe programme (grant agreement No. 101042275, project Stellar-MADE). D.J. is supported by NRC Canada and by an NSERC Discovery Grant. S.V and M.B. acknowledge the support from the European Research Council (ERC) Advanced Grant MOPPEX 833460. F.Me. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon Europe research and innovation program (grant agreement No. 101053020, project Dust2Planets). L.L. acknowledges the support of DGAPA/UNAM PAPIIT grants IN108324 and IN112820, CONAHCyT grant CF-263356 and SECHIHTI grant CBF-2025-I-109.

Appendix

Distribution of line optical depth

Figure A1 shows the distribution of line optical depth $\tau$ for H$^{13}$CO$^{+}$ (a), DCO$^{+}$ (b), and C$^{18}$O (c). Our adopted method for calculating line optical depth is described in § 3.4.

The uncertainty of our derived values of ionisation rates

Figure B1–Figure B3 show the uncertainty of our derived values of ionisation rates $\zeta_{\mathrm{CR}}$ (see § 4.2), by changing $T_{\mathrm{ex}}$ (Figure B1 and Figure B2) and the depth of the structure $L$ (Figure B3).