Early Planet Formation in Embedded Disks (eDisk). II. Limited Dust Settling and Prominent Snow Surfaces in the Edge-on Class I Disk IRAS 04302+2247

While dust disks around optically visible, Class II protostars are found to be vertically thin, when and how dust settles to the midplane are unclear. As part of the Atacama Large Millimeter / submillimeter Array large program, Early Planet Formation in Embedded Disks, we analyze the edge-on, embedded, Class I protostar IRAS 04302 + 2247, also nicknamed the “ Butter ﬂ y Star. ” With a resolution of 0 05 ( 8 au ) , the 1.3 mm continuum shows an asymmetry along the minor axis that is evidence of an optically thick and geometrically thick disk viewed nearly edge-on. There is no evidence of rings and gaps, which could be due to the lack of radial substructure or the highly inclined and optically thick view. With 0 1 ( 16 au ) resolution, we resolve the 2D snow surfaces, i.e., the boundary region between freeze-out and sublimation, for 12 CO J = 2 – 1, 13 CO J = 2 – 1, C 18 O J = 2 – 1, H 2 CO J = 3 0,3 – 2 0,2 , and SO J = 6 5 – 5 4 , and constrain the CO midplane snow line to ∼ 130 au. We ﬁ nd Keplerian rotation around a protostar of 1.6 ± 0.4 M e using C 18 O. Through forward ray-tracing using RADMC-3D, we ﬁ nd that the dust scale height is ∼ 6 au at a radius of 100 au from the central star and is comparable to the gas pressure scale height. The results suggest that the dust of this Class I source has yet to vertically settle signi ﬁ cantly.


INTRODUCTION
The formation of rotationally supported circumstellar disks plays a crucial role in the star and planet formation process.As a consequence of the conservation of angular momentum, much of the material from the larger scale core is channeled to the disk and subsequently accretes onto the protostar itself (e.g.Terebey et al. 1984;Li et al. 2014;Tsukamoto et al. 2022).The reservoir of material in the disk enables the growth of solids and serves as the birthplace of planets (e.g.Testi et al. 2014;Drazkowska et al. 2022;Tu et al. 2022).Nevertheless, the process of dust evolution, from sub-micron-sized particles inherited from the core to planetesimals and planets, requires numerous mechanisms to overcome multiple growth barriers, e.g., the meter-sized barrier (Weidenschilling 1977).One of the most favored mechanisms to overcome the meter-sized barrier is the streaming instability, which can drive rapid growth from pebbles to planetesimals, but it requires comparable densities of the dust and the gas rather than the 1:100 dust-to-gas ratio inherited from the interstellar medium (e.g.Youdin & Goodman 2005;Lesur et al. 2022).One natural process to increase the dust-to-gas ratio is through dust settling (e.g.Gole et al. 2020).
Gaseous disks are vertically extended owing to the vertical pressure support.The balance between the pressure gradient and the vertical gravitational pull sets the gas scale height.In contrast, dust particles, if decoupled from the pressure-supported gas, will inevitably descend to the midplane to form a thin dust layer.Turbulent mixing operates against dust settling by stirring up the dust and prevents the dust from becoming fully settled (e.g.Nakagawa et al. 1986;Dubrulle et al. 1995).While the tendency for settling is well established, the effectiveness of turbulence is not clear and relies on observations for constraints (e.g.Pinte et al. 2016;Ohashi & Kataoka 2019;Villenave et al. 2022).However, observations that can characterize the vertical structure of disks are few in number, since it requires high angular resolution of nearly edge-on disks (Tobin et al. 2010;Lee et al. 2017;Sakai et al. 2017;Lee et al. 2020;Villenave et al. 2020;Michel et al. 2022;Ohashi et al. 2022).
IRAS 04302+2247 (hereafter IRAS 04302) is a Class I (bolometric temperature T bol = 88 K; Ohashi et al. in prep.)protostar, poetically nicknamed the "Butterfly Star" by Lucas & Roche (1997) for its remarkable bipolar reflection nebulae in the near-infrared.Highresolution near-infrared images from the Hubble Space Telescope/NICMOS exhibited a clear dark lane sandwiched between the reflection nebulae and depict a highly inclined system with an obscured central source and bipolar cavity walls that scatter the near-infrared photons (Padgett et al. 1999).A molecular outflow in H 2 was detected in the same direction as the bipolar cavity walls (Lucas & Roche 1998), and the deep absorption silicate feature in the mid-infrared requires significant inclination (Furlan et al. 2008).Indeed, millimeter wavelength observations show an elongated continuum within the near-infrared dark lane, which is evidence of the presence of an edge-on disk (Wolf et al. 2003(Wolf et al. , 2008;;Sheehan & Eisner 2017;van 't Hoff et al. 2020;Villenave et al. 2020).The near edge-on disk orientation facilitates the determination of the geometrical thickness of the dust layer (e.g.Villenave et al. 2020).
Detailed models of IRAS 04302 using scattered light images and the millimeter continuum images, which trace different physical processes and regions of the circumstellar system, have ascertained an inclined system of a disk and envelope (e.g.Wolf et al. 2003;Furlan et al. 2008;Wolf et al. 2008;Eisner 2012;Sheehan & Eisner 2017).Intriguingly, the dust in the envelope is consistent with interstellar medium (ISM) grains (e.g.Lucas & Roche 1997), while the dust in the disk is found to have grown significantly (Wolf et al. 2003;Gräfe et al. 2013;Sheehan & Eisner 2017).Furthermore, Gräfe et al. (2013) suggested that the larger grains in the disk show evidence of radial and vertical decoupling from the small grains.
Recent molecular line observations with ∼ 0.3 to 0.4 achieved by ALMA have begun to resolve the locations where molecules trace the disk surface of IRAS 04302, making the study of its vertical structure possible (van 't Hoff et al. 2020;Podio et al. 2020).van 't Hoff et al. (2020) identified C 17 O in the midplane within 100 au and detected emission in the disk surface layers beyond 100 au which can be explained by freeze-out of CO.In addition, H 2 CO (3 1,2 − 2 1,1 ) mainly originates from the disk surface layers with a large reduction of emission at the midplane where the continuum is located.Podio et al. (2020) also found a similar distribution of emission for 12 CO (2−1), H 2 CO (3 2,1 −2 1,1 ), and CS (5−4).The pattern is consistent with results from thermochemical models that consist of a midplane freeze-out and an elevated molecular layer separated by a snow surface (e.g.Aikawa et al. 2002;Akimkin et al. 2013;Dutrey et al. 2014).Furthermore, different isotopologues of CO trace different densities (van 't Hoff et al. 2020;Podio et al. 2020).
Most of the prior continuum observations have been limited in angular resolution with ∼ 0.2 to 0.5 making it difficult to resolve the vertical structure of the dust (Gräfe et al. 2013;Podio et al. 2020;van 't Hoff et al. 2020).The highest angular resolution of the continuum to date is ∼ 0.06 at λ = 2.1 mm and hints at a flared dust disk (Villenave et al. 2020).The unique view of IRAS 04302 thus serves as a perfect laboratory to study the vertical structure of the dust and gas around a young source in detail.As part of the Early Planet Formation in Embedded Disks (eDisk) program, we present highresolution λ = 1.3 mm continuum (∼ 0.05 or 8 au) and molecular line images (∼ 0.1 or 16 au) obtained from ALMA.
IRAS 04302 is located within the L1536 cloud of the Taurus star-forming region.The whole Taurus starforming region is conventionally assumed to have a distance of 140 pc (Kenyon et al. 1994), but recent parallax measurements found significant depth effects for each cloud.From Gaia, Luhman (2018) and Roccatagliata et al. (2020) found a distance of 161 and 160.3 pc, respectively, for the L1536 cloud.Galli et al. (2018) inferred a distance of 162.7 pc using astrometry from the Very Long Baseline Array.For this paper, we adopt a distance of 160 pc.
The rest of the paper is organized as follows.Section 2 describes the observations and data processing, while Section 3 shows the resulting dust continuum images and molecular line channel maps.We analyze the continuum and line data in more detail in Section 4. We discuss several implications in Section 5 and conclude in Section 6.

OBSERVATIONS
The data are obtained as part of the ALMA Large Program (2019.1.00261.L, PI: N. Ohashi).The details of the survey, including the spectral setup, calibrators, and imaging procedure, are discussed in Ohashi et al. (in prep.).We briefly describe the relevant setup for IRAS 04302.The short baseline data of IRAS 04302 were observed on Dec. 21, 2021 in configuration C-5 with baselines ranging from 15 m to 3.6 km with an onsource integration time of ∼ 35 minutes.The long baseline data were observed on Sept. 30 and Oct. 1 in 2021 with total integration times of ∼ 2.16 hours in configuration C-8 with baselines ranging from 70 m to 11.9 km.The spectral setup was in Band 6 with a representative wavelength of 1.3 mm (225 GHz) for the continuum.The spectral resolution for each detected line is listed in Table 1.
All calibration and imaging tasks utilized the Common Astronomy Software Applications (CASA) package (McMullin et al. 2007) version 6.2.1 and pipeline version 2021.2.0.128.From the pipeline calibrated data, we follow the self-calibration procedure presented in Ohashi et al. (in prep.) which we briefly describe in the following.First, we imaged each execution block separately and aligned the peaks to a common phase center using the fixvis and fixplanets tasks.Second, to adjust for flux calibration uncertainties between each execution block, we scaled the amplitude of the visibilities that were azimuthally binned as a function of uv-distance.We self-calibrated the short-baseline data through three rounds of phase-only calibration.With the self-calibrated short-baseline data, we included the long-baseline data and conducted one round of phaseonly calibration with a solution interval that was the length of each execution block.
We used the tclean task to image the self-calibrated visibilities.The continuum imaging used several Briggs robust weightings from robust=-2 to 2 (Briggs 1995).Smaller robust values correspond to better angular resolution at the expense of increased noise, while larger robust values correspond to better sensitivity albeit with lower angular resolution (e.g., Briggs 1995;Czekala et al. 2021).We show the resulting images in Appendix A. We adopt the image with robust=0.5 as the representative image to compromise between spatial resolution, sensitivity, and image fidelity.
The self-calibration solutions were applied to the measurement set used for the lines and further continuum subtracted using the uvcontsub task.Each line image cube used a robust=0.5 and 2 with the uvtaper set at 2000kλ (or ∼ 0.09 ).The self-calibration and imaging scripts for this source can be found at http: //github.com/jjtobin/edisk.We assume a 10% absolute flux calibration uncertainty, but we only consider the statistical uncertainty for the rest of this paper.The resulting resolution and noise levels for each image are listed in Table 1.
The CLEAN process for lines results in 2D images as a function of frequency ν (i.e., an image cube or channel maps).From the image cube, we can define several 2D quantities to interpret the 3D data.We denote the image cube as I(x, y, ν), where x and y represent the sky coordinates, R.A. and Dec.With a known line transition frequency ν 0 , one can convert from the observed frequency ν to the velocity along the line-of-sight, v, through where c is the speed of light.For the rest of the paper, we express v in km s −1 .The channel width in velocity units ∆v is related to the channel width in frequency units ∆ν through where ∆v and ∆ν are both positive quantities.∆v is shown in Table 1.
Since the image is defined on discrete pixels and spectral channels, we use I i,j,k , where i, j, k are indices, to represent the intensity value at a certain pixel (x i , y j ) and at a certain velocity v k .The integrated intensity image M is defined as (3) The resulting two-dimensional quantity only depends on the sky coordinates and the units are in Jy beam −1 km s −1 .
The peak intensity image P results from taking the peak along the spectrum at each coordinate (x i , y j ) expressed as where "max" represents taking the maximum along the velocity axis with N v points.
The intensity I ν in units of Jy beam −1 can be converted to brightness temperature in kelvins through the Planck function.Suppose the major and minor axes of the beam are θ M and θ m , respectively.The solid angle of the beam is Ω = πθ M θ m /(4 ln 2).Let J ν be the intensity I ν expressed in erg s −1 cm −2 Hz −1 ster −1 which are related by J ν = I ν 10 −23 /Ω.The peak intensity image P can express in brightness temperature with Using Eq. ( 5), we can express P as a brightness temperature.
With the line-of-sight velocities known, one can analyze the velocity structure of the source by extracting a representative velocity at each pixel of the image cube.We use the "peak velocity image," V , which is the velocity that corresponds to the peak of the spectrum at each pixel.The peak velocity map does not rely on an assumption on the profile shape and is shown to be less susceptible to noise compared to other methods (e.g. de Blok et al. 2008;Teague & Foreman-Mackey 2018).
The images, M , P , and V , were created using the CASA task immoment setting the argument moments to 0, 8, and 9, respectively. 1Furthermore, we only consider emission above the 3σ level to avoid "negative" intensities from continuum oversubtraction (see Table 1 for the noise levels).

Continuum
Fig. 1 shows the 1.3 mm continuum image with ro-bust=0.5 and reveals a highly elongated structure that is consistent with past low angular resolution images at millimeter wavelengths (Wolf et al. 2003(Wolf et al. , 2008;;Gräfe et al. 2013).The image has a peak of 1.11 mJy beam −1 with a noise level of σ =14.5 µJy beam −1 .The total flux is 184.15 mJy by integrating the emission above 3σ.The image appears largely symmetric along the major axis, but clearly asymmetric along the minor axis in which the east side is brighter than the west side.The elongated emission is expected from an inclined disklike structure and the kinematic analysis in Section 4 confirms a Keplerian disk.Thus, we will refer to the elongated continuum as simply the (dust) disk.Even with the higher angular resolution compared to previous observations, there is no clear evidence of rings or gaps.
To characterize the continuum image, we fit the disk with a 2D Gaussian using the CASA task imfit.The coordinate center of the 2D Gaussian is one of the free parameters, and we get the best-fit value of (04:33:16.50,+22:53:20.2) in ICRS, which we set as the origin of the image hereafter unless explicitly stated otherwise.We treat the center as the location of the star.The deconvolved full width at half maximum (FWHM) for the major and minor axes are 2.149 ± 0.007 and 0.2385 ± 0.0007 respectively.Assuming a completely flat disk, the ratio between the minor and major axes equals cos i where i is the inclination of the disk (i = 0 • means face-on).With the FWHM from the 2D Gaussian fitting, we derive i ∼ 84 • .Since the disk has a finite vertical thickness, the inclination estimation is a lower limit (see Section 4.1).The position angle (PA) of the major axis of the best-fit 2D Gaussian is 174.77• ± 0.03 • which we adopt as the position angle of the major axis of the system.The total flux from the fitting is 182.6±0.6 mJy (1σ uncertainty).Fig. 2 compares the major and minor axis cuts with the origin set at the center determined from the fitted 2D Gaussian.The cuts are produced by interpolating the image and we also calculate the brightness temperature T b using Eq. ( 5).The brightness temperature is low across the disk with only ∼ 14 K at the peak.For comparison, the peak brightness temperatures at λ = 0.9 mm (ALMA Band 7) and λ = 2.1 mm (ALMA Band 4) are 10 and 6.7 K, respectively (Villenave et al. 2020).The slightly higher peak brightness temperature presented here is likely because the disk is better resolved.The extent of the major axis reaches up to ∼ 2 (320 au) from the center, which is similar to the Band 4 and 7 continuum images from Villenave et al. (2020).The large extent implies a fairly large disk radius, which we constrain in Section 4.1.
To see the asymmetry along the minor axis clearly, we zoom in on the minor axis cut and show a comparison with the beam in the right panel of Fig. 2. The FWHM of the minor axis is resolved by ∼ 3.5 beams.The asymmetry could be due to an intrinsically asymmetric disk or due to a highly inclined axisymmetric disk that is optically thick, has a finite geometrical thickness, and is not seen exactly edge-on.We favor the latter possibility since the asymmetry occurs along the minor axis and is readily consistent with the high inclination and with the direction of the outflow (see Section 3.2) Given that the emission is brighter on the east side, we can infer that the east side is the far side of the disk based on simple expectations of an optically thick disk with decreasing temperature as a function of radius (Lee et al. 2017;Villenave et al. 2020;Ohashi et al. 2022;Takakuwa et al. in prep.)demonstrated through detailed modeling in Sec-tion 4.1).In addition, the optically thinner λ = 2.1 mm (ALMA Band 4) image with similar resolution (∼ 0.06 ; 10 au) does not show a similar asymmetry (Villenave et al. 2020), which is more consistent with our picture than an intrinsically asymmetric disk.
By assuming the emission at ν = 225 GHz comes entirely from the dust thermal emission and is optically thin, one can estimate the total dust mass disk through where S ν = I ν dΩ is the flux density, κ ν is the mass opacity in cm 2 g −1 of dust, D is the distance to the source, T is the temperature in Kelvin, and B ν is the black body radiation using the Planck function.We adopt the opacity of 0.023 cm 2 g −1 of gas from Beckwith et al. (1990) (see also recent evidence from Lin et al. 2021 in support of this prescription and Section 4.1) and assume a dust-to-gas mass ratio of 0.01 to obtain κ ν = 2.3 cm 2 g −1 of dust.We assume T = 20 K which is a commonly adopted value for surveys (e.g. which is optimized at a radius of 50 au (Tobin et al. 2020).With L bol = 0.43L (Ohashi et al. in  The continuum image of IRAS 04302+2247.The white contour marks the 5σ level (see σ in Table 1).The white ellipse in the lower right corner is the beam size (see Table 1) and the length scale is 50 au.The black plus sign marks the center from the best-fit 2D Gaussian.
limit and likely a drastic underestimation given the near edge-on view.

Lines
Fig. 3 shows the integrated intensity images for 12 CO 2-1, 13 CO 2-1, C 18 O 2-1, H 2 CO 3 0,3 -2 0,2 , and SO 6 5 -5 4 with robust=0.5.We also show 13 CO with ro-bust=2 which captures more large-scale emission that is resolved out from robust=0.5.The different molecules trace different spatial scales of the edge-on disk and their images also differ from the continuum image.
In the direction parallel to the disk major axis, which is described by the impact parameter, the extents of C 18 O, H 2 CO, and 12 CO (robust=0.5)appear comparable to the continuum image, while SO clearly spans a smaller range in impact parameter.The 13 CO image with robust=0.5 is the most extended and even more so with robust=2.The large extent of 13 CO with ro-bust=0.5,which is on the order ∼ 4 (640 au) from the center, suggests a gas disk that is larger than the dust disk (see Section 5.2 for more detail).
In the direction parallel to the disk minor axis, which corresponds to the "vertical" direction of an edge-on disk, the more optically thin lines, C 18 O, H 2 CO, and SO, are more confined to regions just to the east and west of the continuum.This suggests that these molecules trace the disk surface and not the outer edges of the disk in the radial direction.In contrast, the more optically thick lines, 12 CO and 13 CO, are far more vertically extended.The emission traced by 13 CO with robust=0.5appears to reach ∼ 2 (320 au) in the vertical direction (see Section 3.3).
A common feature across all images is the lack of emission near the supposed disk midplane and the emission appears to form a V-shaped pattern to the north and south.The same feature was observed in previous lower angular resolution observations in C 17 O 2-1 and H 2 CO 3 1,2 -2 1,1 by van 't Hoff et al. (2020) and in 12 CO 2-1, CS 5-4, H 2 CO 3 1,2 -2 1,1 by Podio et al. (2020).The lack of emission is largely due to freeze-out, especially at larger impact parameters, and the V-shape is a natural result of the snow surface given the typical 2D temperature structure of an irradiated disk (e.g., Aikawa & Herbst 1999;Dutrey et al. 2017;van 't Hoff et al. 2018;Qi et al. 2019;Zhang et al. 2019;Flores et al. 2021).Indeed, an absorption feature due to CO ice is also detected for this source in the infrared (Aikawa et al. 2012).At smaller impact parameters where the continuum could be optically thick, it could also be due to dust extinction (see Section 5.1).Intriguingly, the 13 CO image with ro-bust=0.5(and, similarly, the case with robust=2) not only has the V-shape where the emission diverges, but the emission converges at large impact parameters beyond ∼ 3 (480 au), enclosing a dark cavity and resembling the shape of the number "8" overall.We can infer that in the midplane, 13 CO is frozen out at small radii, but reappears at large radii (see Section 3.3).We discuss the cause for the re-emergence in Section 5.2.
Apart from the major axis, there is a lack of emission directly along the minor axis of the disk for 12 CO, 13 CO, C 18 O, and H 2 CO (it is less clear for SO).This can be  The cuts along the major and minor axis of the continuum disk by interpolating the image.The cut along the major axis is in green and that along the minor axis is in orange.The origin is at the center of the fitted 2D Gaussian (04:33:16.5,+22:53:20.2).The bottom and top axes mark the offset from the origin along the cut in arcsec and au.The positive location for the major axis is along the northern part of the disk, while the positive location for the minor axis is along the eastern part of the disk.The left and right axes mark the intensity in mJy beam −1 and brightness temperature (using the full Planck function) in Kelvin respectively.The line segment to the upper left corner represents the length of the FWHM of the beam.The shaded region is the intensity below 3σ.Right panel: Zoom-in comparison between the minor axis cut (solid black line) and the beam (dotted black line).
explained by optical depth effects (see van 't Hoff et al. 2018 for a visualization) through the following.For a rotating disk seen edge-on, only the material with projected speeds near the systemic velocity, which is along the minor axis, can contribute.At those channels, we trace regions further from the disk where it is colder since there is more material along the line-of-sight (see also Section 3.3 for the channel maps).For optically thicker lines, like 12 CO and 13 CO, much of the emission can even be resolved out and the depression is even more pronounced.Another factor that can decrease the brightness at small impact parameters is beam dilution (Flores et al. 2021) though the high resolution images here are likely less susceptible.
Another common feature seen in Fig. 3 is that the lines are all brighter on the east side compared to the west side, except for SO which appears brighter on the west side.The brighter east side can be interpreted as an inclination effect for a disk with a two-dimensional temperature distribution (Dutrey et al. 2017;Flores et al. 2021).We can infer that the brighter eastern side is the far side of the disk and the western side is the near side.The orientation is consistent with the orientation inferred from the continuum.The opposite behavior of SO, however, is puzzling and it could be due to other reasons, like chemical effects, rather than inclination effects (Sakai et al. 2014).Intriguingly, SO 2 was also detected to be brighter on the same side as SO (Garufi et al. 2022).
Fig. 4 expands upon the C 18 O in Fig. 3 and compares the intensity integrated image, peak intensity image, and the peak velocity image.To distinguish the redshifted and blueshifted halves, the peak velocity images are shown relative to a systemic velocity, v sys , of 5.7 km s −1 (see Section 4.2 for the measurement of v sys ).The velocity gradient of C 18 O shows a clear signature of rotation at the disk surface.The other two optically thin tracers, H 2 CO and SO, are shown in Fig. 5 and Fig. 6.Similar to C 18 O, both tracers also follow the disk surface and show similar velocity features.
From Fig. 7 and Fig. 8, we also see the same blueshifted and redshifted halves for 12 CO and 13 CO near the dust continuum which are similar to the C 18 O, H 2 CO, and SO, but there are additional extensions that do not follow what is expected from rotation.Notably, the extension towards the east side of the 12 CO image at ∼ 3 from the center is blueshifted.The level of blueshift increases with increasing distance from the center which is consistent with a Hubble-type outflow(e.gArce et al. 2007).A blueshifted outflow to the east is also consistent with the orientation of the disk where the east side is the far side.In addition, 12 CO with robust=2 appears more extended than its robust=0.5version, but it is still less extended than 13 CO with ro-bust=2.This is likely because much of 12 CO remains resolved out even with robust=2.

Tracing the CO Snow Line and Snow Surface
The optically thinner tracer C 18 O probes the snow line and the immediate snow surface better.For the convenience of the discussion from this section, we define x as the impact parameter along the disk major axis, where positive x lies in the northern part of the major axis and y as the location along the disk minor axis where positive y lies along the blueshifted side of the jet axis to the east of the disk midplane.The origin, x = 0 and y = 0, corresponds to the center of the fitted 2D Gaussian from Sec. 3.1.
A key feature of the C 18 O is a region near the midplane that clearly lacks emission which can be attributed to freeze-out (e.g.Dutrey et al. 2017;van 't Hoff et al. 2020;Villenave et al. 2022).Given the fine resolution, we can trace the snow surface to ∼ 0.1 .We show selected channel images of C 18 O with robust=0.5 in Fig. 9 and focus only on the northern half of the disk, since the freeze-out zone appears symmetric to the southern half (see Fig. 4).To increase the signal-to-noise ratio and to limit the number of channel images, we averaged every 3 channels and the noise level used for the figure is correspondingly decreased by √ 3. To outline the snow surface, we give a simple prescription (motivated by a similar prescription in Lee et al. 2021): This describes an increasing surface that begins from the snow line R s to some transition radius R v after which the V-shaped snow surface begins to close and ends at R e where the gas re-emerges.H v is the height at R s .The closing of the snow surface at large impact parameters is less clear in C 18 O, but more obvious in 13 CO which we show in Fig. 10 and discuss later.We estimated the parameters to be R s = 0.8 (130 au), and R e = 2.8 (448 au) by eye.Although the disk is not perfectly edgeon, we assume a symmetric outline across y = 0 for simplicity.
The mid-velocity channels (∼ 2 − 4 km s −1 for the northern half) show the iconic V-shaped emission expected for a snow surface.Under the simple expectation that C 18 O should exist from the center to the snow line, the emission should appear from the disk center at high velocities and emerge away from the center with decreasing velocity until the emission begins to concentrate along the disk minor axis at velocities near v sys (see other edge-on sources with Keplerian rotation, e.g., Dutrey et al. 2017, Teague et al. 2020, Flores et al. 2021).Thus, the snow line R s is based on the maximum impact parameter with emission that exists between the east and west surfaces.Note that the location could be an upper limit due to contamination from finite beam averaging of the east and west snow surfaces.Nevertheless, the R s of 0.8 (130 au) appears consistent with previous constraints using C 17 O (van 't Hoff et al. 2020).
While the optically thinner C 18 O probes the snow line and the immediate snow surface at small impact parameters, 13 CO reveals a complete freeze-out zone explained by the following.Fig. 10 shows selected channels for the northern half of the 13 CO with robust=2 for better detection of the larger scale structure.The first three high velocity channels also show the distinct V-shaped snow surface that extends further than C 18 O.At the low velocity channels, the east side and the west side of the emission appear to connect at large impact parameters (∼ 3 ) forming an apparent "cap" to the V-shaped Lin et al.  emission that closes the opening.The difference between C 18 O and 13 CO is likely due to optical depth and sensitivity.Since 13 CO can be detected more easily at lower column densities, we can identify the full spatial extent (or complete) freeze-out zone of CO, while C 18 O can only reveal the partial freeze-out zone.From the cap, it appears that CO is no longer frozen-out on grains at the larger radii even though one may expect that the temperature is lower than the inner radii.We extend this discussion in Section 5.2.
Another intriguing feature is the non-Keplerian, blueshifted feature in the southeast atmosphere of the disk.Fig. 11 shows selected 13 CO channel map with robust=2 focusing on the southern half.At redshifted channels (bottom row of Fig. 11), the disk near the midplane shows the typical Keplerian rotation for an edge-on disk (like that of C 18 O in Fig. 9, but for the blueshifted half).In addition, we see the freeze-out zone and the outer cap that is similar to the northern half (Fig. 10).However, at blueshifted channels (top row of Fig. 11), the southern half is not void of emission as expected from Keplerian rotation, but hosts large extensions to the east.The extension is along the disk minor axis at the most blueshifted channel of Fig. 11 and extends to the south when closer to the system velocity.Intriguingly, the edge of the extension closest to the disk appears to match the eastern edge of the redshifted Keplerian part in shape (the emission in the bottom row of Fig. 11).Thus, it appears that the blue extension is aware of the atmosphere of the southeast part of the disk and forms an interface.With a simple modification to Eq. ( 8), we outline the interface by: where H 0 is the distance from the center along the minor axis, H a is the height at some transition radius R a , and R c is the outer radius of the cap.We find that H 0 = 0.5 (80 au), H a = 0.9 (140 au), R a = 2.2 (350 au), and R c = 3.9 (620 au) by eye.We show the outline symmetric across x = 0 and y = 0 for convenience of discussion.
From the difference in the kinematics within and outside the southwest interface, we can distinguish the disk component and the envelope component.The southeast blue extension outside the interface is connected to even larger distances at ∼ 10 shown in the 13 CO with robust=2 (bottom row of Fig. 8).Given that the extension is closest to the disk at high blueshifted channels (e.g., ∼ 4.7 km s −1 ) and becomes more extended at lower blueshifted channels (e.g., ∼ 5.4 km s −1 ), the nature of the extension can be explained by infalling material from behind the plane-of-sky that lands onto the southern half of the disk.In such a scenario, it would make sense for an interface to form, since the infalling material moving from behind the plane-of-sky has to collide with the disk material moving into the plane-of-sky.The existence of infalling material is not too surprising given evidence in other younger Class 0/I disks (Pineda et al. 2020;Alves et al. 2020;Valdivia-Mena et al. 2022;Garufi et al. 2022) or even late-stage infall onto Class II disks (e.g., Tang et al. 2012;Ginski et al. 2021;Huang et al. 2020Huang et al. , 2021;;Gupta et al. 2023; see also Kuffmeier et al. 2020).
Though the outline of the interface was determined from the kinematic difference in the southeast part of the disk, the southwest part of the outline appears to also separate the disk from a broad extension to the west.Different from the southeast blue extension, the west extension is redshifted as the southern part of the Keplerian disk should be, making it indistinguishable kinematically and thus the outline from Eq. ( 9) may not mark a clear interface.However, given that there is no symmetric counterpart across y = 0 on the east side  at the same channels, it is morphologically distinct from the material within the outline.It is unclear what the nature of the west red extension is.

ANALYSIS
In this section, we analyze the data presented in Section 3 in more detail.Section 4.1 analyzes the continuum image through forward ray-tracing of the dust and provides constraints on the dust scale height and inclination.Section 4.2 analyzes the position-velocity (PV) diagram of C 18 O and measures the stellar mass.

Continuum Forward Ray-Tracing
Although a 2D Gaussian captures the overall features, such as the position angle and the overall shape, certain deviations stand out.Fig. 12a shows the original contin- The black cross marks the center of the image.The black solid contour is the snow surface (from Eq. ( 8)), and the black dashed contour represents the "interface" (from Eq. ( 9); see Section 3.3 for more detail).The black ellipse to the lower right corner is the beam.
uum and the fitted 2D Gaussian, while Fig. 12b shows the residuals, which are defined as the observed image subtracted by the 2D Gaussian.The largest deviation is the significant positive residual extending parallel to the disk major axis that is slightly offset from the center to the east.This corresponds to the asymmetry along the minor axis where the east side is brighter.
In this section, we demonstrate that the asymmetry along the minor axis is due to the inclination effect of an optically thick disk.We use a parameterized disk model and use RADMC-3D2 to conduct the ray-tracing (Dullemond et al. 2012).We refrain from conducting the heating/cooling calculations from RADMC-3D given the large computational cost and complexities regarding the dust opacity spectrum (e.g.Birnstiel et al. 2018).The calculation is beyond the scope of this paper and we leave it to a future paper.The parameterized disk model is a similar version of the disk model from Lin et al. (2021) which is suited for a disk viewed near edgeon. 3The model was applied to a Class 0 edge-on source, HH 212 mms, and successfully reproduced the asymmetry along the minor axis across ALMA Bands 3, 6, and 7.In the following, we briefly describe the key parts of the model and include modifications.
We parameterize the disk using the Toomre Q parameter (Toomre 1964) where c s is the isothermal velocity, Ω k ≡ GM * /R 3 is the Keplerian frequency, R is the cylindrical radius, and Σ is the gas surface density.For a gravitationally stable disk, Q must be greater than a value of order unity (e.g.Kratter & Lodato 2016).The pressure scale height of the gas is From basic arguments of vertical hydrostatic equilibrium, the gas density in the midplane is ρ g,mid = Σ/ √ 2π/H g and thus, when combined with Eq. ( 10), we have where R 0 is a characteristic radius, which we take to be the outer radius of the disk.For illustrative purposes, we assume that Q is a constant in the disk, and introduce a characteristic density ρ 0 ≡ M * /(π √ 2πR 3 0 Q), which is the density at the disk outer edge.
Since the dust disk appears vertically thin, we approximate the temperature with just a vertically isothermal prescription: where T 0 is the temperature at the outer edge of the disk and q specifies the temperature gradient.Note that the whole gas disk should have a vertical temperature gradient (warmer temperature in the atmosphere), which is needed for the existence of the clear snow surface (Fig. 3, 10).However, since the bulk of the dust disk appears to lie below the snow surface and there is no continuum dark lane (such as that found for HH 212 mms), the effect of a vertical temperature gradient is likely marginal, and thus, we only use a vertically isothermal profile for the dust disk.
As a further simplification, we fix q = 0.5 which is expected from passively irradiated disks in radiative equilibrium (e.g.Chiang & Goldreich 1997;D'Alessio et al. 1998).This assumption may not be entirely applicable to embedded protostars, which can have additional accretion heating or warming from the envelope (e.g.Butner et al. 1994;Agurto-Gangas et al. 2019).Accretion heating should lead to a steeper temperature gradient, usually q = 0.75 (Armitage 2015), and dominate the inner regions of the disk (Takakuwa et al. in prep.).Envelope warming prevails in the outer regions and should make the temperature gradient shallow (e.g., q ≤ 0.4 from Whitney et al. 2003).The Class I designation of IRAS 04302 motivates a smaller q, however, van 't Hoff et al. ( 2020) found q = 0.75 based on the location of snow lines of H 2 CO and C 17 O though the resolution is not ideal.We also refrain from fitting q directly, since a single wavelength image of an edge-on disk probes a limited range in radius due to the high optical depth (Lin et al. 2021).Longer wavelength observations are necessary to probe the temperature of the inner regions and using multiwavelength observations that probe different radii will better constrain q.Thus, given the uncertainties, we fix q = 0.5 as a compromise for this paper and leave the exploration of q for a future study. 4in et al. ( 2021) assumed that the dust and the gas are well-coupled and thus the dust also follows the gas in hydrostatic equilibrium (qualitatively, this means the dust scale height is equal to the gas scale height if the disk is vertically isothermal).However, to directly explore the dust scale height independent of what the gas Lin et al.
scale height should be, we parameterize the dust scale height by where the power-law index is the same as that from the gas scale height, i.e., 1.5 − q/2.Eq. ( 14) allows us to easily explore the effects of height with one parameter H 100 .
By assuming that the midplane density of the dust is related to the midplane density of the gas (Eq.( 12)) through a dust-to-gas mass ratio η, the complete dust density as a function of radius and height is where z is the vertical height and ρ d,0 ≡ ρ 0 /η is the midplane dust density at the outer edge of the disk.Instead of prescribing the dust opacity κ ν (in units of cm 2 per gram of dust) explicitly, we use the characteristic optical depth τ 0,ν defined as The definition makes sense because the characteristic length scale along the line-of-sight for an edge-on disk is R 0 .This parameter reflects the fact that opacity and density are degenerate and it is the optical depth (proportional to the product of opacity and density) that controls how an image appears (see Lin et al. 2021 for detailed derivation and for exploration of how τ 0,ν controls the image of an edge-on disk).In other words, τ 0,ν is a free parameter that we can fit from the image.As an initial exploration for this paper, we conduct the parameter search by hand.To limit the parameter space, we fix the position angle to 174.77 • obtained from the 2D Gaussian fit.The free parameters include τ 0,ν , T 0 , i, R 0 and H 100 in addition to the location of the star (δ RA , δ DEC ).The parameters for the best-fit model are listed in Table 2. Fig. 12c shows that the model compares quite well with the observations.The dust model can easily reproduce the shift along the minor axis towards the far side of the disk (towards the east for the case of IRAS 04302) since the disk is optically thick and highly inclined (Villenave et al. 2020;Takakuwa et al. in prep.).The residuals are shown in Fig. 12d and are evidently much lower than that from the simple 2D Gaussian fit (Fig. 12b).
We find that the H 100 is 6 au.The dust scale height from past modeling efforts based on lower resolution mm-images varies in the literature and ranges from ∼ 2 au to 15 au at a radius of 100 au (Wolf et al. 2003(Wolf et al. , 2008;;Gräfe et al. 2013;Sheehan & Eisner 2017) though it depends on the exact prescription of each model.By resolving the asymmetry along the disk minor axis, the new high-resolution image presented here offers a strong constraint on the dust scale height.In addition, the value is consistent with an independent study that modeled another high-resolution image at Band 4 (Villenave et al. 2023).On the other hand, the derived radius of R 0 = 310 au is consistent with past modeling efforts based on lower resolution mm-images in which case the major axis of the disk was well resolved (Wolf et al. 2003;Gräfe et al. 2013).
The inferred inclination of i = 87 • provides the necessary deviation from being perfectly edge-on (i = 90 • ) which would not produce an asymmetry along the minor axis since both halves across the midplane would be perfectly symmetric (e.g.Wolf et al. 2003).The value is also consistent with the lower limit of ∼ 84 • assuming the disk is completely flat (see Section 3.1).It is not surprising that the actual inclination is larger than the inclination inferred just from the ratio between the minor and major axes, or arccos(minor/major).Using the ratio assumes that only the radial extent contributes to the projected length along the minor axis which is indeed the case for a geometrically thin disk.However, for a highly inclined geometrically thick disk, the vertical thickness contributes to the projected width along the minor axis which decreases arccos(minor/major).
The inferred T 0 of 7.5 K appears to be lower than necessary when compared to what is expected from the estimated snow line of CO.The low temperature profile is necessary because the peak brightness temperature is only ∼ 14 K and yet the disk has to be optically thick to produce the minor axis shift of the continuum.Based on the fitted T 0 , the snow line for CO, assuming a freezeout temperature of 20 K, should be at ∼ 44 au (0.275 ).However, this appears inconsistent with the observed location of the snow line which is ∼ 130 au (∼ 0.8 ) from C 18 O (also similar to what was derived in van 't Hoff et al. 2018 from lower angular resolution observations of C 17 O).One possibility is that the dust temperature profile is correct and the observed C 18 O emission beyond the inferred snow line location (of 44 au) is contaminated by emission from the warmer surface layers due to the finite beam.
Another possibility to alleviate the above discrepancy is through scattering.Scattering makes objects appear dimmer, which means the actual temperature should be higher than what is inferred when assuming no scattering (e.g.Birnstiel et al. 2018).Interestingly, radiation transfer calculations for this source including scattering of 100 µm grains infer a temperature of 20 K at 100 au (Gräfe et al. 2013) which is higher than the 13 K at 100 au based on the model prescribed here.Given that scattering only scales the image intensities and does not alter the relative shape of the image much (Lin et al. 2021), the inferred low temperature could be evidence of scattering, but we leave the incorporation of scattering to a future study.
Intriguingly, the outermost contour of the model appears systematically less extended than the observations along the minor axis (Fig. 12c).This is also seen as two lanes of generally positive residuals to the east and west of the disk in Fig. 12d which suggests a more extended upper layer.However, increasing H c to broaden the image along the minor axis leads to even broader widths at the endpoints of the major axis of the disk.Thus, it appears that the dust scale height should not be too flared at the outer radius compared to the inner radius.This is in fact what we would expect from dust settling of a given grain size, where the outer region should be more settled than the inner region because the Stokes number of the grains increases as the density decreases towards larger radii (Dullemond & Dominik 2004).We leave also this possibility for future exploration.
We found that the characteristic optical depth is τ 0,ν = 0.35 which can be related to the opacity.5 From Eq. ( 16) and the definition of ρ 0 from Eq. ( 12), we can explicitly solve for κ ν through Using the best-fit R 0 and τ 0,ν from this section, the M * derived based on the rotation curve of C 18 O (see Section 4.2), the opacity is κ ν = 0.019Qη in units of cm 2 g −1 of gas.If the disk is gravitationally stable, Q should be greater than of order unity.Otherwise, the disk should fragment (Kratter & Lodato 2016).Thus, taking Q = 1 gives a lower limit to κ ν .We note that the lower limit to the opacity is per mass of gas since it is the gas that contributes most of the mass, and that limits the amount of material.However, theoretical dust models calculate dust opacity with respect to the mass of the dust (e.g.Ossenkopf & Henning 1994) and thus we have to assume a η to directly compare the dust opacity calculations to the observationally constrained opacity presented here.By assuming the standard η = 100, we get κ ν = 1.9Q cm 2 g −1 of dust.The uncertainty of κ ν is 0.5 cm 2 g −1 of dust based on error propagation from the uncertainty of M * derived in Section 4.2.We add the caveat that the opacity can vary spatially which is not captured through the model and thus, the value measured here is an effective opacity of the region observable at Band 6.The conventional Beckwith et al. (1990) opacity at λ=1.3 mm is κ ν = 2.3 cm 2 g −1 of dust (also constrained observationally and assumed η = 100) and the opacity based on HH 212 mms is κ ν = 1.33 cm 2 g −1 of dust (Lin et al. 2021).By taking Q = 1, it appears that the lower limit from IRAS 04302 lies right in between the two previous studies as shown in Fig. 13.For completeness, we have included opacity constraints at other wavelengths for HH 212 mms (Lin et al. 2021) and also another commonly adopted dust opacity model from Ossenkopf & Henning (1994) with calculations adopted for low and high densities.The lower limit from IRAS 04302 disfavors the opacity model from Ossenkopf & Henning (1994) and is more consistent with the Beckwith et al. (1990) prescription.
The proximity of the lower limit from IRAS 04302 to the opacity from HH 212 mms is intriguing, given that HH 212 mms is vastly different compared to IRAS 04302 Lin et al. in class, size of the disk, and stellar mass.While HH 212 mms is likely to be marginally gravitationally unstable given the small stellar mass, bright continuum, and early stage (Tobin et al. 2020), IRAS 04302, as a Class I source, is less certain.Even if grains have a universal opacity, the lower limit from IRAS 04302 need not be similar, since from Eq. ( 17), taking Q = 1 is only a lower limit after all and Q can take on any value greater than 1 if the disk is not marginally gravitationally unstable.
If not purely coincidental, a possible physical explanation is that the grains could be similar between these two systems and both systems are marginally gravitationally unstable which fixes Q to a value of order unity (e.g.Lodato 2007;Kratter & Lodato 2016;Xu & Kunz 2021).It may not be too surprising if IRAS 04302 can also be marginally gravitationally unstable given the large disk, an available reservoir of envelope material, and cold midplane temperature.There is growing evidence of other Class 0/I sources that are marginally graviationally unstable (e.g.Kwon et al. 2011;Tobin et al. 2020;Xu 2022).Furthermore, from an evolutionary standpoint, this is in line with evidence of Class II sources with Q that largely falls within 1 to 10 (e.g.Kwon et al. 2015;Cleeves et al. 2016;Booth et al. 2019;Veronesi et al. 2021;Paneque-Carreño et al. 2021;Ueda et al. 2022;Schwarz et al. 2021;Sierra et al. 2021;Yoshida et al. 2022;Lodato et al. 2022).

Deriving the Stellar Mass from Disk Rotation
Given the clear evidence of rotation (e.g., right panel of Fig. 4), we further analyze the rotation curve for IRAS 04302 using the position-velocity (PV) diagram along the major axis of the disk.Fig. 14 shows the PV diagram for C 18 O with robust=0.5.We choose C 18 O since it is optically thinner and only traces the disk as opposed to 12 CO and 13 CO which are more susceptible to surrounding envelope material.Also, C 18 O is better detected than the other optically thin lines.To create the PV diagram, we use the position angle derived from the Gaussian fit of the continuum (see Section 3.1) and we use a slit with a width of ∼ 2 beams to increase the signal-to-noise necessary for the analysis below.
We use the Spectral Line Analysis/Modeling (SLAM) code (Aso & Sai 2023) 6 to extract the rotation curve from the PV diagram (Aso et al. 2015;Sai et al. 2020).
Inferring the rotation properties relies on first assigning pairs of radius and velocity points based on the PV diagram and later fitting the points to a rotation Figure 13.The lower limit to the dust opacity (absorption cross section per gram of dust assuming a dust-to-gas ratio of 0.01) inferred from the IRAS 04302 disk (marked as an orange cross) in comparison to other millimeter dust opacities from the literature.The error bar is the uncertainty associated with the uncertainty from the stellar mass.The filled circles with solid lines are the lower limit to the dust opacity for the HH 212 mms disk from Lin et al. (2021).The corresponding lighter, shaded region is the uncertainty associated with the stellar mass and also the Toomre Q parameter (ranging 1 to 2.5) and the darker, shaded region is the uncertainty from the noise.The open circle is the Beckwith et al. (1990) opacity at 1.3 mm and its line segment represents the opacity index of 1.The open squares are opacities from Ossenkopf & Henning (1994) at 1 and 1.3 mm.curve.Details of SLAM are described in Ohashi et al. (in prep.), but we describe the essential steps and parameters adopted here.
For the first step, we aim to trace the "outer" edge of the PV diagram (the top of the second quadrant and the bottom of the fourth quadrant).We use the 5σ level for each spectrum along the position as the representative pairs of radius and velocity, which corresponds to the "edge" method in SLAM.The reason is as follows.For an edge-on disk, the line-of-sight at a particular impact parameter x (i.e., the position along the major axis) crosses several radii.The spectra is simply the collective emission of material along that line-of-sight each with varying levels of projected velocities (without considering any complications from finite line width).Along the line-of-sight, there is a minimal radius that contributes the maximal velocity and that is the location in planeof-sky which equals the impact parameter x (see e.g., Dutrey et al. 2017 for an illustration).Thus, in the spectra, we would expect that the maximum velocity where we have detection is precisely the representative velocity for the radius that equals (absolute value of) the impact parameter.Complications arise when considering finite line width, temperature effects, inclination, and detection levels, which can be addressed through modeling.However, as a working expectation, we use the 5σ level for each spectrum along the position to fit for the Keplerian rotation (Seifried et al. 2016).To assess how sensitive the parameters are to the chosen level, we also use the 3σ level.
Another common way to extract representative pairs of radius and velocity from the PV diagram is to take the mean of the intensity profile, which corresponds to the "ridge" method in SLAM (Aso et al. 2015;Yen et al. 2017;Sai et al. 2020).This extraction usually underestimates the true stellar mass (e.g.Maret et al. 2020), but we use it to complement the "edge" method described above to assess the systematic uncertainty.Conventionally, there are two ways to take the mean of the intensity profile, either along the velocity axis (i.e., the spectra at a certain impact parameter) or along the position axis (i.e., the profile of the image of a certain channel).From experimentation, we find that using both was necessary to trace the PV diagram.
The noise level used here is assessed in regions of the PV diagram where no emission is expected.We have σ = 0.976 mJy beam −1 which is ∼ √ 2 smaller than the channel map noise level (Table 1) as expected from our adopted slit width.In addition, we avoid fitting the spectra within 0.5 (80 au), since the PV diagram is even qualitatively different from the typical Keplerian rotation curve.The lack of high-velocity emission could be due to the lack of material at inner radii (Dutrey et al. 2017) or dust extinction, but we leave the verification for future exploration and focus on fitting the Keplerian parts in practice.
The next step involves fitting a rotation curve to the data points which we use where r b is a characteristic radius, v b is the rotational velocity at r b , and p is the power-law index of the rotation profile.The sign of v has been adjusted to account for the definition of x in this paper (positive along the northern part of the major axis which is consistent throughout the paper).If the disk is in Keplerian rotation, we should retrieve p = 0.5 and one can infer the stellar mass M * from v b = GM * /r b sin i where i is the inclination.The left panel of Fig. 14 shows the assigned pairs of position (radius) and velocity from the edge method using the 5σ level of the spectra (which corresponds to the intensity in the vertical direction of Fig. 14) at the largest (absolute) velocity with respect to v sys .The inferred rotation curve (plotted as a white curve) follows the outer edge of the PV diagram reasonably well.Considering only statistical uncertainty, we find that p = 0.52 ± 0.02 which verifies that the disk is consistent with Keplerian rotation.The systemic velocity v sys is 5.72 ± 0.02 km s −1 and the stellar mass is M * = 1.65 ± 0.02 M assuming an inclination of i = 87 • derived from the dust continuum (see Section 4.1).
To assess the systematic uncertainty, we compare the edge method with 5σ to the edge method with 3σ and the ridge method.The edge method using data points at the 3σ level (not shown in Fig. 14 for brevity) yielded a stellar mass of M * ∼ 2M (see Table 3 for a comparison of the results).Measuring a larger stellar mass is not too surprising, since adopting a lower threshold adds to the range of the measured velocity which could be due to the line width and it can artificially increase the measured stellar mass.In the other extreme, the extracted points for the ridge method (right panel of Fig. 14) yielded a smaller stellar mass of M * ∼ 1.2M * as expected (e.g.Maret et al. 2020).Given the large spread in measurements depending on the assumed method, we adopt M * = 1.6 ± 0.4 M and v sys = 5.7 ± 0.1 km s −1 .
The measured stellar mass is similar to the adopted mass of 1.7 M from Gräfe et al. (2013) although the detailed derivation was not described in the literature.Otherwise, as far as we know, there are no other published measurements of the stellar mass through dynamical measurements.

DISCUSSION
Figure 14.For both plots, the color scale is the pv diagram of C 18 O and the black contour shows the 5σ level.The left plot includes the fitted points from the "edge" method along the velocity axis, while the right plot includes the fitted points from the "ridge" method along the velocity axis (circles) and position axis (crosses).The best-fit rotation curve, using M * = 1.6M , for each is in solid white line and the best-fit vsys is shown as a horizontal dashed white line.The shaded region is the range of the rotation curve using M * = 1.6 ± 0.4M .The bottom and top axes are the location along the major axis in units of arcseconds and au respectively.The left and right axes show the velocity and that relative to vsys.

Evidence for Dust Extinction
As described in Section 3.2, the lack of line emission along the major axis of the disk at large impact parameters is due to freeze-out which gives the iconic V-shaped emission (van 't Hoff et al. 2020).However, at small impact parameters where we do not expect freeze-out, a depression is shared across all lines and is especially obvious from the moment 0 images of H 2 CO and SO (Fig. 3, 5, 6).
The depression along the innermost parts of the major axis can be explained by dust extinction.The lack of emission for SO and H 2 CO is likely because their snow lines lie well within the τ = 1 surface of the dust, i.e., the location where the optical depth to the observer is 1.In such a scenario, the dust essentially buries the emission behind the τ = 1 surface.For example, from the dust model shown in Section 4.1, the impact parameters where the total optical depth equals 1 and 5 are ∼ 215 au and ∼ 105 au, respectively.H 2 CO has a larger freezeout temperature at ∼ 70 K (Noble et al. 2012) and we would expect the snowline to be at impact parameters much less than 105 au well into the optically thick re-gions of the dust disk.On the other hand, the observed CO snowline from imaging of ∼ 130 au is roughly in the translucent region between the two limits which makes it possible to see the emission from the midplane.
Dust extinction can also explain the asymmetry of the high velocity emission from 12 CO.We show the channel maps focused on the innermost region of the disk in Fig. 15.At high blueshifted channels (top row of Fig. 15), 12 CO initially emerges as a single point in the east side (v = −2.85km s −1 ) until a second point appears on the west side (v = −2.21km s − 1 ).The same is true for the high redshifted channels (bottom row of Fig. 15) where the single point in the east side at the highest velocity channel (v = 14.30km s −1 ) and the second point in the west appears at a lower redshifted channel (v = 13.03km s −1 ).A similar behavior is evident for C 18 O in Fig. 9.With the nearly edge-on view, the emission from the front side of the disk can be seen unobstructed, while the emission from the back side of the disk must travel through the dust disk to reach the observer.
The existence of dust extinction is consistent with the requirement that the dust must be optically thick to produce the continuum asymmetry along the minor axis (see Section 4.1; unless it is due to an intrinsic asymmetry in the density distribution of the disk).It is also not too surprising as other sources also have examples of dust extinction, for example, the rings of HD 163296 (Isella et al. 2018) and DG Tau B (Garufi et al. 2020).
We note that an inclined disk with a two-dimensional temperature structure of a warmer surface and colder midplane could also contribute to the asymmetry in the brightness between the near-and far-sides (e.g.Flores et al. 2021).The brightness asymmetry further away from the major axis of the disk is more likely from the inclination effect.A complete radiation transfer including both dust and gas would be required to identify the separate contributions.

The outer cap of 13 CO
An intriguing part of the 13 CO morphology is the detection at ∼ 4 from the center along the disk major axis even though the molecule is not seen from ∼ 1 to 3 along the disk major axis.We interpret the lack of emission due to freeze-out and the transition to the freeze-out zone extends into the atmosphere resembling the shape of "V."However, we detect 13 CO at ∼ 4 in the form of a cap that closes off the freeze-out zone which means the molecule is no longer frozen-out and somehow "reemerges" at larger radii where the temperature is usually expected to be lower.The cap also exists for 12 CO, but is only visible when robust=2.0 (Fig. 7 bottom row), which suggests that much of the emission along with the cap is mostly resolved out for 12 CO.The C 18 O cap is not evident with robust=0.5 and also not clear with robust=2.0which could be due to the lack of sufficient signal-to-noise for the optically thinner isotopologue.
A natural question is whether the cap belongs to the disk or envelope.Fig. 16 shows the PV diagram along the major axis of 13 CO with robust=2.0with a slit width of ∼ 1 beam.For comparison, we show the Keplerian curve with M * = 1.6M and v sys measured from C 18 O in Section 4.2.At impact parameters within ∼ ±1.5 (240 au), the Keplerian curve follows the outer extent of the PV diagram quite well which is similar to the case of C 18 O (Fig. 14).Within ±[1.5 , 2.8 ], there is a lack of material that follows Keplerian rotation which corresponds to the freeze-out zone.The snow line from C 18 O is 0.8 (Section 3.3) which is less than the inner boundary of 1.5 here due to significant contamination from the warmer surface from larger beam averaging.The cap begins at ∼ ±2.8 and appears to follow the Keplerian rotation curve up to ∼ ±3.9 (620 au).For the southern part (negative x-axis of Fig. 16), the emission stops and we can directly identify the same edge of the emission in the channel maps in Fig. 11.For the norther part (positive x-axis), there appears to be a sharp break in the PV diagram in which case much of the emission appears more redshifted than Keplerian.Given the consistency with Keplerian rotation, we reason that the cap belongs to the Keplerian rotating disk and the sharp deviations from Keplerian rotation at ∼ 3.9 correspond to the edge of the gas disk outside of which is a part of the envelope.
There are two other edge-on sources with an apparent cap, namely 2MASS J16281370-2431391 (so called "Flying Saucer"; Dutrey et al. 2017) and SSTC2D J163131.2-242627(Oph 163131 for short; Villenave et al. 2022) which are both Class II sources.Dutrey et al. (2017) showed that 12 CO also emerges beyond the freeze-out zone.They found that the transition coincided with a change in grain properties and proposed that the behavior was expected if an efficient rise of UV penetration was re-heating the disk.
The disk around Oph 163131 also showed a 12 CO cap (Flores et al. 2021;Villenave et al. 2022).Since the transition from the inner "V"-shaped region to the outer cap region roughly coincided where the millimetercontinuum disk ends and the disk's scattered light stops, Flores et al. (2021) also interpreted the behavior as external UV radiation providing an additional source of heating to the outer part of the disk where dust particles are not present.
In contrast to the two sources with 12 CO caps, the cap of IRAS 04302 is seen in the optically thinner 13 CO which suggests that there is much more material.Another difference is the location of the cap region.Unlike the two sources whose caps begin at the extent of their millimeter-continuum, the continuum disk of IRAS 04302 clearly ends (at ∼ 1.9 ; 310 au) well before 13 CO emerges (at ∼ 3 ; 500 au).In other words, the freeze-out zone extends beyond the millimeter-continuum disk.A potential explanation is that the freeze-out temperature could decrease in the low density region in the outer disk (Harsono et al. 2015).Another possibility is that external UV irradiation still impacts the disk and heats up the outer region, but the UV photons are efficiently blocked by the smaller grains which are invisible at mm-wavelengths.This scenario may suggest significant radial drift of the larger ∼ mm-grains which is not too surprising given the much smaller radius of the dust disk (310 au; see Section 4.1) compared to the radius of the gas disk (620 au).In fact, when modeling the scattered light and lower resolution mm-continuum simultaneously, Gräfe et al. 11.12 11.76 12.39 13.03 13.66 14.30 Figure 15.Selected channel maps of 12 CO in comparison to the continuum.The top row correspond to the high blueshifted channels and the bottom row are the high redshifted channels.The red color map is the 12 CO J=2-1 emission above 3σ (see Table 1), while the black contours mark the 5, 10, and 20σ levels.The underlying blue color map is the continuum, while the grey contour is the 5σ level of the continuum.The black ellipse in the lower right of each image represents the beam for the 12 CO.The text in the upper left corner of each plot is the velocity in km s −1 .
(2013) required one population of large grains with a smaller radius and another population of small grains with a larger radius.Another related possibility is also warming of the outer disk, but from the envelope (e.g.Whitney et al. 2003).

Dust Settling in the Class I stage
One of the most striking features of the IRAS 04302 disk is the shift of the intensity peak along the minor axis of the continuum image which is a tell-tale sign of dust with finite vertical extent, i.e., non-settled dust.This feature exists for several other sources among the eDisk sample, including CB 68 (Kido et al. in prep.), L1527 IRS (van 't Hoff et al. in prep.),IRS 7B (Ohashi et al. in prep.; Takakuwa et al. in prep.), , IRAS 32 (Encalada et al. in prep.),BHR 71 (Gavino et al. in prep.),IRAS 04169+2702 (Han et al. in prep.), and IRAS 16253-2429 (Aso et al. in prep.).
In one extreme, dust settled into an infinitely thin sheet should appear symmetric across the minor axis and for disks with rings, the rings and gaps should not show azimuthal variation (e.g.Pinte et al. 2016;Doi & Kataoka 2021).Several observations of Class II sources show that the dust is predominantly well settled (e.g.Andrews et al. 2018;Long et al. 2018;Villenave et al. 2020;Doi & Kataoka 2021;Liu et al. 2022;Villenave et al. 2023).One of the clearest case is SSTC2D J163131.2-242627(or Oph 163131 for short) whose gaps are resolved even though the disk is near edge-on (Villenave et al. 2022, i ∼ 84 • ).The inferred dust scale height is ≤ 0.5 au at 100 au, which is an order of magnitude smaller than that of IRAS 04302.Furthermore, the significant difference in the vertical extent of the gas and dust also shows that dust is decoupled from the gas over most of the disk volume away from the midplane (e.g.Villenave et al. 2020;Law et al. 2021Law et al. , 2022)).
In the other extreme, the Class 0 source, HH 212 mms, hosts a clear dark lane sandwiched between two bright lanes in the dust continuum at ∼ 1 mm, which is evidence that the dust is elevated high enough to trace the warm surface layers.The dust scale height is ∼ 12 au at a radius of ∼ 36 au and the dust was shown to follow the gas in hydrostatic equilibrium (Lee et al. 2017;Lin et al. 2021).
From Section 4.1, we found that the dust scale height is 6 au at a radius of 100 au.For comparison, the gas pressure scale height from Eq. ( 11) is H g = 5.8±0.7 au at a radius of 100 au after adopting M * = 1.6±0.4Mfrom Section 4.2 and the dust isothermal temperature profile of Eq. ( 13) with the best-fit T 0 (only the stellar mass uncertainty is included here).The effectively equivalent scale heights given the uncertainties suggest that the dust has not separated from the gas vertically.
We caution that there is ambiguity in the midplane temperature, since the temperature derived from dust modeling appears different from the temperature inferred from the freeze-out location of CO.Using the snow line of 130 au (see Section 3.3) and assuming a freeze-out temperature of 20 K with q = 0.5, the temperature at 100 au is ∼ 23 K and results in H g = 7.6 ± 1.0 au.Considering the ambiguity of the temperature profile from the two scenarios, we have ∼ 0.8 ≤ H d /H g ≤∼ 1.We also note that H d inferred from Section 4.1 assumes a mixed, single population of grains.However, if grain growth has occurred, we may expect grains of different sizes to settle at various characteristic heights (e.g.Dubrulle et al. 1995).Nevertheless, the inferred H d represents the characteristic height of the bulk of the material that is responsible for the λ = 1.3 mm emission which is already different from the Class II sources where the dust responsible for the emission at the same wavelength has already settled to a much smaller scale height as mentioned above.At face-value, the non-significant level of dust settling may pose difficulties for the streaming instability to produce planetesimals (Gole et al. 2020) and thus delay planet formation.
Although the dust traced by 1.3 mm continuum is non-settled, the dust in general appears very distinct from the distribution of gas molecules (demonstrated in Fig. 3) and also very distinct from the scattered light images of IRAS 04302.Fig. 17 shows a comparison between the 1.3 mm continuum, 13 CO, and scattered light images from the Hubble Space Telescope (HST) at 1.6 µm (Padgett et al. 1999).We describe the correction for proper motion in Appendix C. Strikingly, each image traces a spatially distinct location.The 1.3 mm continuum appears only near the midplane, while the scattered light only exists in the bipolar cavities. 13CO fills the atmospheric regions of the disk and reaches beyond the radial extent of the 1.3 mm continuum and scattered light.Nevertheless, the gas pressure scale height H g of 6 ∼ 7 au may not be too surprising, since the line emission can typically be at several pressure scale heights above the midplane (e.g.Dullemond & Dominik 2004;Wolff et al. 2021;Flores et al. 2021;Villenave et al. 2022;Law et al. 2022;Paneque-Carreño et al. 2022), and small dust grains are present in the bipolar nebula to scatter optical/IR light.Detailed modeling using the high-angular resolution observations of the molecular lines with the dust could give a more robust view on the level of dust settling.
Another distinction between the gas and mmcontinuum is the radial extent.The edge of the dust disk has a radius of ∼ 310 au (see Section 4.1), while the edge of the gas disk has a radius of ∼ 620 au (see Section 5.2).In light of the disparity in the dust and gas radii, but similarity in the dust and gas scale heights (see Section 4.1), IRAS 04302 demonstrates that radial settling occurs sooner than vertical settling.Nevertheless, proper forward ray-tracing including both the dust and gas will make the disparity more definitive.
IRAS 04302 is formally a Class I source based on the SED (Ohashi et al. in prep.).Although an object with the Class I designation could actually be a Class II source if viewed edge-on, there is additional evidence that IRAS 04302 is indeed younger than formal Class II sources.First, the scattered light image of IRAS 04302 is noticeably irregular which indicates potential interactions with its envelope.In contrast, scattered light images of Class II sources tend to be well-ordered (Villenave et al. 2020).Second, IRAS 04302 has clear evidence of extended 13 CO and 12 CO emission beyond the Keplerian disk surface with kinematics inconsistent with Keplerian rotation which is likely part of the envelope (Fig. 7 bottom row and 8 bottom row).Thus, it is quite clear that this Class I source is a case where there is rela- tively little dust settling amid infall and outflow.Given that most Class II sources appear settled, we speculate that substantial dust settling should happen between the Class I stage and Class II stage.
It is curious whether IRAS 04302 has any radial substructure given its Class I stage.Rings and gaps are ubiquitous around Class II protostars (e.g Andrews et al. 2018;Long et al. 2018) and these structures could be signposts of planets (e.g.Zhang et al. 2018).Gaps from a highly inclined disk like Oph 163131 were resolved (Villenave et al. 2022), but the order of magnitude larger dust scale height of IRAS 04302 can easily obscure the gaps if there exists any.

Stellar Mass-Luminosity Tension
From the rotation curve of C 18 O, we find that the stellar mass is ∼ 1.65M .Given the stellar mass and depending on the age of the protostar, we should expect a luminosity that is greater than ∼ 2L (e.g.Iben 1965;Kippenhahn & Weigert 1994;Hillenbrand & White 2004).However, the bolometric luminosity which is estimated to be ∼ 0.43L (Ohashi et al. in prep.) is much smaller.
The diminished level of luminosity is likely because of the edge-on view.Most of the stellar photons along the line-of-sight are removed by extinction and not replenished by scattering causing an underestimation of the total luminosity (Whitney et al. 2003).Indeed, Gräfe et al. (2013) relied on a larger input stellar luminosity of 5L to explain both the scattered light and mm-continuum of IRAS 04302.In contrast, younger protostars, like L1527 IRS, do not see such a difference between the bolometric and expected protostellar luminosity (Tobin et al. 2012).We speculate that the small envelope size, low density, and wide outflow cavities for IRAS 04302 results in lots of emission not being reprocessed and simply escape along the polar regions.The much larger envelope in size and mass of L1527 IRS (Tobin et al. 2008), on the other hand, could help capture and reprocess the photons.Whether the observationally inferred low bolometric luminosity of IRAS 04302 is consistent with the newly obtained stellar mass remains to be determined quantitatively.In principle, a fully consistent physical modeling including stellar irradiation can constrain the stellar luminosity since the disk temperature is constrained through imaging (e.g.Gräfe et al. 2013;Sheehan & Eisner 2017), but we leave it as a future effort.

CONCLUSION
As part of the ALMA large program, eDisk, we presented high resolution ALMA Band 6 dust continuum and line emission of the nearly edge-on Class I disk IRAS 04302.Our main results are as follows: 1.The dust continuum image has an angular resolution of ∼ 0.05 (∼ 8 au) and shows a nearly edgeon disk with a clear brightness asymmetry along the disk minor axis.By fitting the disk with a 2D Gaussian, we find that the lower limit to the inclination is ∼ 83.6 • using the ratio of the major and minor axis FWHM.Through forward ray-tracing of the dust, we find that the inclination is ∼ 87 • and that the disk needs to be optically thick and geometrically thick to produce minor axis asymmetry.There is no evidence of rings and gaps, which could be due to the lack of radial substructure or because the highly inclined and optically thick view obscures the gaps.
2. We detect 12 CO, 13 CO, C 18 O, H 2 CO, and SO and find that all five exhibit V-shaped integrated intensity images which can be explained by freezeout near the midplane.From C 18 O, we estimate by eye that the CO snow line is located at ∼ 0.8 (130 au).However, the frozen-out midplane only extends to ∼ 2.8 (450 au) after which we detect the optically thicker tracer 13 CO emission out to ∼ 3.9 (620 au) which forms a "cap" of emission closing the V-shaped opening and produces a welldefined "8"-shaped CO depletion region along the disk major axis (see Fig. 3).
3. By fitting the position-velocity diagram of C 18 O along the disk major axis, we find that the disk is in Keplerian rotation and that the stellar mass is 1.6 ± 0.4M (see Fig. 14).The mass is in tension with the low observationally inferred bolometric luminosity of ∼ 0.43L .
4. The optically thick lines, 12 CO and 13 CO, trace significant amounts of complex extended structures outside of the Keplerian rotating disk.We find 12 CO outflow along the blueshifted jet axis to the east which is consistent with the orientation of the continuum disk in which the far side of the disk is also to the east.In addition, we find blueshifted 13 CO emission next to the redshifted part of the Keplerian disk surface which we suggest as material infalling onto the disk (see Fig. 11).
5. Our most important conclusion is that the dust has yet to settle significantly in the Class I IRAS 04302 disk.We find a dust scale height ∼ 6 au at a radius of 100 au, which is comparable to the gas scale at the same radius.This result, coupled with the lack of dust settling in Class 0 disks, such as HH 212 mms, indicates that substantial dust settling should happen between the Class I stage and Class II stage.In addition, the radial extent of the dust disk is likely smaller than the gas disk which suggests that radial drift occurs sooner than vertical settling.tronomy Observatory is a facility of the National Sci-ence Foundation operated under cooperative agreement by Associated Universities, Inc.

APPENDIX
A. CONTINUUM Fig. A1 compares the continuum images for different robust weightings from -1, -0.5, 0, 0.5, and 1.The asymmetry along the disk minor axis is more evident with robust weightings smaller than 1, but robust weightings less than 0 begins to resolve out the large scale major axis.The robust = 0.5 is a good compromise between resolving the minor axis asymmetry and not resolving out the large scale major axis.
B. DUST MODELING UNCERTAINTY Section 4.1 presented model image with parameters that are able to capture most of the features of the data.In this section, we provide a simple verification that the adopted parameters are the local best-fit by varying each individual parameter.The simple approach assumes that the parameters are not too correlated, which can be verified with better parameter space sampling techniques (e.g.Foreman-Mackey et al. 2013).However, the complete exploration of the multi-dimensional parameter space is beyond the scope of this first-look paper.
For interferometric data, the visibility plane is where one should consider the goodness-of-fit between the model and observation to include the effects of finite sampling in the visibility plane and also it is where the true native uncertainty of each measurement resides.However, as an initial assessment and for easier comprehension of the image, we simply compare the model and observation in the image plane and leave the more complete post processing to a future exploration.To assess the goodness-of-fit, we calculate the reduced χ 2 , defined as χ2 , through where i iterates through each pixel of the image in steps of the beam size.O i and M i are the intensities the data and that of the model, respectively, at the ith pixel and N is the total number of the selected pixels.Considering the 7 free parameters (i, R 0 , T 0 , H 100 , τ 0,ν , δ RA , and δ DEC ), we create a series of models by varying each parameter.The range and step size for each parameter are listed in Table 4.
The right column of Fig. B1 shows the χ2 by varying each parameter around the best-fit values shown in Section 4.1.The intensities along the major and minor axes of the model with each varying parameter are also shown in the left and middle column of Fig. B1 as a comparison to the observation to identify the effects of each parameter.We note that while we only show the major and minor axes profiles, χ2 is evaluated across the image and not only along the major and minor axes.
Fig. B2 shows the χ2 in two-dimensions from varying both δ RA and δ DEC , since both simply describes a translation of the image in the plane-of-sky.

C. VLA BAND KA IMAGE AND PROPER MOTION
To account for proper motion, we utilize the Very Large Array (VLA) Ka-band observed in 2015 (PI: John Tobin; project code: 15A-381).Since the VLA Ka-band image has not been published elsewhere, we briefly describe the calibration procedure and resulting images below.
The VLA Ka-band data were observed in both B and A configuration, which have maximum baseline lengths of ∼11 km and ∼36 km, respectively.The B-configuration data were taken on 2015 Feb 16 with a 3 hour execution, and the A-configuration data were taken on 2015 Aug 16 and 2015 Sept 8 with ∼1.5 hour executions.The observations all used 3C84 as the bandpass calibrator, 3C147 as the flux density calibrator, and J0440+2728 as the complex gain calibrator.During the observations, pointing was updated approximately every hour using the source J0403+2600.The correlator was configured for 3-bit continuum mode with 4 GHz basebands centered at 28.97 GHz and 36.796GHz.The bandwith was broken up into 64 spectral windows, each with 64 channels and 128MHz in width.
The data were processed using the scripted VLA calibration pipeline (version 1.3.1) in CASA 4.2.2.We ran the pipeline twice, we used the first run to identify data that required flagging.We then applied the necessary flags to Lin et al. the data and re-ran the pipeline on the edited dataset.Then to prepare for imaging the data, we combined the three measurement sets into a single measurement set using the CASA task concat.
We created the image with the robust set to 0.5 and uvtaper set to 3000kλ.The noise level is σ =6.3 µJy beam −1 .The resolution is 0.13 × 0.13 with a beam position angle of −72 • .The representative frequency is 33 GHz (9.1 mm).Following Section 2, we assume a 10% absolute flux calibration uncertainty, but only consider the statistical uncertainty for the rest of this section.
Fig. C1 (left) shows a continuum image that is centrally peaked and largely elongated along the north and south.The southern part of the disk appears slightly broader and brighter than the northern part.Using imfit from CASA, we fit a 2D Gaussian and obtain a center of (04:33:16.4952, +22:53:20.34) in ICRS.The integrated flux from the fitted 2D Gaussian is 490 ± 30 µJy, while the integrated flux above 3σ is 380 µJy.The FWHM of the deconvolved major axis is 750 ± 50 mas and that of the minor axis is 150 ± 20 mas.The position angle is 175 • ± 1 • which is consistent with the fitted 2D Gaussian for Band 6.The consistent position angle is evidence that the Band Ka image still detects the edge-on dust disk, while the smaller major axis FWHM is expected as the disk is optically thinner at the longer wavelength (e.g.Lin et al. 2021).
Fig. C1 (right) also shows the cuts along the major and minor axes using the position angle derived from Section 3.1.The major axis does not appear symmetric from the origin.For example, the secondary peak at −0.25 corresponds to a trough at +0.25 with values of 51 µJy beam −1 and 21 µJy beam −1 respectively and the difference is ∼ 5σ.Similarly, the secondary peak at +0.35 corresponds to a trough at −0.35 .At face value, the secondary peaks could suggest the existence of substructure and rule out symmetric rings given the non-symmetric locations of the secondary peaks.However, as demonstrated in the case of L1527 IRS, secondary peaks may not correspond to physical substructure when observed with better integration times (Nakatani et al. 2020;Sheehan et al. 2022).
The peak of the continuum is 0.115 mJy beam −1 and, equivalently, the brightness temperature is 8.7 K using the full Planck function.The brightness temperature is much lower than the ∼ 14 K from Band 6 in Section 3.1, but slightly higher than ∼ 6.7 K at Band 4 (Villenave et al. 2020).Although one may expect that, for an edge-on disk, the longer wavelength should trace the inner regions with higher temperature and lead to higher brightness temperature (Lin et al. 2021), the low value is likely because the vertical extent of the disk is unresolved.
By comparing with the Band 6 continuum from this work, we find a proper motion of ∼ (7, −17) mas yr −1 in RA and Dec, respectively.The proper motion is similar to that of L1527 IRS which is in the same Taurus region (Loinard et al. 2002).

Figure 2 .
Figure2.Left panel: The cuts along the major and minor axis of the continuum disk by interpolating the image.The cut along the major axis is in green and that along the minor axis is in orange.The origin is at the center of the fitted 2D Gaussian (04:33:16.5,+22:53:20.2).The bottom and top axes mark the offset from the origin along the cut in arcsec and au.The positive location for the major axis is along the northern part of the disk, while the positive location for the minor axis is along the eastern part of the disk.The left and right axes mark the intensity in mJy beam −1 and brightness temperature (using the full Planck function) in Kelvin respectively.The line segment to the upper left corner represents the length of the FWHM of the beam.The shaded region is the intensity below 3σ.Right panel: Zoom-in comparison between the minor axis cut (solid black line) and the beam (dotted black line).

Figure 3 .
Figure 3.The integrated intensity images for 12 CO (top left), 13 CO with robust=0.5(top center), 13 CO with robust=2 (right), C 18 O (bottom center-left), H2CO 30,3-20,2 (bottom center-right), and SO (bottom right).The continuum is shown in the bottom left.The color scale starts from 0. The horizontal bar represents 100 au which is the same across all panels.The white ellipse in the lower right corner of each plot represents the beam size.

Figure 4 .
Figure 4.The integrated intensity image (left), the peak brightness temperature map (center), and the peak velocity image (right) for C 18 O with robust=0.5.The velocity image is plotted relative to vsys = 5.7 km s −1 (see Section 4.2).The grey contour outlines the 5σ level of the continuum.The horizontal bar represents the 100 au length scale.The ellipse to the lower right is the beam size.

Figure 6 .
Figure 6.Similar to Fig. 4 but for SO with robust=0.5.
Figure 7.Top row: Similar to Fig.4but for 12 CO with robust=0.5.Bottom row: The moment images for 12 CO with robust=2.0.Note that the image covers a larger region than the top row to show the larger scale structure.
Figure 8.Top row: Similar to Fig.4but for 13 CO with robust=0.5.Bottom row: The moment images for 13 CO with robust=2.0.Note that the image covers a larger region than the top row to show the larger scale structure.

Figure 9 .Figure 10 .
Figure 9. Selected channel images of C 18 O with robust=0.5.Each image here is first averaged by 3 channels to increase the signal-to-noise ratio.The red color map is the C 18 O emission and the black contours mark the 3 and 5σ levels.The underlying blue color map is the continuum, while the grey contour is the 5σ level of the continuum.The black ellipse in the lower right of each image represents the beam for C 18 O.The text in the upper left corner of each plot is the velocity in km s −1 .The inferred snow surface is outlined in black.The white cross marks the origin of the image.

Figure 11 .
Figure11.Selected channel maps of 13 CO (robust=2.0)showing emission above 3σ.The images in the top row are blueshifted and those in the bottom row are redshifted with respect to the system velocity.The grey contour is the 5σ level of the continuum.The black cross marks the center of the image.The black solid contour is the snow surface (from Eq. (8)), and the black dashed contour represents the "interface" (from Eq. (9); see Section 3.3 for more detail).The black ellipse to the lower right corner is the beam.

Figure 12 .
Figure 12.Comparisons between the observed continuum and models.The top row shows the 2D Gaussian model, while the bottom row shows the model with radiative transfer (see Section 4.1).The plots in the left column show the model (in blue contours) plotted against the observed continuum (in black contours).The color maps in the right column show the residuals (the observed subtracted by the model) in mJy beam −1 .The solid and dashed contours trace the 3σ and −3σ levels respectively.

Figure 16 .
Figure16.The colormap is the PV diagram of 13 CO (ro-bust=2.0).The black and grey contours show the 3σ and -3σ levels respectively.The white line is the Keplerian rotation curve and the white horizontal dashed line is vsys.The two white vertical line segments marks ±3.9 to denote the sharp deviation from Keplerian rotation (see Section 5.2 for more detail).The bottom and top axes are the location along the major axis in units of arcseconds and au respectively.The left and right axes show the velocity and that relative to vsys.

Figure A1 .
Figure A1.The continuum images for different robust weightings from -1, -0.5, 0, 0.5, and 1 (left to right).The color scale is the flux density in µJy beam −1 .The white contour shows the 5σ region.The white ellipse to the lower right of each image is the corresponding beam size.

Figure B1 .Figure B2 .
Figure B1.Left column: The comparison between the major axis of the observed continuum and that of the models with various parameters.The grey shaded region of the observed cut is the noise uncertainty.Middle column: The same comparison but along the minor axis.Right column: The χ2 distribution for each value of the considered parameter.The color for each data point corresponds to the color used for each model in the major and minor axes cuts in the left and right columns.The parameters from top to bottom are: i, R0, T0, Hc, and τ0,ν .

Table 1 .
Summary of ImagesNote-seeOhashi et al. (in prep.)for the complete spectral setup.
prep.), we have T ∼ 34 K and the dust mass is M dust ∼ 70M ⊕ .

Table 2 .
Adopted parameters for the dust model

Table 3 .
Results from SLAM

Table 4 .
The grid of parameters considered for the dust model The step size of the parameter.Column (7): The final adopted value which is consistent with Table2.