© M. Kuffmeier et al. 2021

1 Introduction

The formation of disks is a natural consequence of the star formation process, and circumstellar disks are important as they are the birthplace of planets (Keppler et al. 2018; Müller et al. 2018; Christiaens et al. 2019; Haffert et al. 2019; Pinte et al. 2018, 2019, 2020; Teague et al. 2018, 2019; Pérez et al. 2020). Observations of Class I disks around young stellar objects (YSOs) demonstrate that disks predominantly form as a by-product of a collapsing prestellar core (e.g., Wolf et al. 2008; Jørgensen et al. 2009; Murillo et al. 2013; Codella et al. 2014; Ohashi et al. 2014; Maury et al. 2019). In using modern telescopes, such as the Atacama Large (sub-)Millimeter Array (ALMA) or the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) instrument at the Very Large Telescope (VLT), it has become possible to obtain more constraints on the properties of disks. For instance, observations show shadows in a significant fraction of so-called transition disks that can be best explained by a misalignment between inner and outer disks (Avenhaus et al. 2014; Marino et al. 2015; Benisty et al. 2017, 2018; Casassus et al. 2018). Some studies show that a (sub)stellar perturber inside (Nixon et al. 2013) or outside the disk (Doğan et al. 2015) can warp and break the disk, hence misaligned systems might be induced by a perturbing object such as a giant planet or a (sub)stellar companion. For instance, several hydrodynamical models show that misaligned systems can form around compact binary systems (Juhász & Facchini 2017; Facchini et al. 2018; Price et al. 2018). Recently, Nealon et al. (2018) as well as Zhu (2019) have also demonstrated that a mildly misaligned planet located in the gap of a disk can lead to a disk breaking up and substantial misalignment in the inner disk with respect to the outer disk. In particular, Gonzalez et al. (2020) and Nealon et al. (2020b) show that a combination of external binary and an inner planet located in the gap between an inner and outer disk is a promising explanation for the specific case of HD 100453.

However, it has also become clear that some disks, especially those around more massive Herbig stars, are associated with tail- or spiral-like structures. One of the most prominent examples of such a structure is AB Aur, which shows the presence of a large disk associated with a tail of a few 1000 AU in size, while the inner disk is misaligned (Nakajima & Golimowski 1995; Grady et al. 1999; Fukagawa et al. 2004). We note that these large-scale structures are located beyond the spiral structures in the inner disk which are interpreted as signs of ongoing planet formation (Boccaletti et al. 2020) or an inner binary companion (Poblete et al. 2020). This is another example of such a misaligned disk around HD 100546, which hosts an extended arm structure around its disks (Grady et al. 2001; Ardila et al. 2007; Walsh et al. 2017). These stars have in common that they are relatively old (a few million years) (van den Ancker et al. 1998; DeWarf et al. 2003), making it unlikely that the associated structures are remnants of the early protostellar collapse phase. Some authors speculate that these structures might be signs of a recent flyby of an external stellar object (Dai et al. 2015; Winter et al. 2018; Cuello et al. 2019). However, Nealon et al. (2020a) find that an external stellar flyby only causes modest and short-lived misalignment at most, hence the authors rule out flyby-induced misalignment.

Instead, the origin of the misalignment in these systems is likely different. Luminosities of some evolved stars are significantly higher than what stellar evolution models predict (Kenyon et al. 1990). Several studies showed that such enhancements in luminosity can be caused by episodic accretion events, where short-term enhancements of orders of magnitude in the accretion rate induce an increase in luminosity due to the scaling of luminosity L ∝ Ṁ (Vorobyov & Basu 2006; Evans et al. 2009; Baraffe et al. 2009; Dunham et al. 2010). In particular, luminosity bursts might be caused by late infall onto the star–disk system (Offner & McKee 2011; Padoan et al. 2014; Jensen & Haugbølle 2018) that lead to disk instabilities triggering accretion bursts (Kuffmeier et al. 2018).

Following up on the possibility of late infall, recent numerical studies by Kuffmeier et al. (2020) (hereafter KGD20) show that disks associated with larger-scale spiral structures can form during encounter events of an existing star with surrounding gas of sufficient angular momentum (see also Vorobyov et al. 2020). The sizes of the resulting disks crucially depend on the impact parameter of the infalling material as it is the major parameter determining the angular momentum of the gas. This process of disk formation is different from the collapse phase leading to the formation of protoplanetary disks commonly found around young protostars. As this process occurs at a later stage, the corresponding disks are referred to as “second-generation disks”.

The question is whether an infall event cannot only lead to the formation of large extended structures (Dullemond et al. 2019), but also cause the formation of a system hosting a misaligned inner and outer disk as found for young forming disks in the simulations by Bate (2018). Following up on KGD20, we investigate the scenario in which a cloudlet of gas approaches an existing star–disk system with different inclination angles. The aim of this paper is to test whether a second-generation disk can form with a different orientation around an existing inner disk.

Section 2 describes the method used in this work. In Sect. 3, we present the results of our numerical models. Based on the results of the models, we compare the results with observations in Sect. 4 and discuss limitations as well as prospects of our study. In Sect. 5, we summarize the results and present the conclusions.

2 Methods

In this study, we analyze the effect of infall onto an already existing primordial disk. We carried out hydrodynamical simulations of a cloudlet of gas that approaches a star–disk system with different infalling angles. In a final step, we post-processed one snapshot of the simulation and produced an image based on synthetic observations of scattered light and thermal emission from dust grains. We use this image to discuss our results in context with observations.

The simulations were carried out with the moving mesh code AREPO¹ (Springel 2010; Pakmor et al. 2016) solving the hydrodynamical equations and assuming the gas to be isothermal with temperature T_gas = 10 K. The target mass for refinement was set to 10²⁴ g (parameters REFINEMENTCRITERION = 1, DEREFINEMENTCRITERION = 1, REFERENCEGASPARTMASS = 10²⁴, and TARGETGASMASSFACTOR = 1 for setting the code units equivalent to cgs units), that is, cells with twice and half the target mass are refined and derefined. The cells were allowed to have a maximum difference in volume of a factor of 5 (MAXVOLUMEDIFF = 5), while the minimum cell volume was set to 10³⁷ cm³ and the maximum cell volume to 4.54 × 10⁵⁰ (MINVOLUME = 10³⁷, MAXVOLUME = 4.54 × 10⁵⁰).

Analogous to the previous study presented in KGD20, we inserted a star with mass M_* = 2.5 M_⊙ modeled as a point mass at the center of our domain and the star was fixed at its location throughout the simulations. We used the accretion recipe introduced in KGD20, which means that the star accreted 90% of the mass located within its accretion radius at every time step. The accretion radius was set to r_acc = 15 au, which is 1.5 the gravitational softening radius of r_grav,soft = 10 AU. The total length of the cubical box is 184 R_cloudlet and we applied periodic boundary conditions. The simulations were stopped early enough so that they were unaffected by the boundary conditions. The cloudlet approached the star with velocity $$v_{\text{cloudlet}}=\left(\begin{array}{c}{v}_{x,\text{cloudlet}}\\ v_{y,\text{cloudlet}}\\ v_{z,\text{cloudlet}}\end{array}\right)=\left(\begin{array}{c}{v}_{\text{i}}={10}^5\text{cm}{\text{s}}^{-1}\\ 0\\ 0\end{array}\right),$$(1)

while the background gas was at rest with respect to the central star. The center of the cloudlet with an initial radius of R_cloudlet = 887 AU was placed at the following location $$r_{\text{init},\text{cloudlet}}=\left(\begin{array}{c}{x}_{\text{cloudlet}}\\ y_{\text{cloudlet}}\\ z_{\text{cloudlet}}\end{array}\right)=\left(\begin{array}{c}-3.22R_{\text{cloudlet}}\\ -b\cos \alpha \\ -b\sin \alpha \end{array}\right),$$(2)

where b = 1774 AU is the impact parameter and α is the infall angle measured with respect to the xy plane of the coordinate system. We note that R_cloudlet = 887 AU and b = 1774 AU correspond to 0.4b_crit and 0.8b_crit, where b_crit is the impact parameter at which a test particle approaching a 2.5 M_⊙ star with velocity v_i would be deflected by 90° (for more details see Sect. 2 in Dullemond et al. 2019).

Adopting a slightly modified version of the description in Klessen & Hennebelle (2010), the cloudlet mass is defined as $$M_{\text{cloudlet}}(R_{\text{cloudlet}})=0.01M_{\odot }{\left(\frac{R_{\text{cloudlet}}}{5000\text{au}}\right)}^{2.3}.$$(3)

The density of the cloudlet ρ_cloudlet is uniform and the density of the background gas is defined as ${\rho }_{\text{bg}}=\frac{1}{800}{\rho }_{\text{cloudlet}}$. To account for turbulence in the cloudlet, we set the velocity of each gas cell according to a Gaussian random distribution using the prescription of Dubinski et al. (1995). The velocities in Fourier space follow a power spectrum of |v_k |² ∝ k⁻⁴, with k being the wavenumber of perturbation. The power index was set in agreement to the velocity dispersion observed in molecular clouds (Larson 1981). The internal velocity dispersion of the cloudlet is 10% of the initial bulk speed by multiplying each velocity by $$A_v=0.1\frac{v_i}{{\sigma }_v},$$(4)

where σ_v is the velocity dispersion of the field before normalization. The setup of the infalling cloudlet is identical to run I1_t_nw used in KGD20, except for the initial displacement of the cloudlet if α≠ 0.

The major difference to the setup in KGD20 is that we inserted a primordial isothermal disk (T = 10 K) with a radius of R_disk = 50 AU around the star at the center following the description presented in Appendix A of Masset & Benítez-Llambay (2016). The disk midplane is located in the xy plane of the coordinate system and has a column density profile of $$\Sigma (r)={\Sigma }_0{\left(\frac{r}{1\text{AU}}\right)}^{-p},$$(5)

where Σ₀ is the amplitude of the column density and p is the power index. Vertically (i.e., in the z direction), the density of the disk drops exponentially, such that the density profile of the disk is defined as $${\rho }_{\text{disk}}(r,z)=\frac{{\Sigma }_0}{\sqrt{2\pi }{H}_0}{\left(\frac{r}{1\text{AU}}\right)}^{-(p+3/2)}{\text{e}}^{-z^2/(2H{(r)}^2)}$$(6)

for radii in the range of r_in = 20 AU ≤ r ≤ R_disk = 50 AU. We note that H(r) is the scale height of the disk given by $$H(r)=\frac{c_{\text{s}}}{\Omega (r)},$$(7)

where the sound speed c_s for an isothermal gas is constant, $$c_{\text{s}}=\sqrt{\frac{k_{\text{B}}T}{\mu m_{\text{H}}}},$$(8)

with a Boltzmann constant k_B, mean molecular weight μ = 2.3, and hydrogen mass m_H. Furthermore, Ω is the orbital frequency of a Keplerian disk given by $$\Omega (r)=\sqrt{\frac{G{M}_{*}}{r^3}},$$(9)

with a gravitational constant G, that is, H(r) scales as H(r) ∝ r^3∕2 for the isothermal disk considered in the simulations. Also, H₀ is the scale height at r = 1 AU. For r< r_in = 20 AU (r_in is twice the gravitational softening length of r_grav,soft = 10 AU), the disk profile was tapered off such that the density dropped according to $${\rho }_{\text{disk},\text{in}}(r,z)={\rho }_{\text{bg}}+\frac{{\rho }_{\text{disk}}(r=r_{\text{in}},z)-{\rho }_{\text{bg}}}{1+{\text{e}}^{-k\left(\frac{r-r_{\text{in}}}{1\text{AU}}+\frac{d_{\text{dec}}}{1\text{AU}}\right)}},$$(10)

where k = 3.5 and d_dec = 2 AU, determining the sharpness of the density drop. For these values of k and d_dec, the density at r = r_in − ϵ = 20 AU − ϵ (where ϵ is an arbitrarily small positive number) is ≈99.9% of ρ_disk,in(r = 20 au). At radius r = r_in − d_dec = 18 AU, the density drops to ≈ 50% of its maximum at r = 20 AU, plus the negligible contribution of the background density, that is, ρ_disk(18 AU, z) = 0.5 × (ρ_disk(20 AU, z) + ρ_bg). To avoid spurious effects from cells in the vicinity of the star, we forced high refinement of the low-density gas located in the inner region within r_in by applying the criterion described in Appendix B of KGD20. Similarly to the treatment for the region inside the disk, the density was tapered off exterior to r > R_disk = 50 AU according to $${\rho }_{\text{disk},\text{ext}}(r,z)={\rho }_{\text{bg}}+\frac{{\rho }_{\text{disk}}(r,z)-{\rho }_{\text{bg}}}{1+{\text{e}}^{k\left(\frac{r-R_{\text{disk}}}{1\text{AU}}-\frac{d_{\text{dec}}}{1\text{AU}}\right)}},$$(11)

such that the density at R_disk + ϵ is ≈ 99.9% of ρ_disk (r = R_disk, z) and at r = R_disk + d_dec = 52 AU, it is ≈ 50% of the density that the disk would have at this radius without the introduction of an outer edge. We used logistic functions for the tapering to ensure that the density profile asymptotically approaches both the values at the inner and outer edge of the disk as well as the minimum density of ρ_bg. To test the effect of the infalling angle of the cloudlet onto the disk, we varied the initial position of the cloudlet angle with the angle α. The tested angles are α = 0°, 35°, 60°, and 90°. The primordial disk is, in all cases, aligned with the xy plane.

In this study, we define Σ₀ = 170 g cm⁻² for the main runs, which is 0.1 times the column density at 1 AU of the minimum mass solar nebula and we performed simulations using p = 1.5. The rotation axis of the inner disk is parallel or antiparallel to the positive z axis of the coordinate system. A parallel and antiparallel configuration corresponds to prograde and retrograde infall. To test the influence of the inner disk mass, we also performed simulations with a very low inner disk mass of Σ₀ =17 g cm⁻² for an infalling angle of α = 90°, and without inclination α = 0° for prograde or retrograde infall. All simulations were carried out for t = 150 kyr of evolution. Moreover, we carried out a comparison run for with α = 90° inclination for a shorter time sequence of t = 100 kyr using Σ₀ = 170 g cm⁻², but p = 1. Figure 1 shows a sketch of the initial setup of our model. The chosen parameters for the altogether eleven runs are summarized in Table 1.

Fig. 1 Sketch of the initial setup. The initial location of the cloud varies depending on α and the chosen impact parameter b. The staris located at the origin of the coordinate system with the cirumstellar disk in the xy plane. The disk rotation has either a right-handed (corresponding to the prograde case) or left-handed orientation (corresponding to the retrograde case). Here, the illustration shows the setup for α = 60°.

3 Results

For illustration purposes, we show the evolutionary sequence on spatial scales of ± 4000 AU, ± 800 AU, and ± 160 AU around the central star for runs 2 and 5 in Figs. 2 and 3. For comparison, the sequences for runs 1, 3, 4, 10, and 11 are shown in Appendix A.

3.1 Outer disk

The figures show that in all cases, the infalling cloudlet leads to the formation of a second-generation disk associated with larger-scale spiral structures for all of our runs. The outer disk extends up to ~ 1000 AU. During the first several 10 kyr of each simulation, the inner disk is already in a stable state that is close to equilibrium, while the outer disk is in its formation phase. In this initial phase of outer disk formation, the second-generation disk evolves on smaller timescales than the inner disk although the orbital and hence dynamical time t_orb of the inner disk is significantly shorter due to the scaling of t_orb ∝ R^3∕2. Clear signatures of the dynamical nature are the arc structures that form during the encounter events and extend to distances beyond 1000 AU in good agreement with the results found in KGD20. In all cases, the outer disk forms with an inclination to the primordial disk equivalent to the initial infall angle α of the individual runs.

To quantify the velocity profile of the outer disk, we computed the azimuthal velocity v_ϕ. The azimuthal velocity was calculated perpendicular to the angular momentum vector of the infalling cloudlet with respect to the central star given by L = M_cloudlet(r_{init, cloudlet} ×v_cloudlet). As expected from the results in KGD20, the velocity structure of the outer second-generation disk is close to Keplerian for cases 1 to 5 (see Fig. 4) to radial distances of several 100 AU. The profiles only marginally differ for the different infall angles α. Comparing the radial extent of the Keplerian profile at different points in time shows that the outer disks become more compact over time and the density drop at the outer edge of the second-generation disk becomes more distinct. However, the azimuthal velocity is still approximately Keplerian up to radial distances of at least 300 AU for all cases at t = 150 kyr. At the inner edge, the azimuthal velocity profile extends further inward at the end of the simulation (t =150 kyr) for the runs with smaller infall angles α. The Keplerian velocity profile continues to small radii in the case of α = 0 and initial parallel orientation of the cloudlet’s angular momentum vector and inner disk (run 1).

Fig. 2 Zoom-in on column densities for run 2 at t = 0, 50, 100, and 150 kyr (from left to right) measured perpendicular to the initial plane of the inner disk with the star at the center. First row: area projection of 8000 AU × 8000 AU. Second row: 1600 AU × 1600 AU. Third row: 320 AU × 320 AU.

Fig. 3 Same as Fig. 2, but for run 5.

Fig. 4 Average azimuthal velocity vϕ profile of the outer disk for run 1 (purple solid line) and run 2 (solid), run 3 (dashed), run 4 (dotted), run 5 (solid) at t = 50 kyr (left panel), t = 100 kyr (middle panel), and t = 150 kyr (right panel). The darker the color of a line, the larger the inclination angle α. We note that vϕ was computed within a vertical extent of |z| < max(6 AU, 0.05r) from the midplane of the disk for 300 logarithmically spaced radial bins in the range from 10 to 104 AU. The gray line is the Keplerian profile for a 2.5 M⊙ star.

Fig. 5 Projection of the inner disk for runs 1 to 5 (top to bottom) at t = 0, 50, 100, and 150 kyr (from left to right) measured perpendicular to the initial plane of the inner disk and azimuthally averaged column density profile at the shown snapshots (rightmost panels). The azimuthally averaged column density in the rightmost panels was computed within a vertical extent of z = ±6 AU and the darkness of the line corresponds to the time of evolution (from light orange representing t = 0 kyr to dark orange at the end of the simulation). We note that the inner part is affected by accretion onto thecentral star and gravitational softening. We therefore exclude the inner region within r = racc + 5 AU = rgrav,soft + 10 au = 20 AU in the azimuthally averaged column density profiles.

3.2 Inner disk

The properties of the inner disk are more affected by the combination of infall angle α and the orientation of the rotation axis of the inner disk. Figure 5 shows face-on projections for runs 1 to 5 at t = 0 kyr, t = 50 kyr, t = 100 kyr, and t = 150 kyr, as well as the azimuthally averaged radial column density profile of the inner disk. After the start of the simulation, the mass located within the accretion radius r = 15 AU accreted onto the central star. As not all of the mass within the accretion radius was immediately accreted, temporary features such as rings and small spirals of low density formed within the accretion radius. However, as the dynamics of the gas within r_acc = 15 AU are stronglyaffected by the numerical parameters of the accretion recipe and gravitational softening, we excludedthis region from the analysis. Accretion also leads to a drop in column density Σ up to a radial distance of ≈20 AU to ≈25 AU from the central star. However, in all cases, the slope of the inner disk remains similar to the initial setup between ≈ 25 AU and ≈ 50 AU during the first t ≈ 50 kyr. For the disks with a lower density (runs 6 to 8), the inner disk starts to align with the outer disk, while the inner disk remains in its original orientation for the most massive inner disk (run 9). We investigate this aspect in more detail in Sect.3.5.

The size of the inner disk remains similar for the case of prograde infall with α = 0, while the disk becomes smaller for the case of retrograde infall with α = 0. The inner disk in the runs with retrograde infall at low α becomes more compressed with a bump in the density profile at r ≈ 25 AU. Due to angular momentum exchange between the outer and inner disk, the inner disk loses angular momentum, and thus the gas moves inward and piles up. Consistently with the transport of angular momentum, more mass accretes onto the central star in the retrograde scenario compared to the prograde scenario (see Sect. 3.3). In contrast, the density profile changes less compared to the initial profile for increasing inclination angle α, as well as for prograde infall. The azimuthal velocity profile of the inner disk shows a similar pattern (see Fig. 6). The velocity profile is Keplerian beyond the initial size of the inner disk for the case of prograde, noninclined infall. For retrograde infall, the radial extent of the Keplerian profile decreases over time for all cases, but the radial extent shrinks more quickly for lower α than for larger α.

Both the density and velocity profiles of the retrograde infall runs 2 to 5 show a clear trend: the higher the infall angle α, the less affected the size is for the inner disk. This trend is expected as the angular momentum of the infalling cloudlet is maximal opposite to the angular momentum of the inner disk in run 2 (retrograde infall without inclination α), and hence brakes the inner disk. An increasing inclination angle α weakens this effect to the point at α = 90° where the angular momentum vectors do not oppose each other in the disk plane. Consistently, the extent of the inner disk is largest at t = 150 kyr in the case of prograde infall, where an outer and inner disk form in the same plane along the same rotational axis.

During the shrinking phase of the inner disk, the transfer of angular momentum also triggers the formation of strong ring features in the inner disk. The ring features are less deep for prograde than for retrograde infall. Again, the analysis shows a correlation with inclination. The larger the inclination angle with respect to the primordial disk plane, the lower the amplitude is of the forming rings. For run 5 (α = 90°), even after t = 150 kyr, the density distribution of the inner disk is still very smooth, while the disk in run 2 shows the strongest ring structures preceding the formation of the smallest disk toward the end of the simulation.

Fig. 6 Similar as Fig. 4, but for the inner disk. We note that vϕ was computed within a vertical extent of |z| < 6 AU from the midplane of the disk for linearly spaced radial bins with a size of Δr = 2 AU.

3.3 Counter- versus co-rotating disks

To provide a better understanding of an infalling event on the inner disk, we analyze the two extreme cases of perfectly prograde and perfectly retrograde infall of runs 1 and 2 in more detail. Generally, run 2 shows the important result that two counter-rotating disks can form as a consequence of late infall events if the infalling material accretes with an opposite angular momentum orientation. Consistently, the infalling cloudlet in run 1 leads to the formation of a second-generation disk with the same rotation axis as the primordial disk. These results are consistent with results from hydrodynamical simulations by other groups, who find the formation of counter-rotating disks in the case of retrograde infall during the early collapse phase (Vorobyov et al. 2015; Bate 2018). Apart from the opposite orientation, the spiral-like features on larger scales as well as the radial extent of the second-generation disk are similar in both scenarios.

On smaller scales, that is on scales of the transition between the outer and inner disk, we can see substantial differences between thetwo extreme scenarios. The transition between the outer and inner disk is smooth in the aligned infall case. In contrast, the outer and inner disk in run 2 are more distinctly separated by a transition zone of ~ 10 AU in width. As shown inFig. 7, the gap is about 30 AU in radial width, while the inner disk has shrunk to about 25 AU in radius for run 2 at t = 150 kyr. The inner disk in run 2 shrinks more than in run 1 during the evolution. In the misaligned cases (runs 3, 4, and 5), the gap is smaller than in the retrograde scenario without inclination, but still more than ≈ 10 AU, while the size of the inner disk shrinks less, as explained above. The features are a consequence of the opposite orientation of the angular momentum vector of the outer disk (see also Wijnen et al. 2016, 2017a). The transition zone in the form of a gap forms where the gas of the forming outer disk interacts with gas of the inner disk (see also Vorobyov et al. 2016). During the interaction, the outer part of the inner disk transfers part of its angular momentum to the inner part of the outer disk. As a consequence, the inner disk shrinks in size, and hence a gap forms in between the inner and outer disk.

Gaps are an important feature of transition disks. In the majority of cases, the gap size is larger than the size of the inner disks. In our models, the deepest and largest gap occurs for retrograde infall, where the gap size is about 30 AU and exceeds the radius of the shrinking inner disk. By comparing simulations of counter-rotating disks with simulations of disks with an embedded giant planet, Vorobyov et al. (2016) show that the resulting gaps are similar. Therefore, infall-induced gaps in counter-rotating disks might be mistakenly interpreted as signs of a planet-bearing disk. Searching for tracers of accretion onto a giant planet can help to reveal the gap origin, such as was done with Hα emission for PDS 70 (Keppler et al. 2018). In the majority of cases in our simulations, the gap region is smaller than the inner disk. We note that the initial size of 50 AU for the inner disk is a somewhat arbitrary choice in our setup and we expect larger ratiosof the gap size to inner disk radius for smaller initial inner disk sizes and/or larger impact parameters b.

Generally, the possibility of angular momentum transport from an inner to outer disk means that material moves radially inward, and hence we expect an increased accretion rate for the counter-rotating case compared to the aligned case as shown theoretically(Lovelace & Chou 1996) and in earlier hydrodynamical simulations (Kuznetsov et al. 1999) for counter-rotating disks. To test whether this is the case in our simulation, we measured the accretion rate as the amount of mass that accretes onto the central sink particle per unit time.

Figure 8 shows Ṁ of run 1 and 2, which are the two extreme cases of prograde and retrograde infall without inclination. In agreement with expectations, the accretion rate increases to higher values in the retrograde infall run after t ≈ 50 kyr, while Ṁ remains smaller in the run with prograde infall. The delay between the visible increase in the accretion rate compared to the first contact of the encountering material with the inner disk can be understood as the time it takes to propagate the angular momentum transfer from the outer part of the inner disk to the vicinity of the star.

The column densities in the outer disk are significantly lower than in the inner disk, as is expected from the significant mass difference between the inner disk and infalling cloudlet mass. For runs 1 to 5, the initial total mass of the infalling cloudlet is less than 25% of the initial inner disk mass. This can be seen particularly well in Fig. 7, which shows the radial profile of column density for runs 1 and 2, that is, the scenario in which the outer disk forms in the same plane as the inner disk. Similar to the inner disk, the column density profile of the outer disk follows a Σ ∝r^−1.5 relation extending from less than 200 AU to about 1000 AU. During the 150 kyr of evolution, the slope remains similar to Σ ∝ r^−1.5, while the overall profile is shifted to increasing column densities. As shown in the velocity profile (Fig. 4), the outer disk becomes more compact and distinct at its outer edge toward the end of the simulation.

Fig. 7 Column density Σ over radius r for prograde infall without an inclination angle in run 1 (green colors) and retrograde infall without an inclination angle in run 2 (orange colors) at t = 50 kyr (dashed lines), t = 100 kyr (dotted lines), and t = 150 kyr (solid lines).The darkness of the line colors indicates the evolution in time from 50 kyr (lightest color) to the final time (darkest color) at t = 150 kyr. The black dashed line corresponds to a theoretical column density profile for Σ0 = 170 g cm−2 and p = 1.5 which extends to 2000 AU.

Fig. 8 Time evolution of the accretion rate for aligned infall run 1 (purple line) and counter rotating infall run 2 (blue line) with infall angle α = 0°.

3.4 Eccentricity

The infall onto the primordial disk is asymmetric. As a consequence, both the outer and inner disk become eccentric as seen in the projection plots shown in Figs. 2, 3, and 5. To quantify this aspect further, the eccentricity e was computed for each grid point according to $$e=\sqrt{1-\frac{2\left|{ϵ}_{\text{soe}}\right|h^2}{{({\text{GM}}_{\text{tot}})}^2}},$$(12)

where ϵ_soe is the specific orbital energy $$ϵ=\frac{E_{\text{pot}}+E_{\text{kin}}}{{\mu }_{\text{red}}}$$(13)

and h is the specific angular momentum $$h=\frac{L}{{\mu }_{\text{red}}},$$(14)

where ${\mu }_{\text{red}}=\frac{M_{*}m}{M_{\text{tot}}}$ is the reduced mass with m being the mass of the individual grid cell and M_tot = M_* + m ≈ M_* = 2.5 M_⊙. In Fig. 9 we show the time evolution of eccentricities in the inner and the outer disk.

In all cases, the forming second-generation outer disk is significantly more eccentric than the initially circular primordial disk. The time evolution of the eccentricities of the outer disk illustrates the sequence of the encounter event. At the beginning of our run, the outer disk has not formed yet, which is why the eccentricity is close to its maximum value at 1. When the cloudlet approaches the central star, part of the gas is accreted by the central star causing a rapid decrease in eccentricity to e ≈ 0.2 from t ≈ 10 kyr to t ≈ 20 kyr. The eccentricity remains at e ≈ 0.2 for about 50 kyr before it rises again up to e ≈ 0.4 at t = 110 kyr. The rise in eccentricity occurs when gas that initially passes by the star at a high speed returns to the scales of the outer disk for a second time. This pattern is similar for all runs, though there are some minor differences between the setups. The runs with intermediate angles, that is 35° and 60°, remain slightly more eccentric between 20 and 70 kyr than the other three runs with either a perpendicular infalling angle or zero inclination angle. After t = 110 kyr, the disk in the runs with intermediate inclination remain eccentric with e ≈ 0.4, while the eccentricity slowly decreases again for the other runs. These differences are minor and not substantial though.

In all cases, the inner disk remains less eccentric than the newly forming second-generation disk. Comparing the eccentricity in the inner disk of the different runs, however, shows differences. The eccentricity rises more rapidly from an initially circular orbit for theruns without any inclination angle. Also, the eccentricity rises more quickly for the run with a low inclination of 35° compared to the runs with higher inclination. Additionally, for the runs with zero or low inclination angles, the maximum eccentricity is highest with e > 0.1. However, after the eccentricity is at its maximum after t ≈ 90 kyr, the inner disk becomes less eccentric again for the low or zero inclination retrograde cases, while it remains above e ≈ 0.1 for the prograde cases. We note that in the case of retrograde infall at a zero inclination angle (run 2), the drastic increase in eccentricity shown at a radial distance between 25 to 30 AU after t = 130 kyr is not associated with the disk anymore, but it mostly traces the low-density gas beyond the shrinking inner disk (see evolution of the disk in Fig. 5). The reduction in eccentricity in the retrograde low inclination runs correlates to the compression of the inner disk. In contrast, e continuously increases for the high inclination cases for both retrograde and prograde infall until the end of the simulation to e ≈ 0.12. Comparing the retrograde and prograde scenarios at identical infall angles, we find that the rise in eccentricity occurs slightly more rapidly for prograde infall than for retrograde infall on average. This is likely also related to the fact that for retrograde infall more mass is accreted onto the central star, while in the prograde scenarios more mass remains in the domain and contributes to the higher eccentricity.

Fig. 9 Time evolution of azimuthally and radially averaged eccentricity in the inner (top panel) and outer disk (bottom panel). For the inner (outer) disk, the average eccentricity between a radial distance of r = 25 (200) AU to 30 (225) AU and a vertical extent of z = ±0.1r over the entire azimuth is shown.

3.5 Gas infall and orientation of inner with respect to outer disk over time

As the cloudlet is not pressure confined, it starts to expand from t = 0 (see Fig. 10 and also the discussion in Dullemond et al. 2019). During the cloudlet encounter with the central star, part of the gas is captured by the central gravitational potential and forms the outer disk. The remaining part is only deflected. The part of the expanding remnant of the cloudlet furthermost from the star continues to expand, escapes the gravitational pull from the central star, and moves away from the star–disk system. The part of the remaining cloudlet that is closer to the central gravitational potential, but not yet located in the second-generation disk, eventually falls toward the star–disk system. Part of the infalling gas at the border of the highly disrupted cloudlet appears as a faint filament oriented almost perpendicular with respect to the disk plane when looking at the system on scales of ~ 1000 AU along the x-axis (see Fig. 11).

As a consequence of the nonzero inclination angle with respect to the disk plane in runs 3, 4, and 5, the outer disk forms with the corresponding inclination to the inner disk, as illustrated in Fig. 11. The outer disk does not show any warp during the time evolution of the eleven runs. As the gas in the primordial inner disk and the newly forming outer disk exchange angular momentum, the orientation of the inner disk can change during the simulation. Therefore, we compared the orientation of the angular momentum vector of the inner part at different times with each other. In particular, we computed the angular momentum vector L₅₀ (t) for all the gas located within a radial distance r of 15 AU < r < 50 AU from the central star. Figure 12 shows the evolution of the angle between the angular momentum vector at each snapshot with respect to the initial orientation of the angular momentum vector L₅₀ (0). The plot shows that the inner part of the gas undergoes mild changes in all cases. Infall with higher inclination angles generally correlates with higher deviation angles from the initial state of the inner disk.

For the main runs, the disks remain close to their initial state by the end of the simulation. Even in the case of 90 degree infall, δ is less than 12° by the end of the simulation at t = 150 kyr. However, the inner disk becomes continuously more tilted for the runs with higher inclination angles, suggesting that they will eventually align with the outer disk. The inner disk aligns with the second-generation outer disk, as the angular momentum of the outer disk exceeds the angular momentum of the inner disk. The comparison run of the infall angle α =90° with ten times less mass in the inner disk (run 8), hence ten times less angular momentum, as well as the run with a more massive inner disk (run 9) confirm this result (see Fig. 13). For the low-density disks, the inner disk starts to align with the outer disk earlier. By the end of the simulation at t = 150 kyr, the inner disk is almost fully aligned with the outer disk. In contrast, the more massive disk in run 9 is still oriented close to its initial orientation at the end of the simulation at t = 100 kyr.

This also suggests that in configurations where the inner disk is massive enough so that it exceeds the angular momentum of the outer disk that is forming, the outer disk would align with the orientation of the inner disk. We did not run models with such configurations, but it remains an important topic of study for the future. Especially, the question of whether the outer disk becomes (temporarily) warped in a setup of continuous infall is intriguing. Comparing the results of such a scenario with existing work on warped disks induced by a perturbing companion (Papaloizou & Terquem 1995; Larwood et al. 1996; Lubow & Ogilvie 2000; Martin et al. 2019) will be particularly relevant given the observations of warped disks around young stars. In our models, the inner disk changes orientation as a whole, but it does not become warped. This is consistent with previous models of small disks with low viscosity (wave-like regime), where the disk bends globally as a consequence of a warp wave that is relatively large compared to the disk (Nixon & Pringle 2010; Nixon & King 2016). In cases of higher viscosity (diffusive regime), warping and even tearing of the disk would happen more easily (see Sect. 4.3).

Fig. 10 Evolution of run 2 seen along the x-axis for 105 kyr showing the column density every 15 kyr. The images show how the cloudlet expands and how gas falls toward the central star forming the disk while the cloudlet dissolves.

Fig. 11 Inner and outer disk of runs 2 to run 5 (from left to right panel) seen from a viewing angle along the x axis at t = 100 kyr (upper panel) and t = 150 kyr (lower panel).

4 Discussion and future prospects

4.1 Limitations of the model

The key question we address in this paper is whether misaligned disks can form as a result of a late infall event. To test this possibility, the model setup was idealized. The central star only acts as a point mass and the infalling material is initially confined in a spherical cloudlet. In our model, the gas is assumed to be isothermal and cold with a temperature of 10 K. In reality, the inner disk has a temperature profile due to heating by the central star. A more realistic treatment of the thermodynamics would change the properties of the inner disk, for example, causing flaring of the disk, but we showed in KGD20 that heatingfrom the central star only marginally affects the formation of the second-generation disk due to its distance to the heating source.

Another simplification of our model is the lack of magnetic fields. Magnetic fields are supposed to play a major role in transferring angular momentum during the protostellar collapse phase, in particular through the launching of outflows from the inner disk. At later times, magnetic fields are likely to be less important for more evolved disks as the field strength dissipates over time. Magnetic fields likely affect the younger second-generation disk by reducing its disk size via magnetic braking. The corresponding transfer of angular momentum would hence lead to stronger shrinking of the inner edge of the outer disk and thereby indirectly affect the size of the inner disk. Nevertheless, given the high angular momentum of an encounter event for sufficient impact parameter b, we do not expect quenching of infall-induced rotational structures when incorporating magnetic fields.

There is one particular aspect regarding the role of magnetic fields that would be interesting to test in future works though. If magnetic braking transfers enough angular momentum from the outer disk, the angular momentum budget of the outer disk could become less than the angular momentum of the inner disk. As a consequence, the inner disk would cause a stronger torque on the outer disk, and the outer disk would align with the inner disk. This is opposite to the results presented here, where the inner disk starts to align with the plane of outer disk due to the excess of angular momentum of the infalling material.

The biggest idealization in our model is certainly the simplified initial condition of a two-phase medium of empty space and an isolated gas cloudlet. As discussed in KGD20, a real giant molecular cloud (GMC) is filamentary (André et al. 2010), and a star likely encounters assemblies of cold gas deviating from spherical symmetry. Nevertheless, observations of luminosity bursts for sources associated with filamentary arms such as FU Orionis or Z CMa (Liu et al. 2016, 2018) strongly suggest the occurrence of late accretion events. Pineda et al. (2020) consistently reported the presence of an impressive 10000 AU gaseous streamer for the binary system Per-emb-2 (IRAS 03292+3039) with the IRAM NOrthern Extended Millimeter Array (NOEMA). The streamer was likely responsible for past accretion bursts as the chemical abundances are inconsistent with the relatively low present-day luminosity.

Arguably, the best candidate for an infall-induced misaligned disk in action is SU Aur. Ginski et al. (2021) recently reported the presence of disk shadows together with extended filamentary arms that are strikingly similar to what is seen in our simulations. Moreover, Huang et al. (2021) show the presence of extended arms around the star–disk system of GM Aur, which the authors interpret as a sign of infalling motions similar to the cases of AB Aur (Nakajima & Golimowski 1995; Grady et al. 1999; Fukagawa et al. 2004) and SU Aur.

Fig. 12 Time evolution of angle δ for the inner (top panel) and outer disk (bottom panel). For the inner disk, δ was measured as the angle of the angular momentum vector at time t with respect to the initial angular momentum vector of the inner disk L50(0). For the outer disk, δ is the angle between the angular momentum vector of gas located within a range from radius r = 200 AU to r = 225 AU and the initial angular momentum vector of the infalling cloudlet.

Fig. 13 Same as the upper panel in Fig. 12, but for runs 5 (dashed line), 8 (dotted line), and 9 (solid line) only, i.e., the runs with infall angle α = 90°, and varying the initial density profile of the inner disk. The initial inner disk mass is Mi,disk ≈ 4Mcloudlet in run 5, Mi,disk ≈ 0.4Mcloudlet in run 8, andMi,disk ≈ 24Mcloudlet in run 9. In all cases, Mcloudlet = 1.87 × 10−4 M⊙.

4.2 Likelihood of late accretion events

The idealized models investigated in this paper demonstrate that second-generation disks can form around an already existing star–disk system of a Herbig star during a late encounter event. Both computational models (Offner & McKee 2011; Padoan et al. 2014; Jensen & Haugbølle 2018) as well as observations of accretion bursts (e.g., Kenyon et al. 1990; Evans et al. 2009) show that stars can undergo accretion events after their initial formation phase. Recently, Alves et al. (2020) reported on a possible candidate for simultaneous star and planet formation. [BHB2007] 1 hosts a Class II disks while also being connected to the larger-scale environment through filamentary arms that strongly indicate ongoing infall of fresh gas and dust. In particular, analyzing magnetohydrodynamical models that successfully reproduce the stellar initial mass function (IMF), Pelkonen et al. (2021) found that ≈ 90% of the mass of intermediate-mass stars stems from gas unbound to the progenitor core. In contrast, low-mass stars tend to accrete relatively more mass during their initial formation phase from the corresponding prestellar core. Their models neglect stellar feedback and stellar radiation likely affects the formation process of higher mass stars (Kuiper et al. 2015, 2016; Rosen & Krumholz 2020). However, a higher probability of late accretion events for intermediate-mass stars than for low-mass stars could explain that tail and arc structures are more commonly observed around Herbig than T Tauri stars (see discussion in Dullemond et al. 2019).

4.3 Observations of misaligned configurations

In this study, the infalling cloudlet interacts with a single star, but misaligned disks, are also observed around binaries. A particularly interesting case is Oph IRS 43, which is a binary star with two circumstellar disks and a surrounding circumbinary disk. The two circumstellar disks are misaligned with respect to each other, and the surrounding circumbinary disk has an inclination of ~ 70° with respect to the binary orbit (Brinch et al. 2016). Another prominent example of strong misalignment (~ 70°) between theorbital plane of the inner binary and outer disk is HD 142527 (Marino et al. 2015). Also, Alves et al. (2019) found a flattened envelope structure that appears to be misaligned with the orbital plane of the close binary [BHB2007] 11. Carrying out star-forming models, Bate (2018) found misalignment of the binary orbital plane and circumbinary disk in one case during an early stage of star formation for a short time period of ~ 30 kyr. Models similar to the ones presented here, but with infall on an existing binary system, would help to understand whether such misaligned configurations can be formed through late infall events.

In recent years, shadows located in the outer disk have been observed for a number of objects, such as HD 142527 (Avenhaus et al. 2014), DoAr 44 (Casassus et al. 2018), HD 135344B (Stolker et al. 2016), HD 100453 (Benisty et al. 2017), or HD 143006 (Benisty et al. 2018). Moreover, Pontoppidan et al. (2020) reported on the variability of the great shadow in the Serpens Molecular Cloud (Pontoppidan & Dullemond 2005), strongly hinting that a protoplanetary disk is the cause of the shadow.

In Fig. 14, we present a multiwavelength composite image for select bands of the observatories VLT (1.66 μm), SOPHIA/HAWC+ (53.0 μm), and ALMA (870 μm). To get this composite image, we post-processed the results of our hydrodynamical runs with the Monte-Carlo (MC) code POLARIS² (Reissl et al. 2016), which can handle the native Voronoi mesh of AREPO. POLARIS self-consistently calculates the dust temperature assuming that the dust component is in equilibrium with the radiation field. For the radiation we utilized the central star and an additional diffuse interstellar radiation field (ISRF) (Reissl et al. 2016, 2020). The star has an effective temperature of T_eff= 9972 K and a radius of R = 1.0 R_⊙ leading to a luminosity of L = 8.2 L_⊙, while the ISRF follows the parametrization of Mathis et al. (1983) with G₀ = 1 which is typical for the Milky Way. For the dust we assumed spherical grains with a composition of silicates and graphite and applied the canonical MRN Milky Way dust with a minimal grain size of a_min = 5 nm and a power-law size distribution n(a) ∝ a^−3.5 (Mathis et al. 1977). However, circumstellar disks are environments that seem to support significant grain growth (Ossenkopf & Henning 1994; Kataoka et al. 2017). Hence, we used an upper grain size of a_max = 10 μm (Kataoka et al. 2017). The dust mass to gas mass ratio was uniformly set to 1% (Bohlin et al. 1978) everywhere in the mesh.

Finally, we calculated intensity maps via a MC approach (Reissl et al. 2016, 2020) for scattered stellar radiation in the near-IR and via a ray-tracing scheme in the submillimeter for the thermal dust emission (see Appendix B in Reissl et al. 2019, for details). Here, we zoomed into the inner 800 AU around the star and assumed a distance to the observer of 140 pc.

The synthetic observations of our models confirm that misalignment between an inner and outer disk can cast shadows on the outer disk, as shown by Marino et al. (2015). Similar to the UV and mid-IR observations presented, for example, in Marino et al. (2015) and Benisty et al. (2017), we note that the shadows of scattered light on the opposite sides of the outer disk are not exactly parallel, but they have a slight offset. We also emphasize that the gaseous stream from the outer disk onto the inner disk is clearly detectable in multiple bands of our synthetic image. In addition, we see that the dust emission in the mid-IR and submillimeter is also decreased by a factor of about three in the shadowed regions, which means that not all of the light is blocked by the inner disk. We expect a larger shadow width in the outer disk for a more realistic treatment of the thermodynamics in the inner disk (see Sect. 4.1) because a nonisothermal inner disk would be flared and more puffed up.

Nevertheless, even for the thin isothermal inner disks in our simulations, we find that the shielding of radiation by the inner disk has a significant impact on the dust temperature distribution of the outer disk. Considering the differences in eccentricities between the inner and outer disk, we also expect that these differences are detectable in CO velocity channel maps. However, a detailed analysis of these effects goes beyond the scope of this paper. We analyze and discuss the observational constraints in more detail in a follow-up study.

For the observed systems, the inner disk is expected to be smaller than assumed in our models, likely of ~ 1 AU to ~ 10 AU in size. The inner disk is tapered off at radii less than ~20 AU to avoid extensive computational costs. First, the hydrodynamical time step scales with radial distance in a Keplerian disk according Δt ∝ r^1.5, which slows down the simulation when resolving the dynamics at small radii. Second, and more relevant in our case, for the assumed density scaling of Σ ∝ r^−p with p = 1.5 or p = 1 evolving smaller radii with the same resolution as used in this study requires modeling many more cells at the high densities in the inner part of the disk, and hence increases the computational costs. However, given that a second-generation disk with similar properties forms in scenarios without (see KGD20) and with the presence of an inner disk (this study), we expect that outer disks can also form around smaller disks than assumed in this work. The results of this study using a fixed cloudlet configuration and a constant impact parameter for all runs predominantly demonstrates that infall with a substantial inclination angle can lead to the formation of misaligned systems. The size of the outer disk can be smaller for lower impact parameters as previously shown in KGD20.

Our results demonstrate that retrograde infall causes the inner disk to shrink over time. Due to conservation of angular momentum, shrinking is more efficient for lower inclination angles, but it still occurs for substantial infall angles. In cases of retrograde infall with a low angle, infall can hence lead to a configuration that is more of a misaligned ring-disk system. Recently, Kraus et al. (2020) reported on such a configuration around GW Orionis (see also Bi et al. 2020). Kraus et al. (2020) attribute the misaligned configuration in GW Orionis to disk tearing (Nixon et al. 2013) induced by the triple star system at the center. The crucial question is, however, whether breaking at a relatively large radius of ≈ 40 AU can actually occur for low viscosities (i.e., low values of the Shakura-Sunyaev parameter α_SS (Shakura & Sunyaev 1973)) that are expected in protoplanetary disks. The formula defined for the bending-wave regime, where α_SS< H∕r suggests a significantly smaller breaking radius than the formula for the diffusive regime that was used to estimate the breaking radius (see Eqs. (10) and (A3) in Nixon et al. 2013, see also Eqs. (2) and (3) in Facchini et al. 2018). Kraus et al. (2020) assumed relatively large values of α_SS ≈ 0.01 to 0.02 for their simulations and applied the formula defined for the diffusive regime, where α_SS exceeds the aspect ratio of scale height to radius (α_SS > H∕r) for their analytical considerations of the disk breaking radius (Eq. (A3) in Nixon et al. 2013, see also Eq. (2) in Facchini et al. 2018). In fact, Smallwood et al. (2021) carried out smoothed particle hydrodynamics (SPH) simulations with varying viscosities in the disk and found no breaking in the bending-wave regime. They suggest instead that the stellar multiple leads to a warped inner disk after a planet carved a gap in the disk.

Our results strongly suggest that configurations of misaligned inner and outer disks can also arise after infall events. Especially, GW Orionis is located in a clustered environment of Lambda Orionis, where the probability of infall or gas-encounter events is expected to be higher than for isolated objects.

Similar to the results found by Wijnen et al. (2017b), the inner disk tends to align with the orientation of the outer disk for misaligned infall with sufficient angular momentum, eventually erasing the initial misalignment between the disks. However, the inner disk can only align with the outer disk if angular momentum can be transferred between the disks. Considering an early onset of planet formation (Harsono et al. 2018), as well as the frequent observation of ring structures in young disks (ALMA Partnership 2015; Andrews et al. 2018), the presence of a deep gap carved by an accreting planet might quench alignment of the disk inside the radius of the planet. We note that as opposed to the scenario of planet-induced misalignment considered by Zhu (2019), the orientation of the inner disk would remain unchanged.

Fig. 14 Multiwavelength RGB composite image of run 3 for typical observational bands of 1.66 μm (blue, VLT H-band), 53.0 μm (green, SOPHIA/HAWC+ A-band), and 870.0 μm (red, ALMA band 7). The disks were observed from a distance of 140 pc under an angle of 45° with respectto the main axis.

4.4 Future goals for modeling

Ultimately, to model the infall event in a more realistic setup in the future requires models that consistently solve the small-scale physical processes on AU scales, while simultaneously accounting for the protostellar environment. Multiscale simulations similar to the ones presented in Kuffmeier et al. (2017, 2019), but for later times of evolution, are therefore the crucial next step. Carrying out simulations for a sample of objects also allows one to obtain statistics on the occurrence rate of encounter events, as well as to constrain the role of infall in possibly warping the disk around binary stars, thus causing misalignment of the binary orbit and disk.

5 Conclusion

Hydrodynamical models with the moving-mesh code AREPO demonstrate that for a star with an already existing disk, infall of a gaseous cloudlet induces the formation of a second-generation disk around the inner primordial disk. In the likely event of new material approaching the inner disk on an inclined trajectory with respect to the midplane of the inner disk, the orientation of the new outer disk is misaligned with the inner disk. As a consequence of the misaligned configuration, the inner disk casts a shadow on the outer disk as shown in synthetic observations, which is consistent with observations of several transition disks. The orientation of the outer disk is governed by the inclination angle of the infalling material, and hence the misalignment between the inner and outer disk can be substantial.

In the case of retrograde infall, the outer disk counter-rotates with respect to the inner disk and it carves a larger gap between the two disks than for prograde infall. The resulting gap can become larger than the size of the inner disk, which is consistent with observations of several transition disks. For low inclination angles, retrograde infall removes angular momentum from the inner disks, thus shrinking the inner disk and enhancing accretion in comparison to prograde infall.

As a by-product of infall, the forming outer disks are associated with characteristic arc and spiral structures on scales of several 100 to 1000 AU. The outer disk forms with high eccentricity of 0.2 to 0.4, mostly independent of the inclination angle of infall. By the end of the simulations at t = 150 kyr, the eccentricity of the outer disk is 0.3 to 0.4 and the disk’s outer edge becomes more sharp, independent of the infall angle. We therefore predict that infall-induced disks have significant eccentricities, which are expected to be detectable in CO channel maps.

In the inner disk, infall triggers modest eccentricity of e = 0.05 to ≈ 0.1 in all cases. Smaller inclination angles of infall cause larger and more variable eccentricity than infall with high inclination.

During the evolution, the inner disk tends to align with the orientation of the outer disk due the higher angular momentum of the infall material on the larger-scale disk. We show that alignment with the outer disk occurs faster for a lower initial mass of the inner disk, consistent with angular momentum conservation.

Furthermore, infall triggers bumps in the density profile of the inner disk, which could possibly support the trapping of dust grains, and hence provide favorable conditions for dust growth and ultimately planet formation. Finally, if misaligned disks are predominantly a consequence of infall, we would expect a higher fraction of misaligned disks and arc-structures around Herbig stars than TTauri stars as recent star formation models indicate a higher probability of late infall for intermediate-mass stars than for low-mass stars.

Acknowledgements

We thank the referee for a thorough review of an earlier version of the manuscript. In particular, we thank the referee for providing constructive feedback and useful suggestions that improved the quality of this work. We thank Sebastian Wolf for useful comments on the manuscript especially on the observational aspects. This research project of M.K. is mainly supported by a research grant of the Independent Research Foundation Denmark (IRFD) (international postdoctoral fellow, project number: 8028-00025B). M.K. acknowledges the support of the DFG Research Unit ‘Transition Disks’ (FOR 2634/1, DU 414/23-1). M.K. also gratefully acknowledges that this project has received funding from the European Union’s Framework Programme for Research and Innovation Horizon 2020 (2014–2020) under the Marie Skłodowska-Curie Grant Agreement No. 897524. F.G. acknowledges support from DFG Schwerpunktprogramm SPP 1992 ‘Diversity of Extrasolar Planets’ grant DU 414/21-1. S.R. acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 138713538 – SFB 881 “The Milky Way System” (sub-projects B1, B2, and B8)), from the Priority Program SPP 1573 “Physics of the Interstellar Medium” (grant numbers KL 1358/18.1, K.L. 1358/19.2), and acknowledges also support from the DFG via the Heidelberg Cluster of Excellence STRUCTURES in the framework of Germany’s Excellence Strategy (grant EXC-2181/1 – 390900948).

Appendix A Evolution sequences of runs 1, 3, 4, 10, and 11

Here we show the evolution of the column density for the run with an infalling angle of α =0°, but prograde infall (run 1; Fig. A.1). Moreover, we show the sequence of retrograde infall for α =35° (run 3; Fig. A.2) and α =60° (run 4; Fig. A.3). Figure A.1 illustrates that the inner disk remains large with respect to the retrograde case in run 2 (Fig. 2). Also, there is no long-lasting gap forming between the inner and outer disk in the prograde case as discussed in the main text. The sequence for infalling angles α =35° and α =60° of both retro- and prograde infall is illustrated in Fig. A.2, Fig. A.3, Fig. A.4, and Fig. A.5. The sequences illustrate that the outer disk forms misaligned with respect to the orientation of the primordial inner disk and that the inner disk becomes smaller for retrograde infall.

Fig. A.1 Same as Fig. 2, but for run 1.

Fig. A.2 Same as Fig. 2, but for run 3.

Fig. A.3 Same as Fig. 2, but for run 4.

Fig. A.4 Same as Fig. 2, but for run 10.

Fig. A.5 Same as Fig. 2, but for run 11.

References

All Tables

All Figures

Fig. 1 Sketch of the initial setup. The initial location of the cloud varies depending on α and the chosen impact parameter b. The staris located at the origin of the coordinate system with the cirumstellar disk in the xy plane. The disk rotation has either a right-handed (corresponding to the prograde case) or left-handed orientation (corresponding to the retrograde case). Here, the illustration shows the setup for α = 60°.

In the text

	Fig. 2 Zoom-in on column densities for run 2 at t = 0, 50, 100, and 150 kyr (from left to right) measured perpendicular to the initial plane of the inner disk with the star at the center. First row: area projection of 8000 AU × 8000 AU. Second row: 1600 AU × 1600 AU. Third row: 320 AU × 320 AU.
In the text

	Fig. 3 Same as Fig. 2, but for run 5.
In the text

Fig. 4 Average azimuthal velocity vϕ profile of the outer disk for run 1 (purple solid line) and run 2 (solid), run 3 (dashed), run 4 (dotted), run 5 (solid) at t = 50 kyr (left panel), t = 100 kyr (middle panel), and t = 150 kyr (right panel). The darker the color of a line, the larger the inclination angle α. We note that vϕ was computed within a vertical extent of |z| < max(6 AU, 0.05r) from the midplane of the disk for 300 logarithmically spaced radial bins in the range from 10 to 104 AU. The gray line is the Keplerian profile for a 2.5 M⊙ star.

In the text

Fig. 5 Projection of the inner disk for runs 1 to 5 (top to bottom) at t = 0, 50, 100, and 150 kyr (from left to right) measured perpendicular to the initial plane of the inner disk and azimuthally averaged column density profile at the shown snapshots (rightmost panels). The azimuthally averaged column density in the rightmost panels was computed within a vertical extent of z = ±6 AU and the darkness of the line corresponds to the time of evolution (from light orange representing t = 0 kyr to dark orange at the end of the simulation). We note that the inner part is affected by accretion onto thecentral star and gravitational softening. We therefore exclude the inner region within r = racc + 5 AU = rgrav,soft + 10 au = 20 AU in the azimuthally averaged column density profiles.

In the text

	Fig. 6 Similar as Fig. 4, but for the inner disk. We note that vϕ was computed within a vertical extent of |z| < 6 AU from the midplane of the disk for linearly spaced radial bins with a size of Δr = 2 AU.
In the text

Fig. 7 Column density Σ over radius r for prograde infall without an inclination angle in run 1 (green colors) and retrograde infall without an inclination angle in run 2 (orange colors) at t = 50 kyr (dashed lines), t = 100 kyr (dotted lines), and t = 150 kyr (solid lines).The darkness of the line colors indicates the evolution in time from 50 kyr (lightest color) to the final time (darkest color) at t = 150 kyr. The black dashed line corresponds to a theoretical column density profile for Σ0 = 170 g cm−2 and p = 1.5 which extends to 2000 AU.

In the text

	Fig. 8 Time evolution of the accretion rate for aligned infall run 1 (purple line) and counter rotating infall run 2 (blue line) with infall angle α = 0°.
In the text

	Fig. 9 Time evolution of azimuthally and radially averaged eccentricity in the inner (top panel) and outer disk (bottom panel). For the inner (outer) disk, the average eccentricity between a radial distance of r = 25 (200) AU to 30 (225) AU and a vertical extent of z = ±0.1r over the entire azimuth is shown.
In the text

	Fig. 10 Evolution of run 2 seen along the x-axis for 105 kyr showing the column density every 15 kyr. The images show how the cloudlet expands and how gas falls toward the central star forming the disk while the cloudlet dissolves.
In the text

	Fig. 11 Inner and outer disk of runs 2 to run 5 (from left to right panel) seen from a viewing angle along the x axis at t = 100 kyr (upper panel) and t = 150 kyr (lower panel).
In the text

Fig. 12 Time evolution of angle δ for the inner (top panel) and outer disk (bottom panel). For the inner disk, δ was measured as the angle of the angular momentum vector at time t with respect to the initial angular momentum vector of the inner disk L50(0). For the outer disk, δ is the angle between the angular momentum vector of gas located within a range from radius r = 200 AU to r = 225 AU and the initial angular momentum vector of the infalling cloudlet.

In the text

	Fig. 13 Same as the upper panel in Fig. 12, but for runs 5 (dashed line), 8 (dotted line), and 9 (solid line) only, i.e., the runs with infall angle α = 90°, and varying the initial density profile of the inner disk. The initial inner disk mass is Mi,disk ≈ 4Mcloudlet in run 5, Mi,disk ≈ 0.4Mcloudlet in run 8, andMi,disk ≈ 24Mcloudlet in run 9. In all cases, Mcloudlet = 1.87 × 10−4 M⊙.
In the text

	Fig. 14 Multiwavelength RGB composite image of run 3 for typical observational bands of 1.66 μm (blue, VLT H-band), 53.0 μm (green, SOPHIA/HAWC+ A-band), and 870.0 μm (red, ALMA band 7). The disks were observed from a distance of 140 pc under an angle of 45° with respectto the main axis.
In the text

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

	Fig. A.2 Same as Fig. 2, but for run 3.
In the text

	Fig. A.3 Same as Fig. 2, but for run 4.
In the text

	Fig. A.4 Same as Fig. 2, but for run 10.
In the text

	Fig. A.5 Same as Fig. 2, but for run 11.
In the text