© The Authors 2025

1 Introduction

Planet formation theory has evolved greatly in recent years, with more complete models now connecting the large range of size scales of dust grain growth, the formation of planetary cores, migration, and gas accretion (for a review, see Drążkowska et al. 2023). Despite the steady advancement, protoplanetary disks, the birthplaces of planets, are still largely treated as being isolated from their environment. After their initial formation, the dust budget is treated as fixed, being converted into planetary cores through internal processes. Issues remain in this framework; in particular, recent observations indicate that the dust mass budget of fully formed protoplanetary disks is insufficient to form the giant planets we observe (Manara et al. 2018).

The mass budget problem has sparked an increased interest in treating the possible onset of formation in young, embedded disks still subject to infall (Drążkowska & Dullemond 2018; Morbidelli et al. 2022; Hühn et al. 2025; Carrera et al. 2025). However, infall may also play a major role in the treatment of older, fully formed disks (Ginski et al. 2021; Huang et al. 2021; Garufi et al. 2022; Pineda et al. 2023; Gupta et al. 2024), to the point where an entirely new, second-generation disk may be formed (Dullemond et al. 2019; Kuffmeier et al. 2020). The assumption of isolated evolution is especially challenged by the increasing number of disk observations, both in scattered light and molecular line emission, indicating direct interactions with their large-scale environment. A recent census of about half of the observable Class II disks in the Taurus star-forming region with VLT/SPHERE reveals hints of environment interactions for about a third of the sample (Garufi et al. 2024). They mostly present in the form of streamers, which are arcs of gas connecting the protoplanetary disk scales with its more diffuse surroundings. Furthermore, some targets show various spiral structures that may be related to these streamers (Calcino et al. 2025). Recent studies of such structures in the AB Aur system reveal that, while dynamical modeling of gravitational instability (Speedie et al. 2024) can explain the observed structure, specific spirals may actually be connected to the larger-scale environment, identifying them as streamers of infalling material. The detection of SO line emission, believed to be a shock tracer, further supports this hypothesis (Speedie et al. 2025).

In the past, studies have been conducted to understand the nature of environment-disk interactions. In particular, the structures around the HL Tau and DG Tau disks have been confirmed as inflowing material associated with an accretion shock (Garufi et al. 2022), and the dynamics of the flow onto DG Tau has been explained by models of cloudlet accretion (Hanawa et al. 2024). Similarly, late infall has been confirmed in SU Aur (Ginski et al. 2021), which also shows a self-shadowing structure that may be related to a misaligned inner disk; a peculiar configuration that may be caused by the observed infall. Recently, there have been an increasing number of detections of misaligned and warped inner disks. Villenave et al. (2024) showed that 75% of the targeted edge-on disks show asymmetries that are likely caused by a tilted inner region. Because streamers have been observed in such systems, infall of material with misaligned angular momentum is a prime candidate to explain their origin (Kuffmeier et al. 2021; Dullemond et al. 2022). In turn, this would indicate that the continuous infall of material, even at later evolutionary stages, is a ubiquitous phenomenon.

This hardening observational evidence of evolved disks interacting with and accreting mass from their environment raises the important question of the origin and nature of these structures. Streamers may arise as a connection from clearly defined mass reservoirs at larger scales down to disk scales (Kuffmeier et al. 2023), modeled on smaller scales in the form of “cloudlets” (Dullemond et al. 2019; Kuffmeier et al. 2021), or from interactions in a system with a stellar companion or a flyby (Dai et al. 2015; Kurtovic et al. 2018; Pfalzner et al. 2024). However, they are detected around a wide range of disk-hosting stars, and the possible rarity of encounters with cloudlets of considerable density or other stars might disfavor such interactions as the underlying reason. Alternatively, they may arise solely as the result of interactions with the diffuse interstellar medium (ISM). Here, material is continuously supplied via Bondi-Hoyle-Lyttleton (BHL) accretion (Krumholz et al. 2006; Padoan et al. 2025; Pelkonen et al. 2025), steadily supplying mass and angular momentum that naturally creates streaming structures without the need for a distinct source. This simple process can reproduce various observed disk properties (Winter et al. 2024), suggesting that it may offer an explanation for the high abundance of observed inflow structures.

Understanding the mechanisms that can create streamer-like structures and shape their morphology, as well as identifying their potential origin, is vital for our understanding of the dynamical and mass-budget evolution of protoplanetary disks, with potential far-reaching consequences for planet formation as a whole. In this work, we investigate two different types of interactions of fully formed protoplanetary disks with their environment using hydrodynamical simulations. First, we probed how a single, well-characterized interaction with a cloudlet of gas can create various streamer-like observational signatures (Section 3.1). We then considered a disk that is embedded in the turbulent interstellar medium to examine structures that arise in the absence of a distinct source of infall (Section 3.2). While signatures of single interactions are not long-lived, those arising from constant interactions with the ISM could persist over longer time scales and would therefore be much more abundantly observable.

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 utilizing the advection algorithm by Masset (2000). The radial extent of the grid was from r_min = 10 AU to r_max = 5000 AU, except the simulation presented in Section 3.1, where r_min = 5 AU. We moved the outer edge further outward to r_max = 50 000 AU for simulations with a turbulent initial condition, described in later sections. In the polar direction, the range was from θ_min = 10° to θ_max = 170°. In the radial and polar direction, we employed an outflow boundary condition, whereas the azimuthal boundary is periodic. There are N_r = 250 cells in the radial, N_θ = 100 cells in the polar and N_φ = 225 cells in the azimuthal direction, resulting in H/∆r = H/(r∆θ) = H/(r∆φ) ≈ 1.3 at r = 5.2AU (≈3 at 100 AU). Here, H is the gas pressure scale height and ∆r, ∆θ, and ∆φ are the resolution in the respective direction. We utilized an isothermal equation of state, P = c²_sρ, where P is the gas pressure, c_s is the isothermal sound speed, and ρ is the gas density. A 1 M_⊙ star was placed at the center of the grid. Additionally, we added artificial viscosity to the simulation, corresponding to an α parameter (Shakura & Sunyaev 1973) of α = 10⁻³. Due to the lack of self-gravity in our simulation, we deliberately omitted indirect terms rooted in the stellar acceleration caused by the gas in the potential.

2.1 Protoplanetary disk setup

We introduced a protoplanetary disk around the star at the center of our simulation, so that the midplane coincides with the equatorial plane. We set the density of the corresponding grid cells to $$\begin{array}{rl}\rho (r,\theta )=& \frac{{\mathrm{\Sigma }}_0}{\sqrt{2\pi }{h}_0{R}_0}{\left(\frac{r_{\mathrm{c}\mathrm{y}\mathrm{l}}(r,\theta )}{R_0}\right)}^{-a-b-1}\\ & \times \exp (-\frac{z(r,\theta {)}^2}{2h_0^2{R}_0^2}{\left(\frac{r_{\mathrm{c}\mathrm{y}\mathrm{l}}(r,\theta )}{R_0}\right)}^{-2b-2})\\ & \times {(1+\exp \left(\frac{r_{\mathrm{c}\mathrm{y}\mathrm{l}}(r,\theta )-R_c}{0.05R_c}\right))}^{-1},\end{array}$$(1)

where a is the surface density slope, b is the flaring index, r_cyl and z are the cylindrical coordinates associated with position (r,θ) in the spherical coordinate system, R_c is the disk cutoff radius, R₀ = 5.2 AU is a reference radius, h₀ is the aspect ratio at that reference radius, and Σ₀ is the surface density at that radius. For the disk itself, we chose Σ₀ = 310.547 g cm⁻², a = 1.5, and R_c = 100 AU, resulting in a disk with a mass of M_d = 0.05 M_⊙. In the interest of numerical stability, we introduced a floor value for the density, ρ_floor = 10⁻²²g cm⁻².

We defined the temperature such that the aspect ratio h = H/r, where H is the pressure scale height, is described by $$h(r,\theta )=h_0{\left(\frac{r_{\mathrm{c}\mathrm{y}\mathrm{l}}(r,\theta )}{R_0}\right)}^b.$$(2)

To represent the case of passive irradiation by a star with solar luminosity, we used h₀ = 0.03799 and b = 0.25. Furthermore, we included a floor temperature of T_floor = 10 K, resulting in the isothermal sound speed being given by $$c_s(r,\theta )={({(h(r,\theta )r_{\mathrm{c}\mathrm{y}\mathrm{l}}(r,\theta )\mathrm{\Omega })}^8+{\left(\frac{k_B{T}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{o}\mathrm{r}}}{\mu }\right)}^4)}^{1/8},$$(3)

where k_B is the Boltzmann constant, Ω is the Keplerian orbital frequency, and μ = 2.3m_p is the mean molecular weight of the gas, with m_p the proton mass.

To take the effects of the gas pressure gradient and the exponential cutoff of the disk into account, we initialized our simulations with an azimuthal velocity that deviates from the Keplerian value, $$\begin{array}{rl}{v}_{\varphi }(r,\theta {)}^2=& {\mathrm{\Omega }}^2{r}_{\mathrm{c}\mathrm{y}\mathrm{l}}^2+c_s^2\\ & \times ((b+1){\left(\frac{z(r,\theta )}{r_{\mathrm{c}\mathrm{y}\mathrm{l}}(r,\theta )h}\right)}^2\\ & -\frac{r_{\mathrm{c}\mathrm{y}\mathrm{l}}(r,\theta )}{0.05R_c(1+\exp (-\frac{r_{\mathrm{c}\mathrm{y}\mathrm{l}}(r,\theta )-R_c}{0.05R_c}))}-a+b-2),\end{array}$$(4)

which can become imaginary, in which case we set v_φ = Ωr_cyl. We set the other two velocity components to zero, v_r = v_θ = 0.

2.2 Setup of the environment

We simulated the late accretion of gas onto Class II disks in two different scenarios. First, we considered a controlled case where the environmental interactions are represented by a single, spherical cloudlet of gas encountering the disk on a hyperbolic orbit. Subsequently, we considered the more realistic case of a disk with a systemic velocity embedded in a turbulent interstellar medium (ISM), utilizing a compressively turbulent initial condition with a given power spectrum and dispersion.

2.2.1 Spherical cloudlet

The center of the initially spherical cloudlet was placed at a position described by the initial distance d₀ and the impact parameter b, as well as the speed at infinity v∞ that sets the critical impact parameter $b_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}=\sqrt{G{M}_{\star }/v_{\mathrm{\infty }}^2}$. The initial position (x₀, y₀) is then described by two values in Cartesian coordinates, $$\nu =\arccos \left(\frac{\frac{b^2}{b_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}-d_0}{d_0\sqrt{1+{\left(\frac{b}{b_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}\right)}^2}}\right),$$(5) $$x_0=d_0\cos (\nu ),$$(6) $$y_0=d_0\sin (\nu ),$$(7)

with ν the true anomaly. The initial velocity is given by $$\begin{gathered} v_{x,0}=v_{\mathrm{\infty }}\sqrt{\frac{2b_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}{d_0}+1} \\ \times \cos (\arctan \left(\frac{\sin (\nu )}{\cos (\nu )+{(1+{\left(\frac{b}{b_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}\right)}^2)}^{-1/2}}\right)-\nu +\frac{\pi }{2}), \end{gathered}$$(8) $$\begin{gathered} v_{y,0}=v_{\mathrm{\infty }}\sqrt{\frac{2b_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}{d_0}+1} \\ \times \sin (\arctan \left(\frac{\sin (\nu )}{\cos (\nu )+{(1+{\left(\frac{b}{b_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}\right)}^2)}^{-1/2}}\right)-\nu +\frac{\pi }{2}), \end{gathered}$$(9)

which was applied to all grid cells that belong to the cloudlet, thereby neglecting the spatial extent of it for the calculation of the orbital velocity. We finally rotated the cloudlet’s position and velocity vectors by 30 degrees around a hypothetical x-axis to achieve an out-of-plane infall initial condition.

The cloudlet was initialized with v∞ = 0.5kms⁻¹, d₀ = 3000 AU, a radius of R_cloud = 500 AU, and a mass of M_cloud = 5 × 10⁻³ M_⊙ = 0.1 M_d, resulting in an initial density of ρ_cloud,0 = 5.67 × 10⁻¹⁸ g cm⁻³. We chose b = 0.5b_crit for the impact parameter. The mass ratio was chosen such that the disk disruption due to the cloudlet is moderate. We selected the other parameter values with the goal to keep the radial domain of the grid and the physical time until the interaction small, and ensure b < b_crit across the entire physical extent of the cloudlet. The mass of the cloudlet is higher than suggested by observational data from giant molecular clouds (Falgarone et al. 2004; Klessen & Hennebelle 2010) for our choice of initial radius, but the lack of pressure support from a warm ISM leads to expansion of the cloudlet before encountering the disk. This initial condition was not intended to directly represent realistic physical conditions. Rather, we started with a more compact cloudlet to keep the inaccuracy incurred by neglecting the physical extent when initializing the orbital velocity small. When initially encountering the disk, the mean density of the cloudlet gas is ~10⁻¹⁹ gcm⁻³, which is comparable to the density of filamentary structures in the Taurus star forming region (Pineda et al. 2010; Palmeirim et al. 2013). This initial expansion does not affect the angular momentum of the cloudlet relative to the disk. Therefore, we do not expect it to substantially affect our results. The total runtime of the cloudlet simulation is 42.7 kyr, after which any arising structures have dissipated.

2.2.2 Turbulent interstellar medium

The turbulent ISM simulations were run without a cloudlet, but with a global density ρ_ISM and turbulent initial velocity field, modeled as a Gaussian random field (Dubinski et al. 1995; Krumholz et al. 2006; Seifried et al. 2013). The density was determined via the systemic infall rate, Ṁ_sys, which is a free parameter, resulting in $${\rho }_{\mathrm{I}\mathrm{S}\mathrm{M}}=\frac{{\dot{M}}_{\mathrm{s}\mathrm{y}\mathrm{s}}}{\pi R_{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{i}}^2{v}_{\mathrm{s}\mathrm{y}\mathrm{s}}},$$(10)

where R_bondi = 2GM_*/c²_s is the Bondi accretion radius of the star, M_* is the stellar mass, c_s is the isothermal sound speed of the ISM, and v_sys is the absolute value of the systemic velocity of the star-disk system as it moves through the ISM. Unless specified otherwise, we adopted Ṁ_sys = 10⁻⁸M_⊙ yr⁻¹, so that ρ_ISM = 1.5 × 10⁻²¹ g cm⁻³ for v_sys = 1 kms⁻¹.

The velocity field was realized by randomly sampling N_wav = 10000 sets of the wave number k_n = |k_n|, the wave propagation unit vector k̃_n, the components of the corresponding velocity u_k, and a phase φ_n. We sampled k logarithmically in the interval [k_min, k_max), k̃ on a spherical surface with uniformly chosen φ_k ∈ [0,2π) and μ_k = cos θ_k ∈ [−1,1), and φ_n ∈ [0,2π). We chose a Gaussian random field with a power spectrum P(k) ~ k⁻⁴ as the initial condition, achieved by sampling the components of u_k from a Gaussian distribution with $${\sigma }_k=\sqrt{\frac{k_{\mathrm{m}\mathrm{i}\mathrm{n}}{k}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{N_{\mathrm{w}\mathrm{a}\mathrm{v}}(k_{\mathrm{m}\mathrm{a}\mathrm{x}}-k_{\mathrm{m}\mathrm{i}\mathrm{n}})}\log \left(\frac{k_{\mathrm{m}\mathrm{a}\mathrm{x}}}{k_{\mathrm{m}\mathrm{i}\mathrm{n}}}\right)}{\sigma }_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}{k}^{-0.5},$$(11)

where σ_turb is the standard deviation of the total velocity field in real space, which is a model parameter. The turbulent velocity field at location r is then given by $$v_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}(r)=\sum _{n=0}^{N_{\mathrm{w}\mathrm{a}\mathrm{v}}}{v}_k\cos (k_n{\hat{k}}_n\cdot r+{\varphi }_n).$$(12)

We added systemic infall, caused in reality by a systemic velocity, u_sys, to mimic the movement of the system through the ISM, so that the total gas velocity is given by $$v(r)=v_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}(r)+v_{\mathrm{s}\mathrm{y}\mathrm{s}}{\hat{v}}_{\mathrm{s}\mathrm{y}\mathrm{s}}(r),$$(13)

where we chose ${\hat{v}}_{\mathrm{s}\mathrm{y}\mathrm{s}}=(0,\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$ as the unit vector for the direction of systemic infall to consider one realization of the expected random alignment between the disk and systemic velocity. The velocity initial condition is therefore described by the two parameters σ_turb and v_sys, with the two parameters k_min and k_max describing the limits of the turbulence power spectrum. We note that the turbulence was only introduced as an initial condition, and we did not apply turbulent forcing to drive the turbulence dynamically. As a result, the simulation is no longer self-consistent after the gas from the outer grid edge has reached the disk. At the highest systemic velocity we considered, v_sys = 1kms⁻¹, this occurs after 237kyr, which is much larger than the employed simulation runtime of 30kyr. Nevertheless, the lack of turbulent forcing results in a decrease of the turbulent velocity dispersion over time. In terms of the turbulence Mach number M, a turbulent velocity diffusion of σ_k = {1,0.5,0.25,0.1} kms⁻¹ corresponds to M = {8,4,3.2,2.9} at the beginning of the simulation, which has decayed down to M = {6,3,2.9,2.8} in the end.

2.3 Synthetic observations

Using the Monte Carlo radiative transfer code RADMC3D (Dullemond et al. 2012), we created synthetic observations of select simulation snapshots. We produced observations of the CO J 2-1 transition line, assuming a constant CO to H₂ mass ratio of 10⁻⁴, using 250 points in frequency space covering a spectral range from −5 to +5km s⁻¹ to create maps of the first moment and the peak brightness temperature, T_b,peak. To calculate T_b,peak, we used the Rayleigh-Jeans approximation. We did not employ any prescriptions for the freeze-out or photo-dissociation of CO. Instead of the isothermal temperature distribution employed in the simulations, we first re-computed the temperature with a separate Monte Carlo simulation with N_phot = 10⁸ photons, where we considered the central star to have a blackbody spectrum of a Solar mass star with Solar luminosity. The wavelength grid consists of 100 wavelengths between λ_min = 0.09 μm and λ_max = 3000 μm.

We also created synthetic scattered-light observations. Here, the required dust opacities were calculated using the optool code (Dominik et al. 2021), where we assumed the DIANA standard dust model composition (Woitke et al. 2016). The grains are composed of 70% Mg pyroxene and carbon, which are mixed with a mass ratio of 87 and 13%, respectively. We assumed a porosity of 25% and mix 15 different grain sizes between 1 μm and 3 μm, sampling an n(a) ∝ a^2.5 power-law size distribution. Despite only including such small grains, we assumed a dust-to-gas volume density ratio of 10⁻², assuming that grains in the infall structures cannot grow until reaching the disk. This assumption breaks down inside the disk itself, which we neglect as we focus on the inflow structures themselves. We considered full polarization-dependent scattering, though we flattened the peak of the scattering phase function to avoid numerical issues by “chopping” at 2°. Using this technique, photons scattered at an angle below this value relative to the forward scattering direction are instead treated as not having interacted at all (Dominik et al. 2021). For the spatial distribution of the dust, we assumed that it perfectly follows that of the gas, with a dust-to-gas ratio of 1%.

The scattering Monte Carlo simulation was performed by tracing N_phot,scat = 10⁹ photons. The image was taken at a wavelength of λ_img = 1.245 μm, which we used to calculate the absolute value of the Qφ polarization, $$Q_{\varphi }=|Q\cos (2\varphi )+U\sin (2\varphi )|,$$(14)

where φ is the azimuthal position on the image plane and Q and U are the polarization Stokes vector components.

3 Results

In this work, we investigated the emergence of streamer-like structures in two different scenarios. First, we considered the case of a single gas cloudlet being accreted onto the disk (see Sect. 2.2.1). The primary advantage of this scenario is that it is a controlled environment with a well-defined encounter, however, we lack realism. Second, we considered the more realistic case of a star-disk system moving through the dense, turbulent ISM (see Sect. 2.2.2), which inherently contains more randomness.

3.1 Cloudlet capture

In the cloudlet capture scenario, we find that the arising structures occur for a fundamentally different reason than in previously considered model setups (e.g., Dullemond et al. 2019; Hanawa et al. 2024). With passive stellar irradiation being the only heating source, the background ISM surrounding the initially spherical cloudlet is cold, so that it lacks pressure support and expands freely. If, on the contrary, the ISM was modeled as warm as in previous models, it would act to confine the cloudlet. Differently from our scenario, the arced structure of a streamer would be the natural outcome of a spherical cloud being deformed as the gas elements move along their orbit.

As a result of this difference, the cloudlet has already expanded considerably at the time when it encounters the disk. This is apparent in the second panel of Fig. 1, showing the time evolution of the surface density in and around the disk, as well as the flow lines of the gas. Indeed, in the given field of view, the original spherical geometry is fully lost, and it can be seen that diffuse gas engulfs the disk fully during the encounter. As a direct consequence, over-shooting mass is accumulated on the side of the disk opposite to the infall direction. This morphology is reminiscent of an accretion tail that is the result of converging flow lines caused by gravitational focussing. A key difference to such an accretion tail lies in the angular momentum the accumulated material has with respect to the star. The hyperbolic orbit of the cloudlet has an impact parameter of b = 0.5b_crit, and angular momentum is not lost in the expansion process; instead, it is reflected by the asymmetry of the density and velocity at the time of the encounter.

The gas that is accumulated in the accretion tail remains bound to the star, because b < b_crit. Due to the angular momentum it has with respect to the disk, it falls back onto it in an arced shape. This “fallback streamer” is supported in its shape by the flow of material that is still inflowing from the original direction, creating a shock at the edge of the arc. The resulting streamer morphology can be seen in the third panel of Fig. 1. Due to being caused by the fallback of material that was accumulated in the accretion tail, we find that the direction of the streamer differs from the original inflow direction.

However, the resulting streamer is short-lived. Only ~8.5 kyr later, the surface density of the fallback streamer has decreased in contrast considerably, as shown in the forth panel of Fig. 1. Subsequently, the distribution of mass in the simulation can no longer be recognized as a streamer, although left-over material serves as evidence of the encounter until it is eventually accreted on timescales beyond the runtime of the simulation.

We investigated the fallback shock streamer using a synthetic polarized scattered-light observation, shown in Fig. 2. We find that it produces a bright, extended, and arced signal at a wavelength of λ = 1.245 μm, so that it would indeed be recognized as a streamer in observations, though only the bottom part close to the disk has a flux >0.1 mJy arcsec⁻². Additionally, the inflowing material causes spiral structures to arise at the disk surface, visible in the scattered light. Such a combination of a streamer and spiral structures is found in systems though to be undergoing infall, like SU Aur (Ginski et al. 2021).

Furthermore, we created synthetic observations of the CO J 2-1 transition. The resulting moment 1 and peak brightness temperature maps are shown in Fig. 3. The streamer can be recognized dynamically, adding a recognizable non-Keplerian component to the moment 1 map, similar to systems with observed streamers like HL Tau and DG Tau (Garufi et al. 2022). Closer to the disk, the moment 0 signal of the shock fallback streamer exceeds 600 mJy arcsec⁻² km s⁻¹, and the arc visible in scattered light is also clearly visible in the peak brightness temperature map. In the case of DG Tau, the channel map morphology was explained using a model of a cloudlet confined by a warm medium (Hanawa et al. 2024); here, we find that the accretion of a more extended gas structure, if it has angular momentum relative to the star, can lead to similar results.

Fig. 1 Gas column density of the cloudlet encounter simulation, for six different points in time. The perspective is face-on. The white arrows show the velocity streamlines, computed as a mass-weighted average along the line of sight. The color scale was chosen to highlight low-density features, thereby saturating the protoplanetary disk.

Fig. 2 Polarized scattered light (λ = 1.245 pm) intensity at t = 17.07 kyr for the cloudlet encounter simulation. The inclination of the camera is i = 40°. The white arrow represents the original orbital velocity of the cloudlet, and the white dashed line is the 0.1 mJy arcsec−2 contour line.

Fig. 3 Synthetic CO line emission (J 2-1) at t = 17.07 kyr for the cloudlet encounter simulation. The left panel shows the moment 1 map, with the gray dashed line showing the 600 mJy arcsec−2 km s−1 moment 0 contour. The right panel shows the peak brightness temperature on a square root scale. The gray and white arrows represent the initial orbital velocity of the cloudlet. The white dashed line in the right panel represents the Tb,peak = 20 K contour line.

3.2 Accretion in a turbulent medium

The emergence of fallback shock streamers in the previous section is the result of the encounter with an initially spherical gas cloudlet. Even though it has evolved to a more extended structure by the time it reaches the star, the initial condition is representative of a considerable over-density in the ISM with a symmetric shape. However, it remains an open question how frequent such encounters with clearly discernible cloudlets of gas are, and how well they describe the conditions of the large-scale environment of a star-forming region. Even if they were to serve as a reasonable approximation, we find that the resulting streamer is short-lived. If encounters of this type were the dominant cause for streamers, they should be much rarer than the observational sample suggests (Garufi et al. 2024; Villenave et al. 2024).

In this section, we offer a way to reconcile this discrepancy. In an effort to model the physical conditions of a star-forming region more realistically, we instead modeled the star-disk system as moving through the turbulent interstellar medium with a given systemic velocity (see Sect. 2.2.2). The resulting total velocity field is shown in Appendix A.

In this scenario, streamers can form in a way qualitatively similar to the cloudlet capture case. Rather than creating a reservoir of gas with angular momentum relative to the disk by overshooting cloudlet material, it arises naturally in the form of a BHL accretion tail (Hoyle & Lyttleton 1939; Bondi 1945). In the absence of turbulence, this accretion tail cannot lead to the formation of arc-shaped structures, as no asymmetries that would create angular momentum with respect to the disk would be present. As a result, the material would fall back on a straight trajectory. In the presence of turbulence, however, such asymmetries can naturally be created randomly. Even if the flow structure was unperturbed by the velocity, angular momentum could result from spontaneous enhancements in density. In a simplified picture, one might imagine the large-scale flow lines to vary in direction between two points in time, thereby creating, on average, a qualitatively similar picture as the extended structure caused by the encounter with an expanded cloudlet.

While this idea can help to understand the possibility of fallback streamers arising, the full turbulent picture is more complex. In the following, we discuss simulation results from multiple simulations with various resulting structures. We investigated how the arising structures change with v_sys, σ_turb, and k_max. We considered the influence of M_sys in Appendix B. The numerical choices for these parameters are shown in Table 1.

In our analysis, we find that the morphology of the emerging infall structures is predominantly dependent on the flux of material into the region where the ISM gas has low-enough velocity to be initially bound to the star, and the resulting interaction. Gas originating from farther out regions is too energetic to remain bound to the disk, so that no fallback occurs and no fallback streamers can emerge. Neglecting the spatial details of the turbulent velocity field, we approximated this region as a sphere with radius R_bound, given by $$R_{\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}=\frac{2G{M}_{\star }}{(v_{\mathrm{s}\mathrm{y}\mathrm{s}}+{\sigma }_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}{)}^2}.$$(15)

For all the discussed simulations, we show the absolute value of the angular momentum flux through the surface of this sphere in Fig. 4. Additionally, the streamer morphology is influenced by the direction of the non-systemic, that is, turbulent gas inflow from the ISM. We describe this quantity by the direction of the linear momentum flux through the spherical surface, shown in Fig. 5. The systemic and other symmetrically out- and inflowing contributions to this linear momentum flux are zero. The resulting quantities, although limited to a spatially averaged view, are used in following sections in the discussion of the arising streamers. In the interest of brevity, we only discuss one snapshot for most simulations, at t = 26.07 kyr.

Fig. 4 Absolute value of the angular momentum flux through a spherical shell with radius Rbound as a function of time. The differently colored lines denote different simulations from Table 1. Unless stated otherwise, the shown simulations have vsys = 0.5 km s−1 and kmin = 2π/50 AU. The round markers show the physical time of the snapshot used for the synthetic observations. The cross markers show the physical time of the snapshots used for the time series in Fig. 8.

3.2.1 High systemic velocity: Weakly bound gas accretion

In the first simulation, we considered a scenario where the disk moves through a highly turbulent medium (σ_turb = 1km s⁻¹) with a large systemic velocity (v_sys = 1 km s⁻¹). Here, the ISM gas is highly energetic, so that only gas in a small volume around the disk remains bound to the star, R_bound = 444 AU.

The consequences are two-fold. First, the angular momentum flux through this small volume is the lowest of all simulations and variable on the smallest timescales (see the blue line in Fig. 4). Therefore, any fallback of material does not produce strong arcs. Second, the average linear momentum flux direction also varies strongly and on short timescales. This means that the inflowing gas does not create distinct density enhancements at specific locations, but is smeared out across the volume, preventing the formation of a long-lived and bright BHL accretion tail, as well as streamers that are distinct from this tail.

These considerations serve as the explanation for the majority of features visible in the synthetic observations shown in Fig. 6. In scattered light, sensitive to small changes in the density, the most distinct non-disk feature is a straight extended structure above the disk. In addition, one arc that is spatially separate from the main accretion tail can be seen, originating from the right, but only faintly visible. We find that arced structures arise at different times as well, but, due to the high time variability, we find that they persist only on timescales of a few kyr. The flux of all of these structures is low and likely not observable.

The CO molecular line emission is not as sensitive to these structures. The moment 1 map only shows the upper structure, which is clearly dynamically separated from the background. We find that this particular structure is not a streamer, but an outflow caused by the interaction with the unbound gas encountering a disk. Apart from a small-scale and short-lived minor infall component on the right side, with no corresponding signal in the scattered light, it is the only visible structure with moment 0 signal beyond 190 mJy arcsec⁻² km s⁻¹, as indicated by the contour line in Fig. 6. The peak brightness temperature of the CO emission of the envelope structures is especially faint; the outflow is barely visible and likely not discernible from the background emission observationally. The infall appears more extended here, but at very low intensity, suggesting that the accretion tail cannot accumulate material due to the frequent and strong velocity fluctuations, and likely could not be detected in observations.

We note that we find substructures in the disk in both scattered light and CO synthetic observations. This is also the case for the scenarios described in subsequent sections, and their discussion is the subject of future work.

Fig. 5 Orientation of the non-systemic momentum flux through a spherical shell with radius Rbound as a function of time. The differently colored lines denote different simulations from Table 1. Unless stated otherwise, the shown simulations have vsys = 0.5 km s−1 and kmin = 2π/50 AU. The round markers show the physical time of the snapshot used for the synthetic observations. The cross markers show the physical time of the snapshots used for the time series in Fig. 8. The left panel shows the azimuthal angle θ, and the right panel shows the polar angle φ.

Fig. 6 Synthetic observations of simulation number 1 in Table 1. The left panel shows the polarized scattered light intensity, as in Fig. 2, but at an inclination of i = 30°. The middle and right panels show the CO emission moment 1 and peak brightness temperature maps, as in Fig. 3, but with i = 30°. The peak brightness temperature values are mapped to colors on a square root scale. Here, the gray dashed line is the 190 mJy arcsec−2 km s−1 moment 0 contour. The white dashed line represents the 0.01 mJy arcsec−2 contour line in the left panel, and the 9K contour line in the right panel. The gray and white arrows represent the systemic inflow direction.

3.2.2 Reduced systemic velocity: Multiple arced streamers

For fallback streamers to arise naturally in a simulation with no initial density enhancement, it is essential that the gas encountering the disk remains bound to the star after the encounter, so that the angular momentum flux is visible as an arc. In simulation number 2, we reduced the systemic velocity by a factor of 2, υ_sys = 0.5kms⁻¹, while keeping the ratio between turbulent velocity dispersion and systemic velocity constant, σ_turb/υ_sys = 1. As a result, we find R_bound = 1775 AU.

We expect that, analogously to the scenario in the previous section, significant variations in the infall direction occur, because the systemic velocity and turbulent velocity dispersion are equal. However, unlike previously, the variations should occur over longer timescales due to the larger R_bound. For the same reason, the angular momentum flux should be higher. Indeed, both expectations are fulfilled in the simulation. The angular momentum flux is up to two orders of magnitude higher than previously (see the purple line in Fig. 4). The direction of the linear momentum flux (see Fig. 5) and the angular momentum flux vary comparatively less strongly and on longer timescales.

Here, we find that many spatially distinct infall structures are visible in the scattered light images and the CO emission. This is shown in Fig. 7. In the scattered light image, the strongest intensity originates from a streamer originating from the top right with a flux of ≈0.01 mJy arcsec⁻², while the left structure is less apparent. In the CO moment 1 map, these two streamers are clearly revealed as dynamically separate from each other and the disk. The structures are visible with >225 mJy arcsec⁻² km s⁻¹ across the majority of their spatial extent.

The moment 1 map also reveals that dynamically, the streamers seen in the scattered light consist of multiple sub-streamers: The left one appears to have two components, whereas the right one has three. In the peak brightness temperature map, only the streamer to the right is visible with two components with T_b,peak > 9 K. This brightness temperature is comparable to values observed in the extended structures of GM Tau (Huang et al. 2021). We find that the big streamer also persists for the longest time, ~15kyr. The visible part of left streamer is only weakly arced, likely because it was created earlier in the simulation. In fact, the minimum of the angular momentum flux at ~15 kyr seen in Fig. 4 seems to separate the formation events leading to the two streamers on either side of the disk, as it also corresponds to a shift in linear momentum direction.

As the structures in this setup are highly time dependent due to the strong turbulence, we show a time series of the CO emission in four simulation snapshots in Fig. 8. At first, only weak non-Keplerian features arise around the disk, and the locations of the peak of the CO emission is not correlated with it. In the second panel, the first streamer on the left has formed, and the formation of the one on the right has started. Both big streamers are visible in the third panel, but a third one at the top also appears. It is distinct at this time, but merges with the right streamer by the time of the snapshot shown in Fig. 7. In the end, as seen in the fourth panel, the remaining structures are narrow, fainter, less extended inflows starting to dissipate. In a more realistic scenario with strong turbulence persisting for longer times, continued interactions may continuously lead to the creation of new, strong streamers as old ones dissipate. We therefore expect that fallback streamers could arise abundantly in this type of environment, leading to an increased detection probability compared to the interaction with a single cloudlet.

This strong turbulence scenario is highly sensitive to the maximal wave number, k_max, of the turbulence power spectrum, that is, the presence or absence of small-scale fluctuations. We explore this sensitivity by reducing this value such that k_max = 2π/(1000AU) in simulation number 3. This impacts the arising structures in two major ways, as can be seen in Fig. 9. First, the apparent number of subcomponents decreases, and the spatial separation of the two big streamers increases. The streamers and the ISM appear less chaotic in the CO emission. Both streamers are clearly visible in the moment 0 map (signal >280 mJy arcsec⁻² km s⁻¹) and the peak brightness temperature of the CO molecular emission (>9 K). Combined, it appears that they form a new disk around the existing one, but it can be seen that they originate from two spatially and dynamically distinct mass reservoirs. The formation of these reservoirs is due to two distinct accretion episodes with distinct linear momentum flux. The pink line in Fig. 5 shows a strong change at ~20kyr that separates the two episodes. In addition, the left mass reservoir is connected to the disk in a second, less arced streamer.

We therefore expect bright streamers with spatially well separated origin, like the connection to extended mass reservoirs, to form in dense (accretion rate Ṁ_sys > 10⁻⁸M_⊙ yr⁻¹) environments where small-scale turbulent modes are suppressed or have dissipated due to some physical process.

Fig. 7 Same as Fig. 6, but for simulation number 2, and the level of the moment 0 contour line was adjusted to 225 mJy arcsec−2 km s−1 to enclose the area of the streamer. Additionally, the color range of the brightness temperature map was adjusted.

Fig. 8 Time series of synthetic CO moment 1 maps as in Fig. 7 for the same simulation. The additional black lines are the 7 K contours of the CO peak brightness temperature.

Fig. 9 Same as Fig. 7, but for simulation number 3. The level of the moment 0 contour is 280 mJy arcsec−2 km s−1. The black line shows the 9 K peak brightness temperature contour.

3.2.3 Reduced turbulence: Single arced streamer

For simulation number 4, we reduced the turbulence amplitude relative to the systemic velocity by a factor of 2, that is, σ_turb/v_sys = 0.5. In this scenario, the large-scale turbulent modes are expected to change the direction of inflowing material to a lesser extent. Indeed, the green line in Fig. 5 shows a smaller peak to peak variation of the linear momentum flux direction. It shows a smaller time variability, in part because the turbulent modes with higher k no longer play a significant role. Analogously, the absolute value of the angular momentum flux exhibits less time variability, and is marginally higher in value compared to simulations 2 and 3 at t > 10kyr. The limiting radius for bound accretion is R_bound = 3155 AU.

The emerging infall structures show two predominant features. One is the previously absent straight inflow originating from the BHL accretion tail, which can now be identified more clearly due to less severe turbulent fluctuations. It is supplied by the flow related to the systemic velocity, which carries no angular momentum relative to the star-disk system. Consequently, the resulting feature is not arced and has a substantially extended morphology, so that the appearance of this accretion tail is less reminiscent of a typical streamer in scattered light and the CO moment maps (see Fig. 10). The intensity of the molecular emission is low, so that the accretion tail is very faint with only some areas with T_b,peak > 9 K, and a distinction from the background in the peak brightness temperature map is difficult.

The second feature is an arc of gas, located at the bottom left edge of the disk. The moment 1 map of the CO emission reveals that it is an inflowing structure with a velocity that is different from the straight accretion tail, as shown in the middle panel of Fig. 10, confirming it as a fallback streamer. Compared to the previously discussed case with higher turbulence, only a single, more strongly arced streamer arises in this scenario. Its spatial origin is located close to the straight accretion tail, which is expected due to the lower amplitude of the turbulent waves. The corresponding CO emission is dim, like the straight accretion tail, but the arc is encompassed by the 200 mJy arcsec⁻² km s⁻¹ moment 0 contour. During earlier times in this simulation, an additional arced streamer is visible for a short time interval, from ~20kyr to ~25kyr. At this time, the qualitative inflow structure is similar to the bottom and left streamers found in simulation 2 (see Fig. 7), but the bottom streamer is spatially and dynamically close to the accretion tail as a result of the lower turbulence. As a result, it has dissipated at the time shown in Fig. 10.

Further reduction of the ISM turbulence in simulation 5 further lowers the amplitude of the turbulent velocity fluctuations. To investigate the structures arising in this case, we chose σ_turb/v_sys = 0.2. Here, the time variability of the angular momentum and linear momentum fluxes becomes negligible, as indicated by the orange line in Figs. 4 and 5. Due to R_bound being the highest out of all considered simulations, R_bound = 4930 AU, the absolute value of the angular momentum flux is also maximal. The resulting structures are shown in Fig. 11.

We find that, under these conditions, almost all gas inflowing onto the disk from the environment is part of the straight accretion tail. Scattered light and CO moment maps also show a single, strongly arced fallback-streamer, which is expected from the high angular momentum flux. Its signal is >235 mJy arcsec⁻² km s⁻¹ in moment 0. Due to the low time variability, both the accretion tail and the arced streamer can be seen more clearly in the peak brightness temperature, with T_b,peak > 9K. While no additional streamers emerge at any time during the simulation, this result indicates that even low levels of turbulence can lead to a high angular momentum flux relative to the disk, forming an arced streamer in the process.

4 Discussion

4.1 Fallback streamers should be common

Previous modeling of streamers that only consider them as originating from clear mass filaments, or as the signatures of the accretion of massive gas cloudlets (see Section 3.1), implies that streamers should be rare. This conclusion is strengthened further by the fact that streamer-like structures do not persist for a large fraction of the disk lifetime: only ~1% of it assuming a conservative disk lifetime of 1 Myr.

Unless such encounters are very frequent to replenish the inflow and the related signatures, they fall short of explaining the high abundance of late-infall streamers suggested by observations (Gupta et al. 2023; Garufi et al. 2024). In this case, the fallback streamers we find in the described simulations may offer an explanation. Here, streamers are a common phenomenon for disk-star systems moving through a dense environment. Even though individual streamers do not last considerably longer than in the case of a cloudlet encounter, their emergence does not depend on a specific encounter, but are the result of smaller asymmetries caused by turbulence. As a result, they should arise continuously, increasing the total time they may be observable. In fact, even a low level of turbulence can create a single arced streamer on top of the straight accretion tail.

The fallback streamers, even though their morphology is similar to a stream of material originating from a mass reservoir, do not hold the same information about the large-scale environment. In particular, the infall direction that could be obtained through dynamical modeling is not related to the spatial location of an over-density in the molecular cloud, and a direct connection may not be found. Even in the simplified case of cloudlet accretion, the fallback streamers origin is different from the origin of the cloudlet initially encountering the disk. For turbulent accretion, especially in a highly turbulent environment, the orientation of the streamers is fully random. Here, only the straight infall directly from the accretion tail can be used to determine the systemic velocity of the disk.

Infall of material via BHL accretion after the initial core collapse is also found in parsec-scale star-formation simulations (Padoan et al. 2025; Pelkonen et al. 2025). Previous studies attribute various protoplanetary disk features to this phenomenon (Hennebelle et al. 2017; Kuffmeier et al. 2018, 2023; Winter et al. 2024), further emphasizing its relevance and the need to more in-depth studies of its implications. In this work, we considered the occurrence and morphology of streamers arising in the BHL accretion scenario, and we explored a multitude of environmental conditions. Such considerations were generally not subject of previous large-scale simulations. We focussed on the environment around a single star with a simplified treatment of the environment, identifying fallback as the source for the streamers in our simulations. In the future, our findings should be compared to the morphology and occurrence of accretion structures around disks arising self-consistently in parsec-scale star-formation simulations. The majority of such simulations do not currently reach far into the Class II phase, though.

Fig. 10 Same as Fig. 6, but for simulation number 4. Here, the level of the gray dashed moment 0 contour line is 200 mJy arcsec−2 km s−1.

Fig. 11 Same as Fig. 6, but for simulation number 5, with the level of the moment 0 contour line being 235 mJy arcsec−2 km s−1.

4.2 Streamer morphology as a constraint for local cloud conditions

The quantity and morphology of streamers depends on the local conditions of the ISM as the system is moving through it, especially on the turbulence strength. The observational detection and characterization of streamers onto Class II disks could therefore be used to determine the conditions of their environment. If multiple, arced streamers are detected simultaneously, the disk environment may be very turbulent, whereas the presence of a bright BHL accretion tail suggests lower levels of turbulence. In either case, these structures would indicate a substantial infall rate onto the disk, Ṁ_sys ~ 10⁻⁸ M_⊙ yr⁻¹.

The detection of a straight inflow signature, which can be identified as the accretion tail caused by the systemic flow, serves as strong evidence for a dense environment that may cause additional inflow signatures. The orientation of the straight accretion tail should then be correlated with the proper motion of the disk’s host star determined by, for example, Gaia. In contrast, the orientation of the fallback streamers caused by turbulence should be uncorrelated. Detecting a single streamer would be indicative of a low turbulence environment, but if no corresponding BHL accretion tail is found, a formation scenario related to the cloudlet case may be the favored explanation, implying that the environmental density is low. It may then be caused by such a single encounter, or originate from a mass reservoir. We note that in environments with moderate levels of turbulence, like in our simulation number 4, the presence of an arced streamer may be missed by observations as emission is faint.

The spatial separation and multiplicity of streamers can provide information on the spatial scales the turbulence operates on. For high separation, but low multiplicity that is >1, it might be indicative that small-scale turbulence may be suppressed by some physical process, whose origin may shed more light on the conditions present in the environment of Class II disks. On the other hand, if a high spatial separation is accompanied by a high multiplicity, small-scale turbulence may be operating.

4.3 Implications for disk dynamics and planet formation

The high abundance of observational signatures of streamers, which can be directly related to asymmetric inflow from a dense, large-scale environment engulfing the disk, has strong implications for the dynamical evolution of Class II protoplanetary disks, and thereby for planet formation. Asymmetric inflow has been found to cause a torque that can facilitate the formation of a vortex (Bae et al. 2015; Kuznetsova et al. 2022). This vortex could act as a trap for solids, enhancing the local solid-to-gas ratio to aid in the formation of planetesimals. Furthermore, the inflowing gas itself may locally enhance the turbulence (Lesur et al. 2015; Hennebelle et al. 2017) and related angular momentum transport in the disk, at multiple locations in some cases, which could cause accretion outbursts and allow for the formation of pressure bumps at multiple locations. Gravitational instability may also be triggered as a result of inflow (Kratter & Lodato 2016; Kuffmeier et al. 2018). We will investigate some of these aspects in future work.

In addition to impacting disk dynamics, the late infall that causes these streamers can alleviate the dust-mass-budget problem in Class II disks. Another possible solution for this problem is the early onset of planet formation during the Class 0/I phase, which can strongly profit from infall (Drążkowska & Dullemond 2018; Morbidelli et al. 2022; Hühn et al. 2025), though it remains unclear whether the dynamical conditions in these young disks are favorable for planetesimal formation (Carrera et al. 2025). If dynamically more quiescent Class II disks, however, can also undergo significant accretion over extended time periods, planet formation could start at these late stages and still have enough material to form planets of the observed masses, perhaps even forming planets in multiple generations.

4.4 Model caveats

We modeled the BHL accretion flow onto a protoplanetary disk moving through the ISM in a simplified way. In addition to the described approximation introduced by not driving the turbulence, but just introducing it as an initial condition, we did not model heating and cooling of the gas. Shock heating at the impact locations, where the inflowing gas streamers hit the disk, may be a significant heating source that is also not captured by the Monte-Carlo temperature calculation we perform during post-processing. However, we do not expect a qualitative change to the streamer morphology by a more detailed temperature treatment, but rather an effect on the disk chemistry and the detailed morphology of potentially arising accretion spirals (van Gelder et al. 2021). Other types of modeling may be more suitable for detailed treatment of the impact zone. Furthermore, we did not consider the influence of magnetic fields, which may affect the shape and the angular momentum of the accreted gas, and we neglect stellar feedback to the accretion flow.

When computing the synthetic images, the freeze-out and photodissociation of CO was neglected. The freeze-out of CO onto interstellar dust is particularly relevant for the cold regions far from stellar heating sources. It could reduce the background emission, which would enhance the contrast of the streamer structures and thereby their observability. Conversely, the photodissociation of CO could reduce emission close to the star, balanced by shielding effects of molecular hydrogen and CO itself. To estimate the strength of this effect, we performed order-of-magnitude calculations of the CO column density, assuming photodissociation as the only reaction. We assumed a moderately high unattenuated UV field strength χ = 10⁵, given in units of the Draine field and at a reference radius of 100 AU, used reaction rates from the UMIST database for astrochemistry (Woodall et al. 2007), and employed the self-shielding prescription from Woitke et al. (2009). We neglected extinction by the ISM dust. In regions with a lower density of 10⁻¹⁹ g cm⁻³, 50% of CO at a distance of 1000 AU may dissociate over a streamer lifetime of 10 kyr, but almost all CO inside 500 AU could be dissociated. In regions with a high density of 10⁻¹⁷ g cm⁻³, only 2% of CO is dissociated close to the star at 200 AU. Therefore, photodissociation may play a significant role for low-density streamers arising in high-turbulence regions, but more detailed modeling would be required to accurately assess how astrochemical processes affect the observability of the streamers we find using CO lines.

5 Conclusions

We performed 3D hydrodynamical simulations to investigate the formation of late-infall streamers in dense environments. We modeled late infall as an encounter with a gas cloudlet on an inclined orbit, and we considered the emergence of streamers while a protoplanetary disks moves through the ISM subject to different levels of turbulence. For the latter, a summary of the resulting streamers in a single channel is shown in Appendix C. We draw the following conclusions.

• The encounter of a Class II disk with an expanding cloudlet results in the formation of a “fallback” streamer, where material that initially overshot remains bound and falls back onto the disk with relative angular momentum.

• The apparent spatial origin of the cloudlet streamer differs from the original origin of the cloudlet. Consequently, the continued inflow from the original direction shocks material that is falling back, enhancing the streamer.

• A cloudlet encounter leads to the formation of single, strongly arced streamer and may be relevant in environments with otherwise low density. It dissipates after ~ 10 kyr. Despite its bright CO emission, the short lifetime makes it unlikely to observe cloudlet streamers unless the frequency of such encounters is high.

• Streamers with different morphologies and multiplicities can arise naturally given a high enough accretion rate (~10⁻⁸M_⊙ yr⁻¹), as bound material accumulated during BHL accretion falls back onto the disk with relative angular momentum, similarly to a cloudlet encounter.

• If the disk has a high systemic velocity, υ_sys = 1 km s⁻¹, a majority of accreted material does not remain bound to the disk, so that no extended streamers or even outflows arise.

• For v_sys = 0.5 km s⁻¹ and a turbulent velocity dispersion of σ_turb = 0.5km s⁻¹, accretion is chaotic, but a substantial amount of gas remains bound. Here, multiple arced streamers with subcomponents emerge, and a BHL accretion tail cannot be seen.

• In the absence of small-scale turbulent modes, two streamers connected to distinct mass reservoirs arise, suggesting that these modes may be suppressed in reality if such streamers are observed.

• Reducing the turbulence to σ_turb = 0.25 km s⁻¹ or σ_turb = 0.1 km s⁻¹ allows for the formation of a single arced streamer and a straight BHL tail. Secondary streamers may arise briefly on the timescale of a few kyr.

Our findings suggest that while late infall with a clear origin, represented by cloudlet encounters, can create streamers, they can also arise naturally as an evolved protoplanetary disk moves through the turbulent ISM. This can reconcile ubiquitous detections of late-infall streamers with their short lifetime. Because their morphology is sensitive to the turbulent conditions of the environment and the accretion rate onto the disk, detailed studies of streamers may shed light on the conditions present in the environment of planet-forming disks, which may help to improve our understanding of planet formation.

Acknowledgements

We thank Jiahan Shi, Jorge Pérez González, and Dmitry Semenov for insightful 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 Turbulent velocity field

Fig. A.1 Velocity in Cartesian z direction, vz, at t = 0 of simulations number 1 and 5. The left column panels show a cylindrical slice at φ = 0, and the right column panels show a Cartesian slice at the simulation midplane, θ = π/2. There are different color scales for the two simulations.

We modeled turbulence as a Gaussian random field, whose strength is defined by σ_turb. Figure A.1 shows the realization of this field for one velocity component, Cartesian v_z, for simulations number 1 and 5 from Table 1, which use σ_turb = 1km s⁻¹ and σ_turb = 0.1 km s⁻¹, respectively. The resulting turbulent velocity fluctuations are less significant for simulation number 5, but not negligible. We note that the depicted velocity component contains the contribution from the systemic velocity of the disk.

Appendix B Low accretion rate

Fig. B.1 Same as Fig. 6, but for simulation number 6, with the level of the moment 0 contour line being 60 mJy arcsec−2 km s−1.

The fallback streamers discussed in the main text emerge in a dense environment, that is, where Ṁ_sys = 10⁻⁸ M_⊙ yr⁻¹. However, some stars may undergo considerably lower accretion rates, down to ~10⁻⁹M_⊙ yr⁻¹ (e.g., Herczeg et al. 2023). If less material is accreted onto the disk, related accretion signatures are fainter, and streamers may become undetectable. Therefore, we investigated a low accretion rate of 10⁻⁹ M_⊙ yr⁻¹ in simulation number 6, with a setup that is otherwise identical to simulation number 2, where multiple fallback streamers are visible.

We find that, while the emerging structures are qualitatively similar dynamically, the streamers can indeed no longer be detected, as shown in Fig. B.1. In the peak brightness temperature map, only a single component of the streamer originating from the top right can be seen; all other streamers and their subcomponents can no longer be detected. Its emission is insignificant (T_b,peak < 5 K) compared to the emission in all previous cases and compared to the disk emission, so that this system would not be detected to be moving through an environment with strong turbulence. Almost no scattered light emission can be seen at flux levels >10⁻³ mJy arcsec⁻², and only a small region of the top right streamer can be seen in the moment 0 map with a faint signal of 60 mJy arcsec⁻² km s⁻¹. We therefore conclude that, in order for fallback streamers to be a plausible formation scenario for observed streamers, the infall rate needs to be of the order of ~10⁻⁸ M_⊙ yr⁻¹.

Appendix C Single velocity channels

Fig. C.1 Intensity of the CO emission in a single velocity channel for all simulations listed in Table 1, at the points in time discussed in Section 3.2. The single channel intensity is computed by averaging five consecutive points in frequency space from the RADMC3D simulation, resulting in a channel width of ∆v = 0.2 km s−1. Unless stated otherwise, the shown simulations have vsys = 0.5 km s−1 and kmin = 2π/50 AU.

Figure C.1 shows a summary of the CO emission of the various streamers we find in this work. The emission in a single velocity channel is shown, chosen so that key features of the morphology are visible.

All Tables

All Figures

	Fig. 1 Gas column density of the cloudlet encounter simulation, for six different points in time. The perspective is face-on. The white arrows show the velocity streamlines, computed as a mass-weighted average along the line of sight. The color scale was chosen to highlight low-density features, thereby saturating the protoplanetary disk.
In the text

	Fig. 2 Polarized scattered light (λ = 1.245 pm) intensity at t = 17.07 kyr for the cloudlet encounter simulation. The inclination of the camera is i = 40°. The white arrow represents the original orbital velocity of the cloudlet, and the white dashed line is the 0.1 mJy arcsec−2 contour line.
In the text

Fig. 3 Synthetic CO line emission (J 2-1) at t = 17.07 kyr for the cloudlet encounter simulation. The left panel shows the moment 1 map, with the gray dashed line showing the 600 mJy arcsec−2 km s−1 moment 0 contour. The right panel shows the peak brightness temperature on a square root scale. The gray and white arrows represent the initial orbital velocity of the cloudlet. The white dashed line in the right panel represents the Tb,peak = 20 K contour line.

In the text

Fig. 4 Absolute value of the angular momentum flux through a spherical shell with radius Rbound as a function of time. The differently colored lines denote different simulations from Table 1. Unless stated otherwise, the shown simulations have vsys = 0.5 km s−1 and kmin = 2π/50 AU. The round markers show the physical time of the snapshot used for the synthetic observations. The cross markers show the physical time of the snapshots used for the time series in Fig. 8.

In the text

Fig. 5 Orientation of the non-systemic momentum flux through a spherical shell with radius Rbound as a function of time. The differently colored lines denote different simulations from Table 1. Unless stated otherwise, the shown simulations have vsys = 0.5 km s−1 and kmin = 2π/50 AU. The round markers show the physical time of the snapshot used for the synthetic observations. The cross markers show the physical time of the snapshots used for the time series in Fig. 8. The left panel shows the azimuthal angle θ, and the right panel shows the polar angle φ.

In the text

Fig. 6 Synthetic observations of simulation number 1 in Table 1. The left panel shows the polarized scattered light intensity, as in Fig. 2, but at an inclination of i = 30°. The middle and right panels show the CO emission moment 1 and peak brightness temperature maps, as in Fig. 3, but with i = 30°. The peak brightness temperature values are mapped to colors on a square root scale. Here, the gray dashed line is the 190 mJy arcsec−2 km s−1 moment 0 contour. The white dashed line represents the 0.01 mJy arcsec−2 contour line in the left panel, and the 9K contour line in the right panel. The gray and white arrows represent the systemic inflow direction.

In the text

	Fig. 7 Same as Fig. 6, but for simulation number 2, and the level of the moment 0 contour line was adjusted to 225 mJy arcsec−2 km s−1 to enclose the area of the streamer. Additionally, the color range of the brightness temperature map was adjusted.
In the text

	Fig. 8 Time series of synthetic CO moment 1 maps as in Fig. 7 for the same simulation. The additional black lines are the 7 K contours of the CO peak brightness temperature.
In the text

	Fig. 9 Same as Fig. 7, but for simulation number 3. The level of the moment 0 contour is 280 mJy arcsec−2 km s−1. The black line shows the 9 K peak brightness temperature contour.
In the text

	Fig. 10 Same as Fig. 6, but for simulation number 4. Here, the level of the gray dashed moment 0 contour line is 200 mJy arcsec−2 km s−1.
In the text

	Fig. 11 Same as Fig. 6, but for simulation number 5, with the level of the moment 0 contour line being 235 mJy arcsec−2 km s−1.
In the text

	Fig. A.1 Velocity in Cartesian z direction, vz, at t = 0 of simulations number 1 and 5. The left column panels show a cylindrical slice at φ = 0, and the right column panels show a Cartesian slice at the simulation midplane, θ = π/2. There are different color scales for the two simulations.
In the text

	Fig. B.1 Same as Fig. 6, but for simulation number 6, with the level of the moment 0 contour line being 60 mJy arcsec−2 km s−1.
In the text

Fig. C.1 Intensity of the CO emission in a single velocity channel for all simulations listed in Table 1, at the points in time discussed in Section 3.2. The single channel intensity is computed by averaging five consecutive points in frequency space from the RADMC3D simulation, resulting in a channel width of ∆v = 0.2 km s−1. Unless stated otherwise, the shown simulations have vsys = 0.5 km s−1 and kmin = 2π/50 AU.

In the text