JWST Imaging of Edge-on Protoplanetary Disks. I. Fully Vertically Mixed 10 μm Grains in the Outer Regions of a 1000 au Disk

Scattered light imaging of protoplanetary disks provides key insights on the geometry and dust properties in the disk surface. Here, we present James Webb Space Telescope (JWST) 2–21 μm images of a 1000 au radius edge-on protoplanetary disk surrounding an 0.4 M ⊙ young star in Taurus, Two Micron All Sky Survey (2MASS) J04202144 + 2813491. These observations represent the longest wavelengths at which a protoplanetary disk is spatially resolved in scattered light. We combine these observations with Hubble Space Telescope optical images and Atacama Large Millimeter/submillimeter Array continuum and CO mapping. We find that the changes in the scattered light disk morphology are remarkably small across a factor of 30 in wavelength, indicating that dust in the disk surface layers is characterized by an almost gray opacity law. Using radiative transfer models, we conclude that grains up to ≳10 μm in size are fully coupled to the gas in this system, whereas grains ≳100 μm are strongly settled toward the midplane. Further analyses of these observations, and similar ones of other edge-on disks, will provide strong empirical constraints on disk dynamics and evolution and grain growth models. In addition, the 7.7 and 12. μm JWST images reveal an X-shaped feature located above the warm molecular layer traced by CO line emission. The highest elevations at which this feature is detectable roughly match the maximal extent of the disk in visible wavelength scattered light as well as of an unusual kinematic signature in CO. We propose that these phenomena could be related to a disk wind entraining small dust grains.


INTRODUCTION
Gas-rich protoplanetary disks are the birth site of planetary systems.How the small (≲ 1 µm) grains these disks inherit from their parent molecular cloud evolve to form km-sized planetesimals and, ultimately, rocky planets and cores of giant planets is hotly debated and multiple scenarios have been proposed (see Drazkowska et al. 2022, for a recent review).As a first step, small grains are generally believed to grow to mm-or cm-sized pebbles by coagulation through gentle hit-and-stick interactions (Blum & Wurm 2008).Spatially unresolved (sub)mm observations of disks surrounding T Tauri stars indicate that this process is largely underway by an age of ∼ 1 Myr (Williams & Cieza 2011).Once they have grown to large enough sizes, drag forces induced by the surrounding gas lead to their settling in the midplane and inward drift.The details of these two phenomena depend sensitively on the gas dynamics, which can be both turbulent and subject to various instabilities (Lesur et al. 2022).Mapping out the degree of settling and radial drift thus informs the physics of the initial stages of planet formation and disk evolution.By modifying the local dust opacity in different regions of the disk, settling has a significant influence on even unresolved observations of protoplanetary disks (e.g., Dullemond & Dominik 2004;D'Alessio et al. 2006) but these effects are generally ambiguous and dependent on assumptions about the disk structure.
The advent of high resolution imaging capabilities over the past decade has provided crucial insight on dust settling.In particular, the Atacama Large Millimeter/submillimeter Array (ALMA) observations of disks have shown that the vertical extent of the (sub)mmemitting pebbles is roughly an order of magnitude smaller than that of the gas and/or µm-sized dust that dominates scattered light (e.g., Pinte et al. 2016;Villenave et al. 2022).This is most readily illustrated in the case of edge-on disks (hereafter EODs), whose vertical extent is directly exposed to the observer.From a practical standpoint, these systems present the advantage that the disk blocks direct starlight, thus alleviating the need for coronagraphic techniques that are otherwise necessary.They are however severely underluminous as a result and, in general, cannot be observed with ground-based adaptive optics instruments.Spacebased observations, with the Hubble Space Telescope (HST) and James Webb Space Telescope (JWST), are therefore necessary to image them at the required resolution.Achieving similarly high angular resolution from the optical to the millimeter regime is key to constraining settling in disks associated with T Tauri stars.Indeed, direct comparisons of HST and millimeter observations of EODs have confirmed that settling is prevalent (Duchêne et al. 2003;Villenave et al. 2020;Wolff et al. 2021).
However, because those observational approaches are primarily sensitive to grains sizes that are several orders of magnitudes apart, they only provide constraints in two dramatically different regimes: µm-sized (and smaller) grains are fully coupled with the gas, whereas mm-sized pebbles are strongly decoupled.These are necessary constraints but it is critical to determine the behavior of intermediate-sized grains (several to a few tens of µm), which are only partially coupled to the gas, to fully disambiguate between disk and dust evolution models.Observations in the mid-infrared regime are crucial to fill this gap, as they are sensitive to such intermediate grains.In addition, thanks to the expected reduced opacity of dust at mid-infrared wavelengths, such observations should naturally probe deeper in the disk than the surface layers that are imaged with HST.Unresolved analyses of a few EODs have suggested that grains up to ∼ 10 µm are present in the upper scattering layers (Pontoppidan et al. 2007;Sturm et al. 2023a).This has been corroborated for the few EODs that have been imaged in the 3-5 µm range to date (Duchêne et al. 2010;Tobin et al. 2010;McCabe et al. 2011), although only two have been resolved at 10 µm (McCabe et al. 2003;Perrin et al. 2006).The successful deployment of JWST now provides a new opportunity to image EODs in scattered light through the near-and mid-infrared ranges.
This paper focuses on the disk surrounding 2MASS J04202144+2813491, hereafter Tau 042021.The object is located on the Western edge of the L1495 complex in the Taurus star forming-region, at a distance of 130 pc (Galli et al. 2019).Seeing-limited images first identified the system as an EOD, with two parallel nebulae representing the disk surfaces where stellar photons are scattered towards the observer (Luhman et al. 2009).Remarkably Tau 042021 is one of the largest disks known in Taurus, 50% larger than the iconic EOD HH 30 (Burrows et al. 1996).Deeper, higher-resolution HST images subsequently revealed a smooth and mostly symmetrical disk with a radius of at least 400 au, as well as a bright collimated bipolar atomic jet (Duchêne et al. 2014).This indicates that the central star is actively accreting, consistent with the detection of an apparent UV excess (see Appendix A).Based on near-infrared spectroscopy, Luhman et al. (2009) estimated a spectral type of M1±2, suggesting an 0.4-0.8M ⊙ for the central star.Andrews et al. (2013) reported the first millimeter observation of the system, with a flux density that places it in the top quartile among other Taurus members with a similar spectral type, suggesting a relatively massive disk, consistent with its large physical size.High-resolution 0.9-1.3mm ALMA observations were subsequently obtained by Villenave et al. (2020).These marginally resolved the disk along the vertical axis, with an apparent vertical size that is several times thinner than the separation between the two disk surfaces revealed in scattered light images.Broadly speaking, the ALMA maps support the presence of strong settling in the Tau 042021 disk.Villenave et al. (2020) further estimate a lower limit to the disk dust mass of 8 × 10 −5 M ⊙ although optical depth effects could boost this mass by up to an order of magnitude.
This paper presents the first results of an ongoing near-to mid-infrared imaging campaign of most known EODs with JWST (GO programs 2562 and 4290 in Cycles 1 and 2, co-PIs: F. Ménard and K. R. Stapelfeldt).Specifically, we present the first 2-21 µm JWST observations of Tau 042021 to further constrain the degree of settling in this disk.We also present CO ALMA maps to locate the warm molecular layer of the disk and to estimate the dynamical mass of the central star.Furthermore, we make use of archival GALEX, Spitzer and Herschel observations to complete the SED of the system and to probe the presence of mid-infrared emission lines.The paper is organized as follows: Section 2 summarizes the new and archival observations used here, and the corresponding observational results are presented in Section 3. In Section 4, we construct several toy models, with varying degrees of settling and dust evolution to try and match the appearance of Tau 042021 across the wavelength range.We then discuss the implications of our findings in Section 5.

JWST
We observed Tau 042021 with JWST as part of GO program 2562, using NIRCam and MIRI in two consecutive visits starting at 2023 January 23 UT 09:00.See Table 1 for more detail including total exposure times per filter.We used the dual channel capability of NIR-Cam to obtain simultaneous F200W and F444W images.With MIRI, we obtained consecutive images with the F770W, F1280W and F2100W filters.With both instruments, we used 4-point dither patterns to improve spatial sampling and mitigate detector artifacts: NIRCam used module B with the STANDARD subpixel dithering (because this target is spatially small, we did not use any large primary dither to fill in the SW inter-chip gaps), while MIRI used the BRIGHTSKY subarray and EXTENDED SOURCE dither pattern.Exposure settings (e.g.readout pattern and number of groups per integration) were chosen to avoid saturation on target, and the number of integrations per position were set to achieve the desired total integration times for high SNR on fainter surrounding nebulosity, based on ETC predictions scaled from the HST images.
We obtained the level-3 reduced and combined i2d data products from the Mikulski Archive for Space Telescopes (MAST).The retrieved files were reduced with pipeline version 1.8.2 (Bushouse et al. 2022) and using CRDS reference file context 1041.These frames are corrected for flat-field, cosmic rays and other outliers, and optical distortion, but not for background.The background in the dither-combined NIRCam images is low and uniform and we simply subtracted the median value across the frame to remove it.This approach fails for the MIRI images due to residual gradients in the background.We therefore retrieved the individual level-2 data products from MAST, generated a background frame by median-combining the exposures at each individual location of the dither patterns, subtracted it from each frame, and manually aligned and median-combined the resulting images to produce the final images of Tau 042021.Finally, we used a 5 ′′ square aperture in all JWST images to measure fluxes of 2.4, 2.6, 2.3, 1.8 and 1.9 mJy in the F200W, F444W, F770W, F1280W and F2100W filters, respectively, with a conservative 10% uncertainty based on comparing these fluxes with those in larger apertures.We previously observed Tau 042021 with HST/ACS as part of Program 12514 (PI: K. R. Stapelfeldt).The observations were first presented in Duchêne et al. (2014) and Stapelfeldt et al. (2014).We obtained consecutive images with the F606W and F814W filters in a single orbit using CRSPLIT=2 to reject cosmic rays.The images presented here are the drz data products available through MAST, which are corrected for flat field and optical distortion.The background level is low and uniform in both filters, and we subtracted the median value over the entire images to correct for it.

ALMA
Continuum maps in ALMA bands 4, 6 and 7 were previously presented in Villenave et al. (2020), to which we refer the reader for details about observations and data reduction.Here we report the CO J = 2 − 1 observations associated with the already published band 6 continuum map (program 2016.1.0771.S, PI: G. Duchêne).The source was observed using a compact array configuration on December 5, 2016 and an extended configuration on October 21, 2016.The spectral set-up was divided into two 1.9 GHz continuum spectral windows, of rest frequency 216.5 GHz and 219.6 GHz, and three spectral windows that were set-up to include the 12 CO, 13 CO, and C 18 O J = 2 -1 transitions at 230.538 GHz, 220.399GHz and 219.56 GHz, respectively.Here, we focus on the 12 CO and 13 CO observations, which were obtained with a sampling of 122 kHz and 61 kHz, respectively.The raw data were calibrated using the CASA pipeline version 4.7.2 (CASA Team et al. 2022).
To maximize the dynamical range, we performed phase self-calibration on the continuum data for the compact configuration, and applied the continuum solution to the 12 CO and 13 CO spectral windows.We extracted the emission lines from the calibrated visibilities after subtracting continuum emission using the uvcontsub function in CASA.We used the tclean task, with a spectral resolution of 0.22 km s −1 and a Briggs robust parameter of 0.5, to produce two sets of line emission maps.For the maps combining both extended and compact configuration observations, we obtain an angular resolution of 0. ′′ 19×0.′′ 37.In addition, we also use a uvtaper to produce a smoothed image of the 12 CO channel map centered at V − V sys = -0.3km s −1 , which achieved an angular resolution of 0. ′′ 74×0.′′ 67.HST-JWST 3-color composite image of Tau 042021; the full field of view spans 18 ′′ .The 0.8, 2 and 7.7 µm images are rendered (using a logarithmic stretch) in the blue, green and red channels, respectively.The orange object to the right of the edge-on disk is a background galaxy.All images from the optical to the mid-infrared reveal a similar morphology, namely two nearly symmetrical nebulae bisected by a dark lane, as expected for EODs.The system's SED has its trough around 20 µm (see Figure 3), and the 21 µm image confirms that the disk is still optically thick to its own radiation out to that wavelength.In the ALMA regime, the disk is either optically thin (and emission is dominated by the denser midplane) or the two scattering surfaces are too close to each other to be fully resolved (possibly due to settling).Villenave et al. (2020) showed that the intrinsic vertical extent is about 0. ′′ 2-0.′′ 3, or 25-45 au, with a tentative identification of separation between layers, favoring the second scenario (see below).
The overall morphology of the disk does not qualitatively change out to 5 µm with both top and bottom nebulae showing increasing curvature as a function of distance from the symmetry axis.In the optical, faint, diffuse scattered light is detected at high elevations (up to ≈ 2 ′′ from the midplane), unlike what is seen in the 2-5 µm regime, despite comparable (or better) dynamic range in the JWST/NIRCam images.The MIRI 7.7 µm and 12.8 µm images reveal a markedly different structure, although the bisecting dark lane is still clearly present.On both sides of the disk a strong symmetrical X-shaped feature is present on top of a structure that is similar to the 4.4 µm image.This will be further explored in Section 3.3.At the longest JWST wavelength, the disk is back to a simple two-nebulae structure, with a significantly reduced radial extent, and a much more centrally condensed peak, which we quantify below.
Of interest to this study is the gradual decline in the thickness of the dark lane, measured as the distance between the peaks of the top and bottom nebulae along the disk's symmetry axis (see Figure 4).To quantify this, we followed a similar method to that presented in Villenave et al. (2020).We first generated vertical brightness profiles, averaging within 0. ′′ 4-0.′′ 5 windows to improve SNR, and fitted sixth to eighth order polynomials to the peak of each nebulae.We repeated the process out to 1. ′′ 5-3 ′′ (depending on filter) on either side of the symmetry axis and then fitted second order polynomials to the spine of each nebulae to evaluate the closest distance between the resulting parabolas, d neb , as well as their peak flux ratio, FWHM and full width at 10% of the peak (FW10%).The latter is useful as an indication of the disk radius.Varying the polynomial degrees used in the method, we estimate uncertainties of about 0. ′′ 02 on d neb .This two-step approach is more reliable than simply using vertical brightness profile along the symmetry axis for wavelengths with significant contribution from the collimated jet.As indicated in Table 2 and shown in Figure 4, the chromaticity of the dark lane thickness is much shallower than for other EODs (Watson & Stapelfeldt 2004;Duchêne et al. 2010;McCabe et al. 2011).We defer discussion of these results to Section 5.1.
Since Tau 042021 is also resolved along the vertical direction in the ALMA observations (Villenave et al. 2020, see also Figure 4), we extracted the profile within 0. ′′ 5 of Closed and open circles represent the adopted and additional photometric datapoints, respectively (see Appendix A).Red diamonds indicate the fluxes measured in the JWST images presented here.Small gray circles represent the Spitzer/IRS spectrum rebinned to a ∆λ/λ = 10% sampling.The observed scatter observed from the optical to the mid-infrared is most likely due to intrinsic variability and/or changes in disk illumination.The dashed blue and red segments are power law fits to the UV/blue and (sub)mm parts of the SED, respectively.
the symmetry axis and noticed that it presents a shoulder characteristic of two distinct components separated by a distance comparable to the projected beam size.We interpret this as an indication that the disk is still optically thick at 890 µm, but with the two disk surfaces not fully separated.We thus fitted a double Gaussian model to this profile and estimated a distance between the two surfaces of 0. ′′ 11±0.′′ 02, roughly ten times thinner than is observed in the optical through the mid-infrared.
In the HST images, the disk presents a significant lateral asymmetry in its outer regions (away from the axis of symmetry), with the top and bottom nebulae extending further in opposite directions from one another.In contrast, in the JWST images, such an asymmetry is strongly suppressed across all wavelengths.There is therefore temporal variability in the illumination from the central star, in line with the findings from Luhman et al. (2009).Unfortunately, there is currently not enough temporal coverage to determine whether this variability is periodic.
The brightness ratio between top and bottom nebula in an EOD is related to the system's inclination for symmetric disk, via the dust scattering phase function and extinction through the disk (Watson et al. 2007).The brightness ratio between the two nebulae is ≈ 0.5 (see Table 2), suggesting a small deviation from the i = 90    Note-The distance between the apexes of the top and bottom nebulae, d neb , is estimated based on fitting a parabola to each spine (see Section 3.1).Their flux ratio is measured along the vertical profile at the location corresponding to this smallest separation.The FWHM of both nebulae and full width at 10% of the top nebula are estimated from the brightness profile along their respective spines.† Dataset with significant contamination by line emission along the collimated jet.
# Dataset with significant contamination by the X feature.
configuration.We note, however, that this ratio has been found to vary and occasionally approach unity (Luhman et al. 2009), which precludes a definitive estimate.In addition, we find a mild increase from the optical to the mid-infrared, from ≈ 0.45 to ≈ 0.6.By comparison, this flux ratio was found to strongly increase and decrease in the cases of HK Tau B and HV Tau C, respectively (Duchêne et al. 2010;McCabe et al. 2011).
The surface brightness profile along the spine of each nebulae is also (indirectly) driven by the scattering phase function since it probes scattering angles that range from ≈ 0 • to ≈ 90 • (along the major axis and at the tips of the each nebula, respectively).As indicated in Table 2, both the top and bottom nebulae have essentially achromatic FWHM, with the fainter (bottom) nebula being broader than the brighter (top) one.This results from an enhanced role of multiple scat-tering for photons escaping through the bottom layer.This is in line with the behavior observed for HV Tau C (Duchêne et al. 2010).Instead, in HK Tau B, the images become increasingly peakier at longer wavelengths (McCabe et al. 2011), suggesting that there is diversity in the dust scattering properties of EODs.In the mid-infrared, the Tau 042021 disk appears to be somewhat more compact as measured by FW10% top , with a decrease of 20-30% relative to the optical and nearinfrared.
Finally, we note the presence of a background galaxy (clearly identified as such in the JWST 2 µm image, see Figure 1) at a projected distance of 6. ′′ 6, slightly interior to the gaseous disk outer edge (see Section 3.4).The galaxy is faintly detected in the HST F814W image but not in the F606W.We thus conclude that the disk must be optically thin at least at λ ≳ 0.8 µm at this location, revealing a low column density of dust in these outer regions.Furthermore, to estimate the proper motion of Tau 042021 between the HST and JWST observing epochs, we evaluated the centroid position of the galaxy in the 0.8 µm and 2 µm images, respectively, and measured the offset from the location of peak surface brightness in the disk (center of the brighter nebula).This yields a proper motion for Tau 042021 of about 28 mas/yr South and 11 mas/yr East, with an uncertainty of about 4 mas/yr based on an estimated ≈1 pixel uncertainty in the disk position.This result is consistent with membership to the L1495 cloud (Galli et al. 2019).Interestingly, the projected separation between Tau 042021 and the background galaxy is declining so that photometric monitoring of the galaxy could be used as a probe of dust extinction in the disk over the next few decades.

The collimated jet
Jet emission in the form of discrete knots is clearly detected in the HST and some of the JWST images.To highlight these features, we applied a 1 ′′ box median filter to these images; the resulting high-passed images are shown in Figure 5.A series of prominent emission knots can be traced up to ≈ 6 ′′ in the HST 0.6 µm image.Jet emission in this filter is typically produced by a combination of Hα and several nearby forbidden lines, in particular [O I], [N II] and [S II] (Ray et al. 2007).Fainter emission is also identified in the same knots in the 0.8 µm image.We speculate that these are due to [O II], [S III] or [Fe II] forbidden lines, which have been seen in some jet-driving T Tauri stars with near-edge-on viewing geometry (Bacciotti et al. 2011;Whelan et al. 2014).Several knots are similarly identified in the JWST 2 µm and 12.8 µm images.In the infrared, H 2 and [Ne II] are among the dominant outflow-related emission line (Lahuis et al. 2007).To confirm that these lines are the likely tracers of the jet knots in the Tau 042021 system, we obtained the Spitzer/IRS spectrum of Tau 042021 from the Combined Atlas of Sources with Spitzer IRS Spectra (CASSIS, 1 Lebouteiller et al. 2011).Although the source is faint (≈ 1 mJy), the detection is solid and the spectrum is of good quality across the entire 5-35 µm range.Of particular interest for this study is the λ ≤ 15 µm regime, which contains the F770W and F1280W bandpasses in which the X-shaped feature is detected (see Section 3.3).As shown in Figure 6, there are multiple emission lines present in the spectrum, with the most prominent being the 12.8 µm [Ne II] line and the H 2 S(1)-S( 7) series (some of these lines were already reported in Sturm et al. 2023b).The former is commonly associated with disk winds and/or outflows in young stars (e.g., Güdel et al. 2010;Pascucci et al. 2020) whereas the latter is generally thought to be due to UV excitation in the uppermost layers of disks and in photodissociative winds (e.g., Nomura & Millar 2005;Lahuis et al. 2007).We also find tentative detections of the H I 6-5, 8-6, 9-7 and 7-6 lines (7.46, 7.50, 11.31 and 12.37 µm, respectively), which are thought to be driven by accretion on the central star (Rigliaco et al. 2015).

A new X-shaped feature
The 7.7 µm and 12.8 µm images reveal a prominent X-shaped feature with an essentially linear shape that is absent in all other images of the system.Such straight features are often a consequence of telescope diffraction and/or detector properties, and JWST/MIRI is known for displaying such structures in some filters (e.g., Gáspár et al. 2021).To assess whether this feature is instrumental or astrophysical, we compare the F770W image of Tau 042021 with that of a bright point source (2MASS J04202192+2813206) observed in the same frame in Figure 7. None of the diffraction features of the hexagonal primary mirror line up with the X-shaped feature, nor does the so-called "cross artefact" that runs along the detector columns and rows.The Xshaped feature associated with Tau 042021 is therefore astrophysical in nature.Given the bandpasses of the F770W and F1280W filters, we interpret the emission in this feature as due to photodissociated H 2 , shock-excited [Ne II] and/or thermal emission from out-of-equilibrium very small grains, either at the disk surface or at the base of disk wind.This can be confirmed with integral field unit JWST observations as these have different spectral characteristics.
Figure 8 shows the 7.7 µm image of Tau 042021 with a stretch that highlights the faint outer substructures.In addition to the X-shaped feature, fainter whiskers are seen extending at similar position angles as the disk surface traced in the 4.4 µm image.To better characterize the location of this feature, we also compare the 7.7 µm image to the molecular warm layer, as determined by the moment zero map of the 12 CO emission (see Section 3.4).It appears that the warm molecular layer is located between the 5-8 µm scattered light disk "surface" and the X-shaped feature itself.
The "arms" of the X-shaped feature appear remarkably linear and are well represented by symmetrical straight lines inclined ≈ 36 • from the disk midplane.The location where the two lines intersect aligns well with the jet features seen in the 12.8 µm image (see Section 3.2) and we therefore assume that it is the location of the central star.We further note that the arms of the feature are resolved along their minor axis, with a FWHM of ∼0.′′ 6 in the 7.7 µm image, i.e., more than twice the instrumental resolution.It is likely that the Xshaped feature is the limb-brightened edge of a biconical structure, whose nature will be discussed in Section 5.2.This feature can be traced up to an elevation of 1. ′′ 5-1.′′ 75 on either side of the midplane, almost as high as the most vertically extended scattered light detected in the HST optical images.

The gaseous disk
The ALMA spectral cubes show a strong detection of the 12 CO and 13 CO emission, which enables us to 1) estimate the dynamical mass of the central star via its (near) Keplerian motion, and 2) probe the full radial extent of the disk.Inspection of the spectral cubes reveals that the northern side of the disk is rotating away from the observer.By comparing the blue-and red-shifted channels, we estimate that V sys = 7.4 km s −1 consistent with the ≈ 7.0 ± 0.5 km s −1 range observed in the L1495 cloud (Galli et al. 2019).This is further confirmation that Tau 042021 is physically associated with that cloud within the Taurus complex.We also note that the 12 CO disk emission is heavily suppressed in the V LSR = 7-8 km s −1 range, consistent with cloud emission (Gold-smith et al. 2008) that is filtered by the interferometer.Apart from this, the moment maps (Figure 9) show the typical morphology of EODs (Dutrey et al. 2017;Louvet et al. 2018), with the two warm molecular layers separated by a cold midplane where CO freezes out on dust grains.
We then constructed the 12 CO and 13 CO positionvelocity diagrams and superimposed representative Keplerian curves (Figure 10).A coarse exploration suggests a stellar mass in the 0.3-0.4M ⊙ range.We note that this estimate is consistent with the lower end of the range derived from the system's spectral type while the higher end of that range is firmly excluded by observations (Luhman et al. 2009).Conversely, the 0.25 M ⊙ estimated by Simon et al. (2019, after scaling to the 130 pc distance used here) is too low.We note that the best-fitting Keplerian model derived by Simon et al. (2019) present significant residuals, which may be due to a systematic bias at high inclination induced by the methodology employed by these authors or to the presence of non-Keplerian motion (see below).While obtaining a precise estimate of the stellar mass is beyond the scope of this study since it does not directly affect our analysis, we conclude that M ⋆ ≈ 0.4 M ⊙ .The moment maps are well defined up to a radius of about 2 ′′ , significantly less than the HST or JWST radial extent but marginally larger than the submillimeter continuum (Villenave et al. 2020).However, faint 12 CO emission is detected out to much larger distances and is best mapped in a few individual channels.On the blueshifted (South) side of the disk, the channel centered at V −V sys = −0.3km s −1 reveals detectable emission in the disk upper layers out to a radius of at least ≈ 8 ′′ (≈ 1000 au, see Figure 11).The fact that the CO emission is most extended at such a low velocity (relative to the systemic velocity) is highly unusual among EODs (e.g., Guilloteau et al. 2016;Louvet et al. 2018;Teague et al. 2020) and suggests non-Keplerian motion.Unfortunately, the corresponding red-shifted channel is devoid of emission due to filtered-out cloud emission.Nonetheless, this measured size is twice as much as in the HST and JWST images, and ≈ 3.5 times more than in the submillimeter continuum emission.This is in line with the typical ratio between gas and dust (submillimeter) component of protoplanetary disks, both edge-on and at lower inclinations (Ansdell et al. 2018;Flores et al. 2021).This channel map shows a notably boxy pattern, essentially flat-topped, suggesting that the warm molecular layer extends over a wide range of radii at a fixed elevation of ∼ 325 au, a height that matches the end points of the X-shaped feature discussed in Section 3.3.

Objectives and modeling setup
With the data presented here, the Tau 042021 disk is now extremely well characterized, with scattered light images spanning a factor of 30 in wavelengths, several thermal emission maps in the (sub)millimeter regime, CO (2-1) spectral cubes, and a well sampled SED that includes a Spitzer/IRS spectrum.While simultaneously modeling all of these observables is the principled methodology to characterizing the disk structure in great detail, it is beyond the scope of this study.Instead, we adopt an exploratory approach that focuses on two key physical aspects, the grain size distribution (indicative of grain growth) and the degree of settling in the disk.To probe these phenomena in an efficient manner, we consider the thickness of the dark lane that separate the two scattered light surfaces (d neb ) as a direct, albeit reductive, observable tracer of dust opacity and disk vertical structure.All model calculations are performed with the MCFOST radiative transfer code (Pinte et al. 2006(Pinte et al. , 2009)).
Concretely, we adopted the following methodology.We first fix the stellar properties based on our current knowledge of the system.Specifically, we adopt a stellar mass of 0.4 M ⊙ (Section 3.4) and an effective temperature of 3700 K based on the system's spectral type and adopting the conversion law from Luhman et al. (2003).We further set the stellar luminosity to L ⋆ = 0.85 L ⊙ based on the median luminosity of Taurus members with the same spectral type (Luhman et al. 2009).Our model also includes a UV excess component that is evident in the short wavelength end of the system's SED.Fitting a power law to the datapoints at λ < 0.5 µm, we find that this excess is characterized by a F ν ∝ ν −2 shape and a total excess of F UV /F ⋆ = 4% (see Figure 3).Finally, we include no foreground extinction, consistent with the GALEX detection of Tau 042021.
To describe the density structure of the gaseous component, we adopt the commonly used tapered edge surface density profile (see, e.g., Williams & Cieza 2011;Andrews 2015): where we set γ = γ exp = 1 and r c = 300 au, coupled with a hard outer cutoff at r out = 400 au.These values were selected based on an initial manual exploration of the parameter space as they produced an acceptable extent along the disk major axis.We assumed a Gaussian vertical density profile, appropriate for a vertically isothermal disk, characterized by a scale height of h 0 = 11 au at r 0 = 100 au and a flaring exponent of β = 1.125 so that h(r) = h 0 ( r r0 ) β .The value for h 0 was computed from Eq. 6 in Andrews (2015), using the midplane temperature of the same initial models.The vertically isothermal assumption is incorrect as the disk surface is superheated but the resulting difference in vertical density profile is most noticeable in the optically thin highest regions of the disk that do not contribute significantly to scattered light images.Besides, the assumption is only employed to set the vertical density structure; the temperature is subsequently computed locally during the radiative transfer calculation.
Finally, we selected an inclination of i = 88 • , which generally leads to peak flux ratios of ≈ 2 between the two disk surfaces, consistent with the observations presented here.While model scattered light images are dependent on all of these parameters, our analysis on the wavelength dependence of the dark lane thickness is not strongly affected as it is a chromatic effect.Thus, fixing these parameters reduces dramatically the dimensionality of the problem while keeping the focus on the topic of interest here, namely the dust properties and especially their spatially-dependent size distribution.Nonetheless, a more thorough exploration of the disk structural parameters will be required in future modeling of the source.As described below, we explore models in which dust particles of all sizes are fully mixed with the gas, as well as some models in which larger grains are settled closer to the midplane.In the absence of direct quantitative constraints on the dust composition in the Tau 042021 disk, we adopted the DIANA mixture, namely a 60/15/25 (volume fraction) mixture of silicates, amorphous carbon and vacuum (Woitke et al. 2016).Instead of the standard Mie theory, applicable for compact spherical grains, this model uses the Distribution of Hollow Spheres method (DHS, Min et al. 2016), which provides absorption and scattering properties that mimic aggregate dust particles.In combination with the selected porosity, this treatment of dust scattering produces stronger forward scattering.In turn, this leads to more centrally peaked images of EODs (since the central part of the top disk surface is seen at very small scattering angles), which matches better with our observations of Tau 042021.
Regarding the grain size distribution, we first consider a "standard" model that is a dn(a) ∝ a −p da power law with p = 3.5 ranging from a min = 0.03 µm to a max = 1 mm, although we will also explore devia-tions from these default values (see Table 3).Furthermore, we consider a "pristine" model, in which we built a grain size distribution that produces an extinction law consistent with that observed in molecular clouds out to λ = 70 µm, albeit with moderately too strong midinfrared silicate features (see Figure 12).We refer the reader to Appendix B for the details of this model.We stress that this model should not be considered as physically realistic but it achieves the desired opacity law, which is the quantity we wish to probe in the Tau 042021 disk.
Since empirical evidence for dust settling of large (millimeter-emitting) dust in protoplanetary disks is now strong (e.g., Villenave et al. 2020), we incorporate a simple settling prescription to explore its qualitative influence on the JWST images of Tau 042021.Rather than adopting a specific physical model, which depends on poorly constrained ingredients such as the level of turbulence in a disk and the microphysics of grain-gas interactions, we adopt a simple analytical prescription, whereby grains with a > a mix have a reduced scale height: h(a) ∝ (a/a mix ) η settl , with η settl < 0. Expectations from hydrodynamical and magnetohydrodynamical turbulent settling suggest η settl ≳ −0.5 (e.g., Dubrulle et al. 1995;Fromang & Nelson 2009), while ALMA observations indicate that h(a ≳ 100µm) ≈ 1 au (Pinte et al. 2016;Villenave et al. 2020).
For each model (i.e., different size distribution and settling parameters), we then adjusted the disk mass to match the observed dark lane thickness at 2 µm.We selected the JWST 2 µm disk image for this as it has the best contrast/resolution combination in our datasets.Having set this last parameter, we then compute images at other wavelengths from the optical to the submillimeter as well as the complete SED, and determine whether it is an acceptable fit to the observations (see Table 3).

Modeling results
Figure 13 presents a gallery of PSF-convolved model images from the optical to the sub-mm.Furthermore, Figure 14 presents the wavelength dependence of the distance between the top and bottom nebulae (referred to as "dark lane thickness" in Table 3), measured in the same manner as in the data, as well as the observed and model SED.
Focusing first on the dark lane thickness metric, the standard model displays a decline of about 25% from the optical to the mid-infrared, similar to, albeit slightly steeper, than in the observations.The disk morphology is reasonably well reproduced, although the images are not quite as centrally peaked as in the observations, Note-The symbols in the last three columns describe the degree to which a given model matches observations: a ✓ symbol represents an acceptable quantitative match, a (✓) indicates a small but significant deviation, while a × symbol points to a large shortcoming of the model.
suggesting that the scattering phase function is not sufficiently forward scattering.The overall SED of the system is also satisfyingly well reproduced, although we note a strong albedo feature around 9 µm that is not seen in the Spitzer observations but is inherent to the dust composition adopted here.The submillimeter regime is where this model shows its largest shortcoming, with the disk much too vertically extended and showing too flared of a structure at large distances from the symmetry axis.This highlights that settling is most likely present in this system.Nonetheless, considering that the model disk mass is only set by the distance between nebulae at 2 µm in our modelling approach, we conclude that this model overall matches surprisingly well with all observations.Conversely, the "pristine dust" model experiences a strongly chromatic dark lane thickness, with a ≳ 50% decline from the optical to the mid-infrared, in line with its much steeper opacity law (Figure 12).In addition, this model underpredicts the SED by more than an order of magnitude in the 2-20 µm range, for lack of effective scatterers in that wavelength regime.Finally, in the submillimeter regime, the model image is much too radially concentrated, as a consequence of the strongly reduced opacity of that dust model at the longest wavelengths.All in all, this model convincingly demonstrates that dust in the Tau 042021 disk is significantly different from that present in molecular cloud.
To try and reconcile the very shallow chromaticity of the dark lane with a fully mixed model, we considered two minor changes to the standard model.In the first one, we adopted a slightly shallower grain size distribution (p = 3.25), while in the other one we increased the minimum grain size to a min = 0.3 µm.Since both of these models produce a flatter extinction law from the optical to the near-infrared (Figure 12), they succeed in matching the dark lane thickness behavior better than the standard model.However, they still produce submillimeter maps that do not match the observed one.Furthermore, the model with increased minimum grain size underpredicts the system's SED at λ ≲ 0.5 µm as all dust grains are too large to scatter such short wavelengths.These models indicate that a modest decrease in the population of small (≲ 0.1 µm) relative to the standard model may yield a satisfying match to the HST and JWST images, although it would still fail in the ALMA regime.
Finally, we considered three settling models in which grains up to a mix = 100 µm, 10 µm or 1 µm are fully mixed, with large grains increasingly settled with η settl = −0.5.These models yield an effective scale height of ≈ 3, 1 and 0.3 au at 100 au for 1 mm grains, respectively, in reasonable agreement with derived scale height in protoplanetary disks (Pinte et al. 2016;Villenave et al. 2022).The a mix = 100 µm model is essentially indistinguishable from the "standard" model in all observations, indicating that this it produces an insufficient degree of settling for mm-emitting grains.For smaller values of a mix , settled models produce much improved 890 µm images, as expected, although the surface brightness along the major axis is not quite flat enough.However, whereas the a mix = 10 µm model shows a dark lane thickness chromaticity that is almost identical to the standard model presented above, i.e., only slightly too steep, the a mix = 1 µm model shows a significant decline at the longest mid-infrared wavelengths due to its significant settling of infrared-scattering grains.Furthermore, the SED for that model shows a deeper trough in the mid-infrared whereas the a mix = 10 µm model  3.All images are shown on a square root stretch except the last column (linear stretch) from 0 to the peak value.matches well with the observed SED, albeit with a modest overprediction in the 30-100 µm regime.
Overall, as summarized in Table 3, we consider this brief model exploration as a strong indication that a model in which grains up to ∼ 10 µm in size are fully mixed with the gas and larger grains (≳ 100 µm) are increasingly settled down to a scale height of ∼ 1 au could simultaneously match the HST, JWST and ALMA continuum observations of Tau 042021.We note that the disk dust masses required to match the dark lane thickness at 2 µm are fairly high, in the ∼ 3 − 10 × 10 −4 M ⊙ range.This is consistent with the lower limit derived independently by Villenave et al. (2020) from ALMA continuum maps.While the dust mass is inversely proportional to the assumed dust opacity and remains therefore somewhat uncertain, the Tau 042021 disk is on the high end of the observed distribution of M d /M ⋆ ratio (e.g., Barenfeld et al. 2016), consistent with it being one the  3).Right: Observed and model SED for the same models.The dip at 9 µm in the model SEDs is due to a strong albedo feature from the silicates adopted in the DIANA mixture.We stress that the models are constructed solely on the basis of matching the dark lane thickness at 2 µm; all other quantities are merely predicted based on having fit that quantity.
largest disks in Taurus (Andrews et al. 2018).We defer additional improvements upon this model, including a different scattering phase function and a more realistic settling prescription, to a further study.

Dust properties and settling
With a 1000 au radius as measured in CO, Tau 042021 is one of the largest EOD in nearby star-forming regions.
Its images across all wavelengths are smooth and largely symmetric and it has no known companion.It is therefore an ideal testbed for models of disk structure and dust evolution.The JWST images we presented here allow us to probe the properties of its dust component, mostly opacity and scattering phase function, as well as the degree of vertical settling in the disk.
In scattered light imaging, the thickness of the dark lane that separates the scattering surfaces is driven by the optical depth along our line of sight to the central star.In turn, this depends on the product of the disk mass and the dust opacity.The very modest drop in dark lane thickness from the optical through the midinfrared indicates that the dust located in the upper layers of the Tau 042021 disk is characterized by a (nearly) gray opacity law over that wavelength range.Indeed, our analysis showed that a model based on a dust population consistent with that found in molecular cloud is strongly inconsistent with observations.This confirms that significant grain growth occurs in protoplanetary disks.Assuming the standard power law size distribution (p = 3.5), the limited dark lane chromaticy in this disk indicates that grains at least as large as 10 µm are present in copious amount in its upper layers.Similar conclusions were reached for other disks based on a va-riety of observation constraints, such as system SED, polarized intensity scattered light imaging or water ice absorption (Pontoppidan et al. 2007;Tazaki et al. 2021a;Franceschi et al. 2023;Sturm et al. 2023a).The JWST images presented here provide a direct probe of the dust opacity in the upper outer layers of the Tau 042021 disk, which complements these other methods.Specifically, the observations presented here demonstrate the abundant presence of 10 µm-sized grains in a region located ∼ 60 au above the midplane at a radius of ∼ 400 au from the central star.Ultimately, combining these different observables, as well as the spectral shape of the midinfrared silicate features, will enable a more granular analysis of the dust properties in the disk upper layers, including allowing for departures from simple power law parametrizations.Indeed, each observational approach is uniquely sensitive to different aspects of the dust component, both spatially and in terms of grain size.
One of the shortcomings of our "standard" model, however, is that it predicts a change in dark lane thickness between the optical and near-infrared ranges that is 10-15% too large.An even shallower dust opacity law seems necessary to account for this, as supported by our models using a shallower size distribution or a larger minimum grain size.This is an indication that the growth of dust particles may introduce a break in the slope of the distribution, as is suggested by numerical models of grain growth and fragmentation (e.g., Birnstiel et al. 2011).However, it is unclear if grains as large as 10 µm can efficiently grow in situ in the uppermost layers of the outer disks, where overall densities are relatively low.Instead, a (turbulent) mechanism may be necessary to loft them up from deeper layers, where they could more easily form (e.g., Tazaki et al. 2021a;Franceschi et al. 2023).
While the fully mixed models we developed can (roughly) reproduce all HST and JWST images of Tau 042021, they produce submillimeter images that do not match the observations.In short, the optical depth through the disk is much too high, and thermally emitted photons scatter off layers that are a few tens of au above the midplane.In our models, we were able to obtain a better match to observations by introducing some settling of large grains, as suggested by other ALMA observations of protoplanetary disks.While we did not consider physical models of dust settling, our parametric approach allows us to reach some general conclusions.To reproduce the images of Tau 042021 through the HST-JWST regime, we find that grains up to ≳ 10 µm must be vertically fully mixed, i.e., turbulence is sufficiently strong to dredge them all the way to the disk surface.Conversely, millimeter-emitting grains (sizes ≳ 100µm) must have a scale height that is roughly one tenth that of the gas and the smallest dust grains.
Let us now consider how this compares to settling model predictions.Pure hydrodynamical simulations (e.g., Dubrulle et al. 1995) indicate that the scale height of the dust component scales as h dust ≈ h gas α α+St , where St is the Stokes number, i.e. the product of the orbital frequency by the stopping time of a grain as a result of the drag force induced by the surrounding gas, and α the standard turbulence parameter (Shakura & Sunyaev 1973).In the large grain limit, St ∝ a, leading to h dust ∝ a −0.5 in the case of weak turbulence (α ≪ St).Therefore, mm-sized dust are expected to have a scale height that is ≈ 10 times smaller than that of 10 µm grains, in line with our findings.Magnetohydrodynamical models, on the other hand, predict shallower dependencies: Fromang & Nelson (2009) find an h dust ∝ St −0.2 approximate relationship, for instance.This may imply that the importance of magnetized transport is reduced in the outer regions of the Tau 042021 disk, which are primarily traced by scattered light images.However, Riols & Lesur (2018) have shown that, while the hydrodynamical solutions may be approximately correct in the context of 1D models, higher dimension simulations reveal more complex behaviors, in which it is nearly impossible to characterize settling through a single scale height quantity.To make matters worse, both types of simulations concur on the fact that the relevant quantity is not directly the size of dust grains but their Stokes number.Because the latter depends also on the local gas density and the material density of the dust particles, the constraint we have obtained from the JWST images presented here can only be indirectly compared to such models.We leave the comparison of more physically-grounded settling models with JWST images to future studies.
While the models we have constructed reproduce the dark lane thickness, they are imperfect when considering the overall morphology of the Tau 042021 disk.In particular, model images appear consistently less centrally peaked than is observed, especially at the longest JWST wavelengths.Although this may be due to a combination of factors, we suspect that a main issue is that the dust prescription we have used leads to a scattering phase function that is not sufficiently forward scattering.The more peaky phase functions could be either due to solid grains with a lower refractive index (i.e., more ice rich) or porous aggregates.Indeed, lower refractive indices of solid grains facilitate efficient transmission of light, leading to more efficient forward scattering (e.g., Figures 10 and 17 in Tazaki et al. 2021b).Porous aggregates are also known as efficient forward scatterers (Min et al. 2012;Tazaki et al. 2019), although the phase function depends on the detailed structure of aggregates, such as fractal dimension (Min et al. 2016;Tazaki et al. 2016;Tobon Valencia et al. 2022).Future modeling of this system should explore such dust properties to improve the match between observed and model scattered light images of Tau 042021 (R. Tazaki, priv. comm.).
Finally, we note that the modest chromaticity of the dark lane thickness in Tau 042021 is in stark contrast with previously-imaged EODs, which become markedly thinner towards longer wavelengths, at least out to 5 µm (Duchêne et al. 2010;McCabe et al. 2011, see also Figure 4).One possible explanation is that this system contains a shallower size distribution than other disks, possibly as a result of more advanced dust evolution.It could also be the case that the Tau 042021 disk is more optically thick than other EODs.Indeed, in the limit of very high optical depths, the chromatic dependence of the opacity law is gradually muted with increasing disk mass as the elevation of τ = 1 surface does not depend much of the total column density in the disk since the density drops rapidly in the disk atmosphere.In other words, if the disk is very optically thick up to ≈ 60au above the midplane at λ < 21 µm and the density drops extremely fast above this layer, the lack of chromaticity could be explained with a wide range of grain sizes.The gradual decline in surface brightness in the optical HST images seems to argue against this scenario but further analyses are needed to definitely rule it out.JWST observations of other EODs will allow us to place Tau 042021 in the larger context.

A possible disk wind
Arguably the most surprising outcome of the JWST observations of Tau 042021 is the X-shaped feature that is prominent at 7.7 and 12.8 µm and resides immediately above the warm molecular layer traced by the ALMA 12 CO mapping, at elevations up to ∼ 225 au.Without follow-up spectroscopic observations, no definitive conclusion can be reached about its nature.Given its large opening angle (∼ 35 • above the midplane), this feature must arise in the highest layers of the disk, or in a wind above the disk, and over a broad range of stellocentric distances.In this context, it is worth recalling two other unusual features of this system: 1) the presence of scattered light in the optical up to ≳ 250 au above the midplane, and 2) the unusually extended CO channel map at a very small velocity relative to the systemic velocity, also marking surfaces that are about ∼ 325 au above the midplane.The similarity in elevation between these various features suggest a common origin.
Since the diffuse optical emission observed in this region is found in two separate filters, it is unlikely to be due to line emission and instead points to a scattered light nature.The absence of this component in even the shortest wavelength JWST image suggests scattering off of small (≲ 1 µm) grains.Conversely, the X-shaped feature is unlikely to be solely due to scattering since it is not detected in the 4.4 µm image.We speculate that it could be due to H 2 photodissociation, or to out-ofequilibrium thermal emission from very small (≲ 0.1µm) grains in the disk atmosphere that are in direct view of the X-and UV-emitting central source.In the H 2 interpretation, we note the lack of detection of the S(0) line (at 28.2 µm) in the IRS spectrum and the absence of the X-shaped feature in the 4.4 µm image, which contains the S( 9) line (at 4.7 µm).In conjunction, these would suggest an excitation temperature of ≲ 100 K (Nomura & Millar 2005;Lahuis et al. 2007).Polycyclic aromatic hydrocarbons (PAHs) are commonly observed in the 3-15 µm range in disks surrounding more massive Herbig AeBe stars (Meeus et al. 2001), and it is plausible that the active accretion from the central source in Tau 042021 can excite such grains even though it is not intrinsically as hot (Siebenmorgen & Heymann 2012).We note that, while the X-shaped feature is prominent in the JWST images presented here, it represents a small fraction of the integrated brightness of the system, so that PAHs could remain undetected in the Spitzer/IRS spectrum.It is also possible that the out-of-equilibrium grains have a more standard composition, in which case their stochastic emission would appears as a sum of warmer-than-expected blackbodies (Draine & Li 2001).Ultimately, the only requirement for this process to occur is that the grains are very small, to the point where the absorption of a single UV photon can overheat them for a short time.If this very small grain interpretation is correct, it could be common among T Tauri systems, and mid-infrared integral field unit spectroscopy of a large sample could confirm this conclusion.
A plausible scenario to explain these various features invokes a photodissociative disk wind (e.g., Owen et al. 2011;Bai et al. 2016).In this picture, the UV emission from the central star (and accretion column) dissociates molecules located at the very top of the warm molecular layer.The excited atoms escape the disk surface and are then placed on ballistic outward trajectories.Recombination of H 2 molecules along them could then produce the X-shaped feature through limb brightening of the conical wind, while the non-Keplerian flow of the wind could possible lead to the odd-looking low-velocity CO channel map.It is also possible that PAHs and other very small grains are entrained along with the gas and get stochastically heated once they reach high enough layers from which they are in direct sight of the central source.Furthermore, if this wind can entrain grains as large as 0.1-1 µm, this would naturally produce a "veil" of optical scattered light well above the disk midplane, as is observed in Tau 042021.Current disk wind models suggest that entrainment of such grains is possible (e.g., Booth & Clarke 2021;Rodenkirch & Dullemond 2022;Franz et al. 2022).In summary, such a scenario can qualitatively reproduce some of the unexpected features observed in this system, but we defer to a later study a physically-grounded wind model.

CONCLUSION
In this study, we have presented new optical through suubmillimeter high-resolution observations of the Tau 042021 edge-on protoplanetary disk.With a ≈ 1000 au radius, this disk is one of the largest disks known in nearby star-forming regions, and it shows a smooth and symmetrical structure.This makes it a perfect prototype to study the disk structure and dust properties in such environments.
We have presented the first JWST observations of Tau 042021, imaging its edge-on protoplanetary disk from 2 to 21 µm, which we interpret in combination with HST images at 0.6 and 0.8 µm.The disk is observed in scattered light throughout the entire wavelength range considered here, and shows very modest levels of chromaticity.In particular, the thickness of the dark lane that separate the two disk surfaces only declines by ∼ 15% from 0.6 to 21 µm, a much shallower decline than has been observed for other edge-on disks (over more limited wavelengths ranges).We also find that up to 21 µm, the vertical extent of the disk is an order of magnitude larger than observed in ALMA continuum maps, a confirmation that the large, millimeter-emitting grains, must be strongly settled.To interpret these observations, we build radiative transfer models that explore the role of the grain size distribution and degree of dust settling in the disk.We conclude that grains up to ≳ 10µm in size are fully coupled to the gas up to the disk surface layers, whereas particles larger than ≳ 100µm must be concentrated ∼ 10 times closer to the midplane.This is the first time observations provide a direct constraint on the degree of settling of grains in the 10 µm size regime via a spatially-resolved analysis of dust opacity.Future studies of other edge-on disks will help place this system in context and identify the underlying physical processes that are responsible for settling and, ultimately, growth, of dust particles in disks.
In addition to analyzing the scattered light properties of the disk, the HST and JWST images reveal multiple knots of emission from the jet launched by the central source.Multiple emission line features are also detected in the Spitzer/IRS spectrum of the system.The ALMA CO observations yield an estimate of the stellar mass (∼ 0.4 M ⊙ ) that is on the low end of the range suggested by past spectroscopic observations.Finally, we have also identified several unexpected in this system.The HST optical images point to the presence of ≲ 1 µm grains located up to 250 au above the disk midplane.The 12 CO channel maps reveal a ∼ 300 auelevation CO layer with unusual kinematic properties, suggesting non-Keplerian motion.Finally, and most intriguingly, the 7.7 and 12.8 µm JWST images reveal a X-shaped feature that is located above the CO warm molecular layer.We propose that these various features can be explained if the disk drives a significant photoevaporative disk wind that entrains small dust grains.Follow-up spectroscopic observations of this system can help identify the nature of the X-shaped feature and its relationship to the other phenomena.To produce the complete SED of the system, we complement the new JWST photometric datapoints presented in this study and the IRS spectrum (Section 3.1) with published and archival data.The complete SED is presented in Figure 3 while Table 4 summarizes the photometry adopted here (see below), as well as the new JWST fluxes.
Tau 042021 was detected by GALEX with its NUV channel (0.18-0.28 µm) in a 2662 s exposure as part of the Medium Imaging Survey.We retrieved the corresponding flux density from the GR6/GR7 source catalog available through the MAST portal2 .The GALEX detection suggests that the source is actively accreting, consistent with the presence of the collimated atomic jet, and that it is not significantly extincted by foreground cloud material.No corresponding FUV (0.13-0.18 µm) observation was taken.In the optical range, we adopted the SDSS DR16 (Ahumada et al. 2020) "Petrossian" magnitudes of Tau 042021.These represent an aperture photometry method with a radius based on the extent of the source; while this method was designed for galaxies, we assume that it operates reasonably well for any smooth, extended source.In the far-infrared, Tau 042021 was observed with Herschel as part of the "Herschel Gould Belt Survey" Guaranteed Time Key Program (André et al. 2010).PACS flux densities are available through the Herschel PACS Point Source Catalog (Herschel Point Source Catalogue Working Group et al. 2020).To complete the SED, we retrieved the level 2.5 SPIRE (250, 350 and 500 µm) data products from the Herschel archive and averaged the individual frames to produce final images in each filter.Tau 042021 is clearly detected -and unresolved -in all filters.Following the SPIRE Handbook, we estimated photometry for the system using Gaussian fitting and applying bandappropriate aperture corrections to obtain the far-infrared flux densities of Tau 042021.The resulting measurements, for which we estimate a 10% relative uncertainty due to a combination of absolute calibration error and uncertainties due to the complex background, are listed in Table 4.We note that a "starless core" was identified by Marsh et al. (2016) at the location of Tau 042021, using the same Herschel data.The SPIRE flux densities estimated by these authors are broadly consistent with our estimates.However, their 160 µm measurement is about 33% lower than that estimated by Herschel Point Source Catalogue Working Group et al. ( 2020), probably as a result of methodological differences.As a result, we decided against adopting their Herschel photometry both to ensure consistency and to avoid an implausible break between 100 and 160 µm.
In some regimes, Tau 042021 is detected in multiple catalogs and surveys, which reveals significant variability in some cases.In the optical, the Pan-STARRS Kron photometry (Chambers et al. 2016) is consitent with the SDSS photometry adopted here to within 20% or so.On the other hand, the UKIDSS photometry (Lawrence et al. 2012) is a factor of ≈ 2 fainter than the 2MASS photometry, and the 3.4 and 4.6 µm WISE flux densities (Cutri & et al. 2012) are about 30% lower than the corresponding Spitzer/IRAC ones (Rebull et al. 2010).Finally, the integrated flux densities measured in the JWST images presented here are 20-60% higher than the Spitzer measurements in the 4-8 µm, whereas the 20 µm JWST flux is consistent with the 24 µm Spitzer one.All in all, the source is affected by a typical variability of up to a factor of ≈ 2 from the optical to the mid-infrared, which is in line with the morphological variability noted by Luhman et al. (2009) in the near-infrared.Here we selected datasets on the basis of producing a smoothly continuous SED but we caution against interpreting any photometric point without taking into consideration this significant variability.
Finally, we performed power law fits to the (sub)millimeter and UV/blue regimes of the SED, which led to spectral indices of α mm = 2.5 (see also Villenave et al. 2020) and α UV = −2.0,defined by F ν ∝ ν α .In particular, at λ ≲ 0.5 µm, the system's SED presents a markedly different spectral index from the rest of the optical regime, suggesting the presence of a significant accretion-driven UV excess and we use the latter index in our modeling setup (Section 4.1).

B. DUST OPACTIY IN MOLECULAR CLOUDS
The extinction law in the interstellar medium (ISM) and molecular clouds has long been used as a probe into the grain size distribution in these environments (Draine 2003, and references therein).In particular, the realization that the extinction law in molecular cloud flattens out in the near-to mid-infrared range revealed the presence of larger grains in these clouds than in the diffuse ISM (e.g., Lutz 1999;Flaherty et al. 2007), suggesting that the initial stages of grain growth occur before the formation of stars and their disks.To determine the extent and timescale of grain growth in disks, it is therefore necessary to develop a model of the dust population of molecular clouds, whose properties can then be compared to observations of protoplanetary disks, such as those presented in the present study.
Unlike the case of the diffuse Galactic ISM, the extinction law in molecular clouds is not a unique curve, as it depends, among other things, on the temperature and column density of the cloud.Furthermore, no single observational method can probe extinction from the short-wavelength optical to the mid-infrared, requiring combining multiple studies that may inherently probe different environments.Fortunately, the agreement between studies is sufficient for our purposes, as we show below.
To construct the extinction law in the molecular clouds, we focus on Spitzer studies.Specifically, Flaherty et al. (2007) and Chapman et al. (2009) provided broadband extinction laws ranging from 3.6 to 24 µm for multiple nearby star-forming regions.In both cases, we computed the weighted average of the individual line of sight extinction laws.We note that the 1 ≤ A K < 2 and A K ≥ 2 curves from Chapman et al. (2009) bracket the 24 µm estimate from Flaherty et al. (2007) and are comparable to the latter at λ ≤ 8 µm.Lim et al. (2015) extend broadband extinction through an infrared dark cloud to longer wavelengths (12, 24 and 70 µm).Again, we computed the average of the various line of sights they studied and anchored their relationship to the 8 µm extinction coefficient from Flaherty et al. (2007). Finally, McClure (2009) used 5-20 µm spectroscopy to sample the extinction law across the silicate feature.To improve signal-to-noise, we resampled their extinction curve to a resolution of R ≈ 15.Altogether, we find that these different extinction laws agree to ≈ 10% across the near-to mid-infrared range (Figure 15).To produce a representative extinction law, we adopt the Flaherty et al. (2007) datapoints from 3.6 to 24 µm, together with the  et al. (2009); the latter corresponds to the lowest value of A24/AK .The dotted orange line indicates the opacity law of the "pristine dust" model that was constructed to match the adopted molecular cloud extinction law (represented by the solid symbols).
Figure 1.HST-JWST 3-color composite image of Tau 042021; the full field of view spans 18 ′′ .The 0.8, 2 and 7.7 µm images are rendered (using a logarithmic stretch) in the blue, green and red channels, respectively.The orange object to the right of the edge-on disk is a background galaxy.

Figure 1
Figure 1 presents a 3-color combination of selected HST and JWST images of Tau 042021, whereas Figure 2 compares all HST and JWST images along with the ALMA band 7 continuum map.All images from the optical to the mid-infrared reveal a similar morphology, namely two nearly symmetrical nebulae bisected by a dark lane, as expected for EODs.The system's SED has its trough around 20 µm (see Figure3), and the 21 µm image confirms that the disk is still optically thick to its own radiation out to that wavelength.In the ALMA regime, the disk is either optically thin (and emission is dominated by the denser midplane) or the

Figure 2 .Figure 3 .
Figure 2. Image gallery with the disk oriented horizontally.All images are shown on same angular scale and a log stretch for all panels except at 21 and 890 µm (for which square root and linear stretches are used, respectively).The corresponding beam is shown in the bottom left of each panel.A faint star is detected from 0.8 to 12.8 µm to the SW of the disk.The bright spike indicated by a ⊕ symbol in the 0.6 µm image is an uncorrected bad pixel "enhanced" due to Fourier shifting of the image. •

Figure 4 .
Figure 4. Left: Vertical cut across the disk from the optical to the mid-infrared, measured within 0. ′′ 25 of the symmetry axis.All curves are normalized to their peak and aligned based on the location of the brighter nebula, except for the 890 µm ALMA profile, which is aligned with the dark lane seen in the scattered light images.Right: Distance between the top and bottom nebulae, normalized to the 0.8 µm value for comparison purposes, as a function of wavelength for Tau 042021 (red stars) and other EODs (black symbols).Data for other systems are from McCabe et al. (2011), Duchêne et al. (2010), and Cotera et al. (2001) andWatson & Stapelfeldt (2004) for HK Tau B, HV Tau C and HH 30, respectively.

Figure 5 .
Figure 5. HST and JWST images of Tau 042021 after application of a 1 ′′ box median filter.The images are rendered with a hard logarithmic stretch to highlight the jet component of the system, leading to saturation in the scattered light disk and X-shaped feature areas.

Figure 6 .
Figure 6.Spitzer IRS spectrum of Tau 042021 with H2, H and [Ne II] lines indicated by the solid blue, dotted red and dashed orange segments, respectively.

Figure 7 .Figure 8 .
Figure 7. MIRI 7.7 µm images of Tau 042021and a bright star in the same frame (left and right panels, respectively).Both images are shown in the native detector orientation and with a hard log stretch to highlight the X-shaped feature and the low-level structure of the PSF.Notice the smaller field-of-view of the right panel.To guide the eye, red dashed segments are drawn at 36 • on either side of the disk midplane and replicated in the right panel.

Figure 9 .
Figure 9. Moment 0 and 1 maps of the 12 CO emission from Tau 042021 (top and bottom, respectively).A 3σ clipping was applied in each datacube prior to computing the moment maps.

Figure 10
Figure10. 12CO and 13 CO position-velocity diagrams (top and bottom, respectively).The dotted, solid and dashed curves correspond to dynamical masses of 0.25, 0.4 and 0.8 M⊙, respectively.The lack of low-velocity 12 CO emission is due to interferometric filtering of fore/background cloud emission.

Figure 12 .
Figure 12.Extinction law of dust models used in this study.The black and gray symbols represent the extinction law observed in molecular clouds (see Appendix B).

Figure 13 .
Figure 13.Model and observed images for various dust grain size distribution and disk structure models.Corresponding parameters are summarized in Table3.All images are shown on a square root stretch except the last column (linear stretch) from 0 to the peak value.

Figure 14 .
Figure14.Left: Dark lane thickness as a function of wavelength (left panel) for the models explored here (see Table3).Right: Observed and model SED for the same models.The dip at 9 µm in the model SEDs is due to a strong albedo feature from the silicates adopted in the DIANA mixture.We stress that the models are constructed solely on the basis of matching the dark lane thickness at 2 µm; all other quantities are merely predicted based on having fit that quantity.

Facilities:
JWST (NIRCam and MIRI), HST (ACS), ALMA APPENDIX A. COMPILING THE SED OF THE SYSTEM

Figure 15 .
Figure 15.Near-to far-infrared extinction law measured in molecular clouds (colored symbols) as well as function forms for the interstellar extinction from Weingartner & Draine (2001) (labeled WD01 in the figure).The Case A and B curves are appropriate for the diffuse ISM and molecular clouds, respectively.Solid circles represent Case B extinction values for standard optical and near-infrared filters.The two sets of blue triangles represent the 1 ≤ AK < 2 and AK ≥ 2 results from Chapmanet al. (2009); the latter corresponds to the lowest value of A24/AK .The dotted orange line indicates the opacity law of the "pristine dust" model that was constructed to match the adopted molecular cloud extinction law (represented by the solid symbols).

Table 2 .
Disk Morphological properties

Table 3 .
Summary of Model Exploration This project has received funding from the European Research Council (ERC) under the European Union's Horizon Europe research and innovation program (grant agreement No. 101053020, project Dust2Planets, PI F.Ménard).GD, KRS, MV, and SGW acknowledge funding support from JWST GO program #2562 provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.We are grateful to the W.M. Keck Institute for Space Studies as the hosts for our week-long team meeting in May 2023.CP acknowledges Australian Research Council funding via FT170100040, DP18010423 and DP220103767.RT acknowledges financial support from a CNES postdoc fellowship.MV research was supported by an appointment to the NASA Postdoctoral Program at the Jet Propulsion Laboratory, administered by Oak Ridge Associated Universities under contract with NASA.Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute.The specific JWST observations analyzed can be accessed via doi:10.17909/rgg7-zk46.This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.0771.S.

Table 4 .
Adopted system SED for Tau 042021