<figure><img src="image_1.png" /></figure>

# JWST Imaging of Edge-on Protoplanetary Disks. II. Appearance of Edge-on Disks with a Tilted Inner Region: Case Study of IRAS04302+2247

Marion Villenave¹, Karl R. Stapelfeldt¹, Gaspard Duchêne²,³, François Ménard³, Schuyler G. Wolff⁴,¹

Marshall D. Perrin⁵, Christophe Pinte³,⁶, Ryo Tazaki³, and Deborah L. Padgett¹

¹Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA; marion.f.villenave@jpl.nasa.gov

²Astronomy Department, University of California, Berkeley, CA 94720, USA

³Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

⁴Department of Astronomy and Steward Observatory, University of Arizona, Tucson, AZ 85721, USA

⁵Space Telescope Science Institute, Baltimore, MD 21218, USA

⁶School of Physics and Astronomy, Monash University, Clayton Vic 3800, Australia

Received 2023 September 15; revised 2023 November 10; accepted 2023 November 13; published 2024 January 17

## Abstract

We present James Webb Space Telescope imaging from 2 to 21 µm of the edge-on protoplanetary disk around the embedded young star IRAS04302+2247. The structure of the source shows two reflection nebulae separated by a dark lane. The source extent is dominated by the extended filamentary envelope at ~4.4 µm and shorter wavelengths, transitioning at 7.7 µm and longer wavelengths to more compact lobes of scattered light from the disk itself. The dark lane thickness does not vary significantly with wavelength, which we interpret as an indication for intermediate-sized (~10 µm) grains in the upper layers of the disk. Intriguingly, we find that the brightest nebula of IRAS04302 switches side between 12.8 and 21 µm. We explore the effect of a tilted inner region on the general appearance of edge-on disks. We find that radiative transfer models of a disk including a tilted inner region can reproduce an inversion in the brightest nebula. In addition, for specific orientations, the model predicts strong lateral asymmetries, which can occur for more than half possible viewing azimuths. A large number of edge-on protoplanetary disks observed in scattered light show such lateral asymmetries (15/20), which suggests that a large fraction of protoplanetary disks might host a tilted inner region. Stellar spots may also induce lateral asymmetries, which are expected to vary over a significantly shorter timescale. Variability studies of edge-on disks would allow us to test the dominant scenario for the origin of these asymmetries.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Planet formation (1241); Radiative transfer (1335); Dust continuum emission (412)*

## 1. Introduction

Exoplanet statistics have shown that most stars have a planet orbiting them (Winn & Fabrycky 2015; Vigan et al. 2021). Because planet formation occurs in the protoplanetary disk phase, studying protoplanetary disk evolution can allow us to better understand planet formation. In particular, it is not yet fully understood how planets can grow from micron-sized dust to kilometer-sized bodies in a few million years, or possibly less. In the current paradigm, high dust concentrations are thought to accelerate grain growth by promoting disk instabilities that lead to planetesimal formation (e.g., streaming instability; Youdin & Lithwick 2007), and subsequently allowing efficient growth via pebble accretion (Lambrechts & Johansen 2012).

Dust vertical settling in the disk is the result of gas drag on dust grains subject to stellar gravity and gas turbulence. This mechanism leads large dust grains to fall into the disk midplane and accumulate there, which is favorable for planet formation. However, this mechanism remains poorly observationally constrained. For relatively large dust particles visible at millimeter wavelengths, it is currently accepted that vertical settling is efficient in the Class II phase, that is typically for disks without a remaining envelope. Although still limited to a

small number of systems, the disk emission of Class II systems at millimeter wavelengths appears geometrically thin, with scale height typically less than 4 au at a radius of 100 au (Pinte et al. 2016; Doi & Kataoka 2021; Liu et al. 2022; Villenave et al. 2022; Pizzati et al. 2023). Yet, no clear observational constraints have been obtained regarding smaller-sized dust particles. In particular, the size of the smallest grains that do experience vertical settling is currently unknown.

Furthermore, limited information is currently available regarding the timescale of vertical settling and its efficiency within younger systems, which still possess part of their primordial envelope. The recent results from the eDISK large program (Ohashi et al. 2023) and two recent independent studies on the embedded disk IRAS04302+2247 (Lin et al. 2023; Villenave et al. 2023), indicated that large dust is not very settled in Class 0 and I systems. This suggests that most vertical settling of large (~millimeter)-sized particles might occur between the Class I and Class II phases.

Highly inclined protoplanetary disks are favorable targets to investigate this mechanism because they allow a direct view of the disk's vertical structure. This paper is part of an ongoing James Webb Space Telescope (JWST) near- to mid-infrared imaging campaign of the largest edge-on disks in nearby star-forming regions (GO programs 2562 and 4290 in Cycles 1 and 2, co-PIs: F. Ménard and K. R. Stapelfeldt). By looking at edge-on disks at near- to mid-infrared wavelengths, the goal of the program is to study vertical dust settling and grain evolution. While the programs mostly include Class II disks, we focus here

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
on the study of the youngest object of the cycle 1 survey, the
embedded IRAS04302+2247 (hereafter IRAS04302).

The system is located in the L1536 cloud in the Taurus star-
forming region (d = 161 ± 3 pc; Galli et al. 2019) and orbits
around a young 1.6 M☉ star (Lin et al. 2023). IRAS04302 has
previously been the focus of a number of studies, using
observations from the optical to the centimeter. The source has
been classified as being in the Class I phase, based on the shape
of its spectral energy distribution (e.g., Kenyon & Hartmann
1995) or its bolometric temperature (Ohashi et al. 2023), and
consistently with its morphology in scattered light (e.g., Lucas
& Roche 1997; Padgett et al. 1999; Wolf et al. 2008). High
angular resolution millimeter continuum observations resolved
a brightness asymmetry along the minor axis, with the east side
being brighter than the west (Lin et al. 2023). Similarly, a
blueshifted outflow and millimeter molecular line enhancement
have also been observed toward the east of the source (Podio
et al. 2020; van’t Hoff et al. 2020). This indicates that the disk
is slightly offset from edge-on, with the eastern side tilted
toward us. Finally, previous studies also found that its
millimeter emission is not significantly concentrated in the
midplane (Lin et al. 2023; Villenave et al. 2023).

In this paper, we present the first mid-infrared observations
of the system, obtained with JWST. The data reduction and
image results are discussed in Sections 2 and 3. In Section 4,
we construct a toy model to explore the impact of a tilted inner
region on the appearance of edge-on disks and compare the
modeling predictions with our JWST observations of
IRAS04302. Then, in Section 5, we discuss the implications
of our modeling results for previous scattered light observa-
tions of a sample of 20 edge-on disks. In Section 5, we also
compare disk appearance for three systems in different
evolutionary stages observed with the JWST, aiming to obtain
constraints on vertical settling and grain properties. Finally, we
summarize the results in Section 6.

## 2. Data Reduction

We observed IRAS04302 with JWST (GO program 2562)
using both NIRCam and MIRI instruments in two consecutive
visits starting on 2023 February 20 UT 01:20. The observations
used five different filters, F200W (2.0 µm), F444W (4.44 µm),
F770W (7.7 µm), F1280W (12.8 µm), and F2100W (21.0 µm).
The NIRCam images were obtained simultaneously with
exposures of 859 s. The MIRI observations were taken
successively with exposures of 633 s at 7.7 µm, 1672 s at
12.8 µm, and 329 s at 21 µm.

For the F200W, F770W, F1280W, and F2100W observa-
tions, we obtained the phase 2 pipeline calibrated data from the
MAST archive, and we re-ran phase 3 using the JWST pipeline
version 1.9.5. In the skymatch step, we used the default
skymethod "global+match" but set subtract=True.
This method first equalizes sky values between the four
exposures taken in each filter. It then finds a minimum sky
value among all input images, and subtracts this background
value from each frame. The average background levels were
0.36, 11.60, 46.70, and 272.27 MJy sr⁻¹ in F200W, F770W,
F1280W, and F2100W, respectively.

In addition, the pipeline-reduced MIRI 12.8 and 21 µm
images show a gradient in the sky brightness between the top
and bottom parts of each frame. Because it appears in each
individual, dithered image, it is likely some instrumental
 artifact. Thus for those two filters, we performed an additional

2D background subtraction step using the Background2D
function in the photutils python package.

For the F444W observations, we obtained the uncalibrated
data from the archive and re-ran each pipeline step using
pipeline version 1.11.1.dev16+gb79a88af. In order to recover
previously saturated pixels, we set suppress_one_group=False in the ramp fit step of phase 1. This allows
us to suppress the ramp fit in the case when only the 0th group
is unsaturated. We use the default pipeline values in phase 2
and 3, except for the skymatch step where we subtract the
background. The average background level in F444W was
0.49 MJy sr⁻¹.

Finally, we resample and align each image to the NIRCam
F200W pixel size and position. We adopt the reproject
python package to resample the images, and align them using a
background galaxy visible in all fields. The alignment does not
require a rotation, only a translation in R.A. and decl. The final
images are presented in Figure 1.

In this work, we also use archival HST NICMOS images,
from program 7418 (PI: D. Padgett). The HST observations
were obtained using filter F205W (2.05 µm). They were
previously published by Padgett et al. (1999). These observa-
tions were taken on 1997 August 19, nearly 25.5 years prior to
the JWST data. For this study, we use the data directly
downloaded from the archive.

## 3. Results

### 3.1. Disk Morphology

We show the five new JWST images of IRAS04302 in
Figure 1. All images present a similar overall shape, displaying
two extended nebulae, separated by a dark lane. This
morphology is typical of scattered light observations of edge-
on sources, indicating that scattering dominates up to 21 µm.
However, we note that the origin of the photons likely varies
between the different wavelengths. Indeed, while at short
wavelengths (e.g., 2 µm) the star dominates, thermal emission
from the hot inner disk progressively becomes important over
the stellar emission as the wavelength increases toward the
mid-infrared. The morphology of the images shows that the
disk is optically thick, even to its own radiation, up to 21 µm.

Besides a globally similar overall shape, the images reveal a
clear change of morphology around 7.7 µm, where the apparent
size of the source gets dramatically smaller. This marks the
limit up to which scattering off the envelope dominates over
scattering off the disk. At short wavelengths (2.0, 4.44 µm), the
scattered light reveals two nebulae surrounded by diffuse
extended emission which includes some filamentary features.
At 2.0 µm, both the east nebula and the south west filament are
visible at distances greater than 20" from the central dark lane.

At 7.7 µm, the source appears much more compact, with
most diffuse emission from shorter-wavelength observations
not detectable. Two nearly vertical lines, along the major axis
direction, are also visible at this wavelength. We note however
that these lines are not physical, but instead are related to the
known "cross artifact" at this wavelength (Gáspár et al. 2021),
also referred to as the cruciform artifact in the MIRI user
documentation.

Finally, at longer wavelengths (12.8, 21 µm) the diffuse
emission disappears totally below the rising thermal back-
ground. At these wavelengths, IRAS04302 shows two smooth
nebulae separated by a dark lane. The absence of diffuse
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Image gallery of the JWST observations of IRAS04302. All images are shown with a log stretch. The bottom left ellipse in each panel indicates the corresponding beam size. All images are shown to the same angular scale, and the square box in the 4.44 µm image represents the extent of the 7.7, 12.8, and 21 µm smaller panels.</figcaption></figure>

<table><caption>**Table 1**
IRAS04302 Morphological Properties</caption><thead><tr><th>λ<br>(μm)</th><th>F<sub>20''×20''</sub><br>(mJy)</th><th>d<sub>neb</sub><br>(arcsec)</th><th>FR<br>W/E</th><th>R<sub>FWHM</sub><br>E (arcsec)</th><th>R<sub>FWHM</sub><br>W (arcsec)</th><th>R<sub>FW10%</sub><br>E (arcsec)</th><th>R<sub>FW10%</sub><br>W (arcsec)</th><th>δ<sub>spines</sub><br>W-E (%)</th></tr></thead><tbody><tr><td>2.00<sup>(E)</sup></td><td>17 ± 2</td><td>1.35 ± 0.10</td><td>0.99 ± 0.01</td><td>2.17 ± 0.25</td><td>1.75 ± 0.43</td><td>3.69 ± 0.13</td><td>3.83 ± 0.17</td><td>12</td></tr><tr><td>4.44<sup>(E)</sup></td><td>32 ± 3</td><td>1.13 ± 0.10</td><td>0.58 ± 0.01</td><td>1.71 ± 0.10</td><td>1.63 ± 0.15</td><td>2.97 ± 0.10</td><td>3.78 ± 0.19</td><td>11</td></tr><tr><td>7.70</td><td>16 ± 2</td><td>1.03 ± 0.10</td><td>0.52 ± 0.01</td><td>1.35 ± 0.14</td><td>1.63 ± 0.10</td><td>2.76 ± 0.10</td><td>3.57 ± 0.10</td><td>7</td></tr><tr><td>12.8</td><td>11 ± 1</td><td>1.10 ± 0.10</td><td>0.78 ± 0.04</td><td>1.30 ± 0.16</td><td>1.43 ± 0.10</td><td>2.69 ± 0.10</td><td>3.12 ± 0.15</td><td>8</td></tr><tr><td>21.0</td><td>101 ± 10</td><td>0.96 ± 0.10</td><td>1.34 ± 0.05</td><td>1.31 ± 0.10</td><td>1.17 ± 0.10</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Note.** $F_{20''×20''}$ are the flux densities integrated over a 20'' × 20'' aperture, $d_{neb}$ the dark lane width, FR the peak flux density ratio of the west over the east nebula, while $R_{FWHM}$ and $R_{FW10%}$ correspond to the radial size of each nebula, respectively measured at 50% or 10% of the peak. $\delta_{spines}$ characterizes the difference in lateral centering between both nebulae. At 21 µm, the signal to noise is not sufficient to measure $R_{FW10%}$. The error bars of $d_{neb}$, $R_{FWHM}$, and $R_{FW10%}$ correspond to the statistical error for the different spine averaging or $0''/1$, whichever is greater. A positive value of $\delta_{spines}$ indicates that the west nebula is more centered toward the south than the east nebula. <sup>(E)</sup> indicates that the images are dominated by the envelope; at these wavelengths our measurements might not trace the disk properties.

emission and smooth shape of the two nebulae suggest that, at these wavelengths, scattering from the disk dominates. In addition, we note that the disk is more centrally peaked and radially smaller at 21 µm than at 12.8 µm, which could be related to dust properties.

For all five filters, we use a 20'' square aperture to measure the overall source flux density. We report the results in the first column of Table 1. Then, we follow the methodology recently highlighted by Duchêne et al. (2023) to measure the structural source's parameters (dark lane width, brightness ratio, size, lateral asymmetries; see also Appendix D of Villenave et al. 2020). We parameterize the vertical position of the two scattering surfaces as a function of radius at each wavelength. This is done up to a radius of 1'' at 21 µm, or 3''/4 at the other wavelengths, and up to 0''/8 vertically from the midplane for all images. To do so, we take minor axis profiles averaged over 0''/4 or 0''/5 along the major axis of the source. For each nebula, we determine the position of the peak by fitting a polynomial of degree 6 or 8. This allows us to obtain four sets of arrays with the peaks coordinates of the east and west nebulae, to which we fit a second-order polynomial to generate the final spines.

We use these polynomial spines to estimate the dark lane thickness (closest distance between the nebulae, $d_{neb}$), the peak

flux ratio of the west over the east nebula (FR), radial full width half maximum ($R_{FWHM}$) and the radial full width at 10% of the peak ($R_{FW10%}$) of each spine. In addition, we characterize the lateral offset between the spines by the parameter $\delta_{spines}$. This parameter quantifies the difference between the median location of each spine, within 10% of their respective peaks. To allow comparison with disks of different sizes, the lateral difference is normalized by the averaged disk size (i.e., mean $R_{FW10%}$, respectively at each wavelength). This final number is reported as $\delta_{spines}$.

We present the average of these values obtained with the four sets of spines and corresponding standard deviation in Table 1. We note that when building the spine without the radial averaging and polynomial fitting $d_{neb}$ is consistently found to be ~0''/1 (~3 pixels) smaller than the values reported in Table 1, except at 2.0 µm where it is 0''/3 broader (highlighting the uncertainty of the methodology when the image appears more messy). Moreover, we caution that for the shortest wavelengths, $\lambda < 7$ µm, the images are likely dominated by the envelope such that our measurements might not trace the disk properties.

Focusing on the wavelengths >7 µm which are mainly sensitive to the disk, Table 1 shows that the west side of the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Comparison of the 2.05 µm HST (left panel) and 2.0 µm unconvolved JWST (middle panel), previously aligned images of IRAS04302. The contours correspond to an intensity level of 5% of the peak (blue is HST, black is JWST), and the white line in the left panel shows the orientation of the normalized cut, displayed in the right panel. A 1D Gaussian was fitted to both cuts (marked in red), and the vertical lines in the right panel indicate the resulting centering of these Gaussians for each data set.</figcaption></figure>

disk is more radially extended than the east side, both at 7.7 and 12.8 µm. In addition, the value of $\delta_{\text{spines}}$ is of order 10%, which illustrates that the two nebulae are not perfectly aligned with each other. The west nebula is found to extend further south than the east nebula at 7.7 and 12.8 µm. Finally, we find that the dark lane thickness does not vary significantly with wavelength (by less than 10% for $\lambda > 4$ µm), which we discuss further in Section 5.2.

Surprisingly, we also notice that the brightest nebula flips with wavelength. At 7.7 and 12.8 µm, the east nebula is brighter than the west nebula, while the opposite is true at 21 µm, where the west nebula is the brightest. We note that Table 1 reports the peak brightness ratio, but this is also true for the integrated brightness of both nebulae. In Section 4, we perform radiative transfer models aiming to reproduce this effect, and we discuss the results in Section 5.1.

## 3.2. Comparison of Two Epochs

IRAS04302 was observed multiple times at ~2 µm (Lucas & Roche 1997; Padgett et al. 1999). In this section, we compare the 2023 JWST observations and the previous 1997 HST 2.05 µm image of Padgett et al. (1999). We aim to look for moving features between the two observations.

To do so, we first aligned the two images. Because no stars were visible in the background of the HST image, we cross-correlated the HST and JWST images to align them. We first resampled the JWST NIRCam image to the coarser HST NICMOS pixel scale. We normalized the JWST and HST images to their respective total counts and, to align them, we then shifted the JWST data to minimize the normalized residuals. We performed this alignment both before and after convolving the JWST observations with an HST point-spread function (PSF), generated using the TinyTim software (Krist et al. 2011).

We display the aligned images and significant contours in Figure 2. The JWST observations reach a signal to noise significantly better than the previous data, such that diffuse emission is seen to a much larger distance. Filaments features that were barely visible with HST are now confirmed with the new image. In addition, we observe some differences in the source's morphology. In particular, around -1.5 eastward of the star, the source is dimmer in the JWST observations than in

the older HST images. This is possibly due to time variation of the density distribution in the nebula and/or variable illumination from the star (e.g., Watson et al. 2007).

In addition, we find that a roundish blob toward the east side of the source (see the insert in Figure 2) seems to move toward the outer regions of the source. To quantify the movement of this blob feature, we estimate the location of the star in the disk-dominated 21 µm observations, and obtain the fluxes along a cut intercepting this star location and the globular feature. We show the fluxes obtained along this cut, normalized to their respective maximum, in the right panel of Figure 2. The cuts show two bright peaks with maximum within ±1", and a small increase of emission, around -3", which corresponds to the globular feature. For each cut, we fitted this small increase by a 1D Gaussian. On the JWST data, this was done twice, both for the HST-convolved image and that with the native resolution. We find that the center of this blob feature moves by 0.10 (resp. 0.16)<fnref n="7" /> from the star between the HST image and JWST convolved data (resp. native JWST data). This motion is pointing directly away from the presumed location of the star.

The relative proper motion of this blob compared to the disk location is +3.9 mas yr⁻¹ in R.A. (resp. +6.2 mas yr⁻¹), and -0.7 mas yr⁻¹ in decl. (resp. -1.1 mas yr⁻¹). Using Gaia and VLBI astrometry, Galli et al. (2019) estimated that the proper motion of the L1536 cloud in which IRAS04302 is located is +10 mas yr⁻¹ in R.A. and -17 mas yr⁻¹ in decl. The proper motion of the blob feature appears unrelated to the proper motion of the disk, which suggests that it is not a background object.

Assuming a distance to the disk of 161 pc (Galli et al. 2019), we estimate that the blob's apparent movement is 16 au (resp. 26 au) in ~25.5 yr. Thus, we infer an apparent velocity of 3.0 km s⁻¹ (resp. 4.8 km s⁻¹) for this feature. This velocity is significantly less than typical velocities of winds or jets (e.g., Pascucci et al. 2023). Yet, this could correspond to a knot of dust entrained in the low-velocity outflow. Indeed, if the dust lies at the outflow cavity edge, its actual velocity could be slightly larger, due to projection effects.

<footnote n="7">⁷ The difference in distance between the convolved and unconvolved data corresponds to ~1 HST pixel. It mostly due to the non-perfectly Gaussian shape of the excess.</footnote>
# 4. Radiative Transfer Modeling

## 4.1. Motivations

The new JWST observations of IRAS04302 show that the flux ratio between the two nebulae is not constant with wavelength (Section 3.1) and, in particular, we found that the brightest nebula switches side between 12.8 and 21 µm. This was an unexpected result. Indeed, at wavelengths between the optical and infrared, we expect forward scattering to be important. This implies that a disk that is slightly offset from edge-on will have one side that is brighter than the other, simply due to its orientation, and independent of the wavelength. Because the MIRI observations of IRAS04302 were performed consecutively within less than one hour, variability likely does not explain this behavior.

Multi-wavelength optical and near-infrared observations have previously been performed in about a dozen edge-on disks, mainly at wavelengths between 0.4 and 2 µm. As expected in the case of axi-symmetric radiative transfer models, the vast majority of systems (e.g., Haro6-5B, LkHa 263C, HV Tau C, PDS144N, HH30, HOPS136, ESO-Hα569, Oph163131, Tau042021; Krist et al. 1998; Chauvin et al. 2002; Stapelfeldt et al. 2003; Perrin et al. 2006; Watson & Stapelfeldt 2007; Fischer et al. 2014; Wolff et al. 2017, 2021; Duchêne et al. 2023) do not show a switch in the brightest nebula.

While the nebulae flux ratio in the edge-on disk HK Tau B varies significantly between 0.6 and 3.8 µm (by a factor ~8; McCabe et al. 2011), to our knowledge, the only other edge-on disk where the brightest nebula switches side between two wavelengths is 2MASS J16281370-2431391 (Flying Saucer), although at different wavelengths than in IRAS04302. Using VLT/ISAAC, Grosso et al. (2003) indeed found that the flux ratio of the two nebulae reverses between 1.3 and 2.2 µm. They suggested that a disk warp might be able to explain the differences by introducing some extra self-absorption in the disk.

Nealon et al. (2019) performed radiative transfer modeling of numerical simulations including a warp, allowing some initial testing of this hypothesis. They showed that, when a warp is included, a true edge-on disk ($i=90°$) observed at 1.6 µm can appear asymmetric, with the brightest nebula depending on the angle of observation of the system (their Figure 6). Here, we aim to extend this study to multiple wavelengths to analyze the impact of a tilted region in the nebulae flux ratio and other morphological properties in edge-on disks.

## 4.2. Model Description

We construct a toy model using the radiative transfer code `mcfost` (Pinte et al. 2006, 2009). The goal of this modeling is to illustrate the effect of a tilted inner disk onto edge-on disk appearance. While we construct our model with some disk and stellar parameters consistent with IRAS04302, we do not aim to reproduce this disk. Instead, we are interested in more general disk properties.

We take a stellar effective temperature and radius of $T_{\text{eff}} = 4500$ K, $R_* = 3.7R_\odot$ (~5$L_\odot$), respectively (Gräfe et al. 2013; Villenave et al. 2023). We assume that the disk is broken into two regions. The inner region extends from 0.1 to 5 au, while the outer region extends from 9 to 300 au. We set the dust masses to $1.6 \times 10^{-6} M_\odot$ and $5 \times 10^{-4} M_\odot$, respectively for the inner and outer disks.

For both regions, the surface density follows a simple power-law distribution, with $\Sigma(r) \propto r^{-1}$, and the vertical extent is parameterized as $H(r) = 41(r/300\ \text{au})^{1.14}$. We assume that dust grains are composed of a mixture of 62.5% of astronomical silicates and 37.5% of graphite (following Gräfe et al. 2013; Villenave et al. 2023), and only include small dust grains which follow a power-law distribution $n(a)da \propto a^{-3.5}da$. We fixed dust sizes between 5 and 30 µm rather than using a broader distribution. This choice aims to mimic the amount of forward scattering observed in IRAS04302, as the shape of the brighter side is defined by the forward scattering off the closer surface. We checked that changing the grains properties does not qualitatively alter the discussion below. Further exploration of the grain properties, including porosity, and aggregate shape will be performed in future work (R. Tazaki et al. 2023, in preparation).

We produce two models, one where the inner region is aligned to the outer disk (“no-tilt” or “axisymmetric” model) and another where the inner region is tilted by 10° compared to the outer disk (“tilted” model). We fix the outer disk’s inclination to 88° (Lin et al. 2023; Villenave et al. 2023), and consider 16 different angles of observations (azimuths) of the system with the tilted inner disk, separated by 22°.5. For different azimuths, the orientation of the outer disk remains the same while the orientation of the tilted inner disk varies. Azimuths 0° and 180° correspond to the orientation where the line of nodes defined by the inner disk and outer disk planes is pointed directly at the observer. A schematic view of the system highlighting four different viewing orientations is shown in Figure 3.

For each model, we produce images at 2.0, 4.44, and 7.7, 12.8, and 21 µm and convolve them by the JWST PSF obtained with the WebbPSF software (Perrin et al. 2014). As for the data (Section 3.1), the models are dominated by scattering at all wavelengths. We then employ the methodology of Section 3.1 to fit the spines of the models, and measure the variation of the flux ratio and lateral asymmetry $\delta_{\text{spines}}$ with the angle of observation (azimuth). We note that for the models, the radial averaging and polynomial fitting did not improve the quality of the spine, so we did not use them.

## 4.3. Results

In Figure 4, we show how the flux ratio and the lateral asymmetry vary with azimuth, for wavelengths between 2 and 21 µm. In this section, we first describe the model features for both the dependence of the flux ratio and the lateral asymmetry with azimuth. Then, in a second phase, we describe the reasons why a tilted inner disk may lead to such variations. Throughout this section the wording “flux reversal” refers to an azimuth where the brightest nebula in the tilted model is opposite to that of the no tilt model. We use “brightest nebula switch” when the brightest nebula switches side between wavelengths.

### 4.3.1. Flux Ratio

In the top panels of Figure 4, we find that while the models with no tilt do not show any clear flux ratio change with azimuth (less than 15% variations, due to radiative transfer noise), the tilted models do show significant variation with the observing azimuth, at all wavelengths. In the case where there is no tilt, the brightest nebula of this model is on the east, implying a flux ratio <1. While this is also true for most
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Schematic geometry of the model, not to scale. The orientation corresponding to an azimuth of 270° allows one to reproduce the flux reversal between 12.8 and 21 µm observed in IRAS04302 data (see Figure 7). In the bottom view, the observer is slightly above the outer disk midplane, the dashed inner region is below that midplane, and the squared inner region is above it.</figcaption></figure>

orientations in the tilted models, some azimuthal range show a flux reversal (flux ratio >1) at 21 µm and for $\lambda < 12.8$ µm. No flux reversal is found for the models at 12.8 µm, although there is some enhancement of the flux ratio within the range where the flux reversal is observed at 21 µm.

At 21 µm, the flux reversal occurs for azimuths between ~210° and 330° (see the schematic orientation in Figure 3), while it is observed for the opposite orientation between 2 and 7.7 µm (azimuths between ~40° and 140°). In both cases, the flux reversal only occurs for about around a third of the possible viewing orientations.

### 4.3.2. Lateral Asymmetries

In addition to generating flux ratio variations with azimuth, introducing a tilted inner region also leads to lateral asymmetries in the disk. The bottom panels of Figure 4 show the azimuthal variation of the parameter $\delta_{spines}$. This parameter, introduced in Section 3.1, quantifies the difference in centering between the spines of both nebulae, where the spines are only considered up to 10% of their respective peaks. By visually inspecting the models, we identify that a value of $\delta_{spines} \ge 5\%$ leads to significant lateral asymmetry, which can be directly detectable by eye.

We find that this threshold is reached for a similar range of orientations at 21 µm and between 2 and 7.7 µm. No clear azimuthal variation is found at 12.8 µm. At 21 µm, this threshold is reached for azimuths between ~0°–40°, 160°–210°, and 310°–360°. On the other hand, at shorter wavelengths the lateral asymmetries are stronger and overtake the threshold for a wider range of orientations (~0°–45°, 130°–230°, 300°–360°). The azimuths where lateral asymmetries are observed in the disk correspond to orientations where the observer is close to perpendicular to the tilted inner region (Figure 3).

Interestingly, we find that, at short wavelengths, important lateral asymmetries occur for a larger range of azimuths than the range where flux reversals were obtained: about 200° versus 100°. More than half of the possible orientations present strong lateral asymmetries at 2 µm, while about a third may lead to flux reversal. At short wavelengths, lateral asymmetries can thus be an easier and more frequent tracer of inner disk misalignment than strong flux variations. We also checked that at smaller inclinations (80° and 85°) this statement remains valid, with lateral asymmetries remaining clearer.

In Figure 5, we illustrate the lateral asymmetries by showing models at 2 and 21 µm observed at azimuth 0° and 180°, found to be associated to large values of $\delta_{spines}$. The spines, within 10% of their respective peaks, are represented on top of the models. From this figure, it is clear that the spines are not perfectly centered on each other and that the lateral asymmetry is stronger at shorter wavelengths. Moreover, we also see that the asymmetry is opposite at 2 µm and at 21 µm. For azimuth 0°, at 2 µm the east nebula is found to the north of the west nebula, while this is the opposite at 21 µm. This illustrates that the cause of the asymmetry, which is the same as for the flux reversal, is not the same at 2 and 21 µm.

### 4.3.3. Origins of the Azimuthal Variations

First, we recall that, because forward scattering is important at infrared wavelengths, a disk that is slightly offset from edge-on will have one side that is brighter than the other, simply due
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Variation of flux ratio (top panels) and lateral asymmetries (bottom panels) with azimuth, for the tilted (solid lines) and axisymmetric (dashed lines) disk models. On the top panels, flux reversal occurs above the horizontal dotted line, and the green region highlights the azimuth where flux is reversed (at 21 µm, left, and at 2 µm, right). The small azimuthal variations (<15%) in the no-tilt models are due to radiative transfer noise. On the bottom panels, the horizontal dotted lines show a value of $\delta_{spines} = \pm 5\%$, corresponding to strong lateral asymmetries. The blue regions show the azimuths where this threshold is surpassed (at 21 µm, left, and at 2 µm, right). The 12.8 µm model is shown on both columns, to scale, for comparison.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** 2.0 and 21 µm images of the tilted models, for azimuths of 0° and 180°. The locations of the spines within 10% of their respective peaks are shown in blue, highlighting clear lateral asymmetries.</figcaption></figure>

to its orientation. By definition, this is the east side in our models.

At short wavelengths (2.0, 4.44, 7.7 µm), stellar emission and very compact hot disk thermal emission are the dominant source of photons in the disk. The tilted inner region, however, affects the propagation of stellar light to the outer disk. In particular, shadows are expected where the inner disk blocks direct starlight (shadowed regions in the top panel of Figure 3). At the same time, the opposite side, which has a direct view to the star, will have more incoming photons and thus it

brightness will be enhanced, relatively. For azimuth 270° as pictured in Figures 3 and 6, this implies that the east nebula is even brighter than for a disk without a tilt. However, for the opposite orientation where the observer stands on top of Figure 3 (azimuth ~90°), the east side will be in the shadow and thus less bright than the west side. This causes the flux reversal observed in the models. For intermediate orientations (azimuth ~0° or ~180°; Figure 3), a similar effect is seen from the side. The sides in the shadow of the inner disk will appear less radially extended than their counterpart with a direct view to the star.

At longer wavelengths, however, the situation is different because thermal emission from the inner disk dominates over the stellar emission and originates from a radially extended region. In fact, the ratio of scattered photons originating from inner disk thermal emission to those originating from the star is more that 2.5 times higher at 21 µm than at 12.8 µm (ratios of 4300 versus 1600). The inner disk thermal emission, which is spatially very different from the central star, is added to stellar light, and scattered through the outer disk. Thus, the outer disk's side toward which the small inner disk is pointing will receive—and thus re-scatter—significantly more photons than the other side and appear brighter. For an azimuth of 270° as pictured in Figures 3 and 6, this implies that the west nebula is brighter than the east nebula and causes the observed flux reversal. At the opposite orientation (azimuth ~90°), the tilt enhances the east side. At intermediate orientations, we expect lateral asymmetries, but on the opposite side as for short wavelengths. In our models, we find that this effect is more important at 21 µm than at 12.8 µm because thermal emission is more important.

In this section, we have demonstrated that a tilted inner region can have a strong effect on the appearance of disks viewed close to edge-on at near- to mid-infrared wavelengths. On the one hand, the tilted inner region can lead to a change of the brightest nebula for about a third of the possible disk
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Schematic view of the photon propagation in the front half a disk for two regimes of wavelengths. The observing azimuth corresponds to 270°, and can reproduce a flux reversal between 12.8 and 21 µm. Forward scattering is more efficient than scattering in other directions as illustrated by the larger arrows.</figcaption></figure>

orientations. Such flux reversal occurs for opposite orientation in the mid-infrared, where thermal emission from the inner disk is important, than in the near-infrared, where stellar photons dominate. Moreover, a tilted inner region can lead to strong lateral asymmetries. This effect is stronger at near-infrared wavelengths and occurs for more than half of the possible viewing orientations. Finally, lateral asymmetries occur over a similar orientation range at near- and mid-infrared, but have opposite directions. While, quantitatively, the presence/degree of reversal and amplitude of lateral asymmetry are dependent on model parameters (e.g., dust properties, tilted inner disk extent and misinclination, viewing geometry), we note that the qualitative behavior is quite general for small misalignments, likely as long as they cast shadows over a broad range of azimuthal angles (e.g., HD139614; Muro-Arena et al. 2020).

## 4.4. Application to IRAS04302

The observations of IRAS04302 show a clear switch in brightest nebula between 12.8 and 21 µm (Section 3.1). In Section 4.3, we demonstrated that the presence of a tilted inner region can produce such a brightest nebula switch. Now, in Figure 7, we compare the data with the tilted and no tilt models observed at an azimuth of 270° (see Figure 3). We find that the tilted model reproduces the switch with a similar amplitude than the data. This would imply that the east side is tilted slightly toward us, and closest to the observer.

This proposed orientation is consistent with the conclusions from Lin et al. (2023). Indeed, using high angular resolution 1.3 mm observations, they resolved an asymmetry in the continuum emission of IRAS04302, finding that the east side is

brighter than the west side. They concluded that the east side traces the backside of the disk at millimeter wavelengths. Besides, the blueshifted outflow and millimeter molecular line enhancement have been observed toward the east in a number of previous studies (Podio et al. 2020; van't Hoff et al. 2020; Lin et al. 2023). These observations are consistent with a moderate tilt and east side of the outer disk tilted toward us.

We note, however, that our toy model does not match every observational feature. On the one hand, the modeled emission is too radially extended at 21 µm compared to the data, suggesting that forward scattering is more important in the data than in the model. Future modeling will be required to infer the dust properties of these systems (R. Tazaki et al. 2023, in preparation).

On the other hand, while we show that the current tilted model observed at an azimuth of 270° is able to reproduce the brightest nebula switch between 12.8 and 21 µm, our analysis in Section 3.1 also found that IRAS04302 has some degree of lateral asymmetry. This suggests that the disk is not viewed perfectly from an azimuth of 270° but instead with some viewing orientation between 270° and 360°.

Finally, it is important to remind that at wavelengths shorter than ~7.7 µm the envelope is important in IRAS04302. No envelope is included in our models, which thus cannot be directly compared with the observations of this system at short wavelengths.

# 5. Discussion

## 5.1. Tilted Inner Regions in Protoplanetary Disks

### 5.1.1. Other Edge-on Disks and General Disk Population

The models presented in Section 4 showed that the impact of a tilted inner region is stronger at short wavelengths, where thermal emission is negligible. We found that for about half of the possible viewing orientations strong asymmetries are expected at 2 µm, while a flux reversal can occur for about a third of the possible orientations. These numbers may vary for different dust properties, inner disk sizes, and tilt geometry. Nevertheless, these predictions encourage us to analyze the diverse optical to near-infrared observations of edge-on protoplanetary disks, looking for any information on the occurrence of misaligned inner regions.

We focus here on the identifications of lateral asymmetries in previous observations. Indeed, even though flux reversal is predicted to occur over a relatively large range of azimuthal orientation, it is less readily detectable. This is because observations over a wide range of wavelengths (ideally from the optical/near-infrared to the mid-infrared) are required to observe the brightest nebula switch at a particular viewing orientation. As highlighted in Section 4.1 there are only 2/12 cases where brightest nebula switch was found between two wavelengths (IRAS40302, Flying Saucer).

Interestingly, a large number of scattered light observations of edge-on disks have revealed lateral asymmetries. To quantify the relative occurrence of lateral asymmetries, we analyzed the list of 22 edge-on disks compiled by Angelo et al. (2023, their Table 5), which gathers known edge-on disks within 150 pc observed at optical or infrared wavelengths. We excluded six systems (IRAS04158+2805, HH390, ISO-Oph31, SSTc2d J162221.0-230402, DG Tau B, CB 26), either because the scattered light observations show a unipolar nebula or because no scattered light observations were previously
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Left panels: JWST 12.8 and 21 µm observations of IRAS04302. Middle panels: tilted disk model at 12.8 and 21 µm, observed with an azimuth of 270°. Right panels: minor axis cut at the center of the disk for the data, tilted, and no tilt models, at 12.8 and 21 µm.</figcaption></figure>

published, and included four additional disks: PDS 144 N (Perrin et al. 2006), Orion114-426 (McCaughrean et al. 1998), Gomez Hamburger (Ruiz et al. 1987; Bujarrabal et al. 2008), and HOPS 136 (Fischer et al. 2014). The final sample thus contained 20 disks, observed at scattered light wavelengths within the range of 0.4–21 µm.

For all disks, we performed a direct measurement of $\delta_{\text{spines}}$. Because of lower dynamic range, for eight disks in the sample (HV Tau C, IRAS04200, HH30, Haro6-5B, Orion114-426, J1607, Flying Saucer, Oph EMM 3), we measured the spine size within 20%, 30%, or 40% of their respective peak, while for the rest of the sample we took the size within 10% of the peaks. As in Section 4.3, we identify disks with strong lateral asymmetries as those having $\delta_{\text{spines}} > 5\%$. We show the disk images, on which the spines are overplotted, in Figure 8, illustrating the prevalence of lateral asymmetries.

We find that 15/20 disks show clear lateral asymmetries, while 5/20 of them do not. Lateral asymmetries at short wavelengths are a strong prediction of our radiative transfer models with a tilted inner region. Thus, the large number of edge-on disks observed with lateral asymmetries suggests that a significant fraction of protoplanetary disks may possess a tilted inner region.

In recent years, a number of protoplanetary disk observations have revealed tilted or warped inner regions. This was discovered mainly by detecting shadows in scattered light (e.g., Benisty et al. 2018; Muro-Arena et al. 2020) or finding gas kinematic signatures (e.g., Pérez et al. 2018). Misalignment evidence is encountered both in young (e.g., Sakai et al. 2018; Yamato et al. 2023) and older protoplanetary disks (e.g., Wolff et al. 2016; Debes et al. 2017), and the broad shadows moving over a several year timescale in the TW Hya disk (Debes et al. 2023) suggests that they are potentially long lived.

Statistical analysis of their occurrence is still limited but several studies of samples of few tens of disks have started to

reveal a large fraction of misaligned systems. For example, Ansdell et al. (2020) analyzed a sample of 24 dipper objects, which are often assumed to host close to edge-on inner disks. They found that the outer disks are consistent with an isotropic distribution of inclinations, suggesting that most of these systems have misalignments. Besides, Bohn et al. (2022) analyzed the inclination of the inner and outer regions of a sample of 20 transition disks observed with VLT/GRAVITY and ALMA. They concluded that six disks were clearly misaligned, five were not, and were inconclusive about the nine remaining systems. The analysis presented here shows that lateral asymmetries are frequent in edge-on protoplanetary disks (15/20). Given that tilted inner disks are also common around other systems, it is realistic that at least some of the asymmetry seen in edge-on disks are induced by a tilt.

Finally, we note that, even though some systems were observed multiple times and at different wavelengths, we limited our analysis to one observation only for each object (single epoch and wavelength), meaning that time variability was not taken into account. Yet, Watson & Stapelfeldt (2007) showed that, for example in HH30, the lateral asymmetry varies rapidly over a 10 yr period. Other examples of variability are that the most recent 2023 observations of Tau042021 shows a symmetric disk (Duchêne et al. 2023), as opposed to the asymmetric classification in Figure 8, while other systems are more stable (e.g., HK Tau B). In the case of HH30, the angle of observation of the system (azimuth) is not expected to change significantly over a 10 yr timescale. If the lateral asymmetry is caused by a tilted inner disk scenario alone then it would imply that the inner disk has to precess very rapidly to explain the rapid changes in lateral differences, which might not be possible. Hence, more complete analysis and modeling would be needed to interpret the origin of the lateral asymmetry of each individual system. We discuss potential alternative
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Gallery of edge-on disks with scattered light observations where we estimate the level of lateral asymmetry. Asymmetric disks have $|δ_{spines}| ≥ 5\%$. The spines within 10% of their respective maxima, used for the measurement of $δ_{spines}$, are shown as blue lines. Their respective centers are marked by a point. A 1" scale bar is also displayed in the bottom right of each panel.</figcaption></figure>

explanations for the lateral asymmetries observed in edge-on disks in Section 5.1.2.

### 5.1.2. Tilt Origin and Alternative Explanations

Several mechanisms are thought to be able to lead to a warped disk. One of the first proposed mechanisms stated that misalignment of the rotation axis of the disk with the magnetic field direction can warp the innermost edge of the disk (AA Tau; Bouvier et al. 1999). Alternatively, binaries misaligned with the disk or an inclined planet orbiting the central star, if massive enough, might also be able to generate misalignments between the inner and outer regions of the disk (Facchini et al. 2013; Nealon et al. 2018; Young et al. 2023). Other scenarios involve the anisotropic accretion or late capture of infalling material with a different angular momentum vector orientation than that of the disk (Dullemond et al. 2019; Kuffmeier et al. 2021). Given the relatively young age of IRAS04302, it is not clear which scenario is most likely.

Alternative explanations for the asymmetries in the edge-on disks optical and near-infrared images have previously been

proposed. For example, hot spots in the stellar surface can allow us to reproduce the asymmetry in HST optical and near-infrared images of HH30 (Stapelfeldt et al. 1999; Cotera et al. 2001). Such a scenario would predict rapid asymmetry variations that may be detectable over timescales of a few days, as in the case of HH30. A hot spot might not lead to lateral asymmetries in the mid-infrared because their temperatures might be too hot, but a colder spot could. A test of this scenario with respect to the inner tilt scenario could include a variability study of the optical/near-infrared asymmetry (a tilted inner region would not rotate as fast as a stellar hot spot), and looking for the presence of a mid-infrared laterally asymmetry counterpart.

Another possibility to generate lateral asymmetries includes the presence of asymmetric envelope or dust cloud material near the source, but detailed models have not been performed to confirm these suggestions. Similarly, one could also imagine that the lateral asymmetries are due to an intrinsic morphological difference between the two sides of the disk (e.g., different flaring between the top/bottom or left/right sides). It
is, however, not clear how such significant morphological
differences could be long lived. For these scenarios, we would
expect the lateral asymmetries to remain similar at all
wavelengths, as opposed to predictions from the tilted models.

Finally, we note that the models with a tilted inner disk predict that lateral asymmetries are more pronounced and opposite between the near-infrared and the mid-infrared. Most disks were previously observed at optical to near-infrared wavelengths, and thus looking for asymmetries in these systems at 21 µm can allow us to test the misalignment hypothesis.

## 5.2. Dark Lane Thickness and Aspect Ratio Variation with Evolutionary Stage

The variation of the dark lane thickness of an edge-on disk is a hint toward its scale height and optical depth (Watson et al. 2007). One of the main goals of JWST GO program 2562 presented here (and upcoming cycle 2, JWST GO program 4290) is to use this measurement to obtain information on the level of dust settling and dust properties across different evolutionary stages. This could allow us to better understand the timescale of dust accumulation of the midplane, and ultimately of planet formation.

Here we aim to look for any trend for the variation of the dark lane thickness and of the disk aspect ratio at near- to mid-infrared wavelengths with evolutionary class. We focus on three edge-on disks, observed with JWST, in different evolutionary stages: the Class I IRAS04302 studied in this work, L1527 (LDN 1527, Class 0/I), and Tau042021 (2MASS J04202144+2813491, Class II). For L1527, we use archival 3.35, 4.44, 7.7, 12.8, and 18 µm JWST observations and determine the disk morphology at the different wavelengths in this disk as in Section 3.1 (see the Appendix). For Tau042021, we use the results obtained by Duchêne et al. (2023) with a similar methodology.

We consider both the dark lane thickness $d_{\text{neb}}$ and the disk aspect ratio of the different observations. The latter is defined as the ratio of $d_{\text{neb}}$ and the averaged disk size at 7.7 and 12.8 µm. The disk size at one wavelength is taken as the smallest value of $R_{\text{FW10\%}}$ reported in Table 1 for IRAS04302, Table 2 (in the Appendix) for L1527, and the value of $\text{FW10\%}_{\text{top}}$ in Table 2 of Duchêne et al. (2023) for Tau042021. The final adopted disk sizes are 2.34 (237 au, at 140 pc; van't Hoff et al. 2023), 2.72 (438 au, at 161 pc; Galli et al. 2019), and 2.63 (342 au, at 130 pc; Galli et al. 2019), respectively for L1527, IRAS04302, and Tau042021.

**Variation of the Dark Lane Thickness.** We show the variation of $d_{\text{neb}}$ with wavelength for these three objects in the top panel of Figure 9. First, we find that, although the disks have different radial extent, the separation of their nebulae appears to be very similar. Moreover, the dark lane thickness does not vary significantly over this wavelength range (typically less than 10%).

Duchêne et al. (2023) presented the first results of the series of observations of edge-on disks from JWST GO program 2562. They focused on the NIRCam and MIRI observations of Tau042021. In this system, they showed that the dark lane thickness does not vary significantly between 2 and 21 µm. By comparing the observations with predictions of a range of radiative transfer models, they concluded that such a small variation indicates that grains have grown to at least 10 µm and that such grains are not significantly settled. Although detailed

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Dark lane thickness (top) or aspect ratio (bottom) as a function of wavelength in L1527 (Class 0/I), IRAS04302 (Class I), and Tau042021 (Class II). The uncertainties correspond to ±0.1 for $d_{\text{neb}}$, propagated into the aspect ratio.</figcaption></figure>

radiative transfer would be required to take optical depth into account, the similar observational result obtained in IRAS04302 and L1527 suggests that these disks have a similar configuration, with some level of grain growth in their upper layers. Consistent results of unsettled ~10 µm dust grains in the outer regions of protoplanetary disks have also been reported in previous studies of other edge-on disks (e.g., Pontoppidan et al. 2007; Sturm et al. 2023).

At other locations than in the disk midplane, Franceschi et al. (2023) showed that the collisional timescale in protoplanetary disks is typically larger than disk lifetimes, preventing efficient dust growth in the upper layers of the disk. Finding intermediate grains in the upper layers of disks could thus suggest that either turbulence in these systems is sufficiently strong to replenish the upper layers with larger grains from the midplane (e.g., Tazaki et al. 2021; Sturm et al. 2023), or that dust particles are able to grow there through processes other than turbulent collision (e.g., coagulation driven by settling; Franceschi et al. 2023).

**Variation of the Aspect Ratio.** The bottom panel of Figure 9 shows the variation of the disk aspect ratio with wavelength, which allows us to reduce potential dependencies of the apparent height with disk radius. We find that the disk aspect ratio at JWST wavelengths tentatively decreases with evolutionary class. At all wavelengths, the young L1527 appears relatively thicker than IRAS04302, which is also thicker than the older Tau042021, although uncertainties are large.

The apparent height of the dust grains observed with JWST is linked to the total disk dust mass and the gas scale height, which itself depends on the disk temperature and stellar mass. Tau042021 and L1527 have a similar stellar mass, but L1527 has both a higher dust mass ($M_{\text{dust}} > 41 M_{\oplus}$) versus $M_{\text{dust}} > 25.5 M_{\oplus}$; Villenave et al. 2020; van't Hoff et al. 2023) and possibly higher disk temperature than Tau042021 (CO freeze out at $R > 350$ au versus $R > 100$ au; Duchêne et al. 2023; van't Hoff et al. 2023), which could explain its higher aspect ratio. On the other hand, the difference in aspect ratio between Tau042021 and IRAS04302 is possibly driven by the
larger dust mass in IRAS04302 ($M_{\text{dust}} > 54.8 M_{\oplus}$ versus $M_{\text{dust}} > 25.5 M_{\oplus}$; Villenave et al. 2020). Altogether the tentative decrease in disk aspect ratio with evolutionary class may be linked to changes in the disk temperature structure, dust mass, or dust properties with age (see, e.g., van't Hoff et al. 2020).

## 6. Conclusions

In this study, we have presented new JWST observations of the Class I IRAS04302+2247, obtained with both NIRCam and MIRI imagers. The source is detected in scattered light at all wavelengths. At 2 and 4.44 µm the disk appears surrounded by a prominent envelope, which gradually disappears at 7.7 µm, leaving clear disk images at 12.8 and 21 µm. We compared 25.5 year older HST observations with the new JWST 2.0 µm image and identified a scattered light blob feature which is moving slowly, directly away from the presumed star location. The analysis of its proper motion indicates that it is likely not a background source, and its low velocity suggests that it corresponds to material entrained by the low-velocity outflow.

We characterized the morphology of IRAS04302 and two other disks of different evolutionary classes (L1527 and Tau042021) and found no significant variation of dark lane thickness with wavelength (typically within ~10%). This is consistent with the presence of relatively large grains (several microns) unsettled in the upper layers of the disk (Duchêne et al. 2023). In addition, the aspect ratio of the disks seems to decrease with evolutionary stage, which could trace evolution of the disk temperature structure, disk mass, and/or grain properties with age.

The characterization of IRAS04302 also revealed that the disk appears asymmetric, with its western nebula more extended to the south than its eastern counterpart. Finally, we find that the brightest nebula switches side between 12.8 and 21 µm.

We performed radiative transfer models introducing a tilted inner region to explore this particular aspect of the observations. We find that a tilted inner region can lead to a switch in the brightest nebula for about a third of the possible viewing orientation, and that such a model is able to reproduce the MIRI observations of IRAS04302. We discuss the viewing orientation of the system based on the modeling.

In addition, for some orientations, the model with a tilted inner disk predicts strong lateral asymmetries in edge-on disks. The intensity and orientation of such asymmetry depend on the wavelength, and are opposite between the near- and the mid-infrared. At short wavelengths (2 µm) about half of all possible viewing orientations are susceptible to create strong lateral asymmetries. We identify 15/20 edge-on disks in the literature with such lateral asymmetries. While stellar spots may also create lateral asymmetries, the frequency of asymmetries suggests that a large fraction of protoplanetary disks might possess misaligned inner regions, consistent with previous observational results. Variability studies of edge-on disks could allow us to test the dominant scenario. Indeed, stellar spots would lead to asymmetric variations over a timescale of several days, while a tilted inner disk might not lead to variation over a year to decade timescale. Observations at longer wavelengths (21 µm), such as enabled by the upcoming JWST GO program 4290, could also allow us to confirm the presence of a tilted inner region in these asymmetric disks. Indeed, the model

predicts that asymmetries would be opposite at 21 µm than at shorter wavelengths.

## Acknowledgments

We thank the anonymous referee for providing a detailed review of our work. We thank Howard Bushouse and Ken MacDonald at STScI for fixing the pixel fitting issues in the JWST pipeline and allowing to recover the saturated pixels in the F444W image. We are grateful to the W. M. Keck Institute for Space Studies as the hosts for our week-long team meeting in May 2023. 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). M.V., K.R.S., G.D., and S.G.W. 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. The research of M.V. was supported by an appointment to the NASA Postdoctoral Program at the NASA Jet Propulsion Laboratory, administered by Oak Ridge Associated Universities, under contract with NASA. R.T. acknowledges financial support from a CNES postdoctoral fellowship. The IRAS04302+2247 data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via 10.17909/kh9r-ym63.

Facilities: JWST, HST.

Software: mcfost (Pinte et al. 2006, 2009), Matplotlib (Hunter 2007), Numpy (Harris et al. 2020).

# Appendix
## JWST Observations of L1527

In Section 5.2, we use JWST NIRCam and MIRI data of L1527 observed as part of DD program 2739 (PI: K. Pontoppidan). We consider observations at five wavelengths, taken with filters F335M (3.35 µm), F444W (4.44 µm), F770W (7.7 µm), F1280W (12.8 µm), and F1800W (18 µm). The MIRI data were observed successively on 2022 September 1, with about 45 minutes exposures per filter. On the other hand, NIRCam observed the source on 2023 September 8 for 1460 s.

We obtained the phase 2 pipeline calibrated data from the MAST archive, and we re-ran phase 3, using the JWST pipeline version 1.9.5. Similarly to the reduction of IRAS04302 presented in Section 2, we set subtract=True in the skymatch step to subtract the background. The average background levels were 0.26, 0.53, 23.0, 85.08, and 170.23 MJy sr⁻¹ at F335W, F444W, F770W, F1280W, and F1800W, respectively.

The final images are shown in Figure 10. The disk is vertically and radially resolved up to 18 µm. Except for the observations at 18 µm all other images show extended scattered light along the minor axis direction. This feature is connected to the spectacular outflow cavity revealed by the NIRCam observations.<fnref n="8" /> Moreover, we see that the disk is not laterally symmetric. Specifically, the top nebula appears to extend toward the left while the bottom nebula is wider to the right. This asymmetry may be due to a tilted inner region as we argue

<footnote n="8">⁸ https://nasa.gov/feature/goddard/2022/nasa-s-webb-catches-fiery-hourglass-as-new-star-forms</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Gallery of the JWST observations of L1527. All images are shown with a log stretch. The bottom left ellipse in each panel indicates the corresponding angular resolution of the observations.</figcaption></figure>

<table><caption>**Table 2**
L1527 Morphological Properties</caption><thead><tr><th>λ (μm)</th><th>F<sub>5"&times;5"</sub> (mJy)</th><th>d<sub>neb</sub> (arcsec)</th><th>FR T/B</th><th>R<sub>FWHM</sub> T</th><th>R<sub>FWHM</sub> B</th><th>R<sub>FW10%</sub> T</th><th>R<sub>FW10%</sub> B</th><th>δ<sub>spines</sub> T-B (%)</th></tr></thead><tbody><tr><td>3.35</td><td>0.3 ± 0.03</td><td>1.09 ± 0.10</td><td>0.41 ± 0.05</td><td>1.18 ± 0.10</td><td>1.81 ± 0.17</td><td>...</td><td>...</td><td>...</td></tr><tr><td>4.44</td><td>8 ± 1</td><td>0.96 ± 0.10</td><td>0.50 ± 0.04</td><td>0.92 ± 0.10</td><td>1.63 ± 0.10</td><td>2.18 ± 0.10</td><td>2.63 ± 0.10</td><td>-9</td></tr><tr><td>7.7</td><td>13 ± 1</td><td>0.95 ± 0.10</td><td>0.64 ± 0.01</td><td>1.18 ± 0.10</td><td>1.42 ± 0.10</td><td>2.22 ± 0.10</td><td>2.45 ± 0.10</td><td>-13</td></tr><tr><td>12.8</td><td>10 ± 1</td><td>0.99 ± 0.10</td><td>0.51 ± 0.03</td><td>1.15 ± 0.10</td><td>1.37 ± 0.10</td><td>2.46 ± 0.11</td><td>2.67 ± 0.13</td><td>-12</td></tr><tr><td>18.0</td><td>72 ± 7</td><td>0.91 ± 0.10</td><td>0.29 ± 0.10</td><td>1.24 ± 0.10</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Note.** Same parameters as in Table 1. A negative value of δ<sub>spines</sub> indicates that the bottom nebula is centered further west than the top nebula. The error bars of d<sub>neb</sub>, R<sub>FWHM</sub>, and R<sub>FW10%</sub> correspond to the statistical error for the different spine averaging or 0<sup>′′</sup> 1, whichever is greater. “T” stands for the top nebula, and “B” for bottom nebula in Figure 10.

in Section 5.1.1. Interestingly, previous millimeter continuum and large-scale <sup>12</sup>CO observations also suggested the presence of a misaligned inner region in this system (Sakai et al. 2018; van’t Hoff et al. 2023).

We measure the flux within a squared aperture of 5". Then, we apply the same methodology as in Section 3.1 to determine the morphological features of the disk, namely dark lane thickness (d<sub>neb</sub>), the peak flux ratio of the bottom over the top nebula (FR), the radial R<sub>FWHM</sub> and R<sub>FW10%</sub> of the nebulae, and δ<sub>spines</sub> to quantify the lateral asymmetry. The results are summarized in Table 2.

## ORCID iDs

Marion Villenave https://orcid.org/0000-0002-8962-448X

Karl R. Stapelfeldt https://orcid.org/0000-0002-2805-7338

Gaspard Duchêne https://orcid.org/0000-0002-5092-6464

François Ménard https://orcid.org/0000-0002-1637-7393

Schuyler G. Wolff https://orcid.org/0000-0002-9977-8255

Marshall D. Perrin https://orcid.org/0000-0002-3191-8151

Christophe Pinte https://orcid.org/0000-0001-5907-5179

Ryo Tazaki https://orcid.org/0000-0003-1451-6836

Deborah L. Padgett https://orcid.org/0000-0001-5334-5107

## References

Angelo, I., Duchene, G., Stapelfeldt, K., et al. 2023, ApJ, 945, 130

Ansdell, M., Gaidos, E., Hedges, C., et al. 2020, MNRAS, 492, 572

Benisty, M., Juhász, A., Facchini, S., et al. 2018, A&A, 619, A171

Bohn, A. J., Benisty, M., Perraut, K., et al. 2022, A&A, 658, A183

Bouvier, J., Chelli, A., Allain, S., et al. 1999, A&A, 349, 619

Bujarrabal, V., Young, K., & Fong, D. 2008, A&A, 483, 839

Chauvin, G., Ménard, F., Fusco, T., et al. 2002, A&A, 394, 949

Cotera, A. S., Whitney, B. A., Young, E., et al. 2001, ApJ, 556, 958

Debes, J., Nealon, R., Alexander, R., et al. 2023, ApJ, 948, 36

Debes, J. H., Poteet, C. A., Jang-Condell, H., et al. 2017, ApJ, 835, 205

Doi, K., & Kataoka, A. 2021, ApJ, 912, 164

Duchêne, G., Menard, F., Stapelfeldt, K., et al. 2023, arXiv:2309.07040

Dullemond, C. P., Küffmeier, M., Goicovic, F., et al. 2019, A&A, 628, A20

Facchini, S., Lodato, G., & Price, D. J. 2013, MNRAS, 433, 2142

Fischer, W. J., Megeath, S. T., Tobin, J. J., et al. 2014, ApJ, 781, 123

Franceschi, R., Birnstiel, T., Henning, T., & Sharma, A. 2023, A&A, 671, A125

Galli, P. A. B., Loinard, L., Bouy, H., et al. 2019, A&A, 630, A137

Gáspár, A., Rieke, G. H., Guillard, P., et al. 2021, PASP, 133, 014504

Gräfe, C., Wolf, S., Guilloteau, S., et al. 2013, A&A, 553, A69

Grosso, N., Alves, J., Wood, K., et al. 2003, ApJ, 586, 296

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hunter, J. D. 2007, CSE, 9, 90

Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117

Krist, J. E., Hook, R. N., & Stoehr, F. 2011, Proc. SPIE, 8127, 81270J

Krist, J. E., Stapelfeldt, K. R., Burrows, C. J., et al. 1998, ApJ, 501, 841

Kuffmeier, M., Dullemond, C. P., Reissl, S., & Goicovic, F. G. 2021, A&A, 656, A161

Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32

Lin, Z.-Y. D., Li, Z.-Y., Tobin, J. J., et al. 2023, ApJ, 951, 9

Liu, Y., Bertrang, G. H. M., Flock, M., et al. 2022, SCPMA, 65, 129511

Lucas, P. W., & Roche, P. F. 1997, MNRAS, 286, 895

McCabe, C., Duchêne, G., Pinte, C., et al. 2011, ApJ, 727, 90

McCaughrean, M. J., Chen, H., Bally, J., et al. 1998, ApJL, 492, L157

Muro-Arena, G. A., Benisty, M., Ginski, C., et al. 2020, A&A, 635, A121

Nealon, R., Dipierro, G., Alexander, R., Martin, R. G., & Nixon, C. 2018, MNRAS, 481, 20

Nealon, R., Pinte, C., Alexander, R., Mentiplay, D., & Dipierro, G. 2019, MNRAS, 484, 4951

Ohashi, N., Tobin, J. J., Jørgensen, J. K., et al. 2023, ApJ, 951, 8

Padgett, D. L., Brandner, W., Stapelfeldt, K. R., et al. 1999, AJ, 117, 1490

Pascucci, I., Cabrit, S., Edwards, S., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 567

Pérez, L. M., Benisty, M., Andrews, S. M., et al. 2018, ApJL, 869, L50

Perrin, M. D., Duchêne, G., Kalas, P., & Graham, J. R. 2006, ApJ, 645, 1272

Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, Proc. SPIE, 9143, 91433X

Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967
Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

Pizzati, E., Rosotti, G. P., & Tabone, B. 2023, MNRAS, 524, 3184

Podio, L., Garufi, A., Codella, C., et al. 2020, A&A, 642, L7

Pontoppidan, K. M., Stapelfeldt, K. R., Blake, G. A., van Dishoeck, E. F., & Dullemond, C. P. 2007, ApJL, 658, L111

Ruiz, M. T., Blanco, V., Maza, J., et al. 1987, ApJL, 316, L21

Sakai, N., Hanawa, T., Zhang, Y., et al. 2018, Natur, 565, 206

Stapelfeldt, K. R., Ménard, F., Watson, A. M., et al. 2003, ApJ, 589, 410

Stapelfeldt, K. R., Watson, A. M., Krist, J. E., et al. 1999, ApJL, 516, L95

Sturm, J. A., McClure, M. K., Law, C. J., et al. 2023, A&A, 677, A17

Tazaki, R., Murakawa, K., Muto, T., Honda, M., & Inoue, A. K. 2021, ApJ, 921, 173

van't Hoff, M. L. R., Harsono, D., Tobin, J. J., et al. 2020, ApJ, 901, 166

van't Hoff, M. L. R., Tobin, J. J., Li, Z.-Y., et al. 2023, ApJ, 951, 10

Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, 651, A72

Villenave, M., Ménard, F., Dent, W. R. F., et al. 2020, A&A, 642, A164

Villenave, M., Podio, L., Duchêne, G., et al. 2023, ApJ, 946, 70

Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11

Watson, A. M., & Stapelfeldt, K. R. 2007, AJ, 133, 845

Watson, A. M., Stapelfeldt, K. R., Wood, K., & Ménard, F. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson, AZ: Univ. Arizona Press), 523

Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409

Wolf, S., Schegerer, A., Beuther, H., Padgett, D. L., & Stapelfeldt, K. R. 2008, ApJL, 674, L101

Wolff, S. G., Duchêne, G., Stapelfeldt, K. R., et al. 2021, AJ, 161, 238

Wolff, S. G., Perrin, M., Millar-Blanchaer, M. A., et al. 2016, ApJL, 818, L15

Wolff, S. G., Perrin, M. D., Stapelfeldt, K., et al. 2017, ApJ, 851, 56

Yamato, Y., Aikawa, Y., Ohashi, N., et al. 2023, ApJ, 951, 11

Youdin, A. N., & Lithwick, Y. 2007, Icar, 192, 588

Young, A. K., Stevenson, S., Nixon, C. J., & Rice, K. 2023, MNRAS, 525, 2616