© The Authors 2026

1 Introduction

In recent years, significant revelations have been reported while trying to answer the question of how exoplanets form (for a review, see Drążkowska et al. 2023) from theoretical and observational perspectives. One of the great challenges in the field is understanding the coagulation of initially small dust grains in protoplanetary disks, the birthplaces of planets. As the grains grow, even though many challenges still remain unexplained, they form the necessary building blocks for planets. Observed substructures such as rings, gaps, and spirals in these disks can significantly affect and shape the growth process (for a review, see Bae et al. 2023). Within the past few decades, many questions have arisen that are yet to be fully answered, concerning, among others, the available dust budget, its distribution through the disk, the timing of the emergence of the substructures, and whether they cause planet formation or are a consequence thereof.

Planet formation theory has classically treated the Class II evolutionary stage of a disk, when most of the material from the natal cloud has disappeared, as the starting point. Consequently, planet formation has often been modeled as isolated from the disk environment. However, recent observations have challenged this picture (e.g., Ginski et al. 2021; Huang et al. 2021; Garufi et al. 2022; Gupta et al. 2024; Speedie et al. 2024), resulting in a large paradigm shift, where models specifically target the environmental interactions and their effect on the disk. Especially, the ongoing interaction of early-stage disks with their natal envelope has been the focus of recent work (e.g., Drążkowska & Dullemond 2018; Morbidelli et al. 2022; Hühn et al. 2025a; Carrera et al. 2025), with the aim to solve the issue of insufficient solid mass being available for forming known giant planets during Class II (Manara et al. 2018, though it is debated, Savvidou & Bitsch 2025) by invoking that the formation of planetesimals, the precursors of planets, might start early, when more solids are available.

Observations indicate a ubiquitous number of disks that also interact actively with their environment at later stages in their evolution, which is evident by large-scale structures, such as arms and tails (commonly refereed to as streamers), which are suggested to result from infall. About 30% of the observed disks in the Taurus star-forming region exhibit these types of structures (Garufi et al. 2024). Thus, it is vital to investigate the nature of these large-scale structures and their effect on the formation of planets. Similar to early-stage infall, which can provide a large number of solids (Hühn et al. 2025a), late accretion might deliver new, potentially chemically distinct material to the disk, affect its kinematics, and introduce substructure that might aid or harm the planet formation process. It might even form an altogether new disk (Kuffmeier et al. 2020; Hühn et al. 2025b).

Streamers can be interpreted as material that is accreted sporadically from the environment (also called late infall). Simulations with simplified models of the capture of a spherical gas cloudlet by a star and disk system (Dullemond et al. 2019) have facilitated the characterization of the structures seen in AB Aur (Calcino et al. 2025b) and DG Tau (Hanawa et al. 2024) as streamers. However, these cloudlet encounters might be rare, and these models may be an oversimplification; in reality, it might be necessary to assume that late infall is a result of Bondi-Hoyle-Lyttleton (BHL) accretion in the turbulent interstellar medium (ISM), for instance, from molecular remnants of the parent giant molecular cloud. Models of this type of interaction are able to produce direct (Hühn & Dullemond 2025) and indirect (Winter et al. 2024) infall characteristics of disks.

Next to the direct evidence provided by streamers, an indirect indicator of interaction of the disk with the environment are warps. With the increasing observational indications of warps and misaligned inner disks (Ansdell et al. 2020; Bohn et al. 2022; Villenave et al. 2024; Kimmig & Villenave 2025; Winter et al. 2025), it is important to understand their origin. Late infall is one of the prime candidates for the creation of these warps (Kuffmeier et al. 2021; Dullemond et al. 2022), further indicating that infall might also be common at later disk stages. Other scenarios for forming a warp include gravitational torques from stellar flybys (e.g., Cuello et al. 2023; Kimmig et al. 2026), inclined binaries (Facchini et al. 2013), or inclined planets (Nealon 2018; Zhu 2019). A warped disk can naturally exhibit additional substructure such as spirals when it is also twisted (Winter et al. 2025), or shadows (Marino et al. 2015; Facchini et al. 2018), which in turn can lead to the formation of spirals (e.g., Su & Bai 2024; Zhang & Zhu 2024; Zhu et al. 2025; Ziampras et al. 2025).

Several systems with streamers and ongoing infall exhibit spiral structures. The spiral pattern can either consist of many flocculent arms, such as in SU Aur (Ginski et al. 2021), or of a few well-defined ones, such as in DR Tau (Mesa et al. 2022; Garufi et al. 2024). Previous work investigating the formation of these structures revealed that infall might be the direct cause of m = 1 spirals (Hennebelle et al. 2017; Calcino et al. 2025a) and might explain at least parts of the flocculent spirals in AB Aur (Calcino et al. 2025b). Gravitational instability (GI) is also theorized to offer an explanation (Speedie et al. 2024). However, GI might be short-lived when irradiation from the central star is taken into account (Rowther et al. 2024), so that it is natural to assume that both effects operate at the same time, with the infall supplying fresh material to keep GI active (Longarini et al. 2025). Well-defined spiral arms have also been shown to be the result of a stellar flyby (Quillen et al. 2005), be caused by planets (Baruteau et al. 2014; Muley et al. 2024), or by gap-edge illumination (Muley et al. 2026).

This raises the question of what distinguishes the signatures caused by different mechanisms and how they can be distinguished. Furthermore, spiral patterns found in the scattered light with the Spectro-Polarimetric High-contrast Exoplanet Research (SPHERE) instrument of the Very Large Telescope (VLT) do not always match the spirals found in CO line observations with the Atacama Large Millimeter/submillimeter Array (ALMA). For example, MWC 758 shows a one-armed spiral in ¹²CO line emission (Winter et al. 2025), while scattered light reveals two spiral arms (Benisty et al. 2015; Ren et al. 2023). In the millimeter continuum, the most prominent features are vortices, although a weak one-armed spiral continues to be visible (Dong et al. 2018). For other disks, there might even be no spiral substructure in the continuum emission, even though other types of observations indicate them, for example, in SU Aur (Ginski et al. 2021). Moreover, there can be discrepancies between the isotopologues of CO (e.g., Huang et al. 2025), as these observations trace different disk layers. Infall might also cause the formation of spiral structures even without direct evidence. For example, the MWC 758 system has no confirmed streamers. To explain its spiral structures, two planet candidates were proposed (Wagner et al. 2019, 2023) instead, but it is unclear whether these planets are connected to the spirals (Boccaletti et al. 2021). A flyby is also an unlikely explanation in the case of MWC 758 (Grady et al. 2013; Reggiani et al. 2018), suggesting the possibility that infall might contribute to the formation of the spirals despite the absence of streamers.

We investigate the spiral structures that arise in select simulation setups of late infall by HD25, under which conditions they occur, how long they are present, and how they affect different layers of the disk. This work acts not only as a comparison to previous similar works using smoothed-particle hydrodynamics (SPH) simulations (Calcino et al. 2025a), but also as a more detailed characterization of the spiral structures that can be caused by late infall in the context of the vertical dependence of disk properties and substructures indicated by the aforementioned observations. A comparison to observed systems also enables insights into the role played by the theorized environmental interactions in forming them and which conditions provide a good match. A broader understanding of the location and the way in which disk substructures are shaped provides more clues on its impact on planet formation, which predominantly occurs close to the midplane, where large grains settle (Youdin & Lithwick 2007), when the streaming instability (Youdin & Goodman 2005) is the mechanism invoked for planetesimal formation (for a review, see Lesur et al. 2023).

2 Methods

We performed 3D hydrodynamical simulations on a spherical grid with logarithmic spacing in the radial direction, and uniform spacing in the polar and azimuthal direction. The computations were performed using the FARGO3D code (Benítez-Llambay & Masset 2016), capable of GPU multiprocessing while using the advection algorithm by (Masset 2000). We followed our approach from Hühn & Dullemond (2025, HD25), and we refer to that paper for a detailed description. In the following, we briefly summarize the main aspects of the simulations.

2.1 Types of accretion models

In HD25, we simulated the late accretion of gas onto a Class II disk, which is included as part of the initial condition, in two different scenarios. The first model is the capture of a cloudlet of gas, initialized in compact form on a hyperbolic orbit. We did not confine the cloudlet in any way. As a result, the cloudlet expands considerably before it encounters the disk, engulfing it in gas rather than creating a gas trail directly onto it. The goal of this more simplistic cloudlet capture setup is to model the encounter of the disk with an interstellar medium (ISM) overdensity in a controlled manner, and not to model the physical conditions in realistic detail. Initializing the cloudlet with a compact size has the advantage of keeping the inaccuracy incurred by treating it as a point source when calculating the initial orbital velocity low.

In the second model, we did not include any initial overdensities. Instead, we model a compressively turbulent background ISM with a systemic velocity with respect to the disk. This scenario represents a more realistic case of turbulent Bondi-Hoyle-Lyttleton (BHL) accretion onto the disk, where overdensities and related streamers emerge naturally as a result of the turbulence. As the disk moves through the ISM, an accretion tail emerges, and turbulent asymmetries result in net angular momentum relative to the star, creating streamers (see HD25 for details). This scenario is more realistic, but offers less control over the accretion process, which is mostly regulated by the infall rate, velocity dispersion, and turbulent power spectrum. We did not drive the turbulence; it was only included as an initial condition.

2.2 Numerical setup

We set the inner edge of the grid to r_min = 5 AU. For the cloudlet capture simulations, we employed r_max = 5000 AU, whereas we set r_max = 50 000 AU for all simulations with a turbulent ISM as initial condition. In the polar direction, the grid extends from θ_min = 10^◦ to θ_max = 170^◦. In the radial and polar direction, we used an outflow boundary condition; the azimuthal boundary was periodic. For the simulations of cloudlet capture, we chose N_r = 366 cells in radial, N_θ = 146 cells in polar, and N_ϕ = 330 cells in azimuthal direction, whereas we used N_r = 486 cells for the turbulent medium simulations. This results in a number of cells per scale height of H/∆r = H/(r∆θ) = H/(r∆ϕ) ≈ 2 at r = 5.2 AU (≈4.2 at 100 AU). Here, H is the gas pressure scale height and ∆r, ∆θ and ∆ϕ describe the resolution in the respective direction.

For the initial condition, we included a protoplanetary disk with a symmetric surface density profile of Σ(r) = Σ₀ (r/R₀)^−a, with R₀ = 5.2 AU and a = 1.5. We applied an exponential cutoff at 100 AU. For each accretion model, we considered two scenarios. First, a high disk mass, Σ₀ = 310.547 g cm⁻² resulting in M_d = 0.05 M_⊙, and second, a low disk mass, Σ₀ = 31.0547 g cm⁻², resulting in M_d = 5 × 10⁻³ M_⊙. The simulations were performed using an isothermal equation of state and a solar-mass star, with an aspect ratio of h(r) = h₀(r/R₀)^f with h₀ = 0.03799 and f = 0.25. We included an α-viscosity of 10⁻³ (Shakura & Sunyaev 1973), but we note that the numerical diffusion introduced by the grid also contributes to this value, and that it could be significant at our employed resolution. The floor density of the simulation was 10⁻²² g cm⁻³.

2.3 Initial conditions for the accretion model

In the cloudlet scenario, we initialized the cloudlet with v_∞ = 0.5 km s⁻¹, d₀ = 3000 AU, R_cloud = 500 AU, and M_cloud = 5 × 10⁻³ M_⊙. This results in an initial density of the cloudlet of ρ_cloud = 5.67 × 10⁻¹⁸ g cm⁻³. The cloud initially expands due to the lack of pressure support, leading to a mean density of ∼10⁻¹⁹ g cm⁻³ at the time of the encounter with the disk, which is comparable to filamentary structures in the Taurus star-forming region (Pineda et al. 2010; Palmeirim et al. 2013). The total simulated time of the cloudlet simulations is 42.7 kyr. We set the impact parameter to $b=0.5b_{\text{crit}}=0.5G{M}_{\odot }/v_{\infty }^2$. Different from the case considered in Hühn & Dullemond (2025, HD25), the orbit of the cloudlet is not inclined with respect to the plane of the disk, leading to in-plane capture. We chose this configuration because we suspect it to be most favorable for purely infall-driven spiral formation. First, the interaction between infalling material and gas is concentrated at the disk edges, so that the perturbation on the disk surface has the highest degree of asymmetry, promoting spiral formation. Second, the contamination of emission from the disk by the arising streamers and the accretion tail is minimal for this orientation. Third, the effect on the disk dynamics, for example, potential warping or tilting, which can also cause spiral patterns (Winter et al. 2025), is also minimal for this configuration, allowing us to separate the impact of the different formation pathways.

For the simulations with a turbulent medium, we considered a systemic infall rate of Ṁ_sys = 10⁻⁸ M_⊙ yr⁻¹ and a systemic velocity of v_sys = 0.5 km s⁻¹, resulting in ρ_ISM = 3 × 10⁻²¹ g cm⁻³. Unlike HD25, we chose the direction of the systemic velocity such that it leads to in-plane streamers, that is, the systemic velocity vector and angular momentum vector of the disk are perpendicular, û_sys = (−1, 0, 0), for the same reason as for the cloudlet capture case. We included small and large scales in the turbulence power spectrum, with k_min = 2π × 10⁻⁵ AU⁻¹ and k_max = 2π/50 AU⁻¹. The turbulent velocity dispersion is σ_turb = 0.5 km s⁻¹. The simulated time of the turbulent medium simulations is 30 kyr.

2.4 Synthetic observations

Using the Monte Carlo radiative transfer code RADMC3D (Dullemond et al. 2012), we created synthetic observations of select simulation snapshots, as in HD25. We produced observations of the ¹²CO and ¹³CO J 2–1 transition line, assuming a constant ¹²CO abundance of 10⁻⁴ and a ¹³CO abundance of 10⁻⁴/70 (Wilson & Rood 1994), using 250 points in frequency space covering a spectral range from −5 km s⁻¹ to 5 km s⁻¹ to create first moment maps. Spatially, we used 400 pixels for each axis, covering a range from −100 AU to 100 AU. As in HD25, we did not employ any prescriptions for the freeze-out or photo-dissociation of the CO isotopologues and recomputed the temperature using RADMC3D instead of using the isothermal temperature distribution from the simulations. We assumed T_gas = T_dust.

The main focus of this work is on synthetic scattered-light observations. Here, the required dust opacities were calculated using the optool code (Dominik et al. 2021), where we assumed the standard dust model composition from the European FP7 project DIANA (Woitke et al. 2016). We assumed a porosity of 25% and mixed 15 different grain sizes between 1 μm and 3 μm, sampling an n(a) ∝ a^−2.5 power-law size distribution. We set the dust-to-gas volume density ratio to 1%, which holds true under the assumption that the grain size distribution on the disk surface, perturbed by the infall, is predominantly given by small grains, because large grains are removed from the surface due to settling, and new small grains are constantly replenished from the ISM. This assumption would not hold for the disk midplane, where the large grains are present and contain a majority of the mass, but these disk regions are not probed by the synthetic scattered light observations. For better numerical treatment of the scattering peak, we applied a chopping of 2^◦ (Dominik et al. 2021), which is shown to be a reasonable assumption in similar conditions of protoplanetary disks (see Kimmig & Villenave 2025, Appendix A). To enable the study of the scattered light emission as a function of time, analyzing many simulation snapshots, we utilized that simple scattering approximation included in RADMC3D, which enables direct integration of the radiative transfer problem by taking advantage of the specific problem geometry. It assumes that all scattered light is starlight, that thermal emission can be neglected, and that multiple scattering does not occur. For all direct depictions of scattered light emission shown in subsequent sections, we compared this approximation to a full Monte-Carlo model and find no qualitative difference.

2.5 Determining residual motions

We characterized the spiral structures found using the synthetic CO line emissions by considering the residual motions, that is, the deviation of the determined moment 1 velocities from a Keplerian disk model. While we could have, in principle, full knowledge of the emission from the unperturbed disk from the initial condition, key parameters that are relevant for an observation determination of an accurate Keplerian model are altered due to the accretion that occurs over the course of the simulation. For example, the apparent emission height, inclination, and position angle are affected. To maintain comparability to observations, we therefore discarded any information obtained from the knowledge of the initial conditions, and fit a Keplerian model purely based on the synthetic observation at the considered time.

The fit was performed using the Monte Carlo disk fitting code discminer (Izquierdo et al. 2021). In order to isolate the disk emission from cloud contamination as much as possible, we restricted the disk model to a maximal size of 100 AU, and excluded emission from the inner 10 AU. In the interest of computational cost, we downsampled the synthetic image by six pixels, and we only considered emission from within 150 AU. Therefore, only emission at a distance 10 AU < r < 150 AU was considered. We performed the fit using 250 walkers that take 10000 steps while fixing the stellar mass to 1 M_⊙, enforcing a centered disk, and fixing the systemic velocity at zero.

Fig. 1 Polarized scattered-light intensity in four snapshots of simulation 1. The solid gray line denotes the 1 mJy arcsec−2 contour, the dashed gray line shows the 0.1 mJy arcsec−2 contour, and the dotted gray line shows the 0.01 mJy arcsec−2 contour. The panel numbering corresponds to the phases of the cloudlet encounter. The direction of the white arrow denotes the original orbital velocity of the cloudlet.

3 Results

We investigated the formation of spiral structures, mainly those on the disk surface which are probed by scattered light observations, in simulations with setups similar to those in HD25. In particular, we considered a cloudlet capturing onto a disk with a mass of 5 × 10⁻² M_⊙ in-plane (simulation 1) and BHL accretion with an accretion rate of 10⁻⁸ M_⊙ yr⁻¹ onto a disk of the same mass (simulation 2). We also consider these scenarios for a disk with a mass of 5 × 10⁻³ M_⊙ (simulations 3 and 4). We put the main emphasis on simulation 1, and explored how a different infall model, or infall onto a lighter disk, changes the picture in subsequent sections.

3.1 Cloudlet capture

The capture of a cloudlet of gas leads to various spiral structures at different times in the simulation. These structures, either spirals that have a low number of arms or that are very flocculent, can be characterized depending on the environmental conditions at that time, and whether they occur in the inner or in the outer disk regions. We show synthetic polarized scattered light observations for four representative snapshots in Fig. 1. The camera angle was chosen such that the disk appears fully face-on. We identify four different phases in our simulation. The first phase is characterized by the initial, main encounter of the cloudlet with the disk, before overshooting gas falls back and streamers form. Here, we find little structure in the outer disk, whereas the inner disk shows several flocculent spiral arms. These flocculent spiral arms persist once the fallback streamer arises, which we call the second phase, but they are considerably different in shape and have fainter emission. At this stage, we also find a strong two-armed spiral, which we refer to as an m = 2 mode, emerging in the outer disk for the first time in the simulation.

We find that another m = 2 spiral structure arises after the fallback streamer has largely dissipated, which we classify as the third phase. This structure is more clearly visible than the one during phase 2. It stretches from inner to the outer disk and is not accompanied by flocculent structures. At this time, the bulk of the mass has already been accreted; this structure is therefore solely caused by the accretion of remnant gas leftover from the cloudlet encounter. Here, the spiral structure is also directly reflected in the contour lines of the emission (see Fig. 1, bottom left panel). Qualitatively, the spiral persists for a long time (∼20 kyr), but it stops extending to the outer disk after ∼10 kyr, leaving behind a structure in the inner disk that appears to be more wound-up (see phase 4 in Fig. 1).

To quantify the strength of the individual spiral modes, we performed an analysis of the Fourier-transformed azimuthal intensity distribution, shown in Fig. 2, where $m=\tilde{\varphi }$ is the Fourier transform of the azimuthal coordinate. The transformation was performed at two distinct radii that we are representative of structure in the inner (r = 25 AU) and outer (r = 70 AU) disk. The classification into the four phases is based on the excitation of the Fourier modes, shown in Fig. 2. This figure confirms the previous findings, where the amplitude of the m = 2 mode is particularly strong in the inner disk starting at phase 3, and the first m = 2 mode with strong amplitude is found in the outer disk during phase 2. The presence of flocculent spirals can be seen in phases 1 and 2, where also higher-order modes are excited. The Fourier analysis suggests that the m = 1 mode is also excited in the outer disk during phase 1. The excitation of this mode is not related to the presence of an m = 1 spiral, but caused by the contamination of the scattered light image from the accretion tail caused by the cloudlet interacting with the disk.

Fig. 2 Azimuthal Fourier decomposition of the polarized scattered-light intensity, normed by the total intensity. $m=\tilde{\varphi }$ is the Fourier transform of the azimuthal coordinate, and the two panels show the transformation for two radii. The dashed white lines separate the phases of the encounter, labeled by the white numbers.

Fig. 3 Polarized scattered-light intensity as a function of azimuth and time for two different radii. The dashed white lines and numbers denote the different encounter phases, as in Fig. 2. The colored dots show the azimuthal peaks in gray when they were not considered further, otherwise by membership to one of the spiral arms. The solid white lines show the linear fit used to determine the pattern speed.

3.1.1 Pattern speed

In an effort to characterize the m = 2 spirals found across the majority of the disk in phases 3 and 4, we calculated the pattern speed, an observable spiral property that is particularly useful when the spirals are caused by a perturber, for example, a planet or another star. We utilized the peak finding algorithm from the scipy Python package to find the location of the spirals in azimuth at the same radii as for the modal analysis above, that is, r ∈ {25, 70} AU. In other words, we found all peaks in two rings at each point in time. The algorithm was applied to the logarithmic intensity values, which were divided by the median initial intensity at each radius. We required that a peak is brighter than the initial median intensity and has a topographic prominence¹ of 0.3. The results are marked in Fig. 3. Afterwards, we used the clustering algorithm DBSCAN from the Python package sklean to separate the two spiral arms. We only considerd peaks that we determined visually to be part of the m = 2 spiral (marked in green and blue in Fig. 3), and discarded all others (marked in gray in Fig. 3).

The pattern speed was then determined through a linear fit of the azimuthal peaks over time, as shown in Fig. 3. We find that the spirals are close to stationary. For the inner disk, the arm marked in blue has a pattern speed of $\dot{\varphi }=0.08{\text{kyr}}^{-1}$, corresponding to a corotation radius r_corot = 1840 AU and the green one has $\dot{\varphi }=0.11{\text{kyr}}^{-1}$ (r_corot = 1460 AU). The pattern speed is especially slow in the outer disk. Here, the blue spiral arm moves with $\dot{\varphi }=0.05{\text{kyr}}^{-1}$ (r_corot = 3150 AU) and the green one with $\dot{\varphi }=0.07{\text{kyr}}^{-1}$ (r_corot = 2040 AU). This result acts as a defining difference to spirals created by other mechanisms, where either a Keplerian pattern speed or one corresponding to a specific perturber are expected, so that the inferred r_corot would be smaller.

Measurements determining spiral pattern speeds have been conducted observationally. For the spiral in the disk around SAO 206462, for example, a pattern speed corresponding to r_corot = 86 AU was measured (Xie et al. 2021). Additionally, in V1247 Ori, the pattern speed corresponds to r_corot = 118 AU. These observed values are much lower than what we find in our simulations. Thus, this can serve as a way to distinguish spirals formed purely through infall from those formed through other mechanisms. We confirmed the low pattern speed of the spirals in our simulations not to be an effect of the output frequency.

Fig. 4 Moment 1 residuals from Keplerian motion of the 12CO line emission. The top row shows three snapshots, taken at the same time as in Fig. 1, where the camera is oriented face-on. The bottom row instead shows images where the camera has an inclination of 30◦. The black lines show polarized scattered light contours, as in Fig. 1. The panel numbers correspond to the evolutionary phase.

3.1.2 Gas kinematics

Substructures found in polarized scattered light trace the density structure of the disk surface layer. In order to obtain information about the kinematical structure and deeper disk layers, we investigated the gas kinematics by creating synthetic observations of ¹²CO and ¹³CO J 2–1 line emission, and considering the first moment of that emission. Due to the fact that the disk is surrounded by gas from the environment, and that we neglect possible photodissociation of CO, the emission of ¹²CO is very optically thick, and the τ = 2/3 surface, where τ is the optical depth, is high above the midplane (see Fig. 6 below). For the ¹³CO emission on the other hand, the τ = 2/3 surface lies deeper inside the disk, where perturbations created by the infall are smaller. We investigated deviations from a purely Keplerian rotation, which would correspond to a fully unperturbed, isolated disk, by computing the residuals of the first moment map using the discminer code. The usage of its Monte-Carlo fitting routines allows for the determination of the apparent disk inclination and position, which is different from the expected orientation given by the camera angle, as possible asymmetries and changes in disk tilt can be introduced during the accretion. The determined residuals are shown in Fig. 4 for ¹²CO emission, and Fig. 5 for ¹³CO emission. The discminer fit parameters can be found in Tables E.1, E.2, and E.3.

The optically thick ¹²CO emission, tracing the upper layers of the disk, shows spiral structures in its residuals. They are most apparent during phases 3 and 4, at the time when the m = 2 spiral is prominent in scattered light. We find that the physical location of the spirals in the disk does not coincide with their location in scattered light. However, the spiral shape matches between the ¹²CO and scattered light emission. The spirals are most prominent when the disk is viewed face-on, suggesting that they are caused by vertical motion induced by the cloudlet capture. In contrast, the residuals we find when the disk is viewed at a 30^◦ angle reveal a more large-scale structure in the outer disk, which we interpret to be a one-armed spiral for phase 2, and a two-armed one for phases 3 and 4.

We find that the emergence of substructures in the CO emission depends strongly on the viewed isotopologue. For ¹³CO, residuals of considerable amplitude are only found during phase 2, where an active streamer is falling onto the disk. The distribution of substructures matches that of the ¹²CO emission, but at a lower amplitude. For the inclined view, a spiral is visible in phase 2, but cannot be found in the later phases. We note that the residuals shown in these later phases are the result of a model limitation in the discminer code, which is related to the characterization of the emission surface and errors introduced by including environment emission in the fit of the Keplerian model.

Fig. 5 Same as Fig. 5, but for 13CO. We note that the value range of the color bar is reduced.

3.1.3 Perturbation depth

The amplitude of the visible structures in the residual motions strongly depends on the isotopologue. This indicates a weak overall influence of the infall on disk kinematics, meaning that strong perturbations only occur in the upper disk layers, where the density is low. This result is expected, as the total mass of the cloudlet is only 10% of the disk mass. In Fig. 6, we characterize the perturbation depth in more detail by calculating the azimuthal coefficient of variation (CV) of the density, ${\sigma }_{\rho }/\overline{\rho }$, over time as a function of height above the midplane. This coefficient is a measure for the strength of the density perturbations at a given height. For the calculation, we used the same two radii as for the calculation of the azimuthal Fourier decomposition. As the residual motions for the two different CO isotopologues suggest, perturbations in the density strongly decrease in layers below three pressure scale heights in the inner disk, and only reach substantial values in the outer disk during the initial cloudlet encounter (phase 1). CV values above 0.1 are not reached below four pressure scale heights in the inner disk after phase 2, when the m = 2 spiral arises. This means that the infall-induced density fluctuations below this height are lower than 10% of the mean density. The effect of a deeper perturbation on the emergence of spiral patterns seems to switch between the inner and outer disk. In the inner disk, the m = 2 spiral is only visible in phases 3 and 4, were the perturbation depth is minimal. In contrast, in the outer disk, the m = 2 spiral disappears for the minimal depth in phase 4, but is strong during the deeply perturbed phases 2 and 3.

Fig. 6 Coefficient of variation $\left({\sigma }_{\rho }/\overline{\rho }\right)$ in azimuth as a function of time and height above the disk midplane. The dashed white lines denote the different phases of the encounter. The black lines represent three (solid), four (dotted), and five (dash-dotted) pressure scale heights. The coefficient is calculated for r = 25 AU and r = 70 AU, shown in the top and bottom panel, respectively. The dashed red, orange, and yellow lines denote the median in time of the τ = 2/3 emission surface of 12CO and 13CO, respectively.

Fig. 7 Volume-weighted averages of the orbital elements of the gas as a function of time and semimajor axis. The left panel shows the eccentricity, while the middle panel shows the argument of periapsis. The vertical dashed lines show the different phases of the encounter. The right panel shows the elliptical orbits for each semimajor axis bin, viewed face-on, using the e and ω values from the other two panels.

3.1.4 Effect on orbital elements

In an effort to shed more light on the mechanism that leads to the formation of the spirals, we calculated the orbital properties of the gas in the directly influenced upper layers. We determined the eccentricity and argument of periastron as a function of semimajor axis, which was split into 100 bins. Values from cells whose semimajor axis values lie within the same bin were averaged, and values were weighted by cell volume. We chose not to use the density or mass of the cell as weight, as this would favor the unperturbed disk midplane. As we are mainly interested in surface-level effects, using the cell volume as weight is a choice that puts emphasis on the upper disk layers. To prevent cells that are unrelated to the disk itself from influencing the calculated average, we only considered cells that are within 4 pressure scale heights from the disk midplane.

The resulting values, as a function of time, are shown in Fig. 7. We find an excitation of eccentricity in the outer disk, up to e ∼ 0.35, once the fallback streamer emerges. Similar to recent results by Calcino et al. (2025a), where twisted elliptical orbits acted as an explanation for the spirals arising in their model of cloudlet capture, we find a shift in the argument of periastron at the locations with excited eccentricity. However, unlike their findings, the resulting elliptical orbits do not align in a way that would promote the formation of low-modal spirals for the majority of the time. The most prominent, but brief, spiral structure that forms as a result of that mechanism is shown in the right panel of Fig. 7. It is an m = 1 spiral in the inner disk, which is not related to any structures visible in scattered light. This indicates that the spirals found in our simulation are a direct consequence of the late infall accretion, and not a consequence of perturbed disk dynamics. We note that, in principle, the induced eccentricity subjects the upper disk layers to a parametric instability, which could drive vertical flow patterns and turbulence (Ogilvie & Barker 2014; Dewberry et al. 2025).

Fig. 8 Same as Fig. 1, but for simulation 2. The white arrows denote the number of spiral arms m determined by eye. For the right panel, the dotted arrow denotes an alternative interpretation where the m = 1 spiral is an m = 2 spiral, but one arm is much fainter.

Fig. 9 Same as Fig. 2, but for simulation 2, where we model a BHL accretion scenario.

3.2 Accretion in a turbulent medium

When late infall occurs in the form of BHL accretion, that is, when gas is accreted from a turbulent environment, the emerging spiral structure is considerably different from that in the cloudlet capture scenario. Only two different kinds of structures arise, which we again used to divide the evolution into two phases, shown in Fig. 8 using two representative snapshots². In phase 1, a faint Bondi-Hoyle accretion tail is present, but no separate streamer has formed yet. Here, an m = 2 spiral can be seen, especially in the outer disk. The distinct spiral arms are marked in Fig. 8. Once a streamer has formed, which we classify as phase 2, the inner disk m = 2 spiral can no longer be seen due to the perturbations at the disk surface, and the scattered light from the streamer arms make it difficult to discern clear spiral arms by eye. A spiral structure with a low modal number remains in the outer disk, although the exact dominating mode is also hard to discern by eye: It could be a one-armed spiral, or a two-armed spiral where one arm is considerably dimmer than the other.

Analogously to the analysis of the cloudlet capture case, we analyzed the strength of the individual spiral modes in Fourier space, shown in Fig. 9. Here, we identified the two evolutionary phases mentioned above. The initial m = 2 spiral in the outer disk during phase 1 is less strong in the inner disk, and frequently disappears. The inner and outer disk exhibit a strong m = 1 mode, although we suspect this to be caused by the contamination of the scattered-light image by the Bondi-Hoyle tail. During phase 2, when the arced streamer is present, the Fourier decomposition is dominated by the m = 1 mode in the inner disk, and the m = 2 mode in the outer disk. Since the streamer appears fainter further from the star, we conclude that the spiral structure present in the outer disk at that time is likely to be truly dominated by an m = 2 mode, and not a result of contamination. However, similar to phase 1, we suspect contamination by the direct signature of the late infall to play a significant role for the inner disk. To examine this further, we explored whether this m = 1 mode may indeed be a disk-inherent structure in the subsequent sections.

A fundamental difference of this accretion scenario compared to cloudlet capture is the continued excitation of higher order modes. This is caused by the sustained presence of a streamer, unlike the cloudlet capture scenario, where the streamer only persist for a short time, and eventually dissipates, leading to accretion of the gas remnant. As a result, the Fourier decomposition of the azimuthal intensity profile calculated for the BHL accretion case is more prone to contamination by the streamer.

Fig. 10 Moment 1 residuals from Keplerian motion of the 12CO and 13CO isotopologues for simulation 2. Top: same as the top row of Fig. 4. Bottom: same as the top row of Fig. 5. The arrows denote the number of spiral arms m determined by eye.

3.2.1 Gas kinematics

We further investigated how the gas kinematics are affected in this scenario. Analogously to Section 3.1.2, we computed the residual motion found in the moment 1 map of the ¹²CO and ¹³CO line emission (see Fig. 10). For phase 1, we find that the innermost disk region exhibits an m = 1 spiral in the gas kinematics, which is not seen in the Fourier decomposition of the scattered light (see Fig. 9) due to the perturbations introduced by the ongoing infall. It is visible in the emission of both isotopologues, indicating that it is not only the result of direct streamer emission. It is instead a result of a kinematic perturbation induced by the streamer. This is apparent because the perturbation in the ¹³CO emission spatially matches a perturbation in the ¹²CO emission, originating from the upper layers, and this perturbation is connected to the direct emission of the streamer. We also find an m = 1 spiral in the outer disk regions, which spatially matches one of the arms found in the scattered light (see Fig. 8).

We find the most interesting kinematic structure for phase 2. Here, we see strong residual motions in the inner disk that coincide with the location of the brightest scattered-light structures. The inner disk spirals exhibit an m = 2 pattern in the CO emission, which is especially visible for ¹³CO. The outer disk regions show residual motions that are akin to an m = 1 spiral. These findings are the direct opposite of the apparent scattered-light structures, where the Fourier decomposition shows a strong m = 1 mode is present in the inner disk, and a strong m = 2 mode is present in the outer disk (see Fig. 9, where distinct spiral arms are marked with arrows). A determination of the number of spiral arms by eye suggests that there is an m = 1 spiral in the outer disk, or an m = 2 spiral where one arm is very dim, which is matched by the kinematic analysis. Considering the CO lines emission, it is easier to discern observable features caused by the streamer from those exhibited by the disk itself due to jumps in the velocity residuals. The emission of molecular lines is also not directly affected by shadowing (but it is affected indirectly via changes in the temperature). This allows us to conclude that streamer contamination and shadowing are likely explanations for the differences between structures seen in the CO emission and the dominant modes of the Fourier-decomposition of the scattered light.

In contrast to the cloudlet capture case, we find no meaningful structure in the residual motion of the CO isotopologues when the disk was viewed at an inclination of 30^◦, which is shown in Appendix A. However, more careful modeling of the disk motion that is less susceptible to issues related to environmental emission and non-smooth changes in the emission height might allow for the discovery of residual motions in the more inclined case. Including photodissociation would remove emission from the environment and could further improve the analysis of inclined disks. This is subject to future work.

3.2.2 Perturbation depth

The depth toward the midplane of the perturbations induced by the infall in this scenario is overall comparable to what we find in the cloudlet capture case (Section 3.1.3), see Fig. A.1. As expected, the deepest perturbation occurs during phase 2, when material is actively accreting in the form of a streamer, reaching down to z ≈ 3H in the outer disk. Differently from the cloudlet capture case, the perturbation depth does not vary greatly in time, because the accretion does not stop during the simulation. The minimal height above the midplane for which the density CV is higher than 0.1 is shallower than during the cloudlet impact. This indicates that the disk density structure is not noticeably affected by the infall beyond the upper surface layers.

3.2.3 Orbital gas elements

We analyzed the orbital elements in the same way as in the cloudlet capture case. The accretion of material leads to an excitation of the eccentricity in the surface layers, and a simultaneous shift in the argument of periastron, shown in the left and middle panel in Fig. 11. There are several key differences to the excitation of eccentricity in the cloudlet capture scenario. First, the maximal excitation is slightly lower, reaching e ∼ 0.25, present during phase 2 when the streamer is active. Second, the semimajor axes of the excited gas are limited to a smaller range of values, close to 60–70 AU. Lastly and most importantly, the alignment of the elliptical orbits caused by the shift in the argument of periastron creates an m=1 spiral pattern, shown in the right panel in Fig. 11, similar to the findings by Calcino et al. (2025a).

Even though we find a similar alignment of the gas orbits, we find the resulting pattern to be unrelated to the spiral patterns that arise in the scattered-light synthetic observation. While the strong m = 1 mode in the inner disk (see Fig. 9) may initially suggest the presence of a spiral created through this mechanism, its shape and location do not match the scattered-light images found in Fig. 8. Furthermore, the ellipse pattern is already present during phase 1, where the m = 1 component of the Fourier decomposition is weak in the inner and outer disk, and where no visual match to the pattern exists. The disk eccentricity might, in principle, also cause the spiral patterns found in the residual motions, motivated by the idea that gas in an eccentric disk exhibits vertical motion due to a breathing mode (Ragusa et al. 2024). Indeed, in the CO residual motion of the two considered isotopologues, various m = 1 patterns arise, but they are of a different shape and at a different location. We therefore conclude that, while the orbits are affected in the same way as found by Calcino et al. (2025a), the spiral patterns we find in the scattered light and residual motions are a direct consequence of the late infall, rather than an alignment of eccentric orbits. The disturbances at the top layers of the disk due to the infall act to conceal any structures that the orbital shift might create, and the lower disk layers traced by the less abundant isotopologues are not sufficiently perturbed.

Fig. 11 Same as Fig. 7, but for simulation 2.

3.3 Lower disk mass

In previous sections, we found that the induced spiral structures are mainly caused by the disturbances in velocity and density in the upper layers of the disk. The reason for the shallow influence of the infall on the disk is the comparatively low mass of accreted material compared to the disk mass. We therefore expect the mass of the disk itself to have a major impact on our results, as a lower-mass disk would be subject to deeper perturbations. In this subsection, we investigated the two infall scenarios, a cloudlet capture scenario and gas accretion from a turbulent medium, with a reduced disk mass by a factor of 10. We summarize the most important differences in the following. A more detailed analysis is given in Appendix B.

As expected, the disks are more deeply perturbed when the disk mass is lower, which is more pronounced for the turbulent accretion scenarios. This is shown in Fig. 12, where we show the lowest value of z/H for which the CV of the density is at least 0.5. The figure can also act to quantify the results discussed in previous sections: The perturbations by BHL accretion are deeper at later times, but shallower than the perturbation by the cloudlet capture during the time when the cloudlet encounters the disk. This result generally holds true for the lower disk mass case, but here, the BHL perturbation is deeper than the cloudlet capture one at earlier times. Specifically, it occurs at t ∼ 15 kyr in the inner disk and t ∼ 20 kyr in the outer disk for the lower disk mass, but at t ∼ 20 kyr in the inner and t ∼ 25 kyr in the outer disk for the regular disk mass cases. For the lower disk mass, the deepest perturbation is 2.5 pressure scale heights for BHL accretion for the inner and outer disk, although this occurs at different times for the different disk regions. In the case of cloudlet capture, the disk with lower mass is perturbed down to 3 pressure scale heights in the inner, and 2.5 pressure scale heights in the outer disk.

A natural consequence of the change in perturbation depth is a corresponding change to the arising spiral structures in the polarized scattered-light intensity, especially for the low spiral modes. We provide a comparison for all four simulations in Fig. 13. Least affected by the change in disk mass is the relative amplitude of the m = 1 mode in the inner disk. Here, the amplitude is diminished compared to the regular mass disk around t ∼ 20 kyr for the BHL accretion case, but enhanced around the same time for the cloudlet capture case. More significant differences occur in the outer disk. Here, the amplitude of the m = 1 mode is strongly enhanced for the cloudlet capture case at t ∼ 20 kyr, and again at t ∼ 30 kyr. The increase is less severe, but present for the turbulent medium accretion simulation, at t ∼ 15 kyr.

We find that for the m = 2 mode, a decrease in disk mass is detrimental to the relative amplitude. This is especially true for the cloudlet capture case, where the mode disappears fully in the inner disk. In the outer disk, while not disappearing fully, it is strongly decreased for the cloudlet capture case and the scenario where gas is accreted from the turbulent medium.

Fig. 12 Perturbation depth in units of the pressure scale height as a function of time for all four simulations. The left panel shows the depth for the inner disk (r = 25 AU), and the right panel shows the depth for the outer disk (r = 70 AU). The dashed lines represent simulations with regular disk mass, solid lines ones with low disk mass. The blue lines describe cloudlet capture simulations, whereas the red lines describe accretion in a turbulent medium.

Fig. 13 Normalized Fourier-transformed intensity of the m = 1 (left panels) and m = 2 (right panels) modes of the polarized scattered-light intensity of all four simulations. The top panels show the intensity in the inner disk (r = 25 AU), whereas the bottom panels show the intensity in the outer disk (r = 70 AU).

4 Discussion

4.1 Spiral formation mechanism

In this work, we show that the two considered infall mechanisms, that is, cloudlet capture and BHL accretion, can cause spiral structures, and that the m = 2 mode dominates the arising pattern. This result deviates from the results of previous studies, which mainly found the m = 1 mode to dominate in the case of asymmetric accretion (Hennebelle et al. 2017; Calcino et al. 2025a). A more detailed analysis of the physical origin of the spiral structures is subject of future work, but in the following, we discuss the main potential spiral formation mechanism in our simulations, and the main differences to previous work.

Without consideration for secondary mechanisms, the effects of arbitrary perturbations to any protoplanetary disk quantity on its dynamics can be studied by performing a decomposition into a series, for example, azimuthal Fourier modes or vertical Hermite modes (see, e.g., Zhang et al. 2025). A perturbation whose decomposition is dominated by an m = 1 Fourier mode would then be associated with creating an m = 1 spiral pattern. To first order, this is the reason for the creation of this pattern in previous works, where the infall is highly asymmetric. The infall originates from one direction and interacts with the side of the disk facing that direction, either fully in the case of extended infall, or in a spatially confined way when the infall takes the form of a streamer. Other mechanisms can create m = 1 spirals, for example, a disk warp (Winter et al. 2025) or the excitation of disk eccentricity and a twist (Calcino et al. 2025a). However, we find that these mechanisms do not cause the spirals arising in our simulations.

No such strong asymmetry in the interaction with the infalling material exists in our simulations. This is a consequence of the fallback of material in the cloudlet and BHL case (see HD25). In the simulations of cloudlet capture, the expanded cloudlet interacts with the disk predominantly from one side during the initial encounter, and the newly bound, initially overshooting material falling back interacts with the disk predominantly from the other side. This double-sided interaction has a strong m = 2 mode in its Fourier decomposition, leading to m = 2 spirals. These spirals are not identifiable in synthetic observations during the entire interaction. The reason for this is that the disk surface is perturbed by the accretion, in the form of sound- and shockwaves. This is especially relevant during the initial interaction between the cloudlet and the disk, and during the period where the fallback streamer is active. After these strong interactions have subsided, the spirals are more clearly visible; they are still caused by font- and backside interactions with the leftover cloudlet gas, but the accretion rate is lower and corresponding surface perturbations are less relevant. The two-sided interaction being the main driver of the m = 2 spiral morphology also offers an explanation for the result that the spiral arms of the cloudlet-induced m = 2 spirals are almost stationary: The spirals are launched from the interaction sites between disk and infall, which do not change drastically over the course of the simulations. This matches the previous findings by Calcino et al. (2025a) and Hennebelle et al. (2017).

The mechanism forming m = 2 spirals in the BHL accretion simulations is similar in concept, but the details are different. The two-sided nature of the interaction is given by the systemic flow on one side, and the fallback of material from the accretion tail on the other side. While a turbulent medium is necessary to create fallback streamers (see HD25), turbulence is not a required ingredient to form m = 2 spirals, as the two-sided interaction arises due to the fallback of material in the accretion tail. When turbulence is included in the consideration, the m = 2 spiral formation mechanism does not change, because the deviations it introduces in the density and velocity fields of the systemic ISM flow are small and do not have a significant impact on the fallback mechanism. Therefore, this scenario is best compared to the ASNR2 case of Hennebelle et al. (2017), where inflow is asymmetric but nonrotating. While not discussed in the original work, this simulation also shows a weak m = 2 spiral pattern. In our study, we did not perform a pattern speed analysis for the BHL accretion simulations, because the surface is too strongly perturbed in the scattered light to allow for a proper characterization of the spiral location. However, given that the spirals presumably form through the same mechanism, we expect the spirals to also be almost stationary in the BHL accretion case.

In summary, the mechanism forming m = 2 spirals in our simulations is the two-sided interaction of infalling material with the disk, caused by the fallback of newly bound material. The essence of this mechanism is the same for both accretion mechanisms, although their observational implications are different. In the cloudlet capture case, the main interaction leads to stronger perturbations, but the corresponding spiral patterns are difficult to discover due to the corresponding surface perturbations. The interaction with the leftover gas after the main encounter is favorable for the detection of spirals, because the accretion occurs without creating a fallback streamer, but it is limited in time. On the other hand, BHL accretion creates more persistent interactions, solving the issue of a short interaction time compared to the disk lifetime. At the same time, the persistent interactions lead to the continuous creation of fallback streamers, and the corresponding constant surface perturbation make the detection of spirals more challenging in observations.

4.2 Implications for the observational interpretation of spirals

We find a large variety of spirals arising as the result of surface-level disturbances solely caused by late infall. In particular, we find that prominent m = 2 spirals can arise in scattered light, and m = 1 spirals in the CO residual motions, which are unrelated to other mechanisms commonly invoked to form spirals. The infall that causes these structures also does not necessarily need to be visible in the form of a streamer. The most prominent scattered-light spiral structure arises in Fig. 1, at a time when only a faint remnant of the original cloudlet is accreted after the main infall event, with an active fallback streamer, has passed. With this in mind, special care should be taken when inferring disk conditions or planet candidates from the existence of spiral structures alone. When no additional information, such as a considerable pattern speed or perturbations in deeper layers, is available that can exclude a late-infall event causing spirals, an observed spiral might be mistaken for evidence of a perturber, GI, or a warped disk that might not reflect the physical reality. For example, a warped and twisted disk was recently found to be a possible explanation for the spirals in the ¹²CO line emission in MWC 758 (Winter et al. 2025); however, none of our spiral structure, neither in scattered light, nor in CO, is related to a warp structure (see Appendix D).

Deep perturbations could act as evidence for the formation of spirals through other mechanisms, because late infall does not drive them directly unless a lot of mass is accreted (see Appendix B). It remains to be investigated to which extent this statement holds true when considering that the infall could trigger secondary mechanism, like the Rossby wave instability (Bae et al. 2015; Kuznetsova et al. 2022), or the excitation of turbulence (e.g., Fleming & Stone 2003). It is currently unclear whether the Rossby wave instability might be triggered by the infall when mass is not directly deposited at the midplane, as 3D studies of it assume hydrostatic equilibrium (e.g., Meheut et al. 2012), which does not hold in this case.

Based on the patterns we see in our simulations, we predict that structures like in SU Aur (Ginski et al. 2021) should be most common for disks under the influence of late infall: A clearly visible streamer with flocculent spiral structures. Especially in the inner disk, we find that the perturbation caused by streamers prevent the formation of spirals with strong low-numbered Fourier modes, so that we would not expect to see such spirals at the same time as a streamer. This conclusion holds for observations of scattered light; even in the abundant CO isotopologues, low-m spirals can exist during active infall streamers. In general, the spiral morphology in scattered light and CO line emission could be used to infer information about the environmental conditions causing streamers of material, by distinguishing between accretion from a turbulent ISM or accretion caused by a single encounter with a cloudlet-like overdensity.

The spiral patterns that we discover in the CO line emission kinematic residuals are the clearest when the disk is viewed almost face-on, with only a small deviation from a perfect 90^◦ angle being introduced by the infall. This indicates that the motions causing the spiral patterns have a strong vertical component. We find that, during the accretion of the cloudlet remnant, the spirals are not visible when the disk is viewed at an inclination of i = 30^◦. They are only visible while the streamer is active in the cloudlet capture case. However, this could be related to difficulties in fitting a Keplerian model to the emission, caused by the abundant emission of the surrounding material. More detailed modeling, or the consideration of the photodissociation of CO which reduces emission in the surrounding cloud, could allow for the detection of spirals also in the inclined case. If this is indeed the case, we would expect the perturbations at the surface level, caused by the turbulent flow, to differ between the front and back emission surfaces, which could act as another way to distinguish infall-related perturbations from those created by other mechanisms.

As discussed in the previous sections, there are two key factors that allow us to distinguish the spirals caused by late infall from spirals caused by other mechanisms: The pattern speed and the dependence of the structures on the height above the mid-plane. In some cases, the shape of the spirals themselves can also aid in the distinction. For example, depending on the cooling rate, GI causes flocculent structure in the outer, and spirals with fewer arms in the inner disk (Rowther et al. 2024). This is the direct opposite of infall-induced spirals, where flocculent structure is found in the inner disk, and m = 2 spirals are found in the outer disk.

4.3 Implications for disk dynamics and planet formation

The presence of spiral substructures in the disk raises the question whether they could aid planet formation. On the one hand, infall generally delivers new solids that can aid, or even rejuvenate, planet formation. In addition, spirals can accumulate dust and create density enhancements that would be beneficial, for example for the formation of planetesimals or for retaining solids for pebble accretion. For this to be effective, trapping needs to occur at the disk midplane, but we find that the midplane is not directly affected by the infall; the spiral structure does not penetrate deeply enough into more massive disks. For these disks, that is, for M_d ≳ 0.05 M_⊙, we therefore conclude that the induced spirals do not influence planet formation neither positively, nor negatively.

This finding does not come without a model-based caveat. The presence of a strongly turbulent disk surface, excited by the infall, could drive weak turbulence in the midplane (e.g., Fleming & Stone 2003), which would be missed in our simulations. Numerical viscosity introduced by the misalignment between the gas velocity and the cell surfaces is not a significant problem due to our choice of a spherical grid, but numerical diffusion might play a significant role. Even if this was fully negligible, low levels of turbulence cannot be resolved in our simulations because the eddy scale, $\sqrt{\alpha }H$, is smaller than the cell size. Resolving the excitation of turbulence with α = 10⁻² requires a cell size of 0.1H, or 10 cells per scale height, which is not reached inside the disk, but even turbulence below this value could influence planet formation. All in all, while our results indicate that the midplane is not directly perturbed by the infall, our setup does not allow us to make a definitive statement about whether planet formation could be aided by the secondary excitement of turbulence via the turbulent upper layers.

This result strongly depends on the ratio between the accreted mass through infall and the mass of the already present proto-planetary disk. We find that, for a lower initial disk mass but the same accreted mass, the perturbations reach deeper into the disk (see Appendix B). When the infall perturbs the disk deeply enough to affect the midplane structure, dust traps can be created, which affects planet formation. This would be more likely to occur at later stages of the disk, when most of the mass has been depleted. Alternatively, a higher infall rate onto a disk of unchanged mass would also affect layers close to the mid-plane, but might lead to the disk and surrounding cloud being interpreted as a Class I source. This suggests that late infall may be especially important to consider toward the end of the protoplanetary disk lifetime.

4.4 Model caveats

Our simulations only considered the evolution of the gas in the disk and its surroundings, and we made no attempt to model the transport or evolution of any solids. This was motivated by the computational cost that would be introduced by integrating additional components, especially when multiple grain sizes and grain growth are to be treated. For the large-scale streamers themselves, it is reasonable to assume that the dust is small and therefore well coupled, so that an inclusion of the dust is not very impactful. While assuming perfectly mixed and small grains is also suitable for the consideration of scattered light and CO line emission images at sufficient emission heights, modeling of the dust size distribution, including coagulation, and dust-gas interactions is necessary for a more detailed investigation of the impact that late infall has on dust transport, dynamics, and planet formation. Even though we only find surface-level perturbations in the gas, the growth, settling, drift, and vertical mixing could, in principle, connect the perturbed surface with the seemingly unperturbed midplane and impact planet formation directly.

We modeled the BHL accretion by introducing compressible turbulence to the large-scale environment around the disk as an initial condition, but without dynamically driving the turbulence. As a result, the turbulence decays over time in our model, whereas in reality, the strength of the turbulence at a given time is a reflection of the environmental conditions in the region the disk moves through. The decay is more significant for the smaller scales, as scales smaller than the cell size cannot be resolved, and with our choice of spherical coordinates with log-radial spacing, the cell size increases with increasing distance from the star. Flow that originates from these further-out regions, reaching the disk at later times, has a turbulent power spectrum that is suppressed at small scales (higher k values). We do not expect that this decay has a significant impact on our results, but we note that the emergence of spirals that have low-m Fourier modes with a high amplitude could be related to this artificial suppression of small-scale turbulence at that time. As we utilized a simple isothermal temperature model in our simulations, we also did not treat the heating that would occur as a result of the shocks and turbulent compression in the environment. The former is likely to be result in substantial heating as streamers hit the disk surface (van Gelder et al. 2021). In the future, it should be investigated how this heating affects the spiral structures.

We did not consider the effects of freeze-out and photo-dissociation when computing the CO abundances used for the creation of the synthetic CO line emission images. While we do not expect these effects to have a significant impact on the abundances inside the disk and at its surface, the resulting optical depth of the cloud surrounding the disk is likely overestimated as a result. As a result, it has a stronger effect on the discminer Keplerian fit, and fits performed on a real observation might have more success with discovering residual motions without introducing emission-height related model residuals. As we only resolve the disk with a few cells per pressure scale height in the vertical direction, we might also incur inaccuracies in the effective height above the midplane a particular CO isotopologue is emitting from.

5 Conclusions

We performed 3D hydrodynamics simulations and radiative transfer modeling to investigate the emergence of spiral structures as a result of late infall. For our investigation, we used the simulation setups leading to streamers from Hühn & Dullemond (2025), considering the capture of a gaseous cloudlet and the accretion of gas from a turbulent medium. We draw the following conclusions:

• In the cloudlet capture scenario and considering the scattered-light emission, the inner disk is strongly perturbed when the cloudlet hits the disk and a streamer forms, leading to a flocculent structure. In contrast, m = 2 spirals can form in the outer disk at this time, which only occurs in the inner disk after the initial streamer has subsided. This scenario generally favors the formation of well-defined spiral arms after the infall has largely subsided;

• The pattern speed of the m = 2 spiral during the remnant accretion after cloudlet capture is low, ranging between $\dot{\varphi }=0.05{\text{kyr}}^{-1}$ (r_corot = 3150 AU) and $\dot{\varphi }=0.11{\text{kyr}}^{-1}$ (r_corot = 1460 AU), depending on the considered radius and spiral arm;

• In the scenario in which the disk accretes mass from a turbulent medium, where the infall rate fluctuates but does not fully subside, the formation of surface-level m = 2 spirals, visible in scattered light, is generally not favored. Here, the disk shows more flocculent structures during the simulation, but before the first streamer arises, an m = 2 spiral can form;

• Spiral structures can also be seen in ¹²CO line emission through residuals from Keplerian motion, but they do not match the scattered-light emission patterns. They are best seen when the disk is viewed face-on, suggesting that they are caused by vertical motion. The amplitude of their individual angular modes also differs from the scattered-light case;

• The perturbations induced by the infall are mostly limited to the upper layers of the disks, to about four pressure scale heights. This is directly related to the smaller residual motions found in the emission of the more optically thin ¹³CO;

• While changes to the gas orbital elements can be induced for some time in the upper layers, the m = 1 spiral pattern induced for the turbulent medium accretion case is not visible in scattered light or in molecular line emission;

• Perturbations reach deeper into a disk when its mass is lower compared to the infall, which is expected. This affects the spirals that emerge in scattered light in many aspects; most notably, it disfavors the emergence of an m = 2 spiral during remnant accretion.

Our results demonstrate that the observed spiral structures do not necessarily need to be induced by a clear perturber, even when the spiral arms are very well defined. Instead, they emerge naturally as a result of infall. Despite this, they can be discerned from spirals formed through other mechanisms through their low pattern speed. More generally, the spiral perturbations do not reach the midplane of the disk, unless its mass is low compared to the accreted mass, which is the most relevant location for planet formation. Spirals formed solely through infall, without considering any additional mechanisms they might trigger, are therefore not expected to greatly affect planet formation beyond the delivery of new material, especially at the beginning of the Class II stage.

Acknowledgements

We thank the anonymous referee for suggestions that helped improve the manuscript. We thank Haochang Jiang and Josh Calcino for stimulating discussions. L.-A.H. acknowledges funding by the DFG via the Heidelberg Cluster of Excellence STRUCTURES in the framework of Germany’s Excellence Strategy (grant EXC-2181/1 – 390900948). 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 grants INST 35/1134-1 FUGG, 35/1597-1 FUGG and 37/935-1 FUGG.

References

Appendix A Additional figures for turbulent medium accretion (high disk mass)

Fig. A.1 Same as Fig. 6, but for simulation 2.

Fig. A.2 Residual motion found in CO isotopologue line emission for simulation 2 at an inclination of i = 30◦. Left: Same as bottom row of Fig. 4. Right: Same as bottom row of Fig. 5.

In this Appendix, we present additional figures for simulation 2, which were not discussed in Section 3.2. In Fig. A.1, the depth of the perturbations induced by the BHL accretion is shown. Similar to the cloudlet capture case, the perturbations do not reach deeper than three pressure scale heights in the inner disk, and only does so very briefly for the outer disk. This allows us to conclude that, for a disk mass of M_d = 0.05 M_⊙, any discovered substructure is limited to the surface for both infall mechanisms, even when assuming initial conditions or infall rates that are on the higher end of realistic values.

In Fig. A.2, we present the residuals of the moment 1 line emission of ¹²CO and ¹³CO for an inclination of i = 30^◦. Unlike the face-on case, no considerable substructure can be seen. We note that, due to the optical thickness of the surrounding turbulent medium and the strong disturbance of the disk kinematics in the upper disk layers, fitting a Keplerian model becomes increasingly difficult within the capabilities of the discminer code, creating severe unphysical residuals, not allowing any conclusions to be made for an inclination of i = 30^◦. In addition, a straight line of residuals can be seen at x = 0. This is caused by missing emission in the 0 km s⁻¹ velocity channel in the synthetic observations, which is the result of the turbulent medium surrounding the disk. Due to the simulation setup, it has no velocity toward the line of sight, so that emission in that channel is dominated by the background emission. A different camera angle would shift this strong residual to a different channel, but would not fix the issue. In reality, the freeze-out of CO might make the line emission less optically thick than we assume (see Section 4). For ¹³CO, the surrounding medium does not obscure the disk as strongly.

Appendix B Additional figures for simulations with low disk mass

In this Appendix, we present additional Figures for simulations number 3 and 4, which have a lower disk mass, making them more susceptible to deeper perturbations. In Fig. B.1, we show the scattered light synthetic images of a selection of the snapshots considered for the regular disk mass simulations. The visible structures match the findings presented in Fig. 13: For the cloudlet capture, m = 2 spirals disappear in the inner disk and become weaker in the outer regions, whereas the m = 1 mode becomes stronger, especially in the outer disk. For the BHL accretion, the m = 2 mode actually becomes more strongly visible during phase 1, but structures mostly disappear during phase 2. In general, the observed structure is comprised less of flocculent spiral arms. To test to what extent the change in visible structures is related to the change in the scattering height, we also created synthetic scattered light observations with a dust-to-gas ratio higher by a factor of 10 compared to the regular value. Even though this is an unrealistic value, it allows us to reproduce the solid mass assumed in simulations 1 and 2. The results are shown in Fig. B.2, showing that the change in emission height does not account for all differences found in the lower disk mass case, but it can account for some changes in visible structure. For the cloudlet case, the spiral structure seems more dominated by the m = 2 mode for the higher dust-to-gas ratio during phase 2, as is the case for simulation 1, and some spiral structure is recovered in the inner disk. For the BHL accretion case, the m = 2 mode during phase 1 disappears with the higher dust-to-gas ratio, which also holds true for the flocculent variations during phase 2.

Fig. B.1 Synthetic scattered light observations of structures found in the simulations with a lower disk mass. Left: Same as zoomed-in area of panels 2 and 3 in Fig. 1, but for simulation 3. Right: Same as zoomed-in areas in Fig. 8, but for simulation 4.

Fig. B.2 Same as Fig. B.1, but with an artificially high dust-to-gas ratio of 0.1, reproducing the dust mass used for simulations 1 and 2.

Fig. B.3 Same as the left panel of Fig. 7, but for simulation 3 (left) and simulation 4 (right).

In Fig. B.3, we show the eccentricity of the gas for the cases of lower disk mass. The locations of excitation in time and semimajor axis do not change significantly, though for the cloudlet case, the excitation of eccentricity is stronger. This is an expected result, as the momentum of the infalling material has a greater impact on the total momentum of the gas. For the BHL case, the excitation is at a similar level, or lower at a few points in time, but the inner disk becomes more excited.

Fig. B.4 Same as Fig. 4 (left) and Fig. 5 (right), but for simulation 3.

Fig. B.5 Same as Fig. 4 (left) and Fig. 5 (right), but for simulation 4.

In Fig. B.4, the residuals of the gas motion as determined from the first moment map of the CO isotopologue line emission are shown for the cloudlet capture case. In the face-on images, m = 1 spiral structure dominates for phases 3 and 4, and they can be seen with more clarity than for the regular disk mass case. For the panel showing the phase 2 snapshot, the structure is less obvious and spirals are more flocculent. Similar to the high disk mass case, the issue of fitting a Keplerian model remains and is further amplified by the stronger perturbations that the disk with lower mass is subject to, creating further severe unphysical residuals in the i = 30^◦ images. We show the same kind of images for the BHL accretion case in Fig. B.5, and find an analogously increased prominence of low-mode spirals. Most notably, the m = 1 spiral structure is very clear in this case, especially during phase 1.

The depth of the perturbations induced in the lower-mass disks are shown in Fig. B.6, providing more context on the analysis based on a threshold value shown in Fig. 12. In line with expectations, perturbations protrude further toward the midplane. In the cloudlet case, even the midplane can be weakly affected directly at the time of the main encounter, but they retract to the surface for the later accretion stages, affecting only layers at z ⪆ 3H. For BHL accretion, deeper layers are affected only at later times, but for significantly longer compared to the cloudlet case. Perturbations occur close to the midplane during the majority of phase 2.

Fig. B.6 Same as Fig. 6, but for simulation 3 (top) and simulation 4 (bottom).

Appendix C Inclined infall

Fig. C.1 Left: Same as zoomed-in areas of Fig. 1. Right: Same as zoomed-in areas of Fig. 8.

In this Appendix, we present synthetic scattered light observations of two simulations where the infall direction is inclined relative to the disk plane. They are shown in Fig. C.1. The cloudlet capture simulation was run with the same parameters as simulation 1, but the position vector and initial velocity of the cloudlet were rotated by 30^◦ around the x-axis. The BHL accretion simulation was run with the same parameters as simulation 2, but the systemic velocity unit vector is ${\widehat{v}}_{\text{sys}}=\left(0,\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}\right)$. In the cloudlet case, m = 2 spirals arise during the same phases as in the case of in-plane cloudlet capture, but they are less prominent during phase 3. We suspect that this is because the infall perturbations are more symmetrically spread across the disk surface, diminishing the effect. In the limit of perfectly face-on accretion event, we would not expect that spirals form directly. In this case of inclined infall, the disk is also more susceptible to other spiral formation mechanisms, like the excitation of eccentricity (Calcino et al. 2025a), or disk warping (Winter et al. 2025), which could contribute to the visible spiral patterns. A detailed analysis of this scenario is outside the scope of this work. The accretion tail found in the BHL accretion simulation has an orientation relative to the observer that is unfavorable for the study of spiral patterns, though infall streamers are more easily observed (see HD25).

Appendix D Warp structure

Previous work has shown that a warped and twisted disk can create spiral structures, and that this mechanism can explain aspects of the spiral structure in MWC 758 (Winter et al. 2025). While our simulation setups was not chosen to promote the formation of a disk warp, we nevertheless investigated whether a warp could be the cause of the substructures we find in this Appendix. Figures D.1, D.2, D.3 and D.4 show three parameters that are commonly considered for the characterization of a warped disk structure for all four simulations. These parameters are the disk tilt β, $$\beta (r)=\arccos \left(l_z(r)\right),$$(D.1)

the twist γ, $$\gamma (r)={\arctan }_2\left(l_y(r),l_x(r)\right),$$(D.2)

and the warp strength ψ, $$\psi (r)=|\psi (r)|=r\left|\frac{\partial l(r)}{\partial r}\right|.$$(D.3)

Fig. D.1 Warp properties as a function of radial distance and time, where a darker color denotes a later snapshot. Left: Disk tilt in the reference frame described by the total angular momentum, β′. Middle: Disk twist in the same reference frame, γ′. Right: Warp strength ψ′ in that reference frame. All values are for simulation 1.

Fig. D.2 Same as Fig. D.1, but for simulation 2.

Fig. D.3 Same as Fig. D.1, but for simulation 3. We note that the range of values on the ordinate is increased compared to that Figure.

Fig. D.4 Same as Fig. D.1, but for simulation 4. We note that the range of values on the ordinate is increased compared to that Figure.

Here, l is the unit vector of the cell volume-weighted average angular momentum in a spherical shell at radius r with a thickness given by the grid resolution, ∆r. They are calculated in a different reference frame than that one given by our grid configuration; it is defined as the frame where the total disk angular momentum coincides with the z direction, so that the transformation is given by $$l^{\prime }=R_x\left(-{\varphi }_x\right)R_z\left(-{\varphi }_z\right)l,$$(D.4)

where R_i are rotation matrices about axis i, and the rotation angles are given by $${\varphi }_x=\arccos \left(l_{tot,\text{Z}}\right),$$(D.5) $${\varphi }_z=\arccos \left(-\frac{l_{tot,\text{y}}}{\sqrt{1-l_{tot,\text{z}}^2}}\right),$$(D.6)

with l_tot being the unit vector of the total disk angular momentum, integrated within 150 AU. The initial condition of the cloudlet capture simulations lead to a capture in the plane of the disk. We therefore do not expect a strong warp to occur, which we confirm with this analysis. In the regular disk mass case, the maximal tilt in the bulk of the disk is 0.1^◦. Here, there is no considerable twist, that is, γ does not change significantly with radius, and the maximal strength is ψ ∼ 10⁻². The warp is not much more significant for the lower disk mass case, where the maximal tilt is 0.8^◦, there is no twist, and the strength only reaches ψ ∼ 10⁻¹ in the outer disk.

For the simulations of BHL accretion, the systemic velocity is also in the disk plane, which does not promote warping of the disk, but the turbulent nature of the environment could, in principle, lead to more significant warping. However, we also do not find this to be the case. For the regular disk mass case, the tilt reaches β ∼ 0.4^◦ and the strength reaches ψ ∼ 10⁻². For the low disk mass case, the strength is of the same order of magnitude, but lasts longer, and the maximum tilt reaches close to β ∼ 1^◦. Neither of these simulations exhibit considerable twisting. We therefore conclude that warping is not the source of our CO line emission substructures, as expected. It remains to be investigated how out-of-plane accretion could cause more significant warp effects, which was found in previous works about late infall (e.g., Kuffmeier et al. 2021), and how that might cause the creation of substructures.

Appendix E Discminer fit parameters

All Tables

All Figures

Fig. 1 Polarized scattered-light intensity in four snapshots of simulation 1. The solid gray line denotes the 1 mJy arcsec−2 contour, the dashed gray line shows the 0.1 mJy arcsec−2 contour, and the dotted gray line shows the 0.01 mJy arcsec−2 contour. The panel numbering corresponds to the phases of the cloudlet encounter. The direction of the white arrow denotes the original orbital velocity of the cloudlet.

In the text

	Fig. 2 Azimuthal Fourier decomposition of the polarized scattered-light intensity, normed by the total intensity. $m=\tilde{\varphi }$ is the Fourier transform of the azimuthal coordinate, and the two panels show the transformation for two radii. The dashed white lines separate the phases of the encounter, labeled by the white numbers.
In the text

Fig. 3 Polarized scattered-light intensity as a function of azimuth and time for two different radii. The dashed white lines and numbers denote the different encounter phases, as in Fig. 2. The colored dots show the azimuthal peaks in gray when they were not considered further, otherwise by membership to one of the spiral arms. The solid white lines show the linear fit used to determine the pattern speed.

In the text

	Fig. 4 Moment 1 residuals from Keplerian motion of the 12CO line emission. The top row shows three snapshots, taken at the same time as in Fig. 1, where the camera is oriented face-on. The bottom row instead shows images where the camera has an inclination of 30◦. The black lines show polarized scattered light contours, as in Fig. 1. The panel numbers correspond to the evolutionary phase.
In the text

	Fig. 5 Same as Fig. 5, but for 13CO. We note that the value range of the color bar is reduced.
In the text

Fig. 6 Coefficient of variation $\left({\sigma }_{\rho }/\overline{\rho }\right)$ in azimuth as a function of time and height above the disk midplane. The dashed white lines denote the different phases of the encounter. The black lines represent three (solid), four (dotted), and five (dash-dotted) pressure scale heights. The coefficient is calculated for r = 25 AU and r = 70 AU, shown in the top and bottom panel, respectively. The dashed red, orange, and yellow lines denote the median in time of the τ = 2/3 emission surface of 12CO and 13CO, respectively.

In the text

Fig. 7 Volume-weighted averages of the orbital elements of the gas as a function of time and semimajor axis. The left panel shows the eccentricity, while the middle panel shows the argument of periapsis. The vertical dashed lines show the different phases of the encounter. The right panel shows the elliptical orbits for each semimajor axis bin, viewed face-on, using the e and ω values from the other two panels.

In the text

	Fig. 8 Same as Fig. 1, but for simulation 2. The white arrows denote the number of spiral arms m determined by eye. For the right panel, the dotted arrow denotes an alternative interpretation where the m = 1 spiral is an m = 2 spiral, but one arm is much fainter.
In the text

	Fig. 9 Same as Fig. 2, but for simulation 2, where we model a BHL accretion scenario.
In the text

	Fig. 10 Moment 1 residuals from Keplerian motion of the 12CO and 13CO isotopologues for simulation 2. Top: same as the top row of Fig. 4. Bottom: same as the top row of Fig. 5. The arrows denote the number of spiral arms m determined by eye.
In the text

	Fig. 11 Same as Fig. 7, but for simulation 2.
In the text

Fig. 12 Perturbation depth in units of the pressure scale height as a function of time for all four simulations. The left panel shows the depth for the inner disk (r = 25 AU), and the right panel shows the depth for the outer disk (r = 70 AU). The dashed lines represent simulations with regular disk mass, solid lines ones with low disk mass. The blue lines describe cloudlet capture simulations, whereas the red lines describe accretion in a turbulent medium.

In the text

	Fig. 13 Normalized Fourier-transformed intensity of the m = 1 (left panels) and m = 2 (right panels) modes of the polarized scattered-light intensity of all four simulations. The top panels show the intensity in the inner disk (r = 25 AU), whereas the bottom panels show the intensity in the outer disk (r = 70 AU).
In the text

	Fig. A.1 Same as Fig. 6, but for simulation 2.
In the text

	Fig. A.2 Residual motion found in CO isotopologue line emission for simulation 2 at an inclination of i = 30◦. Left: Same as bottom row of Fig. 4. Right: Same as bottom row of Fig. 5.
In the text

	Fig. B.1 Synthetic scattered light observations of structures found in the simulations with a lower disk mass. Left: Same as zoomed-in area of panels 2 and 3 in Fig. 1, but for simulation 3. Right: Same as zoomed-in areas in Fig. 8, but for simulation 4.
In the text

	Fig. B.2 Same as Fig. B.1, but with an artificially high dust-to-gas ratio of 0.1, reproducing the dust mass used for simulations 1 and 2.
In the text

	Fig. B.3 Same as the left panel of Fig. 7, but for simulation 3 (left) and simulation 4 (right).
In the text

	Fig. B.4 Same as Fig. 4 (left) and Fig. 5 (right), but for simulation 3.
In the text

	Fig. B.5 Same as Fig. 4 (left) and Fig. 5 (right), but for simulation 4.
In the text

	Fig. B.6 Same as Fig. 6, but for simulation 3 (top) and simulation 4 (bottom).
In the text

	Fig. C.1 Left: Same as zoomed-in areas of Fig. 1. Right: Same as zoomed-in areas of Fig. 8.
In the text

	Fig. D.1 Warp properties as a function of radial distance and time, where a darker color denotes a later snapshot. Left: Disk tilt in the reference frame described by the total angular momentum, β′. Middle: Disk twist in the same reference frame, γ′. Right: Warp strength ψ′ in that reference frame. All values are for simulation 1.
In the text

	Fig. D.2 Same as Fig. D.1, but for simulation 2.
In the text

	Fig. D.3 Same as Fig. D.1, but for simulation 3. We note that the range of values on the ordinate is increased compared to that Figure.
In the text

	Fig. D.4 Same as Fig. D.1, but for simulation 4. We note that the range of values on the ordinate is increased compared to that Figure.
In the text