© The Authors 2024

1 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), and HL 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 main-sequence 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.

This observational evidence motivated recent studies of dust dynamics and growth in numerical hydrodynamics models that simulate the transition from a collapsing prestellar core to the star and disk formation stage self-consistently. Zhao et al. (2018) suggest that dust growth and the removal of the small dust grains from the Mathis, Rumpl, and Nordsiek (MRN) dust size distribution changes the magnetic resistivity significantly and promotes disk formation in a magnetized environment. Magnetohydrody-namics simulations of disk formation by Tsukamoto et al. (2020) demonstrate that protoplanetary disks form when dust physics is included in the calculation of resistivity regardless of the dust size distribution. Bate & Lorén-Aguilar (2017) and Koga & Machida (2023) found that dust grains with sizes greater than 10 µm can dynamically decouple from gas already during the cloud core collapse. Interestingly, dust grains with sizes ≫1 µm 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 smoothed-particle 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 precom-puted 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, 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 par-allelization. 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.

2 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 Coar-ray 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.

2.1 Gas component

The dynamics of gas is followed by solving the equations of continuity and momentum $$\frac{\partial {\rho }_{\text{g}}}{\partial t}+\nabla \cdot \left({\rho }_{\text{g}}\mathit{v}\right)=0,$$(1) $$\frac{\partial }{\partial t}\left({\rho }_{\text{g}}v\right)+\nabla \cdot \left({\rho }_{\text{g}}v\otimes v+\mathbb{I}P\right)=-{\rho }_{\text{g}}\nabla \text{Φ}-{\rho }_{\text{d},\text{gr}}f,$$(2)

where ρ_g is the volume gas density, υ the gas velocity, Φ the gravitational potential, P the gas pressure, 𝕀? 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 $$P_k=c_{\text{s},0}^2{\rho }_{\text{g}}^{{\gamma }_k}{\displaystyle \prod _{i=1}^{k-1}{\rho }_{\text{c},i}^{{\gamma }_i-{\gamma }_{i+1}}},\text{ for }{\rho }_{\text{c},k-1}\le {\rho }_{\text{g}}<{\rho }_{\text{c},k}\text{, }$$(3)

where $c_{\text{s},0}=\sqrt{ℛT_{\text{in}}/\mu }$ is the initial isothermal sound speed, T_in the initial gas temperature, 𝓡 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 piece-wise 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_{\text{s}}^2{\rho }_{\text{g}}^{{\gamma }_1}$. 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.

2.2 Dust component

For the dust component we use the model originally developed in Vorobyov et al. (2018) with modifications to account for the back-reaction 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⁻³ µ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⁻³–10⁻⁴ limits in the bulk of the disk and approach 10⁻² 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⁻² (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 prestel-lar 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: $$N(a)=C\cdot a^{-\text{p}}$$(4)

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 $$\frac{\partial {\rho }_{\text{d},\text{sm}}}{\partial t}+\nabla \cdot \left({\rho }_{\text{d},\text{sm}}\mathit{v}\right)=-S\left(a_{\max }\right),$$(5) $$\frac{\partial {\rho }_{\text{d},\text{gr}}}{\partial t}+\nabla \cdot \left({\rho }_{\text{d},\text{gr}}u\right)=S\left(a_{\max }\right),$$(6) $$\frac{\partial }{\partial t}\left({\rho }_{\text{d},\text{gr}}\mathit{u}\right)+\nabla \cdot \left({\rho }_{\text{d},\text{gr}}\mathit{u}\otimes \mathit{u}\right)=-{\rho }_{\text{d},\text{gr}}\nabla \text{Φ}+{\rho }_{\text{d},\text{gr}}f+S\left(a_{\max }\right)\mathit{v},$$(7)

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)υ 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 $$f=\frac{1}{2m_{\text{d}}}{C}_{\text{D}}\sigma {\rho }_{\text{g}}(\mathit{v}-\mathit{u})|\mathit{v}-\mathit{u}|,$$(8)

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_{\text{d}}=8v_{\text{th}}/3|\mathit{v}-\mathit{u}|$ (Weidenschilling 1977), where υ_th is the mean thermal velocity of gas. The drag force then takes a simple form $$f=\frac{\mathit{v}-\mathit{u}}{t_{\text{stop }}},$$(9)

where t_stop is the stopping time expressed as $$t_{\text{stop }}=\frac{a{\rho }_{\text{s}}}{{\rho }_{\text{g}}{v}_{\text{th }}},$$(10)

where ρ_s = 3.0 g cm⁻³ 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 $\sqrt{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.

2.3 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 $$\frac{\partial a_{\max }}{\partial t}+(\mathit{u}\cdot \nabla )a_{\max }=\mathcal{D},$$(11)

where the rate of dust growth due to collisions and coagulation is computed using a monodisperse model (Birnstiel et al. 2012) $$\mathcal{D}=\frac{{\rho }_{\text{d}}{u}_{\text{rel}}}{{\rho }_{\text{s}}}.$$(12)

This rate includes the total volume density of dust ρ_d = ρ_d,sm + ρ_d,gr, dust material density ρ_s, and relative velocity of grain-to-grain collisions defined as $u_{\text{rel}}={\left(u_{\text{br}}^2+u_{\text{turb}}^2\right)}^{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) $$u_{\text{turb }}=\sqrt{\frac{3\alpha }{\text{St}+{\text{St}}^{-1}}{c}_{\text{s}}}\text{, }$$(13)

where c_s is the adiabatic speed of sound, a is the parameter describing the strength of turbulence in the disk, and St is the Stokes number defined as St = t_stopΩ_ĸ. Here, Ω_K is the Kep-lerian 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 Ω_ĸ but that worked poorly in disks perturbed by spiral density waves and mass infall from the envelope. In this paper, we set Ω_ĸ 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 Ω_ĸ. This, however, requires determining the disk-to-envelope interface on the fly and is subject to future code improvements.

The value of a 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⁻⁴ (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 gravitoturbu-lence (Riols & Latter 2018). If gravitoturbulence modifies the value of a, 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 $$S\left(a_{\max }\right)=-\frac{\text{Δ}{\rho }_{\text{d},\text{sm}}}{\text{Δ}t},$$(14)

where $\text{Δ}{\rho }_{\text{d},\text{sm}}={\rho }_{\text{d},\text{sm}}^{n+1}-{\rho }_{\text{d},\text{sm}}^n$ 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 $$\Delta {\rho }_{\text{d},\text{sm}}={\rho }_{\text{d},\text{sm}}^{n+1}-{\rho }_{\text{d},\text{sm}}^n=\frac{{\rho }_{\text{d},\text{gr}}^n{\displaystyle {\int }_{a_{\min }}^{a_{*}}{a}^{3-\text{p}}}\text{d}a-{\rho }_{\text{d},\text{sm}}^n{\displaystyle {\int }_{a_{*}}^{a_{\max }^{n+1}}{a}^{3-\text{p}}}\text{d}a}{{\displaystyle {\int }_{a_{\min }}^{a_{\max }^{n+1}}{a}^{3-\text{p}}}\text{d}a},$$(15)

where the indices n and n +1 denote the current and next hydro-dynamic 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 $$a_{\text{frag }}=\frac{{\rho }_{\text{g}}{u}_{\text{frag }}^2}{3\alpha {\rho }_{\text{s}}{v}_{\text{th}}{\text{Ω}}_{\text{K}}}\text{, }$$(16)

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 Ω_ĸ 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.

2.4 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 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.

2.5 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–Lax– van 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 poly-tropic 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 $$\frac{\partial \mathit{Q}}{\partial t}+\nabla \cdot \mathit{F}(\mathit{Q})=0,$$(17)

where Q is the vector of conservative variables (ρ_g, ρ_d,sm,ρ_gυ, E) and F(Q) is the flux of conservative variables (ρ_gυ,ρ_d,smυ,ρ_gυ ⊗ υ + 𝕀?P, (E + P)υ). 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 piece-wise 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 $$\frac{\partial {\mathit{Q}}_d}{\partial t}+\nabla \cdot \mathit{F}\left({\mathit{Q}}_d\right)=0,$$(18)

where Q_d is the vector of conservative variables (ρ_d,gr,ρ_d,gra_max,ρ_d,gru) and F(Q_d) is the flux of conservative variables (ρ_d,gru,ρ_d,gra_maxu, ρ_d,gru ⊗ 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 $$\text{Φ}\left(x_i,y_j,z_k\right)=-G{\displaystyle \sum _{i^{\prime },j^{\prime },k^{\prime }}\frac{M\left(x_{i^{\prime }},y_{j^{\prime }},z_{k^{\prime }}\right)}{\sqrt{{\left(x_i-x_{i^{\prime }}\right)}^2+{\left(y_j-y_{j^{\prime }}\right)}^2+{\left(z_k-z_{k^{\prime }}\right)}^2}}},$$(19)

where $M\left(x_{i^{\prime }},y_{j^{\prime }},z_{k^{\prime }}\right)$ is the mass of gas and dust in a numerical cell (i′, j′, k′) with coordinates $\left(x_{i^{\prime }},y_{j^{\prime }},z_{k^{\prime }}\right)$, 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) $${\mathit{v}}^{n+1}=\epsilon D+E,{\mathit{u}}^{n+1}=-D+E,$$(20) $$E=\frac{C_1}{\epsilon +1},D=\frac{C_2}{\epsilon +1}\exp \left(-\frac{\epsilon +1}{t_{\text{stop }}}\text{Δ}t\right),$$(21) $$C_1={\mathit{v}}^n+\epsilon {\mathit{u}}^n,C_2={\mathit{v}}^n-{\mathit{u}}^n,$$(22)

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 $$\frac{\partial \left({\rho }_{\text{d},\text{gr}}{a}_{\text{max}}\right)}{\partial t}+\nabla \cdot \left({\rho }_{\text{d},\text{gr}}{a}_{\max }u\right)={\rho }_{\text{d},\text{gr}}\mathcal{D}+a_{\max }S\left(a_{\max }\right).$$(23)

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 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 𝒟 is set to zero during the collapse and disk formation stage.

2.6 Hybrid parallelization approach

ngFEOSAD is built using a combination of Coarray Fortran and OpenMP. Coarray Fortran extends modern Fortran, specifically Intel Fortran in our case, to distributed memory systems using the Single Program Multiple Data (SPMD) model, in which all multiprocessors execute the same program but on different data. This extension is realised through the addition of a datatype to arrays called a codimension and with a small number of additional directives. It is built on the Message Passing Interface (MPI) standard but only uses one-sided communication as opposed to pure MPI which is two-sided. This allows a multiprocessor to directly access the memory of another multiprocessor without requiring explicit cooperation from the target multiprocessor, something called Direct Memory Access (DMA). This capability simplifies the parallel programming process and can lead to more readable code because it reduces the need for complex coordination between multiprocessors. OpenMP, on the other hand, is used for the shared memory parallelism possible with multiprocessors and works within each distributed Coarray image, requiring only a further small number of additional directives. When combined, Coarray Fortran and OpenMP enable a highly optimized, parallel implementation of our key technique – the Cartesian nested grid.

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³ 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 step-by-step solution procedure and more details on the parallel realization using the Coarray Fortran technique can be found in Appendix A.

2.7 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³ 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.

3 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.

3.1 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.57 law. For a gravitationally unstable disk that settles in a state with Q ~ 1, the n_g ≃ r⁻³ 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 $$H_{\text{g}}=\frac{{\text{Σ}}_{\text{g}}}{\sqrt{2\pi }{\rho }_{\text{g}}(z=0)},$$(24)

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 x-coordinate 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. Postpocess-ing 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).

Fig. 1 Number density for different time instances in the inner 100×100 au2 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.

3.2 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 of its rarefied density and low turbulent velocity, see Eq. (12). We acknowledge that our dust growth model is designed for proto-planetary 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 amax 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.

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 proto-planetary 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 $${\text{Σ}}_{\text{MMSN }}=1700{\left(\frac{r}{1.0\text{\hspace{0.17em}au}}\right)}^{-1.5}\left[ \text{g} {\text{cm}}^{-2}\right].$$(25)

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 pro-tostellar 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 $${\xi }_{\text{d}2 \text{g}}(z=0)=\frac{{\rho }_{\text{d},\text{gr}}(z=0)+{\rho }_{\text{d},\text{sm}}(z=0)}{{\rho }_{\text{g}}(z=0)},$$(26)

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 nonhomo-geneous, 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 $${\xi }_{\text{d}2 \text{g}}={10}^{-2.7}{n}_{\text{g}}^{+0.08}\left[ {\text{cm}}^{-3}\right].$$(27)

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. (2018, 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.

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 Disk vertical scale height Hg (blue circles) and the gas number density in the midplane ngas (red crosses) as a function of distance in both directions from the disk center along the x-coordinate.

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

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 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 Total dust-to-gas mass ratio ξd2g as a function of the gas number density ng. 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 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.

3.3 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_{\text{peb}}^0=0.5$ mm and St ≥ St₀ = 0.01. More specifically, we calculate the minimum pebble size as $$a_{\text{peb}}^{\min }=\max \left(a_{\text{peb}}^0,a_{\text{St}}^0\right),$$(28)

where $a_{\text{St}}^0$ is the size of dust grains that have the Stokes number equal to St₀ for the local disk conditions, $a_{\text{St}}^0=a_{\text{max}}{\text{St}}_0/\text{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_{\text{peb}}^{\min }$ (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 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⁻³ (as in the current work) to 10⁻⁴ 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₀ 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 $$\frac{H_{\text{d}}}{H_{\text{g}}}\simeq \frac{{\text{Σ}}_{\text{d},\text{gr}}}{{\text{Σ}}_{\text{g}}}\times {\left({\xi }_{\text{d}2 \text{g}}(z=0)-\frac{{\text{Σ}}_{\text{d},\text{sm}}}{{\text{Σ}}_{\text{g}}}\right)}^{-1}=\frac{{\text{Σ}}_{\text{d},\text{gr}}}{{\text{Σ}}_{\text{g}}}\frac{{\rho }_{\text{g}}}{{\rho }_{\text{d},\text{gr}}}.$$(29)

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, υ_ϕ > 2υ_r,

• the gas rotational velocity must be faster than twice the vertical infall velocity, υ_ϕ > 2υ_z,

• gas must not be thermally supported, $\frac{\rho v_{\varphi }^2}{2}>2P$,

• the gas number density must be higher than 10⁹ cm⁻³.

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.

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.

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

3.4 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 disk. Dust settling in the low-rotation model is slightly more efficient, likely because the spiral density waves are weaker in this model. The difference, however, is marginal.

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

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

4 Model caveats and future developments

This is our first step toward fully three-dimensional numerical hydrodynamics modeling of the formation of protoplanetary gas-dust 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. 2018, 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⁻³ (adopted here) may have a strong impact on the dust and pebble distribution in the disk. Moreover, protoplane-tary 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³ when the code is updated to use the GPUs for the gravitational potential solver. The current resolution with N = 64³ 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.

5 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⁻³ –10⁻¹ 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 as compared to the corresponding characteristics in the MMSN. Low Stokes, along with a nonsteady character of the spiral pattern, hinderstrong dustaccumulationinthespiral 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 non-homogeneous and patchy character. Interestingly, pebbles are absent in the inner 10 au because of the low St < 10⁻² 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.

Acknowledgements

We thank the anonymous referee for constructive and useful comments and suggestions. This work was supported by the FWF project I4311-N27 (E.I.V., J.M.) and RFBR project 19-51-14002 (I.K. and V.E.). Simulations were performed on the Vienna Scientific Cluster (VSC) (https://vsc.ac.at/) and on the Narval Cluster provided by Calcul Quebec (https://www.calculquebec.ca/) and the Digital Research Alliance of Canada (https://alliancecan.ca/).

Appendix A Solution procedure and parallel realization in Coarray Fortran

In this section, we provide essential details of the solution procedure and parallelization technique employed in ngFEOSAD. For a parallel implementation of the solution procedure on nested grids, we use a geometric decomposition in which each grid is placed on a separate image of the Coarray Fortran (CAF) program (see Figure A.1). Below we present the main procedures of the ngFEOSAD computational cycle.

• 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 (Car-ray 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 self-gravity 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.

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.

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, υ_L, υ_R, ω_L, ω_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: $$Q_{\text{L}}=\left(\begin{array}{c}{\rho }_{\text{L}}\\ {\rho }_{\text{L}}{u}_{\text{L}}\\ {\rho }_{\text{L}}{v}_{\text{L}}\\ {\rho }_{\text{L}}{w}_{\text{L}}\end{array}\right),Q_{\text{R}}=\left(\begin{array}{c}{\rho }_{\text{R}}\\ {\rho }_{\text{R}}{u}_{\text{R}}\\ {\rho }_{\text{R}}{v}_{\text{R}}\\ {\rho }_{\text{R}}{w}_{\text{R}}\end{array}\right)\text{, }$$(B.1) $$F_{\text{L}}=\left(\begin{array}{c}{\rho }_{\text{L}}{u}_{\text{L}}\\ {\rho }_{\text{L}}{u}_{\text{L}}{u}_{\text{L}}\\ {\rho }_{\text{L}}{v}_{\text{L}}{u}_{\text{L}}\\ {\rho }_{\text{L}}{w}_{\text{L}}{u}_{\text{L}}\end{array}\right),F_{\text{R}}=\left(\begin{array}{c}{\rho }_{\text{R}}{u}_{\text{R}}\\ {\rho }_{\text{R}}{u}_{\text{R}}{u}_{\text{R}}\\ {\rho }_{\text{R}}{v}_{\text{R}}{u}_{\text{R}}\\ {\rho }_{\text{R}}{w}_{\text{R}}{u}_{\text{R}}\end{array}\right)$$(B.2)

• 2.

  We calculate the exact value of the dust transport rate according to the following equation: $$S=\frac{\sqrt{{\rho }_{\text{L}}}{u}_{\text{L}}+\sqrt{{\rho }_{\text{R}}}{u}_{\text{R}}}{\sqrt{{\rho }_{\text{L}}}+\sqrt{{\rho }_{\text{R}}}},$$(B.3)

• 3.

  We write down a one-wave HLL scheme in the following basic form: $$F^{1\text{ WHLL }}=\{\begin{array}{ll}{F}_{\text{L}},\hfill & 0\le S,\hfill \\ F_{\text{R}},\hfill & 0\ge S\hfill \end{array}.$$(B.4)

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.

Appendix C Numerical tests

In this section, we provide the results of extensive testing of the main constituent parts of the ngFEOSAD code.

Appendix C.1 Tests on the numerical scheme for gas dynamics

Verification of the numerical method for solving hydrodynamic equations for the gas component was carried out on a set of classical one-dimensional Godunov tests, tests that consider the development of physical instabilities of the Kelvin-Helmholtz and Rayleigh-Taylor types, as well as on the Sedov test of a point explosion. Such a test suite is aimed at checking the method for correctly reproducing the characteristic gas flows and has either analytical solutions or a general estimate of the nature of the gas flow during the growth of instabilities. All presented tests in this section were performed on a single Cartesian grid with equidistant cell spacing.

Appendix C.1.1 The Godunov tests for gas dynamics

Here, we consider a set of Godunov tests in the form of the Sod problem with a sonic point, the Einfeldt problem with the formation of a strong rarefaction region, and the Lax problem with the formation of strong shock waves. The initial configuration for these three test problems is given in Table C.1, where x₀ is the initial position of the boundary between the two adjacent states (the subscripts L and R denote the left-hand side and right-hand side states, respectively). For computational experiments, 200 computational cells were used.

The main goal of the Sod problem is to test the ability of the numerical method to simultaneously reproduce a weak shock wave, a contact discontinuity, and a rarefaction wave. Unlike the standard realizations of this test problem, our initial setup is intentionally complicated by the presence of a nonzero velocity at the left state, which forms initially a supersonic gas flow. Roe-type methods, when reproducing such a solution, may have a feature in the form of an entropy trace at the point of the discontinuity of the rarefaction wave. In our version of the method (see Figure C.1), all waves are reproduced correctly, including the absence of the entropy trace. The use of piece-wise parabolic reconstruction of physical variables makes it possible to reproduce a weak shock wave with a numerical dissipation of only two cells.

The second Godunov test, the Einfeldt test, assesses the ability of a numerical method to reproduce a rarefaction region caused by gas expansion in the opposite directions. The main problem of the test is the nonphysical growth of the internal energy at the boundary between the two initial states, with a fairly good quality of reproduction of the pressure, density, and velocity distributions. We also found it to be the case in our numerical method – the internal energy grows nonphysically in the rarefaction region (see Figure C.2), but such growth is limited in amplitude, rather smooth, and does not have oscillations, which speaks in favor of the chosen numerical technique.

Fig. C.1 Sod test problem: panel (a) displays the gas density, panel (b) – pressure, panel(c) – velocity, and panel (d) – internal energy. The symbols present the numerical data, while the lines correspond to the analytical solution.

Fig. C.2 Einfeldt test problem: density (a), pressure (b), velocity (c), and internal energy (d). The symbols present the numerical data, while the lines correspond to the analytical solution.

In the Lax test, the ability of the code to accommodate a significant pressure drop (5 orders of magnitude) is assessed. The key feature of the flow is the low-amplitude precursor wave that can be seen in the internal energy ahead of the shock wave. Our computational experiments (see Figure C.3) indicate that such a wave is reproduced correctly without significant dissipation. We also note the absence of unphysical oscillations on the shock wave and the contact discontinuity, which may occur in high order schemes.

Fig. C.3 Lax test problem: density (a), pressure (b), velocity (c), and internal energy (d). The symbols present the numerical data, while the lines correspond to the analytical solution.

Appendix C.1.2 The Kelvin-Helmholtz and Rayle¡gh-Taylor instabilities

An important property of numerical methods is their ability not to suppress physical instabilities when they are numerically resolved. To verify this property of ngFEOSAD, we study the development of Kelvin-Helmholtz and Rayleigh-Taylor instabilities. The Rayleigh-Taylor instability develops at the boundary of two gases with different densities in the common gravitational field if the heavier component is located on top of the lighter one. The Kelvin-Helmholtz instability develops at the boundary of two gases with different densities moving relative to each other. Both types of instabilities lead to the development of nonlinear hydrodynamic turbulence.

Kelvin-Helmholtz instability. We consider a square domain [−0.5; 0.5]². The initial density profile and the x-velocity component are given as: $${\rho }_0(x,y)=\{\begin{array}{ll}1,\hfill & \left|y\right|>0.25+0.01(1+\cos (8\pi x))\hfill \\ 2,\hfill & \left|y\right|\le 0.25+0.01(1+\cos (8\pi x))\hfill \end{array}$$(C.1) $$v_{x,0}(x,y)=\{\begin{array}{ll}-0.5,\hfill & \left|y\right|>0.25+0.01(1+\cos (8\pi x))\hfill \\ 0.5,\hfill & \left|y\right|\le 0.25+0.01(1+\cos (8\pi x))\hfill \end{array}$$(C.2)

The gas pressure is set equal to p = 2.5 and the adiabatic index is γ = 1.4. For the development of instability, a sinusoidal perturbation of the gas density and velocity is introduced. The development of the Kelvin-Helmholtz instability is shown in Figure C.4. The characteristic time for the development of the Kelvin-Helmholtz instability can be estimated as: $${\tau }_{\text{KH}}=\frac{\lambda \left({\rho }_{\text{in }}+{\rho }_{\text{out }}\right)}{v_{\text{rel }}\sqrt{{\rho }_{\text{in }}{\rho }_{\text{out }}}},$$(C.3)

where λ = 1/4 is the inverse frequency of the sinusoidal disturbance, υ_rel = 1 the relative velocity of motion between gases of different densities, ρ_in = 2 the density of the inner denser region, ρ_out = 1 the density of the outer rarefied region. For such parameters, the characteristic development time of the Kelvin-Helmholtz instability is τ_KH ≈ 0.53, which agrees well with the results of computational experiments shown in Figure C.4. Indeed, the formation of turbulent eddies occurs at the boundary of two gases at t ≈ 0.5.

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.

Rayleigh-Taylor instability. We consider the area [−0.25; 0.25] × [0; 1.5] with the following parameters: $${\rho }_0(x,y)=\{\begin{array}{ll}1,\hfill & y>0.75\hfill \\ 2,\hfill & y\le 0.75,\hfill \end{array}$$(C.4) $$p=0.15-\rho \cdot g\cdot \left|y\right|\text{, }$$(C.5)

where p is the hydrostatic equilibrium pressure, g = 0.1 the gravitational acceleration, y the vertical coordinate, and γ = 1.4 the adiabatic exponent. The vertical component of velocity is perturbed according to the following rule: υ_y,0(x, y) = A(y − O.75)[l + cos(2πx)][l + cos(2πy)], where $$A(y)=\{\begin{array}{ll}{10}^{-2},\hfill & y\le 0.01\hfill \\ 0,\hfill & y>0.01\hfill \end{array}.$$(C.6)

The development of the Rayleigh-Taylor instability is shown in Figure C.5. The growth of the perturbation amplitude of the Rayleigh-Taylor instability is determined as: $$f(t)=0.01\mathrm{\hspace{0.17em}exp\hspace{0.17em}}(t\sqrt{Ag})\simeq 0.01\mathrm{\hspace{0.17em}exp\hspace{0.17em}}(0.25t),$$(C.7)

where A = 2/3 is Atwooď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.

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.

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]³ the adiabatic exponent γ = 5/3, the initial density in the domain ρ₀ = 1, and the initial pressure ρ₀ = 10⁻⁵. At time t = 0, the internal energy e₀ = 0.6 is released. The explosion area is limited by the radius r₀ = 0.01. The computational grid with 100³ cells is used. The density and momentum profiles at t = 0.05 are shown in Figure C.6.

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.

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; Stoy-anovskaya 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_st0p = ρ_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 code reproduces well the position of the shock and rarefaction waves. The shocks are spread over just 2–3 grid zones, which can be expected for the 3rd order accurate parabolic interpolation scheme used in ngFEOSAD.

Appendix C.2.2 Dusty wave

The dusty wave problem tests the propagation and damping of linear waves in the mixture of gas and dust. The advantage of this test is that it has analytical solutions for both weak and strong coupling between gas and dust. The initial setup consists of a sinusoidal wave of the following form $${\rho }_0=\widetilde{{\rho }_{\text{g}}}+\delta \sin (kx),{\rho }_{\text{d},0}=\widetilde{{\rho }_{\text{d}}}+\delta \sin (kx),$$(C.8) $$v_0=\delta \sin (kx),u_0=\delta \sin (kx).$$(C.9)

Fig. C.7 Sod shock tube problem for a strongly coupled gas-dust mixture. The top and middle panels present solutions for the gas component, while the bottom rows show solutions for the dust component. The red lines present the analytical solution and the filled circles are the numerical solution. The green line the upper left panel presents the analytic solution for the pure gas for comparison.

Fig. C.8 Solution of the dusty wave problem for weak (top row) and strong (bottom row) coupling between gas and dust. The left and right columns present the gas and dust velocities, respectively. The black dashed curves are the numerical solutions, while the red line shows the analytical solution.

We adopt isothermal gas with $c_s=1,\widetilde{{\rho }_{\text{g}}}=1,k=1$, and δ = 0.01 throughout this section and vary K and $\widetilde{{\rho }_{\text{d}}}$. All simulations were done with 200 grid zones on unit distance and periodic boundary conditions. The Courant number is set to 0.2. The analytical solution is taken from Laibe & Price (2011).

Figure C.8 presents the gas and dust velocities at t = 0.5 for K = 1 and K = 1000, namely, for the cases of weak and strong coupling. The dust-to-gas ratio is $\widetilde{{\rho }_{\text{d}}}/\widetilde{{\rho }_{\text{g}}}=1.0$. We achieve a reasonably good agreement with the analytical solution for both the gas and dust velocities, and also in the limits of weak and strong coupling between gas and dust. A small mismatch between the analytic and numerical solution is typical for this test (Laibe & Price 2012; Lorén-Aguilar & Bate 2014). Finally, we note that both tests, the Sod shock tube and the dusty wave, took the back reaction of dust on gas into account.

Appendix C.2.3 Sod shock tube for pure dust

Finally, we test the numerical solver on pure dust without gas. We consider the initial setup similar to that shown in Table C.1 for the Sod test for gas, noting only that dust is pressureless. The key in this test is the correct reproduction of the dust front propagation velocity, low dissipation of the discontinuity, as well as the reproduction of the delta function at the discontinuity front as a function of the dust density. Figure C.9 demonstrates that all features of the solution are reproduced correctly and with low dissipation.

Fig. C.9 Sod problem test for pure dust: density (a) and velocity (b). The lines present the analytic solution while the symbols correspond to the numerical data.

Appendix C.3 Gravitational potential solver

Unlike many other grid-based codes that solve the Poisson equation to find the gravitational potential, ngFEOSAD employs the convolution theorem to compute the triple sum (Eq. 19). A detailed explanation of the method, its extension to nested grids, and comparison with a more common method of conjugate gradients, as well as the pertinent tests can be found in Vorobyov et al. (2023b). The solver is characterized by a global second-order convergence, except for the nested grid interfaces where the convergence is linear. The mean errors on the standard tests – an oblate ellipsoid and a wide-separation binary – stay below 0.05% and 0.5%, respectively.

The advantage of our original method is that it can be easily parallelized on nested grids and does not require finding the boundary potential on the coarsest grid via the multipole expansion. The method can be significantly accelerated using the graphic processing units, although in the current work we do not utilize this feature.

Appendix C.4 Global collapse problems

In this section, we consider the test problems that are most relevant to modeling the gravitational collapse of prestellar cores. These problems are inclusive in the sense that they test the code performance on nested grids for both the gravitational potential and hydrodynamics solvers.

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.

Appendix C.4.1 Collapse of a cold nonrotating gaseous sphere

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⁴ 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² should be conserved and equal to (Norman et al. 1980) $$M(J)=M_0\left[1-\left(1-\frac{J}{{\text{Ω}}_0{R}_0^2}\right)\right].$$(C.10)

Here, Ω is the angular velocity of the collapsing cloud at a distance r from the z-axis, and M₀, Ω₀, and R₀ 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 gravitational 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 R0 = 10⁴ au, n₀ = 5 × 10⁴ cm^–3, and Ω₀ = 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³.

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.

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 nonconservation is only a minor fraction of the entire cloud mass, namely, less than 0.04%.

Appendix C.4.3 Protostellar accretion in spherical collapse problems

The accretion rate of a spherically collapsing cloud on to a pro-tostar has an analytic solution of the following form (Shu 1977; Kratter et al. 2010) $$\dot{M}=\frac{c_{\text{s}}^3}{G}\times \{\begin{array}{ll}0.975,\hfill & (A=2)\hfill \\ {(2A)}^{3/2}/\pi ,\hfill & (A\gg 2)\hfill \end{array}$$(C.11)

for the initial gas density distribution defined as $$\rho (r)=\frac{A{c}_{\text{s}}^2}{4\pi G{r}^2}.$$(C.12)

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⁵ 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 th computational cell within the sink sphere exceeds ρ_crit = m_H µ n_crit, where n_crit = 10¹³ 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 $$\frac{d{\rho }_i}{dt}=-b\left({\rho }_i-{\rho }_{\text{crit }}\right),$$(C.13)

which has an analytical solution of the following form $${\rho }_i={\rho }_{0,i}{\exp }^{-b\text{Δ}t}+{\rho }_{\text{crit }}\left(1.0-{\exp }^{-b\text{Δ}t}\right).$$(C.14)

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 $$\dot{M}={\displaystyle \sum _i\frac{{\rho }_{0,i}-{\rho }_i}{\text{Δ}t}}\times \text{Δ}{V}_i,$$(C.15)

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 $\dot{M}$is also found in other collapse simulations (e.g., Vorobyov & Basu 2005). Soon the mass accretion rate drops and stabilizes at ≈ 2.5 × 10^–5M_⊙ 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. 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.

References

All Tables

All Figures

	Fig. 1 Number density for different time instances in the inner 100×100 au2 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.
In the text

	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.
In the text

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

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

In the text

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.

In the text

	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.
In the text

	Fig. 7 Total dust-to-gas mass ratio ξd2g as a function of the gas number density ng. 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.
In the text

	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.
In the text

	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.
In the text

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

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

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

	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.
In the text

	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.
In the text

	Fig. C.1 Sod test problem: panel (a) displays the gas density, panel (b) – pressure, panel(c) – velocity, and panel (d) – internal energy. The symbols present the numerical data, while the lines correspond to the analytical solution.
In the text

	Fig. C.2 Einfeldt test problem: density (a), pressure (b), velocity (c), and internal energy (d). The symbols present the numerical data, while the lines correspond to the analytical solution.
In the text

	Fig. C.3 Lax test problem: density (a), pressure (b), velocity (c), and internal energy (d). The symbols present the numerical data, while the lines correspond to the analytical solution.
In the text

	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.
In the text

	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.
In the text

	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.
In the text

	Fig. C.7 Sod shock tube problem for a strongly coupled gas-dust mixture. The top and middle panels present solutions for the gas component, while the bottom rows show solutions for the dust component. The red lines present the analytical solution and the filled circles are the numerical solution. The green line the upper left panel presents the analytic solution for the pure gas for comparison.
In the text

	Fig. C.8 Solution of the dusty wave problem for weak (top row) and strong (bottom row) coupling between gas and dust. The left and right columns present the gas and dust velocities, respectively. The black dashed curves are the numerical solutions, while the red line shows the analytical solution.
In the text

	Fig. C.9 Sod problem test for pure dust: density (a) and velocity (b). The lines present the analytic solution while the symbols correspond to the numerical data.
In the text

	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.
In the text

	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.
In the text

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.

In the text