© The Authors 2023

1 Introduction

With a typical lifetime of a few million years (Haisch et al. 2001; van der Marel & Mulders 2021), protoplanetary disks are known to rapidly accrete their gas and dust content onto the central pre-main-sequence star, with a typical accretion rate of ~10⁻⁹−10⁻⁸ M_⊙ yr⁻¹ (e.g., Hartmann et al. 2016). The accretion phenomenon is ultimately controlled by angular momentum transport and outflow mass-loss processes. Such processes shape the disk structure and global evolution as well as its dispersal, hence being of crucial importance in understanding the first steps of planet formation (Armitage 2011, 2019). Particularly, the disk turbulence shapes the global density distribution within which planets forms, and it strongly impacts the evolution of dust particles which are the building blocks of planets (Dubrulle et al. 1995; Ormel & Cuzzi 2007; Youdin & Lithwick 2007; Birnstiel et al. 2010). Despite its fundamental role, the nature of turbulence in protoplanetary disks has not been well constrained yet. Pure hydrodynamic-driven mechanisms have been suggested for generating turbulence, such as vertical shear instability (VSI; e.g., Urpin & Brandenburg 1998; Nelson et al. 2013; Lin & Youdin 2015; Manger et al. 2020; Flock et al. 2020) or baroclinic instabilities (e.g., Klahr & Bodenheimer 2003; Raettig et al. 2013). Nevertheless, theoretical studies show that the resulting turbulence level is typically too weak to explain the observed accretion rates of protoplanetary disks (see the review of Lesur et al. 2022). Other mechanisms such as gravitational instability (GI; e.g., Lin & Pringle 1987; Lodato & Rice 2004; Vorobyov & Basu 2009) can provide significant turbulence, but only in the early stages of disk evolution. Consequently, magnetohydro-dynamic (MHD)-driven mechanisms such as MHD winds (e.g., Blandford & Payne 1982; Suzuki & Inutsuka 2009; Bai et al. 2016; Bai 2016) and magnetorotational instability (MRI; e.g., Balbus & Hawley 1991, 1998; Hawley et al. 1995) are currently thought to be the main candidates for driving disk accretion.

The MRI-driven accretion can be substantially modified and even suppressed at some locations in the disk by three nonideal MHD effects: Ohmic resistivity, the Hall effect, and ambipo-lar diffusion. They arise due to the weak level of ionization in the disk (e.g., Gammie 1996; Fleming et al. 2000; Sano & Stone 2002a,b; Fleming & Stone 2003; Inutsuka & Sano 2005; Ilgner & Nelson 2006; Turner et al. 2007, 2010; Turner & Sano 2008; Perez-Becker & Chiang 2011) or because of the high drift between the ions and electrons in regions of low gas densities and strong magnetic field strengths (e.g., Bai & Stone 2011). In general, how much and where the MRI activity is suppressed is a complex problem that significantly depends on the dust and gas properties, the magnetic field strength, as well as the complex ionization chemistry (e.g., Bai & Goodman 2009; Bai 2011a; Delage et al. 2022). The direct consequence is that a magnetically dead zone arises, and it is characterized by a low gas accretion rate (e.g., Dzyurkevich et al. 2010). It causes a steep increase in the disk turbulence at the transition from the dead zone to the MRI active region, the so-called dead zone outer edge. This location has been previously hypothesized to be a prime location for dust particle trapping, hence potentially explaining some of the observed substructures in pro-toplanetary disks (e.g., Regály et al. 2012; Flock et al. 2015; Pinilla et al. 2016).

To further investigate the potential dust trapping power of the dead zone outer edge, one needs to build a model whose outputs can be compared to current dust continuum and gas observations of million-year-old disks. Such a model thus necessarily requires a time-dependent framework where nonideal MHD calculations are self-consistently combined with gas and dust evolution on million-year timescales. Some theoretical works have studied the detailed behavior of MRI active and nonactive regions by performing 3D local shearing box or global simulations (e.g., Dzyurkevich et al. 2010; Turner et al. 2010; Bai & Stone 2011; Bai 2011a; Okuzumi & Hirose 2011; Flock et al. 2011, 2015). However, these studies did not implement a full treatment for dust evolution (dynamics and grain growth processes combined), as well as could not evolve the protoplanetary disk over millions of years. Conversely, some papers used the Shakura-Sunyaev α-disk model (Shakura & Sunyaev 1973), wherein the quantity α encodes the disk turbulence level in order to make possible, on million-year timescales, the implementation of gas and dust evolution altogether with an educated guess for the MRI-induced α parameter (e.g., Pinilla et al. 2016). Crucially, though, ad hoc prescriptions of α have been adopted, without accounting for the detailed physics of the MRI.

A more consistent approach to access the dead zone outer edge as a potential location for dust trapping thus has yet to be developed since the interplay between dust evolution and MRI-driven accretion has been poorly understood. Particularly, it is still unclear how dust evolution impacts the MRI-driven accretion in protoplanetary disks on million-year timescales. Indeed, most previous works have not implemented the feedback of dust evolution on the ionization state of the disk, which is crucial to accurately describe the MRI activity. A possible way to investigate such interplay is by using a “trade-off” model that combines a 1D viscous disk model and nonideal MHD calculations: the viscosity parameter α is determined by the MRI-driven turbulence accounting for the nonideal MHD effects, with a careful modeling of the gas ionization degree. Such a model allows one to capture the essence of the MRI in a noncomputationally expensive way, which makes the coupling with 1D gas and dust evolution models (including growth, fragmentation, settling, and radial drift) feasible on million-year timescales.

Okuzumi & Hirose (2012) investigated the impact of dust evolution on MRI-driven turbulence by coupling the dust coagulation equation, an analytic model for the disk ionization including charged grains, and an empirical model for α based on nonideal MHD simulations. Nevertheless, their work modeled the dust as a two-population phase instead of considering the full dust size distribution, neglected dust radial drift, and most importantly did not include ambipolar diffusion which is fundamental to accurately describe the MRI activity in the outer part of protoplanetary disks.

Delage et al. (2022) put forward a trade-off model specifically designed to study the MRI-driven accretion in the outer protoplanetary disk (r ≳ 1 au, where r is the distance from the central star). Their model allowed them to compute a self-consistent MRI-induced α parameter given stellar, gas, and dust properties in the framework of viscously driven accretion, accounting for both Ohmic resistivity and ambipolar diffusion with a careful modeling of the gas ionization degree. In their paper, they provide some insights regarding the potential impact of gas and stellar evolution on the MRI activity, particularly the dead zone outer edge. They show that the MRI-driven turbulence becomes stronger when the total disk gas mass decreases due to a higher gas ionization degree. Additionally, they indirectly show that a higher stellar X-ray luminosity leads to stronger MRI-driven turbulence overall due to a higher stellar X-ray ion-ization rate. However, their study assumes a fixed grain size across the whole protoplanetary disk, and they did not investigate the interplay between dust evolution and MRI-driven accretion over millions of years.

The aim of this present paper is thus twofold: (1) to understand how the implementation of a dust size distribution impacts the MRI activity, especially the steady-state MRI-driven accretion described in Delage et al. (2022); and (2) to present a pilot study providing a proof of concept that dust evolution alone has substantial impact on the MRI-driven turbulence, particularly the dead zone outer edge. Such a study solely focusing on the effect of dust evolution is a necessary first step toward a better understanding of the MRI–dust coevolution. Here we note that the accretion driven by magnetic disk winds is ignored to focus our efforts on the MRI.

The layout of the paper is as follows. In Sect. 2, we present the disk model employed. In Sect. 3, we describe the numerical implementation of the various simulations making use of our disk model. In Sect. 4, we present the results that investigate the impact of a dust size distribution as well as dust evolution on the MRI activity. In Sect. 5, we discuss the implications of our results. Finally, Sect. 6 summarizes our conclusions.

2 Disk model

We consider a central star of mass M_⋆ and bolometric luminosity L_⋆. Additionally, we assume that the envelope has dispersed to reveal a gravitationally stable and viscously accreting disk of total disk gas mass M_disk. The protoplanetary disk is considered geometrically thin, so that the vertical and radial dimensions can be decoupled into a 1+1D (r, z) framework, where each radial grid-point contains an independent vertical grid. Furthermore, we assume the disk to be axisymmetric and symmetric about the midplane.

The local mass and angular momentum transport are assumed to be solely controlled by the MRI and hydrodynamic instabilities based on the model of Delage et al. (2022), where the disk turbulence level is encoded into the viscosity parameter α. In their 1+1D model (r–z plane), the solution for vertical stratification of gas and dust is required in order to get an approximate 1D description of vertically layered accretion within the Shakura-Sunyaev viscous α-disk model. The main output of their model is thus the effective turbulent parameter $\overline{\alpha }$, defined as the pressure-weighted vertical average of the local turbulent parameter α: $$\overline{\alpha }\left(r\right)=\frac{{\displaystyle {\int }_{-\infty }^{+\infty }\alpha \left(r,z\right)P_{\text{gas}}\left(r,z\right)\text{\hspace{0.17em}d}z}}{{\displaystyle {\int }_{-\infty }^{+\infty }{P}_{\text{gas}}\left(r,z\right)\text{\hspace{0.17em}d}z}},$$(1)

where z is the height from the midplane, and P_gas is the isothermal gas pressure.

In Sects. 2.1 and 2.2, we present the disk properties for the gas and the dust, respectively. From these, we can obtain the effective turbulent parameter, $\overline{\alpha }$, by employing the MRI-driven turbulence model described in Sect. 2.3.

2.1 Gas

The gas is assumed vertically isothermal, with a radial temperature profile set by passively absorbing stellar irradiation $$T\left(r\right)={\left[T_{1\text{\hspace{0.17em}au}}^4{\left(\frac{r}{1\text{\hspace{0.17em}au}}\right)}^{-2}\left(\frac{L_{\star }}{L_{\odot }}\right)+T_{\text{bkg}}^4\right]}^{\frac{1}{4}},$$(2)

where T_{1 au} = 280 K is the gas temperature at r = 1 au for L_⋆ = L_⊙, and T_bkg = 10 K is the background gas temperature corresponding to the primordial temperature of the cloud prior to the collapse. We note that the choice made for the gas temperature model does not have a significant impact on the MRI-driven turbulence, since its temperature dependence is weak (see Appendix D.3 of Delage et al. 2022).

Assuming hydrostatic equilibrium in the vertical direction gives the gas volume density profile $${\rho }_{\text{gas}}\left(r,z\right)=\frac{{\displaystyle {\sum }_{\text{gas}}\left(r\right)}}{\sqrt{2\pi }{H}_{\text{gas}}\left(r\right)}\exp \left(-\frac{z^2}{2H_{\text{gas}}^2\left(r\right)}\right),$$(3)

where Σ_gas is the gas surface density chosen as explained in Sect. 3, and H_gas = c_s/Ω_K is the disk gas scale height with the isothermal sound speed $c_{\text{s}}=\sqrt{k_{\text{B}}T/\mu m_{\text{H}}}$ and the Keplerian angular velocity ${\text{Ω}}_{\text{K}}=\sqrt{G{M}_{\star }/r^3}$. Here k_B is the Boltzmann constant, µ = 2.34 is the mean molecular weight (assuming solar abundances), m_H is the atomic mass of hydrogen, and G is the gravitational constant. Finally, the total number density of gas particles is defined as n_gas = ρ_gas /m_neutral, with m_neutral = µm_H the mean molecular mass.

2.2 Dust particles

Observations of protoplanetary disks at different wavelengths suggest that the dust particles can significantly decouple from the gas, depending on the grain size (see e.g., Andrews et al. 2016; van Boekel et al. 2017; Huang et al. 2018). Consequently, including a dust model is required.

We considered that the dust phase consists of a distribution of dust particles with different sizes. Each grain is assumed to be a perfect compact sphere of intrinsic volume density ρ_bulk = 1.4 g cm⁻³ (consistent with the solar abundance when H₂O ice is included in grains; see Pollack et al. 1994). For a grain of size a, the corresponding mass is then $m\left(a\right)=\frac{4}{3}\pi {\rho }_{\text{bulk}}{a}^3$.

Most of the transport and dynamics of dust particles in pro-toplanetary disks are regulated by the interactions between them and the gas. A way to quantify the importance of the drag forces on the dynamics of a dust particle (hence the level of coupling between that dust particle and the gas) is by its Stokes number, defined as the dimensionless version of the stopping time of that particle. Near the midplane, the Stokes number of a dust particle of size a is given by $$\text{St}\left(r,a\right)=\frac{\pi }{2}\frac{a{\rho }_{\text{bulk}}}{{\displaystyle {\sum }_{\text{gas}}\left(r\right)}}\cdot \{\begin{array}{ll}1\hfill & {\lambda }_{\text{mfp}}/a\ge 4/9\hfill \\ \frac{4}{9}\frac{a}{{\lambda }_{\text{mfp}}}\hfill & {\lambda }_{\text{mfp}}/a<4/9\hfill \end{array},$$(4)

where ${\lambda }_{\text{mfp}}={\left(n_{\text{gas}}{\sigma }_{{\text{H}}_2}\right)}^{-1}$ is the mean free path of gas particles, with ${\sigma }_{{\text{H}}_2}=2\times {10}^{-15}{\text{cm}}^2$ the molecular cross-section for H₂ (e.g., Brauer et al. 2008; Birnstiel et al. 2010).

2.2.1 Dust settling

Dust models predict that grains tend to settle more efficiently toward the midplane as they grow in size (Dubrulle et al. 1995). Additionally, dust settling appears to be at play by ALMA observations of edge-on disks (see, e.g., Villenave et al. 2019, 2020). As a result, we can expect the number of dust particles to drop significantly above a dust scale height, H_dust, that can be much smaller than the gas scale height, H_gas. Assuming dust stirring to be induced by the MRI-driven turbulence, we can relate H_dust, H_gas, St, and $\overline{\alpha }$ by (Dubrulle et al. 1995; Youdin & Lithwick 2007; Yang et al. 2018) $$H_{\text{dust}}\left(r,a\right)=H_{\text{gas}}\left(r\right)\sqrt{\frac{\overline{\alpha }\left(r\right)}{\overline{\alpha }\left(r\right)+\text{St}\left(r,a\right)}.}$$(5)

This expression is given for each grain species of size a. We assume that St/D_gas is independent of z, where D_gas is the gas diffusion coefficient. We also implicitly approximate D_gas by the vertically integrated gas kinematic viscosity $$\overline{\nu }=\overline{\alpha }{c}_{\text{s}}{H}_{\text{gas}}.$$(6)

In the vertical direction, we then assume that the dust volume density profile follows a Gaussian distribution. For each grain species of size a, it is described by $${\rho }_{\text{dust}}\left(r,z,a\right)=\frac{{\displaystyle {\sum }_{\text{dust}}\left(r,a\right)}}{\sqrt{2\pi }{H}_{\text{dust}}\left(r,a\right)}\exp \left(-\frac{z^2}{2H_{\text{dust}}^2\left(r,a\right)}\right),$$(7)

where Σ_dust is the dust surface density for each grain species of size a that we describe in Sect. 2.2.2. The total dust volume density (accounting for all grain species) is thus defined as $${\rho }_{\text{dust,\hspace{0.17em}tot}}\left(r,z\right)={\displaystyle \sum _a{\rho }_{\text{dust}}\left(r,z,a\right).}$$(8)

Finally, the number density for each grain species of size a is given by $$n_{\text{dust}}\left(r,z,a\right)=\frac{{\rho }_{\text{dust}}\left(r,z,a\right)}{m\left(a\right)},$$(9)

with m(a) being the corresponding grain mass. It follows that the total number density (accounting for all grain species) is determined by $$n_{\text{dust,tot}}\left(r,z\right)={\displaystyle \sum _a{n}_{\text{dust}}\left(r,z,a\right).},$$(10)

2.2.2 Dust surface density

The remaining quantity to complete the description of the dust phase is the surface density for each grain species of size a (Σ_dust(r, a)). In this paper, it is determined either by assuming that the dust size distribution follows a fixed power-law, or is directly obtained from a dust evolution model.

Power-law dust size distribution

It is determined by the three parameters a_min, a_dist,Max, and p_dist,Exp; respectively the distribution minimum grain size, maximum grain size, and exponent. In the grain size range [a, a + da], it reads as follows: $${n^{\prime }}_{\text{dust}}\left(a\right)\text{d}a\propto \{\begin{array}{ll}{a}^{p_{\text{dist,Exp}}}\text{d}a\hfill & \text{if\hspace{0.17em}}{a}_{\min }\le a\le a_{\text{dist,Max}}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array},$$(11)

where ${n^{\prime }}_{\text{dust}}\left(a\right)$ refers to the dust number density per grain size, and is a distribution function over a. It differs from n_dust(a) above, which is already the quantity integrated over the bin size around a.

The dust surface density for each grain species of size a, Σ_dust(r, a), then follows from the conservation of the total dust mass: The quantity Σ_dust,tot(r) = Σ_a Σ_dust(r, a) must be equal to f_dg,tot(r)Σ_gas(r), where f_dg,tot(r) is the vertically integrated total dust-to-gas mass ratio (accounting for all grain species) at radius r. We note that f_dg,tot is a free parameter for the models of this present paper using the power-law dust size distribution (Models I–III). We assume this quantity to be radially constant equal to 10⁻².

Dust evolution model

The dust size distribution can substantially differ from a power-law because dust particles collide with each other leading to coagulation or fragmentation, and are transported across the disk by different mechanisms such as thermal Brownian motion, vertical stirring and settling, turbulent mixing, and drift (azimuthal and radial). As a result, we need a dust evolution model that includes all these mechanisms. For this purpose, we used the code DustPy¹ (Stammler & Birnstiel 2022), which can simulate the advection of gas and dust, along with the growth and fragmentation of multiple grain species, based on the model of Birnstiel et al. (2010). Below, we summarize the main ingredients in the model.

Following Nakagawa et al. (1986) and Takeuchi & Lin (2002), the dust radial velocity for each grain species of size a is given by: $${\upsilon }_{\text{r,dust}}\left(r,a\right)=\frac{{\upsilon }_{\text{gas}}\left(r\right)}{1+{\text{St}}^2\left(r,a\right)}-\frac{2\text{St}\left(r,a\right)}{1+{\text{St}}^2\left(r,a\right)}\eta \left(r\right){\upsilon }_{\text{K}}\left(r\right),$$(12)

with $\eta =-\frac{1}{2}{\left(H_{\text{gas}}/r\right)}^2$ (d lnP_gas,imd/d ln r) being the gas pressure support parameter (e.g., Birnstiel et al. 2016), v_K = rΩ_K the Keplerian orbital velocity, and P_gas,mid the midplane isothermal gas pressure $P_{\text{gas}}={\rho }_{\text{gas}}{c}_{\text{s}}^2$. Moreover, υ_gas corresponds to the gas viscous velocity, and is determined by $${\upsilon }_{\text{gas}}=-\frac{3}{{\displaystyle {\sum }_{\text{gas}}\sqrt{r}}}\frac{\partial }{\partial r}\left({\displaystyle {\sum }_{\text{gas}}\overline{\nu }\sqrt{r}}\right).$$(13)

In addition to advection, each dust particle of size a diffuses according to the concentration gradient, with a dust radial diffusivity approximated by (Youdin & Lithwick 2007) $$D_{\text{dust}}\left(r,a\right)=\frac{\overline{\nu }\left(r\right)}{1+{\text{St}}^2\left(r,a\right)}.$$(14)

Given that both the dust radial and vertical diffusion coefficients are determined by Eq. (14), where the latter results in Eq. (5), we implicitly assume the turbulence to be isotropic.

The transport for each grain of size a is described by the following 1D radial advection-diffusion equation, for the dust surface density Σ_dust(r, a) (Birnstiel et al. 2010): $$\frac{\partial {\sum }_{\text{dust}}}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left\{r\times \left[{\displaystyle {\sum }_{\text{dust}}{\upsilon }_{\text{r,dust}}-{\displaystyle {\sum }_{\text{gas}}{D}_{\text{dust}}}\frac{\partial }{\partial r}}\left(\frac{{\sum }_{\text{dust}}}{{\sum }_{\text{gas}}}\right)\right]\right\}=0.$$(15)

There are no sink terms included in Eq. (15) because we ignore the potential loss of dust particles through wind entrain-ment from, for example, internal photoevaporative winds or MHD disk winds (see e.g., Gárate et al. 2021; Rodenkirch & Dullemond 2022).

In order for the dust evolution model to be complete, grain coagulation (growth) and fragmentation needs to be included. Indeed, each grain of the dust distribution can be transported across the protoplanetary disk, and their size can also evolve through sticking and fragmentation (Birnstiel et al. 2010). All of these processes are included in DustPy by solving the Smoluchowski coagulation equation (Smoluchowski 1916), simultaneously with the transport of the grains (Eq. (15)).

We note that all dust-related quantities (and $\overline{\alpha }$; see Sect. 2.3) become time-dependent when a dust evolution model is employed. Furthermore, f_dg,tot(r) = Σ_dust,tot(r)/Σ_gas(r) is no longer a free parameter for the models of this present paper using the dust evolution model (Models IV–VI). Instead, it is determined by the dust evolution calculations, through solving for Σ_dust(r, a) at each time step. The initial f_dg,tot (required to start the dust evolution model) is chosen to be radially constant equal to 10⁻².

2.3 MRI-driven turbulence model

In this paper, we use and improve the MRI-driven turbulence model of Delage et al. (2022). The main output of this model is an effective radial profile for the MRI-induced viscosity parameter, $\overline{\alpha }$, as shown by the flowchart in Fig. 1. The next paragraph summarizes the main ingredients.

Their model is a 1+1D global magnetically driven disk accretion model built to study the outer region of class II proto-planetary disks (r ≳ 1 au), which accretes viscously solely due to the MRI and hydrodynamic instabilities. It has the advantage to capture the essence of the MRI-driven accretion, without resorting to computationally expensive 3D global nonideal MHD simulations. It includes the key following physical processes: (1) disk heating by passively absorbing stellar irradiation; (2) dust settling; (3) nonthermal ionization from stellar X-rays and galactic cosmic rays as well as the decay of short- and long-lived radionuclides; and (4) ionization chemistry based on a semi-analytical chemical model that captures the charge state of the disk dust-gas mixture, hence carefully modeling the gas ionization degree. In order to know where the MRI can operate in the disk, the general methodology is to compute the magnetic diffusivities of the nonideal MHD effects as well as their corresponding Elsasser numbers from the ionization chemistry, and apply a set of conditions for sustaining active MRI derived from 3D numerical simulations. These conditions account for the suppression of the MRI by Ohmic resistivity and ambipo-lar diffusion, but ignore for now the role of the Hall effect. In the MRI-dead zones (where the MRI is suppressed), it is further assumed that the gas can still accrete due to a small constant hydrodynamic turbulent parameter α_hydro, induced by hydrodynamic instabilities, such as VSI (e.g., Flock et al. 2020; Barraza-Alfaro et al. 2021). For given stellar, gas and dust properties, the Shakura-Sunyaev viscosity parameter, α, can thus be determined self-consistently under the framework of viscously driven accretion from detailed considerations of the MRI with Ohmic resistivity and ambipolar diffusion. It is computed both as a function of radius and height, eventually leading to the effective turbulent parameter ᾱ, which is the key output quantity for further coupling with ID gas and dust evolution models. In this work, the rms magnetic field strength (B) is numerically constrained by our MRI-driven turbulence model, and chosen such that the MRI activity is at the maximal efficiency as permitted by the two nonideal MHD effects considered. In other words, В is found such that it maximizes the accretion rate in the MRI active region at any radii (Sect. 3.2 of Delage et al. 2022). We note that the accretion driven by magnetic disk winds is ignored to solely focus on the roles of the MRI and hydrodynamic instabilities.

In Delage et al. (2022), the authors assumed a mono-disperse dust distribution of fixed size to describe the dust phase. To improve their model with a dust size distribution, we now define the following three quantities, at any locations (r,z) in the protoplanetary disk: the representative grain size a_dust,rep, the representative grain cross-section σ_dust,rep, and the representative grain mass m_dust,rep. They read as follows: $$a_{\text{dust,rep}}\left(r,z\right)=\frac{1}{n_{\text{dust,tot}}}{\displaystyle \sum _{a}a{n}_{\text{dust}}\left(r,z,a\right),}$$(16) $${\sigma }_{\text{dust,rep}}\left(r,z\right)=\frac{1}{n_{\text{dust,tot}}}{\displaystyle \sum _a\pi a^2{n}_{\text{dust}}\left(r,z,a\right),}$$(17) and $$\begin{array}{l}{m}_{\text{dust,rep}}\left(r,z\right)=\frac{1}{{\rho }_{\text{dust,tot}}}{\displaystyle \sum _a\frac{4}{3}\pi {\rho }_{\text{bulk}}{a}^3{\rho }_{\text{dust}}\left(r,z,a\right).}\\ \end{array}$$(18)

In practice, ρ_dust, n_dust, a_dust, σ_dust, m_dust defined in Delage et al. (2022) must now be replaced by ρ_dust,tot, n_dust,tot, a_dust,rep, σ_dust,rep, m_dust,rep, respectively. Indeed, these new five dust quantities encapsulate all the necessary information to successfully implement a dust size distribution in their model, hence making possible the coupling between dust evolution and MRI-driven turbulence calculations. We note that the definition of m_dust,rep does not really matter because it only intervenes through the dust Hall parameter, which is actually quasi-independent of it since the grain mass is always much larger than m_neutral for any grains of size a ≥ 0.1 µm (see Eq. (21) of Wardle 2007). Furthermore, we note that the representative grain size can be thought of the total grain size per unit volume called “C_tot” in Okuzumi et al. (2011a,b), Ormel & Okuzumi (2013), whereas the representative grain cross-section can be seen as their total grain surface area called “A_tot” We can write a_dust,rep = C_tot/n_{dust, tot} and σ_dust,rep = A_tot/n_dust,tot. For the ionization chemistry, what actually matters are the quantities C_tot and A_tot. We note that the dust mainly affects the ionization chemistry (hence the MRI-driven turbulence) by A_tot, and weakly by C_tot.

Except for the grain size (here we assume a dust distribution with different sizes rather than a mono-disperse distribution) and the vertically integrated total dust-to-gas-mass ratio f_dg,tot (see Sect. 2.2.2), the parameters used to run our MRI-driven turbulence model are either taken from Table 1 of Delage et al. (2022), if not explicitly mentioned, or as follows: M_⋆ = 1 M_⊙, L_⋆ = 2 L_⊙, M_disk = 0.05 M_⋆, α_hydro = 10⁻⁴, and s_e = 0.6. We note that the electron sticking coefficient, s_e, is chosen equal to 0.6 (instead of 0.3 as in their paper). This updated value is more compatible with the detailed calculations conducted by Bai (2011a).

Fig. 1 Flowchart of the MRI-driven disk accretion model presented in Delage et al. (2022). This model captures the essence of the MRI-driven turbulence in a 1+1D framework, accounting for the following: stellar properties (gray symbols), disk gas properties (blue symbols), disk dust properties (red symbols), nonthermal ionization sources (yellow symbols), ionization chemistry modeling the gas ionization degree (green symbols), disk magnetization properties (powder blue symbols), and nonideal MHD calculations (dark pink symbols). The main output of the model is an effective radial profile for the MRI-induced viscosity parameter, ᾱ (Eq. (1)). In this paper, we improve the dust phase modeling with a dust size distribution, either by assuming a fixed power-law distribution with different properties or the outputs from dust evolution obtained with DustPy (see Table 1 and text in Sect. 2.3 for further explanations).

3 Simulation setups

Table 1 summarizes all the models considered in the present paper. Below we describe them in details, and explain the numerical implementation.

3.1 Models I–III

Models I–III investigate how a fixed power-law dust size distribution with varying parameters impacts the MRI activity, especially the steady-state MRI-driven accretion described in Delage et al. (2022). It is determined by the three parameters a_min, a_dist,Max, and p_dist,Exp (see Eq. (11)), and its effect is studied by varying them one at a time. We refer to such a dust distribution as “Power-law” in Table 1. Here we note that providing the fragmentation velocity, υ_frag, is irrelevant because there is no dust evolution in this case.

For these models, Σ_gas is computed alongside $\overline{\alpha }$ using our MRI-driven turbulence model (Sect. 2.3) in the following mode: we seek for a steady-state MRI-driven accretion for the gas (the gas accretion rate is radially constant, implying that all regions of the disk accrete at the same rate), corresponding to the given dust and stellar properties. Through an iterative process (see Sect. 4.2 of Delage et al. 2022), $\overline{\alpha }$ and Σ_gas are computed together in order for the gas to follow such a regime. We refer to it as “Steady-State 1 (SS1)” in Table 1.

3.2 Models IV–VI

Models IV–VI investigate the impact of dust evolution on the MRI-driven turbulence (there is no gas evolution here). We particularly monitor how the effective MRI-induced turbulent parameter, $\overline{\alpha }$, varies as a function of various dust evolution snapshots. To do so we partially couple the dust evolution model employed with DustPy to our MRI-driven turbulence model. We emphasize that this coupling is not self-consistent yet, since we have not treated the evolution of dust and $\overline{\alpha }$ simultaneously: At each dust evolution snapshot, the new radial profile of $\overline{\alpha }$ derived from our MRI-driven turbulence model is not reinjected into the dust evolution model. This implies that the feedback from a change in the turbulence level due to dust evolution is not accounted for the next steps of the dust evolution calculations (although, a change in $\overline{\alpha }$ is expected to impact both the gas and the dust). We justify our choice by reminding that the framework of this present work aims to isolate the effect of dust evolution on the MRI-driven accretion. In this context, our goal is only to investigate whether dust evolution can change the MRI-induced α over time at all, as well as where in the protoplanetary disk substantial changes can occur, if any. Besides, we note that a full self-consistent coupling between MRI-driven turbulence calculations, gas and dust evolution will be addressed in a future paper focusing on the potential dust trapping power of the dead zone outer edge.

In this paper, the methodology for partially coupling DustPy with our MRI-driven turbulence model is the following: At each disk radius, we assumed the grain size distribution to initially follow a MRN-like distribution of interstellar grains (Mathis et al. 1977), with a_min = 0.55 µm, a_dist,Max = 1 µm and p_dist,Exp = −3.5. We then let the dust phase to evolve until 5 Myr (referred as “DustPy” in Table 1). For various dust evolution snapshots, we then computed the five dust quantities (Eqs. (8), (10), (16), (17), (18)), which are used as inputs into our MRI-driven turbulence model. The corresponding ᾱ can thus be computed for each dust evolution snapshot. Doing so means that we employ our MRI-driven turbulence model following its other mode where the steady-state accretion assumption is relaxed: Computing, on the fly, the effective MRI-induced turbulent parameter, $\overline{\alpha }$, corresponding to any provided gas (e.g., Σ_gas(r), dust (e.g., Σ_dust(r, a)) and stellar (e.g., the stellar X-ray luminosity L_XR) properties.

Since gas evolution is turned off, we need to provide the gas surface density, Σ_gas, that is used for both the dust evolution and the MRI-driven turbulence model. It is fixed to an input profile, which follows what we refer to as either “Steady-State 2 (SS2)” or “LBP” in Table 1 (see the profiles in Panel b of Figs. 4 and 6, respectively). On the one hand, SS2 means that Σ_gas follows the gas surface density profile obtained for the steady-state MRI-driven accretion solution corresponding to Model II, with the same MRN-like grain size distribution described in the previous paragraph, and after applying Rayleigh adjustment to it (see Appendix A). On the other hand, LBP means that follows the classical radial profile of a power law combined with an exponential cutoff (self-similar solution, Lynden-Bell & Pringle 1974) with a total disk gas mass M_disk = 0.05 M_⋆, and a critical radius R_c = 80 au. Additionally, we need to provide the $\overline{\alpha }$ used to perform the dust evolution. We fixed it to an input profile (see at t = 0 yr in Panel a of Figs. 4, 6 or 8) obtained by our MRI-driven turbulence model, assuming the same MRN-like grain size distribution, and the same gas surface density profile (following either condition SS2 or LBP) as described above.

In the light of the recent laboratory experiments on icy particles (Gundlach et al. 2018; Musiolik & Wurm 2019; Steinpilz et al. 2019), we assume the fragmentation velocity, υ_frag, to be radially constant equal to 1 m s⁻¹ in Models IV and V. For completeness, we also adopt the higher value of 10 m s⁻¹ (Wada et al. 2011; Gundlach et al. 2011; Gundlach & Blum 2015) in Model VI.

3.3 Numerical implementation

In all our simulations, the radial grid is computed from r_min = 1 au to r_max = 200 au, with N_r = 256 cells logarithmically spaced. For every radial grid-point r є [r_min, r_max], the corresponding vertical grid is computed from the disk midplane (z = 0) to z_max(r) = 5 H_gas(r), with N_z = 512 cells linearly spaced. We note that the vertical grid is only used to run our 1+1D MRI-driven turbulence model.

Regarding the grain size distribution, we always consider a logarithmic grid of grain species whose size range from a_min (free parameter depending on the model chosen) to 250 cm (fixed across the models). Furthermore, we consider seven mass bins per mass decade (choice based on the work of Ohtsuki et al. 1990; Drazkowska et al. 2014). For example, the total number of mass bins becomes N_m = 141 for a_min = 0.55 µm.

When the dust evolution model is employed, we particularly need to set the dust outer boundary condition. In Model IV, we impose at each time step a constant power-law on the dust surface density, leading to an inflow of dust particles over time consistent with the steady-state MRI-driven accretion assumption for the gas. Indeed, assuming that the gas is in steady-state MRI-driven accretion means that there is a gas and dust reservoir outside the simulation domain due to viscous spreading. In Models V and VI, we impose a floor value on the dust surface density, which essentially prevents the inflow of dust particles from the disk outer regions.

In the dust evolution code DustPy, there are δ parameters that control the turbulent collision velocities, vertical stirring, and radial diffusion of dust particles. Similar to the documentation, we use the symbols δ_turb, δ_vert, and δ_rad, respectively, for these physical processes. It is up to the user to decide whether such parameters are independent of each other, and depend or not on the α parameter that regulates the gas viscous evolution. In our simulations, we assume the mixing of dust particles to be driven by the MRI and hydrodynamic instabilities captured in our viscous parameter $\overline{\alpha }$. In other words, we assume that ${\delta }_{\text{turb}}={\delta }_{\text{vert}}={\delta }_{\text{rad}}=\overline{\alpha }$. Finally, we use DustPy version 0.5.8.

4 Results

4.1 The effect of a dust size distribution on the steady-state MRI-driven accretion

To study the effect of a dust size distribution on the steady-state MRI-driven accretion investigated in Delage et al. (2022), we run a set of simulations where we assume the dust size distribution to follow a fixed power-law described by a_min, a_dist,Max, and p_diSt,Exp (see Eq. (11)). In the following, we present the impact of a variation in a_dist,Max, while fixing a_min = 0.1 µm and p_dist,Exp = −3.5. This set of simulations corresponds to Model I, and the results are presented in Fig. 2.

Figures 2a, e and f show that a higher adist,Max leads to stronger MRI-driven turbulence overall, a higher gas accretion rate, and a more compact dead zone. When larger grain sizes are included in the dust distribution, the overall ionization level becomes higher (Fig. 2c), leading to enough charged particles in the gas phase for the magnetic field to couple with, mainly due to two reasons: (1) Dust settling becomes more important, which locally leads to an increase in the dust-to-gas mass ratio at the midplane, hence promoting the ionization power of radionu-clides (dominating the ionization process in the inner regions of the dead zone). (2) The total grain surface area, Atot, decreases (i.e., grains can less efficiently adsorb free electrons or ions onto their surfaces), resulting in the gas-phase recombination more easily dominating the recombination process over grain surface adsorption. A_tot decreases because a fraction of the smaller sizes in the dust distribution is replaced by larger sizes, implying that charged particles in the gas phase encounter, per unit volume and on average, less small dust particles.

Consequently, the MRI can operate with stronger magnetic field strengths on average (Fig. 2d). Finding stronger B overall for increasing adist,Max means that ambipolar diffusion becomes less stringent, allowing for MRI-driven turbulence with stronger magnetic field strengths. Furthermore, the MRI can operate closer to the central star, leading the dead zone outer edge to be almost located as twice as close for a_dist,Max = 10 µm compared to a_dist,Max = 1 µm. Interestingly, we notice that for the given model parameters the dead zone outer edge is within 1050 au, and even 10–20 au for a_dist,Max ≥ 100 µm. The situation in which the midplane would be almost entirely MRI-dead (gray dashed line in Fig. 2a), obtained assuming a mono-disperse dust distribution of fixed size a_mono = 0.1 µm, is thus greatly mitigated when a dust distribution with different sizes is taken into account. We find that the presence of micron-sized particles in the dust size distribution, on top of the submicron-sized particles, prevents the dead zone from extending up to ~100 au for our choice of the magnetic field strength and configuration.

Figure 2a shows that the solid colored lines lie within the dashed gray and black lines, representing the two limiting scenarios for α in which all the dust would be either in the form of grains of size a_mono = 0.1 µm or a_mono = 1 cm, respectively. Since the minimum grain size of Model I is a_min = 0.1 µm and its maximum grain size is a_dist,Max = 1 cm, it is expected that the α obtained from such dust distributions with different grain sizes can neither represent a less MRI active scenario than the dashed gray line, nor a more MRI active scenario than the dashed black line. It is worth mentioning that the models with growth (solid colored lines) converge toward the grain-free scenario for the model parameters considered, corresponding to the dashed black line. The dashed black line mimics the grain-free scenario because the recombination process occurs in the gas phase rather than onto the grains surface when all grains are of size 1 cm. Indeed, the dashed black line of Fig. 2a is identical to the red and purple lines of Fig. 8a of Delage et al. (2022) as well as their blue and yellow lines of Fig. 9a, implying that this is the grain-free steady-state solution because it is independent of the dust properties considered.

From Fig. 2, we can infer that the presence of larger sizes in the dust distribution (due to dust growth) substantially impacts the MRI-driven turbulence. Particularly, we expect dust growth to have a major positive impact on the MRI activity in regions where grain surface absorption is the main process for recombination (this recombination regime is highly sensitive to the dust properties), whereas such impact is expected to be weak in regions where the recombination process is dominated by gas-phase recombination (this recombination regime is weakly dependent of the dust properties). Figure 2a even shows that the overall effective turbulence level in the inner regions of the dead zone ($\overline{\alpha }$ for r ≲ 10 au) is noticeably different depending on a_dist,Max. This suggests that dust growth is able to change the activity in the MRI active layer sitting above the dead zone such that the effective turbulence level in the dead zone increases. Furthermore, it seems that dust growth does not need to be very efficient to generate a significant boost in the MRI-driven turbulence. Indeed, the positive feedback obtained from the presence of larger grain sizes on the MRI activity is getting less noticeable once the maximum grain size is larger than 100 µm.

For example, the gas accretion rate and the dead zone outer edge do not change as much for a_dist,Max ≥ 100 µm as 1 µm ≤ adist,Max ≤ 100 µm (Figs. 2e and f).

Here, we have only presented the effect of the variation in a_dist,Max on the steady-state MRI-driven accretion. In Appendix B, we also show the results for the variation in the two remaining parameters a_min (Model II) and p_dist,Exp (Model III), respectively. The main conclusion is that increasing any of these three parameters leads to a decrease in the total grain surface area, A_tot; hence stronger MRI-driven turbulence overall, a higher accretion rate, and a more compact dead zone. In other words, we expect that any changes occurring in the dust size distribution due to the evolution of the dust phase should impact the MRI-driven turbulence mainly through the quantity A_tot.

Fig. 2 Impact on the steady-state MRI-driven accretion when varying the maximum grain size, adist,Max, of the fixed power-law dust size distribution (Model I). Going from solid light-colored to dark-colored lines, adist,Max spans the range from 1 µm to 1 cm. The panels show the steady-state radial profiles of several key quantities (see also Delage et al. 2022, for their definition), for the model parameters M⋆ = 1 M⊙, L⋆ = 2 L⊙, Mdisk = 0.05 M⋆, fdg,tot = 10−2, αhydro = 10−4, amin = 0.1 µm, and pdist,Exp = −3.5. Panel a: pressure-weighted vertically integrated turbulent parameter, $\overline{\alpha }$. Panel b: gas surface density, Σgas. Panel c: midplane total ionization rate, ζ. Panel d: midplane optimal rms magnetic field strength, B. Panel e: (constant) gas accretion rate, ${\dot{M}}_{\text{acc,gas}}$. Panel f : midplane radial dead zone outer edge location, RDZ. For comparison, the dashed gray and black lines in panel a display the steady-state quantity $\overline{\alpha }$ obtained assuming the limiting case of a mono-disperse population of dust with size amono = 0.1 µm and amono = 1 cm (equivalent to the grain-free scenario), respectively.

4.2 The effect of dust evolution on the MRI-driven turbulence

The dust size distribution is expected to deviate from a power-law as the dust phase evolves, even in the gaps and pressure maxima (see e.g., Andama et al. 2022). Consequently, we run two dust evolution simulations that we partially couple to our MRI-driven turbulence model, relaxing the assumption of the steady-state accretion, in order to investigate the effect of dust evolution on the MRI activity. These two differ in how the gas surface density profile was chosen (either following the condition SS2 or LBP, as explained in Sect. 3.2), as well as in the dust outer boundary condition we adopt.

4.2.1 SS2 condition for the gas surface density profile

This simulation assumes that the gas surface density profile is stationary, following the steady-state MRI-driven accretion solution corresponding to Model II with a MRN-like dust size distribution (a_min = 0.55 µm, a_dist,Max = 1 µm and p_dist,Exp = −3.5; see Appendix B.1), and after applying Rayleigh adjustment to it (see Appendix A). The corresponding $\overline{\alpha }$ computed with our MRI-driven turbulence model is used to perform the dust evolution. From the dust phase perspective, therefore, the dead zone is invariant over time, with its outer edge always located at ~27 au.

Figure 3 shows that the dust can grow up to sizes larger than millimeters inside the dead zone (r ≲ 27 au), whereas the maximum grain size reaches an upper limit of ~10 µm outside of it, with an abrupt transition at the dead zone outer edge. This is because the maximum grain size can be limited by either drift when the drift timescale exceeds the growth timescale, or fragmentation when the collision velocity between dust particles exceeds the material fragmentation velocity υ_frag (Brauer et al. 2008; Birnstiel et al. 2009). For each limit, the time-dependent maximum Stokes number reachable by a dust particle is determined by $${\text{St}}_{\text{drift}}={\left|\frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}{P}_{\text{gas,mid}}}{\text{d\hspace{0.17em}}\ln r}\right|}^{-1}\left(\frac{{\upsilon }_{\text{K}}}{c_{\text{s}}}\right)f_{\text{dg,tot}},$$(19)

and $${\text{St}}_{\text{frag}}=\frac{1}{3\overline{\alpha }}{\left(\frac{{\upsilon }_{\text{frag}}}{c_{\text{s}}}\right)}^2.$$(20)

The $\overline{\alpha }$ used to perform the dust evolution is shown at t = 0 yr in Fig. 4a. It has a mean value of ~1.7 × 10⁻⁴ in the dead zone, and ~3 × 10⁻³ in the MRI active region, with a sharp increase near the dead zone outer edge. Additionally, Fig. 3 indicates that the maximum grain size is set by the fragmentation barrier (dashed cyan line), everywhere in the disk for the entire dust evolution simulation. Since St_frag is inversely proportional to α, and the turbulent collision velocity between grains is proportional to it with $\text{Δ}{\upsilon }_{\text{turb}}\approx \sqrt{3\overline{\alpha }/\text{St}{c}_{\text{s}}}$ (for St ≪ 1, Ormel & Cuzzi 2007), dust particles can grow into larger sizes in the dead zone compared to the MRI active region.

The turbulence level used to perform the dust evolution is so low in the regions within 30 au that grain coagulation is effective during the first stages of dust evolution (from 0 yr to 0.1 Myr), even leading to a depletion in the submicron-sized particles in some parts of the dead zone (top row of Fig. 3). We can better appreciate this depletion by looking at Figs. 5a and b. These panels show the temporal evolution of the midplane representative grain size (tracing the smallest sizes of the dust distribution), and the midplane equivalent size of the representative grain cross-section (relevant size of the dust distribution involved in the ionization chemistry), respectively. There is a growth wave propagating inside-out from the initial time until 0.01 Myr, which is attributed to the location where the dust size distribution becomes skewed toward larger sizes because submicron-sized particles grow much quicker than fragmentation can replenish them, as shown by Fig. 5c. This panel displays the temporal evolution of the midplane equivalent size of the representative grain mass (tracing the largest sizes of the dust distribution) in the protoplanetary disk. Since this quantity increases quickly in the dead zone within 0.01 Myr, it implies that larger sizes indeed become present in the dust size distribution. The effective growth of submicron-sized particles into larger sizes can also be seen by comparing Figs. 5d and e. While the midplane total number dust density (n_dust,tot) decreases within 0.1 Myr, the midplane total dust density (ρ_dust,tot) increases. Since the total dust content almost remains constant within this period of time, a decrease in n_dust,tot is primarily related to a decrease in the relative proportion of small dust particles in the dust distributions that have grown into larger sizes (hence the increase of ρ_dust,tot)·

The first direct consequence of the initial effective grain growth in the dead zone, within 0.1 Myr of dust evolution, is that more dust particles can settle toward the midplane. Settling causes an increase in the local dust-to-gas mass ratio at the midplane (ρ_gas is constant here, while ρ_dust,tot increases as shown by Fig. 5d). Consequently, the ionization power of radionuclides increases, leading the midplane total ionization rate to increase for r ≲ 30 au (Fig. 4c). The second direct consequence is that ambipolar diffusion becomes weaker (particularly where grain surface adsorption dominates the recombination process), since the total grain surface area, A_tot, decreases when the dust grow (see Sect. 4.1) or when the dust distribution is skewed toward larger grain sizes (see Appendix B.2). In this regard, Fig. 4d shows that the MRI is allowed by ambipolar diffusion to have stronger magnetic field strengths for r ≲ 30 au, as dust evolves from 0 yr to 0.1 Myr, particularly near the dead zone outer edge.

From these two consequences, we can understand the temporal evolution of $\overline{\alpha }$ for r ≲ 30 au, within the first stages of dust evolution (from 0 yr to 0.1 Myr): A higher ionization level combined with a less stringent ambipolar diffusion result in the MRI activity being able to operate closer to the central star (Fig. 4f shows that the dead zone outer edge deceases within 0.1 Myr), with stronger turbulence generated for any regions within 30 au (Fig. 4a). Interestingly, we notice that ᾱ first varies in the very inner regions of the dead zone due to an increase in the radionuclide ionization rate from dust settling, followed by a variation in the outer regions of the dead zone due to a decrease in the total grain surface area, A_tot, from dust growth. Looking at Fig. 3, we find that dust particles reach their maximum sizes between 0.01 Myr and 0.1 Myr, which coincides with the rimescale over which the MRI activity has substantially changed (see Figs. 4a and f). This indicates that the timescale over which dust evolution significantly impacts the MRI-driven turbulence is a timescale of local dust growth (see e.g., Eq. (30) of Birnstiel et al. 2016, for the analytic formula).

On another note, we notice that ᾱ continues to increase while R_DZ decreases between 0.1 Myr and 1 Myr (Figs. 4a and f). Since the maximum grain sizes have already been reached within 0.1 Myr (Fig. 3), the explanation for such behaviors no longer lies in dust growth alone. Figure 5d shows that ρ_dust,tot significantly decreases in the regions within ~30au, between 0.1 Myr and 1 Myr. We can infer that the dust content decreases in these regions, due to radial drift which is faster for particles that have grown to larger sizes. When the dust is removed from the disk, the dust-to-gas mass ratio decreases, which allows the MRI to operate closer to the central star with stronger activity (see Sect. 6.3 of Delage et al. 2022). Indeed, although a decrease in the dust-to-gas mass ratio implies a lower mid-plane total ionization rate for r ≲ 30 au due to less radionuclides (Fig. 4c), the dust is far less efficient at sweeping up free electrons and ions from the gas phase due to a decrease in the total grain surface area, hence resulting in a net increase in the gas ionization degree.

Despite some dust growth within the first stages of dust evolution for r ≳ 30 au, we notice that it is not enough to make a real impact on the temporal evolution of the MRI-induced effective turbulent parameter, the midplane total ionization rate or the optimal rms magnetic field strength (Figs. 4a, c, and d). It is because the dominant recombination process in these regions is gas-phase recombination, which is weakly dependent on the dust properties. Consequently, it is expected not to see much change in terms of MRI activity for r ≳ 30 au, within 1 Myr.

Finally, we need to focus on the late stages of dust evolution (t > 1 Myr), where all the five dust quantities used to couple DustPy with our MRI-driven accretion model reach a quasi-steady-state (Fig. 5). In the present model, we assumed that the gas surface density is fixed to a steady-state profile (Fig. 4b). While the steady-state profile describes the inner inward accreting region (r ≲ R_t) of the viscously evolving disk, it cannot capture the outer viscously expanding region (r ≳ R_t). Here R_t corresponds to the transition radius at which the gas motions changes from inward to outward in a viscously evolving disk (e.g., Hartmann et al. 1998). We thus expect the presence of an outer region outside our simulation domain (r ≳ 200 au), which can feed the disk (1 ≲ r ≲ 200 au) of submicron- and micron-sized dust particles. Here we aimed to mimic this situation by choosing the dust outer boundary condition such that there is an inflow of (small) dust particles through the outer boundary of the domain. Consequently, the fact that the dust reaches a quasi-steady-state for t > 1 Myr appears to be directly linked to such a choice. It leads the MRI-driven turbulence and the dead zone outer edge to be roughly constant from 1 Myr all the way until 5 Myr of dust evolution. Particularly, the apparent temporal dead zone stability whilst the dust evolves is thus a mere artifact of our choice for the dust boundary condition in this case. It is worth noting that the radioactive decay of²⁶ Al has been ignored in the present paper, which means that the radionuclide ionization rate for t > 1 Myr is expected to be smaller than what we assumed (see Eq. (15) of Delage et al. 2022). However, it is not too bothersome because (Delage et al. 2022, see, e.g., their Fig. C.l) showed that the ionization from radionuclides only dominates the total ionization rate deep within the dead zone (radially and vertically). Consequently, the location of the dead zone outer edge should not significantly vary, in the cases considered, if we account for the decrease over time in the radionuclide ionization rate due to the decay of²⁶ Al. This statement also holds for the time evolution for t > 1 Myr seen in the next section.

Although the coupling between dust evolution and our MRI-driven turbulence model is only partial in this paper (see Sect. 3.2), it is clear that dust evolution has significant impact on the MRI activity: dust settling, grain coagulation, and fragmentation drive the change in the gas ionization degree, hence on ᾱ at each radius r, on a timescale of local dust growth. In the specific case of Model IV, it is thus expected that dust evolution drives the gas away from the assumed steady-state MRI-driven accretion (condition SS2 for the gas surface density profile), in the regions within ~30 au, by substantially changing the effective turbulent parameter $\overline{\alpha }$ within 0.1 Myr. Indeed, any changes in ᾱ due to dust evolution directly modifies the 1D advection-diffusion equation for the gas. We expect this change to carry on for the next steps of the dust evolution calculations because it depends on $\overline{\alpha }$ and the gas properties. Consequently, we find that the steady-state MRI-driven accretion that we impose on the gas is actually not physically consistent for r ≲ 30 au.

Fig. 3 Temporal evolution of the dust surface density distribution per logarithmic bin of grain size σ(r, a) (see Eq. (C.1)), for Model IV. In each panel, the solid white line shows a Stokes number of unity (radial drift reaches its maximal efficiency), the dotted green line shows the drift limit (Eq. (19)), and the dashed cyan line shows the fragmentation limit (Eq. (20)). The horizontal solid dark gray lines show a grain size of 1 mm.

Fig. 4 Impact of dust evolution on the MRI-driven turbulence, for Model IV. The panels show the temporal evolution of the same quantities as in Fig. 2, except for the gas accretion rate. Also, the gas surface density is now fixed to the displayed input profile in panel b. We emphasize that these quantities do not describe steady-state MRI-driven accretion (unlike Fig. 2), since they are recalculated at each dust evolution snapshot, through partial coupling between the ID radial dust evolution model employed and our MRI-driven turbulence model. On this note, the corresponding temporal evolution of the five dust quantities used to perform such a coupling are shown in Fig. 5. Here we note that, in panel f, the dead zone outer edge coincides at t = 0 yr, t = 100 yr and t = 1000 yr.

Fig. 5 Temporal evolution of the five dust quantities used to couple the ID radial dust evolution model employed with our MRI-driven turbulence model, for Model IV. These quantities are the representative grain size adust,rep (Eq. (16)) in panel a, the representative grain cross-section σdust,rep (Eq. (17)), the representative grain mass mdust,rep (Eq. (18)), the total dust density ρdust,tot (Eq. (8)) in panel d, and the total dust number density ndust,tot (Eq. (10)) in panel e. Instead of displaying σdust,rep and mdust,rep, we show the equivalent grain size $a_{\text{dust,}\sigma }{\left({\sigma }_{\text{dust,rep}}/\pi \right)}^{\frac{1}{2}}$ corresponding to σdust,rep (panel b), and the equivalent grain size $a_{\text{dust,m}}={\left(3m_{\text{dust,rep}}/4\pi {\rho }_{\text{bulk}}\right)}^{\frac{1}{3}}$ corresponding to mdust,rep (panel c). These equivalent grain sizes better indicate the dominant size when considering the grain cross-section or the grain mass, each of which is important for the MRI calculations (see Sect. 2.3). The panels particularly show the midplane radial profiles. We note that adust,rep and adust,m traces the smallest and largest sizes of the dust distribution, respectively, whereas adust,σ is the relevant size of the dust distribution involved in the ionization chemistry.

Fig. 6 Same as in Fig. 4, except for Model V. The corresponding temporal evolution of the five dust quantities is shown in Fig. 7. The gas surface density (panel b) is fixed to the input profile corresponding to the classical self-similar solution (Lynden-Bell & Pringle 1974), with a total disk gas mass Mdisk = 0.05 M⋆, and a critical radius Rc = 80 au. As a result, the gas is no longer in a steady-state MRI-driven accretion at t = 0 yr, unlike Model IV. Here we note that, in panel f, the dead zone outer edge coincides at t = 0 yr and t = 100 yr.

4.2.2 LBP condition for the gas surface density profile

Here the goal is to see whether the results drawn in the previous section still hold with a different set of assumptions for the gas surface density profile and the dust outer boundary condition. In this case, we assume that the gas surface density profile follows the classical self-similar solution with a total disk gas mass M_disk = 0.05 M_⋆, and a critical radius R_c = 80 au (Fig. 6b). Here we emphasize again that the gas surface density profile does not evolve with time, so that we can solely focus on the effect of dust evolution on the MRI-driven turbulence. The dust outer boundary condition is chosen such that there is no inflow of small dust particles in the simulation domain, unlike Model IV. Additionally, the $\overline{\alpha }$ used to perform the dust evolution is now derived from the gas surface density profile mentioned above by our MRI-driven turbulence model (see Fig. 6a at t = 0 yr). From the dust perspective, it means that the dead zone is still invariant over time, but with its outer edge now always located at ~40 au.

Similar to what we find in the previous section, dust growth is fragmentation-limited everywhere in the disk and throughout its whole evolution, with dust particles growing into larger sizes in the dead zone (r ⪝ 40 au) compared to the MRI active region (see Fig. C.1). The main difference, though, is that the particles do not reach sizes as large in the dead zone (there are less grains of size a ≥ 1 mm than Model IV), but do reach larger sizes in the MRI active region (grains can be as large as ~100 µm), resulting in a much smoother transition in the dust surface density per logarithmic bin of grain size (σ(r, a)) at the dead zone outer edge. This can be explained by the fact that the $\overline{\alpha }$ used to perform the dust evolution in Model V is on average higher in the dead zone and lower in the MRI active region than the one used in Model IV, with a much smoother transition at the dead zone outer edge.

Another common feature between the results of this case (Model V) and Model IV is the effective grain coagulation in the regions within 40 au, during the first stages of dust evolution (from 0yr to 0.1 Myr). Specifically, it results in a similar growth wave attributed to the depletion of small dust particles that grow into larger sizes (Figs. 7a–c). In the same fashion as in Sect. 4.2.1, effective grain coagulation leads to: (1) more settling of dust particles toward the midplane (ρ_dust,tot increases as shown in Fig. 7d), allowing for a higher midplane ionization rate in the regions within 40 au (Fig. 6c); (2) less stringent ambipolar diffusion, implying that the MRI can operate with stronger magnetic field strengths (B increases as shown by Fig. 6d). Consequently, the MRI activity is able to operate closer to the central star as the dust phase evolves (Fig. 6f shows that R_DZ decreases over time), with stronger turbulence generated in the regions within 40 au (Fig. 6a). Similar to Model IV, $\overline{\alpha }$ (hence the dead zone outer edge) undergoes significant variation within 0.1 Myr, which corresponds to the timescale over which dust particles have reached their maximum size (Fig. C.1). This again suggests that the timescale over which dust evolution impacts the MRI-driven turbulence is determined by the timescale of local dust growth.

Comparing Figs. 6a, f with Figs. 4a, f, respectively, we notice that the temporal evolution of ᾱ and R_DZ between 0.1 Myr and 1 Myr is less pronounced in Model V compared to Model IV. Indeed, the dust particles within 40 au grow into smaller sizes compared to Model IV, which makes their radial drift slower. As a result, they can be present in the disk for a longer period of time, meaning that the removal of the dust content is delayed: Fig. 7d shows that ρ_dust,tot still increases between 0.1 Myr and 1 Myr and only starts decreasing from 1 Myr, while it decreases from 0.1 Myr in the case of Model IV as shown by Fig. 5d.

In regions of the disk beyond ~40 au (MRI active region from the dust perspective in Model V), Fig. 6 shows that, within 1 Myr, the effective MRI-induced turbulent parameter, the midplane total ionization rate, and the optimal rms magnetic field strength vary a bit more than in Model IV. We explain this behavior by noticing that the dust particles can grow into larger sizes in the MRI active region of Model V, especially near the dead zone outer edge (r ~ 40 au). This region marks the transition for the recombination process between grain surface adsorption and gas-phase, and is therefore more sensitive to the dust properties compared to the outer regions. For r ≳ 80 au, though, we can see that the MRI activity is pretty much steady within 1 Myr. This is expected since R_c = 80 au corresponds to the critical radius where the gas surface density profile drops exponentially, and therefore where the gas-phase is the main channel for recombination due to a high gas ionization degree.

To complete the comparison between Model IV and Model V, we now need to discuss the temporal evolution of the MRI-driven turbulence during the late stages of dust evolution (t > 1 Myr). In the previous section, we saw that $\overline{\alpha }$ and R_DZ were roughly steady after 1 Myr of dust evolution, since the five dust quantities used to couple DustPy with our MRI-driven accretion model reach a quasi-steady-state due to the steady influx of dust from the outer boundary. In Model V, though, ρ_dust,tot and n_dust,tot keep decreasing from 1 Myr all the way until 5 Myr (Figs. 7d and e), since the dust content is gradually removed by radial drift which is no longer compensated for by an inflow of small particles in the simulation domain as in Model IV. This implies that $\overline{\alpha }$ continues to increase while R_DZ decreases between 1 Myr and 3 Myr because the MRI can operate closer to the central star with stronger activity for decreasing dust-to-gas mass ratio (as explained in the previous section). Nonetheless, a salient result of Model V is the temporal evolution of the MRI activity between 3 Myr and 5 Myr. Although the dust keeps being removed from the disk, $\overline{\alpha }$ and R_DZ become stationary (Figs. 6a and f). It thus suggests that the MRI activity in the disk becomes weakly dependent on the dust properties when the dust content drops below a certain threshold, which is caused by a lack of dust particles to efficiently sweep up free electrons and ions from the gas phase. In other words, gas-phase recombination dominates in most of the protoplanetary disk after 3 Myr, and the dust has no longer a significant impact on the ionization chemistry. The disk dust-gas mixture thus behaves as a grain-free plasma after 3 Myr, and the dead zone outer edge becomes stationary because the MRI activity evolution becomes primarily controlled by the gas which is not evolving here. This result suggests that the dead zone may potentially be able to survive the protoplanetary disk evolution over a few million years when the MRI is the main driver for the disk accretion. We further discuss this idea in Sect. 5.2.

The results of this section emphasize that dust evolution has significant impact on the MRI activity, regardless of the assumptions made for the gas surface density profile or the dust outer boundary condition. Particularly, the MRI-induced ᾱ undergoes substantial change within the timescale over which the dust particles grow.

Fig. 7 Same as in Fig. 5, except for Model V.

5 Discussion

5.1 The effect of the fragmentation velocity

Laboratory experiments of particle collisions are crucial in understanding the growth of grains from interstellar medium micron-sized dust to millimeter- and centimeter-sized pebbles in protoplanetary disks. First, they demonstrate the potential outcome (sticking, fragmentation, bouncing or mass transfer) after grain collisions with a given initial relative velocity (see e.g., Windmark et al. 2012; Birnstiel et al. 2016). Second, they help to constrain the velocity threshold (fragmentation velocity υ_frag) for either effective growth or destructive collisions (e.g., Blum & Wurm 2000, 2008; Kimura et al. 2015, 2020). In the classical picture, it has been commonly thought that amorphous water-ice particles are stickier than silicates (Wada et al. 2011; Gundlach et al. 2011; Gundlach & Blum 2015). Theoretical models employing dust evolution thus usually adopt a fragmentation velocity of 10m s⁻¹ for ice particles. However, recent laboratory experiments show that this threshold velocity sensitively depends on the composition and temperature of the colliding particles. They found that water-ice particles may be as fragile as silicates, or even more so, resulting in the fragmentation velocity potentially being as low as 1 ms⁻¹ (Gundlach et al. 2018; Musiolik & Wurm 2019; Steinpilz et al. 2019). In Models IV and V, we therefore have experimented with υ_frag = l m s⁻¹. For completeness, Model VI investigates how this parameter impacts our results by taking a value of 10m s⁻¹. We note that, except for υ_frag, Model VI has the same setup as Model V.

For larger values of the fragmentation velocity (υ_frag), grain collisions are less destructive which improves the efficiency of coagulation. The fragmentation barrier thus becomes less stringent, allowing grains to grow into larger particles than Model V, everywhere in the disk (Fig. C.2). Within 0.1 Myr of dust evolution, the grains can grow into such large sizes that their Stokes numbers reach values close to unity. This results in the dust content being quickly removed from the disk between 0.1 Myr and 1 Myr due to effective radial drift, hence the quick increase in the MRI activity (Fig. 8). Indeed, Fig. 9d shows that ρ_dust,tot first increases from 0 yr to 0.1 Myr (caused by dust growth), then plummets by two orders of magnitude between 0.1 Myr and l Myr. As a result, the dust-to-gas mass ratio significantly decreases within l Myr (ρ_gas is stationary), implying that the drift barrier gradually becomes more stringent because St_drift decreases when f_dg,tot decreases (see Eq. (19)). The main consequence is that dust growth is no longer solely fragmentation-limited, unlike Model V. Instead, the temporal evolution of the dust surface density distribution per logarithmic bin of grain size (Fig. C.2) shows that it transitions from being fragmentation-limited to drift-limited in the entire protoplanetary disk, in a million-year timescale. The larger grains thus radially drift inward before they can collide and replenish the smaller ones, resulting in the depletion of the latter at a wide range of radii within the dead zone as seen by the dust particles (r ≲ 40 au), from 1 Myr. The lack of small dust particles leads the total dust number density to plummet from 0.1 Myr to t > 1 Myr by at least two order of magnitudes (Fig. 9e), and the grain cross-section representative to be skewed toward particles as large as 100 µm in the disk inner regions (Fig. 9b).

Since the rapid drift of the larger particles significantly reduces the dust content in the protoplanetary disk, as well as prevents the smaller ones to be replenished efficiently, the total grain surface area is significantly reduced after 0.1 Myr. Consequently, gas-phase recombination becomes the main channel after 0.1 Myr, which makes the MRI activity weakly dependent on the dust properties. Particularly, it leads the effective MRI-induced turbulent parameter, the dead zone outer edge, and the optimal rms magnetic field strength to become stationary (Figs. 8a, d, and f) because there is no gas evolution accounted for here. This needs to be put in the context of the stationary temporal evolution of the MRI-driven turbulence between 3 Myr and 5 Myr seen in Sect. 4.2.2. In the present model, such a stationary behavior occurs much earlier than in Model V because the dust content is removed on a much shorter timescale due to more effective radial drift. In other words, the disk dust-gas mixture behaves as a grain-free plasma much faster for the present model. Once the treatment of dust and MRI calculations is done simultaneously, we thus expect dust evolution to have a less long-term impact on the MRI activity in the regions of the disk with a higher fragmentation velocity.

Fig. 8 Same as in Fig. 6, except for Model VI. The corresponding temporal evolution of the five dust quantities is shown in Fig. 9. Here we note that, in panel f, the dead zone outer edge coincides at t = 0 yr and t = 100 yr.

Fig. 9 Same as in Fig. 5, except for Model VI.

5.2 The potential long-lived state of the dead zone in protoplanetary disks

One of our salient results is that, once the full self-consistent treatment of gas and dust evolution with MRI calculations is done, we expect the temporal evolution of the MRI-driven turbulence to be controlled first by dust evolution, then gas evolution. Indeed, we saw in Sects. 4.2.1 and 4.2.2 that the MRI-induced $\overline{\alpha }$ changes on a timescale of local dust growth, which is significantly shorter than the viscous evolution timescale in general. As long as there are enough dust particles in the disk to dominate the recombination process for the ionization chemistry, we thus expect the MRI activity evolution to be controlled by dust evolution. Once that is no longer the case, the MRI activity evolution is then expected to be controlled by gas evolution and occurs on a viscous evolution timescale, since the dust would no longer have a significant feedback on the ionization chemistry and the disk dust-gas mixture would behave as a grain-free plasma.

The timescale marking the transition from “dust-dominated” to “gas-dominated” MRI-driven turbulence depends on the physical properties of the protoplanetary disk. For instance, if the disk initially has a low dust content, the transition is expected to happen earlier in the disk lifetime. Conversely, if the disk initially has a high dust content or if it has pressure bumps (nonsmooth disk), this transition is expected to occur over a longer timescale. Regarding the latter, it is commonly accepted that pressure bumps are a possible explanation for the observed disk substructures, and can be created by various mechanisms such as embedded massive planets (see e.g., Zhu et al. 2012; Pinilla et al. 2012, 2015; Dong et al. 2015; Pyerin et al. 2021), dust particle growth by condensation near the ice lines (e.g., Stammler et al. 2017), or magnetic disk winds (e.g., Suriano et al. 2017, 2018, 2019; Hu et al. 2022). If there are pressure bumps in the disk, the main mechanism removing the dust (i.e., radial drift) is not as efficient because a substantial amount of particles can be trapped there. As a result, pressure bumps allow the dust to be in the disk for a longer period of time, hence delaying when the transition from dust-dominated to gas-dominated MRI-driven turbulence occurs.

In Sects. 4.2.2 and 5.1, we saw that the dead zone shrinks over time due to dust evolution (R_DZ decreases), and eventually becomes stationary after 3 Myr of disk evolution for Model V, and 0.1 Myr for Model VI. Unlike Model IV (Sect. 4.2.1), these models have no inflow of small particles feeding the outer boundary of the simulation domain. The stationary nature of the dead zone comes from the fact that the disk dust-gas mixture eventually behaves as a grain-free plasma after some time of evolution, where the MRI activity evolution is primarily controlled by the gas which is not evolving here. Indeed, Figs. 6f and 8f show that the midplane dead zone radial extent is larger than ~ 10 au, at any time and for both Models V and VI. Consequently, our results show that dust evolution alone does not lead to a complete reactivation of the dead zone in protoplanetary disks. As long as there are no mechanisms that can efficiently ionize the gas in the regions of the disk where the dead zone sits at and that the gas dispersal in those regions occurs in timescales of a few million years, the dead zone may potentially be able to survive the protoplanetary disk evolution over a few millions years when the MRI is the main driver for the disk accretion. This supports that a disk evolution model including X-ray photoevap-orative dispersal and a dead zone in the inner regions is a feasible idea in order to successfully explain the main observable properties of transitions disks such as extended gaps and high accretion rates (see Gárate et al. 2021).

5.3 The effect of the dust distribution minimum grain size

In the models where we partially coupled dust evolution with MRI calculations (Models IV, V and VI), we used a dust distribution minimum grain size, a_min, of 0.55 µm. This quantity is a free parameter in dust evolution models. From the modeling perspective, the choice of a_min should not matter because the dust particles quickly forget their initial grain size due to coagulation and fragmentation. As a result, running a pure dust evolution simulation with a_min = 0.1 µm or a_min = 0.55 µm should not substantially change the dust density output.

However, this is no longer true if MRI calculations are now combined with the dust evolution model. As we saw in Appendix B.1, taking either a_min = 0.1 µm or a_min = 0.55 µm for the distribution minimum grain size leads to appreciably different outcomes in terms of the MRI activity (particularly, the dead zone morphology and the location of its outer edge). Indeed, if a larger grain size is used for the distribution minimum size, the representative grain size involved in the ionization chemistry is skewed toward larger sizes, implying that the total grain surface area, A_tot, decreases and the overall MRI activity increases (see Figs. B.1a and f). In the context of a full coupling between dust evolution and MRI calculations, the MRI-driven turbulence would be sensitive to the dust distribution and its minimum grain size at any point in time (except if it is gas-dominated as discussed in the previous section). The choice for the adopted value of amin is thus crucial.

Tazaki & Dominik (2022) recently investigated the effect of monomer size and composition on scattering polarization of dust particles by using an exact light scattering technique. By comparing their simulations to observations, they estimated the monomer radius of dust particles to be no greater than 0.4 µm for several protoplanetary disks. They even found that a minimum grain size of 0.1–0.2 µm appears to explain the recent polarimet-ric observations of the disk around HD 142527. Nevertheless, they have not excluded the possibility that the monomers could actually be much smaller than 0.1 µm. If that were the case, polycyclic aromatic hydrocarbon (PAH), representing the smallest end of a grain size distribution, would need to be considered in the MRI calculations. Counter-intuitively, though, it has been shown by various works that including PAHs in the dust size distribution reduces ambipolar diffusion, hence enhancing the overall MRI activity (e.g., Bai 2011b; Zhao et al. 2016; Marchand et al. 2020). Consequently, if one wants to accurately describe the MRI-driven turbulence in protoplanetary disks, further studies are required with the aim to provide realistic constraints on the minimum grain size of the dust distribution.

6 Summary and conclusions

In this pilot study, we provide an important step toward a better understanding of the MRI-dust coevolution in protoplanetary disks, with the aim to present a proof of concept that dust evolution ultimately plays a crucial role in the MRI activity. To this end, we divided our analysis into two parts: First, we studied how a fixed power-law dust size distribution with varying parameters impacts the MRI activity, especially the steady-state MRI-driven accretion described in Delage et al. (2022), by employing and improving their 1+1D MRI-driven turbulence model. Second, we relaxed the steady-state accretion assumption in this newly improved turbulence model, and partially coupled it to the dust evolution code DustPy. Doing so has allowed us to unveil, for the first time, some insights about how the evolution of dust (dynamics and grain growth processes combined) and MRI-driven accretion are intertwined on million-year timescales, from a more sophisticated modeling of the gas ionization degree. Our key results can be summarized as follows:

1. Dust coagulation and settling lead to a higher gas ioniza-tion degree (the recombination rate onto grains decreases), resulting in stronger MRI-driven turbulence as well as a more compact dead zone. On the other hand, fragmentation has an opposite effect because it replenishes the disk in small dust particles, which are very efficient at sweeping up free electrons and ions from the gas phase (the recombination rate onto grains increases). Since the dust content of the protoplanetary disk decreases over millions of years of evolution due to radial drift, the MRI-driven turbulence overall becomes stronger and the dead zone more compact until the disk dust-gas mixture eventually behaves as a grain-free plasma;

2. The MRI activity evolution (hence the temporal evolution of the MRI-induced α parameter) is controlled by dust evolution and occurs on a timescale of local dust growth, as long as there are enough dust particles in the disk to dominate the recombination process for the ionization chemistry. Once that is no longer the case, the MRI activity evolution is expected to be controlled by gas evolution and occurs on a viscous evolution timescale;

3. Dust evolution alone does not lead to a complete reactivation of the dead zone since the dust eventually has no significant impact on the ionization chemistry when the disk dust-gas mixture behaves as a grain-free plasma. Such result suggests that the dead zone may potentially be able to survive the pro-toplanetary disk evolution over a few million years when the MRI is the main driver for the disk accretion, as long as there are no mechanisms that can efficiently ionize the gas in the regions of the disk where the dead zone sits at and that the gas dispersal in those regions occurs in timescales of a few million years;

4. For typical T-Tauri stars, the dead zone outer edge is expected to be located roughly between 10 au and 50 au during the protoplanetary disk lifetime for our choice of the magnetic field strength and configuration;

5. The MRI activity evolution in protoplanetary disks is expected to be crucially sensitive to the choice made for the minimum grain size of the dust distribution, especially in the early stages of the disk lifetime when the dust has significant feedback on the ionization chemistry. Further studies focusing on constraining such a minimum grain size are thus fundamental.

The evolution on million-year timescales of the MRI activity in protoplanetary disks is a complex problem that significantly depends on the dust properties and how it evolves (this study) as well as the gas and stellar properties. A comprehensive approach to investigate the potential dust trapping power of the dead zone outer edge thus requires a time-dependent framework where MRI calculations are self-consistently combined with gas, dust, and stellar evolution on million-year timescales. Armed with our new framework combining our MRI-driven turbulence model and dust evolution (DustPy), we aim to achieve such a self-consistent model as our next step.

Acknowledgements

This work made extensive use of the Astropy (Astropy Collaboration 2013), Matplotlib (Hunter 2007), Numpy (Harris et al. 2020), Scipy (Virtanen et al. 2020) software packages. T.N.D., P.P. and M.G. acknowledge support provided by the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research. C.-C.Y. is grateful for the support from NASA via the Astrophysics Theory Program (grant number 80NSSC21K0141), NASA via the Emerging Worlds program (grant number 80NSSC20K0347), and NASA via the Theoretical and Computational Astrophysics Networks program (grant number 80NSSC21K0497). S.O. is supported by JSPS KAKENHI Grant Numbers JP18H05438, JP19K03926, JP20H01948, and 20H00182. M.F. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 757957). T.B. and S.M.S. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 714769 and funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grants 361140270 and 325594231. This research was also supported by the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP) which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311.

Appendix A Rayleigh adjustment for a nonuniform radial grid

Appendix A.1 Context

As explained in Delage et al. (2022), a gas surface density profile resulting from steady-state accretion of our MRI-driven turbulence model (Sect. 2.3) often leads to a striking mathematical discontinuity located at the dead zone outer edge. This discontinuity arises from the on/off criteria for active MRI used, inducing a steep change in the local turbulent parameter α at each transition between the dead zone and the MRI active region. As proposed by Yang & Menou (2010), though, such a steep gas transition would not happen because the gas would not be stable and would rearrange itself due to turbulent diffusion on a dynamical timescale (Rayleigh adjustment process). In Model IV, we chose Σ_gas to follow the gas surface density profile obtained for the steady-state MRI-driven accretion solution corresponding to Model II, with MRN-like grain size distribution (a_min = 0.55 µm, a_dist,Max = 1 µm, and p_dist,Exp = − 3. 5). Consequently, we need to apply Rayleigh adjustment in order to smooth the density gradient at the dead zone outer edge, and thus avoid any potential physically inconsistent state while running the dust evolution model employed. After applying Rayleigh adjustment, the new Σ_gas obtained is the gas surface density profile that we refer as “Steady-State 2 (SS2)” in Table 1.

Below we revisit the simple algorithm put forward by Yang & Menou (2010) to implement Rayleigh adjustment. Particularly, we describe a new and more robust method that can be applied to any Rayleigh unstable gas surface densities, for a nonuniform radial grid.

Appendix A.2 Method

A radially one-dimensional diffusion equation in polar coordinates (r, ϕ) reads $$\frac{\partial {\sum }_{\text{gas}}}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(r𝒟\frac{\partial {\sum }_{\text{gas}}}{\partial r}\right),$$(A.1)

where Σ_gas is the gas surface density and Ɗ is the diffusion coefficient. Eq. (A.1) represents a conservation law $$\frac{\partial {\sum }_{\text{gas}}}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left(rℱ\right)=0,$$(A.2)

where the flux is defined by $ℱ\left({\displaystyle \sum {}_{\text{gas}}\left(r,t\right);r,t}\right)\equiv -𝒟\frac{\partial {\displaystyle \sum {}_{\text{gas}}}}{\partial r}$. Integrating Eq. (A.2) over a concentric ring from r = r₁ to r = r₂ and dividing the result by the area of the ring gives $$\frac{\partial }{\partial t}\left(\frac{2}{r_2^2-r_1^2}{\displaystyle {\int }_{r_1}^{r_2}{\sum }_{\text{gas}}r\text{d}r}\right)+\frac{2}{r_2^2-r_1^2}\left[r_2ℱ\left(r_2,t\right)-r_1ℱ\left(r_1,t\right)\right]=0.$$(A.3)

Defining the cell average as $$Q\left(t\right)\equiv \frac{2}{r_2^2-r_1^2}{\displaystyle {\int }_{r_1}^{r_2}{\sum }_{\text{gas}}r\text{d}r}$$(A.4)

and integrating Eq. (A.3) from t = t₁ to t = t₂ gives $$Q\left(t_2\right)=Q\left(t_1\right)-\frac{2\text{Δ}t}{r_2^2-r_1^2}\left[r_{2}F\left(r_2\right)-r_{1}F\left(r_1\right)\right],$$(A.5)

where Δt = t₂ – t₁ and $ℱ\left(r\right)\equiv \frac{1}{\text{Δ}t}{\displaystyle {\int }_{t_1}^{t_2}ℱ\left(r,t\right)\text{d}t}$.

Consider a nonuniform grid with a mapping from the index space ρ to the physical space r, that is, r = r(ρ) = r_ρ. We adopt the convention that the cell edges are indexed by integers and the cell centers by half-integers. Eq. (A.5) then gives the Godunov method for the (n + 1)-th time step $$Q_{j+1/2}^{n+1}=Q_{j+1/2}^n-\frac{2\text{Δ}t}{r_{j+1}^2-r_j^2}\left(r_{j+1}{F}_{j+1}^n-r_j{F}_j^n\right),$$(A.6)

where $$Q_{j+1/2}^n\equiv \frac{2}{r_{j+1}^2-r_j^2}{\displaystyle {\int }_{r_j}^{r_{j+1}}{\sum }_{\text{gas}}\left(r,t_n\right)r\text{d}r},$$(A.7)

and $$F_j^n\equiv \frac{1}{\text{Δ}t}{\displaystyle {\int }_{t_n}^{t_{n+1}}ℱ\left(r_j,t\right)\text{\hspace{0.17em}d}t},$$(A.8)

with Δt = t_n+1 – t_n.

To proceed, one could consider a Riemann problem and use the solution to evaluate Eq. (A.8). Instead, we approximate it by assuming Σ_gas(r,t) is nearly unchanged over t ∈ [t_n, t_n+1] and hence $$F_j^n\approx ℱ\left(r_j,t_n\right)=-{𝒟\frac{\partial {\sum }_{\text{gas}}}{\partial r}|}_{r=r_j}={-𝒟\frac{\text{d}\rho }{\text{d}r}\frac{\partial {\sum }_{\text{gas}}}{\partial \rho }|}_{r=r_j}.$$(A.9)

It can then be discretized using central differences: $$F_j^n\approx -𝒟\left(r_j\right){\rho }^{\prime }\left(r_j\right)\left(Q_{j+1/2}^n-Q_{j-1/2}^n\right),$$(A.10)

where ρ(r) is the inverse function of r(ρ). Therefore, Eq. (A.6) becomes $$Q_{j+1/2}^{n+1}=Q_{j+1/2}^n+\frac{2\text{Δ}t}{r_{j+1}^2-r_j^2}\left[r_{j+1}{\rho }^{\prime }\left(r_{j+1}\right)𝒟\left(r_{j+1}\right)\left(Q_{j+3/2}^n-Q_{j+1/2}^n\right)-r_j{\rho }^{\prime }\left(r_j\right)𝒟\left(r_j\right)\left(Q_{j+1/2}^n-Q_{j-1/2}^n\right)\right].$$(A.11)

If the diffusion coefficient Ɗ is a constant, it can be simplified as $$Q_{j+1/2}^{n+1}=Q_{j+1/2}^n+\frac{2𝒟\text{Δ}t}{r_{j+1}^2-r_j^2}\left[r_{j+1}{\rho }^{\prime }\left(r_{j+1}\right)\left(Q_{j+3/2}^n-Q_{j+1/2}^n\right)-r_j{\rho }^{\prime }\left(r_j\right)\left(Q_{j+1/2}^n-Q_{j-1/2}^n\right)\right].$$(A.12)

Next, we find the stability condition for the time step Δt. Assume that $$Q_{j+1/2}^n\sim e^{ikj},$$(A.13)

where k is any wavenumber in the index space p. Substituting this into Eq. (A.12), a von Neumann stability analysis gives $$\begin{array}{l}{Q}_{j+1/2}^{n+1}\sim e^{ikj}\left\{1+\frac{2𝒟\text{Δ}t}{r_{j+1}^2-r_j^2}\left[r_{j+1}{\rho }^{\prime }\left(r_{j+1}\right)\left(e^{+ik}-1\right)-r_j{\rho }^{\prime }\left(r_j\right)\left(1-e^{-ik}\right)\right]\right\}\hfill \\ =e^{ikj}\left[1+\left(\alpha +\beta \right)\left(e^{+ik}-1\right)-\left(\alpha -\beta \right)\left(1-e^{-ik}\right)\right]\hfill \\ =e^{ikj}\left[1-2\alpha \left(1-\cos k\right)+2i\beta \sin k\right]\equiv A{e}^{ikj},\hfill \end{array}$$(A.14)

where $$\begin{array}{l}\alpha \equiv \frac{𝒟\text{Δ}t\left[r_{j+1}{\rho }^{\prime }\left(r_{j+1}\right)+r_j{\rho }^{\prime }\left(r_j\right)\right]}{r_{j+1}^2-r_j^2},\hfill \\ \beta \equiv \frac{𝒟\text{Δ}t\left[r_{j+1}{\rho }^{\prime }\left(r_{j+1}\right)-r_j{\rho }^{\prime }\left(r_j\right)\right]}{r_{j+1}^2-r_j^2},\hfill \end{array}$$(A.15)

and $$A\equiv 1-2\alpha \left(1-\cos k\right)+2i\beta \sin k.$$(A.16)

We note that α > 0 and β < α if r(ρ) is a strictly monotonically increasing or decreasing function of ρ, which should be always the case. To be stable, |A| = (1 – 2α + 2α cos k)² + 4β² sin² k ≤ 1. The extrema of |A| occur at sin k = 0 or cos k = α(2α – 1)/2(α² –β²). The first leads to (1 – 4α)² ≤ 1, or 0 ≤ α ≤ 1/2. The second results in (α – 2β²)² ≥ 0, which is always true. Therefore, the algorithm of Eq. (A.12) can be optimized by adopting $$𝒟\text{Δ}t\equiv {\left\{2\underset{j}{\max }\left[\frac{r_{j+1}{\rho }^{\prime }\left(r_{j+1}\right)+r_j{\rho }^{\prime }\left(r_j\right)}{r_{j+1}^2-r_j^2}\right]\right\}}^{-1}.$$(A.17)

To conclude, a general outline of the Rayleigh adjustment procedure would be as follows: Firstly, the stability of the gas surface density profile, Σ_gas, should be evaluated using the stability criterion derived in Yang & Menou (2010, either their Eq. (3) for a simple one, or a generalization shown in their Sect. 4). In the present paper, we chose to follow their Eq. (3). Secondly, if the gas surface density profile is unstable (which is the case in the steady-state accretion of our MRI-driven turbulence model), the unstable profile should be used as the initial condition and Eq. (A.12) should be iterated with ƊΔt following Eq. (A.17) until the profile has relaxed to marginal stability. Lastly, the relaxed gas surface density profile should be used to run the 1D dust evolution model (similar to what is done in this paper), or to resume the next time step of the 1D gas and dust evolution model (in this case Rayleigh adjustment needs to be applied at each time step when solving for the gas).

Appendix B Dependence of the steady-state MRI-driven accretion on the power-law dust size distribution

In this appendix we perform a follow-up on the effect of a fixed power-law dust size distribution on the steady-state MRI-driven accretion. Particularly, we investigate the impact of a variation in the dust distribution minimum grain size (Appendix B.1) as well as its exponent (Appendix B.2).

Appendix B.1 Variation in the distribution minimum grain size

We investigate the impact of a variation in a_min, while fixing a_dist,Max= 1 µm and p_dist,Exp = −3.5. This set of simulations corresponds to Model II, and the results are presented in Fig. B.1.

A higher amin implies stronger MRI-driven turbulence overall, a higher gas accretion rate, and a more compact dead zone (Figs. B.1a, B.1e and B.1f). Unlike what we saw in Sect. 4.1, though, we notice that a change in a_min primarily impacts the dead zone outer edge location, which can be almost twice closer from the central star for a_min = 0.55 µm compared to a_min = 0.1 µm. The gas accretion rate and the optimal r.m.s. magnetic field strength appear, indeed, to be weakly dependent on a_min (Figs. B.1d and B.1e).

When a larger grain size is used for the dust distribution minimum size, the gas ionization degree becomes higher closer to the central star (Fig. B.1c), since gas-phase recombination can dominate the recombination process from smaller radial distances. For increasing value of a_min, the total grain surface area, A_tot, decreases and the dust becomes less efficient at sweeping up free electrons and ions from the gas phase at smaller radial distances, hence the MRI being able to operate more easily.

It is important to note that the distribution minimum grain size, a_min, is a free parameter in dust models with a distribution of sizes. From Fig. B.1, it can be seen that a_min is crucial in determining the MRI-driven turbulence. We discuss the implications in Sect. 5.3.

Fig. B.1 Impact on the steady-state MRI-driven accretion when varying the minimum grain size, amin, of the fixed power-law dust size distribution (Model II). Going from solid light-colored to dark-colored lines, amin spans the range from 0.1 µm to 0.55 µm. Panels a–f show the same quantities as in Fig. 2, but for the model parameters M⋆ = 1M⊙, L⋆ = 2 L⊙, Mdisk = 0.05 M⋆, fdg,tot = 10−2, αhydro = 10−4, adist,Max= 1 µm, and pdist,Exp= −3.5. For comparison, the dashed gray and black lines in panel a now display the steady-state quantity $\overline{\alpha }$ obtained assuming the limiting case of a mono-disperse population of dust with size amono = 0.1 µm and amono = 1 µm, respectively.

Appendix B.2 Variation in the power-law exponent

We investigate the impact of a variation in p_dist,Expwhile fixing a_min = 0.1 µm and a_dist,Max = 1 µm. This set of simulations corresponds to Model III, and the results are presented in Fig. B.2.

The power-law exponent (p_dist,Exp) controls the relative proportion of the smaller and larger particles in the grain size distribution: A smaller or more negative value means that the dust distribution is skewed toward the smaller sizes, whereas a larger or more positive value means that it is skewed toward the larger sizes. Looking at Figs. B.2a, B.2e, B.2f, a higher p_dist,Expresults in stronger MRI-driven turbulence overall, a higher gas accretion rate, and a more compact dead zone. In the same spirit as the previous section, we notice that p_dist,Exp mainly changes the dead zone outer edge (although the gas accretion rate and the optimal r.m.s. magnetic field strength vary more with p_dist,Exp than a_min). The dead zone outer edge is, indeed, located at ~63 au for p_dist,Exp= −4.5, whereas it is located at ~27 au for p_dist,Exp= 0.25.

The gas ionization degree becomes higher closer to the central star when the dust size distribution is skewed toward the larger sizes (Fig. B.2c). For increasing value of p_dist,Exp, free electrons and ions are less likely to encounter, per unit volume and on average, the smaller particles of the dust size distribution. Instead, they primarily interact with the larger ones. The total grain surface area (A_tot) thus decreases, allowing for the gas-phase recombination to dominate closer to the central star, hence the MRI being able to operate more easily.

Fig. B.2 Impact on the steady-state MRI-driven accretion when varying the exponent, pdis,Exp, of the fixed power-law dust size distribution (Model III). Going from solid light-colored to dark-colored lines, pdist,Expspans the range from −4.5 to 0.25. Panels a–f show the same quantities as in Fig. 2, but for the model parameters M⋆ = 1M⊙, L⋆ = 2 L⊙, Mdisk = 0.05 M⋆, fdg,tot = 10−2, αhydro = 10−4, amin = 0.1 µm, and adist,Max = 1 µm. For comparison, the dashed gray and black lines in panel a now display the steady-state quantity α obtained assuming the limiting case of a mono-disperse population of dust with size amono = 0.1 µm and amono = 1 µm, respectively.

Since the ionization chemistry is primarily controlled by the smaller sizes of the dust distribution, the MRI activity crucially depends on their relative proportion. Particularly, our results suggest that the MRI activity substantially increases in the regions of the disk where the dust size distribution is cut-off at micron-sized particles rather than submicron-sized particles. This can happen when submicron-sized particles get depleted enough due to effective grain coagulation and less frequent fragmentation.

Appendix C Temporal evolution of the dust surface density distribution for Models V and VI

We define the vertically integrated dust surface density distribution per logarithmic bin of grain size, σ, as (Birnstiel et al. 2010) $$\sigma \left(r,a\right)={\displaystyle {\int }_{-\infty }^{+\infty }{n^{\prime }}_{\text{dust}}\left(r,z,a\right)m\left(a\right)adz.}$$(C.1)

Defining σ(r,a) as in Eq. (C.1) makes it a grid-independent dust density unlike the mass integrated over each numerical bin (Σ_dust(r, a)). This way, all plots of σ(r,a) are meaningful without the knowledge of the size grid that we used. Below we show the temporal evolution of σ(r,a) for Models V and VI. The one for Model IV is shown in Fig. 3.

Fig. C.1 Same as in Fig. 3, except for Model V.

Fig. C.2 Same as in Fig. 3, except for Model VI.

References

All Tables

All Figures

Fig. 1 Flowchart of the MRI-driven disk accretion model presented in Delage et al. (2022). This model captures the essence of the MRI-driven turbulence in a 1+1D framework, accounting for the following: stellar properties (gray symbols), disk gas properties (blue symbols), disk dust properties (red symbols), nonthermal ionization sources (yellow symbols), ionization chemistry modeling the gas ionization degree (green symbols), disk magnetization properties (powder blue symbols), and nonideal MHD calculations (dark pink symbols). The main output of the model is an effective radial profile for the MRI-induced viscosity parameter, ᾱ (Eq. (1)). In this paper, we improve the dust phase modeling with a dust size distribution, either by assuming a fixed power-law distribution with different properties or the outputs from dust evolution obtained with DustPy (see Table 1 and text in Sect. 2.3 for further explanations).

In the text

Fig. 2 Impact on the steady-state MRI-driven accretion when varying the maximum grain size, adist,Max, of the fixed power-law dust size distribution (Model I). Going from solid light-colored to dark-colored lines, adist,Max spans the range from 1 µm to 1 cm. The panels show the steady-state radial profiles of several key quantities (see also Delage et al. 2022, for their definition), for the model parameters M⋆ = 1 M⊙, L⋆ = 2 L⊙, Mdisk = 0.05 M⋆, fdg,tot = 10−2, αhydro = 10−4, amin = 0.1 µm, and pdist,Exp = −3.5. Panel a: pressure-weighted vertically integrated turbulent parameter, $\overline{\alpha }$. Panel b: gas surface density, Σgas. Panel c: midplane total ionization rate, ζ. Panel d: midplane optimal rms magnetic field strength, B. Panel e: (constant) gas accretion rate, ${\dot{M}}_{\text{acc,gas}}$. Panel f : midplane radial dead zone outer edge location, RDZ. For comparison, the dashed gray and black lines in panel a display the steady-state quantity $\overline{\alpha }$ obtained assuming the limiting case of a mono-disperse population of dust with size amono = 0.1 µm and amono = 1 cm (equivalent to the grain-free scenario), respectively.

In the text

Fig. 3 Temporal evolution of the dust surface density distribution per logarithmic bin of grain size σ(r, a) (see Eq. (C.1)), for Model IV. In each panel, the solid white line shows a Stokes number of unity (radial drift reaches its maximal efficiency), the dotted green line shows the drift limit (Eq. (19)), and the dashed cyan line shows the fragmentation limit (Eq. (20)). The horizontal solid dark gray lines show a grain size of 1 mm.

In the text

Fig. 4 Impact of dust evolution on the MRI-driven turbulence, for Model IV. The panels show the temporal evolution of the same quantities as in Fig. 2, except for the gas accretion rate. Also, the gas surface density is now fixed to the displayed input profile in panel b. We emphasize that these quantities do not describe steady-state MRI-driven accretion (unlike Fig. 2), since they are recalculated at each dust evolution snapshot, through partial coupling between the ID radial dust evolution model employed and our MRI-driven turbulence model. On this note, the corresponding temporal evolution of the five dust quantities used to perform such a coupling are shown in Fig. 5. Here we note that, in panel f, the dead zone outer edge coincides at t = 0 yr, t = 100 yr and t = 1000 yr.

In the text

Fig. 5 Temporal evolution of the five dust quantities used to couple the ID radial dust evolution model employed with our MRI-driven turbulence model, for Model IV. These quantities are the representative grain size adust,rep (Eq. (16)) in panel a, the representative grain cross-section σdust,rep (Eq. (17)), the representative grain mass mdust,rep (Eq. (18)), the total dust density ρdust,tot (Eq. (8)) in panel d, and the total dust number density ndust,tot (Eq. (10)) in panel e. Instead of displaying σdust,rep and mdust,rep, we show the equivalent grain size $a_{\text{dust,}\sigma }{\left({\sigma }_{\text{dust,rep}}/\pi \right)}^{\frac{1}{2}}$ corresponding to σdust,rep (panel b), and the equivalent grain size $a_{\text{dust,m}}={\left(3m_{\text{dust,rep}}/4\pi {\rho }_{\text{bulk}}\right)}^{\frac{1}{3}}$ corresponding to mdust,rep (panel c). These equivalent grain sizes better indicate the dominant size when considering the grain cross-section or the grain mass, each of which is important for the MRI calculations (see Sect. 2.3). The panels particularly show the midplane radial profiles. We note that adust,rep and adust,m traces the smallest and largest sizes of the dust distribution, respectively, whereas adust,σ is the relevant size of the dust distribution involved in the ionization chemistry.

In the text

Fig. 6 Same as in Fig. 4, except for Model V. The corresponding temporal evolution of the five dust quantities is shown in Fig. 7. The gas surface density (panel b) is fixed to the input profile corresponding to the classical self-similar solution (Lynden-Bell & Pringle 1974), with a total disk gas mass Mdisk = 0.05 M⋆, and a critical radius Rc = 80 au. As a result, the gas is no longer in a steady-state MRI-driven accretion at t = 0 yr, unlike Model IV. Here we note that, in panel f, the dead zone outer edge coincides at t = 0 yr and t = 100 yr.

In the text

	Fig. 7 Same as in Fig. 5, except for Model V.
In the text

	Fig. 8 Same as in Fig. 6, except for Model VI. The corresponding temporal evolution of the five dust quantities is shown in Fig. 9. Here we note that, in panel f, the dead zone outer edge coincides at t = 0 yr and t = 100 yr.
In the text

	Fig. 9 Same as in Fig. 5, except for Model VI.
In the text

Fig. B.1 Impact on the steady-state MRI-driven accretion when varying the minimum grain size, amin, of the fixed power-law dust size distribution (Model II). Going from solid light-colored to dark-colored lines, amin spans the range from 0.1 µm to 0.55 µm. Panels a–f show the same quantities as in Fig. 2, but for the model parameters M⋆ = 1M⊙, L⋆ = 2 L⊙, Mdisk = 0.05 M⋆, fdg,tot = 10−2, αhydro = 10−4, adist,Max= 1 µm, and pdist,Exp= −3.5. For comparison, the dashed gray and black lines in panel a now display the steady-state quantity $\overline{\alpha }$ obtained assuming the limiting case of a mono-disperse population of dust with size amono = 0.1 µm and amono = 1 µm, respectively.

In the text

Fig. B.2 Impact on the steady-state MRI-driven accretion when varying the exponent, pdis,Exp, of the fixed power-law dust size distribution (Model III). Going from solid light-colored to dark-colored lines, pdist,Expspans the range from −4.5 to 0.25. Panels a–f show the same quantities as in Fig. 2, but for the model parameters M⋆ = 1M⊙, L⋆ = 2 L⊙, Mdisk = 0.05 M⋆, fdg,tot = 10−2, αhydro = 10−4, amin = 0.1 µm, and adist,Max = 1 µm. For comparison, the dashed gray and black lines in panel a now display the steady-state quantity α obtained assuming the limiting case of a mono-disperse population of dust with size amono = 0.1 µm and amono = 1 µm, respectively.

In the text

	Fig. C.1 Same as in Fig. 3, except for Model V.
In the text

	Fig. C.2 Same as in Fig. 3, except for Model VI.
In the text