Clumpy accretion in pre-main-sequence stars as a source of perturbations in circumstellar disks

The development of perturbations in the circumstellar disks of pre-main-sequence stars caused by clumpy accretion was investigated. Here we perform 3D hydrodynamical smoothed particle hydrodynamics simulations of disks perturbed by a recent clump accretion event. These simulations are further explored by radiative transfer calculations to quantify the observational appearance of such disks. It was shown that the density waves in the disks were formed at the fall of the clump. After several revolutions they can transform into spirals and ring structures. Their images in millimeter wavelengths are very similar to those observed with Atacama Large Millimeter/submillimeter Array in some protoplanetary disks. We assume that clumpy accretion may be the source of such structures.


INTRODUCTION
After formation of a young star and its circumstellar disk due to the gravitation collapse of the protostellar cloud the process of accretion of matter from the nearest surrounding the star can continue and have a form of the clumpy accretion. Graham (1992) was probably the first who used this term. He tried to explain by such a way the observations of the strong extinction events observed in some young variables. Later this type of accretion has been considered by Hartmann & Kenyon (1996) for explanation of the FUOR's phenomenon. This idea is quite popular up to now (Zhu et al. 2010;Bae et al. 2013;Hartmann & Bae 2018). The mechanism of formation of chondrule as a result of the clumpy accretion was discussed Tanaka et al. (1998). Obviously, the fall of the clump on the disk should cause disturbances at the place of the fall. It is interesting to trace how this disturbance will develop and what structures on the disk can be caused by it.
The detection of various structures on images of protoplanetary disks is one of the most interesting results obtained with the ALMA interferometer (see, e.g. Cazzoletti et al. 2018;Huang et al. 2018a,b;Pérez et al. 2018). Ring and spiral structures are most often observed. Less commonly, structures resembling highly elongated vortices are observed. A number of papers are devoted to theoretical studies of the formation of such structures. Their formation is associated with perturbations in the disks caused by the motion of the forming planets (e.g., Ruge et al. 2013;van der Marel et al. 2015;Dong et al. 2015;Demidova & Shevchenko 2016;Jin et al. 2016;Dong et al. 2018), the development of various kinds of instabilities in the disks (e.g., Banzatti et al. 2015;Zhang et al. 2015;Birnstiel et al. 2015;Okuzumi et al. 2016;Johansen et al. 2009;Bai & Stone 2014;Takahashi & Inutsuka 2014;Dullemond & Penzlin 2018), a large scale vertical magnetic field (Suriano et al. 2018) or with the destruction of large bodies in collisions (Demidova & Grinin 2019;Nayakshin et al. 2020). In all these papers, the source of disturbance is in the disk itself.
In our paper we discuss an alternative scenario for the formation of the observed structure. We investigate in the first time the dynamical response of circumstellar disk on the perturbation associated with the clumpy accretion events. Using the hydrodynamic simulations we calculate the disk images at 1 mm and discuss the results in the context with interferometric observations of the protoplanetary disks with ALMA.

INITIAL CONDITION
A model consists of a young star of solar mass (M * = M ⊙ ) embedded in a gas disk with total mass is M disk = 0.01M ⊙ . At the beginning of simulation the disk matter was distributed azimuthally symmetrical within the radii r in = 0.5 and r out = 50 AU. The initial density distribution of the disk is where Σ 0 is arbitrary scale parameter, which is determined by disk mass. Hydrostatic scale height is H(r) = κT mid (r)r 3 GM * µmH , where κ, G and m H are the Boltzmann constant, the gravitational constant and the mass of a hydrogen atom and µ = 2.35 is the mean molecular weight (Dutrey et al. 1994). Following Chiang & Goldreich (1997) we determine the law of midplane temperature distribution T mid (r) = 4 γ 4 R * r T * , where γ = 0.05 (Dullemond & Dominik 2004). The calculations were performed in the local thermodynamic equilibrium approximation P (r, z, t) = c 2 (r)ρ(r, z, t), where P is local pressure at the moment t and c is a sound speed. It was assumed that in the vertical direction along z the disk is isothermal. The temperature of the star was assumed to be T * = 5780 K and star radius is R * = R ⊙ . The disk relaxed during 600 years, and then the remnant of the fallen gas clump was added into it.

Impulse approximation to the clump-disk collision
When a clump falls onto a disk, part of its kinetic energy is converted into thermal energy. An infrared spot should appear on the image of the disk where the clump fell. The thermal relaxation time in protoplanetary disks at the distance ≥ 20 AU from the star is much less than the local orbital period (Malygin et al. 2017). Therefore, the remnant of the fallen clump quickly comes to thermodynamic equilibrium with the matter of the disk. It participates in the Keplerian motion of the matter, while maintaining the residual velocity component orthogonal to the plane of the disk.
At the initial moment of time, we assume that the remnant of the clump has already reached thermodynamic equilibrium with the matter of the disk. The remnant was generated as a density perturbation on the disk in the form of a disk segment bounded by radii R 0 and step dR and distributed over the azimuthal angle φ with the axis of symmetry along negative part of x-axis (φ = 30 • for all models). The density of matter in the disturbance exceeded the local density of the disk by a factor of K = Σ cl Σ d , where Σ cl and Σ d are local surface density of remnant and disk respectively. The remnant moves prograde. The particle velocity of the remnant is equal to V (R) = L · V k (R), where V k (R) is Keplerian velocity at a given distance from the star and L is parameter of the problem. The velocity vector had a residual inclination to the disk plane sin(I) = V z (R) V (R) (Fig. 1). The residual angle of inclination of the remnant depends on the initial angle of the fall and the amount of kinetic energy that is spent on heating the disk in the region of fall. We considered a number of the possible options for the value and inclination of the velocity vector. The parameters of all calculated models are given in the Table 1.

METHODS
The evolution of the remnant in gas disk was simulated by the SPH method (smooth particle hydrodynamics). The calculations were performed using the code Gadget-2 (Springel et al. 2001;Springel 2005) modified by us (Demidova 2016). In total, from 5 · 10 5 to 2.5 · 10 6 particles of the gas disk and from 5 · 10 3 to 2 · 10 5 particles of the perturbation were involved in the simulations. The calculations took into account the self-gravity of the disk.
The simulated region was divided by 200 × 30 × 90 cells in spherical coordinates (R, θ, φ), in which the average values of the SPH particle density were determined. We assume that dust particles with a size of 1, 10, 100 microns and 1 mm are well mixed with the gas, and are distributed according to the law dn(s) ds ∝ s −3.5 , where n is the concentration and s is the size of the dust grain (Dohnanyi 1969). The dust to gas mass ratio in the disk was 0.01 as the average in interstellar medium. The dust opacity is calculated using Mie theory for Magnesium-iron silicates (Dorschner et al. 1995). RADMC-3D (Dullemond et al. 2012) code was used for the 3-D radiative transfer calculations. The emergence of the remnant of the clump of matter in the disk leads to the propagation of density waves in horizontal and vertical directions. The strongest perturbations arise at large values of the parameters K, L, and I, as expected. Due to the differential rotation of the Keplerian disk the remnant stretches and transform during the time. A local increase in the surface density leads to the appearance of corresponding large-scale inhomogeneities in the disk images.  Note-The column of "Structures" lists the types of observed asymmetries in the order of their appearance in the disk images. The column of "Lifetimes" is a long-lived structures lifetime.

Perturbations at large radii
For the models discussed here, the initial position of the disturbing remnant of the clump was set equal to R 0 = 20 AU with a step of dR = 5 AU. Calculations have shown that the parameter L, which characterizes the kinetic energy of the falling clump, has the greatest influence on the type of disturbance in the disk. Therefore, we will sequentially discuss the three energy regimes considered in our models.

Sub Keplerian perturbations
In this case, when the parameter is L = 0.8, due to the differential rotation of the Keplerian disk the remnant stretches and transform to a piece of arc reembling a cyclonic vortex and then turn into a spiral during one revolution around the star (∼ 125 years). Since the matter of the disk is involved in its motion the spiral dips and thickening of the density are visible in the disk (Fig. 3). In the central part of the disk, the spiral splits into two branches. The spiral structure quickly (∼ 300 years) twists into an asymmetric ring, which, scattering in the disk, retains the asymmetry until the end of the calculations (600 years). This asymmetry rotates with the disk.
The wave also passes in the vertical direction along the disk, which extends inside and to the edge of the disk. Perturbation twists the central plane of the disk. The maximum distortion of the disk plane occurs near the axis y, but does not coincide with its position. The inner parts of the disk inclines relative to the periphery. Over time, the radius of the outer boundary of the inclined area increases. It distorts 30 AU at the time of 600 years (Fig. 2). In this case, the tilt of the disk plane relative to the initial position is approximately equal to 0.2 • at I = 5 • and 0.9 • in the case of I = 30 • . An increase in the angle of I does not affect the speed of propagation of global perturbation to the   edge of the disk (Fig. 4). Increasing the initial density (parameter K) of the clump increases the vertical distortion of the disk.
The perturbations described above show themselves on the images of protoplanetary disks. The form of asymmetric structures in images corresponds to the perturbations of the density of the disk matter. However, on the periphery of the disk is visible a shadow from the matter above the disk plane. The asymmetric ring-shaped structure on the images has a horseshoe-shaped form (Fig. 5). Calculations have shown that the minimum value of K in which this structure   is visible on the images, equal to 3 (about 0.1 of Jupiter mass). An increase in the parameter K to 5 increases the brightness of the structure and their lifetime, but its horseshoe form is preserved.  For this class of models of the phase of the disintegration of the clump remnant during the first orbit similar to the previous one (Fig. 6). A piece of arc is converted to the one-hand spiral. It twisted into a symmetric ring-shaped structure during the next few revolutions (∼ 500 years).

Keplerian perturbations
An increase of the radius of distortion of the central plane of the disk relative to the initial position is faster than in the previous case. It reaches 40 AU at the time 600 years. The maximum inclination angle is 0.2 • in the case I = 5 • Figure 13. The same as in Fig. 5 for two models with 2.5 · 10 6 particles. The models parameters K = 3, I = 5 • , L = 1 (left) and K = 3, I = 30 • , L = 1.2 (right) at the time 1000 years. and 0.8 • at I = 30 • . The direction of maximum distortion in the vertical direction is also close to the axis y, but does not coincide with it.
The visible asymmetry on the images of the disk are also appropriate for surface density (Fig. 7). In this case, the propagation of the density wave gives multi-lane image of the protoplanetary disk at a certain point in time (right image of Fig. 7). An increase in the angle of I affects the image of the ring-shaped structure. Two symmetric weakly pronounced spirals are visible on the images instead of the ring, if the value of I ≥ 20 (Fig. 8).

Super Keplerian perturbations
In this case, the clump matter motion in the protoplanetary disk causes severe density perturbations and significantly distorts the disk in the vertical direction. As in previous cases, during one convolution of the clump, the vortex-like structure is stretched into a spiral, which is converted into two spirals during the next revolution (Fig. 9). Each of the spirals is logarithmic, they are shifted by phase relative to each other by 180 • . The form of the spirals weakly depends on the angle of inclination of I.
Disk distortion in the vertical direction differs from the models described above for this case. The periphery of the disk is distorted, and the inner parts of the disk deforms weaker. The waves are still propagating along the disk in a vertical direction at the time of 600 years as seen from Fig. 10. In this case, the distortion of the inner region depends on the inclination angle and has the opposite character for I < 20 and I ≥ 20.
It become seen noticeable asymmetry of the spirals on the periphery of the disk in images with an increase in the angle of inclination of I (Fig. 11).
For this class of models, the case K = 1 (corresponds to the mass of the remnant ∼ 12M ⊕ ) was considered. In this case, the perturbation is weaker; however, the two-arm spiral can also be identified in the disk images. It is more pronounced at a larger inclination angle I (Fig. 12).

Long-term dynamics
For the models described above, the number of SPH particles was 5 · 10 5 . These models have a lower resolution compared to the models, the calculations of which involved 2.5 · 10 6 particles. However, the calculations showed that at the initial phases of clump destruction, the images of disks obtained on the basis of models with a small number of particles show the same structures as for more accurate models. However, to study the long-term dynamics of the fallen clump remnant, higher resolution is required. We have calculated two limiting cases that correspond to the parameters that cause the minimum (K = 3, I = 5, L = 1) and maximum (K = 3, I = 30, L = 1.2) disturbance in the disk. Over time, density waves scattered and the all structure settles down to the plane of the disk. The ring of the first model stretches along the radius and loses brightness mixing with the matter of the protoplanetary disk with time. It is faintly noticeable at 1000 years after the fall of the clump (Fig. 13 left). Than it is not visible against the background of the disk matter. For the second model spiral waves are still noticeable after 1000 years (Fig. 13 right), but after ∼ 2000 years they completely disappear. An asymmetric ring structure can be seen on the disk by the time of 2000 years (Fig. 14).
The charateristic time of the disk dynamic relaxation after the fall of the clump also depends on the place of its fall. For example, with R 0 = 30 AU and the same clump parameters as in the previous model, the lifetime of spirals and ring structures generated by its fall increases to 4 × 10 3 years. At R 0 = 50 AU, the characteristic relaxation time of disturbances on the disk is even longer: ∼ 10 4 years.
On the Fig. 15 it is seen that the plane of the disk settles to its original position over time. However, even 2000 years after the fall of the clump, a slight inclination of the disk plane near the axis y remains. For the first model, it is equal to ∼ 0.14 • , and for the second, it is ∼ 0.72 • .

Perturbations at small radii
In this class of models the perturbation was located near the star at the distance R 0 = 10 AU with the step of dR = 2 AU. The clump had parameters φ = 30 • , K = 3 and I = 30 • . The parameter L was varied. The mass of the clump was about 13M ⊕ .   16 shows images of the disk for three models with the L parameter 0.8, 1 and 1.2 after one period (36.5 years). Bright dense structures and areas of shadow, which are caused by matter rising above the plane of the disk, are visible on the disk. But all structures are scattered after next period for models with L ≤ 1. For case L = 1.2 waves propagate along the disk, which in time from 250 to 500 years can be seen in the images as a multi-lane structure (Fig. 17).
In the vertical direction, the disk is distorted in all considered cases. The maximum distortion is achieved near the y axis as for models distant from the star. The Fig. 18 shows the average values of z along the y-axis for two points in time 365 and 1460 years. One can see in all cases, the perturbation propagates outward from the inner part of the disk, tilting its central plane. With an increase in L, the final inclination of the disk increases, but remains within 0.5 • .

CONCLUSION
Calculations have shown that at the initial stages of the remnant of the clump disintegration, the structures that are visible in the images of the protoplanetary disk are similar for the entire set of the models. However, the shape of the final long-lived structure primarily depends on the kinetic energy of the falling clump.  During the first revolution of the center of the remnant (at the initial moment of time), an arc-like structure resembling a vortex is visible in the disk image. Similar structures are observed, for example, in objects HD 135344B (Cazzoletti et al. 2018) and HD 143006 (Pérez et al. 2018). At the next stage of evolution, the image shows a tightly wound spiral, as, for example, in the case of object HD 163296 (Huang et al. 2018b). In the case when the residual velocity of the remnant does not exceed the Keplerian velocity, a ring is a long-lived structure, which can also be asymmetric. In this case, at a certain moment in time, the passage of a wave over the disk can give a multi-lane structure in the image. The ring-shaped structure is visible in the images of a number of objects, for example, HD 169142 (Fedele et al. 2017), HD 97048 (van der Plas et al. 2017), RU Lup, Elias 24, AS 209, GW Lup (Huang et al. 2018b). In the case of a high kinetic energy of a clump, a two-armed spiral appears on the disk image, which, after several thousand years, transforms into an asymmetric ring. Two-arm spirals were obtained on images of objects Elias 27, IM Lup, WaOph 6 (Huang et al. 2018a). The median age of these sources is ∼ 1 Myr (Huang et al. 2018b), and the youngest sources have estimated ages about a few tenths of Myr (Huang et al. 2018a).
Since the velocity vector of the remnant has a residual inclination relative to the plane of the disk, it is distorted. However, over time, the inclination of the plane of the disk relative its original position decreases. The tilt angle at the end of calculations does not exceed 1 • . Probably, for a more noticeable change in the inclination of the disk, a large mass of the clump remnant is required. In this work, we were looking for the minimum mass of the remnant, which develops into a structure visible in the image of the protoplanetary disk. It turned out that in the case of a high-energy fall, the minimum mass of the remnant of the clump is ∼ 10M ⊕ . If the residual velocity of the remnant does not exceed Kepler one, then the minimum mass is ∼ 0.1M J .
Ring-shaped structures formed from the material of the remnant of the fallen clump and the protoplanetary disk are long-living structures and can exist for more than 3000 years. Thus with a sufficiently large mass of the falling clump, the evolution of the disturbance can lead to the formation of a planet on an inclined orbit. It should be noted that Toomre parameter in the disks under consideration is equal to ∼ 40 at a distance of 20 AU. Consequently, in a more massive disk, the fall of a clump several times denser than the material of the disk can trigger the process of gravitational collapse and the formation of a planet.
The fall of the clump near the star can cause not only a FUOR flare, but also a strong increase in circumstellar extinction, leading to a deep and prolonged weakening of the optical brightness of the star and an increase in its infrared radiation. Such events were observed in three young objects: V1184 Tau (Grinin et al. 2009), RW Aur (Shenavrin et al. 2015;Koutoulaki et al. 2019) and AA Tau (Covey et al. 2021).
Let us make a rough estimate of the additional mass of circumstellar matter, which can be added to the mass of the disk during the lifetime of a star due to episodic falls of clumps. Suppose that the average age of stars with protoplanetary disks 10 6 years (Huang et al. 2018b) and the average lifetime of one disturbance 3 × 10 3 years. This will give the probability of observing one event P 1 ∼ 3 × 10 −3 . The probability close to unity will be obtained when ∼ 3 × 10 3 clumps fall during the disk lifetime (∼ 1 Myr). If the clump mass of 0.1M J is necessary to create a strong disturbance (see above), the disk mass will increase during this time by 0.03M ⊙ at an average accretion rate on the disk of 3 × 10 −8 M ⊙ yr −1 . Such an increase in mass is not critical for a typical disk mass of 0.01 − 0.2M ⊙ and an accretion rate onto a young star of ∼ 10 −7 − 10 −8 M ⊙ yr −1 .
Thus, the mechanism of the clumpy accretion in protoplanetary disks can explain the formation of the main types of structures identified in images of protoplanetary disks. In addition, one clump falling at an angle onto a protoplanetary disk can produce multi-lane bright ring-shaped structures. So far, we have demonstrated here the fundamental possibility of obtaining the observed structures on protoplanetary disks in the model of the clumpy accretion. Appreciate the complexity of this process, it needs a more detailed consideration, taking into account the thermal regime in the perturbed region of the disk.