Misaligned disks induced by infall

Arc- and tail-like structures associated with disks around Herbig stars can be a consequence of infall events occurring after the initial collapse phase of a forming star. An encounter event of gas with a star can lead to the formation of a second-generation disk after the initial protostellar collapse phase. Additionally, observations of shadows in disks can be well described by a configuration of misaligned inner and outer disk, such that the inner disk casts a shadow on the outer disk. Carrying out altogether eleven 3D hydrodynamical models with the moving mesh code AREPO, we test whether a late encounter of an existing star-disk system with a cloudlet of gas can lead to the formation of an outer disk that is misaligned with respect to the primordial inner disk. Our models demonstrate that a second-generation disk with large misalignment with respect to an existing primordial disk can easily form if the infall angle is large. The second-generation outer disk is more eccentric, though the asymmetric infall also triggers eccentricity of the inner disk of $e\approx 0.05$ to $0.1$. Retrograde infall can lead to the formation of counter-rotating disks and enhanced accretion. As the angular momentum of the inner disk is reduced, the inner disk shrinks and a gap forms between the two disks. The resulting misaligned disk system can survive for $\sim 100$ kyr or longer without aligning each other even for low primordial disk masses given an infall mass of $\sim 10^{-4}$ M$_{\odot}$. A synthetic image reveals shadows in the outer disk similar to the ones observed in multiple transition disks that are caused by the misaligned inner disk. We conclude that late inclined infall onto a star-disk system leads to the formation of a misaligned outer disk. Infall might therefore be responsible for observations of shadows in at least some transition disks.


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 Müller et al. 2018;Christiaens et al. 2019;Haffert et al. 2019;Pinte et al. 2018Pinte et al. , 2019Pinte et al. , 2020Teague et al. 2018Teague et al. , 2019Pé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). 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 an inner and an outer disks (Avenhaus et al. 2014;Marino et al. 2015;Benisty et al. 2017Benisty et al. , 2018Casassus et al. 2018). Against the background of studies showing that a (sub-)stellar perturber inside (Nixon et al. 2013) or outside the disk (Dogan et al. 2015) can warp and break the disk, 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 Marie Skłodowska-Curie Action global fellow around compact binary systems Facchini et al. 2018;Price et al. 2018). Recently, Nealon et al. (2018) as well as Zhu (2019) also demonstrated that a mildly misaligned planet located in the gap of a disk can lead to break-up of the disk and large misalignment of the inner disk with respect to the outer disk. In particular, Gonzalez et al. (2020) and Nealon et al. (2020b) showed that a combination of external binary and an inner planet located in the gap between 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 tailor spiral-like structures. One of the most prominent example 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). 1 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 Myr) (van den Ancker et al. 1998;De-Warf et al. 2003), making it unlikely that the associated structures are remnants of the early protostellar collapse phase. Some authors speculated that these structures might be signs of a re-arXiv:2110.04309v1 [astro-ph.SR] 8 Oct 2021 A&A proofs: manuscript no. aanda Fig. 1. Sketch of the initial setup. The initial location of the cloud varies depending on α and the chosen impact parameter b. The star is located at the origin of the coordinate system with the cirumstellar disk in the xy-plane. The disk rotation has either right-handed (corresponding to the prograde case) or left-handed orientation (corresponding to the retrograde case). Here, the illustration shows the setup for α = 60 • . cent 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 at most modest and short-lived misalignment, hence the authors rule out flybyinduced misalignment.
Instead, the origin of the misalignment in these systems is likely different. Luminosities of some evolved stars are significantly higher than stellar evolution models predict (Kenyon et al. 1990). Several studies showed that such enhancements in luminosity can be caused by episodic accretion events, where shortterm 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 & Mc-Kee 2011;Padoan et al. 2014;) 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 can not 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 different orientation around an existing inner disk.
Section 2 describes the method used in this work. In section 3, we present the results of our numerical models. Based on the results of the models, we compare the results with observations in section 4 and discuss limitations as well as prospects of our study. In Section 5, we summarize the results and present the conclusions.

Methods
In this study, we analyze the effect of infall onto an already existing primordial disk. We carry out hydrodynamical simulations of a cloudlet of gas that approaches a star-disk system with different infalling angles. In a final step, we post-process one snapshot of the simulation and produce 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 are carried out with the moving mesh code arepo 2 (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 is set to 10 24 g (parameters RefinementCriterion = 1, Derefine-mentCriterion = 1, ReferenceGasPartMass = 10 24 and Tar-getGasMassFactor = 1 for setting the code units equivalent to cgs-units), i.e., cells with twice/half the target mass are refined/derefined. The cells are allowed to have a maximum difference in volume of a factor of 5 (MaxVolumeDiff = 5), while the minimum cell volume is set to 10 37 cm 3 and the maximum cell volume to 4.54 × 10 50 (MinVolume = 10 37 , MaxVolume = 4.54 × 10 50 ).
Analogous to the previous study presented in KGD20, we insert a star with mass M * = 2.5M modeled as a point mass at the center of our domain and the star is fixed at its location throughout the simulations. We use the accretion recipe introduced in KGD20, which means that the star accretes 90% of the mass that is located within its accretion radius at every timestep. The accretion radius is 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 apply periodic boundary conditions. The simulations are stopped early enough that they are unaffected by the boundary conditions. The cloudlet approaches the star with velocity while the background gas is at rest with respect to the central star. The center of the cloudlet with initial radius R cloudlet = 887 AU is placed at location where b = 1774 AU is the impact parameter and α is the infall angle measured with respect to the xy-plane of the coordinate system. 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.5M star with velocity v i would be deflected by 90 • (for more details see section 2 in (Dullemond et al. 2019)). Adopting a slightly modified version of the description in Klessen & Hennebelle (2010), the cloudlet mass is defined as  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.
The density of the cloudlet ρ cloudlet is uniform and the density of the background gas is defined as ρ bg = 1 800 ρ 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 | 2 ∝ k −4 , with wavenumber of perturbation k. The power index is 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 where σ v is the velocity dispersion of 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 insert a primordial isothermal disk (T = 10K) with radius 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 where Σ 0 is the amplitude of the column density and p is the power index. Vertically (i.e., in z-direction), the density of the disk drops exponentially such that the density profile of the disk is defined as for radii in the range of r in = 20 AU ≤ r ≤ R disk = 50 AU. H(r) is the scale height of the disk given by where the sound speed c s for an isothermal gas is constant, with Boltzmann constant k B , mean molecular weight µ = 2.3 and hydrogen mass m H . Ω is the orbital frequency of a Keplerian Article number, page 3 of 20 A&A proofs: manuscript no. aanda disk given by with gravitational constant G, i.e., H(r) scales as H(r) ∝ r 3/2 for the isothermal disk considered in the simulations. H 0 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 is tapered off such that the density drops according to 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 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, i.e., ρ disk (18 AU, z) = 0.5 × (ρ disk (20 AU, z) + ρ bg ). To avoid spurious effects from cells in the vicinity of the star, we force 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 is tapered off exterior to r > R disk = 50 AU according to 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 use logistic functions for the tapering to ensure that the density profile asymptotically approaches both the values at the inner/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 vary the initial position of the cloudlet angle with the angle α. The tested angles are α = 0 • , 35 • , 60 • and 90 • . The primordial disk is on all cases aligned with the xyplane.
In this study, we define Σ 0 = 170 g cm −2 for the main runs, which is 0.1 times the column density at 1 AU of the minimum mass solar nebula and we perform simulations using p = 1.5. The rotation axis of the inner disk is parallel or anti-parallel to the positive z-axis of the coordinate system. A parallel/antiparallel configuration corresponds to prograde/retrograde infall. To test the influence of the inner disk mass, we also perform simulations with very low inner disk mass of Σ 0 = 17 g cm −2 for an infalling angle of α = 90 • , and without inclination α = 0 • for prograde or retrograde infall. All simulations are carried out for t = 150 kyr of evolution. Moreover, we carry out a comparison run for with α = 90 • inclination for a shorter time sequence of t = 100 kyr using Σ 0 = 170 g cm −2 , but p = 1. The chosen parameters for the altogether eleven runs are summarized in table 1.

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 run 2 and 5 in Fig. 2 and Fig. 3. (For comparison, the sequences for run 1, 3, 4, 10 and 11 are shown in appendix A.)

Outer disk
The figures show that in all cases the infalling cloudlet leads to the formation of a second-generation disk associated with largerscale 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 time scales 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 compute the azimuthal velocity v φ . The azimuthal velocity is 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 still is 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 inwards 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).

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. Fig. 5 shows face-on projections for run 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, mass located within the accretion radius r = 15 AU accretes onto the central star. As not all of the mass within the accretion radius is immediately accreted, temporary features such as rings and small spirals of low-density form within the accretion radius. However, as the dynamics of the gas within r acc = 15 AU are strongly affected by the numerical parameters of the accretion recipe and gravitational softening, we exclude this 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 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 section 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 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 section 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, non-inclined 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 is the size of 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 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 r [au] 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 is the amplitude 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 towards the end of the simulation.

Counter-vs 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 run 1 and 2 in more detail. Generally, run 2 shows the important result that two counterrotating 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 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, i.e. on scales of the transition between outer and inner disk, we can see substantial differences between the two extreme scenarios. The transition between outer and inner disk is smooth in the aligned infall case. In contrast, outer and inner disk in run 2 are more distinctly separated by a transition zone of ∼ 10 AU in width. As shown in Fig. 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 (run 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. 2016Wijnen et al. , 2017a. The transition zone in 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 interac-tion, 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 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) showed 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 done with Hα emission for PDS 70 ). In the majority of cases in our simulations, the gap region is smaller than the inner disk. Note that the initial size of 50 AU for the inner disk is a somewhat arbitrary choice in our setup and we expect larger ratios of 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 inner to outer disk means that material moves radially inwards, 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 measure 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, that are the two extreme cases of prograde and retrograde infall without inclination. In agreement with expectation, 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 expected from the significant mass difference of 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  5. Projection of the inner disk for run 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 is 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). Note that the inner part is affected by accretion onto the central star and gravitational softening. We therefore exclude the inner region within r = r acc + 5 AU = r grav,soft + 10 au = 20 AU in the azimuthally averaged column density profiles. Fig. 7, which shows the radial profile of column density for run 1 and 2, i.e., 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.

Eccentricity
The infall onto the primordial disk is asymmetric. As a consequence both outer and inner disk become eccentric as seen in the projection plots shown in Fig. 2, Fig. 3  ig. 6. Similar as Fig. 4, but for the inner disk. v φ is computed within a vertical extent of |z| < 6 AU from the midplane of the disk for linearly spaced radial bins of size ∆r = 2 AU.  this aspect further, the eccentricity e is computed for each grid point according to where is the specific orbital energy h is the specific angular momentum where µ red = M * m M tot is the reduced mass with m being the mass of the individual grid cell and M tot = M * + m ≈ M * = 2.5M . 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 with 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, i.e., 35 • and 60 • remain slightly more eccentric between 20 to 70 kyr than the other three runs with either 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.
The inner disk remains in all cases less eccentric than the newly forming second-generation disk. Comparing the eccen- tricity in the inner disk of the different runs, however, shows differences. The eccentricity rises more rapidly from an initially circular orbit for the runs without any inclination angle. Also, the eccentricity rises more quickly for the run with 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. Note that in the case of retrograde infall at zero inclination angle (run 2), the drastic increase in eccentricity shown at radial distance between 25 to 30 AU after t = 130 kyr is not associated with the disk anymore, but 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 with the compression of the inner disk. In contrast, e continuously increases for the high inclination cases for both retro-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 on average slightly more rapidly for prograde infall than retrograde infall. 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.
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 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 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 non-zero inclination angle with respect to the disk plane in the 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 compare the orientation of the angular momentum vector of the inner part at different times with each other. In particular, we compute the angular momentum vector L 50 (t) for all the gas located within a radial distance r of 15 AU < r < 50 AU from the central star. Fig. 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 50 (0). The plot shows that the inner part of the gas undergoes mild changes in all cases. Consistently, 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 degrees 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 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 that it exceeds the angular momentum of the forming outer disk, 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 study for the future. Especially, the question 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 vations of warped disks around young stars. In our models, the inner disk changes orientation as a whole, but does not become warped. This is consistent with previous models of small disks with low viscosity (wave-like regime), where the disk globally bends 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 section 4.3).

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 is 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, e.g., causing flaring of the disk, but we showed in KGD20 that heating from the central star only marginally affects the formation of the second-generation disk due to its distance to the heating source. ig. 12. Time evolution of angle δ for inner (top panel) and outer disk (bottom panel). For the inner disk, δ is measured as the angle of the angular momentum vector at time t with respect to the initial angular momentum vector of the inner disk L 50 (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.
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 secondgeneration 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 can 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 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(Liu et al. , 2018 strongly suggest the occurrence of late accretion events. Consistently, Pineda et al. (2020) 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) showed 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.

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; 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 its 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 Article number, page 11 of 20 A&A proofs: manuscript no. aanda magnetohydrodynamical models that successfully reproduce the stellar initial mass function (IMF) Pelkonen et al. (2020) find 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(Kuiper et al. , 2016Rosen & 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).

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, as well as 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 orbital plane of inner binary and outer disk is HD 142527 (Marino et al. 2015). Also, Alves et al. (2019) find a flattened envelope structure that appears to be misaligned to the orbital plane of the close binary [BHB2007] 11. Carrying out star-forming models, Bate (2018) finds misalignment of 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 in understanding whether such misaligned configurations can be formed through late infall events.
In the last 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 . Moreover, Pontoppidan et al. (2020) reported on the variability of the great shadow in the Serpens Molecular Cloud (Pontoppidan & Dullemond 2005) strongly hinting at a protoplanetary disk as the cause of the shadow.
In Fig. 14 we present a multi-wavelength 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-process the results of our hydrodynamical runs with the Monte-Carlo (MC) code polaris 3 (Reissl et al. 2016) that 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 utilize the central star and an additional diffuse interstellar radiation field (ISRF) (Reissl et al. 2016(Reissl et al. , 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 0 = 1 typical for the Milky Way. For the dust we assume spherical grains with a composition of silicates and graphite and apply 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, Rumpl, & Nordsieck 1977). How-3 http://www1.astrophysik.uni-kiel.de/~polaris/ ever, circumstellar disks are environments that seem to support significant grain growth (Ossenkopf & Henning 1994;Kataoka et al. 2017). Hence, we use an upper grain size of a max = 10 µm (Kataoka et al. 2017). The dust mass to gas mass ratio is uniformly set to 1 % (Bohlin et al. 1978) everywhere in the mesh. Finally, we calculate intensity maps via a MC approach (Reissl et al. 2016(Reissl et al. , 2020 for scattered stellar radiation in the near-IR and via a ray-tracing scheme in the sub-mm for the thermal dust emission (see Appendix B in Reissl et al. 2019, for details). Here, we zoom into the inner 800 AU around the star and assume a distance to the observer of 140 pc.

AU
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 e.g. 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 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 sub-mm 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 section 4.1) because a non-isothermal 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 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 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 still occurs for substantial infall angles. In cases of retrograde infall with 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 equation 10 and A3 in Nixon et al. (2013), see also equations 2 and 3 in Facchini et al. (2018)). Kraus et al. (2020) assume relatively large values of α SS ≈ 0.01 to 0.02 for their simulations and apply 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 (equation A3 in Nixon et al. (2013), see also equation 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 bendingwave 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 et al. 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. Note that different from the scenario of planet-induced misalignment considered by Zhu (2019), the orientation of the inner disk would remain unchanged.

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 account for the protostellar environment. Multi-scale simulations similar to the ones presented in Kuffmeier et al. (2017Kuffmeier et al. ( , 2019, but for later times of evolution are therefore the crucial next step. Carrying out simulations for a sample of objects also allows 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.

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 to 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 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 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 that shrinks the inner disk and enhances 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. Appendix A: Evolution sequences of run 1, 3, 4, 10 and 11 Here we show the evolution of the column density for the run with 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 longlasting gap forming between 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.
Article number, page 15 of 20 A&A proofs: manuscript no. aanda Article number, page 20 of 20