Numerical modeling of gravitational instability outcomes in multiphase circumstellar discs

To suggest consistent route toward planetesimal and planetary core formation in circumstellar discs we study gravitational instability outcomes in massive multiphase (gas-collisionless bodies) disc. Such unstable massive disc can be formed together with protostar in molecular cloud collapse with increased ratio of solids to gas density. In our calculations we found regimes when low-massive solid bodies subdisc drastically affect global structure formation in the disc, whose mass is constituted mainly by the gas. We demonstrated also that solitary areas of high gas density can concentrate solids, producing multiphase clumps, which can be considered as a cradle of large bodies formation.


Introduction
Circumsolar disc evolution goes via growth of solid bodies -from nanometer dust grains to kilometersized planets. Formation of large bodies (planetesimals and planet embryos) from meter-sized boulders -major bottleneck in the planet formation process. To suggest consistent route towards planetesimals and planetary cores formation in circumstellar discs we study gravitational instability outcomes in massive multiphase disc. Such unstable massive discs can be formed together with protostar in molecular cloud collapse with increased ratio of solids to gas density. Possible mechanisms of rapid gas cooling (thermodynamic processes) and enrichment the gas with solid phase is discussed in [3].

Mathematical model of circumstellar disc at the stage of clump formation 2.1 Basic equations
Computational experiments reported in the paper were carried out within a quasi-3D model of the disc. It means that we neglected the vertical motion of matter and considered dynamics of the disc where its entire mass is concentrated inside the equatorial plane of the system. Gravitational potential of the disc is calculated using the 3D Laplace equation. Quasi-3D approximation is appropriate for the short-term (some turnovers) processes in the two-phase disc with central body [1]. Gas component is described by equations of gas dynamics: Equations of gas dynamics include surface quantities that are obtained from volume quantities by integration with respect to vertical coordinate z: σ par,gas = +∞ −∞ ρ par,gas dz; p * = +∞ −∞ pdz. Here, v = (v x , v y ) is the two-component gas velocity, and p * is the surface gas pressure. T * = p * σ , S * = ln T * σ γ * −1 are the quantities similar to gas temperature and entropy. γ * is a 2D version of γ, which is related with constant ratio of specific heats as γ * = 3 − 2/γ. Φ -is the gravitational potential in which the motion occurs.
In our model, a solid-phase subdisc is represented by 1-10 m solids moving in such a manner that two solids collide with a frequency not exceeding one event per turnover round the protostar. Then, the dynamics of subdisc of primary solids can be described by the Vlasov equation neglecting collisions of solids at the times of several turnovers: where a = −∇Φ, a -the acceleration of particles in external and self-consistent field, u = (u r , u φ ) is the two-component velocity of particles, f = f (t, r, u) is the function of particle velocity distribution u at a point with r coordinate of the disc. This function is related with the surface density of particles as σ par = f dudz. Note that radius of solids does not appear explicitly in the Vlasov equation; however, the description of subdisc of primary solids as a collisionless system implies a certain size range of these solids.
Φ is the gravitational potential, which is the sum of potential of the motionless central body and potential of the disc, Φ = Φ 1 + Φ 2 , Φ 1 = −M c /r, M c is the mass of central body. Φ 2 is the potential of self-consistent gravitational field, which is found as a solution of mixed problem for the Laplace ∂Φ 2 ∂z | z=0 = 2π(σ par + σ gas ). As was noted above, the gaseous and solid-phase subdiscs at a protostellar stage are considered to be thin discs. The equations are written with dimensionless variables. The basic dimensional quantities are G -the gravitational constant and R 0 = 10 au, M 0 = M = 2 · 10 30 kg, which are typical size and mass of the system.

Boundary and initial conditions
At zero time, surface temperature and density of the disc are specified. In calculations presented in the paper, the density of gas and subdisc of primary solids was taken as a Mackloren disc of mass M par,gas and radius R: σ par,gas (r) = Gas temperature at zero time is specified as T * (r) ∼ σ(r), using a given T 0 , which is the temperature in the 'center' of the disc.
Initial velocities of solids are specified as a sum of regular and chaotic components, u = u' + u", where u' is regular velocity, and u" is chaotic one. Gas velocity and regular velocity of particles are determined from the equilibrium of centrifugal and centripetal gravitational forces:

Numerical methods and setups
To follow the disc evolution we used PIC-SPH code Sombrero. In this code we implemented PICmethod to solve Vlassov equation, we used iterative combined method which includes Fourier transformation and sweep to solve Poisson equation. Gas component of the disc was simulated by SPHmethod [2]. We used cubic spline for 2D space as an interpolation kernel. The smoothing length was varied in space and time according to the formula h i = 2 m i 10 −3 + σ i , which was implemented through arithmetic kernel averaging. As a form of energy equation we used entropy formulation together with summation interpolant to calculate surface density. To prevent particle interpenetrations we applied standard artificial viscosity with α = 1, β = 1, introducing it to the equation of motion only. For nearest neighbor search we used linked-list algorithm with efficient method of particles assortment on each time step.

Results of modeling
Protoplanetary discs can experience different type of gravitational and convective instabilities that can produce structures of different patterns. Solitary areas of high density can be considered as a cradle for large body and planetary embryo formation possibly providing conducive environment for sticking mechanisms and local gravitational collapses. In the paper we want to underline 'mutual-influence' effect of gas and solids moving under the common gravitational field in the disc. First -and most surprising one -effect shows that low-massive solid bodies subdisc can drastically affect structure formation in the disc, whose mass is constituted mainly by the gas. Second -that gaseous clumps even before their collapsing capture solid bodies by their gravitational field. The Jeans length forms several criteria of local and global instability development in discs and can be used for instability outcome prediction. For gas and collisionless medium it can be estimated as Gσ gas . Jeans length for hybrid system, when gas and solid move under the common gravitational field was estimated in [3]. It was found that effective Jeans length is a nonlinear combination of gas and solid bodies quantities: This formula means in particular that not most massive, but most unstable component of the system defines effective lambda -Jeans length for hybrid system. In molecular cloud 98% of mass is presented by gas and only 2% by dust. Due to sedimentation of dust in the equatorial plane the ratio of 07005-p.3

Instabilities and Structures in Proto-Planetary Disks
EPJ Web of Conferences Figure 2. Global surface density of gas (first line), and angular velocity of individual gas clumps in clump coordinate system (second line). T 0 = 0.026, γ * = 5/3 for the first column, T 0 = 0.01, γ * = 5/3 for the second column, T 0 = 0.026, γ * = 11/10 for the third column. The disc is resolved by 640000 SPH particles, 40000000 PIC particles, the gravitational field is calculated using 400 × 512 × 400 grid cells. solids to gas surface density can vary in the range of 0.01-0.1. To consider regime of clump formation, caused by the presence of solids in the disc, we modeled initially unstable massive subdisc of gas and boulders. We reproduced the dynamics of some turnovers of the disc with mass M = 0.55M 0 and radius R = 2R 0 = 20 au rotating round the central body with mass M c = 0.45M 0 . For quasi-3D-model specifying the ratio of surface densities and mass of entire disc, we determined the masses of gas and subdisc of primary solids, M gas = 0.52M 0 , M par = 0.03M 0 . Initial temperature was T * (r) ∼ σ(r) with T 0 = 0.026 in the center, γ * = 5/3.
Let us compare results of two calculations (see Fig.1), where we varied parameter of solids only. In the first case we used velocity dispersion of solids v d = 0.01 and in second case -v d = 0.2 (in dimensionless parameters c s ∼ 0.1). And in first case we see individual clumps formation, when in second -spiral sleeves only. We underline it again that for the same initial parameters of gas -massive component of the system -we got two different outcomes of gravitational instability development. Varying parameters of the system (gas temperature in the center and γ * ), we demonstrated that all gaseous clumps on the early stage of their existence are rotating objects (vortices) with pressure and density maximum in the center and with solid-bodies rotation around the center independently of their size and formation time (see Fig.2). Gravitational potential of such vortices capture solid bodies and produce epicyclical trajectories for some of them. We found that the efficiency of such solid's capturing for non-collapsing gaseous objects strictly depends on initial velocity dispersion of solids: the less initial dispersion is, the more portion of solids is kept during the movement of the clump.