Impact of dust size distribution including large dust grains on magnetic resistivity: an analytical approach

This paper investigates the impact of dust size distribution on magnetic resistivity. In particular, we focus on its impact when the maximum dust size significantly increases from sub-micron. The first half of the paper describes our calculation method for magnetic resistivity based on the model of \citet{1987ApJ...320..803D} and shows that the method reproduces the results of a more realistic chemical reaction network calculations reasonably well. Then, we describe the results of the resistivity calculations for dust distributions with large maximum dust grains. Our results show that resistivity tends to decrease with dust growth, which is particularly true when the dust size power exponent $q$ is $q=2.5$. On the other hand, the decrease is less pronounced when the dust size power exponent $q$ is $q=3.5$, i.e., when the small dust is also responsible for the dust cross-section. Our results suggest that detailed dust coagulation and fragmentation processes play a vital role in the magnetic resistivities in protostar formation.

The degree of impact of the non-ideal effect depends on the magnetic resistivities, which are determined by the amount of charged particles, and hence ionization chemistry. In the ionization chemistry, dust grains absorb the ions and electrons and affect their abundance. The adsorption efficiency depends on the total cross section of dust grains. Furthermore, a large population of charged small dust grains ( 10 nm) can contribute to the conductivities (Zhao et al. 2016).
The dust size distribution in the previous studies is often assumed to be that of the interstellar medium (ISM) such as MRN size distribution (Mathis et al. 1977) or submicron sized dust grains. However, the distribution may change through dust coagulation during protostar formation, particularly in the protoplanetary disks. By assuming that the relative velocity among the dust is determined by turbulence, the growth timescale t growth of dust grains in the disk is calculated as (Ormel & Cuzzi 2007), t growth = 1.6 × 10 3 α where ng, α, ρmat, a d , f , cs, M * denotes the gas number density, viscous α value, material density and size of the dust grains, dust-to-gas mass ratio, sound velocity, and mass of the central protostar, respectively. fX means fX = ( f X ). Here we assume turbulent velocity and timescale of largest eddy to be ∆vL = √ αcs and tL = Ω −1 , respectively.
Thus, the dust growth timescale is about 100 times smaller than the age of Class 0/I young stellar objects (YSOs), and dust growth may proceed even in the early evolution of circumstellar disks. Actually, recent 3D simulation by Tsukamoto et al. (2021b) shows that dust growth in the disk (and reflux of large dust to the envelope). The question we address in this paper is how the change of dust size distribution caused by the dust growth affects resistivities. Care should be taken in applying the usual chemicalreaction-network calculations to calculate the resistivities with large dust grains because the mean grain charge Z of 1 µm is typically (Draine & Sutin 1987) Z ∼ −20a d,10 µm T10K, when abundances of ions and electron is much larger than that of dust grains. Since dust grains with different charges need to be treated as different chemical species in chemical reaction network calculations, a vast number of charged dust species should be considered, which is computationally demanding. Therefore, analytical models of ionization chemistry such as Draine & Sutin (1987); Okuzumi (2009);Tsukamoto et al. (2021a) are more suitable for the calculation of resistivities with large dust grains. Thus, we adopt this approach in this paper. This paper is organized as follows. In §2, we describe our analytical model which is based on Draine & Sutin (1987). In §3.1, we validate the analytical model by comparing it with chemical reaction calculations. In §3.2, we investigate the magnetic resistivities with large dust grains. Finally, the results are summarized and discussed in §4.

Equilibrium of ionization recombination reaction
We start from equations for chemical equilibrium in the gas phase.
, are average ion velocity, gas-phase recombination rate coefficient, ionization rate. Here ni = k n (k) i and ng = k n (k) g are the total number density of ion and neutral respectively.J i(e) (I, Z) is the effective cross sections normalized by σ d (I) = π(a d (I)) 2 between dust grains and ion (electron). a d is the dust radius. (A) and A denotes the average over dust size (I) and charge (Z), respectively. Draine & Sutin (1987) derive the approximation formula for the effective cross section of charged particles as, where ν = ν(Z, q i(e) ) = Ze/q i(e) where q i(e) are charge of ion or electron. τ = τ (I) = a d (I)kBT /e 2 is the normalized temperature and e is the elementary charge. These formula are correct within few % for τ > 10 −3 (thus, the particle size of a d > 1.67nm(T /(10K)) −1 which is enough for our purpose).
UsingJ(τ, ν), we obtain, for singly charged ions and (5 for electrons. The charge neutrality condition is given as The governing equations for dust charging is the detailed balance equation for dust grains which is given as where we define The final governing equation is number density conservation for each I n d (I) = Z n d (I, Z).
We assume that ng, ζ, s i(e) , a d (I), β, and n d (I) are known. Our purpose is to obtain ni, ne, and n d (I, Z). In the following three subsections, we describe the procedure for computing these quantities.

High τ case
The equations in the previous subsection hold as long as the dust charge remains at Z = ±1, 0. However, as τ becomes large (i.e., the temperature of the gas increases or the size of the dust increases), the typical dust charge becomes Z ∼ −kBT a d /e 2 = −τ ≪ −1 and very small. Thus, the strategy in the previous section of solving the detailed balancing equations in sequence is not useful for Z ≪ 1 (or τ ≫ 1) because we have to consider a large number of detailed balance equations. On the other hand, for large τ , the dust charge distribution can be treated as a continuous distribution, known to become Gaussian distribution. Then Z , J i , J e can be obtained analytically (Draine & Sutin 1987;Okuzumi 2009). For τ ≫ 1, we can approximateJi(I, Z) andJe(I, Z) of equation (4) and (5) assuming Z < 0, for singly charged ion and for electron. For τ ≫ 1, the solution of detailed balance equation (i.e., equation (7)) is given as (Draine & Sutin 1987;Okuzumi 2009), where The dimensionless parameter ψ is the solution of the equation of siniuiJi(I, Z ) = seneueJe(I, Z ), and hence, ψ is a function of ǫ. Equations (24) to (27) are derived from the detailed balance equation (7) with the assumption of n d (I, Z + 1) ∼ n d (I, Z) + ∂n d (I, Z)/∂Z and Je(I, Z + 1) ∼Je(I, Z ) (for the detail, see Okuzumi 2009). Using these equations, we can calculate Z , J i(τ, Z) , and J e(τ, Z) for high temperature case as,

Number density of ions and electrons
In the previous two sections we have obtained Z , J i , and J i for each I in the high and low τ limits.
Here we explicitly write the variables of these quantities. By summing these up for I, we can calculate Z , J i , and J i as a function of ǫ.
Then ni and ne are obtained from equation (3) as a function of ǫ, 1 + 4βζng 1 + 4βζng The charge neutrality condition becomes, Equations (35) is a nonlinear algebraic equation for ǫ, and we solve this equation with Newton-Raphson method.

Conductivity and magnetic resistivity
Using the ni, ne, and n d (I, Z) obtained in the previous sections, the conductivity is calculated as follows (Wardle 2007), where σO,H,P are the Ohmic, Hall, and Pedersen conductivities, respectively, of the charged species.
is the product of the cyclotron frequency and the collision frequency with the neutral gas. The subscript s denotes the charged species. Here ns and qs are the number density and charge of the species s. B and c are the magnetic field strength and speed of light, respectively. γs = σv s/(ms + mg) and σv s is the collisional momentum transfer rate between species s and the neutrals. mg is the mean mass of the gas. The momentum transfer rate between neutral and charged species was calculated using the equations described in Pinto & Galli (2008). The conductivities of dust grains are separately calculated from low temperature and high temperature dust size distribution (equations (14) and (23)) and then summed up. This treatment is necessary to correctly calculate the Pedersen conductivity to which both the positively and negatively charged dust grains positively contribute. This method adds the conductivity of the dust in duplicate, but we confirmed that this does not cause an error because the contribution of dust at higher temperature is small (see §3.1).
The Ohmic, Hall, and ambipolar resistivities are calcu-lated as

Chemical reaction network calculation
We perform chemical reaction network calculations to compare with the analytical model above.
In the chemical reaction network calculations, we consider ion species H3O + , OH + , H2O + and neutral species H, H2, He, CO, O2, Mg, O, C, HCO, H2O, OH, N, Fe. We also consider neutral and singly charged dust grains, G 0 , G − , G + . We consider cosmic-ray ionization, gas-phase and dust-surface recombination, and ion-neutral reactions. We also considered the indirect ionization by high-energy photons emitted by direct cosmic-ray ionization (described as CRPHOT in the UMIST database). The initial abundance and reaction rates are taken from the UMIST2012 database (McElroy et al. 2013). We neglect grain-grain collisional neutralization so that the chemical network calculations are consistent with the analytical model. The chemical reaction network is solved using the CVODE package (Hindmarsh et al. 2005). We calculate the conductivities using the abundances of charged species in the equilibrium state.

RESULTS
In this section, we compare the analytical calculation with chemical network calculations and previous studies to justify the analytical calculations in §3.1. In §3.2, we investigate the impact of the dust size distribution with large dust grains on magnetic resistivity using the analytical calculation.

Validation of the analytic model
In this subsection, we compare the analytic model with the chemical reaction network calculation. Here, we assume that the temperature is T = 10(1 + γT (ng/nc) (γ−1) ) K, where γ = 7/5 and nc = 2.6 × 10 10 cm −3 , the magnetic field is 0.2n 1/2 g, cm −3 µG (i.e., assuming flux freezing; see e.g., Nakano et al. 2002), the dust internal density is ρmat = 2 g cm −3 , the dust-to-gas mass ratio is f = 0.01, and the cosmic ray ionization rate of ξCR = 10 −17 s −1 except for the calculations presented in figure A1. We assume that the dominant ion is HCO + and adopt its recombination rate of β = 2.4 × 10 −7 (T /300) −0.69 and its mean molecular weight of µI = 29 for ion in the analytic model. ion e -G + G Gion analytic eanalytic Ganalytic G analytic G + analytic Figure 1. Fractional abundance of ion, electron and positively, negatively, and neutral dust grains from the analytic calculation (dotted lines). The fractional abundance of total ions, electrons, and charged and neutral dusts from chemical reaction calculations are also plotted with solid lines. The orange and blue-shaded regions represent the regions within a factor of three of the electron and ion abundances from chemical reaction network, respectively. The dust size is assumed to be constant of a d = 0.1 µm in the top panel. The dust size distribution is assumed to be MRN size distribution in the bottom panel.
MRN dust size distribution (Mathis et al. 1977), in which the dust size is assumed to be where q = 3.5, amin = 5 nm, and amax = 250 is a constant for normalization. µg = 2.34, µH = 1.4, and ρ d is the dust mass density. dn d da d is the number of dust grains whose sizes are between a d and a d + da d per hydrogen nucleus.
The top panel of figure 1 shows the fractional abundance with a d = 0.1 µm. The ion abundance of the analytic model is in good agreement with the chemical reaction network. On the other hand, our model tends to overestimate the electron abundance in high density region. The figure shows that ion and electron abundance difference between analytic calculation and chemical reaction network is within a factor of three. The negatively charged dust grains G − and neutral dust grans G 0 are dominant in low ( 10 10 cm −3 ) and high-density regions ( 10 10 cm −3 ), respectively, and are in a good agreement between the analytic model and chemical reaction network. Although the abundance of G + and G 0 slightly different between the two calculations around 10 8 cm −3 , this does not cause errors for conductivities.
Bottom panel of figure 1 shows the fractional abundance with MRN size distribution. Even with the size distribution, the similar trend is seen as in the case with a d = 0.1 µm, and the analytic model and chemical reaction network are in a good agreement. Figure 2 shows the conductivities from the chemical network calculation and the analytic calculation of this work with a d = 0.1 µm. The figure shows that electrons dominate Ohmic conductivity, so there is about a factor of three discrepancy over the entire region due to differences in the electron abundance. On the other hand, for the Hall and Pedersen conductivities, the deviation is much smaller than for the Ohmic conductivity because ions dominate them. As a result, the resulting error for Ohmic resistivity is also about a factor of three, and for Hall and ambipolar resistivity, the error is smaller than that apart from ηH of the very low and high density region as shown in figure 4. In this figure, we also plot the conductivities obtained by the method of Okuzumi (2009) which is valid in τ ≫ 1. Since τ < 1 almost everywhere in this plot, we can see that the difference becomes larger when the dust charge affects the ion/electron abundance. Figure 3 shows the conductivities from the chemical network calculation and the analytic model with MRN size distribution. In the low-density region of ng 10 11 cm −3 , Ohmic conductivity is determined by that of electrons as with a d = 0.1 µm, so there is up to about a factor of three discrepancy due to differences in the abundance of electrons. In the high-density region of ng 10 11 cm −3 , Ohmic conductivity is determined by that of dust grains, and the discrepancy becomes much smaller. For Hall conductivity, the deviation is sufficiently small apart from the low-density region of ng 10 5 which is due to simplified ion treatment and high-density region of ng 10 14 which is due to the difference of electron abundance. The deviation is sufficiently small for Pedersen conductivity because they are mainly determined by small dust grains. As a result, the resulting error for resistivities is also within a factor of three for MRN size distribution apart from the low density region of ng < 10 −6 cm −3 where non-ideal MHD effect is not important as shown in figure 5.

Average dust charge from low to high temperature
The results in the previous subsection only confirm that the analytical model is consistent with chemical reaction network calculations for small τ (i.e., small dust and low temperature). In this subsection, we will further check that our model in the high-temperature regime is consistent with Draine & Sutin (1987). Figure 6 shows the mean dust grain charge Z as a function of normalized temperature τ . The cyan dashed line shows the analytic formula of Draine & Sutin (1987) which is given as where τ0 = 8/(πµ)(me/mp) and µ = (sene/ni) 2 (mi/mp).
The figure shows Z obeys Z ∝ τ and our analytic model well reproduces the analytic formula of Draine & Sutin (1987). This means that our model can correctly calculate the ionization state of dust grains not only for small τ but also large τ .

Impact of dust size on magnetic resistivity
In this subsection, we investigate the impact of dust growth on the magnetic resistivities. We assume power law dust size distribution (equation (43)) and vary amin, amax, and q as parameters.
The important quantity is total cross-section Stot of the  dust grains, which determines the absorption efficiency of ions and electrons. If 4 > q > 3 and amax ≫ amin, Stot ∝ amin −q+3 /amax −q+4 and both small and large dust grains affects the total cross section. In this case, it is expected that the impact of dust growth would be less significantly on the absorption efficiency of ions and electrons, and hence resistivities. On the other hand, if q < 3 and amax ≫ amin, Stot ∝ amax −1 and only large dust grains determines the total surface area. Thus, it is expected that dust growth would significantly change the resistivities for q < 3. It is pointed out that q ∼ 2.5 when the coagulation process dominates while q ∼ 3.5 when the disruption process dominates (Miyake & Nakagawa 1993). Thus, we investigate the resistivities with q = 2.5 and q = 3.5.
On the other hand, a large population of charged small grains ( 10 nm) can be responsible for the conductivities (Zhao et al. 2016, see also figure 3), and decreases the resistivities at envelope and disk (Zhao et al. 2018a;Koga et al. 2019;Marchand et al. 2020;Tsukamoto et al. 2020). Thus, whether such small grains exist or not would also be important. In this subsection, we investigate the resistivities with amin = 5 nm and amin = 0.1 µm. Figure 7 shows the resistivities with amin = 5 nm and q = 2.5 for different amax. The figure shows that ηO tends to decrease with increasing amax aside from amax = 0.25 µm to amax = 2.5 µm in high density region. The increase of ηO and ηA from amax = 0.25 µm to amax = 2.5 µm in the high-density region is due to the decrease of dust conductivity. Then, it can be seen that resistivities converge to the power law of the form of ηO,H,A ∝ n p g with respec- In this plot, we set a d = 1 µm and ng = 10 4 cm −3 . The cyan dashed line shows the fitting formula of Draine & Sutin (1987) (equation (44)). The red dotted line shows the Z = ψτ , i.e., the mean dust charge of the high-τ case.
tive constant power exponent p as the dust size increases to amax = 2.5 × 10 3 µm. They are power laws determined by the chemical equilibrium of cosmic-ray ionization and gas-phase recombination. ηH becomes positive (shown with dashed lines) in amax > 250 µm almost entire density region. This is because the relative velocity between ions and electrons determines the Hall current. Figure 8 shows the resistivities with amin = 5 nm and q = 3.5 for different amax. The figure shows that ηO tends to decrease with increasing amax in the low-density region, which is due to the decrease of the total cross-section. On the other hand, ηO increases in the high-density region, which is due to the decrease of dust conductivity. ηH and ηA also increase in the high-density region, which is also due to the decrease of dust conductivity there. The difference between q = 2.5 (dotted lines) and q = 3.5 (solid lines) is striking. If we compare ηO and ηA at 10 15 cm −3 between figure 7 and 8 , the difference is more than 10 3 times greater, which may significantly affect the disk evolution. Figure 9 shows the resistivities with amin = 0.1 µm and q = 2.5 for different amax. By removing small dust grains, resistivity behavior is simplified because the dust grains themselves are no longer responsible for conductivity and serve only as adsorber of ions and electrons. The figure shows that ηO tends to decrease with increasing amax in the entire density region. It can be seen that resistivities converge to the single power law as the dust size increases, which is the same as the figure 7. Again ηH becomes positive in amax > 250 µm almost entire density region, which is also consistent with the figure 7. Figure 10 shows the resistivities with amin = 0.1 µm and q = 3.5 for different amax. Similar to the result of figure 9, ηO tends to decrease with increasing amax. However, the decrease is less pronounced because of the contribution of the small dust to the cross-section. ηA increases in the low-density region of ng 10 10 cm −3 and converges to the single power law as the dust size increases. On the other hand, it decreases with increasing amax in the high density region of ng 10 10 cm −3 . However, again the decrease is less pronounced. ηH tends to decrease with increasing amax.

DISCUSSION
In this paper, we investigate the impact of dust size distribution with large dust grains on the magnetic resistivities using the analytic method based on Draine & Sutin (1987). Our test results show that the analytic model can correctly calculate the ionization state from small τ (low temperature or small dust) to large τ (high temperature or large dust). Therefore, the method is applicable to a dust size distribution that simultaneously contains small dust with τ < 1 and large dust with τ > 1, and can be used over a broader class of dust size distribution than previous study which uses the Gaussian charge distribution such as Okuzumi (2009); Tsukamoto et al. (2021a).
The calculation results with large dust grains show that the resistivity tends to decrease with dust growth. This is particularly true when the dust size power exponent q is q = 2.5 (i.e., in the case the coagulation process dominates in the dust size evolution, and only large dust grains are responsible for the dust cross-section). On the other hand, the decrease is less pronounced when the dust size power exponent q is q = 3.5, (i.e., in the case the disruption process dominates in the dust size evolution, and the small dust grains are also responsible for the dust cross-section). Our results suggest that detailed dust coagulation and fragmentation processes play a crucial role to investigate the impact of non-ideal effects, in particular in the high density region of 10 10 cm −3 .
Recently, Marchand et al. (2021) proposed a similar method that also can be used to calculate magnetic resistivity analytically. Our numerical tests showed that our method seems to be more robust and applicable over a wide parameter range (e.g., when the dust grains are highly depleted). When we implemented and tested their algorithm, we found that it did not converge in the limit where the total dust charge goes zero. This is because their method uses ψ as the basic variable and use equation (A.3) and (A.4) of Marchand et al. (2021) to determine ǫ and ni. However their equation (A.4) becomes singular when ǫ = 1 and Z = 0 i.e., gas phase recombination determines the ionization state. That would be a reason why the solution does not converge.
Our model does not include charge neutralization due to grain-grain collisions which many previous studies have included (e.g., Umebayashi & Nakano 1990;Tsukamoto et al. 2015b;Marchand et al. 2016). One might find this to be a flaw in our model. However, we would argue here that inclusion of charge neutralization by grain-grain collisions is debatable and may not necessarily describe a realistic dust charge state in dense region. This is because, when (sub-)micron-sized small dust particles collide, the dust particles tends to coalesce and grow rather than bounce because the collisional velocity of small dust grains tends to be much smaller than their bouncing velocity ] a max =2.5 10 -1 µm a max =2.5 10 0 µm a max =2.5 10 1 µm a max =2.5 10 2 µm a max =2.5 10 3 µm Figure 7. η O , η H , and η A with a min = 5 nm and q = 2.5. The black, blue, red, cyan, orange lines show the results of amax = 2.5 × 10 −1 µm, 2.5 × 10 0 µm, 2.5 × 10 1 µm, 2.5 × 10 2 µm, 2.5 × 10 3 µm, respectively.
where m th = 1.1 × 10 −15 g (Weidling et al. 2012). On the other hand, the relative velocity of the dust (assuming Brownian motion) is given as By solving the inequality of ∆v stick > ∆vBrown, we can conclude that the silicate dust grain with a d > 1.2(T /100K) 3/4 nm tends to stick rather than bounce, which is much smaller than the minimum size of MRN size distribution (amin = 5 nm). Note also icy or porous dust grains are even more sticky (Wada et al. 2009). Hence, grain-grain collisions lead to dust growth and a change in the dust size distribution rather than bounce.
Therefore, the approximation that ignores charge neutralization due to grain-grain collisions is not necessarily a flaw in our model, but rather a difference in the approximation of how we view the dust collision process. Note also that neglecting grain-grain neutralization does not cause the change on the resistivities when we consider the large dust (e.g., 1 µm), which is the main subject of current and our subsequent studies. This is because the charged dust grains determines the resistivities only when there are sufficient sub-micron dust grains ( 100 nm).
On the other hand, the comparison with Guillet et al. (2020) is difficult because their dust size distribution changes as density increases and they also include non-thermal dust drift due to ambipolar diffusion. However, the following consistent trends can be observed. Our figure 9 shows that smaller power exponent q (meaning that small dust aggregates are less abundant) causes significant decreases of ηO and increase of ηA around intermediate density (ng ∼ 10 9 cm −3 ). On the other hand, figure 9 of Guillet et al. (2020) shows that larger VAD (causing the removal of small dust) results in the decrease of ηO and the increase of ηA around the intermediate density. These points are consistent.
The method described in this paper is less computationally expensive than conventional chemical reaction network calculations and easily converges to the solutions because it only requires to perform one-dimensional Newton-Raphson method twice (in determining ψ and ǫ). Our method typically requires only a few (typically 1-4 times) iterations for each Newton-Raphson calculation to obtain a result. Therefore, it can be easily used in 3D simulations with negligible computational costs. Upon request, we will provide a sample implementation of our method to the readers.
We plan to use the analytic model in this paper in 3D MHD simulations of disk formation and evolution which incorporates dust growth.

ACKNOWLEDGMENTS
This work is supported by JSPS KAKENHI grant number 18H05437, 18K13581, 18K03703.
APPENDIX A: COMPARISON WITH Marchand et al. (2021) For comparison with the previous study by Marchand et al. (2021), we plot the number density of ions and electrons and − Z calculated from our analytic model in figure A1. In this figure, we assume the parameters of Marchand et al. (2021). This figure can be directly compared to figure 2 of Marchand et al. (2021). The results in this figure are in good agreement with their results. Quantitatively, we confirmed that ni/ne converges to Θ ≡ se(mi/me) 1/2 = 107 in ng 10 12 cm −3 , which is also in good agreement with their results.