© The Authors 2025

1 Introduction

Understanding the physical structure of protoplanetary disks is of critical importance to understanding planet formation. Historically, protoplanetary disks have been assumed to be axisymmetric, where all material is distributed along the same plane. However, a growing number of observations are revealing the presence of asymmetric features, such as spirals, shadows, and other lateral asymmetries (e.g., reviews by Benisty et al. 2023; Bae et al. 2023). Some of these asymmetries can be indications of a misaligned inner region or warp casting a shadow onto the outer part of the disk. Misalignments or warps might in fact be fairly common in protoplanetary disks. A local misalignment (i.e., broken disk) or warp (i.e., smooth misalignment) will locally perturb the planet formation process and can be a cause of the diversity of planets.

Warps can be formed in many different scenarios. One scenario is a bound companion object, either a planet (Facchini et al. 2014) or a binary companion that is inclined with respect to the disk plane. This could either be an outer companion (Zanazzi & Lai 2018a), a system with a circumbinary disk (Facchini et al. 2013; Lodato & Facchini 2013; Zanazzi & Lai 2018b; Deng & Ogilvie 2022; Rabago et al. 2023; Young et al. 2023), or any combination of the two. This companion object exerts a gravitational torque on the disk. This torque causes the concentric rings within the disk to precess. Because the torque strongly depends on the distance to the gravitational perturber, the concentric rings precess on different timescales. This produces the warp shape. The same principle applies to a flyby of a massive object inclined with respect to the disk (Picogna & Marzari 2014; Dai et al. 2015; Xiang-Gruess 2016; Cuello et al. 2019, 2023, Kimmig et al., in prep.). Such flybys can be common, especially in the early phases of star and disk formation (Pfalzner & Govind 2021). A further formation scenario is the late infall of material, which can add angular momentum in a different direction, thus distorting the disk (Dullemond et al. 2019; Ginski et al. 2021; Kuffmeier et al. 2021; Krieger et al. 2024). Finally, a change in the rotation axis during the formation process due to a magnetic field is also thought to cause misalignments (AA Tau; Bouvier et al. 1999).

The diversity of formation scenarios already hints at the variety of appearances of misalignments in observations. Additionally, warps are highly dynamic (Papaloizou & Pringle 1983; Papaloizou & Lin 1995). Once the warp is formed through one of the formation scenarios, internal torques come into play. These internal torques arise due to the misalignment of neighboring orbits, which leads to pressure gradients in the disk. The pressure gradients then cause resonant motions of the gas, called sloshing motions (Dullemond et al. 2022), and vertical breathing modes (Held & Ogilvie 2024). These motions exert a torque on the different orbital planes within the disk, and hence change the orientation of the orbits. Thus, the net effect of the internal torque and pressure gradients is the evolution of the warp shape.

In typical conditions of a protoplanetary disk (especially low viscosity), the warp travels as a wave through the disk, often referred to as the wave-like regime (Lubow & Ogilvie 2000; Gammie et al. 2000; Ogilvie & Latter 2013a,b). Its wave speed is related to the sound speed in the disk (Nixon & King 2016). The wave is dampened over time, mainly by viscosity, on a damping timescale that can be estimated by t_damp ≈ 1/(αΩ_K), with Ω_K as Kepler frequency evaluated at the outer regime of the disk (Lubow & Ogilvie 2000). Here, α is the Shakura-Sunyaev viscosity parameter (Shakura & Sunyaev 1973). This timescale typically ranges from 10³ to 10⁵ years (Foucart & Lai 2014; Poon et al. 2021; Rowther et al. 2022; Kimmig & Dullemond 2024). During their evolution, warps may also break into separate disks under certain conditions (Nixon et al. 2013; Doğan et al. 2023; Young et al. 2023; Rabago et al. 2024).

Because of their dynamic evolution, warps can take on various shapes leading to a plethora of disk misalignment observations. Asymmetric observations linked to shadows due to warps are recurrent in scattered light imaging of low- to midinclination protoplanetary disks (e.g., Garufi et al. 2018, 2022; Benisty et al. 2023). Low-inclination systems can show narrow dark lanes (Benisty et al. 2017; Pinilla et al. 2018; Keppler et al. 2020), broad dark regions (Garufi et al. 2016; Debes et al. 2017; Muro-Arena et al. 2020; Kraus et al. 2020; Hashimoto et al. 2024), or a combination of both (Stolker et al. 2016; Benisty et al. 2018). Furthermore, the sensitivity and spectral resolution of the Atacama Large Millimeter Array (ALMA) are now revealing a larger number of warped disks, for example via gas kinematic observations (e.g., Casassus et al. 2015; Loomis et al. 2017; van der Plas et al. 2017; Pérez et al. 2018; Sakai et al. 2019; Yamato et al. 2023).

In parallel, an increasing number of studies also aim to directly compare the orientation of the inner and the outer regions of protoplanetary disks. Francis & van der Marel (2020) compared the inclination and position angle of the inner disk (≲20 au) and main ring of 14 transition disks observed at high angular resolution with ALMA, finding that 8 of them are significantly misaligned. Ansdell et al. (2020) instead analyzed a sample of 24 dipper stars, for which the inner disk is believed to have an inclination close to edge-on, and found that the outer disks follow a random distribution of orientations, indicative of significant misalignment for these systems. Bohn et al. (2022) compared the inclination of the outer disk, based on ALMA observations, with that of the inner disk, obtained using VLT/Gravity (Very Large Telescope) for a subsample of 20 transition disks. They identify 6 of the 20 disks that have clear evidence of misalignments, all of which show shadows on their outer disks at scattered light wavelengths, compatible with the derived inner disk geometry. Kluska et al. (2020) use reconstructed images of near-infrared interferometry observations to evaluate the symmetry of the innermost region of 15 protoplanetary disks around Herbig AeBe stars. They find evidence of non-axisymmetric emission that can be linked to warps in 5 of these 15 disks.

A key to interpreting the plethora of misalignment observations are models for identifying signatures of warped disks. Such models are able to provide kinematic insight (Casassus et al. 2015; Young et al. 2022; Zuleta et al. 2024) or chemical signatures (Young et al. 2021) in disks containing warps. Shadows can be investigated in synthetic observations in scattered light (Nealon et al. 2019; Juhász & Facchini 2017). Facchini et al. (2018) investigated shadow signatures of broken disks in intermediate inclination disks and were able to reproduce both dark lanes for strong misalignments, and broad regions for low to moderate misalignments. The location of the observed shadows can be used to constrain the geometry of the inner, misaligned disk (Min et al. 2017). The shadows can be variable, for example for precessing broken disks (Nealon et al. 2020).

Asymmetries can also occur in disks observed edge-on (Villenave et al. 2024). Edge-on disks typically appear in scattered light observations as two nebulae separated by a dark line (e.g., Burrows et al. 1996; Stapelfeldt et al. 1998). Those two nebulae correspond to the surface regions of the disk, while the dark line occurs because the material close to the midplane obscures the light from the host star. Thus, the characteristics of the dark line strongly depend on the properties and vertical distribution of the dust contained in the disk (Watson et al. 2007). On the other hand, the surfaces of the disk are less dense, and therefore multiple scattering events allow photons to escape the disk and reach the observer, thus leading to the two bright nebulae. In a perfectly planar, unwarped disk, viewed perfectly edge-on, the nebulae would appear entirely symmetric. Even though observing such a disk with a slight inclination variation would lead to one of the nebulae being brighter due to forward scattering, the extent of both nebulae along their major axis would still be symmetric. However, observations show that a large fraction of edge-on protoplanetary disks at scattered light wavelengths display a shift or lateral asymmetry of their nebulae (15 of the 20 disks analyzed by Villenave et al. 2024, see also Berghea et al. 2024). A misalignment of the inner disk region is proposed to explain this asymmetry, as reproduced in some models (Facchini et al. 2018; Nealon et al. 2019; Villenave et al. 2024).

In this work, we specifically target the scattered light appearance of edge-on and almost edge-on disks containing a warp. We aim to investigate how the asymmetry of the two nebulae changes depending on the warp morphology and to find limits of observability. In addition, we also evaluated the warp geometry that can lead to a switch in the brightest nebula with wavelength, which has been observed in a few systems (Grosso et al. 2003; Villenave et al. 2024). We present our models in Sect. 2 and the results in Sect. 3. In Sect. 4, we link our results to observations by evaluating the misalignment between the outer regions of edge-on disks and jets presumably linked to the inner region. We conclude our work in Sect. 5.

2 Methods

To study the appearance of warped edge-on protoplanetary disks, we performed radiative transfer simulations using RADMC-3D (Dullemond et al. 2012). We created most models in a wavelength of λ = 0.8 μm, but extended it to λ = 2 μm, 12 μm, and 21 μm for some parameter sets. Images in this near- to mid-infrared regime are often called scattered light images as scattered light usually dominates the appearance of disks in these wavelengths. In the longer, mid-infrared regime wavelengths thermal emission can contribute to the radiation from disks, and we include the thermal emission in all of our models. However, even with the thermal emission, the scattered light scattered light dominates the outcome of the images. The images are scaled by 2.28 au per pixel, and for the evaluation using spine fitting, we assumed a distance of 140 pc.

2.1 Model setup

We set up an initial gas density structure, where the surface density profile follows a power law with an exponential cut-off both at the inner and outer disk edge, according to $$\mathrm{\Sigma }(r)={\mathrm{\Sigma }}_0{\left(\frac{r}{R_0}\right)}^{-p}\left(\frac{1}{1+\exp \left(\frac{r-r_{\text{out}}}{0.05r_{\text{out\hspace{0.17em}}}}\right)}\right)\left(\frac{1}{1+\exp \left(\frac{r_{\text{in}}-r}{0.05r_{\text{in}}}\right)}\right),$$(1)

with r_in and r_out as inner and outer disk edge, respectively, and Σ₀ as surface density at the reference radius R₀, which we set to R₀ = 5.2 au. According to typical temperature structures in protoplanetary disks, the vertical extent can be described by $$h(r)=\frac{H_{\mathrm{p}}}{r}=h_0{\left(\frac{r}{R_0}\right)}^{i_{\mathrm{f}\mathrm{l}}},$$(2)

where h₀ is the aspect ratio at the reference radius R₀ and i_fl is the flaring index. To parametrize the warp, we implemented a tilt profile inspired by Martin et al. (2019) of $$i(r)=i_{max}[\frac{1}{2}\tanh \left(\frac{r-r_{\text{warp\hspace{0.17em}}}}{r_{\text{width\hspace{0.17em}}}}\right)-\frac{1}{2}],$$(3)

with i_max as maximum warp tilt, r_warp the warp location, and r_width is the width of the warp transition. This setup is identical to Kimmig & Dullemond (2024), with adjusted parameters.

In all our models, the disk ranges from r_in = 2 au to r_out = 200 au, while the grid extends from 0.5 au to 260 au in the radial direction. We used a full spherical grid (i.e., the vertical extent is set from 0 to π). For the resolution, we chose 250 grid cells in the radial direction, 100 cells in the azimuthal direction, and 250 cells in the vertical direction. We set the aspect ratio to h₀ = 0.1 and use a flaring disk geometry with i_fl = 0.1. Similar parameters have been obtained from the modeling of scattered light images of several disks (e.g., Muro-Arena et al. 2018; Sturm et al. 2023, see also Angelo et al. 2023).

We chose for the star’s properties a mass of M_* = 2.5 M_⊙, a radius of R_* = 1.8 R_⊙, and an effective temperature of T = 10 500 K. For the gas distribution, we set a slope of p = 1 and Σ₀ = 34.3 g/cm² at the reference radius of R = 5.2 au, resulting in a gas disk mass of M_gas = 0.01 M_*.

To investigate the effect of differently warped disks, we varied the warp tilt i_max from 0° to 10° in steps of 1°. For the location of the warp, we set a fiducial value of r_warp = 20 au, but varied the location of the warp for a warp tilt of 3° to r_warp ∈ {5 au, 50 au, 100 au}. To get a smooth warp transition at all locations, we scaled the warp width with the location of the warp using r_width = 0.25 r_warp. We can measure the warp strength in each case by evaluating the maximum displacement of the disk midplane in terms of the pressure scale height at the location of the warp. Equation (3) shows that the tilt angle between the disk midplane at the warp location r_warp with respect to the inner (as well as outer) disk midplane is i_max/2. We can then calculate the displacement w from the height of the midplane z_mid(r_warp) with respect to the midplane of the inner (or outer) disk $$w=\frac{z_{\text{mid\hspace{0.17em}}}\left(r_{\text{warp\hspace{0.17em}}}\right)}{H_{\mathrm{p}}\left(r_{\text{warp\hspace{0.17em}}}\right)}=\frac{r_{\text{warp\hspace{0.17em}}}\tan \left(\frac{i_{\text{max\hspace{0.17em}}}}{2}\right)}{H_{\mathrm{p}}\left(r_{\text{warp\hspace{0.17em}}}\right)}.$$(4)

The numbers of w for the cases in our investigation are given in Table 1 and Fig. 1 shows vertical cross sections through the gas density for three example cases with different warp tilts (for r_warp = 20 au).

Fig. 1 Vertical cross section of the gas density setup for an unwarped disk, a warp of 3°, and a warp of 10°. The dash-dotted line indicates the disk midplane.

2.2 Radiative transfer

For the radiative transfer simulations, we need to include dust in our models. We assumed small dust particles that are perfectly coupled to the gas and well-mixed, which means that the structure of the dust perfectly follows the structure of the gas disk. Then, we chose a dust-to-gas ratio of 10⁻³, leading to a total dust mass of M_dust = 10⁻⁵ M_*. We note that this dust-to-gas ratio is slightly lower than typically assumed because small dust particles likely only make up a fraction of the total dust mass in protoplanetary disks. Despite recent studies on the behavior of larger dust particles in warped disks (Cuello et al. 2019; Aly et al. 2021, 2024), many aspects remain unexplored, which is why we decided to only consider a single grain size of a = 10 μm. This is the simplest possible description, in contrast to a dust size distribution, as such distributions are not yet well constrained for protoplanetary disks. The assumption that dust of this size is well coupled to the gas is supported by observations, at least in some systems (Pontoppidan et al. 2007; Sturm et al. 2023; Duchêne et al. 2024; Villenave et al. 2024). We set a composition of carbon and silicate (pyroxene with 70% magnesium and 30% iron) with a mass ratio of carbon to pyroxene of 0.87/0.13, and a porosity of 25%, according to the DIANA standard dust model (Woitke et al. 2016). We created the corresponding opacity using optool (Dominik et al. 2021). The resulting scattering and absorption opacities are shown in Fig. A.1.

To make sure that the qualitative results still hold for a different choice of dust grain sizes, we explored grain sizes of a = 0.5 μm and 1 μm, as well as a dust size distribution in Appendix C. Overall, we find similar results for the lateral asymmetries and flux ratios in all of these models.

The scattering phase functions in our setup are very forward-peaked, especially for wavelengths up to a few micron. This can lead to very expensive simulations. However, we can use a technique to save computation time, called chopping. In this technique, we chop off a few degrees of the scattering phase function and treat rays within that cone as if they were not interacting with the scattering particle. This technique helps to avoid bright spots that appear in the radiative transfer models if not enough scattering photons are used. We carefully checked the applicability of this technique in our case in Appendix A. For the Monte Carlo temperature calculations, we use 10⁸ photon packages. For the near-/mid-infrared radiative transfer images, we used 10⁹ scattering photon packages for all images at the wavelengths of 0.8 μm, 2 μm, and 21 μm. We found that at 12 μm, the same amount of photon packages leads to noisy images, which is why we used twice as many photon packages for this wavelength.

Because warped disks are not axisymmetric, we need to synthetically observe the disk from different directions. For that, we define an angle φ, so that φ = 0° looks at the warp tilt from the side, as illustrated in Fig. 2. We can exploit some symmetries in the disk, as the disks in our models are untwisted. Therefore, we only need to consider a quarter of the φ-space if the outer disk is viewed perfectly edge-on (inclination of exactly i = 90°). The remaining regime of φ can then be achieved by appropriate rotations of the images. In some models, we chose to look at the outer disk slightly inclined (i = 85°). In this case, we need to cover half of the azimuthal regime, as a part of the symmetry is broken.

Fig. 2 Definition of the azimuthal φ-angles. The colors represent the height above the x–y-plane: purple is in plane, blue-green above, and red-yellow below. If the outer disk is viewed perfectly edge-on, the quadrant of 270°–0° is identical to the quadrant 0°–90°, only flipped upside down. The sector 90°–270° is identical to the points of sector 270°–90° appropriately rotated and flipped. We chose this definition of φ to be consistent with Villenave et al. (2024).

Fig. 3 Schematic representation of the spine fitting. The light circles show the location of the maxima after a high-order polynomial fit of the minor-axis profiles. The solid lines represent the final spines after fitting the maxima by a polynomial of order 2. The diamond represent the locations where the flux reaches 10% with respect to the maximum of the spine, and the squares are the center of the spine (i.e., the mid-point of these two points). The lateral asymmetry is the ratio between the distance of the centers to averaged disk size.

2.3 Analysis

Images of edge-on protoplanetary disks in the mid- to near-infrared wavelength regime typically show two nebulae separated by an elongated dark lane. In the presence of a tilted or warped inner region, the illumination of the upper disk surfaces does not follow spherical symmetry and depends on the viewing orientation. This can lead to changes in the flux ratio between the two nebulae and/or to lateral asymmetries.

To quantify the lateral asymmetry and flux ratio in relation to different warp parameters, we follow the methodology employed by Duchêne et al. (2024) and following works on JWST observations. This method identifies a spine for each of the two scattering nebulae¹. To do so, the method first estimates the location of the two maxima along the minor-axis profiles for every location along the nebulae. The minor-axis profiles, averaged to a resolution of 0.1″ or 0.2″ (for a disk model at a distance of 140 pc), are fitted by a polynomial function of order 6 or 8, which is used to assess the location of the two maxima. Those are then fitted by a polynomial function of second order to determine the two spines. On both spines, we then identify the two points on each side of its maximum where the flux reaches 10% with respect to the maximum. We define the mid-point of these two points as the center of the spine, and define the distance between them as the disk size. Finally, following Villenave et al. (2024), we define the lateral asymmetry δ_spine as the distance between the center of both spines in the direction along the nebulae’s major axis, divided by the averaged disk size. The procedure described above is repeated multiple times for different radial extents of the scattering nebulae, between 0.9″ and 1.6″ with steps of 0.1", and for other parameters (averaged resolution, polynomial order) in order to estimate an error for the measurement of the lateral asymmetry. A schematic representation of the spine fitting is shown in Fig. 3.

In addition, independent of the spine fitting, we estimate the integrated flux ratio between the top and bottom nebulae. We do this by taking the ratio of the total flux included in the top nebula divided by the total flux of the bottom nebula. We find that this definition is the inverse of that of Villenave et al. (2024). Therefore, the flux ratio in our unwarped models is greater than one for a disk inclined at 85°, in contrast to theirs being smaller than one.

Fig. 4 Near-infrared radiative transfer models of disks with different amplitude warps, defined by imax, observed at λ = 0.8 μm at φ = 0° (top two rows). The bottom row shows models with a varying warp location or in other words a varying size of the misaligned inner region. For the color, we use a logarithmic scale.

3 Results

Protoplanetary disks in near-/mid-infrared show, if they are viewed edge-on (at an observer inclination close to 90°), two nebulae separated by a dark lane. In the following section (Sect. 3.1), we present the near-infrared images² at λ = 0.8 μm for different warp morphologies, as well as our analysis of the lateral asymmetry and the flux ratio between the two nebulae. In Sect. 3.2, we discuss images at multiple wavelengths. In Appendix B, we additionally show the models of an unwarped reference case, both perfectly edge-on and slightly inclined (85°).

3.1 Near-infrared images at a single wavelength

We present near-infrared images at λ = 0.8 μm of different warp amplitudes (characterized by i_max) in Fig. 4, top two rows. The bottom row in this figure shows models with different locations of the warp, or in other words, different sizes of the misaligned inner region. All disks in this figure are viewed perfectly edge-on and at an azimuthal angle of φ = 0°, meaning that the warp tilt is viewed from the side (see Fig. 2). For this viewing angle, the lateral asymmetry is greatest, as the shadowed regions of both nebulae are the most opposite.

The orientation of the warped disk with respect to the observer has a large effect on the visibility of the asymmetry. In Fig. 5, we show how this picture changes for different azimuthal orientations in the model with a misalignment between inner and outer disk of i_max = 3°. The top row shows the perfectly edge-on view (inclination of 90°), while the bottom row presents a disk inclined by 85°. The figure shows that the lateral asymmetry vanishes for azimuthal viewing angles shifted by 90° from zero, as the misaligned inner region is viewed face-on in these cases (see Fig. 2, bottom row).

To quantify the results of the near-infrared images, we evaluate the lateral asymmetry, as well as the flux ratio between the two nebulae as described in Sect. 2.3. Figures 6 and 7 present both quantities for our models with varying warp amplitude and warp location, respectively. In all those models, the outer disk is viewed perfectly edge-on. By eye, we define a lateral asymmetry of ≥10% to be visible and note that this is a somewhat arbitrary definition, but motivated by previously published work (see Villenave et al. 2024). The pink shaded region in Figs. 6 and 7 indicates the region below this visible threshold.

We find that the stronger the warp (higher i_max), the stronger the lateral asymmetry. This result is expected, since stronger misalignments can cast larger shadows. In Fig. 4, we see that the lateral asymmetry of a warped disk can be visible in scattered light with a misalignment of only 2°. This means that even slight warps can in principle be observable. However, we note that this is only the case for the best azimuthal viewing angle. For viewing angles within 30° around φ = 270° (and φ = 90°, respectively), even a stronger warp of 10° misalignment cannot produce a significant lateral asymmetry (see Fig. 6).

The difference in flux between the two nebulae, on the other hand, is greatest for a viewing angle of φ = 270° and, for an edge-on view on the outer disk, equals one if the warp is viewed directly from the side (φ = 0°). In observations, the flux ratio by itself is not a unique indication of a warp or misalignment, since slightly inclined planar disks also exhibit a difference in flux between the two nebulae. However, the variation of the flux ratio with wavelength can be a strong indicator of warps, which we explore in Sect. 3.2.

When changing the location of the warp (bottom row in Fig. 4), we can see a similar effect: the further out the warp, the stronger the asymmetry in the near-infrared images. However, this only holds true up to a certain point, as a warp location of 100 au shows slightly less asymmetry than a warp location of 50 au, see bottom row of Figs. 4 and 7. We suspect this to occur because of the way we set up the warp shape, as we scale the width of the warp with the warp location. This description is the most natural because it is coherent in logarithmic space. However, considered linearly, the warp transition is wider for warps located further out, leading to a less steep warp curve where more light can pass. This reduces the size of the shadowed region and causes less lateral asymmetry.

Additionally, we find the flux ratio for the case of r_warp = 100 au to be smaller than for some warps closer in at all azimuthal angles. This could again be linked to the setup of the warp: because the warp is close to the outer edge, the outer disk is small compared to the other models, which could influence both the lateral asymmetry and the flux ratio.

An interesting fact is that the image of the warped disk at φ = 270° with the outer disk viewed perfectly edge-on (top left image in Fig. 5) looks very similar to an unwarped disk that is inclined with respect to the observer, see Fig. B.1, bottom left panel. For observations, this means that it can be difficult to distinguish between an inclined planar disk and a perfectly edge-on warped system. We find that in our models, the near-infrared image of the inclined planar disk is more strongly peaked toward the star. However, this is unlikely to serve as distinction for observations, since the disk parameters are usually not known precisely.

Fig. 5 Near-infrared radiative transfer models (λ = 0.8 μm) of a disk with tilt imax = 3° observed from different azimuthal angles φ (color in log-scale). The top row shows a perfectly edge-on view, while the bottom row shows the disk observed at an inclination of 85°. We note that for symmetry reasons, we do not need to simulate φ = 90°. The image shown in this figure is the simulation of φ = 270° flipped upside down.

Fig. 6 Lateral asymmetry (top) and flux ratios (bottom) for warped disk models of different warp tilts imax (at λ = 0.8 μm). The outer disk is viewed perfectly edge-on (meaning an inclination of 90°). For symmetry reasons, we only need to model a quarter of the azimuthal regime in this configuration. The curves can be extended to the remaining regime by flipping the curves appropriately. The pink shaded region in the top panel corresponds to a lateral asymmetry of <10%, which is what we assume not to be visible by eye.

Fig. 7 Same as Fig. 6, but for different warp locations rwarp for a warp amplitude characterized by imax = 3°. Again due to symmetry reasons, the curves can be extended to the whole azimuthal regime. As in Fig. 6, the pink shaded region highlights the lateral asymmetry we assume not to be visible by eye.

3.2 Near-/mid-infrared radiative transfer models in multiple wavelengths

A few observations of protoplanetary disks show a swap of the brightest nebula with wavelength: In shorter wavelengths one nebula appears brighter, while the other nebula appears brighter in longer wavelengths. This phenomenon is rare, so far it has only been observed in three disks: IRAS04302+2247, HH30, and 2MASS J16281370-2431391 (Flying Saucer). The swap does not occur at the same wavelength for the observed disks. In IRAS04302 +2247, the swap occurs between λ = 12.8 μm and 21 μm (Villenave et al. 2024), while the switch for the Flying Saucer occurs between λ = 1.3 μm and 2.2 μm (Grosso et al. 2003). Finally, Tazaki et al. (2025) identified a change in the brightest nebula of HH30 between 7.7 μm and 12.8 μm. They mentioned intrinsic flux variability between epochs as a potential origin for a wavelength swap. However, with their observations being observed consecutively within 3h, they ruled out this possibility for this system. On the other hand, Villenave et al. (2024) were able to recreate this phenomenon in a broken disk with a misaligned inner region. In this work, we aim to investigate if a continuously warped disk exhibits a similar behavior.

For that, we present radiative transfer simulations in four different wavelengths: λ = 0.8 μm, 2 μm, 12 μm, and 21 μm. We tested the multi-wavelength behavior for two different warp tilts i_max = 3° and 10°. In these models, the outer disk is viewed at an inclination of 85°, because the phenomenon only occurs for slightly inclined disks, as discussed later in this section. We indeed find a swap in brightest nebula in specific cases. Fig. 8 shows an example of such a case, in especially, a warp tilt of 10° and an azimuthal viewing angle of φ = 60°. For the shortest wavelength, λ = 0.8 μm, the bottom nebula is brighter, while for λ = 21 μm, the top nebula clearly appears brighter.

The phenomenon of the switch between brightest nebula can be explained by the shadow from the misaligned inner region. First, we need to understand that in a general disk that is slightly inclined with respect to the line of sight, one of the nebulae appears brighter than the other. This is because forward scattering is dominant over backward scattering (see, e.g., Benisty et al. 2023, their Fig. 2), which makes the nebula closest to the observer brighter. Now, an inclined warped disk could be oriented in such a way that this brighter nebula is actually shadowed by the misalignment. In this case, the nebula actually appears darker. However, in longer wavelengths, the disk could be able to emit more thermal radiation, and therefore the shadow would be less pronounced. Thus, the nebula can again appear brighter. To test the hypothesis, we performed additional radiative transfer simulations with the scattering switched off (option scattering_mode_max=0 in RADMC-3D). We indeed find almost no emission for λ = 0.8 μm and 2 μm, indicating that the light in the scattering simulations indeed results purely from scattered light. In contrast, we find thermal emission in these test simulations in the wavelengths λ = 12 μm and 21 μm, where it is more significant for the latter. However, we want to point out that the scattered light still dominates the appearance of the disks at both of these wavelengths. In summary, the fact that the thermal emission becomes more significant for longer wavelengths means that the shadows in the scattered light become less pronounced, which enables the swap in brightest nebula with wavelength.

Evaluating the flux ratios quantitatively in Fig. 9, we see that for a φ-regime of between 30° and 90° (and due to symmetry also up to 120°), a flux ratio swap between wavelengths can be possible. This only occurs for the stronger warp of i_max = 10° in our models. For the moderate warp of i_max = 3°, the variation of the flux ratio with azimuthal angle is less strong at all wavelengths, so that the ratio never drops below one.

Concentrating on the disk containing a stronger warp, we see that the shape of the flux ratio curve differs with wavelength. Generally, for shorter wavelengths the curve is steeper and there is a stronger variation of the flux ratios with azimuth. Interestingly, a wavelength of λ = 2 μm shows more variation than λ = 0.8 μm. This could be due to a difference in the chopped scattering opacity and, consequently, a difference in albedo (see Fig. A.1). The flux ratio curve for λ = 21 μm is shallowest, leading to a flux ratio always above one. At the shorter wavelengths, the flux ratio drops below one for orientations where the misalignment is shadowing the otherwise brighter nebula. We see that this swap of brightest nebula can occur at different wavelengths, depending on the azimuthal orientation of the disk. Recall that we evaluate the flux ratio by using the integrated flux in our models. We additionally checked the flux ratios when considering only the peak intensity and found that the curves are qualitatively the same, but slightly noisier and shifted, so that only the curve of λ = 2 μm drops below one. However, the swap still occurs, and with slightly different parameters, such as a different inclination of the outer disk, a similar result can be achieved in peak intensity flux ratios.

In our models, we find that the azimuthal curve at 21 μm follows the same qualitative shape as the shorter wavelengths, with their maximum at φ = 270° and their minimum at φ = 90° (see Fig. 9). On the other hand, Villenave et al. (2024) find that the 21 μm curve is flipped with respect to the shorter wavelengths. This difference could lie in the different opacities and geometry of the models. Specifically, Villenave et al. (2024) investigated a broken disk, while we investigate a continuously warped disk. In a broken disk, the inner edge of the outer disk and the outer edge of the inner disk could thermally radiate onto the outer disk in longer wavelengths, and therefore leading to the inversion of the curve. Because in the models of this work, there are no gap edges, we expect this thermal emission to be less strong.

At this point, we want to note that we expect the phenomenon of the swap in brightest nebula with wavelength to occur rarely, because a very specific orientation of the observed disk is required. In especially, the outer disk needs to be slightly inclined with respect to the observer, and the azimuthal viewing orientation of the system is restricted to a specific regime with respect to the orientation of the misalignment. This expectation connects to the fact that, as of now, only three disks have been observed to exhibit this phenomenon. However, if observed, a change in the brightest nebula with wavelength can be a strong indicator of a misalignment between inner and outer disk.

Fig. 8 Near-/mid-infrared radiative transfer models of a disk with tilt imax = 10° at φ = 60° with an inclination of 85° observed in different wavelengths (color in log-scale). Considering the integrated flux, the bottom nebula is brighter than the top at shorter wavelengths, but at λ = 21 μm, the top nebula is brighter.

Fig. 9 Flux ratios for models of warped disks with a warp tilt of 10° (solid lines) and of 3° (light dashed lines) in four different wavelengths. The outer disk is oriented with an inclination of 85°. In this case, we need to calculate the models for half of the φ-regime, the other half is symmetric. The dashed black line indicates the line where the flux ratio equals one (where both nebulae have the same brightness). We note that the y-axis is plotted in log-scale.

3.3 Caveats

The findings in this work provide valuable insight, although they should be considered within their limitations, which we discuss in this section.

Most importantly, the results significantly depend on the disk model. The vertical thickness and flaring of the disk influence the extent of the shadows. In thinner disks with a lower aspect ratio, we expect a stronger lateral asymmetry, because a larger part of the outer disk’s surface can be hidden behind the misalignment. In thicker disks, on the other hand, the outer disk can be puffed up enough to bulge out of the shadow. A similar argument applies when considering differently flared disks, which occurs for different temperature structures. With more flaring, the surface can bulge out more, leading to less lateral asymmetry. In the same line of arguments, the size of the disk, in especially with respect to the location of the warp, can also influence the lateral asymmetry. Additionally, the warp parameterization has a strong impact on the shadowed region. As it is not well constrained yet, we studied one of the simplest warp shapes in this work: a warp that is only tilted in one direction. In general, there can be more complicated shapes, for example twisted disks or multiply warped or twisted disks, leading to more complex shadow signatures. Since the overall parameter space for disk and warp shapes is extremely large, a further exploration would go beyond the scope of this work.

Another factor influencing the results of the near-/mid-infrared images is the scattering opacity, which depends on the dust properties and total dust mass (Watson et al. 2007). A higher dust mass would lead to a higher opacity, and we therefore expect the dark lane between the nebulae to be wider, possibly influencing the lateral asymmetry. Different dust sizes or a size distribution can also strongly affect the opacity. In our models, we only considered small dust grains that are perfectly coupled with the gas. Although it is likely that larger dust particles are present in protoplanetary disks, we do not expect them to have a strong effect, as they are likely settled to the midplane (e.g., Villenave et al. 2020). We additionally recall our use of the so-called chopping mechanism of 2° to calculate the scattering matrix due to the forward peaked phase functions at small wavelengths. Although we thoroughly tested the applicability of this mechanism to our case in Appendix A, it remains an approximation. We find that we slightly overestimate the total flux by a few percent when using the chopping mechanism. Since we are mainly interested in the lateral asymmetry and flux ratios, this does not have a large effect on the results in this work.

We also note that the near-/mid-infrared images presented in this study are not convolved to a specific instrument resolution, in order to keep this study general as convolution may depend on specific telescope settings. As the images are fairly smooth, we do not expect a convolution to significantly change the results.

Lastly, we mention that our definition of the lateral asymmetry is somewhat arbitrary. We lean the definition on previously published work, but it can be defined differently, which could have an impact on the results. Our fitting method strongly depends on the parameters chosen, namely the radial range of the spines, the radial and vertical averaging of the minor-axis profiles, the order of the polynomial function fitting the minor-axis profiles, the assumed location for the center of the dark lane. In this work, we explored how varying these parameters will impact the results, reaching error bars between a few percent and tens of percents (see Figs. 6 and 7). However, due to the limited parameter space explored, the relative comparison between models is likely more robust than the absolute values of the lateral asymmetries.

4 Occurrence of warps in protoplanetary disks

In Sect. 3, we investigated the level of lateral asymmetry that can be expected for different warp amplitudes and locations, under our fiducial assumptions for dust opacities and warp geometry. Here, we aim to link these predictions to direct observations of edge-on protoplanetary disks. We compare the outer disk and jet orientation of four previously published edge-on disks observed at scattered light wavelengths. Our goal is to look for misalignment between the very inner regions of the disks, where the jet is launched, and the outer disk seen in scattered light or at millimeter wavelengths. Then, we compare apparent misalignments with the observed lateral asymmetry to assess the potential geometry of the systems.

We selected the disks based on the recent study of Villenave et al. (2024), who found that 15 out of the 20 systems they analyzed display strong lateral asymmetries. Seven sources of this sample display clear jet features, among which we identified four sources with a published value for the position angle (PA) of their jets from high angular resolution observations in the literature (with the other three disks, HV Tau C, Tau042021, and ESO Hα-574 having a jet whose PA is not explicitly characterized in the literature; Stapelfeldt et al. 2003; Duchêne et al. 2014, 2024). We report the values for the jet position angle in Table 2. We note, however, that the published jet PA of HH48, reported by Stapelfeldt et al. (2014), suffered from image distortion in the HST Advanced Camera instrument (K. Stapelfeldt, priv. comm.). In Table 2, we report our visual estimate of the jet PA of the distortion corrected image, leaving a more precise measurement for future work. All disks have been observed and well characterized at millimeter (with ALMA) or scattered light wavelengths (with HST), so we additionally report the outer disk position angle in Table 2. If no misalignment is present, we expect the position angle of the jet to be exactly perpendicular to that of the disk’s major axis. The scattered light images together with the corresponding jet orientation are displayed in Fig. 10.

All four systems show hints of small levels of misalignment of their inner disk. In the case of HH 48, we find a misalignment between jet and disk bottom nebula of about 2°, while no misalignment with the top nebula is identified. HH48 displays an apparent misalignment between its top and bottom nebulae, hinting at a distorted morphology. A somewhat similar misalignment can also be found in some of our models due to illumination effects, such as in Fig. 8, top two panels. However, for a good comparison, the current observational constraints on the dust and gas morphology of HH48 are not sufficient. In ESO Hα-569, the jet visible in the scattered light image does not have a published PA, but is thought to be launching the HH 919 jet, which could be consistent with either being aligned, or misaligned to the outer disk by up to 10°. Finally, even though the jet and disk of HH30 and Haro6-5B are consistent with no misalignment, the jet of both sources wiggles at large distances. This can be explained in case of a varying orientation of the accreted material, suggesting a precession of the inner disk (Terquem et al. 1999; of order of ~1° for HH30, Louvet et al. 2018). Alternatively, interaction with a non-uniform environment may lead jets to bend. However, regular wiggles as observed in these systems are not necessarily expected in that case.

Following the methodology previously described (Sect. 2.3, Villenave et al. 2024), we estimated the lateral asymmetry in these disks, as represented on Fig. 10 and indicated in Table 2. The strongest lateral asymmetry is found around Haro6-5B (10–30%), while ESO Hα-569 is barely asymmetric (3–5%). HH48 and HH30 show intermediate levels of asymmetries of order 10%.

HH48 shows the most convincing case for a distorted morphology and hints at a misalignment between the jet and the disk plane of up to 2°. Further, the disk presents a significant lateral asymmetry. Within the modeling presented in Sect. 3, we found that lateral asymmetry around 10% could be obtained for a warped disk with a misalignment between inner and outer disk of i_max = 1–2° if viewed with the optimal orientation (φ ~ 360°), or for stronger warps, up to i_max = 10°, for less favorable viewing angles (φ → 270°), which would be consistent with the apparent jet/disk misalignment observed. Both the significant lateral asymmetry and the locally aligned but wiggling jets of HH30 and Haro6-5B suggest that these systems could also possess warps, that are, however, weak or not viewed in a favorable orientation. Lastly, ESO Hα 569 could be consistent with either no misalignment or an unfavorable warp orientation.

Finally, we also mention that a local change in the disk scale height, leading to asymmetric shadowing of the outer disk, or variation in the stellar illumination (e.g., spots; Stapelfeldt et al. 1999) could also lead to lateral asymmetries. Dedicated imaging campaigns would be needed to determine the timescale of the variability of the lateral asymmetry and disentangle the mechanisms.

Fig. 10 Scattered light observations of four edge-on disks in log-scale, where we estimate the level of lateral asymmetry (blue spines) and indicate the estimated direction of the jet (white lines). A 0.5″ horizontal line is indicated in the bottom left corner of the images.

5 Conclusions

In this work, we present various near-/mid-infrared images from radiative transfer simulations of edge-on protoplanetary disks containing a warp. Edge-on disks in scattered light typically show two nebulae separated by a dark lane. These nebulae correspond to the two surfaces of the disk, whereas the dark lane corresponds to the optically thick midplane. Warped disks contain a misaligned inner region casting a shadow onto the outer disk, leading to visible asymmetries between the two nebulae. We quantified this asymmetry by fitting a spine curve for each nebula individually, and refer to the relative distance between the nebula centers as lateral asymmetry. We investigated the visibility of the lateral asymmetry for different warp parameters, especially the warp amplitude and the location of the steepest warp transition, or in other words, different sizes of the misaligned inner region. Additionally, we evaluated the flux ratio between the two nebulae. Because warped disks are not axisymmetric, we performed radiative transfer simulations for images viewed from different azimuthal angles. Most of our simulations are performed at a wavelength of λ = 0.8 μm, but we extended the simulations to λ = 2 μm, 12 μm, and 21 μm in specific cases.

We find that a warp can be visible in edge-on disks for warp tilts as low as i_max = 2° for optimal azimuthal viewing angles. The larger the warp tilt, the stronger the lateral asymmetry. Additionally, the location of the warp within the disk influences the asymmetry as well. If the warp is located further out, the asymmetry is stronger. This means that in edge-on disks, warps close to the inner edge of the disk are harder to detect. However, our models show a point where the asymmetry decreases for even further warp locations. We suspect this to occur due to the warp transition width that scales with radius. This leads to a shallower warp curve, and consequently more light can pass, causing smaller shadows and, therefore, less lateral asymmetry.

The flux ratio between the two nebulae can also indicate a warp in specific cases. In an unwarped and perfectly edge-on disk, the two nebulae have the same brightness. Due to the shadows in warped disks, the flux ratio can vary depending on the azimuthal angle. By itself, the flux ratio is not a strong indicator of a warp, since the relative brightness of the two nebulae also varies in an unwarped disk when it is not oriented perfectly edge-on. Nevertheless, we find that in some rare occasions, the flux ratio can be an indication of warps in observations after all, when observing the disk in multiple wavelengths. For this scenario, the outer edge of the warped disk should be slightly inclined (not perfectly edge-on; we chose an inclination of 85° in our models). If the disk is additionally oriented azimuthally, such that the nebula that would be brighter if unwarped is now shadowed by the misalignment, this nebula appears darker in the near-infrared. At slightly longer wavelength (12–21 μm in our models), however, the disk’s thermal emission is large enough so that this nebula can appear brighter, as the shadow is less prominent, leading to a swap in brightest nebula with wavelength. We confirm this swap in brightest nebula with wavelength if the warp is strong enough, in our models i_max = 10°. For small warps (i_max = 3°), we cannot reproduce this switch even at the best viewing angle.

In observations, such a swap is found in three edge-on disks so far: IRAS04302+2247, HH30, and Flying Saucer. As expected, this phenomenon is rare, as it requires a very specific orientation of the disk with respect to the observer. However, when it occurs, it is a strong indicator for a misalignment. In observations, the inclination of an observed disk is often inferred from the flux ratio. If the observed disk contains a misalignment, it could lead to a mismatch of inclinations inferred from different wavelengths. Such a mismatch is indeed found for example in Tazaki et al. (2025). We want to point out that observed disk inclinations usually have large uncertainties. However, if a swap with wavelength is observed, it can be a good indication that the disk is inclined and warped.

Additionally, we analyzed the relative position angle between optical jets and the outer disk plane in four previously published edge-on protoplanetary disks. We specifically targeted disks showing some level of lateral asymmetry in previous scattered light observations. In our investigation, we find that all systems show hints of slight misalignments. We use our models to make a rough estimation of the amplitude of misalignment in those disks. For HH48, we find that, depending on azimuthal orientation, a warp between 1–10° could be present, consistent with the measurement of jet and bottom disk surface misalignment (≲2°). In the other systems, the misalignments of a few degrees could be consistent with the apparent wiggling of the jet and the lateral asymmetries in the nebulae. In summary, the weak or lack of misalignments between the jets and outer disks could indicate that warps in protoplanetary disks are typically only of a few degrees. However, more observed systems are needed in order to draw strong conclusions.

Lateral asymmetries in edge-on observations can point to the presence of a warp, and with additional information about the vertical structure of the disk, it can be possible to make estimations about the warp parameters. As always, these results must be interpreted within their limitations, which we elaborated in Sect. 3.3. Although it is hard to infer the exact warp parameters from observations at this point, this study takes the first step in quantifying the asymmetry and trying to match the warp strength to the observed asymmetry. However, the parameter space of different warp shapes is too large to cover in one single work, and we therefore have to leave further explorations for future studies.

Acknowledgements

We thank Cornelis Dullemond, Giovanni Rosotti, Karl Stapelfeldt, Stefano Facchini, and Giuseppe Lodato for their valuable insight and support for this project. We sincerely appreciate the thoughtful and detailed feedback from the referee, which has enhanced the clarity and depth of our work. We also want to thank Philipp Weber and Thomas Rometsch for their numerical and contentual help. We also thank Ryo Tazaki, François Ménard, Gaspard Duchêne and the participants of the edge-on disk workshop (founded by the ERC grant agreement No. 101053020, project Dust2Planets, PI F. Ménard) which provided valuable feedback to this project. We remember the legacy of Willy Kley, who passed away in 2021. CK acknowledges funding from the DFG research group FOR 2634 “Planet Formation Witnesses and Probes: Transition Disks” under grant DU 414/23-2 and KL 650/29-1, 650/29-2, 650/30-1. MV and CK acknowledge support from the European Union (ERC Starting Grant DiscEvol, project number 101039651, PI: G. Rosotti). CK acknowledges support from Fondazione Cariplo, grant No. 2022-1217. Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. Additionally, we acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant INST 37/935-1 FUGG. Plots in this paper were made with the Python (Van Rossum & Drake 2009) library matplotlib (Hunter 2007). We acknowledge also the usage of numpy (Harris et al. 2020), astropy (Astropy Collaboration 2022), and scipy (Virtanen et al. 2020).

Appendix A On the forward scattering peak

For this work, we needed to perform a large amount of simulation runs, as each viewing angle φ in each wavelength requires its own radiative transfer simulation. To keep the computation time feasible for this study, each simulations should not take too long to run. One of the main factors for computation time in radiative transfer simulations is the amount of scattering photon packages used. The more photon packages are used, the closer the result gets to reality, but also the longer the simulation takes.

A commonly used technique to get away with fewer photon packages, but still end up with a smooth result, is called angular chopping (see, e.g., documentation of optool, Dominik et al. 2021). This technique makes use of the fact that light that is scattered within only a few degrees around the forward scattering direction, can be approximated as unscattered. This means that these photon packages, that would be scattered within this forward peaked regime, can be treated as if they did not interact with the scattering particle. Internally, optool uses the chopping angle (option -chop in optool) to set the values of the scattering matrix within the cone of the set angle to the value at the cone edge. The tool compensates for the chopping in order to ensure energy conservation by adjusting the scattering cross section, which leads to a slightly altered opacity for wavelengths with a significant forward peak in the phase function. Fig. A.1 shows the comparison of the scattering opacity with and without chopping.

Fig. A.1 Scattering opacity (mass-weighted) with chopping of 2° (solid) and without chopping (dash-dotted) for the DIANA standard model, as used in this work. The mass-weighted absorption opacity is indicated with the orange dotted line. The lower panel shows the albedo η = κscat/(κabs + κscat). The gray dotted lines indicate the wavelengths investigated in this work.

Clearly, chopping is especially important if the scattering phase function is extremely forward peaked, meaning that most photons are scattered within only a few degrees from the original direction. For short wavelengths, where the forward scattering peak in the phase function is extreme, Fig. A.1 shows that the chopping mechanism reduces the opacity in this regime, while it does not influence the longer wavelengths, where the phase function is not as strongly forward-peaked. Without the chopping mechanism, radiative transfer codes can have difficulty to deal with such an extreme forward peak, leading to bright spots in the resulting image if too few scattering photons are used. These bright spots can indeed be avoided by using an enormous amount of scattering photons. However, this would lead to long computation times.³

Fig. A.2 Phase function (i.e., Z11-component of the scattering matrices) for the DIANA standard model in the wavelength λ = 0.8μm without the chopping mechanism and with chopping by 2°, using varying internal resolutions (nang) in optool. The bottom panel shows the same curves as the top panel, but with the y-axis linear and stretched so that the difference for larger scattering angles becomes visible. Except for the cut by 2° at the forward peak, the phase function is nearly the same for the case without chopping and high internal resolution (green line) and with chopping for the default internal resolution (blue solid line).

By looking at the phase function for our dust setup, we see that at the wavelength of 0.8 μm, the phase function indeed is extremely forward peaked. In Fig. A.2 we plot the Z₁₁-component of the scattering matrix, whose integral determines the total scattering cross section σ_scat with (e.g., Dullemond et al. 2012) $${\sigma }_{\text{scat\hspace{0.17em}}}={\int }_0^{\pi }\sin \theta \mathrm{d}\theta {\int }_0^{2\pi }\mathrm{d}\varphi Z_{11}(\theta ,\varphi ).$$(A.1)

Since we consider three additional wavelengths in the main part of the paper, we additionally checked the scattering phase function for them. At 2 μm, we also find a strong forward peak, but it is less strong than for λ = 0.8 μm by several orders of magnitude. For the longer wavelengths of 12 μm and 21 μm, there is only a very weak forward-scattering peak, which means that a consideration at these wavelengths does not require the chopping mechanism. However, the chopping mechanism also does not have a large effect on the results at these wavelengths. For simplicity, we therefore treat all simulations in the main part using the same opacity setup with chopping of 2°.

Interestingly, we noticed that for such extreme forward scattering peaks as in our case, the default internal angular resolution⁴ of optool is not enough to accurately determine the scattering phase function. We initially noticed a discrepancy between the phase functions with and without chopping using the default angular resolution n_ang = 180 of optool, see yellow and blue lines in Fig. A.2. We found that this difference in the scattering phase function leads to a relative difference in intensity in the images of about 30%. We suspected the reason for this difference between the phase functions with and without chopping to be linked to the extreme forward scattering peak. The tool optool normalizes the scattering matrix to produce the correct integral. As the forward-scattering peak has a strong contribution to the integral, a slight inaccuracy within the peak could lead to strong inaccuracies for larger scattering angles. Therefore, a higher internal angular resolution ensures an accurate calculation of the forward scattering peak, leading to a more accurate determination of the normalized values for high scattering angles. On the other hand, if the extreme forward peak is cut off using the chopping mechanism, the normalization can be more accurate for high scattering angles in this case as well. Using a higher internal angular resolution for the scattering phase function gets rid of this difference, as shown in Fig. A.2 in the green line, which almost perfectly aligns with the blue line in the regime of higher scattering angles. Changing the internal angular resolution for the phase function with the chopping mechanism does indeed not significantly change the phase function (dotted line), as the normalization integral is sensitive to slight deviations when the peak is chopped. We therefore conclude that it is not necessary to use the higher angular resolution for simulations when using a chopped phase function.

Fig. A.3 Comparison of the RADMC-3D images between cases with chopping by 2° and without chopping for 3 · 1010 scattering photons (top row) and 108 scattering photons (bottom row). These models are calculated at the wavelength λ = 0.8 μm. The right column shows the relative error in percent. Here, we ignored intensity values of the images that are smaller than 10−17 in order to avoid artifacts.

To ensure that the chopping does not change the correctness of the result, we ran comparison simulations using different amounts of scattering photons in a test case. Because we performed this test in an early stage of this project, our test case has slightly different disk parameters. We used a warped disk with a disk tilt of i_max = 10° at a warp location of r_warp = 20 au. The range of the disk and the coordinate system are the same as for the simulations in the main part of this work, but the aspect ratio of the disk is h₀ = 0.05 at the reference radius R₀ = 5.2 au and the flaring index is i_fl = 0.25. Additionally, we used a higher dust mass than in our main models (i.e., M_dust = 10⁻⁴ M_*), which leads to a different shape and separation width of the top and bottom nebulae.

In total, we ran four simulations for this test: one set using n_{scat_phot} = 3 · 10¹⁰ scattering photons, and one set with n_{scat_phot} = 10⁸, each set both without chopping and with chopping of 2°. For the simulations without chopping, we used an internal angular resolution of the phase functions of n_ang = 720, while for the simulations with chopping, we kept the default resolution of n_ang = 180, since we were aiming at performing the final simulations using this default resolution.

Fig. A.4 Relative error between the model closest to reality (nphot_scat = 3 · 1010 without chopping, top left panel in Fig. A.3) and the computationally least expensive model (nphot_scat = 108 with chopping by 2°, bottom middle panel in Fig. A.3).

Fig. A.3 shows the image results for all four simulations, and in addition the relative error of the resulting intensities $\delta =\left|\frac{I_{\text{without chopping\hspace{0.17em}}}-I_{\text{with chopping\hspace{0.17em}}}}{I_{\text{without chopping\hspace{0.17em}}}}\right|$ for both sets. To avoid artifacts in the residuals due to very small values in the images, we ignored the areas with intensity values lower than 10⁻¹⁷. The bottom left panel of Fig. A.3 shows the bright spots if too few scattering photon packages are used in the case without the chopping mechanism. The bright spots especially occur in the transition regimes between the optically thick and optically thin part of the disk. They are unphysical artifacts due to the strong forward scattering peak at this wavelength. As expected, there are two ways to avoid these bright spots. One way is to use more scattering photons (top left panel), which is physically closer to reality but has the problem of long computation time. The other way is to use the chopping mechanism (middle column), where we achieve a similarly good result at a lower computational cost.

In Fig. A.4 we show the relative error between the model closest to reality and the trade-off model with low computational cost. The median of the relative error is 2.7%, and we therefore conclude that it is safe to use a chopping of the phase function of 2° in order to save computation time for our final models.

Appendix B Unwarped disk

As a reference case for the models of warped disks, we investigate a planar, untilted disk in near-/mid-infrared. As described in Sect. 2, our disk contains dust of grain size a = 10 μm. Because this model is perfectly axisymmetric, we do not need to consider different viewing angles.

Fig. B.1 Synthetic near-/mid-infrared radiative transfer images in log-scale in different wavelengths of an unwarped disk. The left column shows the disk observed perfectly edge-on, the right column viewed at 85°.

For an observer viewing the disk perfectly edge-on (inclination of 90°), top and bottom nebulae are perfectly symmetric across multiple wavelengths, as shown in Fig. B.1, top row. Here, nebulae have the same shape and brightness. Looking at the disk at a slight inclination of 85°, the top nebula appears brighter in all wavelengths (see the bottom row of Fig. B.1). This occurs mainly due to forward scattering, but backward scattered light from the far side of the disk contributes as well.

Appendix C Testing smaller grain sizes

In this section we explore the lateral asymmetry and flux ratio in models with smaller dust grain sizes in the example of a misalignment of i_max = 3° and a warp location of r_warp = 20 au. In the main part, we focus on a single dust grain size of a = 10 μm. Here, we extend our investigations to a = 1 μm and 0.5 μm, as well as a dust size distribution for comparison. For the dust size distribution, we chose dust sizes from a = 0.1 μm − 10 μm, distributed in a commonly assumed power-law distribution according to n(a)da ∝ a^−3.5da. All other parameters remain the same as in the main part of this work.

Fig. C.1 Radiative transfer simulations of the model with a misalignment of 3° at a location of 20 au, viewed from φ = 0°, for different dust sizes. Images are shown for λ = 0.8 μm.

We find that the appearance of the disk nebulae in the radiative transfer images differs between the models of different dust sizes, as shown in Fig. C.1. In especially, the shape of the nebulae is less rounded and the separation between the nebulae is wider for the smaller dust grain sizes. This is because of the different opacities due to the different models of dust grain sizes. However, when comparing the lateral asymmetry and flux ratio, we find the same trend as in our main models, as shown in Fig. C.2. This good agreement between the different assumptions on the dust indicates that the results of this work in terms of lateral asymmetry and flux ratio can be applied in a more general way.

We further checked for the different dust size models, if the swap in flux ratio with wavelength (see Sect. 3.2) can be reproduced for smaller grains. Fig. C.3 shows that all three dust size models can reproduce the swap in the brightest nebula with wavelength. In particular, for all alternative dust models explored in this section, the flux ratio swap occurs at shorter wavelength than in the case of our fiducial model with a = 10μm, which is noticeable at 12μm instead of 21 μm. We show in Fig. C.4 the images at the different wavelengths for a dust size of a = 1 μm.

We note that for our chosen fiducial dust size of a = 10 μm, the silicate feature in the absorption opacity at roughly λ = 10 μm is weak. However, the other models shown here have a stronger silicate feature because of the smaller grains. Checking the appearance of the images for wavelengths around λ = 9–12 μm, we found that the silicate feature does not have a strong effect on the lateral asymmetry.

Fig. C.2 Same as Fig. 6 but for a misalignment of imax = 3° and different assumptions on the dust size in the models. All data points are for radiative transfer models at the wavelength λ = 0.8 μm.

Fig. C.3 Flux ratio as a function of wavelengths for models of warped disks with a warp tilt of 10° and different dust properties. The models are simulated with φ = 60°, and the y-axis is plotted in log-scale.

Fig. C.4 Same as Fig. 8 but for the models with a dust grain size of a = 1 μm instead of 10 μm.

References

All Tables

All Figures

	Fig. 1 Vertical cross section of the gas density setup for an unwarped disk, a warp of 3°, and a warp of 10°. The dash-dotted line indicates the disk midplane.
In the text

Fig. 2 Definition of the azimuthal φ-angles. The colors represent the height above the x–y-plane: purple is in plane, blue-green above, and red-yellow below. If the outer disk is viewed perfectly edge-on, the quadrant of 270°–0° is identical to the quadrant 0°–90°, only flipped upside down. The sector 90°–270° is identical to the points of sector 270°–90° appropriately rotated and flipped. We chose this definition of φ to be consistent with Villenave et al. (2024).

In the text

Fig. 3 Schematic representation of the spine fitting. The light circles show the location of the maxima after a high-order polynomial fit of the minor-axis profiles. The solid lines represent the final spines after fitting the maxima by a polynomial of order 2. The diamond represent the locations where the flux reaches 10% with respect to the maximum of the spine, and the squares are the center of the spine (i.e., the mid-point of these two points). The lateral asymmetry is the ratio between the distance of the centers to averaged disk size.

In the text

	Fig. 4 Near-infrared radiative transfer models of disks with different amplitude warps, defined by imax, observed at λ = 0.8 μm at φ = 0° (top two rows). The bottom row shows models with a varying warp location or in other words a varying size of the misaligned inner region. For the color, we use a logarithmic scale.
In the text

Fig. 5 Near-infrared radiative transfer models (λ = 0.8 μm) of a disk with tilt imax = 3° observed from different azimuthal angles φ (color in log-scale). The top row shows a perfectly edge-on view, while the bottom row shows the disk observed at an inclination of 85°. We note that for symmetry reasons, we do not need to simulate φ = 90°. The image shown in this figure is the simulation of φ = 270° flipped upside down.

In the text

Fig. 6 Lateral asymmetry (top) and flux ratios (bottom) for warped disk models of different warp tilts imax (at λ = 0.8 μm). The outer disk is viewed perfectly edge-on (meaning an inclination of 90°). For symmetry reasons, we only need to model a quarter of the azimuthal regime in this configuration. The curves can be extended to the remaining regime by flipping the curves appropriately. The pink shaded region in the top panel corresponds to a lateral asymmetry of <10%, which is what we assume not to be visible by eye.

In the text

	Fig. 7 Same as Fig. 6, but for different warp locations rwarp for a warp amplitude characterized by imax = 3°. Again due to symmetry reasons, the curves can be extended to the whole azimuthal regime. As in Fig. 6, the pink shaded region highlights the lateral asymmetry we assume not to be visible by eye.
In the text

	Fig. 8 Near-/mid-infrared radiative transfer models of a disk with tilt imax = 10° at φ = 60° with an inclination of 85° observed in different wavelengths (color in log-scale). Considering the integrated flux, the bottom nebula is brighter than the top at shorter wavelengths, but at λ = 21 μm, the top nebula is brighter.
In the text

Fig. 9 Flux ratios for models of warped disks with a warp tilt of 10° (solid lines) and of 3° (light dashed lines) in four different wavelengths. The outer disk is oriented with an inclination of 85°. In this case, we need to calculate the models for half of the φ-regime, the other half is symmetric. The dashed black line indicates the line where the flux ratio equals one (where both nebulae have the same brightness). We note that the y-axis is plotted in log-scale.

In the text

	Fig. 10 Scattered light observations of four edge-on disks in log-scale, where we estimate the level of lateral asymmetry (blue spines) and indicate the estimated direction of the jet (white lines). A 0.5″ horizontal line is indicated in the bottom left corner of the images.
In the text

	Fig. A.1 Scattering opacity (mass-weighted) with chopping of 2° (solid) and without chopping (dash-dotted) for the DIANA standard model, as used in this work. The mass-weighted absorption opacity is indicated with the orange dotted line. The lower panel shows the albedo η = κscat/(κabs + κscat). The gray dotted lines indicate the wavelengths investigated in this work.
In the text

Fig. A.2 Phase function (i.e., Z11-component of the scattering matrices) for the DIANA standard model in the wavelength λ = 0.8μm without the chopping mechanism and with chopping by 2°, using varying internal resolutions (nang) in optool. The bottom panel shows the same curves as the top panel, but with the y-axis linear and stretched so that the difference for larger scattering angles becomes visible. Except for the cut by 2° at the forward peak, the phase function is nearly the same for the case without chopping and high internal resolution (green line) and with chopping for the default internal resolution (blue solid line).

In the text

	Fig. A.3 Comparison of the RADMC-3D images between cases with chopping by 2° and without chopping for 3 · 1010 scattering photons (top row) and 108 scattering photons (bottom row). These models are calculated at the wavelength λ = 0.8 μm. The right column shows the relative error in percent. Here, we ignored intensity values of the images that are smaller than 10−17 in order to avoid artifacts.
In the text

	Fig. A.4 Relative error between the model closest to reality (nphot_scat = 3 · 1010 without chopping, top left panel in Fig. A.3) and the computationally least expensive model (nphot_scat = 108 with chopping by 2°, bottom middle panel in Fig. A.3).
In the text

	Fig. B.1 Synthetic near-/mid-infrared radiative transfer images in log-scale in different wavelengths of an unwarped disk. The left column shows the disk observed perfectly edge-on, the right column viewed at 85°.
In the text

	Fig. C.1 Radiative transfer simulations of the model with a misalignment of 3° at a location of 20 au, viewed from φ = 0°, for different dust sizes. Images are shown for λ = 0.8 μm.
In the text

	Fig. C.2 Same as Fig. 6 but for a misalignment of imax = 3° and different assumptions on the dust size in the models. All data points are for radiative transfer models at the wavelength λ = 0.8 μm.
In the text

	Fig. C.3 Flux ratio as a function of wavelengths for models of warped disks with a warp tilt of 10° and different dust properties. The models are simulated with φ = 60°, and the y-axis is plotted in log-scale.
In the text

	Fig. C.4 Same as Fig. 8 but for the models with a dust grain size of a = 1 μm instead of 10 μm.
In the text