Evolving grain-size distributions embedded in gas flows

ABSTRACT

We present a numerical approach for accurately evolving a dust grain-size distribution undergoing number-conserving (such as sputtering) and/or mass-conserving (such as shattering) processes. As typically observed interstellar dust distributions follow a power law, our method adopts a power-law discretization and uses both the grain mass and number densities in each bin to determine the power-law parameters. This power-law method is complementary to piecewise-constant and linear methods in the literature. We find that the power-law method surpasses the other two approaches, especially for small bin numbers. In the sputtering tests, the relative error in the total grain mass remains below 0.01 per cent independent of the number of bins N, while the other methods only achieve this for N > 50 or higher. Likewise, the shattering test shows that the method also produces small relative errors in the total grain numbers while conserving mass. Not only does the power-law method conserve the global distribution properties, it also preserves the inter-bin characteristics so that the shape of the distribution is recovered to a high degree. This does not always happen for the constant and linear methods, especially not for small bin numbers. Implementing the power-law method in a hydrodynamical code thus minimizes the numerical cost while maintaining high accuracy. The method is not limited to dust grain distributions, but can also be applied to the evolution of any distribution function, such as a cosmic ray distribution affected by synchrotron radiation or inverse-Compton scattering.

hydrodynamics, plasmas, methods: numerical, dust

1 INTRODUCTION

Within the interstellar medium (ISM), dust grains are an important ingredient as they lock up substantial fractions of the heavy elements (Draine et al. 2007), produce the dominant contribution to the opacity for radiation upward of the Lyman limit (Draine & Lee 1984), contribute a significant part of gas heating through photoelectric heating (Bakes & Tielens 1994), and provide the surface on to which chemical elements can accrete and react (Garrod & Herbst 2006). They also make up a significant fraction of the mass of the ISM with a canonical value of about 1 per cent. This value comes from observational constraints by fitting extinction, scattering or polarization of background stellar radiation and IR dust emission (e.g. Knapp & Kerr 1974; Jura 1979). These restrictions further imply that the dust grains follow a power-law size distribution given by

$$\begin{eqnarray*} \frac{\mathrm{ d}n(a)}{\mathrm{ d}a} \propto a^{-3.5}, \end{eqnarray*}$$

(1)

where a is the grain radius and |$\frac{\mathrm{ d} n(a)}{\mathrm{ d}a} \mathrm{ d}a$| is the number density of grains with radii in the range [a, a + da] (Mathis, Rumpl & Nordsieck 1977; Weingartner & Draine 2001).

Although the typical grain-size distribution follows the Mathis et al. (1977) (MRN) distribution, local variations are expected as the grain-size distribution is set by balancing dust production, growth, and destruction processes. It is thought that dust grains are produced in the upper atmospheres of asymptotic giant branch (AGB) stars (Maercker et al. 2018) as well as formed (Nozawa et al. 2003) and destroyed (Kirchschlager et al. 2019) in supernova remnants. In the dense, quiescent regions of molecular clouds, grains primarily grow in size due to coagulation and mantle accretion (Jones & Williams 1985; Liffamn & Clayton 1989; Ossenkopf 1993; Inoue 2003; Ormel et al. 2009; Asano et al. 2013; Ysard et al. 2016; Jones et al. 2017). In contrast, in star-forming regions, shocks mainly accommodate the destruction of dust grain sizes due to sputtering, shattering, and vaporization (Tielens et al. 1994; Jones, Tielens & Hollenbach 1996; Flower & Pineau des Forêts 2003; Guillet, Pineau des Forêts & Jones 2007; Guillet, Jones & Pineau des Forêts 2009; Hirashita & Yan 2009; Guillet, Pineau des Forêts & Jones 2011; Anderl et al. 2013; Van Loo et al. 2013). In protoplanetary discs, both growth and fragmentation of dust grains take place (Brauer, Dullemond & Henning 2008; Birnstiel, Dullemond & Brauer 2010; Birnstiel et al. 2018; Dullemond et al. 2018; Homma & Nakamoto 2018; Tamfal et al. 2018).

As these grain processes affect the overall grain distribution, they can have significant effects on the dynamics of e.g. the ISM and protoplanetary discs. For example, in the outflows of young stellar objects (YSOs), dust grains are important charge and current carriers and therefore determine the structure of C-type shocks (van Loo et al. 2009), while, in protoplanetary discs, different dust grain sizes influence the growth and structure of the magnetorotational instability (Salmeron & Wardle 2008). Furthermore, grain-processing leads to observational signatures. In typical ISM conditions, silicon is adsorbed on to dust grains. However, gas phase SiO is detected in the clumpy structure of YSO outflows due to the shock-induced sputtering releasing silicon into the gas phase (e.g. Martin-Pintado, Bachiller & Fuente 1992; Mikami et al. 1992). Thus, it is essential to accurately follow the evolution of the dust grain distribution to model both the dynamics and emission signatures.

Previous studies of the dust processing have used different methods. In the simplest approach only a few dust species with specified radii, i.e. typically, one or two grain species representing small and/or large dust grains, are evolved (e.g. Draine, Roberge & Dalgarno 1983; Van Loo et al. 2013; Hirashita 2015). While this is appropriate to capture number-conserving grain processes such as sputtering and mantle accretion, it does not adequately model mass-conserving processes such as shattering¹ and coagulation. A more rigorous approach is to follow the dust grain distribution using a discrete grain-size distribution. Here the distribution is updated, by considering number conservation or mass conservation, to redistribute the grains across the bins in a way appropriate to capture the modelled grain process (e.g. Mizuno, Markiewicz & Voelk 1988; Liffamn & Clayton 1989; Jones et al. 1994; Jones et al. 1996; McKinnon et al. 2018). As most of these studies focus on coagulation and shattering, i.e. the mass-conserving grain processes, these models only follow the total dust mass and not the total grain number. In order to also model sputtering and mantle accretion, McKinnon et al. (2018, hereafter McK18) modified the discrete distribution approach by using a piecewise-linear discretization in each grain-size bin to conserve both mass and particle numbers. They show that this technique is second-order accurate in the number of size bins. However, to achieve an accuracy of the order of 1 per cent in both mass and number conservation, it is necessary to use about 50–100 bins. Thus, including such a scheme into a numerical hydrodynamics (HD) or magnetohydrodynamics (MHD) code significantly increases the computational cost.

In this paper, we propose a different approach to the discretization within a size bin in order to minimize the numerical cost while maintaining high accuracy. As the typical observed dust distribution follows a power law, i.e. the MRN distribution, we will adopt a power-law discretization. Section 2 outlines the method describing the determination of the power-law coefficient and index, and the redistribution of the dust grains for both number-conserving and mass-conserving processes. In Section 3, the method is applied to a number of test problems and the results are discussed. Then Section 4 describes the modifications needed to the HD and MHD equations in order that these reflect the power-law distribution of the dust grains, and, finally, our conclusions are given in Section 5.

2 NUMERICAL METHOD FOR GRAIN DISTRIBUTION EVOLUTION

Here we will describe the numerical methods to evolve a grain power-law distribution. First we consider the construction of a piecewise power-law distribution and the formulation of the power-law coefficients and indices. As we compare the power-law method to piecewise-linear and piecewise-constant discretizations in Section 3, we also present these formulations and explain how they relate to one another. Finally, the routines to redistribute mass and number density of grains across the distribution bins are discussed for grain processes that conserve the grain numbers or total grain mass.

2.1 Discrete power-law distribution

Although dust grains are generally irregularly shaped (e.g. Draine 2003), for the purpose of this paper we assume they are spherical. This significantly simplifies the treatment of the dust grains as the grain distribution depends only on the grain radius a. Furthermore, we assume that the range of grain radii is limited to the range [a_min, a_max]. This range is then divided up logarithmically with a spacing determined by

$$\begin{eqnarray*} \Delta a = \frac{\log (a_{\rm max}/a_{\rm min})}{N}, \end{eqnarray*}$$

(2)

where N is the number of bins. This means that the edges of bin i are effectively

$$\begin{eqnarray*} a_{i} &=& a_{\rm min}\mathrm{e}^{i \Delta a} \nonumber\\ a_{i+1} &=& a_{\rm min}\mathrm{e}^{(i+1)\Delta a}, \end{eqnarray*}$$

(3)

where i = 0, 1, ..., N − 1. Now we assume that the differential grain-size number density distribution in bin i has a power-law shape, i.e.

$$\begin{eqnarray*} \left.\frac{\partial {n(a, t)}}{\partial {a}}\right|_i = A_i a^{-\alpha _{i}}, \end{eqnarray*}$$

(4)

where |$\frac{\partial {n(a,t)}}{\partial {a}} \mathrm{ d}a$| is the number density of grains in a size interval [a, a + da] at time t. Note that, as grain processes change the distribution function, the power-law coefficient A_i and index α_i are implicitly time-dependent.

To determine the power-law coefficient and index, we use the grains’ bin-averaged number density, n(t), and mass density, ρ(t), which are followed according to the redistribution routine in Section 2.2. The bin-averaged number density in bin i at time t is given by the integral

$$\begin{eqnarray*} n_i(t) = \int \nolimits _{a_i}^{a_{i+1}} \frac{\partial {n(a, t)}}{\partial {a}}\Big|_{i}\, \mathrm{ d}a = \left\lbrace \begin{array}{l{@}{\quad}l} \frac{A_i}{1 - \alpha _i}\left(a^{1-\alpha _i}_{i+1} - a^{1-\alpha _i}_{i}\right) & \alpha _i \ne 1 \\ A_{i} \log (a_{i+1}/a_i) & \alpha _i = 1 \end{array}\right. . \nonumber\\ \end{eqnarray*}$$

(5)

It is actually convenient to use equation (3) to rewrite this expression as

$$\begin{eqnarray*} n_i(t) = A_i a^{1 - \alpha _i}_{i+1/2} \Delta a\ \mathcal {F}\left((1 - \alpha _i) \frac{\Delta a}{2}\right), \end{eqnarray*}$$

(6)

where |$\mathcal {F}(x) = \sinh (x)/x$| and applies to all values of α_i, and a_i+1/2 = a_mine^{(i + 1/2)Δa}. However, this does not uniquely determine A_i and α_i. Therefore we need a second expression provided by the grains’ bin-averaged mass density

$$\begin{eqnarray*} \rho _i(t) = \int \nolimits _{a_i}^{a_{i+1}}{m(a)\left.\frac{\partial {n(a, t)}}{\partial {a}}\right|_i \mathrm{ d}a}, \end{eqnarray*}$$

(7)

where

$$\begin{eqnarray*} m(a) = \frac{4\pi \rho _g}{3} a^3 \end{eqnarray*}$$

(8)

is the grain mass with ρ_g the density of the grain material. Similarly to equation (6) we obtain

$$\begin{eqnarray*} \rho _i(t) = \frac{4\pi \rho _\mathrm{ g}}{3} A_i a^{4 - \alpha _i}_{i+1/2} \Delta a\ \mathcal {F}\left((4 - \alpha _i) \frac{\Delta a}{2}\right). \end{eqnarray*}$$

(9)

By combining equations (6) and (9) and recognizing that the average grain mass in bin i is m_i = ρ_i/n_i, we derive an expression that solely depends on the power-law index α_i, i.e.

$$\begin{eqnarray*} \frac{m_i}{m(a_{i+1/2})} \mathcal {F}\left((1 - \alpha _i) \frac{\Delta a}{2}\right) - \mathcal {F}\left((4 - \alpha _i) \frac{\Delta a}{2}\right) = 0 . \end{eqnarray*}$$

(10)

This expression needs to be solved numerically using a root-finding algorithm, for example the Newton–Raphson method. As |$\mathcal {F}$| in equation (10) is a monotonic function, only a few iterations are needed to find the solution, especially if the initial guess is close to the root. Once α_i is determined, the value of A_i can be directly calculated from equation (6) or equation (9).

The power-law description can be compared to methods previously used in the literature such as the piecewise-constant and piecewise-linear ones (e.g. Mizuno et al. 1988; Jones et al. 1994; Jones et al. 1996; Liffamn & Clayton 1989; McK18). The piecewise-constant discretization takes on a constant value for the distribution in bin i according to

$$\begin{eqnarray*} \left.\frac{\partial {n(a, t)}}{\partial {a}}\right|_i = \frac{n_i(t)}{(a_{i+1} - a_{i})}, \end{eqnarray*}$$

(11)

where we have chosen here that the constant reflects the total number density of grains in the bin. Alternatively, one can choose the distribution to reflect the mass density of the grains in the bin. A clear disadvantage is that this method only accurately describes the number density or the mass density, but not both. The piecewise-linear method of McK18 fixes this by assuming a linear distribution around the bin’s mid-point a_c,i = (a_i + a_{i + 1})/2

$$\begin{eqnarray*} \left.\frac{\partial {n(a, t)}}{\partial {a}}\right|_i = \frac{n_i(t)}{(a_{i+1} - a_{i})} + s_{i}(t)(a - a_{c,i}), \end{eqnarray*}$$

(12)

where the slope s_i(t) is chosen so that the mass density in the bin is equal to ρ_i. Note, however, that the linear distribution can become negative and non-physical if the slope is too steep. This is remedied by imposing a slope limiter ensuring positivity of the distribution function and conserving grain mass density (see Section 3.2.1 of McK18). Unfortunately, this also implies that the grain numbers within the bin are not conserved. The piecewise-constant and piecewise-linear methods can be considered to be first and second-order approximations to the power-law, respectively. The accuracy depends on the bin size and on the distribution that needs to be modelled. For example, if the distribution is flat within the bin, all three methods give identical results as s_i = 0 in the piecewise-linear method and α_i = 0 in the power-law method.

In describing the methods we have implicitly assumed that the grain distribution fills an entire bin. This does not necessarily need to be true, especially near the limits of the distribution r_min and r_max (where a_min ≤ r_min < r_max ≤ a_max). Then it is possible that the distribution is skewed towards one of the bin edges. In the piecewise-linear method of McK18, this causes the distribution function to become negative within the bin and slope limiting is required. On the other hand, the power-law method always guarantees positivity, but, unfortunately, a skewed distribution produces a large power-law index resulting, numerically, in a floating point error. To take into account the possibility that bin i is only partly populated, we set |$a^*_i = \max (a_i, r_{\rm min})$| and |$a^*_{i+1} = \min (a_{i+1}, r_{\rm max})$| as bin limits. Furthermore, we need to take |$\Delta a^* = \log (a^*_{i+1}/a^*_i)$| and |$a^*_{i+1/2} = a^*_i e^{\Delta a^*/2}$| to determine A_i and α_i. This small modification avoids floating point errors and conserves both grain mass and numbers. Note that this also means we are not restricted to logarithmically uniform bin widths, but can have randomly sized bins.

While using Δa* in the root finding algorithm for equation (10) improves the conservation properties of the distribution function, it also highlights a concern when Δa* becomes small, i.e. the root finding algorithm does not find a unique solution for α_i, or finds no solution at all. This is due to the shape of the function |$\mathcal {F}$|⁠. For small values of x, and thus Δa*, the function reduces to |$\mathcal {F}(x) \approx 1 + \frac{x^2}{3\!}$| and is numerically a constant for x < 1.5 × 10⁻³ (as the second term falls below the machine precision of 10⁻⁷). Given that Δa* is of the same order as x, we are limited to using the root finding algorithm for values of Δa* > 5 × 10⁻³. For values below this limit, we opt to use α_i = 0 and thus assume that the distribution is constant within the bin.

2.2 Redistribution of grain numbers and mass

The evolution of an advected dust grain-size distribution can be expressed as (Tsai & Mathews 1995)

$$\begin{eqnarray*} \frac{\partial }{\partial t}\left(\frac{\partial n(a,t)}{\partial a}\right) + \nabla \cdot \left(\frac{\partial n(a,t)}{\partial a} {\bf v}\right) &=& - \frac{\partial }{\partial a}\left(\frac{\partial n(a,t)}{\partial a} \frac{\mathrm{ d}a}{\mathrm{ d}t}\right) \nonumber\\ &&+\, S(a, t), \end{eqnarray*}$$

(13)

where da/dt is the rate of change of the grain radius and S(a, t) a source or sink of grains. Note that, if the right-hand side terms are equal to zero, this just represents changes in the distribution due to advection. Therefore, physical processes that affect the grain-size distribution are described by the terms on the right-hand side of the expression. The first term represents processes that increase or decrease the grain radius and conserve the total grain numbers, e.g. sputtering and mantle accretion, while production and destruction processes are included in the second one. Of the latter processes we only focus on the ones that conserve mass, such as shattering and coagulation, as other processes like supernova dust production are straightforward to implement. In the following we use sputtering and shattering, which are relevant to C-type shocks, to illustrate the methods for the grain distribution evolution.

2.2.1 Number-conserving processes

In the ISM, the impact of neutral particles and ions on dust grains releases, or sputters, grain material such as Si, Mg, and O, at a rate (Tielens et al. 1994)

$$\begin{eqnarray*} \frac{\mathrm{ d} N_{ s}}{\mathrm{ d}t} = 2\pi a^2 n_{ p} u_{ p} Y_{\rm s}(u_{ p}), \end{eqnarray*}$$

(14)

where N_s is the number of sputtered species, n_p the number density of impact particles, u_p the relative speed between the impacting particles and the grains and Y_s(u_p) the sputtering yield for species s integrated over all impact angles and evaluated for an impact velocity of u_p. Then the change in rate of the grain radius can be determined using the mass of the sputtered species m_s and is given by

$$\begin{eqnarray*} \frac{\mathrm{ d}a}{\mathrm{ d}t} = - \sum _{s,p}{\frac{m_s n_p}{2 \rho _g} u_p Y_{\rm s}(u_p)}. \end{eqnarray*}$$

(15)

Note that da/dt does not explicitly depend on the grain radius, but it does have a weak dependence through the sputtering yield. This is because the relative speed between the impinging neutral or ion species and the grains is grain radius-dependent (e.g. van Loo et al. 2009) and, for small grains, the projectile particle may be able to pass through the grain, therefore reducing the sputtering yield at these sizes (Bocchio, Jones & Slavin 2014).

When da/dt (or |$\dot{a}$|⁠) is a constant the time evolution of the grain distribution simply reduces to

$$\begin{eqnarray*} \frac{\partial n(a,t+\Delta t)}{\partial a} = \frac{\partial n(a-\dot{a}\Delta t, t)}{\partial a}, \end{eqnarray*}$$

(16)

and it is possible to split the effect of number-conserving processes from equation (13). Then the number density distribution at time t + Δt for bin i is given by

$$\begin{eqnarray*} n_i(t+\Delta t) &=& \int \nolimits _{a_i}^{a_{i+1}}{\left.\frac{\partial {n(a, t+\Delta t)}}{\partial {a}}\right|_i \,\mathrm{ d}a} \nonumber\\ &=& \sum _{j=0}^{N-1}{\int \nolimits _{a^j_i}^{a^{j+1}_{i+1}}{\left.\frac{\partial {n(a-\dot{a}\Delta t, t)}}{\partial {a}}\right|_j \,\mathrm{ d}a}}, \end{eqnarray*}$$

(17)

where |$a^j_{i} = \max [a_i, a_j+\dot{a}\Delta t]$| and |$a^{j+1}_{i+1} = \min [a_{i+1}, a_{j+1} + \dot{a}\Delta t]$|⁠. Thus, to determine the evolved number density in bin i, it is only necessary to determine from which size bin j the dust grains now residing in bin i come from. This can be done simply by calculating the position of the edges of bin j at time t + Δt, i.e. |$[a_j + \dot{a} \Delta t, a_{j+1} + \dot{a} \Delta t]$|⁠, and establishing which bins overlap with bin i. The contribution to the number density is then worked out analytically using A_j and α_j describing the power-law distribution in bin j at time t. Similarly, we can determine the updated volume density using

$$\begin{eqnarray*} \rho _i(t+\Delta t) &=& \int \nolimits _{a_i}^{a_{i+1}}{m(a)\left.\frac{\partial {n(a, t+\Delta t)}}{\partial {a}}\right|_i \,\mathrm{ d}a},\nonumber\\ &=& \sum _{j=0}^{N-1}{\int \nolimits _{a^j_i}^{a^{j+1}_{i+1}}{m(a)\left.\frac{\partial {n(a-\dot{a}\Delta t, t)}}{\partial {a}}\right|_j \,\mathrm{ d}a}}. \end{eqnarray*}$$

(18)

From the updated n_i(t + Δt) and ρ_i(t + Δt) we can solve for the power-law coefficient A_i and the index α_i at time t + Δt to find the discrete distribution function |$\frac{\partial n(a, t+\Delta t)}{\partial a}|_i$|⁠.

2.2.2 Mass-conserving processes

In contrast to sputtering, shattering due to grain–grain collisions conserves the total mass of the grain distribution but not their total number. Above a threshold impact velocity, some volume fraction of the grains involved will fragment into smaller dust grains which themselves follow a size distribution, that is |$\frac{\partial N_{\rm frag}}{\partial a} \propto a^{-3.3}$| (Jones et al. 1996). The evolution of the grain distribution is then described as (e.g. Jones et al. 1994)

$$\begin{eqnarray*} S(a,t) &=& -\frac{\partial n}{\partial a} \int \nolimits _{a_{\rm min}}^{a_{\rm max}}{\mathrm{ d}\mathrm{ }a_1 \frac{\partial n}{\partial a_1} \alpha (a, a_1)} + \frac{1}{2}\int \nolimits _{a_{\rm min}}^{a_{\rm max}}{\mathrm{ d}a_{1}} \frac{\partial n}{\partial a_1} \nonumber\\ &&\times \,\int \nolimits _{a_{\rm min}}^{a_{\rm max}}{\mathrm{ d}a}_{2} \frac{\partial n}{\partial a_2} \alpha (a_1, a_2) \frac{\partial N_{\rm frag}}{\partial a}(a, a_1, a_2), \end{eqnarray*}$$

(19)

where α(a₁, a₂) = π(a₁ + a₂)²v_rel(a₁, a₂), when multiplied by the grain number density is the collision frequency of grains with radius a₁ and a₂ above a threshold velocity for shattering and, otherwise, is equal to zero. |$\frac{\partial N_{\rm frag}(a, a_1, a_2) }{\partial a} \,\mathrm{ d}a$| is the number of grains with radii in the range [a, a + da] produced by interactions of grains with radius a₁ and a₂. Note that the first term describes the removal of dust grains from the interval and the second term the contribution to it due to fragmentation and requires integration over the entire grain-size distribution.

For the purpose of evolving a discrete distribution, equation (19) needs to be integrated over the different bins. Hence, the change of the number density as a function of time for bin i is given by

$$\begin{eqnarray*} S_i(t) &=& \int \nolimits _{a_i}^{a_{i+1}}{\mathrm{ d}a S(a,t)} \nonumber\\ &=& -\sum _{j=0}^{N-1}{\pi v_{\rm rel}(\langle a\rangle _i,\langle a\rangle _j)\ n_i n_j \left\lbrace \langle a^2\rangle _i + 2\langle a\rangle _i\langle a\rangle _j + \langle a^2\rangle _j\right\rbrace }\nonumber\\ &&+\, \frac{1}{2}\sum _{j=0}^{N-1}{\sum _{k=0}^{N-1}{\pi v_{\rm rel}(\langle a\rangle _j,\langle a\rangle _k)\ n_j n_k N_{{\rm frag},i}^{j,k}}}\nonumber\\ && \times \, \lbrace \langle a^2\rangle _j + 2\langle a\rangle _j\langle a\rangle _k + \langle a^2\rangle _k\rbrace , \end{eqnarray*}$$

(20)

where

$$\begin{eqnarray*} \langle a^l \rangle _i = \frac{1}{n_i}\int \nolimits _{a_i}^{a_{i+1}}{\mathrm{ d}a\ a^l \frac{\partial n}{\partial a}}, \end{eqnarray*}$$

(21)

where l is an integer, and

$$\begin{eqnarray*} N_{{\rm frag},i}^{j,k} = \int \nolimits _{a_i}^{a_{i+1}}{\mathrm{ d}a \frac{\partial N_{\rm frag}}{\partial a} (a, \langle a\rangle _j, \langle a\rangle _k)} \end{eqnarray*}$$

(22)

is the number of grains with sizes between [a_i, a_{i + 1}] due to fragmentation by collisions of grains within bin j and k. Here, we assumed that the distribution of grain fragments is the same for all grains in bin j and k, i.e. |$\frac{\partial N_{\rm frag}}{\partial a} (a, a_1, a_2) = \frac{\partial N_{\rm frag}}{\partial a} (a, \langle a\rangle _j, \langle a\rangle _k)$|⁠. If we know the analytic form of the size distribution of fragments and its dependence on the projectile and target radii, a more accurate version of equation (20) can be derived. Furthermore, it is presumed that all grains in a size bin have the same velocity and, thus, the relative velocity between two bins is also constant. Using equation (20), the number density in bin i at time t + Δt is then

$$\begin{eqnarray*} n_i(t+\Delta t) = n_i(t) + S_i(t) \Delta t. \end{eqnarray*}$$

(23)

Likewise, the volume density can be updated using

$$\begin{eqnarray*} \rho _i(t+\Delta t) = \rho _i(t) + S^{\prime }_i(t) \Delta t, \end{eqnarray*}$$

(24)

where |$S^{\prime }_i(t)$| can be derived from multiplying equation (19) with m(a) and then discretizing the integrals. We then find

$$\begin{eqnarray*} S^{\prime }_i(t) &=& -\frac{4\pi \rho _g}{3}\langle a^3\rangle _i \sum _{j=0}^{N-1}{\pi v_{\rm rel}(\langle a\rangle _i,\langle a\rangle _j)\ n_i n_j}\nonumber\\ &&\times \,\lbrace \langle a^2\rangle _i + 2\langle a\rangle _i\langle a\rangle _j + \langle a^2\rangle _j\rbrace \nonumber\\ &&+\, \frac{1}{2}\sum _{j=0}^{N-1}{\sum _{k=0}^{N-1}{\pi v_{\rm rel}(\langle a\rangle _j,\langle a\rangle _k)\ n_j n_k m_{{\rm frag},i}^{j,k}}}\nonumber\\ &&\times \, \lbrace \langle a^2\rangle _j + 2\langle a\rangle _j\langle a\rangle _k + \langle a^2\rangle _k\rbrace . \end{eqnarray*}$$

(25)

Here the mass transferred to bin i due to collisions of grains within bins j and k is given by

$$\begin{eqnarray*} m_{{\rm frag},i}^{j,k} = \int \nolimits _{a_i}^{a_{i+1}}{\mathrm{ d}a\ m(a) \frac{\partial N_{\rm frag}}{\partial a}(a, \langle a\rangle _j, \langle a\rangle _k)}. \end{eqnarray*}$$

(26)

This means that the radius of the fragmented grains is not taken into account, but only an appropriate mass for all grains within a bin is assumed. Note that this assumption must also be reflected in the mass-loss term, i.e. the first term on the left-hand side, as otherwise a systematic discrepancy arises between the mass-loss due to fragmentation and the redistributed mass across the distribution. Again, such simplifications are not needed when we know the analytic expression of the fragment distribution in terms of the radii of the colliding grains. Equations (20) and (25) are analogous to the expressions of other authors who have used either a piecewise-constant or linear description for the discrete distribution function (e.g. Mizuno et al. 1988; Jones et al. 1994; Jones et al. 1996; Hirashita & Yan 2009, McK18).

3 TESTS AND RESULTS

To test the power-law description of the grain distribution, we apply the methods of Section 2 to the test problems outlined in McK18. As these tests have analytical solutions, this allows a direct analysis of the performance of the method, but also a direct comparison with both the piecewise-constant and linear methods studied in McK18. Note that these tests do not necessarily represent physical or realistic situations.

3.1 Sputtering of a boxcar distribution

Here we will test the convergence of the error in the total grain mass depending on the number of size bins used. McK18 show that the piecewise-linear method exhibits a 1/N² scaling of the convergence and, thus, second-order behaviour, which is an improvement of the piecewise-constant method that is only first order.

The initial distribution is taken to be a boxcar function

$$\begin{eqnarray*} \frac{\partial n}{\partial a} = \left\lbrace \begin{array}{l{@}{\quad}l} 1\, {\rm cm^{-4}}& \text{if} a_L \le a \le a_R\\ 0\, {\rm cm^{-4}}& \text{otherwise.} \end{array}\right. \end{eqnarray*}$$

(27)

where [a_L, a_R] = [a_min(a_max/a_min)^3/8, a_min(a_max/a_min)^1/2] and a_min and a_max are set to 0.001 µm and a_max = 1 µm, respectively. Contrary to McK18, who adopt a grain growth rate, we take a constant grain sputtering rate of |$\dot{a} = - 2.4\times 10^{-7} {\rm cm Gyr^{-1}}$| applied for a time of t = 5 Gyr in 100 equal time steps. A constant sputtering rate is used here to ensure that the test is analogous to that of McK18. In reality the sputtering rate is size-dependent via the sputtering yield (e.g. Bocchio et al. 2014). Grains which are sputtered to a size smaller than a_min are assumed to be too small to participate in further sputtering and are removed from the distribution. As sputtering only affects the grain mass, the final distribution is still a boxcar function between the limits |$[a_L + \dot{a}t, a_R+\dot{a}t]$|⁠.

Fig. 1 shows the fractional error in the total grain mass as a function of bin number (from N = 8 to N = 2048) for the piecewise-constant, piecewise-linear and power-law methods. Both the piecewise-constant and linear methods show their expected first- and second-order dependence, respectively, on bin size and the latter method outperforms the former. However, the power-law method surpasses both of these with an accuracy below 0.1 per cent over all bin numbers (N = 8−2048). Especially for a small number of bins we find that the accuracy of the power-law method is more than 4 orders of magnitude better than the two other methods. The linear method only reaches this accuracy for N = 512 and the piecewise-constant for N = 2048.

Figure 1.

Fractional error of the total grain mass as functions of number of bins, N, for an initial boxcar distribution affected by sputtering. The boxes show the results for the power-law (blue), piecewise-linear (green), and piecewise-constant method (red). The dashed lines show a 1/N scaling (red) and 1/N² scaling (green).

It is pertinent to understand where these differences in the fractional error between the methods come from. In principle all methods should describe the distribution equally well as, for example, the piecewise-linear method should reduce to the piecewise-constant method (see Section 2.1). Furthermore, for the power-law method, the power-law index is set to α = 0 for N = 2048 as the root-finding algorithm breaks down (see Section 2.1). This implies that it also reduces to the piecewise-constant method, yet it produces a result that is nearly two orders of magnitude better than the piecewise-constant method. The only difference is the treatment of the distribution edges. As the distribution evolves due to sputtering, it moves across bins, but does not necessarily continue to cover an entire bin at the distribution limits. However, the piecewise-constant method dictates that the grains are uniformly distributed in a bin and, likewise, the linear method uses slope-limiting to distribute the grains across the entire bin. This causes the discrete distribution to be smeared out at its edges (see Fig.   2). Only the power-law method follows the distribution edges and takes them into account when determining the distribution function in the bin. Modifying the piecewise-constant and linear methods to follow the distribution limits, as in the power-law method, leads to an improved accuracy, with the relative errors in the mass below 10⁻⁴ for all bin sizes. Note that the treatment of the distribution edges in the power-law method also produces the variations seen in the relative error as a function of the bin number.

Figure 2.

Initial boxcar distribution (solid black) with N = 128 evolved by applying a sputtering rate of |$\dot{a} = - 2.4\times 10^{-7} {\rm cm \ Gyr^{-1}}$| for 5 Gyr using the piecewise-constant (red), piecewise-linear (green), and power-law method (blue) and compared to the analytic distribution (black dotted).

3.2 Sputtering of a MRN distribution

While the boxcar distribution of the previous section shows that it is important to carefully treat the edges of the distribution, it is not representative of realistic grain-size distributions. In the ISM the size distribution for silicate and carbon grains is given by a power-law (see equation 1 Mathis et al. 1977). Therefore, in this test, the three methods are tested on their ability to follow the evolution of a power-law distribution affected by sputtering.

We initialise each bin between [a_min, a_max] with the number and mass density calculated using |$\frac{\partial {n(a)}}{\partial {a}}|_i = a^{-3.5}$|⁠. We assume the same a_min, a_max, sputtering rate and evolution time as in the previous section. While we already minimize the errors occurring at the distribution edges by completely filling the full grain-size range, we further use the modified piecewise-constant and linear methods as described in the previous section (that is, we track the distribution limits). As a result, the distribution is not affected by the issues arising when the edge of the distribution falls within a bin, as in the boxcar test, and all the differences are due to the ability of each method to describe the underlying grain-size distribution.

Fig. 3 shows the fractional error in the total grain mass at the final time for the different methods. While, for the boxcar distribution, the modified piecewise-constant and linear methods have relative errors of the order 10⁻⁵ for all bin sizes, this is no longer true. Especially, the modified piecewise-constant method shows a linear behaviour in the fractional error with large errors at small bin numbers, i.e. |$\gt\! 10{{\ \rm per\ cent}}$| for N ≤ 32. The modified piecewise-linear method is significantly better but still performs poorly at small bin numbers, i.e. N < 16. Only the power-law method consistently produces errors smaller than 10⁻³ for all bin sizes. However, note that, for large values of N, the modified piecewise-linear method is better than the power-law method as the latter reduces to the modified piecewise-constant method when the bin size becomes small (see Section 2.1).

Figure 3.

Fractional error of the numerical grain mass as functions of number of bins, N, for an initial MRN distribution affected by sputtering. The squares show the results for the power-law (blue) and the modified piecewise-linear (green) and piecewise-constant (red) methods.

The discrete distribution function for the evolved MRN distribution is shown in Fig. 4 for N = 8 and 128. From the figure it is clear that the power-law method describes the power-law distribution very well. The modified linear method does capture the analytic solution well at the small grain radii, but less so at the large radii where it needs to apply slope limiting. However, the modified linear method is always better than the modified piecewise-constant one when describing a power-law distribution and, as N increases, the modified piecewise-linear method converges to the power-law method. Eventually, the modified piecewise-constant method will also converge but only at much larger values of N. This is expected as a power-law distribution, such as in equation (4), can be approximated to second order as

$$\begin{eqnarray*} \left.\frac{\partial n(a)}{\partial a}\right|_i = A_i a_0^{-\alpha _i} - \alpha _i A_i a_0^{-\alpha _i - 1} (a - a_0), \end{eqnarray*}$$

(28)

where a₀ is a grain size in the interval [a_i, a_{i + 1}]. The piecewise-linear and piecewise-constant discretizations are expressed similarly and, thus, eventually converge as the bin size decreases. Note that the convergence is quicker for shallower power laws.

Figure 4.

MRN distribution evolved by applying a sputtering rate of |$\dot{a} = - 2.4\times 10^{-7} {\rm cm \ Gyr^{-1}}$| for 5 Gyr using the modified piecewise-constant (red) and piecewise-linear (green) methods and the power-law method (blue) and compared to the analytic distribution (black dotted) for N = 8 (top) and N = 128 (bottom) bins.

Although all methods are able to conserve the total grain mass and reproduce the final distribution to a high degree for N = 128 (see Figs 3 and 4), it is also useful to evaluate the distribution function at specific grain radii. For processes such as sputtering, as the grains shrink, the number of grains does not change. Fig. 5 shows the grain number normalised to their initial value for grains with an initial grain radius of 0.01 or 0.5 μm for the three methods over a time range of 5 Gyr. While the power-law method maintains a constant grain number for both grain radii, both the modified piecewise-constant and linear methods show errors of the order of 10–15 per cent. These errors do not remain constant but vary significantly over time with large discontinuities when the grains move from one bin to another. Note that the 0.01 μm grains move through many more bins than the 0.5 μm ones before they reach a_min and are removed from the model around 3.8 Gyr. Thus, the power-law model does not only preserve global properties of the distribution, but also the inter-bin characteristics, unlike the modified piecewise-constant and linear methods.

Figure 5.

Evolution of the number density as function of time for grains with a initial radius of 0.01|$\rm \mu m$| (dotted) and 0.5 |$\rm \mu m$| (solid) for the modified piecewise-constant (red) and piecewise-linear (green) methods and the power-law (blue) distributions with N = 128. The number densities are normalised to the analytic values.

3.3 Grain shattering

The previous tests dealt with grain sputtering, a process that conserves the total number of grains in the distribution while the total mass of grains in the distribution changes. Here we look at grain shattering, in which the total mass of grains is conserved but the number of grains is altered significantly due to the production of many small fragments. When two grains of differing sizes, and different velocities, collide at a relative velocity exceeding a threshold value, some portion of the grains are fragmented. These smaller fragments can themselves be treated as spherical grains that follow a power-law grain-size distribution.

Here we carry out the same shattering test as in McK18. Only the collision between large grains (≥0.1 μm) causes fragmentation, with both grains completely destroyed, and, for simplicity, the fragments are distributed across the full size range of [a_min, a_max], where |$a_{\rm min} = 0.001\,{\rm \mu m}$| and |$a_{\rm max} = 1 {\rm \,\mu m}$|⁠, following a size-distribution ∝ a^−3.3. Note that, as fragments can be larger than the fragmenting grains, this model does not consider only shattering but also some degree of grain growth. The collision velocity between the grains is set to 3 kms⁻¹. With these assumptions we have

$$\begin{eqnarray*} \frac{\partial N_{\rm frag}}{\partial a}(a, a_1, a_2) = \frac{0.7\left(a_1^3 + a_2^3\right)}{a_{\rm max}^{0.7} - a_{\rm min}^{0.7}} a^{-3.3}, \end{eqnarray*}$$

(29)

resulting in

$$\begin{eqnarray*} N_{{\rm frag},i}^{j,k} = \frac{0.7}{2.3} (\langle a^3\rangle _j + \langle a^3\rangle _k) \frac{a_i^{-2.3} - a_{i+1}^{-2.3}}{a_{\rm max}^{0.7} - a_{\rm min}^{0.7}}, \end{eqnarray*}$$

(30)

and

$$\begin{eqnarray*} m_{{\rm frag},i}^{j,k} = (\langle m\rangle _j + \langle m\rangle _k) \frac{a_{i+1}^{0.7} - a_{i}^{0.7}}{a_{\rm max}^{0.7} - a_{\rm min}^{0.7}} , \end{eqnarray*}$$

(31)

where 〈m〉 = 4|$\pi$|ρ_g/3〈a³〉. Note that the latter expression is the same as in McK18 (their equation 63) and, to compare results, we will use the same expression. However, as the fragmentation distribution, equation (29) has a simple analytic expression, it is possible to derive more accurate, exact versions for equations (20) and (25). We will not give those expressions here, but we will use those to produce quasi-analytic solutions (using N = 256 bins) and to evaluate the approximations made in |$N_{{\rm frag},i}^{j,k}$| and |$m_{{\rm frag},i}^{j,k}$|⁠. To differentiate between the two power-law methods we will refer to the former as the default method and the latter as the exact method.We use the same initial conditions as in McK18, that is a lognormal distribution represented by a piecewise discrete distribution over N = 8 bins and assume that only the largest grains contribute to the shattering as v_rel(a_i, a_j) = 3 kms⁻¹ if a_i ≥ 0.1 μm and a_j ≥ 0.1 μm and is equal to zero otherwise. Fig. 6 shows the grain distribution due to shattering for a time of 150 Myr. At this time a reasonable amount of large grains have been shattered so that, while the initial distribution is only slightly modified at the large radii end of the distribution, the small grains follow the a^−3.3 distribution resulting from the fragmentation. The piecewise-linear method of McK18 describes this distribution reasonably well, especially if we evaluate the distribution at the geometric mid-point of the bins. Furthermore, quantitatively, the piecewise-linear routine produces a relative error in the total number density of about 10 per cent and conserves the total mass density exactly. However, a closer inspection of the distribution shows that the distribution is not adequately described, particularly at the bin edges. This is due to the slope limiting which needed to be performed at the small grain sizes (<0.1 μm) in order to ensure positivity of the distribution (while conserving mass). At the same time we see that the distribution of the large grains remains uniform within the bin (that is, a zero slope). This is because, in the method of McK18, the average grain size does not change if a bin loses mass, but only when it gains mass. While it is not a significant problem for this specific test where shattering is treated alone, reproducing the distribution shape becomes important when number-conserving grain processes are also considered.

Figure 6.

Distribution of grains due to shattering after 150 Myr for the piecewise-linear method (green) and the power-law method. The dashed blue line uses the default method with an approximated mass deposited in a bin, while the dotted red line uses the exact method using an exact mass calculation. The quasi-analytic solution is given by the black solid line and the initial condition by the dotted line.

The power-law method does describe the grain distribution across the full range of grain sizes more accurately with only minor deviations from the analytic solution at the large grain sizes. The exact method performs slightly better than the default one in reproducing the analytic distribution. Both methods have a relative error below machine precision for the mass density and below 2 per cent for the number density. Although the exact method does describe the distribution better than the approximate method, the errors are similar as a discrete distribution with N = 8 cannot adequately model the break in the analytic distribution at a = 0.1 μm. This break does not coincide with a bin edge, but instead falls in the middle of a bin. The benefits of the exact method become clear if we check the relative error in number density in each bin. The total number density is dominated by a single bin with the smallest grain radius, thus errors at the larger grain radii are not quantified by the relative error in the total grain number density. In the bin with the largest radii the relative error is below 0.1 per cent for the exact method, 8 per cent for the approximate method and reaches 25 per cent for the linear method. The accuracy at which the exact power-law method can reproduce the distribution in a bin reflects in improved performance at longer evolution times. If we run the shattering test for a longer time, for example up to t = 1  Gyr, the relative error in the total number density increases to |$\approx 10{{\ \rm per\ cent}}$| for the default method, but only to |$\approx 5{{\ \rm per\ cent}}$| for the exact method. Thus, it is crucial that the shattering integrals include as much information as possible to minimize the effect of error on the redistribution of fragments across the grain sizes, especially if modelling both mass- and number-conserving processes.

3.4 Combined grain sputtering and shattering

While Sections 3.2 and 3.3 show that the power-law method performs well for number-conserving processes and mass-conserving processes individually, these processes often arise in combination. Therefore, we study here the combined effect of sputtering and shattering of grains on an initial MRN distribution.

In this test we will model the dust grain evolution as it occurs within a C-type shock front moving through a medium of n_H = 10⁶ cm⁻³ with a dust-to-gas ratio of 0.01. In this situation the dynamics of the grains is determined by the balance of Lorentz forces and collision forces with neutral particles. This results in an effective velocity difference between small grains that are coupled to the magnetic field and move with ions and electrons and the large grains moving with speeds close to that of the neutrals. Guillet et al. (2007) show that the grain radius at which this transition occurs is between ∼7.5 × 10⁻⁶ and 2.5 × 10⁻⁵ cm depending on the density of the gas. Here we assume a discontinuous transition at

$$\begin{eqnarray*} a_{t} = a_{\rm min}\left(\frac{a_{\rm max}}{a_{\rm min}}\right)^{\frac{3}{4}} \approx 1.17\times 10^{-7}{\rm cm }, \end{eqnarray*}$$

(32)

where a_min and a_max are the same distribution limits as used in all previous tests, and that the velocity difference is 15 kms⁻¹. (Note that the transition is always on a bin edge.) Hence, only the small grains, a < a_t, will experience non-thermal sputtering due to neutral species, while shattering is due to collisions of small grains with large grains. For simplicity, we apply the same shattering procedure as in the previous section, that is both grains completely shatter and the fragments are distributed across the full range of grain sizes. Using equation (15) we can estimate the rate at which the grain radius decreases, i.e. we find |$\frac{\mathrm{ d}a}{\mathrm{ d}t}\approx -10^{-12}~{\rm cm~s^{-1}}$|⁠. We evolve the distribution for 10⁶ s. As there is no analytic solution for this problem, we assess the results using the converged solution for the distribution function as the bin number increases. We find that both the linear and power-law methods converge to the same solution.

Fig. 7 shows the grain distribution for the piecewise-linear and the power-law methods using N = 8 bins. Comparing the results with the initial MRN distribution we find that sputtering changes the slope of the distribution towards the small grains while shattering changes it for the large grains, an effect previously noted by e.g. Bocchio et al. (2014), Bocchio et al. (2016), Kirchschlager et al. (2019). Both the power-law and piecewise-linear methods are close to the converged distribution function, although the linear method is affected by slope limiting to ensure positivity of the distribution function. Slope limiting conserves mass, but not numbers, and this is reflected in the error relative to the converged distribution function. The linear method has a relative error in the total mass of only 2 per cent but it is 7 per cent in the total numbers. In comparison, the relative errors for the power-law method are below 1 per cent even for N = 8 bins. However, the piecewise-linear method does converge quickly and achieves the same accuracy with N = 32 bins.

Figure 7.

Final distribution after sputtering and shattering are applied in combination for a time of 10⁶ s with the sputtering rate −10⁻¹² cm s⁻¹ for the piecewise-linear (green) and power-law (blue) methods (N = 8). The initial MRN distribution is also shown as well as the converged solution (dashed magenta, N = 256).

The results of this test depend on the relative strength of the sputtering and the shattering. Therefore, as the linear and power-law methods perform differently for number-conserving and mass-conserving processes, we also perform the test with a sputtering rate an order of magnitude larger and smaller. Fig. 8 shows the grain distributions for these two additional models. For the higher sputtering rate, we find that the evolution of the distribution is dominated by sputtering. The sputtering removes more small grains from the distribution compared to the model with the default rate. Hence, less projectiles are available to shatter the large grains and, consequently, the distribution function at large grain radii does not evolve as much. Both the piecewise-linear and power-law methods with N = 8 bins are close to the converged solution with the relative error in the total mass about 2 per cent for the linear method and 0.1 per cent for the power-law method. However, the error in the total numbers is up to 20 per cent for the linear method while it is less than 1 per cent for the power-law method. The linear method achieves the same performance as the power-law method for N = 64 bins. For the lower sputtering rate, the evolution is mainly due to shattering as the number of large grains drops significantly. Sputtering does not remove many grains at the small radii so more projectile grains can collide with the large grains and shatter them. Again both methods reproduce the converged distribution very well, and this also is revealed in the relative errors. The relative error in the total grain number for the linear method (which is always the largest error) is only 3 per cent.

Figure 8.

Same as Fig. 7 but with a sputtering rate of −10⁻¹¹ cm s⁻¹ (top) and −10⁻¹³ cm s⁻¹ (bottom).

This test shows that the power-law method maintains its high level of accuracy for N = 8 bins even when combining number-conserving and mass-conserving processes. To achieve the same level with the piecewise-linear method, one needs to model the distribution with more bins, i.e. N > 32. Additionally, as with the previous tests, the power-law method is also able to produce the correct shape of the final distribution with only N = 8 bins, something which the linear method has been unable to achieve across all tests. This further enforces the usefulness of the power-law method for following the evolution of a grain-size distribution when limited computational resources are required or necessary.

4 IMPLEMENTATION IN HYDRODYNAMICAL CODES

Although the focus of this paper is on following the evolution of a grain-size distribution due to grain processing, our goal is to use this discrete power-law prescription in hydrodynamical simulations. To incorporate the grain-size distribution into, for example, a multifluid MHD code such as the one of Van Loo et al. (2013) the equations for the grain fluids need to be altered.

The starting point for this modification are the dynamical equations for each individual (pressureless) grain fluid which, in a weakly ionized plasma, are given by the continuity equations and the reduced momentum equation

$$\begin{eqnarray*} &&{\frac{\partial n}{\partial t} + \nabla \cdot \left(n {\bf v}\right) = S(a,t) + S_{\rm sputt}(a,t),}\nonumber\\ &&{\frac{\partial \rho }{\partial t} + \nabla \cdot \left(\rho {\bf v}\right) = S^{\prime }(a,t) + S^{\prime }_{\rm sputt}(a,t),}\nonumber\\ &&{\alpha \rho \left({\bf E} + {\bf v} \times {\bf B}\right) + \rho \rho _{\mathrm{ n}} K_{\mathrm{ gn}} ({\bf v}_\mathbf{ n} - {\bf v}) = 0,} \end{eqnarray*}$$

(33)

where S(a, t) is the grain shattering loss term given by equation (19) for the number density, S′(a, t) the shattering loss term for the mass density, S_sputt(a, t) and |$S^{\prime }_{\rm sputt}(a,t)$| the corresponding terms for grain sputtering, α = Ze/m the charge-to-mass ratio of the grains, E and B the electric field and magnetic field of the medium, and ρ_n and v_n the neutral mass density and velocity. Also, K_gn is the collision coefficient between the grain and the neutrals given by (Draine 1986)

$$\begin{eqnarray*} K_{\mathrm{ gn}} = \frac{8}{3} \frac{\pi a^2}{m_\mathrm{ n} + m}\left(\frac{2k_\mathrm{ B}T_\mathrm{ n}}{\pi m_\mathrm{ n}}\right)^{1/2} \left(1 + \frac{9\pi ({\bf v}_\mathrm{ n} - {\bf v})^2}{128k_\mathrm{ B}T_\mathrm{ n}}\right)^{1/2}, \end{eqnarray*}$$

(34)

with T_n the temperature of the neutral gas and m_n the mass of a neutral particle. Since we are modelling grains in a weakly ionized plasma, the collision frequency between grains and charged particles is negligible and we need only consider grain-neutral collisions. As n and ρ in equation (33) are effectively |$\frac{\partial n}{\partial a}\mathrm{ d}a$| and |$m(a) \frac{\partial n}{\partial a}\mathrm{ d}a$|⁠, we can find the governing equations by integrating the above equations over the range of radii for a given bin. For bin i, these then become

$$\begin{eqnarray*} &&{\frac{\partial n_{i}}{\partial t} + \nabla \cdot \left(n_{i} \bar{{\bf v}}\right) = S_i + S_{i,{\rm sputt}},}\nonumber\\ &&{\frac{\partial \rho _{i}}{\partial t} + \nabla \cdot \left(\rho _{i} \bar{{\bf v}}\right) = S^{\prime }_i + S^{\prime }_{i, {\rm sputt}},}\nonumber\\ &&{\langle Z\rangle _i e n_i \left({\bf E} + \bar{{\bf v}} \times {\bf B}\right) + n_i \rho _{\mathrm{ n}} K^*_{\mathrm{ gn}} ({\bf v}_\mathrm{ n} - \bar{{\bf v}}) = 0,} \end{eqnarray*}$$

(35)

where S_i and |$S^{\prime }_i$| are given by equations (20) and (25), S_{i, sputt} and |$S^{\prime }_{i,{\rm sputt}}$| represent the sputtering losses in bin i and 〈Z〉_ie the average grain charge. Furthermore,

$$\begin{eqnarray*} K^*_{\mathrm{ gn}} = \frac{8}{3} \pi \langle a^2\rangle _i \left(\frac{2k_\mathrm{ B}T_\mathrm{ n}}{\pi m_\mathrm{ n}}\right)^{1/2} \left(1 + \frac{9\pi ({\bf v}_\mathrm{ n} - \bar{{\bf v}})^2}{128k_\mathrm{ B}T_\mathrm{ n}}\right)^{1/2} \end{eqnarray*}$$

(36)

is the mean specific collision coefficient between neutrals and grains in bin i.

Note that, in order to derive these expressions, the only assumption we have made is that all grains within a bin have the same velocity. This is similar to the premise made for the collision frequency in the shattering process (see Section 2.2.2). However, from the reduced momentum equation, it is clear that the grain velocity depends on the grain radius through the Hall parameter (the ratio of the gas-grain collision frequency to the gyrofrequency) β = ZeB/mρ_nK_gn ∝ a⁻¹. For small Hall parameters, that is β < 0.1, the grains move with the neutrals, while, for β > 2, they move with the electrons and ions. Thus, there is only a small range of β, or grain radius, for which grains have a velocity in between and only in the bin where this transition occurs can some error in the dynamics be expected. In a subsequent paper studying grain processing in C-type shocks, we will analyse this further.

To deal with the grain processing due to shattering and sputtering in a hydrodynamical code, the routines described in Section 2 are used. Shattering can be included as a source function during the advection update, while we have the option to do the same for the sputtering or to operator split. The latter is preferred as the method described in Section 2.2.1 already gives the updated grain distribution function. Optimally, the operator split is done using Strang splitting with half a time-step before the advection update and half of a time-step after it.

5 DISCUSSION AND CONCLUSIONS

In this paper, we present a numerical method to follow the evolution of a dust grain-size distribution undergoing grain processes which either conserve grain mass or grain numbers. Guided by observations of typical ISM dust distributions, our method uses a power-law prescription to specify the distribution within a bin. Using the number and mass density of grains within a bin, the coefficient and index of the power law can be uniquely determined. We also explicitly track the grain size limits of the distribution. Furthermore, we describe the methods to evolve the discrete power-law distribution due to number-conserving or mass-conserving grain-processes and illustrate this with grain sputtering and grain shattering. The power-law method is complementary to the methods employing either a discrete piecewise-constant or piecewise-linear distribution (e.g. Mizuno et al. 1988; Jones et al. 1996; McKinnon et al. 2018).

The tests performed here show that the power-law method significantly outperforms both the piecewise-constant and linear methods for following the evolution of the distribution function, especially when the distribution is covered by a small number of bins. The main reason is, of course, that the discrete power-law method is naturally suited to modelling a continuous power-law distribution as often occurs in the ISM. The linear and the constant method only provide a second order and a first order approximation, respectively. In part it is also because we follow the pper limits and take th of the distributionem into account when deriving the distribution properties. This is important when considering number-conserving processes and the full radius range of a bin is not filled. For these processes, both the piecewise-constant and linear methods then diffuse the distribution limits. By implementing the same technique as in the power-law method, the relative errors can be reduced in the other two.

The power-law method is more effective for treating mass-conserving processes than the other methods, with the best results occurring when more information of the physical processes (in our case grain shattering) is included when evaluating the integrals. All methods conserve mass to machine precision accuracy, but the number density of the grains is better reproduced with the power-law method. While uncertainties are expected when modelling physical processes, it is best to avoid numerically induced ones. As mass- and number-conserving processes are often modelled together, the combined shattering and sputtering test demonstrated that the power-law method will provide the best results especially for small bin numbers (that is N = 8). For larger bin numbers both the power-law and piecewise-linear methods produce similar results.

The aim of this work is to provide an efficient numerical method that describes the evolution of a dust grain distribution due to advection and grain physics accurately in large-scale simulations. To avoid a large demand in numerical resources, it is beneficial to cover the grain distribution with a minimal number of bins. As our power-law method produces very small errors even for N = 8, it is perfect to include this approach in a numerical hydrodynamics code. One drawback is that operations such as |$\rm pow()$|⁠, |$\rm \log ()$| and |$\rm \sinh ()$| are considerably more CPU expensive than linear operations. This is especially important when finding the root of equation (10). Using standard algorithms, the power-law method using N = 8 is only as fast as the linear method with N = 128. However, one can use alternative algorithms and approximations for these operations so that the CPU cost of the power-law method is only 1.5 times that of the piecewise-linear method for the same number of bins. Thus, the power-law method not only provides a more accurate, but also viable, alternative to the piecewise-linear method. The implementation of the grain-processing methods in a hydrodynamics code is straightforward and only needs minimal alterations to the equations for a single dust grain fluid. The main assumption is a constant velocity for all the grains within a bin. In a weakly ionised plasma this assumption is likely not to be restrictive. In a subsequent paper we will further investigate this assumption when modelling C-type shocks with a full dust grain distribution.

Although the focus has been on the evolution of dust grain-size distributions, the application of this method is not limited to dust alone. The power-law method can be used to follow the evolution of any distribution function, especially those which exhibit a power-law distribution. One possible field of application is, for example, the energy loss of a cosmic ray distribution due to synchrotron radiation or inverse-Compton scattering.

ACKNOWLEDGEMENTS

We thank the anonymous referee for their helpful comments that have improved the paper. We thank Tom Hartquist for the useful discussions on grain physics and Paola Caselli for her continued interest and support. RS thanks the Science and Technology Facilities Council (STFC) and the Center for Astrochemical Studies at the Max Planck Institute for Extraterrestrial Physics for funding this PhD project. SVL is supported by an STFC consolidated grant. The data presented in this paper are available at https://doi.org/10.5518/751.

Footnotes

1

We note that shattering is not always strictly mass-conserving, since fragments may be produced which are smaller than the minimum grain radius limit of the distribution.