Dust growth and pebble formation in the initial stages of protoplanetary disk evolution

,


Introduction
Protoplanetary disks form during the collapse of slowly rotating prestellar cloud cores.These objects are mainly composed of molecular hydrogen and helium with a small admixture of interstellar sub-micron-sized dust grains.During the initial stages of evolution, the disk is surrounded by a massive gas-dust envelope, which makes it difficult to directly observe the processes that take place in the disk vicinity.Yet, this stage may be crucial for the subsequent processes of planet formation, because here the initial stages of dust growth from sub-µm grains to millimeter-sized aggregates should occur.The accelerated formation of protoplanets is indirectly supported by the detection of ring-and gap-like structures in young protoplanetary disks, such as GY 91 (Sheehan & Eisner 2018), IRS 63 (Segura-Cox et al. 2020), andHL Tau (ALMA Partnership 2015), all of which are no older than 0.5 Myr.These structures can presumably be formed by the gravitational torques from yet unseen planets with a mass equal to or greater than Neptune (Zhu et al. 2012;Dipierro et al. 2015;Dong et al. 2018).
Furthermore, recent estimates of the dust mass in young Class 0 and I disks and also in more evolved T Tauri disks revealed that the dust mass systematically decreases with the age of the system (Tobin et al. 2020;Tychoniec et al. 2020), in agreement with earlier predictions (Greaves & Rice 2011).In particular, the available mass of dust in Class 0 and I disks is comparable to that contained in exoplanets around mainsequence stars, while the dust mass in T Tauri disks is insufficient to explain the bulk mass distribution of exoplanets.This again accentuates the importance of the earliest stages of disk formation in the planet formation paradigm.
have indeed been inferred observationally in protostellar clouds (Galametz et al. 2019).The physical mechanism to explain their existence has not been fully understood.Using smoothedparticle hydrodynamics simulations, Bate (2022) found that dust growth in collapsing clouds occurs only in the inner few hundreds of astronomical units, and in particular in the first hydrostatic core where dust can quickly grow in excess of 100 µm.On the other hand, large grains can be entrained with the disk wind into the protostellar envelope as was recently found by Tsukamoto et al. (2021) and Koga & Machida (2023).
Efficient dust growth in the disk formation phase has been reported in the thin-disk hydrodynamics simulations of Vorobyov et al. (2018), who demonstrated that, depending on the value of the turbulent α parameter and the threshold velocity of dust fragmentation, dust can grow from sub-µm grains to millimeter-or even centimeter-sized aggregates already in the embedded phase of disk evolution.Fast growth agrees with analytic expectation, according to which the dust-size doubling time is inverse proportional proportional to the dust-to-gas ratio and the local Keplerian angular velocity (Birnstiel et al. 2016).Using three-dimensional (3D) magnetohydrodynamics simulations, Lebreuilly et al. (2020) and Koga et al. (2022) found efficient settling and the drift of grains with sizes greater than a few tens of microns in the first core and in the newly formed disk, producing dust-to-gas ratios a factor of several greater than the initial value.Marchand et al. (2023) employed a precomputed dust growth model and show that grain sizes can reach ≥ 100 µm in an inner protoplanetary disk 1 kyr after its formation.Tsukamoto et al. (2023) found that dust growth can influence the disk evolution by changing the coupling between gas and the magnetic field, making the gas density profiles steeper and the disk mass smaller.
The fast dust growth solves only part of the planet formation problem.Grown dust must avoid a catastrophic inward radial drift caused by the mismatch between the angular velocities of gas and dust in a Keplerian disk (Weidenschilling 1977).This can be achieved if the disk pressure has local maxima or strongly varying spatial gradients (e.g., Dr ążkowska et al. 2016;Flock et al. 2015).Water ice lines can also help to reduce the inward radial drift (Pinilla et al. 2016;Molyarova et al. 2021;Lau et al. 2022).Spiral arms and gaseous clumps in gravitationally unstable disks can also collect dust particles if they grow sufficiently large to be characterized by the Stokes number ∼0.1-1.0 (e.g., Rice et al. 2004;Baehr & Klahr 2019;Deng et al. 2021;Baehr et al. 2022), though the efficiency of this process may be weakened by the transitory nature of the spiral arms and tidal destruction of gaseous clumps in young protoplanetary disks (Vorobyov et al. 2018;Vorobyov & Elbakyan 2019).
In our previous works, we have studied dust dynamics and growth in young protoplanetary disks using two-dimensional thin-disk approximation (Vorobyov et al. 2018(Vorobyov et al. , 2023a;;Vorobyov & Elbakyan 2019;Elbakyan et al. 2020;Molyarova et al. 2021).These studies have brought important insights in the spatial distribution of pebbles, efficiency of dust accumulation, composition of ice mantles, and dust growth efficiency in the disk midplane.Here, we present our first results of dust dynamics and growth obtained in fully three-dimensional hydrodynamic simulations of a gas-dust disk in the initial stages of its formation.We developed an original code based on numerical techniques that are different from our previous thin-disk model, which allows us to test our previous conclusions and obtain new insights into the three-dimensional distribution of dust and gas in young protoplanetary disks.
The paper is organized as follows.In Sect.2, we describe the basic equations entering our model of a gas-dust disk.We also provide a brief explanation of the original code and the novel Coarray Fortran methodology for distributed memory parallelization.In Sect.3, we present our main results.The caveats and future developments are briefly discussed in Sect. 4. The main conclusions are summarized in Sect. 5.In Appendices A-C, the details of the numerical solution and an extensive suite of test problems are provided.

Model description
In this section, we present the basic equations and the solution methods, which constitute the foundation of the Formation and Evolution Of Stars and Disks on nested grids (ngFEOSAD) code.ngFEOSAD computes the co-evolution of gas and dust during the gravitational collapse of a rotating prestellar core including the initial stages of protoplanetary disk formation using a fully three-dimensional approach.We employed Coarray Fortran for the distributed memory parallelism, which is much simpler in coding than the usual message-passing interface (MPI), while being only slightly inferior to the MPI in terms of the speed.This allows focusing on the solution of scientific problems rather than on the intricacies of parallel programming.Thanks to using fast-Fourier transforms to solve the gravitational potential on nested meshes, ngFEOSAD is fast and environmentally friendly, requiring only a few modern computer nodes to run the models.In particular, we used four nodes with three tasks per node and 16 cores per task, amounting to 192 cores in total on AMD EPYC 7713 processors.For a typical run-time of 15 days, this corresponds roughly to 70 000 core hours.We adopt a barotropic equation of state, which includes the initial isothermal collapse stage followed by the star and disk formation stages (see Eq. ( 3)).We note that the forming protostar is not fully resolved in our simulations with limited spatial resolution (≥0.14 au in the central regions).Yet we do not introduce the sink cell as the results may be sensitive to the details of the sink cell implementation (Machida et al. 2014;Vorobyov et al. 2019;Hennebelle et al. 2020), which requires a careful consideration.Our approach is valid for the very early stages of evolution, when the central object may still be a first hydrostatic core (∼5 au) or in the early stages of collapse toward a protostar.

Gas component
The dynamics of gas is followed by solving the equations of continuity and momentum where ρ g is the volume gas density, u the gas velocity, Φ the gravitational potential, P the gas pressure, I the unit tensor, and f the friction force per unit dust mass between gas and dust.
The gas pressure is related to the gas density via the following barotropic equation A202, page 2 of 20 Vorobyov, E. I., et al.: A&A, 683, A202 (2024) Table 1.Parameters of the barotropic relation.
1.6667 2.5 × 10 12 3 1.4 1.5 × 10 15 4 1.1 -Notes.The coefficients γ i and n c,i were taken from Machida & Nakamura (2015).The number density n c,i is related to the volume density ρ c,i as ρ c,i = µm H n c,i , where m H is the mass of hydrogen atom.
where c s,0 = RT in /µ is the initial isothermal sound speed, T in the initial gas temperature, R the universal gas constant, and µ = 2.33 the mean molecular weight.The values of the indices i and k distinguish the four individual components of the piecewise form shown in Table 1.The first component corresponds to the initial isothermal stage of cloud core collapse, while the other three components account for the postcollapse nonisothermal evolution with various ratios of the specific heat γ i(k) (Masunaga & Inutsuka 2000;Machida & Nakamura 2015).We note that when k = 1 the product term is unity by definition, and the pressure reduces to P 1 = c 2 s ρ γ 1 g .We also formally set ρ c,0 to zero.The gas pressure P is set equal to P k , where the index k is defined by the condition that ρ g falls in between the corresponding critical densities, ρ c,k−1 ≤ ρ g < ρ c,k .

Dust component
For the dust component we use the model originally developed in Vorobyov et al. (2018) with modifications to account for the backreaction on gas in Stoyanovskaya et al. (2018).In this model, dust is considered as a pressureless fluid (see Vorobyov et al. 2022, for justification) with two populations: small dust that is dynamically linked to gas and grown dust that can dynamically decouple from gas.In particular, small dust are grains with a size between a min = 5 × 10 −3 µm and a * = 1 µm and grown dust are aggregates ranging in size from a * to a max .The latter is the maximum value, which is calculated using a simple dust growth model and can be variable in space and time.The value of a * that separates small dust from grown dust is chosen so as to make small dust dynamically coupled to gas.In the gas-dust disk simulations, dynamics of small dust is often considered in the terminal velocity approximation (e.g., Laibe & Price 2014;Lin & Youdin 2017).However, as was shown in Vorobyov et al. (2022, see Appendix D), the drift timescales of µm-sizes dust in a protoplanetary disk is comparable to its typical age.Hence the full coupling approximation is justified on short evolutionary timescales considered in this work.We also note that the small dust-to-gas mass ratio never exceeds 0.01 in our model and is often much lower.In particular, the corresponding local values are in the 10 −3 -10 −4 limits in the bulk of the disk and approach 10 −2 from below only in the disk outermost regions and in the disk upper atmosphere (see also Sect. 3.3).We note, however, that the settling instability may develop at the dust-to-gas mass ratios as low as 10 −2 (Squire & Hopkins 2018) and neglecting the dynamics of small dust and its backreaction on gas may not be justified in this case.Initially, all dust in a collapsing prestellar cloud is in the form of small dust, which can grow and turn into grown dust as the collapse proceeds, and the disk forms and evolves.It is assumed that dust in both populations is distributed over size according to a simple power law: where C is a normalization constant and p = 3.5.We note that the power index p is kept constant in this study, but may also vary if the form of this variation can be inferred from a more sophisticated approach that involves solving for the Smoluchowski equation for multiple dust bins (e.g., Dr ążkowska et al. 2019).
The dynamics of both dust populations is followed by solving the equations of continuity and momentum where ρ d,sm and ρ d,gr are the volume densities of small and grown dust, respectively, u the velocity of grown dust, and S (a max ) the conversion rate between small and grown dust populations.We note that the presence of the S (a max )u term in the last equation and the absence of the mirrored term for the small dust (recall that we assume full coupling for small dust and do not solve the corresponding dynamics equation) formally means that the dust momentum is not conserved.This should be a reasonable approximation for as long as the mass of small dust component is much smaller than that of grown dust, which is generally fulfilled as soon as the disk forms and dust grows (see Figs. 9 and 11).
The drag force per unit dust mass f can be written as where σ is the dust grain cross section, m d the mass of a dust grain, and C D the dimensionless friction parameter.In this work, we consider only the Epstein drag regime, in which case the friction parameter can be written as C d = 8v th /3|u − u| (Weidenschilling 1977), where v th is the mean thermal velocity of gas.The drag force then takes a simple form where t stop is the stopping time expressed as where ρ s = 3.0 g cm −3 is the material density of dust grains.Since grown dust in our model has a spectrum of sizes from a * to a max , the size of dust grains a in Eq. ( 10) should be weighted over this spectrum.However, we use a max when calculating the stopping time.Our choice is justified as follows.Since we use only two dust size bins, the span between a * to a max may become as large as several orders of magnitude during the disk evolution and the average size of dust grains √ a * a max becomes much smaller than a max .However, we are interested in the dynamics of dust grains that are the main dust mass carriers.For the chosen slope p = 3.5, large grains near a max mostly determine the dust mass (Birnstiel et al. 2016).Therefore, we use the maximum size of dust grains a max when calculating the values of t stop .A more consistent approach requires introducing multiple narrower bins for grown dust and setting the mean dust size for each bin, but this is outside the scope of the current work.

Dust growth model
The ngFEOSAD code includes a simple model for dust growth.Initially, in the prestelar cloud, all dust is in the form of small dust grains with a maximum size a max = a * .As the cloud begins to collapse and especially during the disk formation stage the maximum size of dust grains is allowed to increase.We use the following equation describing the time evolution of a max where the rate of dust growth due to collisions and coagulation is computed using a monodisperse model (Birnstiel et al. 2012) This rate includes the total volume density of dust ρ d = ρ d,sm + ρ d,gr , dust material density ρ s , and relative velocity of grainto-grain collisions defined as u rel = (u 2 br + u 2 turb ) 1/2 , where u br and u turb account for the Brownian and turbulence-induced local motion of dust grains, respectively.In particular, u turb has the biggest contribution to dust growth in our model and is expressed as (Ormel & Cuzzi 2007) where c s is the adiabatic speed of sound, α is the parameter describing the strength of turbulence in the disk, and St is the Stokes number defined as St = t stop Ω K .Here, Ω K is the Keplerian angular velocity.We emphasize here that this definition of the Stokes number crucially depends on the assumption of a Keplerian disk.However, young disks may deviate from the Keplerian pattern of rotation.We tried to use the local angular velocity as a proxy of Ω K but that worked poorly in disks perturbed by spiral density waves and mass infall from the envelope.
In this paper, we set Ω K equal to that of a 0.05 M ⊙ central object.
Considering that we limit our runs to the very early stages of star and disk formation, this should introduce only small errors to u turb .We also note that the Stokes number in the infalling envelope should be normalized to the free-fall time rather than to Ω K .This, however, requires determining the disk-to-envelope interface on the fly and is subject to future code improvements.The value of α is a free parameter in our model and it describes the strength of the magnetorotational instability (MRI) in the disk.Since the MRI can only be directly studied in focused high-resolution magnetohydrodynamics simulations (e.g., Bai 2013), we use a parametric approach and set α = 0.001 as a constant of space and time.This value is also consistent with observations (see Rosotti 2023).We note, however, that the MRI may be fully suppressed with α ≤ 10 −4 (Bai 2013;Dullemond et al. 2022) and this may have important consequences for dust dynamics and growth (Vorobyov et al. 2018).On the other hand, spiral density waves that develop in young protoplanetary disks may stir up grown dust via the process known as gravitoturbulence (Riols & Latter 2018).If gravitoturbulence modifies the value of α, which in this case may become anisotropic (Baehr & Zhu 2021), the effect on dust growth and pebble formation can be substantial (Vorobyov et al. 2023a).
As the maximum size of dust grains a max begins to increase, part of the small dust population is converted to grown dust.This process is described by the conversion rate S (a max ), which can be expressed as where ∆ρ d,sm = ρ n+1 d,sm − ρ n d,sm is the mass of small dust (per unit volume) converted to grown dust during one hydrodynamic time step ∆t.The expression for ∆ρ d,sm is derived in Molyarova et al. (2021) and Vorobyov et al. (2022) assuming the conservation of dust mass and continuous dust size distribution across a * .It can be written as follows where the indices n and n + 1 denote the current and next hydrodynamic steps of integration, respectively.We note that the drift of grown dust with respect to small dust (recall that small dust is dynamically coupled with gas) can result in the development of a discontinuity in the distribution function N(a) at the interface a * between the small and grown dust populations.However, this may occur only if the drift timescales are shorter than the dust growth timescales, which is usually not realized in our models (see Fig. 8 in Vorobyov et al. 2022).
The dust growth in our model is terminated at the so-called fragmentation barrier (Birnstiel et al. 2012).We note that another important limiting mechanism -the drift barrier -is considered self-consistently via the solution of the dust dynamics equations.The maximum size up to which dust particles are allowed to grow due to the fragmentation barrier is defined as where u frag is a threshold value for the relative velocity of colliding dust grains, above which grain-to-grain collisions lead to destruction rather than to growth and Ω K is the Keplerian angular velocity.When a max exceeds a frag , the growth rate D is set equal to zero and a max is set equal to a frag .We note that if the local conditions in the disk change so that a frag drops below the current value of a max (e.g., when temperature increases or gas density decreases), then we also set a max = a frag .This effectively implies that part of the grown dust with sizes in the [a frag : a max ] range is shattered via collisions and the process of dust conversion reverses, namely, grown dust is now converted to small dust.A more sophisticated approach that takes into account the growth and fragmentation probabilities, such as in Akimkin et al. (2020), will be considered in a follow-up study.We note that the full solution of the Smoluchowski equation (e.g., Dr ążkowska et al. 2019) and is computationally expensive for three-dimensional numerical simulations.

Nested grid layout
Numerical simulations that follow self-consistently the transition from a collapsing cloud core to the stage of star and disk formation have to deal with vastly changing spatial scales.This can be achieved by using SPH codes (e.g., Bate 2022) or specifically designed grid-based codes.These latter use either the spherical coordinate system with radial and angular refinement toward the coordinate center and disk midplane (e.g., Meyer et al. 2017;Hosokawa et al. 2016;Oliva & Kuiper 2023), or adaptive meshes (Fromang et al. 2006;Hennebelle et al. 2020), or nested Cartesian grids (e.g., Machida et al. 2014;Tomida et al. 2017).The ngFEOSAD code adopted the latter approach.In comparison to spherical coordinate systems it does not suffer from the problem A202, page 4 of 20 Vorobyov, E. I., et al.: A&A, 683, A202 (2024) of singularity at the coordinate center and along the rotation axis, and is therefore more suitable for studying the formation of binary and multiple stars, which can freely move around and through the coordinate center.The angular momentum conservation test with the known analytic solution (Norman et al. 1980) demonstrate good performance, see Appendix C.4.2.We use 12 nested grids with the linear size of the outermost m = 1 grid equal to 0.09 pc.The number of grid cells per coordinate direction of each nested grid is N = 64.The effective numerical resolution in the inner 4.5 au is 0.14 au and it remains at a sub-au level up to 36 au.The forming disk in our runs is usually contained within 60-70 au.The vertical extent of the disk (from the midplane) is usually resolved by 6-7 grid zones up to a radial distance of 10 au and slightly deteriorates at larger radial distances.No mirror symmetry with respect to the z = 0 plane is imposed.

Numerical solution method
To solve the equation of hydrodynamics for the gas component, we employ the Godunov method with a combination of HLLC (Harten-Lax-van Leer contact) and HLL (Harten-Laxvan Leer) Riemann solvers (Toro 2019).We note that the gravity and friction forces are considered separately after the Godunov step.We also consider the equation for the total energy E in the test problems but in the actual simulations we employ the polytropic equation of state, which allows us to eliminate the energy equation.A more accurate approach requires also solving for the radiation transfer equation, which is planned as a future update to the ngFEOSAD code.
The equations of hydrodynamics for the gas component are written in the vector form where Q is the vector of conservative variables (ρ g , ρ d,sm , ρ g u, E) and F(Q) is the flux of conservative variables (ρ g u, ρ d,sm u, ρ g u ⊗ u + IP, (E + P)u).We note that here we also consider the continuity equation for small dust because it is dynamically coupled to gas.At each level of the nested grid, we use the classical finite-volume Godunov scheme to solve Eq. ( 17), in which we place the values of conservative variables at the center of the cells.The fluxes of conservative quantities at the interfaces between the cells are the solution to the Riemann problem.In the case of an interface between the nested grids, the solution to the Riemann problem is averaged for larger cells.For smaller cells, the averaging is not required.For the isothermal gas the HLLC scheme is replaced with the HLL scheme to avoid the singularity in the values of gas pressure between the limiting characteristics.The time step is calculated using the Courant-Friedrichs-Lewy condition with the Courant number set equal to 0.3.To increase the order of accuracy of the numerical method on smooth solutions and reduce dissipation at discontinuities, we use a piece-wise parabolic reconstruction of physical variables with integrals along the sound characteristics as parameters of the Riemann problem.A detailed description of the piecewise parabolic reconstruction, integration along characteristics, the choice of variables for reconstructions, and the rationale for using sound characteristics can be found in Kulikov & Vorobyov (2016); Kulikov et al. (2019); Kulikov (2022).
The dynamics equations of the grown dust component can be written in the vector form as follows where Q d is the vector of conservative variables (ρ d,gr , ρ d,gr a max , ρ d,gr u) and F(Q d ) is the flux of conservative variables (ρ d,gr u, ρ d,gr a max u, ρ d,gr u ⊗ u).We note that we introduced here an additional variable, namely, the product of ρ d,gr and a max .The reason for that is explained later in this section.The equations for the dust component are also solved by the Godunov method, similar to that used to solve the hydrodynamic equations.The key difference from the method of solving the gas hydrodynamic equations is the use of a one-wave HLL Riemann solver (details can be found in Appendix B) owing to the specifics of the pressureless dust fluid.
Once the Godunov steps for the gas and dust dynamics are accomplished, we update the gas and dust momenta due to the gravitational acceleration ∇Φ.The gradient of the potential is discretized using the standard eight-point stencil for the central difference scheme.The gravitational potential Φ is computed using the following triple sum where M(x i ′ , y j ′ , z k ′ ) is the mass of gas and dust in a numerical cell (i ′ , j ′ , k ′ ) with coordinates (x i ′ , y j ′ , z k ′ ), and G is the gravitational constant.To accelerate the summation over all numerical cells we employ the convolution theorem (Binney et al. 1987).
The extension of the convolution method to nested grids can be found in Vorobyov et al. (2023b), see also Appendix C.3.The update of the gas and dust velocities due to the friction force f between dust and gas is computed separately after the Godunov and gravity substeps using the following scheme (Lorén-Aguilar & Bate 2015b; Stoyanovskaya et al. 2018) where, ε = ρ d,gr /ρ g is the grown dust-to-gas ratio.Equations ( 20)-( 22) were derived assuming constant values of t stop and ε during a hydrodynamical time step ∆t, in which case an analytic solution for the updated gas and dust velocities is possible.The scheme performs well on the standard Sod shock tube and dusty wave test problems (see Appendix C.2).The numerical solution of Eq. ( 11) for the time evolution of the maximum dust size a max requires certain manipulations with its original form.In particular, the solution is split into two steps.First, the growth rate D is set to zero and the input from dust advection is computed.Then, the second term on the left-hand side is set to zero and the input from dust growth is considered.However, the second term on the left-hand side is not a divergence of a flux and, hence, cannot be easily incorporated into the Godunov conservative numerical scheme.
Nevertheless, Eq. ( 11) can still be recast in the conservative form with the help of the continuity Eq. ( 6).The resulting equation reads as follows The left-hand side is the continuity equation for the product ρ d,gr a max and we can use the same numerical method as for the A202, page 5 of 20 Vorobyov, E. I., et al.: A&A, 683, A202 (2024) grown dust density to solve it.The product ρ d,gr a max can then be updated to account for the combined source term shown on the right-hand side of Eq. ( 23) using the values of ρ d,gr and a max from the previous time step.The updated value of a max is finally calculated using the updated value of ρ d,gr obtained from Eq. ( 6).We checked and confirmed that a max remains unchanged when the growth term D is set to zero during the collapse and disk formation stage.In our design, each Coarray image, or separate program instance, performs calculations for a different subset of the domain, in our case different levels of the nested mesh.This setup means that computational tasks are efficiently spread across multiple multiprocessors given the balanced workload between images.It also means that memory transfers between distributed memory spaces are limited to only those necessary between nested meshes, something minimized by our mathematical methods for example in solving the gravitational potential (Vorobyov et al. 2023b), reducing one of the main computational overheads associated with distributed memory systems.Shared memory parallelism using OpenMP, providing each parallel thread with a single dedicated physical multiprocessor core, is then used to accelerate local computation.
Our models use 12 nested grids, each with a resolution of 64 3 corresponding to an effective minimum grid size of 0.14 au, deployed on the Vienna Scientific Cluster (VSC-5) which has two sockets per node with each socket connected to a multiprocessor with 64 physical cores.In our case, shared memory parallelism acceleration within a nested grid saturates beyond 32 cores, and so we assign two images per socket, resulting in four images per node and three nodes in total.This results in a very efficient, load-balanced and hardware-optimized domain decomposition.In such a configuration, simulations of the disk formation and evolution to the 18 kyr stage take only 12 days to perform and require only 384 cores, a very compact and still accelerated result.Description of the stepby-step solution procedure and more details on the parallel realization using the Coarray Fortran technique can be found in Appendix A.

Initial and boundary conditions
For the initial conditions we consider a Bonnor-Ebert sphere with mass M cl = 0.8 M ⊙ and radius R cl = 8 × 10 3 au.The initial temperature is set equal to 10 K.The initial dust-to-gas mass ratio is 0.01.To initiate the collapse we give the cloud an initial positive gas and dust density perturbation of 30%, namely, we increase the corresponding densities by 30%.The ratio of rotational to gravitational energy is β = 0.5% and the ratio of thermal to gravitational energy is 70%.Initially, most dust is in the form of small dust grains with a maximum size of 1 µm.For computational reasons, we also add a small amount of grown dust, about 10% of the total dust mass budget, which has a maximum size of a max = 2 µm.The cloud is initially submerged into an external medium with negligible density and velocity.The gas and dust dynamics is calculated for the entire computational domain, but the gas and dust velocities in the external medium are reset to zero after every time step.With this procedure, some mass and angular momentum may leak across the boundary between the Bonnor-Ebert sphere and the external medium.The magnitude of this effect does not exceed a few per cent of the total mass and angular momentum budget initially contained in the Bonnor-Ebert sphere and affects only the outer regions of the collapsing envelope.

Results
In this section, we provide the results of our numerical simulations, focusing on the efficiency of dust growth, pebble formation, and dust settling in a young protoplanetary disk.

Fiducial model. Gas distribution
Figure 1 presents the evolution of the gas number density n g in the fiducial model in the disk midplane for 18 kyr after disk formation.A distinct two-armed spiral pattern forms after ≈5.0 kyr.The calculation of the Toomre Q parameter confirms that the disk is gravitationally unstable with Q ≈ 1.The gas disk grows with time owing to accretion from the collapsing cloud and reaches a radius of ≈60 au by t = 18 kyr.The highest density is found in the disk center where the protostar is forming, though our numerical resolution is insufficient to follow the second Larson collapse to stellar densities.The distribution of n g as a function of distance from the disk center follows the n g ≃ r −2.57law.For a gravitationally unstable disk that settles in a state with Q ∼ 1, the n g ≃ r −3 law is expected.This mismatch is likely caused by the nonsteady character of the young disk and also by deviations from Keplerian rotation owing to a notable contribution from the massive disk to the stellar gravity.
Figure 2 displays the gas volume density distribution at the vertical cuts taken through the x − z plane at y = 0. Interestingly, as the disk forms and grows in size, the vertical shape of the disk changes.At t = 1.5 and 5.5 kyr the disk surface is rather smooth and the disk thickness monotonically increases with distance from the disk center.However, already at t ≥ 11.5 kyr the disk surface attains a wave-like shape.A comparison of Figs. 1 and 2 suggests that the disk bulges at the position of the spiral arms and shrinks in the inter-armed regions.
To investigate how the disk thickness behaves with radial distance from the disk center, we calculated the vertical scale height of the gas disk as where Σ g is the surface density of the gas disk (the details on how we calculate this quantity are provided later in this section) and ρ g (z = 0) is the gas volume density in the disk midplane.This equation holds in the assumption of the local vertical hydrostatic equilibrium, which may not be fulfilled in a gravitationally unstable disk.Moreover, Riemann solvers may require higher vertical resolution than currently achieved in our modeling to maintain the vertical equilibrium and dust settling may exacerbate the problem further.Therefore, our calculations of H g should be considered with caution, only as a first but still useful approximation.
Figure 3 presents the radial profile of H g taken along the xcoordinate with z = 0 and y = 0 at a time instance of t = 15 kyr.For convenience, the corresponding distribution of the gas number density n g is also shown.The forming star occupies the inner few astronomical units and H g is almost constant there.The value of H g starts growing as the protostar transits to the disk at r > 2-3 au.At r = 10 au, the vertical scale height is H g ≈ 1.3 au and at r = 50 au H g ≈ 3 au.Excluding the inner several au, H g scales as r 0.6 .The local peaks in n g that are seen at r ≈ 8 au and 30 au are the density enhancements in the spiral arms that cross the y = 0 line at these locations.The values of H g also have local maxima at the positions of the spiral arms.This analysis confirms that the disk bulges at the position of the spiral arms, at least in the polytropic approximation.This would make such spirals easier to detect observationally thanks to illumination by the stellar radiation.On the other hand, the more distant segment of the spiral arm may be hidden behind the closer segment, making the spiral structure more difficult to discern.Postpocessing with, for example, RADMC-3D is required to assess the detectability of the spirals in this case (e.g., Dullemond et al. 2012;Meyer et al. 2019;Cadman et al. 2020).

Fiducial model. Dust distribution and size
We focus now on the spatial distributions of the dust maximum size a max and the total dust-to-gas mass ratio in the disk midplane and across the disk vertical extent.Figure 4 gives a detailed view of the dust maximum size at t = 18 kyr after disk formation.The left-hand side subpanel in the first row and also the third row indicate that a max in the disk midplane increases sharply from ≈1 µm at the disk-envelope interface to a few millimeters in the disk outer regions, and then a max grows slowly to centimeters in the inner several tens of astronomical units.The maximum values can reach decimeters at about 1 au, but we note that this region is occupied by the forming star and should be considered with care.The dominance of millimeter-and centimeter-sized grains imply that the spiral pattern may be detectable already at this yearly stage (Cadman et al. 2020).We do not find any notable dust growth in the bulk of the infalling envelope because A202, page 7 of 20  of its rarefied density and low turbulent velocity, see Eq. ( 12).We acknowledge that our dust growth model is designed for protoplanetary disks and not for infalling envelopes, for which certain modifications, for example, to the Stokes number and turbulent α parameter are required (Commerçon et al. 2023).This is beyond the scope of the current paper.
To estimate the efficiency of dust growth, we calculate the ratio of the fragmentation barrier a frag to the maximum dust size a max in the disk midplane.The right-hand side subpanel in the first row of Fig. 4 demonstrates that dust growth has almost saturated in the inner disk where the ratio approaches unity (note that the scale bar is in log units).There is also a region of the lower spiral arm where dust growth is saturated but this is due to low a frag in this rarefied region.In general, however, the spiral arms are the disk regions where dust has not yet reached its upper values determined by collisional fragmentation and there is a potential for further dust growth there.
The distribution of a max in the vertical cut taken along y = 0 (second row in Fig. 4) is characterized by higher values closer to the disk midplane.This can be explained not only by dust settling but also by a more favorable dust growth environment in the disk midplane characterized by higher dust densities.
A202, page 8 of 20 Finally, the bottom row in Fig. 4 presents the Stokes number as a function of radial distance from the forming protostar in the disk midplane.The Stokes number gradually increases outward.Interestingly, St reaches a value of 0.01 only beyond 10 au, which has important consequences for the spatial distribution of pebbles as discussed later in the text.Another important feature is that St never reaches values close to 1.0, as often assumed in dust dynamics simulations with a fixed dust size without growth (e.g., Rice et al. 2004;Boss et al. 2020).These works considered preset isolated protoplanetary disks, more appropriate for the T Tauri phase, and found efficient dust accumulation in the spiral arms for dust grains with St ∼ 1.0 or dust size 1-10 m.Our modeling indicates that such values of the Stokes number and dust size are not easily realizable in more realistic young protoplanetary disks formed via gravitational collapse, at least in the early 20 kyr of evolution.
This difference can be understood if we compare the gas density and temperature of our young disk with those of the Minimum Mass Solar Nebular (MMSN).Figure 5 presents the radial profiles of the gas surface density and temperature at 18 kyr.The surface density profile of the MMSN is taken from (Hayashi 1981) and has the following form Clearly, Σ g is about an order of magnitude higher than Σ MMSN .
The bottom panel indicates that the water snow line in our disk, taken to be at a radial position where the gas temperature equals 150 K, is beyond 10 au.This also indicates that our disk is warmer than that of the MMSN.Finally, we note that the protostellar mass at t = 18 kyr is lower by a factor of 10 than the Sun.All these factors taken together lower the Stokes number for cm-sized dust particles below 0.1 in a young protoplanetary disk.
Figure 6 presents the total dust-to-gas mass ratio at t = 18 kyr.The top panel displays the corresponding values calculated in the disk midplane as The value of ξ d2g (z = 0) throughout most of the disk midplane is elevated and can become as high as 0.1.At the same time, the spatial distribution of ξ d2g (z = 0) is highly spatially nonhomogeneous, which was also found to be the case for gravitationally unstable disks in our earlier thin-disk simulations (Vorobyov & Elbakyan 2019).The highest values of ξ d2g (z = 0) are found in between 1.0 and 10 au, and also in the outer disk parts.The latter regions, however, are characterized by quite low values of the gas and dust densities, and may be unfavorable from the point of view of planet formation.The middle panel in Fig. 6 presents the dust-to-gas ratio calculated across the disk vertical extent at y = 0.The plot reveals a strong vertical segregation in the dust-to-gas ratio, with the highest values found near the midplane and lowest in the upper/lower disk atmosphere.The disk midplane is enhanced in dust relative to the fiducial value of ξ d2g = 0.01, while the disk atmosphere is strongly depleted in dust.These variations are caused by dust growth and settling toward the disk midplane.Interestingly, there is a mild enhancement in ξ d2g just above/below the disk surface.This enhancement is caused by grown dust falling in faster than gas in the envelope (recall that the prestellar core in our setup contained a small amount of grown dust with a max = 2 µm, see Sect.2.7).Dust is pressureless and falls in on a free-fall time scale.Gas falls in slower because it has substantial pressure support against gravity, with the ratio of thermal to gravitational energy set equal to 70% in the prestellar core, much higher than the corresponding ratio for rotational energy, 0.5%.
The bottom panel presents the values of ξ d2g (z = 0) as a function of the radial distance from the disk center.The plot confirms that the dust-to-gas ratio in the disk midplane is enhanced through most of the disk extent, except for the outermost regions where its value may drop below 0.01.The peak values approaching ξ d2g (z = 0) = 0.1 are found at radial distances from several to ten astronomical units.This dust concentration is also revealed as a dust ring in the top panel of Fig. 6.
To probe if spiral arms can efficiently accumulate dust in our model, we take a disk snapshot at t = 15 kyr with a well developed two-armed spiral pattern.Spiral arms are the disk regions with the enhanced gas density, and we may expect a positive correlation between ξ d2g and n g in the disk midplane if the accumulation of dust is present in the arms.We cut a radial region between 30 and 50 au to remove disk regions where spiral arms are less prominent.Figure 7 displays the corresponding relation in the disk midplane and the red line plots the power-law best fit to the data of the form The positive correlation takes place but is rather weak.The lack of strong dust concentration toward spiral arms was also reported in Vorobyov et al. (2018Vorobyov et al. ( , 2022)), and hence this phenomenon is not an artifact of the thin-disk simulations as was suggested in Boss et al. (2020).We think that the dust concentration in spiral arms is hindered by a rather low Stokes number in young   disks.The nonsteady character of the spiral pattern, with the arms changing its form on timescales of a few kyr, may also work against strong dust concentrations.Our modeling emphasizes the need for realistic disk formation and dust growth models to consider dust accumulation in spiral arms.

Fiducial model. Pebbles
It is interesting to consider the distribution and time evolution of pebbles in the disk, since they may play a crucial role in the growth of protoplanetary cores (Lambrechts & Johansen 2012;Ida et al. 2016;Johansen & Lambrechts 2017).Following Vorobyov et al. (2023a), we define pebbles as dust grains with radius greater than or equal to a 0 peb = 0.5 mm and St ≥ St 0 = 0.01.More specifically, we calculate the minimum pebble size as where a 0 St is the size of dust grains that have the Stokes number equal to St 0 for the local disk conditions, a 0 St = a max St 0 /St.We then search the disk locations where these conditions on the dust size and Stokes number are fulfilled.Noting that pebbles share the same size distribution with grown dust, the mass of pebbles in each computational cell can be calculated knowing the mass of grown dust, the slope of the dust size distribution (fixed at p = 3.5 in our work), the minimum and maximum sizes of grown dust, a * and a max , respectively, and the minimum size of pebbles, a min peb (see Vorobyov et al. 2023a, for details).Figure 8 shows the spatial distribution of pebbles at t = 18 kyr.Clearly, pebbles are distributed unevenly across the disk midplane.Such a patchy spatial distribution of pebbles is also A202, page 10 of 20 Vorobyov, E. I., et al.: A&A, 683, A202 (2024) Fig. 9. Time evolution of the total gas, total dust, grown and small dust, and pebble masses in the disk.Note that the dust and pebble masses are scaled up by a factor of 100.
found at earlier evolution times.There are regions where pebbles are entirely absent, and most interestingly, in the inner 10 au.Dust can grow to centimeters there but the Stokes number happens to be lower than 0.01 because of the rising gas temperature and density in the innermost disk regions, see Eq. ( 10).This finding is in agreement with our previous work using thin-disk simulations (Vorobyov et al. 2023a), but we note that decreasing the turbulent α from 10 −3 (as in the current work) to 10 −4 may have a strong impact on the pebble distribution in the disk.We also found that pebbles are concentrated to the midplane, which is a consequence of dust settling.
Figure 9 presents the time evolution of the total pebble mass in the disk as it forms and evolves.The total masses of gas and dust (both small and grown) in the disk are also shown for comparison.We note here that we do not distinguish between the disk and the forming protostar, so that the mass of gas and dust include both.On the other hand, pebbles are absent in the inner several astronomical units and the mass of pebbles is that contained in the disk only.We confirm with 3D simulations that dust growth proceeds very fast, a conclusion also made in our earlier thin-disk simulations (Vorobyov et al. 2018).The first generation of pebbles appears after 2.5 kyr from the disk formation instance and the mass of pebbles reaches 31.7 M ⊕ by the end of our simulations.It can be noted that the mass of pebbles is not growing steadily but there are episodes when its total mass decreases.This is caused by the varying local disk conditions in which pebbles find themselves as they drift through the disk.In particular, the Stokes number may decrease below the critical value of 0.01 because of the locally rising gas density and/or temperature.Also, the maximum size of dust grains may decrease if the local fragmentation barrier decreases.Overall, however, the pebble mass increases with time and appears to saturate at ≈31.7 M ⊕ .Given the uncertainties with the definition of the Stokes number in young disks, we also investigated how the pebble mass would change if St 0 is varied by a factor of two.The corresponding mass of pebbles at t = 18 kyr decreased/increased approximately by the same factor.Finally, we note that the mass of small dust is much lower than that of grown dust throughout most of the considered disk evolution period, and by the end of simulations becomes lower than that of pebbles.
To assess the efficiency of grown dust settling, we plot the radial distribution of the following ratio This equation is derived assuming Gaussian distributions of gas and dust volume densities in the vertical direction, which is an acceptable approximation, and assuming that small dust has the same vertical scale height as that of gas.To calculate the disk surface density, we adopted the disk tracking conditions outlined in Joos et al. (2012).In particular, we use the following criteria to determine if a particular computational cell belongs to the disk: -the gas rotational velocity must be faster than twice the radial infall velocity, v ϕ > 2v r , -the gas rotational velocity must be faster than twice the vertical infall velocity, v ϕ > 2v z , -gas must not be thermally supported, ρv 2 ϕ 2 > 2P, -the gas number density must be higher than 10 9 cm −3 .If any of the first three conditions fail, a particular grid cell is not qualified as belonging to the disk.The forth condition must always be fulfilled.
Figure 10 displays the ratio of the scale heights H d /H g at t = 18 kyr.Grown dust settling is present across the entire disk extent, being strongest in the inner 10 au, where H d /H g ≃ 0.2.At larger distances, however, dust settling weakens to H d /H g ≃ 0.4-0.5, most likely owing to the perturbing influence of spiral density waves, as was shown in Riols et al. (2020).Our values (0.4-0.5) are still lower than those quoted by these authors (≥0.6), which can be attributed to different simulation techniques.We model the development of spiral arms in global disk formation simulations, while Riols et al. (2020) do that for the closed-box simulations.The regions with H d /H g > 1.0 are caused by strong perturbations in the flow at the envelope to disk transition, which is a transitory effect.

Effects of lower β
We have also considered the collapse of a cloud core with half the rotation velocity as compared to the fiducial model.The disk in this model forms about 3 kyr later, but in general the disk evolution in both models follows a similar path.In Figs.11 and 12 we show the time evolution of the disk and pebble masses, as well as the radial distribution of the H d /H g ratio.The mass of pebbles is a bit lower at the same time instance after the disk formation, t = 18 kyr, and is equal to 18.3 M ⊕ .The time evolution of the total pebble mass is also nonmonotonic, which is a consequence of the spatially patchy distribution of pebbles across the A202, page 11 of 20

Model caveats and future developments
This is our first step toward fully three-dimensional numerical hydrodynamics modeling of the formation of protoplanetary gasdust disks, which builds upon our previous extensive studies of long-term gas-dust disk evolution using the thin-disk approximation (e.g., Vorobyov et al. 2018Vorobyov et al. , 2023a;;Molyarova et al. 2021).Further developments of the ngFEOSAD code are needed to improve the realism of simulations.The barotropic approach is to be replaced with more realistic radiative transfer simulations (e.g., Flock et al. 2017) and magnetic fields have to be introduced to better describe the disk dynamics and flow in the envelope and the disk (e.g., Tsukamoto et al. 2020;Lebreuilly et al. 2020;Koga & Machida 2023).The simplified calculations of the thermal balance may affect the dust growth rate and the fragmentation velocity in our models, which in turn may have an impact on the spatial distribution and mass of pebbles in the disk.The lack of magnetic field may affect the dynamics of both gas and dust, with nonmagnetized disks being more massive and extended compared to their magnetized counterpart.Dust growth is sensitive to ice mantles (Okuzumi et al. 2016) and the dynamics of volatiles is to be considered together with dust evolution to accurately capture this effect (Molyarova et al. 2021).A careful consideration to the sink cell model as a proxy of the forming star has to be given to avoid possible numerical artifacts (Machida et al. 2014;Vorobyov et al. 2019;Hennebelle et al. 2020, see also Appendix C.4).The dust fragmentation barrier, dust settling, and dust growth depend on the α parameter, and choosing a value other than 10 −3 (adopted here) may have a strong impact on the dust and pebble distribution in the disk.Moreover, protoplanetary disks are unlikely to have spatially and temporally constant values of the turbulent α parameter and a more accurate model that considers, for example, luminosity bursts triggered by the magnetorotational instability is to be developed (e.g., Bae et al. 2014;Vorobyov et al. 2020).
Better dust growth models that include more dust size bins and finally transit to the solution of Smoluchowski equation are needed (e.g., Dr ążkowska et al. 2019), although this might be very challenging in 3D simulations from the numerical point of view.The use of two bins forced us to consider the dynamics of dust with the maximum size because they are the main dust mass carriers.However, the dynamics of dust grains with smaller size is also important as they may substantially contribute to the dust opacity and, hence, to the disk thermal balance.Finally, we plan to perform higher resolution simulations with the base number of grid cells N = 128 3 when the code is updated to use the GPUs for the gravitational potential solver.The current resolution with N = 64 3 is insufficient to fully resolve the vertical extent of the gas disk, not to mention that of the grown dust given that dust settling is present, which may affect our conclusions on the dust settling efficiency.

Conclusions
We have developed an original three-dimensional numerical hydrodynamics code on nested meshes, ngFEOSAD, which computes the coevolution of gas and dust in a protoplanetary disk, which in turn forms self-consistently via the gravitational collapse of a rotating cloud.We employed a novel Coarray Fortran, combined with OpenMP, parallelization scheme, allowing us to efficiently use compact high-performance computing resources.The dust dynamics and growth, including pebble formation, were studied in detail during the initial 18 kyr of disk evolution.Our results can be summarized as follows.
The disk becomes gravitationally unstable and develops a spiral pattern already after several thousand years from its formation instance.The radial distribution of the gas volume density in the disk midplane is approximately proportional to r −2.5 , which is typical of gravitationally unstable disks.In the vertical extent, the disk bulges at the position of the spiral arms.
The dust-to-gas mass ratio in the disk midplane is highly nonhomogeneous throughout the disk extent, in agreement with our earlier thin-disk simulations (Vorobyov & Elbakyan 2019) and Lebreuilly et al. (2020).The dust-to-gas mass ratio also shows prominent vertical stratification, indicative of dust vertical settling.The ratio of the vertical scale heights of the dust and gas disks was found to be in the 0.2-0.5 range.Such a moderate dust settling can be attributed to perturbing action of the spiral arms, as was already noted by Riols et al. (2020).
Dust grows fast in a young protoplanetary disk and reaches sizes on the order of 1-10 mm already after 15 kyr.However, the Stokes number remains rather low, St ∼ 10 −3 −10 −1 throughout the disk extent, with the largest values closer to the disk outer edge.These low St are caused by systematically higher densities and temperatures in young disks around subsolar mass stars A202, page 12 of 20 Vorobyov, E. I., et al.: A&A, 683, A202 (2024) as compared to the corresponding characteristics in the MMSN.Low Stokes, along with a nonsteady character of the spiral pattern, hinder strong dust accumulation in the spiral arms, although in general the dust-to-gas ratio in the midplane is enhanced by a factor of several compared to the initial 1:100 value.Our modeling demonstrates the importance of coupled disk formation and dust growth calculations in realistically assessing the efficiency of spiral arms as dust traps in comparison to models of isolated disks with a fixed dust size or Stokes number (e.g., Rice et al. 2004;Boss et al. 2020).
The spatial distribution of pebbles exhibits a highly nonhomogeneous and patchy character.Interestingly, pebbles are absent in the inner 10 au because of the low St < 10 −2 in these regions.Pebbles are also found to be strongly concentrated to the disk midplane.The total mass of pebbles in the disk increases with time and reaches a few tens of Earth masses by 18 kyr.
Variations in the initial rotation rate of the collapsing cloud were found to have minor effect on our conclusions.We note that when calculating dust growth we assumed a spatially and temporally constant value of the turbulent α parameter.We plan to explore the effects of this simplification in follow-up works.

Nested Meshes
Level 1. Projection.This first step implements the averaging of the conserved hydrodynamic quantities (i.e., densities and momenta for both gas and dust) from a finer to a coarser grid to ensure their conservation on nested grids.The recalculation scheme is illustrated in Figure A.2.The conservative values are averaged over four (or eight in the three-dimensional case) cells using conservation laws.An example of the CAF code with one-sided communications is shown, which transfers part of the array to a coarser grid.We note that the recalculation of the conservative values is carried out sequentially from a finer to a coarser grid level (not in parallel mode).While this removes the CAF (distributed memory) parallelization, it does not affect OpenMP (shared memory) parallelization.In any way, this procedure is computationally inexpensive.2. Boundary.Here, the boundary conditions are implemented by adding one outer layer of cells to each level of the nested grids.For the outer coarsest grid, the boundary conditions are determined from the physics of the problem.On all higher-level grids, the boundary values are taken from the corresponding cells on the underlying coarser grid.The gravitational potential is linearly interpolated to improve the accuracy (see Vorobyov et al. 2023b).This step is implemented in a parallel mode for both CAF and OpenMP.3. Primitive.This step finds the physical variables (e.g., densities and velocities of both gas and dust) from conservative variables.This procedure is implemented in parallel (Carray Fortran and OpenMP) for each level of the nested grids without data exchanges between the grids.4. Courant.Here, the time step is calculated, which is the same for all nested grid levels.Since each grid level is located in its own CAF image, we use the CAF reduction procedure call co_max to find the maximum wave speed over all grid levels, which is then used to determine the minimum time step applied to all nested grids. 5. Reconstruction.This step computes the piece-wise parabolic reconstruction of physical variables in all cells and in all coordinate directions at each grid level.Both CAF and OpenMP parallelizations are employed.6. Godunov.Here we solve the Riemann problem for each level of the nested grids.At the boundary with a coarser grid, four cells participate in the solution of the Riemann problem with the same neighboring cell on the coarser grid.This step employs both CAF and OpenMP parallelization.7. Gravity.This step accounts for the effects of disk selfgravity and the gravity of the star(s) (if formed).We note that in the present work we do not introduce a sink particles as proxy for the star but plan to do that in the future.The self-gravity of both gas and dust is considered.The detailed explanation of the applied method can be found in (Vorobyov et al. 2023b).The details of the CAF-OpenMP realization will be presented elsewhere.8. Friction.The computation cycle is completed with the calculation of the friction force between the dust and gas and its effect on the gas and dust momenta.This step employs the full power of both CAF and OpenMP parallelization.

Appendix B: Riemann solver for grown dust
To find the fluxes of conservative variables for the gas component (see Eq. 17), we use a combination of the HLLC and HLL Riemann solvers, the description of which can be found in Toro (2019).To find the fluxes for the pressureless grown dust (see Eq. 18), we use a modified HLL solver.To describe the scheme for solving the Riemann problem for grown dust, we consider a configuration of the left (L) and right (R) cells, in which the values of ρ L , ρ R are the dust densities, u L , u R are the dust longitudinal velocity, v L , v R , w L , w R are two transverse dust propagation velocities.The one-wave HLL algorithm for grown dust is then as follows.
1. We form vectors of conservative variables and their fluxes: 2. We calculate the exact value of the dust transport rate according to the following equation: A202, page 14 of 20 Vorobyov, E. I., et al.: A&A, 683, A202 (2024) 3. We write down a one-wave HLL scheme in the following basic form: Similarly to the solution of the hydrodynamic equations for gas, a piece-wise parabolic representation of physical variables is used to increase the order of accuracy on smooth solutions and to reduce dissipation at discontinuities.In this case, instead of the piece-wise constant values of the physical quantities q L and q R , where q is understood as the density and velocity functions, an integral over the parabola constructed along the characteristic S is considered.Kelvin-Helmholtz instability can be estimated as: Rayleigh-Taylor instability is determined as: where A = 2/3 is Atwood's number.For example, at time t = 13, the amplitude is f (t = 13) ≃ 0.25, which is in good agreement with the model results.Indeed, the initial position of the interface between different gas states is located at y = 0.75.At t = 13, the density perturbation propagates approximately dy = ±[0.25 : 0.3] from the initial position of the interface between the two states.
Appendix C.1.3.The Sedov point explosion test Traditionally, the Sedov test assesses the ability of a numerical method to reproduce flows with high Mach numbers, which may also occur in star formation.We set the computational domain [−0.5; 0.5] 3 , the adiabatic exponent γ = 5/3, the initial density in the domain ρ 0 = 1, and the initial pressure p 0 = 10 −5 .At time t = 0, the internal energy e 0 = 0.6 is released.The explosion area is limited by the radius r 0 = 0.01.The computational grid with 100 3 cells is used.The In our model test, the sound speed of the background medium is negligibly small, so the Mach number reaches M ∼ 2000.ngFEOSAD quite well reproduces the position of the shock wave, as well as the density and momentum profiles.
Appendix C.2. Tests on the dust dynamics scheme and friction force In this section, we tested the dust solver on the problems that consider a pure dust fluid and a mixture of gas and dust, including the friction force between the two components.We have considered two numerical schemes to compute the momentum exchange between the gas and dust components, including the backreaction of dust on gas.The first method makes use of an analytic solution for the updated gas and dust velocities as described in (see, e.g., Lorén-Aguilar & Bate 2015a; Stoyanovskaya et al. 2018).The second method employs a fully implicit scheme to compute the updated gas and dust velocities as presented in Stoyanovskaya et al. (2018).Both schemes perform well on the standard tests, such as the dusty wave and the Sod shock tube problems.Here, we present the results for the method that uses the analytic solution, which is also implemented in the ngFEOSAD code.
Appendix C.2.1.Sod shock tube for a mixture of gas and dust This test is used to assess the ability of a numerical algorithm to accurately track the position of shock waves and contact discontinuities when for a mixture of gas and pressureless dust.Initial conditions involve two discontinuous states, with a hot dense gas on the left and cold rarified gas on the right.In particular, we set the pressure and density of gas at x ∈ [0.0.5] to 1.0, while at x ∈ [0.5.1.0]the gas pressure is 0.1 and gas density is 0.125.The velocity of the γ=1.4 gas is initially zero everywhere.The dust component has the same initial distribution as that of gas, but the dust pressure is set to zero.The dust-to-gas ratio is, therefore, equal to unity everywhere.The analytic solution for the gas and dust mixture is only known in the limit of short stopping times compared to the time of shock wave propagation.We used the SPLASH code (Price 2007) to generate the analytic solution.We follow the terminology of Laibe & Price (2012) and introduce the drag parameter K and define the stopping time as t stop = ρ g /K, where ρ g = 1 is the gas density in left half of the shock tube.
Figure C.7 present the test results for K = 5000, namely, for a strongly coupled gas-dust mixture.The first and second rows present the resulting gas distributions at t=0.24, while the bottom row shows those of dust.The numerical resolution is 200 grid zones and the Courant number is set to 0.2.Clearly, the As a first test, we consider the collapse of a cold gaseous sphere.
We set a computational domain with a size of 0.15 pc in each direction and introduce a homogeneous gaseous sphere with radius 0.0652 pc.Initially, the number density of the sphere is set equal to 10 4 cm −3 .The rest of the computational domain is filled with rarefied gas of negligible density.The gas temperature is set to 0.1 K to imitate the cold sphere collapse.The mass of the cloud is 0.52 M ⊙ .We have chosen m = 9 nested grids with the number of grid cells in each direction N = 64, which corresponds to a effective minimal resolution of 1.83 au.
The results are shown in Fig. C.10.The code reproduces well the analytic solution in the inner regions, only slightly overestimating the analytic value at 99.6% of the free fall time.In the outer regions, an anomalous peak develops, which is likely caused by degrading spatial resolution at large distances from the coordinate center on nested grids.
Appendix C.4.2.Collapse of a rotating sphere In the absence of any mechanisms for angular momentum redistribution, such as during the gravitational collapse of an axisymmetric cloud, the mass with specific angular momentum less than or equal to J = Ωr 2 should be conserved and equal to (Norman et al. 1980) (C.10) Here, Ω is the angular velocity of the collapsing cloud at a distance r from the z-axis, and M 0 , Ω 0 , and R 0 are the initial mass, angular velocity, and radius of the cloud.We use this property to assess the code's ability to conserve angular momentum during the gravitational collapse.
The setup consists of a homogeneous gaseous sphere with the initial mass 1.4 M ⊙ and the ratio of rotational to gravita- tional energy set equal to 1.0%, which is typical for prestellar clouds (Caselli et al. 2002).The initial radius of the cloud, the number density, and the angular velocity are R 0 = 10 4 au, n 0 = 5 × 10 4 cm −3 , and Ω 0 = 4 × 10 −14 s −1 , respectively.Since we use the Cartesian grid, the cloud is initially submerged into a low-density medium.This causes the outer cloud layers to expand rather than collapse.This is undesirable as the analytic solution stands only for the pure collapse when the spherical layers fall in but do not mix with each other.To reduce this effect we set the cloud isothermal temperature to 0.1 K.We also found that the analytic solution holds only on regular grids.Therefore, we use one nested grid with numerical resolution 256 3 .
The test results are presented in Figure C.11 at t = 96.7% and t = 99.3% that of the free-fall time, namely, when the central density has increased more than by three and four orders of magnitude, respectively.The angular momentum for the bulk of the collapsing cloud is well conserved.The horizontal arrangement of the model data at low J, which can be seen even at the onset of collapse, is caused by insufficient numerical resolution and is not a signature of angular momentum nonconservation.A higher resolution would make these data points to fit better the red line.The deviation from the analytic solution occurs only at the smallest J near the rotation axis and high above and below the midplane of the collapsing cloud.These deviations are caused by small diffusion of the outer cloud layers into the external rarefied medium, which is difficult to avoid.In any case, the mass of the cloud material that shows the angular momentum A202, page 19 of 20 Vorobyov, E. I., et al.: A&A, 683, A202 (2024) nonconservation is only a minor fraction of the entire cloud mass, namely, less than 0.04%.When A = 2, the gas density distribution turns into a marginally stable singular isothermal sphere (Shu 1977).When A > 2, the system is out of virial equilibrium and begins to collapse.
In practice, we initialize the initial gas density distribution that is characteristic of the Bonnor-Ebert sphere with temperature 10 K, central density n BE,0 = 4.2 × 10 5 cm −3 and ratio of the central to the outer densities set equal to 20.The resulting radius of the Bonnor-Ebert sphere is ≈ 8000 au and the mass is 0.9 M ⊙ .The initial configuration is given a positive density perturbation of 30% to initiate the gravitational collapse.
The accretion on the protostar is realized by setting a sink sphere in the center of the collapsing cloud with an accretion radius of 3 au.If the gas density in the ith computational cell within the sink sphere exceeds ρ crit = m H µ n crit , where n crit = 10 13 cm −3 , the gas density is artificially decreased and the excess mass is transferred to the star positioned in the center of the sink sphere.The following equation describes the process Here, ρ 0,i is the gas volume density at the current time (at which the condition on the volume density is checked), ∆t is the time step, and b is the rate of mass evacuation from the sink sphere set equal to the local Keplerian angular velocity Ω K .This definition of the sink sphere is similar to that of Machida et al. (2014) but we use an analytic solution to this problem.The resulting protostellar mass accretion rate is defined as where ∆V i is the volume of the ith computational cell and the summation is performed over all cells within the sink sphere for which ρ 0,i > ρ crit .The top panel in Figure C.12 presents the number density distributions of the collapsing cloud at different times.We used m = 12 nested grids with N=64 grid cells in each coordinate direction.As the collapse proceeds, the density distribution outside of the central plateau acquires the form similar to Eq. (C.12), as can be seen from the red line showing the density profile given by Eq. (C.12) with A = 5.0.The bottom panel displays the resulting mass accretion rate on to the star vs. time.The star forms at about t ≈ 0.1589 Myr, which is manifested by a sharp increase in the mass accretion rate about this time instance.Such a peak in Ṁ is also found in other collapse simulations (e.g., Vorobyov & Basu 2005).Soon the mass accretion rate drops and stabilizes at ≈ 2.5 × 10 −5 M ⊙ yr −1 , which is only slightly higher than the analytically predicted value of 1.5 × 10 −5 M ⊙ yr −1 for A = 5.We attribute difference to deviations of the Bonnor-Ebert sphere from the form given by Eq. (C.12).

Fig. 1 .
Fig. 1.Number density for different time instances in the inner 100×100 au 2 computational domain.The time is shown since the instance of disk formation.The inserts show the radial distribution of the gas number density for each time instance on the log scale.

Fig. 2 .
Fig. 2. Vertical cuts showing the gas number density in the x − z at y = 0.The same time instances as in Fig. 1 are considered.

Fig. 3 .
Fig. 3. Disk vertical scale height H g (blue circles) and the gas number density in the midplane n gas (red crosses) as a function of distance in both directions from the disk center along the x-coordinate.

Fig. 4 .
Fig.4.Spatial and radial distribution of the maximum dust size a max and the Stokes number at t = 18 kyr.The top row consists of two subpanels: the left-hand side subpanel showing a max and right-hand side one presenting the ratio of the fragmentation barrier a frag to the actual maximum size of dust grains a max , both in the disk midplane.The second row is the vertical cut taken across the disk at y = 0, third and bottom rows are a max and St, respectively, as a function of radial distance in the disk midplane in all directions.The scale bars for a max are in log cm, while the color bar for a frag /a max is in log dimensionless units.

Fig. 5 .
Fig. 5. Comparison of the gas surface density and temperature in our model with those of the MMSN.The blue dots in the top and bottom panels are the model gas surface densities and temperatures at t = 18 kyr, respectively.The red line in the top panel shows the surface density profile for the MMSN.The dashed horizontal line in the bottom panel marks the gas temperature of 150 K, taken to be the water sublimation temperature.

Fig. 6 .
Fig.6.Total dust-to-gas mass ratio at t = 18 kyr.The top panel shows the two-dimensional distribution in the disk midplane, while the middle panel presents the data across the vertical cut taken at y = 0.The bottom panel is the dust-to-gas ratio as a function of radial distance in the disk midplane in all directions.The color bars in the top and middle panels are in log units.

Fig. 7 .
Fig. 7. Total dust-to-gas mass ratio ξ d2g as a function of the gas number density n g .All data are for the disk midplane and for a radial annulus 30-50 au capturing the most prominent spiral arm.The red line is the power-law best fit to the model data.

Fig. 8 .
Fig. 8. Spatial distribution of pebbles in the disk at t = 18.1 kyr.The top panel presents the number density of pebbles in the disk midplane, while the bottom panel shows the corresponding values across the disk vertical cut taken at y = 0.

Fig. 10 .
Fig. 10.Dust settling efficiency as a function of radial distance.Shown is the ratio of the grown dust-to-gas scale heights H d /H g at t = 18.1 kyr.

Fig. 11 .
Fig. 11.Similar to Fig. 9 but for the model with half the rotation velocity of the cloud core.

Fig. 12 .
Fig. 12.Similar to Fig. 10 but for the model with half the rotation velocity of the cloud core.

Fig
Fig. A.1.Decomposition of the computational domain.Each level of nested grids is placed on a separate CAF image.As an example, nested grids with three levels of refinement are shown.The number of coarray images is equal to the total number of nested grids.
Fig. A.2. Averaging of conservative hydrodynamic quantities from a finer to a coarser grid.The colors represent different values of conserved quantities.Only inner part of the coarser grid is affected by the averaging.
Fig. C.4.Development of the Kelvin-Helmholtz instability.Panel (a) presents the initial state, panel (b)t = 0.2, panel (c)t = 0.5, and panel (d)t = 0.6.The color bar shows the gas density in dimensionless units.

Fig
Fig. C.5.Development of the Rayleigh-Taylor instability.Panel (a) presents the initial state, panel (b)t = 72, panel (c)t = 10, and panel (d)t = 13.The color bar shows the gas density in dimensionless units.

Fig
Fig. C.6.Density (left) and momentum (right) in the Sedov blast wave problem.The solid lines represent the exact solution, while the symbols correspond to the model data.

Fig
Fig. C.10.Gravitational collapse of a cold gas sphere.The colored circles correspond to different times in terms of the free-fall time.The green lines present the corresponding analytic solutions.

Fig
Fig. C.11. Mass of the cloud with specific angular momentum less or equal to J. The top panel presents the M(J) − J relation at the onset of gravitational collapse, while the bottom panel present the corresponding relation at different times after the onset of collapse.The red line is the analytic solution and the symbols present the numerical data.
Appendix C.4.3.Protostellar accretion in spherical collapse problemsThe accretion rate of a spherically collapsing cloud on to a protostar has an analytic solution of the following form(Shu 1977 (ρ i − ρ crit ) , (C.13) which has an analytical solution of the following form ρ i = ρ 0,i exp −b ∆t +ρ crit 1.0 − exp −b ∆t .(C.14)

Fig
Fig. C.12. Collapse of a gaseous sphere followed by star formation and protostellar accretion.The top panel shows the radial distributions of the gas number density at different time instances during the collapse and after the formation of the star.The red line is the fit to Eq. (C.12) with A = 5.The bottom panel displays the resulting mass accretion on to the star.The time is counted from the onset of simulations.