© The Authors 2023

1 Introduction

Life, as we know it, is largely formed from the elements C, H, O, N, and S (CHONS), which, in protoplanetary disks, are carried mostly as ices or organics on dust grains (Henning & Semenov 2013; Furuya et al. 2013; Krijt et al. 2023). These grains are accreted by newborn planets, either as they form or during a later delivery by planetesimals residing in warm or cold debris reservoirs (asteroid and Kuiper Belt analogs) in more mature systems. Such organic-rich icy grains are formed in molecular clouds and are funneled into the planet-forming regions of disks via an infalling protostellar envelope during the early stages of star formation (e.g., Visser et al. 2009; Hincelin et al. 2013; Drozdovskaya et al. 2016). During infall, some of the ices sublimate during episodic, protostellar mass accretion events (Kim et al. 2012), and they could dissociate to their atomic constituents in the gas phase before reforming and freezing out again in the disk, partly resetting the disk’s ice chemistry relative to the cloud chemistry (Ballering et al. 2021). The relative amount of different ice species inherited by the disk from the cloud versus reset during this process is observationally still poorly constrained (Öberg & Bergin 2021).

Within disks, the ice species are expected to be vertically and radially segregated by 3D “snowlines”, which is the point where species sublimate from the ice phase. Grain growth and gravitational settling produce a vertical dust distribution that dictates where the stellar radiation is absorbed and scattered, setting a strong vertical temperature gradient. A corresponding chemical gradient arises between the disk midplane, where planetesimals form from large, icy grains, and the disk’s upper layers, where UV radiation photodesorbs the ice from the grain surface at a vertical snowsurface (Dominik et al. 2005; Kamp et al. 2018). Similarly, a radial temperature gradient gives rise to snowlines at certain distances from the central star, where inherited or reformed ices sublimate. Disk snowlines gave rise to Solar System planets, moons, and small bodies with different compositions (Bergin et al. 2015) and may have also set the bulk atmospheric properties of the planets (Öberg et al. 2011b; Madhusudhan 2019). The distribution of major CHONS ices during the initial stages of planet formation in protoplanetary disks can be inferred from comets, asteroids, and icy moons in the Solar System (e.g., with Rosetta, Hayabusa 2, OSIRIS-REx, Juice) or in situ snapshots of ices at different stages of planet formation. The latter, in particular, is critical in order to understand the variety of processes and planetary outcomes during planet formation (Krijt et al. 2023).

Individual ice species are detectable in the laboratory through fundamental vibrational modes seen in the infrared between 1 and 140 µm (e.g., Rocha et al. 2022, and references therein). These laboratory spectra can be used to identify astro-physical ices in different environments. Ices have been detected in absorption against the hot continuum provided by stars located behind cold, dense molecular clouds in protostellar envelopes, in the cold outer regions of edge-on protoplanetary disks (Boogert et al. 2015), and in emission from the midplanes of a select number of disks around T Tauri and Herbig Ae/Be stars from 44–63 µm with the Infrared Space Observatory (ISO) and the Herschel Space Observatory (Chiang et al. 2001; McClure et al. 2012, 2015; Min et al. 2016). Ice emission requires full radiative transfer modeling to be fully interpreted; however, for pure absorption spectroscopy, it has been standard practice to calculate logarithmic ice optical depths from spectroscopic data using the following: $$\tau =-\text{ln\hspace{0.17em}}\left(\frac{F_v}{F_{v\text{,cont}}}\right),$$(1)

where F_v is the observed flux and F_v,cont is the flux of the fitted continuum. Then lab-measured band strengths of the ice features were used to convert the observed optical depths to column densities by $$N=\frac{{\displaystyle \int {\tau }_v\text{d}v}}{A},$$(2)

where v is wavenumber in cm⁻¹ and A is the band strength in cm molecule⁻¹. If the ice features could be assumed to come from the same absorbing region, then relative ice abundances could be derived from column densities in a model-independent way. In the past, most of the available ice spectra consisted of spatially unresolved spectroscopy of point sources; therefore, the simplifying assumptions of pure ice absorptions within the same beam seemed appropriate for ices seen toward background stars and protostars. In this manner, ice abundances relative to H₂O have been derived in a few molecular clouds using lines of sights toward background stars, many protostars, and a minimal number of edge-on disks (EODs) using the Very Large Telescope (VLT) or the Spitzer Space Telescope.

This assumption was known to be incorrect for ices in any source with an extended physical structure, such as a disk. However, recent work has suggested that this traditional assumption of pure ice absorption may not hold in most astrophysical environments relevant to star formation. When icy mantles grow to micrometer sizes, they are seen through a combination of both absorption and scattering, which modifies their feature profiles and may strongly impact how well column densities translate into relative abundances, even toward targets with no source structure such as stars behind molecular clouds (Dartois et al. 2022). This effect would be expected to be enhanced toward sources with an internal structure, such as protostars or disks.

For disks in particular, the observations themselves are challenging. First, ices are concentrated in the midplane, which only becomes optically thin in the far-infrared. In the absence of space-based far-infrared telescopes, the only current option to observe ices is to probe them in the mid-infrared, although this samples ices in the disk atmosphere rather than in the disk mid-planes, which is optically thick in the infrared. At these wavelengths, ices have been detected by broadband, coronagraphic ground-based scattered light imaging of sufficiently bright, face-on disks, for example the Herbig Ae stars HD 142527 and HD 100546 (Honda et al. 2009, 2016), and through spectroscopy of edge-on systems (i > 70°). In the latter case, the optically thick midplane acts as a coronagraph to block thermal emission from the central star and hot dust in the disk’s midplane at 1–5 au, reducing the contrast enough to detect the signatures of cold ices beyond 20 au in absorption against the faint light from the star or inner dust rim scattered through the disk’s atmosphere (Pontoppidan et al. 2007; Sturm et al. 2023c). This technique allows the measurement of multiple ice species simultaneously, which is critical for determining CHONS abundances.

Previous spectroscopic studies of edge-on disks were limited by their technical capabilities and sample selection. Unresolved, single-slit ground-based spectroscopy detected water ice in a select number of disks, including CRBR 2422.8-3423, HV Tau C, HK Tau B, a couple of Orion silhouette disks, and the disks around the Herbig Ae/Be stars PDS 144 N and PDS 453 (Thi et al. 2002; Terada et al. 2007; Terada & Tokunaga 2012, 2017). However, detection of the main carbon dioxide (CO₂) and carbon monoxide (CO) ice features are limited by absorption from Earth’s atmosphere and the sky background at 4–5 µm. Additional observations from space of edge-on disks ASR41, 2MASS J1628137-243139, HV Tau, HK Tau, UY Aur, and IRAS 04301 with AKARI (Aikawa et al. 2012) and CRBR 2422.8-3423 with Spitzer IRS (Pontoppidan et al. 2005) confirmed the ground-based water detections and added CO₂ detections for the latter five sources.

However, the last two systems demonstrate the second difficulty with ice measurements in disks: target selection. Without ancillary data, it can be difficult to disentangle whether the observed ice bands originate in the disk itself or in other foreground material (e.g., in a surrounding cloud or a protostellar envelope). Selecting samples of mature, envelope-free disks should solve this problem, but historically it has been difficult to identify edge-on disks a priori. Extinction, including self-extinction by their outer disks and scattering from the inner edge, can result in edge-on disks being mistaken for younger protostars due to their rising shape in the mid-infrared when identified purely based on spectral energy distributions (SEDs; Robitaille et al. 2006; McClure et al. 2010, Fig. 3). Only recently have small samples of non-embedded edge-on disks been compiled and resolved at multiple wavelengths, using optical scattered light and resolved millimeter imaging (Villenave et al. 2020). For these validated disks, the sensitivity, wavelength coverage, spectral resolution, and spatial resolution of the James Webb Space Telescope’s (JWST) near- and mid-infrared spectrographs, NIRSpec and MIRI, now allow the 2D spatial distribution of ices in such disks to be mapped with a 0/1 to 0/3 pixel scale.

As part of the JWST Early Release Science Program, Ice Age, we observed the nearly edge-on (i = ~83°; Sturm et al. 2023b) protoplanetary disk HH 48 NE around a K7 star in the Chamaeleon I star-forming region (185 pc; Gaia Collaboration 2021) that has an estimated foreground-cloud visual extinction of only A_V = 5. A flared disk of 1.3–1.7″ diameter (120– 160 au radius) was seen around HH 48 NE by the Hubble Space Telescope (HST; Stapelfeldt et al. 2014) and by the Atacama Large (sub-)Millimeter Array (ALMA; Villenave et al. 2020; Sturm et al. 2023b), a size that is well suited to the JWST/NIRSpec instrument’s field of view. In the current study, we present the first results of these JWST/NIRSpec observations of the ices in the outer regions of HH 48 NE.

This paper describes the first complete gas and ice inventory of an edge-on protoplanetary disk observed with JWST/NIRSpec in a model-independent way. Due to the richness and novelty of these data, new science can be gleaned directly from the data prior to radiative transfer fitting. We discuss the interpretation of disk-integrated and spatially resolved ice features, finding that multiple scattering from many regions in the disk makes the traditional metrics of optical depths and column densities meaningless. This highlights the ultimate necessity of radiative transfer for determining ice abundances in protoplanetary disks.

The paper is organized as follows: We first present the observations in Sect. 2. Then, in Sect. 3, we provide details on the detected gas and ice features in the disk-integrated spectrum, along with the spatial variations of the ice feature depth. In Sect. 4, we discuss the implications of the observations and pose unresolved questions for future research. Finally, Sect. 5 summarizes our findings and provides concluding remarks.

2 Observations

The protoplanetary disk around HH 48 NE was observed with NIRSpec onboard JWST on 2022 July 19, as part of the Director’s Discretionary-Early Release Science (DD-ERS) program “Ice Age: Chemical evolution of ices during star formation” (ID 1309, PI: McClure). The spectroscopic data were taken with a standard 4-point dither pattern for a total integration of 2451 s, with the pointing centered on 11^h04^m23.56^s, −77° 18′07.26″. The disperser-filter combination used was G395H-F290LP, which covers wavelengths from 2.87–4.08 µm and 4.19–5.25 µm at a resolving power of R ~2700. Because the target is extended, target acquisition was not carried out, and HH 48 NE was blind-pointed into the NIRSpec IFU aperture for spatially resolved imaging spectroscopy across the disk and inner outflow extent. A second pointing was observed at an offset of 3 arcmin southeast of HH 48 NE to subtract the astrophysical background from the target flux.

The data were processed through the November 2022 development version of the JWST NIRSpec pipeline (2022_2a; 1.8.3.dev26+g24aa9bld), which included a first fix for a pronounced “ringing” throughput structure seen in the previous version of the calibrated three dimensional (3D) IFU cube product for NIRSpec. Calibration reference file database version 11.16.16 was used, which included the updated onboard flat field and throughput calibrations for absolute flux calibration accuracy estimate on the order of ~10%. The standard steps in the JWST pipeline were carried out to process the data from the 3D ramp format to the cosmic ray corrected slope image. Further processing of the 2D slope image for WCS, flat fielding, and flux calibration were also done using standard steps in the “Level 2” data pipeline calwebb_spec2. To build the calibrated 2D IFU slice images into the 3D datacube, the “Stage 3” pipeline was run step-wise, and intermediate products were investigated for accuracy. The four dithers were combined using the “drizzle” algorithm with equal weighting and full pixel regions used. The final pipeline processed product presented here was built into 3D with the outlier bad pixel rejection step turned off, as running this over-corrected and removed target flux. As a result, clusters of bad pixels are seen in some velocity channels of the IFU cube that are flagged by hand when necessary based on the 2D spatial information.

Fig. 1 Overview of the different spatial components of the HH 48 NE system using a composition of observations. The JWST/NIRSpec disk continuum at 4 µm is shown in red, the jet is shown in blue using HST/ACS observations of the [O I] line at 0.6 µm (Stapelfeldt et al. 2014), and the vertically extended disk wind is shown in green using the H2 0–0 S(9) transition. The region used to extract a disk-integrated spectrum is shown by the dashed rectangle. The continuum in this region is shown on a logarithmic scale in Fig. 4. The direction toward HH 48 SW is marked with a white arrow.

3 Results

Figure 1 presents the spatially resolved continuum image of HH 48 NE at 4 µm together with the extended H₂ 0–0 S(9) emission at 4.7 µm and the jet traced by [O I] at 0.6 µm with HST (Stapelfeldt et al. 2014). The cube is rotated by −10° using scipy cubic spline interpolation to align the major axis of the disk horizontally for visualization purposes. The total field of view (FOV) of the NIRSpec instrument is 3″ × 3″, which means that HH 48 SW (located at 2/3; Sturm et al. 2023b) lies outside of the field, but some hints of gas emission from HH 48 SW are visible toward the southwest (see Fig. 1). The HH 48 NE protoplanetary disk appears spatially resolved in the NIRSpec spectral cube, with a major axis of l·76 (296 au) and a minor axis of 1″.Ί (204 au) at 4 µm determined using a 5 σ clip. The spatial resolution of the observations varies from 0/09 at 3 µm to 0/15 at 5 µm (17–28 au), which implies that the spaxels (0.1 ″ × 0.1″) are slightly smaller than the size of the effective PSF at the longest wavelengths. The continuum is asymmetric in the lower surface, extending further toward the west than to the east (cf. the same continuum on a log scale in the top left panel of Fig. 4). This pattern is consistent with the optical HST observations (Stapelfeldt et al. 2014; Sturm et al. 2023b) and may be a result of gravitational interaction between the two sources.

A 1D, disk-integrated spectrum is extracted within a rectangle of 1.″6 × 1.″1 around the disk, centered on 11^h04^m23.38^s, −77° 18′06.46″ (see Fig. 1) by summing up the flux over the spaxels with S /N > 5. We estimate the noise level of the disk-integrated spectrum in a flat spectral region without gas and ice features (3.85–4 µm), which results in a local standard deviation of σ = 0.027 mJy, or a S/N of ~110 on the continuum. We present the disk-integrated NIRSpec spectrum of HH 48 NE in the top panel of Fig. 2. The spectrum shows a fairly flat continuum of ~3 mJy, which is consistent with the broadband Spitzer/IRAC fluxes calculated using a Gaussian decomposition of the two binary components (HH 48 NE and SW, separated 2.″3 from each other) by Dunham et al. (2016). This spectrum provides the total throughput of the stellar energy through the disk and has the highest achievable S/N on a spectrum in this data set for ice analysis. Several gas emission lines and ice absorption features are detected, which we present in the subsequent sections.

Fig. 2 Overview of the disk-integrated JWST/NIRSpec spectrum of HH 48 NE. Top panel: disk-integrated JWST/NIRSpec spectrum with the main gas lines labeled. The molecular hydrogen O and S series lines are marked in cyan. Bracket, Humphreys, and Pfund recombination series of atomic hydrogen are marked in green, red, and purple, respectively. The CO R (4.3–4.7 µm) and P (4.7–5.2 µm) branches of ro-vibrational transitions are marked in orange. The two literature Spitzer/IRAC data points are shown in dark red for comparison. The wavelength range 4.08–4.19 µm is not covered in the used spectral setup. Middle panel: disk-integrated JWST/NIRSpec spectrum with the gas lines removed in black. The procedure for removing the CO lines is described in detail in Appendix A. The positions of the main vibrational modes of the more abundant ice species are labeled, with the continua used to calculate their optical depths shown as dashed lines with corresponding colors. Bottom panel: optical depth of the ice features in the disk-integrated spectrum shown in the middle panel. The two different colors represent optical depths for two different continuum fits, one for the H2O and CO2 features in blue, and one for the weaker ice features neglecting the CO2 and H2O features in orange. Both corresponding continua are shown as a dashed line in the middle panel.

3.1 Gas inventory

Ro-vibrational molecular emission lines and hydrogen recombination lines are ubiquitously found in protoplanetary disks, arising from the warm inner disk regions (Pontoppidan et al. 2014). In edge–on disks, these – usually spatially unresolved – lines can be observed through photon scattering in a similar fashion as the thermal dust radiation. Detected gas emission lines in HH 48 NE are labeled in the upper panel of Fig. 2. Atomic hydrogen recombination lines are detected from the Bracket (n = 4), Pfund (n = 5), and Humphreys (n = 6) series. These lines are spectrally unresolved and have a constant ratio with respect to the spatially resolved scattered continuum. This indicates that the lines have a similar region of origin as the continuum emission in the inner disk, and the spatial extent of the emission is a result of photon-scattering in the outer disk. In this spectral setting, there are 7 molecular hydrogen lines from the υ = 1−0 (O; ∆J = −2) and v = 0–0 (S; ∆J = +2) series. These lines have a slightly larger spatial extent than the disk continuum (see Fig. 1), and likely result from a vertically extended disk wind or outflow cavity.

Additionally, we detect many ro-vibrational lines from the CO R (∆J = 1) and P (∆J = −1) branches, up to J = 48 (E_up = 10 056 K) and from several vibrationally excited states. These lines are broader than expected based on the instrumental resolution, which means that they may trace a warm, highvelocity component in the inner disk (υ ~ 150 km s⁻¹) or the wind-launching region. The estimated line width is consistent with high resolution VLT/CRIRES observations of similar systems (see e.g., Herczeg et al. 2011; Brown et al. 2013; Banzatti & Pontoppidan 2015) that have a range of full width at half maximum (FWHM) of 50–200 km s⁻¹ which corresponds to an emitting radius of ~0.05 au. The spatial extent of the CO lines is similar to the continuum in the upper surface, but the contribution from the lower surface is missing (see Fig. A.1). The origin of this phenomenon is unknown, but it could result from a difference in the vertical height of the CO-emitting layer compared to the continuum which could block part of the emission.

We removed the CO lines from the spectrum to accommodate a better fit of the continuum and identification of the ice absorption features, which are the main focus of this paper. For this purpose, we fit a simple LTE slab model to the lines with a constant temperature (2070 K) and column density (1.4 × 10⁻¹⁹ cm⁻²), and subtracted it from the spectrum (see Appendix A). Since the spatial extent of the CO emission is similar to the continuum in the upper surface, with a constant flux ratio of 10 (see Fig. A.1), we assume that the propagation of photons from CO molecules is similar to the continuum emission that originates in the central regions of the disk. This implies that we detect the same CO gas component throughout the disk, and its spatial extent is the product of scattering to the outer disk. Therefore, we scale the best-fit CO spectrum with a constant scaling factor and subtract it from every spaxel that contains dust emission. Additionally, we flag the hydrogen lines and present the gas line subtracted spectrum in the middle panel of Fig. 2.

3.2 Ice inventory

The middle panel of Fig. 2 clearly reveals broad absorption features due to different ice species, after subtraction of the CO and other gas phase lines. We fit the continuum around the main ice absorption features with a polynomial, using the wavelength ranges and fit orders specified in Table 1 and show the result of the continuum fit in the middle panel of Fig. 2. Since the spectrum does not include the full blue wing of the water ice feature, we fit a separate linear continuum for the water (H₂O) feature that is extrapolated down to 2.87 µm. The ice spectrum is then converted to an ice optical depth¹ using Eq. (1). The optical depth indicates the fraction of light at each wavelength that is absorbed by the icy dust grains, on a logarithmic scale, which is linearly correlated with the ice column density encountered along the light path neglecting radiative transfer effects. The resulting spectrum is shown in blue in the bottom panel of Fig. 2. The integrated optical depths for the detected features are given in Table 1, with the same band strengths as in McClure et al. (2023) for consistency. These values cannot be used to calculate column densities directly due to the flux contributions of multiple light paths through the disk, as we will show in Sect. 4.2. We would like to note that the exact choice for the continuum baseline can result in an uncertainty of up to 10% in integrated optical depth.

We describe the major observed ice features in the next subsections. The shape and depth of ice features can be significantly altered by radiative transfer and grain growth effects, especially in scattering-dominated edge-on disks (see e.g., Dartois et al. 2022), which complicates the analysis of the ice bands. A detailed study of the ice composition, environment, and abundance requires radiative transfer modeling is therefore left for a forthcoming paper.

3.2.1 H₂O and NH₃: 3 µm

The 3.0 µm absorption feature is primarily attributed to the H₂O OH stretching mode and is similar to previous spectra observed toward background stars with JWST. The feature exhibits a well-defined peak with an optical depth of 0.9 and a relatively smooth shape, displaying apeak at 3.0 µm and an absorption wing on the red side. The NIRSpec filters used do not capture the blue wing of the feature below 2.8 µm. Analyzing the red wing is complicated by a potential weak emission feature of polycyclic aromatic hydrocarbons (PAH) at 3.28 µm (Draine & Lee 1984). A comparison with PAH features at MIRI wavelengths, that are strongly correlated with this feature (Draine & Lee 1984), is required to determine its actual contribution to the continuum. Currently, a continuum is fitted avoiding potential weak absorption features within the OH stretch mode of water using the wavelength ranges specified in Table 1, illustrated as the orange line in the middle panel of Fig. 2.

Other features of the 3 µm feature can be ascribed to a combination of additional species and radiative transfer effects². The narrow feature centered at 2.93 µm is attributed to the ammonia NH stretch mode (Dartois et al. 2002). Its position and shape are uncertain since it is located right on the blue edge of the filter, which further complicates the exact location of the continuum. The feature at 3.2 µm is similar to those reported for the silhouetted disk 216-0939 in the M43 region and the edge-on disk around the Herbig star PDS 453 (Terada & Tokunaga 2012, 2017), and is interpreted as a result of large-particle-size (~ 1 µm) crystallized water ice absorption. It is known that both larger grains and a higher degree of crystallization can cause a shift from the typical wavelength of this feature (3.1 µm) to 3.2 µm (Smith et al. 1989). The overall strength compared to the peak at 3 µm is a factor of 2–3 lower than in the former sources, which suggests that the fraction of crystallized H₂O ice is lower in HH 48 NE. The broad feature at ~3.5 µm is observed in multiple background and protostellar sources, correlating with the water ice absorption (e.g., Brooke et al. 1996, 1999; Dartois et al. 2002; Shimonishi et al. 2016). This absorption feature is often attributed to ammonia hydrates, implying an OH stretching mode absorption shifted by interaction with the nitrogen atom of ammonia embedded in water ice. The appearance and optical depth of this mode is consistent with the detection of the NH stretch at around 2.97 µm for ammonia. Some potential additional contributions by CH stretching modes have also been presumed in this range, based on laboratory experiments using H atom irradiated carbon grains (Mennella 2010). Although the distinct, narrow peak of CH₃OH at 3.54 µm is absent (see Fig. 9 in McClure et al. 2023), a small contribution of CH₃OH cannot be excluded based solely on this feature and will be constrained better with the combined MIRI and NIRSpec data. Unfortunately, the less abundant water isotopologue, HDO, with a feature at 4.07 µm is not within the observed wavelength range, due to the NIRSpec G395H detector gap, and is therefore excluded from the analysis.

3.2.2 CO₂ and ¹³CO₂: 4.2–4.4 µm

The carbon dioxide (CO₂) absorption feature at 4.26 µm is detected with a similar peak optical depth as the H₂O feature. The peak position of the observed CO₂ stretching mode band is shifted compared to pure CO₂ (4.27 µm). This shift can be attributed to CO₂ in polar ice environments or by ices exposed to energetic processing at low temperatures (Ioppolo et al. 2013, 2022). The main isotopologue has a red wing in emission (see red dashed line in the middle panel of Fig. 2), which is likely caused by a phase scattering effect at a high degree of disk inclination involving the scattering opacity of large (>1 µm) grains (Dartois et al. 2022). The ¹³CO₂ isotopologue absorption band is located in this scattering wing, so for that reason we fitted the continuum for this feature separately with a local fit specified in Table 1. A local fit introduces further uncertainty in the derived column density, but it is necessary here as the main CO₂ feature contributes to the continuum of the ¹³CO₂ feature.

As for ¹²CO₂, the peak of the ¹³CO₂ feature is shifted relative to the lab-measured peak of 4.38 µm of pure ¹³CO₂ to 4.39 µm. Boogert et al. (2000) showed that this peak shift can be ascribed to absorption by ¹³CO₂ ice mixed with molecules with a large dipole moment, for example, H₂O and CH₃OH, at high concentrations. The derived ¹²CO₂/¹³CO₂ ratio is low, with a ratio in peak optical depth of 14 and a ratio in the integrated optical depth of 22. Assuming that both features arise from the same region in the disk, and using the band strengths reported in Gerakines et al. (1995) and Gerakines & Hudson (2015) for ¹²CO₂ and ¹³CO₂, respectively, we find a molecular ratio of 14 compared to the local interstellar medium (ISM) ratio of 77 (Wilson 1999) and the local dark cloud ratio of 69–87 (McClure et al. 2023). This low ¹²CO₂/¹³CO₂ ratio indicates that the main isotopologue may be saturated³ at certain positions along the line of sight, a point to which we return in Sect. 4.2.

3.2.3 CO, OCN⁻, and OCS: 4.4–5 µm

The right hand side of the bottom panel of Fig. 2 shows the spectrum in the range 4.32–5.24 µm, where features attributed to CO at 4.67 µm and ¹³CO₂ at 4.39 µm are detected. The main CO stretch feature at 4.67 µm is detected at a peak optical depth of 0.2. The feature shows signs of a potential red wing in emission, very similar to that observed for the CO₂ feature. The absorption feature appears asymmetric and the total observed FWHM of the feature (12.4 cm⁻¹) is also broader than expected from pure CO (3 cm⁻¹; Palumbo et al. 2006). This may indicate the presence of an additional component of a CO mixture with CO₂ or H₂O, similar to Pontoppidan et al. (2003) and Öberg et al. (2011a). Meanwhile, ¹³CO is not detected, but the 3σ upper limit obtained (τ < 0.03) results in a ¹²C/¹³C ratio lower limit of 6, which is consistent with the observed ratio of 14 for CO₂.

The observed CO ro-vibrational gas lines produce in some places a quasi-continuum by line overlap. We have corrected for this effect as much as possible by subtracting a slab model (see Appendix A), but the remaining systematic error of this method is hard to determine. Residual weak ice absorption should therefore be treated with care. The feature at 4.61 µm is consistent with lab measurements of the cyanate anion OCN⁻ (e.g., Grim & Greenberg 1987; Hudson et al. 2001). This feature is observed in a variety of sources, including low- and high-mass YSOs (e.g., Pendleton et al. 1999; van Broekhuizen et al. 2004; Boogert et al. 2015) and it is consistent with the OCN⁻ band observed in laboratory analogs after UV-photolysis, ion bombardment, and thermal processing of various icy mixtures (e.g., Palumbo et al. 2000; van Broekhuizen et al. 2004). The feature centered at 4.9 µm is consistent with OCS ice. A comparison with laboratory data suggests that the feature arises from OCS embedded in mixed ices and in regions of the disk at different temperatures (Palumbo et al. 1995, 1997; Ferrante et al. 2008; Garozzo et al. 2010; Boogert et al. 2022).

3.3 Spatial variation in ice features

Since the disk is resolved in both vertical and the radial directions, we can map the strength of the ice features spatially for the first time in a protoplanetary disk. In Fig. 3 we present three normalized spectra across the disk along the major axis (middle panel) and the minor axis (right panel). The overall structure of the spectrum is similar at all locations of the disk. Along both the major and the minor axis, there is no big difference in the ice features and in the CO ro-vibrational lines. Spectra extracted further away from the central source along the major axis show a steeper slope in the water ice feature at ~3.3 µm, which may indicate that the relative amount of PAH emission increases in the outer disk. The spectrum extracted from the lower surface shows an increased slope of the continuum, likely as a result of the radiative transfer effects for the disk that has an inclination of 83°.

We fit the same continua as for the disk-integrated spectrum (see Table 1) in each spaxel with a continuum signal >5σ, estimated between 3.85–4 µm. We subsequently fit a Gaussian for each absorption feature, leaving amplitude, position, and FWHM as free parameters to determine the “peak optical depth” or ice feature amplitude in every spaxel. The resulting maps for the strong ice features (H₂O, CO₂, and CO) are presented in Fig. 4. We note that the parameters (r, z) denote projected horizontal and vertical coordinates on the sky rather than intrinsic properties of the disk. Since the disk is not exactly edge-on, but rather inclined at ~83° (Sturm et al. 2023b), the upper surface of the disk (z > 0) is different from the lower surface (z < 0) in terms of absolute flux.

The peak optical depth maps consistently show a butterfly pattern, with the optical depth increasing toward the midplane in the vertical direction and outwards in the radial direction. The magnitude of the variations is only minor (<50%) considering the multiple orders of magnitude difference in the continuum flux (blue panel in the top left of Fig. 4) and dust column density along a pencil beam⁴ (note that the color bar is truncated at a peak optical depth of 0.6 for H₂O and CO₂). There are some potential east-west asymmetries in the disk’s lower surface, that could be mainly attributed to a higher continuum flux on the west side, which means that we can trace the ice features further out on that side (see discussion in beginning of Sect. 3). Within the range where continuum emission is detected toward the east, there is no significant difference between the east and west side in any of the ice features. Overall, vertical and radial cuts through the disk follow similar trends everywhere in the disk, but vary independent from each other. In the subsequent sections, we will describe the variations in optical depth with the height and radius of the disk independently.

Fig. 3 Spatial variations in the spectrum along the major and minor axes of the disk. Left: continuum flux of HH 48 NE at 4 µm on a logarithmic scale. The colored squares mark the regions used to extract the spectra in corresponding colors in the middle and right panel. Middle: three spectra extracted along the major axis of the disk at radial distances of −70, 0 and 70 au with respect to the center. Right: three spectra extracted along the minor axis of the disk at heights of −25, 50 and 80 au with respect to the dark lane.

3.3.1 Radial variations

To highlight the radial trend in the ice optical depth, we took the median of the peak optical depth along the major axis of the disk, see bottom left panel of Fig. 4. The 1 σ variations on this radial profile across the vertical extent of the disk are marked as a shaded region around the median for reference. Surprisingly, the shape is similar throughout the vertical extent of the disk. The overall trend is that the optical depth increases with radius with an increase by 20% at 100 au with respect to the center. H₂O, CO₂ and CO show identical spatial profiles within the uncertainty, with CO a constant factor of ~5 weaker relative to H₂O and CO₂. The ice absorption is in general symmetric in both sides of the disk. The slight asymmetry in the CO profile could be a result of the low S/N and the low number of pixels per bin at larger distances (see top right panel). The radial profile cannot be directly translated into a radially resolved column density profile due to complex radiation transfer, including the fact that because of the edge-on geometry of the disk all sight-lines sample the full radial extent of the disk. We will return to this point in Sect. 4.3.

3.3.2 Vertical variations

To highlight the vertical variations in the ice optical depth, we took the median of the peak optical depth along the minor axis of the disk (see bottom right panel of Fig. 4). The 1 σ variations on this profile across the radial extent of the disk are marked as a shaded region around the median. We note that the shape is similar throughout the radial extent of the disk. The overall trend in the three major ice species is that the vertical profiles peak around the dark lane and decrease quicker on the upper surface than on the lower surface. We note that the origin is taken as the dark lane, which coincides with the highest peak optical depths, but the stellar position could likely be slightly above the dark lane due to projection effects. The three major ice species (H₂O, CO₂ and CO) have identical profiles within the uncertainty from −70 au in the lower surface till 80 au in the upper surface. We will return to the interpretation of this result in Sect. 4.3.

Fig. 4 Overview of the variations in the relative strength of the three strongest ice features, H2O, CO2, and CO. Top left panel: continuum flux of HH 48 NE at 4 µm on a logarithmic scale. The white line denotes the 5σ detection limit on the continuum, determined in a featureless part of the spectrum, namely the 3.85–4 µm region. This field of view is the same as the dashed rectangle in Fig. 1. Other panels in the top row: peak optical depth per spaxel of H2O, CO2, and CO ice, respectively, determined using a Gaussian fit to the optical depth spectrum of each spaxel. Bottom left panel: median combined radial profile of the peak optical depth with the 1σ spread around this line at different disk heights as a shaded region. The black line underneath shows the radial profile of the continuum at 4 µm at z = 50 au. Bottom right panel: median combined vertical profile of the peak optical depth with 1σ spread around this line at different disk radii as a shaded region. The black line underneath shows the vertical profile of the continuum at 4 µm at r = 0 au. We note that z = 0 is chosen as the dark lane between the upper and lower surface: in this frame, the star is located at z ~10 au due to projection effects.

4 Discussion

Observations of ices in protoplanetary disks are crucial to understanding the chemistry that sets the composition of planetesimals and planetary atmospheres, and we show that with JWST we can trace the ice directly in late-stage edge-on disks. Traditionally, the observed ices were believed to come from the regions that were simultaneously cold enough for particular ice species to freeze out but also sufficiently high up in the optically thin disk atmosphere to absorb stellar light scattered into the outer disk. Under these assumptions, the column densities of these ices would have provided relative measurements of CHONS-rich ice abundances in the disk layers just above the planet-forming midplane. However, the observations presented in this work reveal a low ¹²CO₂/¹³CO₂ ratio, suggesting that the ¹²CO₂ feature is saturated, even though the minimum flux is significantly above JWST/NIRSpec’s sensitivity limit. This raises the question of whether ratios between column densities, calculated from the observed optical depths using Eq. (2), are at all meaningful proxies for relative ice abundances. Additionally, contrary to our expectations, only minor spatial variations are observed in the strength of the ice absorption features. This leads to a second question of which disk region(s) the observations are really sampling. In the subsequent sections, we explore which fundamental constraints on the answers to these questions can be extracted from such spatially resolved observations in a model-independent fashion.

4.1 Where in the disk do these ice features originate?

The radiative transfer in edge-on disks is complicated, as only very specific light paths reach the observatory. In contrast to the pencil beam scenario, e.g. background star observations behind a dark cloud, all of the observed photons come from a similar region in the warm inner disk, before they are scattered into our line of sight (see Sturm et al. 2023c, for a deconstruction of the continuum). Figure 5 illustrates examples of the photon propagation that could explain the observed absorption features and spatial variations. Most of the photons emitted by the star and the warm inner disk are absorbed along the way or scattered out of the line-of-sight. Therefore, the observed ice absorption features may represent both the optically thin regions above and beyond the optically thick disk midplane (see right panel Fig. 5) and potentially some part of the disk midplane itself, depending on the number of times each photon is scattered. Since the disk is axisymmetric, this means that all of the light that we observe has traveled through a similar radial extent of the disk, independent of the projected disk radius at which the light appears on the 2D JWST/NIRSpec array. The slight increase of optical depth with radius (seen in the bottom left panel of Fig. 4) can be understood as a geometric effect: the further away we look from the star, the longer the light has traveled on average through the disk which effectively increases the ice column density along the light path (path a + b > c in the left panel of Fig. 5) and the more likely the light has encountered ice along the way.

The light that we observe at every height in the disk took a complex path with one or more scattering events in-between (right panel of Fig. 5). The closer the light appears to come from the outer disk midplane (the dark lane at z = 0 au), the greater the chance that the light has scattered through cold parts of the disk. This explains why we see a peak in the optical depth values at a height of ~ 10 au above the midplane (bottom right panel of Fig. 4). Since the disk is not exactly edge-on, but rather inclined at ~83°, light that has traveled through the lower surface has a longer path through the cold midplane compared to light that traveled through the upper surface. This results in slightly elevated optical depths on the lower surface compared to the upper surface, as seen in the bottom right panel of Fig. 4. A final important result from Fig. 4 is that there are no clear vertical cut-offs – or snowlines – visible in the ice feature depths, indicating that ices may extend up to 100 au in the disk, as discussed in Sect. 4.3.

Fig. 5 Cartoon of the physical origin of the radial (left) and vertical (right) profiles of the ice optical depths. Left: since we observe ice absorption against the warm inner regions of a (close to) edge-on disk, all lines of sight trace ice along the full radial extent of the disk. This results in a largely constant optical depth profile along the radial direction. Since the path of the green photon trace is longer in the ice-rich outer disk (a + b) than the direct path (c), the traced ice column densities increase with radial separation of the observation. Right: direct lines of sight (blue) toward the star and warm inner disk are blocked by the optically thick disk midplane (white region). The observed light follows therefore a specific path through the optically thin disk regions, including one or multiple scattering events. The flux in a JWST spaxel is composed of photons with different paths (orange and yellow), that have not necessarily crossed equal lengths in the cold parts of the disk. This means that photon paths that cross large parts of the ice-rich regions of the disk can be saturated locally, even though the optical depth of the combined contributions is a defined, small number (see inset). The optical depth increases toward the midplane, because the light paths through the cold disk layers are more likely compared to paths higher up in the disk (green). Since the scattered light from the lower surface is tilted toward the cold midplane, and similarly, light from the upper surface is tilted away from the cold midplane, the profile is asymmetric with deeper ice features in the lower surface. Inset: qualitative schematic of the different flux contributions in the spectrum. The spectrum always has a contribution from light that passes through little ice (ambient light continuum; orange line) and light that passes through the cold regions of the disk, usually requiring multiple scattering events along the way (yellow line). The ratio of the contributions determines the optical depth in the observations if the feature is saturated. Photons with a single scattering event are less likely to reach regions with ice in the disk than photons with multiple scattering events (see also Fig. 3 in Sturm et al. 2023c).

4.2 Optical depths do not trace total ice column in outer disk

Our data analysis suggests that ice feature optical depths, and by extension formal column densities, in protoplanetary disks, cannot be used to derive relative ice abundances in a straight-forward way. We have illustrated in the previous section (Sect. 4.1) that the observed absorption features are likely not representative of the full extent of the disk, but trace only optically thin regions above the midplane and in the outer disk. In addition, multiple light paths are combined in one line of sight, which implies that each spaxel contains multiple light contributions with different levels of ice absorption (see the inset panel of Fig. 5 for an illustration of this concept). This means that there may be an ambient-light continuum underneath the ice features with much less ice absorption that fills in the ice absorption features, which has important consequences for the data analysis. CO₂, H₂O, and potentially CO are locally saturated along the light path, indicating that light scattered through the cold parts in the disk is saturated in the ice absorption features (yellow line in the inset panel of Fig. 5). The ice optical depth is then fully determined by the fraction of light scattered through the icy regions in the disk, relative to the light scattered through the ice-free regions (i.e., the ratio of orange to yellow flux in the inset panel of Fig. 5); the measured optical depths do not represent the total column density of ice along the light path.

Increasing the amount of ice in the cold regions of the disk will not have a significant effect on the measured optical depth: the depth of the ice feature in the inset of Fig. 5 is determined primarily by the level of the orange line rather than the depth of the yellow line. This effect was already predicted using radiative transfer models of HH 48 NE in Sturm et al. (2023c), where we studied the sensitivity of the ice features on the physical disk parameters. We used a model tailored to HH 48 NE with constant ice abundances throughout the disk (taken as 80, 22, and 99 ppm for H₂O, CO₂, and CO respectively), assuming that all molecules are frozen out below their desorption temperature and below a UV threshold of A_v = 1.5. Varying the abundances in the model, we find that similar peak optical depth values of ~0.9–1.2 for H₂O and CO₂ in HH 48 NE models require abundances of 2100 times the assumed ISM values: high enough to saturate these ice absorption features locally in the disk. An accurate column density can only be determined if we know what fraction of the light is scattered through regions without ice in the disk; i.e. we must determine the ambient light flux indicated by the orange line in the inset of Fig. 5, which yields the zero point for the flux of the yellow line. The radiative transfer models of HH 48 NE show that photons with a single scattering event – and therefore less ice absorption than photons with multiple scattering events – account for ~1/3 of the flux at NIRSpec wavelengths, which is comparable with the ambient light contribution observed (see Fig. 3 in Sturm et al. 2023c). The fraction of ambient scattered light is impossible to know without modeling the radiative transfer through the disk. Therefore, column densities derived using Eq. (2) should be interpreted as lower limits.

An important implication of this is that naive ice column density ratios, calculated by multiplying the observed integrated optical depth in source integrated spectra with the band strength, do not necessarily correspond to similar relative ice abundance ratios in the disk. The region where molecules are frozen out is different for every molecule, and therefore there are also differences between the ratio of ambient scattered light to scattered light that encounters locally high column densities of each ice species. Since the ambient light continuum is similar for the CO₂ and H₂O bands – due to their proximity in wavelength, and their similar snowline location in the disk – this results in comparable optical depths in the two features (~0.9), even though their actual relative abundance may be very different. The ratio of light that reaches the ices near the midplane will be even more different for volatile species like CO.

With these caveats in mind, we calculate formal column densities, presenting them as lower limits, to enable comparison with literature observations. The lower limits on the column of H₂O, CO₂, and CO are >1.9 × 10¹⁸ cm⁻², > 2.1 × 10¹⁷ cm⁻² and >2.3 × 10¹⁷ cm⁻², respectively. The detection of ¹³CO₂ with a column density of >1.3 × 10¹⁶ cm⁻² is crucial, as this allows us to determine the actual CO₂ column density, using ISM assumption for the ¹²C/¹³C isotope ratio. Since the main isotopologue is locally saturated and we know that ¹²CO₂ and ¹³CO₂ should be located in the same region in the disk, we can correct for the ice-free ambient-light continuum for the ¹³CO₂ feature by removing a constant flux (orange line in the inset panel of Fig. 5) underneath the peak of the CO₂ feature. Removing the ice-free ambient-light continuum results in an increase of the peak optical depth of CO₂ from 0.06 to 0.010 and an increase in the integrated optical depth by a factor of 1.6. The estimated ¹³CO₂ column density using this approach is 2.1 × 10¹⁶ cm⁻². If we assume the ISM ¹²C/¹³C ratio of 77 (Wilson 1999), then we can find an estimate of the actual ¹²CO₂ column density of 1.6 × 10¹⁸ cm⁻².

We note that this estimated column density is one order of magnitude higher than the lower limit found from multiplication of the observed integrated optical depth of the ¹²CO₂ feature with the band strength, illustrating the importance of detecting the less abundant isotopologues with JWST. Inferring molecular abundances from the main isotopologue ice features, as traditionally done in the literature, may drastically underestimate the amount of ice present and should be avoided. Therefore, it is critical to observe either less optically thick, rare isotopologues like HDO, ¹³CO₂ and ¹³CO or a disk of sufficiently low mass in order to determine the correct order of magnitude for the column densities of ices in the comet-forming regions of disks.

Fig. 6 Spatial variations in the CO ice peak optical depth (colorscale) compared to ALMA CO J = 2−1 observations presented in Sturm et al. (2023b; contours). The beam of the ALMA observations (0/51 × 0.″31) is shown in the bottom left corner. The JWST/NIRSpec observations have an effective resolution of 28 au or 1.5 spaxels.

4.3 Excess ice in disk atmosphere?

Assuming that the ice absorption takes place in the optically thin disk outer regions, the vertical snow surfaces caused by thermal desorption and photo-desorption should produce a noticeable drop in the peak optical depth with disk height (z in Fig. 4), as explained in Sect. 4.1. However, the three most prominent ice features H₂O, CO₂, and CO extend vertically over 125 au with only a minor decrease in strength. This is especially curious for CO ice, which is thought to sublimate thermally at ~20 K in protoplanetary disks (Collings et al. 2003; Pinte et al. 2018; Minissale et al. 2022, and references therein) and is therefore limited to about one hydrostatic scale height, which lies at ~ z/r 0.15 (Pinte et al. 2018, Law et al. 2021 and see Sturm et al. 2023c). The observed CO ice is even seen at similar maximum height as the CO J = 2−1 gas observations reported in Sturm et al. (2023b; see Fig. 6), which corresponds to temperatures of >50 K.

There are four potential scenarios that could result in a similar resolved profile for the three main ice features. Here, we describe the possibilities and identify the one that is the most viable, given what is known about this particular source.

First, the extra absorption provided by a foreground cloud or envelope with high ice column densities could explain the relatively constant peak optical depths. In that case, the whole disk would function as the light source, and all the emission that is scattered from the disk would pass through a similar cloud, which would give a relatively similar ice absorption at each spaxel across the disk. However, this option is unlikely given the measured visual extinction toward stars in the neighborhood (A_V ~3) by Gaia (Gaia Collaboration 2021) and the low visual extinction derived from the constraints on the disk geometry and stellar properties (A_V ~5, Sturm et al. 2023b). The integrated optical depths of the ice features (Table 1) are a factor of 5 above the observed trend with AV toward background stars (Boogert et al. 2015), which means that a potential foreground cloud contributes <20% of the total depth to the ice absorption. Additionally, we find that the absorption features of CO₂ and H₂O are saturated at a flux >0, which is only possible if absorption in the disk dominates. If a foreground cloud dominates the absorption, the absorption features would saturate at 0 due to the lack of ice-free ambient-light continuum. The observed optical depths at z = 125 au are therefore likely not due to a foreground cloud.

Second, the observed profile could be a result of radiative transfer in the disk. If most of the ice absorption happens in regions relatively close to the star, and the observed features are scattered through the disk in a similar way as the continuum, the vertical distribution could be uniform. This is, however, not in line with our current knowledge of the temperature structure of the disk (Sturm et al. 2023b), because the CO snowline is estimated to be at r ~ 50 au, comparable to observed snowlines in other systems (Qi et al. 2013; van’t Hoff et al. 2020; Podio et al. 2020; Zhang et al. 2021a) and the scattering of the ice band should therefore be inherently different than the continuum light from the warm inner disk and the CO₂ and H₂O features. Alternatively, light could have previously been scattered down below the CO snow surface, before scattering into the disk atmosphere and then again into the line-of-sight toward us (yellow line in the right panel of Fig. 5). However, radiative transfer models including such anisotropic scattering show layered disk ice segregation as a result of temperature gradients and photo-desorption (Ballering et al. 2021; Arabhavi et al. 2022; Sturm et al. 2023c). Even though photon packages may have encountered multiple scattering events along the light path (see also Sturm et al. 2023c), there is no reason to assume that scattering to deeper layers is more efficient than scattering through the upper layers alone, which would still result in noticeable vertical differences.

Third, a fraction of the CO ice could be trapped in ice with a much higher desorption temperature, and therefore survive at much higher temperatures than that expected of a pure CO ice. Laboratory experiments of CO and CO₂ desorption with H₂O show that CO deposited on top of a H₂O ice base layer can diffuse into the H₂O matrix, remaining trapped due to structural changes of the H₂O ice while heated, and effectively co-desorbs at the crystallization temperature of H₂O (150 K; e.g., Collings et al. 2003; Palumbo et al. 2006; Fayolle et al. 2011; Lauck et al. 2015). Similar diffusion is also seen for CO deposited onto CO₂ (Palumbo et al. 2006; Simon et al. 2019), where the diffusion into CO₂ is more efficient than into H₂O, and CO can also be trapped in a CH₃OH ice matrix (Ligterink et al. 2018). Entrapment of part of the CO in ices with a higher desorption temperature would still be in line with observations of gas-phase CO that show that most of the CO desorbs at ~20 K, given that the features are saturated (see the discussion in Sect. 4.2). The hypothesis that CO is trapped in the ice matrix is supported by the width of the CO feature, which is broader than that of pure CO, and also by shifts in the peak wavelengths of the ¹²CO₂ and ¹³CO₂ features, which can occur due to a H₂O or CO rich environment. Such profile signatures were seen by McClure et al. (2023) in the CO and CO₂ ice toward background stars in the same star-forming region, and ultimately found to originate in mixed ices, after comparison with laboratory data. Due to the complicated and extended disk geometry, a detailed retrieval of the ice mixtures from the profiles must be done in conjunction with a full radiative transfer model, which we reserve for a future paper.

The trapping of CO in mixed ices would not explain the lack of snow surfaces for CO₂ and H₂O, though, unless there is a mechanism for those ices to remain trapped in a less volatile material on the grain surface. However, chemical modeling including vertical turbulent mixing and diffusion shows that icy grains can reach the observable disk surface in turbulent disks, effectively mixing ices upward faster than they can be photo-desorbed in the upper layers (Semenov et al. 2006, 2008; Furuya et al. 2013; Furuya & Aikawa 2014; Woitke et al. 2022). The grain sizes of the icy grains inferred from the H₂O and CO₂ features (>1 µm; see Sect. 3.2), that are thought to trace mainly the disk atmosphere, could be a direct result of mixing ice rich micron-sized dust grains to the disk upper layers. While the timescales for thermal desorption are short, this could potentially explain the lack of a clear CO₂ and H₂O snowsurface. Regardless of the mechanism, if it is confirmed that a fraction of the ice survives up to the disk surface, this would have major consequences on the disk chemistry and models of planet formation, as we discuss in Sect. 4.4 below.

4.4 Implications for the outer disk chemistry

Using the initial ice radiative transfer models of HH 48 NE presented in Sturm et al. (2023c) and described in more detail in Sect. 4.2, a first order estimate of the molecular abundance can be derived using the relation between the peak optical depth of the absorption feature and the molecular abundance. We find that abundances should be significantly higher than the ISM abundance (taken as 80, 22, and 99 ppm for H₂O, CO₂, and CO respectively) in order to reproduce the observed optical depths of the three main ice species and the level of saturation traced with the ¹³CO₂ ice observation. Additionally, the observed ice features are stronger than anticipated based on current chemical models (e.g. Ballering et al. 2021; Arabhavi et al. 2022). This could mean that ices are efficiently inherited from the prestellar cloud stage (Ballering et al. 2021). However, the results suggest that there is even more ice (>2×) than can be explained by inheritance from the pre-stellar phase alone.

If more ices are frozen out at the disk midplane than expected, this could impact the interpretation of the CO, C, and H₂O gas depletion seen in the outer regions of some disks with ALMA (e.g., van Zadelhoff et al. 2001; Miotello et al. 2017; Sturm et al. 2023a; Du et al. 2017; van Dishoeck et al. 2021). Efficient conversion from volatile gas to the ice phase has long been suggested as an explanation for the CO and H₂O depletion in the outer disk gas, requiring a combination of freezeout, chemical conversion, grain growth, settling, radial drift, and vertical mixing (Bosman et al. 2018; Schwarz et al. 2018; Krijt et al. 2018, 2020). The non-detection of CH₃OH in our data suggests that, if chemical conversion is the primary depletion mechanism, CO is more likely converted into CO₂ than in more complex molecules like CH₃OH.

The co-existence of CO ice in the emitting layer of CO gas due to entrapment in less volatile ices, as discussed in Sect. 4.3, could be a another mechanism contributing to the low CO abundances measured in the gas phase of a large population of T Tauri stars. However, lab measurements suggest that entrapment could only result in CO gas depletion by less than a factor of 2 (Simon et al. 2019), which is insufficient to explain the 2 orders of magnitude depletion observed in older systems like TW Hya and DL Tau (Kama et al. 2016; Sturm et al. 2022); in several T Tauri systems, the inner disk is also depleted in atomic C, which is still best explained by efficient trapping of icy dust grains at the edges of disk gaps (McClure 2019; McClure et al. 2020).

Determining which of these mechanisms is operating in HH 48 NE will be essential for procedures that measure disk gas masses using volatiles like CO, as the CO gas mass is frequently used as one proxy for the bulk H₂ gas mass in disk surveys, in the absence of better tracers like HD (Bergin et al. 2013; McClure et al. 2016; Sturm et al. 2023a). As one of the most fundamental properties of disks around young stars, good constraints on the total gas mass are crucial to understanding disk evolution and the potential of any given disk to form planets.

Finally, the observed low ¹²CO₂/¹³CO₂ ratio of 14 compared to the ISM ratio of 77 (Wilson 1999), is interpreted in Sect. 4.2 to be likely a result of local saturation of the main feature. However, isotopic ratios in protoplanetary disks can be used as tracers of chemical fractionation processes during planet formation that are reflected in observed planetary atmospheres (see e.g., Zhang et al. 2021b; Yoshida et al. 2022). If future radiative transfer modeling suggests that the ¹²CO₂/¹³CO₂ isotopic ratio is genuinely lower than in the ISM, rather than deriving from local saturation effects, then this could indicate fraction-ation or isotope-selective destruction processes at work in this disk (Boogert et al. 2000). In that event, edge-on disks could represent an exciting new avenue to explore solid-state fraction-ation in a way that is comparable to sample return missions to Solar System bodies.

4.5 Comparison with previous observations

In this paper, we present the first detection of ¹³CO₂ ice in a Class II disk, along with the first firm localization of CO ice to the disk, and the tentative identification of trace ice species OCN⁻, OCS, and NH₃. While CO ice was previously detected toward a disk line of sight with Very Large Telescope (VLT; Thi et al. 2002), it was later shown that the majority of the CO ice must be present in intervening cloud material (Pontoppidan et al. 2005). With AKARI, CO and ¹³CO ice were only seen toward younger, embedded sources (Aikawa et al. 2012).

A qualitative comparison with previous observations reveals that the near-IR spectrum of HH 48 NE is remarkably similar to IRAS 04302, an edge-on Class I protostar, but bears less resemblance to the other sources observed by AKARI (Aikawa et al. 2012). This includes the absorption signatures near 2.9 and 3.1 µm, seen in HH 48 NE and IRAS 04302, but not in other systems, and the extent of the scattering versus hydrate contribution in the 3 µm red wing. The red scattering wing and relative optical depths of the CO₂ and CO features are also similar. This similarity is somewhat surprising given that IRAS 04302 has an additional envelope component (Wolf et al. 2003) that is not seen in HH 48 NE (Sturm et al. 2023b). It may suggest that the envelope of IRAS 04302 is tenuous enough not to contribute to the near-infrared ice features, or it may indicate that the envelope ices in Class I sources are similar to disk ices in Class II sources.

4.6 Future prospects and outlook

Clearly, observations of additional edge-on systems are needed to understand the relative importance of intrinsic disk properties vs. viewing geometry in setting the near-IR spectrum. Including more sources in the sample will help to constrain which properties are inherent to protoplanetary disks, and which should be attributed to source-specific properties like age, envelope material, and dynamical interaction with a binary component. Future comparison of HH 48 NE’s near-infrared with the comparable NIRSpec spectrum from the neighboring edge-on protostar Ced 110 IRS 4, also observed as part of the Ice Age JWST ERS program, will enable us to quantify the similarity of ices between the Class I and Class II stages.

Additionally, many ice species discussed in this work have corresponding absorption features in the JWST/MIRI wavelength range, 5–28 µm. Longer wavelengths may trace different regions in the disk due to differences in scattering and there could be a significant contribution from direct thermal dust emission at longer wavelengths (see Fig. A.1. in Sturm et al. 2023c). However, a future comparison between NIRSpec and MIRI observations will help to constrain the strength of the potential PAH emission feature at 3.28 µm, since that feature correlates with the features at 6 and 8 µm. Additionally, it could confirm our conclusions about the presence of NH₃ and lack of CH₃OH, as both of these species have stronger absorption features at 9.0 and 9.7 µm, respectively.

In this work, we have clearly demonstrated that neglecting the role of radiative transfer in the interpretation of ice observations toward sources with substantial physical structure, like edge-on protoplanetary disks, can lead to incorrect abundances that could be off by an order of magnitude, in the case of CO₂. Due to the radiative transfer, absorption feature shapes and depths are affected, which complicates determination of abundance ratios and ice environments. In forthcoming papers, we will study the degree to which the effect of the radiative transfer on the feature shape can be disentangled from the effects of the chemical environment of the ices, based on the results of this paper. Additionally, we will use these observations to improve our already existing ice model (Sturm et al. 2023c) to study the spatial variations in ice optical depth and molecular abundances in more detail and derive an ice to silicate ratio using the additional JWST/MIRI spectra.

5 Conclusions

We have presented the first spatially resolved observations of the main ice species in a mature, envelope-free edge-on Class II protoplanetary disk, HH 48 NE, using JWST/NIRSpec. These results demonstrate the novel capabilities of JWST for revealing the ice inventories of disks. After analyzing the observed ice features and spatially resolved trends, we can conclude the following:

• A variety of ice features are detected including the major ice components H₂O, CO₂, and CO and multiple weaker signatures from less abundant ices, NH₃, OCN⁻, and OCS. Additionally, we detected the isotopologue ¹³CO₂ ice absorption feature for the first time in a protoplanetary disk.

• The combination of an integrated absorption ratio of ¹²CO₂/¹³CO₂ = 14, in contrast to the ISM value of 77, plus similar peak optical depths of H₂O and CO₂ indicates that the main ¹²CO₂ and H₂O features are saturated, despite an excellent S/N at the bottom of both features. Their saturation implies that the spectrum has a strong contribution from an ice-free ambient scattered-light continuum, which reduces the apparent optical depth, despite the high ice column densities encountered by light paths through the cold disk regions. Using the ¹³CO₂ absorption feature, we show that the naive calculation of column densities from the observed integrated optical depth of the ¹²CO₂ feature results in an underestimation of the total ice column in the disk by more than one order of magnitude. Detecting the less abundant isotopologues with JWST is therefore crucial for measuring molecular abundance ratios in the ice.

• The shape of the crystallized water ice feature and the red shoulder of the CO₂ feature indicate that a significant fraction of the probed dust grains along the photon path are larger than ~1 µm.

• Radial variations in ice composition and abundance are hard to trace in edge-on disks due to the complex radiative transfer effects, assuming that all light originates from the warm inner disk and travels through the entire radial extent at all observed projected radii on the sky. A slight increase in feature depth with projected radial separation from the star is interpreted as a result of an increased light path length through the cold regions in the outer disk.

• Vertical variations in ice optical depth are also small and H₂O, CO₂, and CO reveal similar distributions, although these molecules have significantly different desorption temperatures. While a contribution from ices in a foreground cloud and complicated radiative transfer scattering making the ice features appear to come from higher layers cannot be excluded entirely, we find that CO ice likely survives higher up in the disk than expected based on their desorption temperature of ~20 K. This may be due to the entrapment of CO in less volatile ices such as H₂O and CO₂.

• The ice features are deeper than anticipated from radiative transfer modeling assuming typical interstellar ice abundances with respect to hydrogen. The high column densities of ice observed in HH 48 NE suggest that we see a result of gaseous CO and H₂O conversion to ices, confirming the impact of earlier proposed mechanisms to explain the observed low CO gas fluxes. If it is verified, by modeling of the ice profiles, that part of the CO ice is trapped in CO₂ and H₂O ice and survives high up in the disk, then CO diffusion into the ice could also contribute to the CO depletion seen in the gas phase.

Understanding the gas/ice ratio of important volatiles is crucial in our understanding of the chemistry in planet-forming regions and the physical properties of protoplanetary disks such as the gas mass. Many of the features of these resolved ice observations and their implications on the chemistry and physics of protoplanetary disks can be readily understood in a model-independent way, as demonstrated in this work. Future comparison over the full wavelength coverage of JWST, including MIRI data, will further constrain some of the ice abundances, as well as the ice/rock ratio. However, the spatial origin of some observed signatures remains uncertain and will require advanced radiative transfer modeling, as well as a larger sample of edge-on disk observations, to fully extract relevant information about the birthplace of planets and comets.

Acknowledgements

We would like to thank the anonymous referee for suggestions that improved the manuscript. Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA), by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101019751 MOLDISK). M.K.M. acknowledges financial support from the Dutch Research Council (NWO; grant VI.Veni.192.241). D.H. is supported by Center for Informatics and Computation in Astronomy (CICA) grant and grant number 110J0353I9 from the Ministry of Education of Taiwan. D.H. acknowledges support from the National Technology and Science Council of Taiwan through grant number 111B3005191. E.D. and J.A.N. acknowledge support from French Programme National ‘Physique et Chimie du Milieu Interstellaire’ (PCMI) of the CNRS/INSU with the INC/INP, co-funded by the CEA and the CNES. M.N.D. acknowledges the Swiss National Science Foundation (SNSF) Ambizione grant no. 180079, the Center for Space and Habitability (CSH) Fellowship, and the IAU Gruber Foundation Fellowship. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). S.I., E.F.vD., and H.L. acknowledge support from the Danish National Research Foundation through the Center of Excellence “InterCat” (Grant agreement no.: DNRF150). M.A.C. was funded by NASA’s Fundamental Laboratory Research work package and NSF grant AST-2009253.

Appendix A CO ro-vibrational gas line fitting

The HH 48 NE observations include a forest of ro-vibrational CO lines from 4.3 µm to 5.28 µm that are spectrally overlapping and therefore create a quasi-continuum underneath the lines. This results in uncertainty on the continuum for the analysis of the ice absorption bands. In an attempt to remove the gas lines, we followed a similar approach as in Grant et al. (2023) by fitting a LTE slab model to the observed spectrum with a constant temperature, column density and emitting area. First, we estimated the continuum by selecting regions in the brightest CO lines and selecting regions with the least amount of line overlap and ice absorption features (see the purple dots and gray areas in Fig. A.2, respectively). Next, we fitted LTE slab models to the lines in the regions marked with a gray horizontal line in the middle panel of Fig. A.2 using a grid of models varying the column density, excitation temperature and linewidth. Since we observe an unknown small fraction of the CO line flux scattered through the outer disk, we normalize the fitted spectrum to the observed spectrum minimizing the reduced χ² value.

The best-fitting model has an excitation temperature of 2070 K, a column density of 1.4 × 10¹⁹ cm⁻², with a line FWHM of 150 km s⁻¹. We note that these values result in a model that fits the CO emission lines well and improves the quality of the continuum, but that these values should not be trusted to imply physical quantities or regions of origin. The radiative transfer can have a significant impact on the line ratios over the spectral region that is used for the line fitting. We can conclude that the emitting area has to be warm (>1500 K) to fit the lines with high E_up and that the lines are spectrally resolved indicating high velocities.

The continuum-subtracted CO emission follows a similar spatial distribution as the continuum at similar wavelengths (see Fig. A.1), which indicates that the CO likely comes from the obscured inner disk or jet launching region and is scattered through the outer disk in our line of sight in a similar way as the continuum. The CO gas emission is less prominent in the lower surface, which indicates that the emitting region may be higher up than the emitting region of the continuum. We used the ratio of the CO lines to the continuum (right panel of Fig. A.1) to scale the best-fitting CO spectrum and remove the lines from the continuum in every individual spaxel.

Fig. A.1 Comparison between the spatial distribution of the CO emission and the continuum emission. Left: spatial distribution of the median continuum flux in the region with CO gas emission (4.36–5.2 µm). Middle: spatial distribution of the median of the peak fluxes of the CO gas emission lines. Right: ratio of the left two panels showing a constant ratio in the upper surface, but a higher ratio in the lower surface due to missing CO flux.

Fig. A.2 Overview of the CO gas line fitting procedure following the steps described in the text. Top panel: observed source integrated spectrum of HH 48 NE with selected troughs in between the emission lines in purple. Regions used to fit the continuum underneath the lines are marked in gray. Middle panel: continuum subtracted CO emission spectrum, corrected for the optical depth in the CO ice feature (blue region). The best fitting model is shown in red, regions taken into account in the fitting are marked with a gray horizontal line. Bottom panel: residual spectrum of subtracting the model from the data. The standard deviation over this region is 0.05 mJy, compared to 0.27 mJy in the CO line-free continuum at ~4 µm (see Fig. 2).

References

All Tables

All Figures

Fig. 1 Overview of the different spatial components of the HH 48 NE system using a composition of observations. The JWST/NIRSpec disk continuum at 4 µm is shown in red, the jet is shown in blue using HST/ACS observations of the [O I] line at 0.6 µm (Stapelfeldt et al. 2014), and the vertically extended disk wind is shown in green using the H2 0–0 S(9) transition. The region used to extract a disk-integrated spectrum is shown by the dashed rectangle. The continuum in this region is shown on a logarithmic scale in Fig. 4. The direction toward HH 48 SW is marked with a white arrow.

In the text

Fig. 2 Overview of the disk-integrated JWST/NIRSpec spectrum of HH 48 NE. Top panel: disk-integrated JWST/NIRSpec spectrum with the main gas lines labeled. The molecular hydrogen O and S series lines are marked in cyan. Bracket, Humphreys, and Pfund recombination series of atomic hydrogen are marked in green, red, and purple, respectively. The CO R (4.3–4.7 µm) and P (4.7–5.2 µm) branches of ro-vibrational transitions are marked in orange. The two literature Spitzer/IRAC data points are shown in dark red for comparison. The wavelength range 4.08–4.19 µm is not covered in the used spectral setup. Middle panel: disk-integrated JWST/NIRSpec spectrum with the gas lines removed in black. The procedure for removing the CO lines is described in detail in Appendix A. The positions of the main vibrational modes of the more abundant ice species are labeled, with the continua used to calculate their optical depths shown as dashed lines with corresponding colors. Bottom panel: optical depth of the ice features in the disk-integrated spectrum shown in the middle panel. The two different colors represent optical depths for two different continuum fits, one for the H2O and CO2 features in blue, and one for the weaker ice features neglecting the CO2 and H2O features in orange. Both corresponding continua are shown as a dashed line in the middle panel.

In the text

Fig. 3 Spatial variations in the spectrum along the major and minor axes of the disk. Left: continuum flux of HH 48 NE at 4 µm on a logarithmic scale. The colored squares mark the regions used to extract the spectra in corresponding colors in the middle and right panel. Middle: three spectra extracted along the major axis of the disk at radial distances of −70, 0 and 70 au with respect to the center. Right: three spectra extracted along the minor axis of the disk at heights of −25, 50 and 80 au with respect to the dark lane.

In the text

Fig. 4 Overview of the variations in the relative strength of the three strongest ice features, H2O, CO2, and CO. Top left panel: continuum flux of HH 48 NE at 4 µm on a logarithmic scale. The white line denotes the 5σ detection limit on the continuum, determined in a featureless part of the spectrum, namely the 3.85–4 µm region. This field of view is the same as the dashed rectangle in Fig. 1. Other panels in the top row: peak optical depth per spaxel of H2O, CO2, and CO ice, respectively, determined using a Gaussian fit to the optical depth spectrum of each spaxel. Bottom left panel: median combined radial profile of the peak optical depth with the 1σ spread around this line at different disk heights as a shaded region. The black line underneath shows the radial profile of the continuum at 4 µm at z = 50 au. Bottom right panel: median combined vertical profile of the peak optical depth with 1σ spread around this line at different disk radii as a shaded region. The black line underneath shows the vertical profile of the continuum at 4 µm at r = 0 au. We note that z = 0 is chosen as the dark lane between the upper and lower surface: in this frame, the star is located at z ~10 au due to projection effects.

In the text

Fig. 5 Cartoon of the physical origin of the radial (left) and vertical (right) profiles of the ice optical depths. Left: since we observe ice absorption against the warm inner regions of a (close to) edge-on disk, all lines of sight trace ice along the full radial extent of the disk. This results in a largely constant optical depth profile along the radial direction. Since the path of the green photon trace is longer in the ice-rich outer disk (a + b) than the direct path (c), the traced ice column densities increase with radial separation of the observation. Right: direct lines of sight (blue) toward the star and warm inner disk are blocked by the optically thick disk midplane (white region). The observed light follows therefore a specific path through the optically thin disk regions, including one or multiple scattering events. The flux in a JWST spaxel is composed of photons with different paths (orange and yellow), that have not necessarily crossed equal lengths in the cold parts of the disk. This means that photon paths that cross large parts of the ice-rich regions of the disk can be saturated locally, even though the optical depth of the combined contributions is a defined, small number (see inset). The optical depth increases toward the midplane, because the light paths through the cold disk layers are more likely compared to paths higher up in the disk (green). Since the scattered light from the lower surface is tilted toward the cold midplane, and similarly, light from the upper surface is tilted away from the cold midplane, the profile is asymmetric with deeper ice features in the lower surface. Inset: qualitative schematic of the different flux contributions in the spectrum. The spectrum always has a contribution from light that passes through little ice (ambient light continuum; orange line) and light that passes through the cold regions of the disk, usually requiring multiple scattering events along the way (yellow line). The ratio of the contributions determines the optical depth in the observations if the feature is saturated. Photons with a single scattering event are less likely to reach regions with ice in the disk than photons with multiple scattering events (see also Fig. 3 in Sturm et al. 2023c).

In the text

	Fig. 6 Spatial variations in the CO ice peak optical depth (colorscale) compared to ALMA CO J = 2−1 observations presented in Sturm et al. (2023b; contours). The beam of the ALMA observations (0/51 × 0.″31) is shown in the bottom left corner. The JWST/NIRSpec observations have an effective resolution of 28 au or 1.5 spaxels.
In the text

Fig. A.1 Comparison between the spatial distribution of the CO emission and the continuum emission. Left: spatial distribution of the median continuum flux in the region with CO gas emission (4.36–5.2 µm). Middle: spatial distribution of the median of the peak fluxes of the CO gas emission lines. Right: ratio of the left two panels showing a constant ratio in the upper surface, but a higher ratio in the lower surface due to missing CO flux.

In the text

Fig. A.2 Overview of the CO gas line fitting procedure following the steps described in the text. Top panel: observed source integrated spectrum of HH 48 NE with selected troughs in between the emission lines in purple. Regions used to fit the continuum underneath the lines are marked in gray. Middle panel: continuum subtracted CO emission spectrum, corrected for the optical depth in the CO ice feature (blue region). The best fitting model is shown in red, regions taken into account in the fitting are marked with a gray horizontal line. Bottom panel: residual spectrum of subtracting the model from the data. The standard deviation over this region is 0.05 mJy, compared to 0.27 mJy in the CO line-free continuum at ~4 µm (see Fig. 2).

In the text