Evolution of dust porosity through coagulation and shattering in the interstellar medium

The evolution of grain size distribution and grain porosity by coagulation is formulated by O09 and O12. Although their main target was dust evolution in protoplanetary discs, the generality of their formulation allows us to apply it to interstellar dust. We also include shattering in this paper. The evolution of grain size distribution and grain porosity by coagulation and shattering is described by the distribution functions of grain mass ($m$) and grain volume ($V$)¹¹1We use the upper case $V$ for volumes and the lower case $v$ for velocities throughout this paper. at time $t$, and is governed by the 2D Smoluchowski equation. However, as noted by O09, solving the 2D Smolchowski equation requires a huge computational expense. Therefore, we concentrate on the moment equations for $V$, but we still solve the full distribution function for $m$.

To treat porous grains, it is convenient to introduce the following two types of grain radius: characteristic radius ($a_{\mathrm{ch}}$) and mass-equivalent radius ($a_{m}$). The characteristic radius is related to the volume as $V=(4\pi/3)a_{\mathrm{ch}}^{3}$, while the mass-equivalent radius is related to the grain mass as $m=(4\pi/3)a_{m}^{3}s$, where $s$ is the bulk material density. For compact grains, $m=sV$, so that $a_{\mathrm{ch}}=a_{m}$. For porous grains, $a_{\mathrm{ch}}>a_{m}$. Using these two grain radii, the filling factor, $\phi_{m}$, is expressed as $\phi_{m}\equiv(a_{m}/a_{\mathrm{ch}})^{3}(\leq 1)$ (see also Ormel et al., 2009), while the porosity is $1-\phi_{m}$. Since it is necessary to specify $s$, we simply adopt the same value as used in HM20 ($s=2.24$ g cm^-3 taken from graphite), but adopting $s=2$–3.5 g cm^-3 appropriate for interstellar grains (e.g. Weingartner & Draine, 2001b) does not change our conclusion significantly. We vary some parameters relevant for porosity, which are more important in this paper (Section 2.5). Thus, we simply use the same parameters as in HM20 for the material properties (specifically $s$ and $Q_{\mathrm{D}}^{\star}$ introduced later) that do not directly affect the porosity evolution.

Before presenting the basic equations, we need to introduce some quantities. We use the grain mass distribution instead of the grain size distribution since there is ambiguity in the grain size for porous grains. We denote the distribution function of $m$ and $V$ at time $t$ as $f(m,\,V,\,t)$. The grain mass distribution at time $t$, $\tilde{n}(m,\,t)$ is defined such that $\tilde{n}(m,\,t)\,\mathrm{d}m$ is the number density of grains whose mass is between $m$ and $m+\mathrm{d}m$, and is related to the above distribution function as

We also introduce the first moment of $f$ for $V$ as

This quantity ($\bar{V}$) is the mean volume of grains with mass $m$. We define the following two quantities, $\rho(m,\,t)\equiv m\tilde{n}(m,\,t)$ and $\psi(m,\,t)\equiv\bar{V}(m,\,t)\tilde{n}(m,\,t)$, which represent the grain mass distribution weighted for the grain mass and volume, respectively (the latter can also be regarded as the volume distribution function). The integration of $\rho(m,\,t)$ for $m$ is related to the dust-to-gas ratio, $\mathcal{D}$, as

where $\mu_{\mathrm{H}}=1.4$ is the gas mass per hydrogen, $m_{\mathrm{H}}$ is the hydrogen atom mass, and $n_{\mathrm{H}}$ is the hydrogen number density.

The Smolchowski equation contains the collision kernel, which is the product of the collision cross-section and the relative speed of colliding pair. The collision kernel for colliding grains with masses $m_{1}$ and $m_{2}$ is denoted as $K_{m_{1},m_{2}}$, which is evaluated in Section 2.2.

It is straightforward to derive the following general equations applicable for both coagulation and shattering based on equations (6) and (7) of O09:

where $\tilde{\theta}(m;\,m_{1},\,m_{2})$ describes the distribution function of grains produced from $m_{1}$ in the collision between grains with masses $m_{1}$ and $m_{2}$, and $(V_{1+2})_{m_{1},m_{2}}^{m}$ is the volume of the newly produced grain with mass $m$ in the above collision. Since the expression is not exactly the same as that used by O09 and O12, we explain the derivation of the above equations in Appendix A. In particular, in our formulation, we count the collisional products originating from $m_{1}$ and $m_{2}$ separately. O09, in contrast, considered the collisional pair $m_{1}$ and $m_{2}$ only once. This is why O09 has a factor 1/2 in front of the second term on the right-hand side in equations (4) and (5).

We aim at solving equations (4) and (5) for $\rho$ and $\psi$. Since $\rho/m=\psi/\bar{V}$ is always satisfied, $\bar{V}$ is not an independent quantity. Thus, these two equations are closed if we give $K$, $V_{1+2}$ and $\bar{\theta}$. In the following subsections, we explain how to determine these three quantities.

In computing the grain size distribution, we discretize the entire grain radius range ($a=3\times 10^{-4}$–10 $\micron$) into 128 grid points with logarithmically equal spacing. We set ${\rho}_{\mathrm{d}}(m,\,t)=0$ at the maximum and minimum grain radii for the boundary conditions. We use the algorithm described in appendix B of Hirashita & Aoyama (2019) to solve the discretized equations.

The collision kernel is determined by the product of the grain cross-section and the relative velocity. The cross-section ($\sigma_{1,2}$) in the collision of grains with masses $m_{1}$ and $m_{2}$, referred to as grain 1 and 2, respectively, is estimated as $\sigma_{1,2}=\pi(a_{\mathrm{ch1}}+a_{\mathrm{ch2}})^{2}$, where the characteristic radii of grains 1 and 2 are $a_{\mathrm{ch1}}$ and $a_{\mathrm{ch2}}$, respectively (see the beginning of Section 2.1 for the definition of the characteristic radius; i.e. $a_{\mathrm{ch}}=a_{m}\phi_{m}^{-1/3}$).

The motion of dust grains is assumed to be induced by turbulence (e.g. Kusaka et al., 1970; Voelk et al., 1980). We basically adopt the same simple model, originally taken from Ormel et al. (2009), for the grain velocity $v_{\mathrm{gr}}$ as adopted in our previous papers:

where $\mathcal{M}$ is the Mach number of the largest-eddy velocity, and $T_{\mathrm{gas}}$ is the gas temperature. This is similar to equation (18) of Hirashita & Aoyama (2019), but modified to include the filling factor $\phi_{m}$ (see appendix A of Ormel et al. 2009). The derivation of the above equation does not strictly hold for highly supersonic regime, but we practically use $\mathcal{M}$ here as an adjusting parameter for the grain velocity. Under fixed ISM conditions and grain material properties, the velocity is reduced to a function of $m$. Since the filling factor $\phi_{m}$ changes as a function of time, $v_{\mathrm{gr}}(m)$ also varies with time. In considering the collision rate between grains 1 and 2 with $v_{\mathrm{gr}}(m_{1})=v_{1}$ and $v_{\mathrm{gr}}(m_{2})=v_{2}$, we estimate the relative velocity $v_{1,2}$ by

where $\mu=\cos\theta$ ($\theta$ is an angle between the two grain velocities) is randomly chosen between $-1$ and 1 in every calculation of the collision kernel (Hirashita & Li, 2013).

Using the above quantities, the collision kernel is estimated for the collision between grains 1 and 2 as

Note that the collision kernel depends on the filling factor of the larger grain as $\propto\phi_{m}^{-1/3}$. Thus, the grain–grain collision rate is enhanced if the grains become more porous (fluffy).

The treatment of collisional products is the key for the evolution of porosity. The porosity is strongly affected by the production of vacuum volume in a grain after grain–grain sticking (coagulation) and the compaction associated with compression in grain–grain collisions. However, these processes are strongly nonlinear and dependent on the actual grain shapes. Moreover, some assumptions in O09 and O12 are not applicable to interstellar dust. In particular, their formulation is based on single-sized ($\sim 0.1~{}\micron$) monomers, since they are interested in protoplanetary discs, where grains grow into much larger sizes than interstellar grains. In the ISM, in contrast, we also consider the evolution of grains much smaller than the ‘typical’ grain size ($\sim 0.1~{}\micron$). To make the problem more complicated, shattering also occurs, producing a large number of small grains. Because of the above difficulty and the difference from O09 and O12, we take the following approach that simplifies the treatment of porosity evolution but still preserves the essence of physical processes regarding grain–grain sticking and compression.

To investigate the effects on observed dust properties, we also calculate extinction curves. We examine two representative grain materials: silicate and carbonaceous species, and investigate how the porosity evolution driven by coagulation and shattering affects the extinction properties of these materials. For silicate, we use the optical constants of astronomical silicate adopted by Weingartner & Draine (2001b). For carbonaceous dust, we adopt amorphous carbon (amC) from ‘ACAR’ in Zubko et al. (1996). Graphite is another possible carbonaceous material, and we confirmed that it produces similar results to amC except at wavelengths where the 2175 Å feature is prominent. We also found that the wavelength where the feature peaks shifts with the porosity, but such a shift is not observed in the Milky Way extinction curves. Thus, special care should be taken of the modelling of the 2175 Å feature, e.g. by treating the 2175 Å carrier such as graphite and PAHs (e.g. Li & Draine, 2001) as a component separated from porous grain species (Voshchinnikov et al., 2006), which is out of the scope of this paper. Thus, we adopt amC in this paper, noting that, as mentioned above, amC and graphite produce similar results at wavelengths not affected by the 2175 Å feature.

The optical properties of dust is calculated using the effective medium theory (EMT), which averages the dielectric permittivity using a mixing rule. We adopt the Bruggeman mixing rule. This mixing rule as well as the Garnett mixing rule gives reasonable extinctions as long as the grains do not have substructures (e.g. monomers) larger than the wavelength (Voshchinnikov et al., 2005). As shown later, the porous grains are formed by coagulation of grains smaller than $\sim 0.1~{}\micron$ and we are interested in wavelengths longer than 0.1 $\micron$. Thus, the above condition is satisfied in this paper. Even if we model the inhomogeneity using multi-layered spheres, the difference in the extinction is expected to be less than $\sim 10$ per cent (Voshchinnikov et al., 2005; Shen et al., 2008). Using the Bruggeman mixing rule, the averaged dielectric permittivity $\bar{\varepsilon}$ is obtained from

where $\varepsilon_{2}$ is the dielectric permittivity of the material (note that the first material is assumed to be vacuum so $\varepsilon_{1}=1$).²²2We adopt Gaussian-cgs units. We assume that each dust grain is a sphere with radius $a_{\mathrm{ch}}$ and refractive index $\bar{m}=\sqrt{\bar{\varepsilon}}$. The cross-section $C_{\mathrm{ext},m}$, which is a function of $m$ under a given $\phi_{m}$ at each epoch, is calculated by the Mie theory (Bohren & Huffman, 1983).

The extinction at wavelength $\lambda$ in units of magnitude ($A_{\lambda}$) can be calculated using the grain size distribution as

where $L$ is the path length. We present the extinction in the following two ways: $A_{\lambda}/N_{\mathrm{H}}$ and $A_{\lambda}/A_{V}$ (the $V$ band wavelength corresponds to $\lambda^{-1}=1.8\,\mu{\rm m}^{-1}$ and $N_{\mathrm{H}}=n_{\mathrm{H}}L$ is the column density of hydrogen nuclei). The first quantity indicates the extinction per hydrogen (thus, it is proportional to the dust abundance), while the second is useful to show the shape of extinction curve. In both quantities, the path length $L$ is canceled out.

Coagulation and shattering are governed by the same form of equation (equations 4 and 5). These two processes occur in different ISM phases. We assume that coagulation occurs in the dense clouds of which the physical conditions are described by $n_{\mathrm{H}}=10^{3}$ cm^-3 and $T_{\mathrm{gas}}=10$ K (typical of molecular clouds); and that shattering takes place in the diffuse ISM characterized by $n_{\mathrm{H}}=0.3$ cm^-3 and $T_{\mathrm{gas}}=10^{4}$ K. These values are similar to the ones adopted for coagulation and shattering by Hirashita & Yan (2009). The density affects the grain collision time-scale, which scales with $(n_{\mathrm{H}}v_{\mathrm{gr}})^{-1}\propto n_{\mathrm{H}}^{-3/4}$ (see equation 6); that is, the adopted duration for shattering/coagulation is degenerate with the gas density. Therefore, we fix the above physical conditions for the gas and examine the grain size distributions at various times (i.e. for various durations of the process), keeping in mind that the same value of $n_{\mathrm{H}}^{3/4}t$ gives similar results.

For the normalization of the grain velocities adjusted by $\mathcal{M}$ in equation (6), we apply $\mathcal{M}=3$ for shattering and $\mathcal{M}=1$ for coagulation (HM20). These values of $\mathcal{M}$ are adopted to obtain a similar level of grain velocities to those calculated for magnetized turbulence by Yan et al. (2004).

We not only treat shattering and coagulation separately, but also examine some cases where these two processes occur at the same time. The coexistence of shattering and coagulation is a reasonable approximation for a volume of the ISM wide enough to sample a statistically significant amount of both ISM phases and/or for a time much longer than the mixing time-scale of the two ISM phases ($\sim 10^{7}$ yr; e.g. McKee 1989). Since coagulation and shattering occur in the media with different densities, it is convenient to define the grain mass and volume distributions per hydrogen number density as $\tilde{\rho}\equiv\rho/n_{\mathrm{H}}$ and $\tilde{\psi}\equiv\psi/n_{\mathrm{H}}$, respectively. In this case, we calculate the evolution of $\tilde{\rho}$ and $\tilde{\psi}$ by

where $\eta_{\mathrm{dense}}$, which is treated as a free parameter, is the dense gas fraction determining the weight of each process, the subscript ‘tot’ indicates the total changing rates of $\tilde{\rho}$ and $\tilde{\psi}$, and the subscripts ‘coag’ and ‘shat’ mean the changes caused by coagulation and shattering, respectively. For the two terms on the right-hand side of the above equations, we apply equations (4) and (5). The change by coagulation and that by shattering are mixed with a ratio of $\eta_{\mathrm{dense}}:(1-\eta_{\mathrm{dense}})$; in other words, $\eta_{\mathrm{dense}}$ determines the fraction of time dust spends in the dense ISM. In calculating shattering and coagulation, we use $\rho$ and $\psi$, but when we sum up the contributions from coagulation and shattering, we divide them by $n_{\mathrm{H}}$ and obtain $\tilde{\rho}$ and $\tilde{\psi}$. Here we implicitly neglect gas which hosts neither coagulation nor shattering. Existence of such a gas component effectively lowers the efficiencies of both shattering and coagulation (or makes the time-scales of these two processes longer) equally and does not affect the relative roles of these two processes.

We examine the evolution of grain size distribution and filling factor by coagulation in the dense ISM (the pure coagulation case). We show the results in Fig. 1. Coagulation continuously forms large grains, depleting small grains. The bump in the grain size distribution at large grain radii is prominent after coagulation takes place significantly. The height of the bump is almost constant, reflecting the mass conservation in coagulation. Micron-sized grains are formed on a time-scale of 100 Myr, which is consistent with Hirashita & Voshchinnikov (2014, see their figure 4a). Note that this time-scale scales with the assumed density roughly as $\propto n_{\mathrm{H}}^{-3/4}$ as mentioned in Section 2.5 (recall that we adopt $n_{\mathrm{H}}=10^{3}$ cm^-3 for the dense ISM).

We also show the filling factor $\phi_{m}$ in Fig. 1. We observe that the filling factor steadily decreases up to $t\sim 100$ Myr because coagulation creates porosity. The decrease of filling factor is governed by coagulation of small grains, which are, however, effectively depleted by coagulation. As the abundance of small grains becomes lower, the decrease of the filling factor is saturated. For large ($a\gtrsim 0.1~{}\micron$) grains, the porosity is determined by the assumption of maximum compaction regulated by $\epsilon_{V}$ in equation (14). The filling factor at large grain radii roughly approaches $\phi_{m}\sim 1/(1+\epsilon_{V})\simeq 0.67$ (recall that $\epsilon_{V}=0.5$). In fact, the value is slightly larger than 0.67, because of our treatment of equation (14): There are still cases where $V_{\mathrm{1,comp}}+V_{\mathrm{2,comp}}<V_{1+2}<(1+\epsilon_{V})(V_{\mathrm{1,comp}}+V_{\mathrm{2,comp}})$. In this case, grains whose filling factors are between $1/(1+\epsilon_{V})$ and 1 form, contributing to raising the averaged $\phi_{m}$ above $1/(1+\epsilon_{V})$ at large grain radii.

We also examine the dependence on the parameters relevant for coagulation. First, we present the effect of $E_{\mathrm{roll}}$, which regulates compaction. As indicated in Section 2.3.1, $E_{\mathrm{roll}}$ is determined by the critical displacement $\xi_{\mathrm{crit}}$ in our model. In Fig. 2a, we show the results for $\xi_{\mathrm{crit}}=2$, 10, and 30 Å (note that $E_{\mathrm{roll}}$ is proportional to $\xi_{\mathrm{crit}}$). We only show the results at $t=100$ Myr since the effect of $\xi_{\mathrm{crit}}$ is qualitatively similar at all ages. We observe that the filling factor at $a\sim 0.1~{}\micron$ is affected by the change of $\xi_{\mathrm{crit}}$ (or $E_{\mathrm{roll}}$), with larger $\xi_{\mathrm{crit}}$ showing smaller $\phi_{m}$. This is because compaction is less efficient in the case of larger $\xi_{\mathrm{crit}}$. In contrast, the filling factor at $a\lesssim 0.01~{}\micron$ is not affected by $\xi_{\mathrm{crit}}$ because compaction is not important in that grain radius range. The filling factor also converges to the same value determined roughly by $1/(1+\epsilon_{V})$ at the largest grain radii as mentioned above. The grain size distribution, on the other hand, is insensitive to $\xi_{\mathrm{crit}}$. Since the grain abundance is dominated by large grains, small grains collide predominantly with large grains. In this situation, the collisional cross-section is governed by large grains which are almost compact. Thus, the difference in the porosity at submicron radii does not affect the grain size distribution.

Grain compaction is also affected by the number of contacts $n_{\mathrm{c}}$ (equation 12). In Fig. 2b, we show the results for $n_{\mathrm{c}}=10$, 30, and 100 at $t=100$ Myr. Naturally, the effect of $n_{\mathrm{c}}$ appears at grain radii where compaction is important. As $n_{\mathrm{c}}$ becomes larger, the filling factor is kept lower against compaction, so that the increase of $\phi_{m}$ becomes shallower at large $a_{m}$. The grain size distribution is insensitive also to the change of $n_{\mathrm{c}}$.

We also examine the dependence on $\epsilon_{V}$, which regulates the maximum compaction. As we argued above, the filling factor roughly converges to $1/(1+\epsilon_{V})$ at large grain radii, where grain velocities are high enough for significant compaction. As shown above, the minimum value of $\phi_{m}$ is around 0.3–0.5. Thus, $\epsilon_{V}$ should be smaller than $\sim 1$; otherwise, the filling factor at large grain radii decreases below the minimum value, which means that we add artificial (or contradictory) porosity to the grains in compaction. Thus, we examine $\epsilon_{V}=0$, 0.5, and 1 at $t=100$ Myr in Fig. 2c. We confirm that the effect of $\epsilon_{V}$ appears at large grain radii ($a\gtrsim 0.2~{}\micron$). As mentioned above, $\phi_{m}$ at large grain radii is roughly $1/(1+\epsilon_{V})$. The grain size distribution is affected by $\epsilon_{V}$ since larger porosity makes the collision kernel larger by a factor of $\phi_{m}^{-1/3}$ (Section 2.2), which increases the grain–grain collision rate.

We show the evolution of grain size distribution and filling factor by shattering in the diffuse ISM. Since shattering does not create new porosity in grains, we start with $\phi_{m}=0.5$ for all grains. (If we choose $\phi_{m}=1$ as the initial condition, the filling factor is always 1.) We show the results in Fig. 3a. Shattering continuously converts large grains to small grains, creating grains down to the minimum grain radius $a=3$ Å. After $t\sim 100$ Myr, the grain abundance is dominated by small grains. The filling factor increases, especially at small grain radii, because shattered fragments are compact. Since we assume that the remnants have the same filling factor as the original grains, the filling factors of the largest grains, which are dominated by shattered remnants, do not change. Recall that the largest fragment size is $(0.02)^{1/3}\simeq 0.27$ times the original grain size (Section 2.3.2). Since the maximum grain radius in the initial condition is 0.25 $\micron$, the largest grain radius where $\phi_{m}$ is raised by fragments is $a\simeq 0.27\times 0.25~{}\micron\simeq 0.068~{}\micron$. Thus, the porosity is only changed from its initial value at $a<0.068~{}\micron$. All the grains with $a<a_{\mathrm{min,ini}}=0.001~{}\micron$ are newly formed by shattering; thus, they always have $\phi_{m}=1$.

As explained in Section 2.5.2, the shattered remnants could suffer compaction. In order to examine this effect, we also show the case with compaction of remnants; that is, we adopt equation (24) instead of the first condition in equation (19) for the remnant volume. In Fig. 3b, we show the evolution of grain size distribution and filling factor with including remnant compaction. The grain size distributions are little affected by the treatment of the remnant volume, while the filling factor at large grain radii increases if we take compaction of remnants into account. Because shattering never creates porosity, the filling factors of all grains converge to unity on a time-scale of depleting large grains by shattering (100–300 Myr).

If we consider the grain size distribution in a region large enough to include both the dense and diffuse ISM (or if we consider time-scales much longer than the phase-exchange time; Section 2.5), the effect of coagulation and that of shattering coexist. In our one-zone model, we could simulate this coexisting effect by simultaneously treating coagulation and shattering at each time-step with a weight of $\eta_{\mathrm{dense}}:(1-\eta_{\mathrm{dense}})$ as formulated in equations (22) and (23). We adopt $\eta_{\mathrm{dense}}=0.5$ first as a fiducial value (that is, the grains spend half of their times in the dense ISM). In Fig. 4a, we show the evolution of grain size distribution and filling factor.

We observe in Fig. 4a that, if both coagulation and shattering are present, the grain size distribution maintains a power-law-like shape. It is also interesting to note that the slope is similar to that of the MRN distribution. It has been shown that the slope of grain size distribution converges to a value similar to the MRN distribution if efficient fragmentation and/or coagulation occur (e.g. Dohnanyi, 1969; Williams & Wetherill, 1994; Tanaka et al., 1996; Kobayashi & Tanaka, 2010). Coagulation is slightly ‘stronger’ than shattering, so that large (submicron) grains gradually increase. The grain size distribution and filling factor reach an equilibrium around $t\sim 100$ Myr. The filling factor changes in a way similar to the pure coagulation case (Fig. 1) but the decrease of the filling factor proceeds more if both coagulation and shattering are present because small grains which coagulate to form porosity are continuously supplied. The filling factor at large grain radii converges to $\sim 1/(1+\epsilon_{V})$ as explained in Section 3.1.

To examine the balance between coagulation and shattering, we also show the results for different values of $\eta_{\mathrm{dense}}$ (0.1 and 0.9) in Fig. 4. Naturally, a higher abundance of large grains is obtained for larger $\eta_{\mathrm{dense}}$. For $\eta_{\mathrm{dense}}=0.9$, the result is similar to that in the pure coagulation case (Fig. 1), except that small grains continue to be produced by shattering. The case of $\eta_{\mathrm{dense}}=0.1$ has the following two interesting properties: (i) The evolution of grain size distribution is not monotonic. Indeed, large grains are continuously converted into small grains up to $t\sim 100$ Myr but large grains increase (‘re-form’) after that. (ii) The filling factor becomes very small at $t>100$ Myr. This is because the production of small grains by shattering continuously activates coagulation, which creates porosity. The increased porosity makes the geometrical cross-section larger, and enhances coagulation. Moreover, since the velocities of porous grains are reduced, compaction does not occur as efficiently as in the pure coagulation case. Thus, efficient shattering together with coagulation makes the interstellar dust highly porous and further activates coagulation. This interesting interplay is further discussed in Section 5.

The dependences on the parameters related to coagulation ($\xi_{\mathrm{crit}}$, $n_{\mathrm{c}}$ and $\epsilon_{V}$) are similar to those shown in Section 3.1 (Fig. 2). On the other hand, compaction of shattered remnants (Section 2.5.2) could decrease the porosity if $\eta_{\mathrm{dense}}$ is low. In Fig. 5, we show the results with compaction of shattered remnants for $\eta_{\mathrm{dense}}=0.1$. This figure should be compared with Fig. 4b. We observe that the filling factor becomes higher at $t\sim 300$ Myr at large grain radii, although the results are almost unchanged at younger ages. Compaction of remnants also makes coagulation at later epochs less efficient, producing less large grains at $t\sim 300$ Myr. However, the filling factor is still as small as 0.1–0.2 at $a\sim 0.1~{}\micron$ at $t\sim 300$ Myr; thus, the conclusion that efficient shattering together with coagulation helps to increase porosity is robust.

In Fig. 6, we present the extinction curves corresponding to the grain size distributions and the filling factors for the pure coagulation case (Fig. 1). As expected from the evolution of grain size distribution, the extinction curve becomes flatter as coagulation proceeds. Since a large fraction of the grains become bigger than $\sim 0.1~{}\micron$ after coagulation, the extinction per hydrogen declines at wavelengths shorter than $\sim 2\pi\times 0.1~{}\micron$ for both materials. The effect of porosity, which appears in the difference between $A_{\lambda}$ (thick lines) and $A_{\lambda,1}$ (thin lines) or in the ratio $A_{\lambda}/A_{\lambda,1}$ shown in the bottom window, is clear. The extinction is 10–20 per cent higher in the porous case than in the non-porous case in a large part of wavelengths for both materials.

There are some differences between the two materials. For silicate, porous grains have smaller extinction than compact ones in a certain wavelength range (e.g. $1/\lambda\sim 0.5$–2 $\micron^{-1}$ at $t=10$ Myr) and this range shifts to longer wavelengths with age. Voshchinnikov et al. (2006) and Shen et al. (2008) showed that the porosity could both increase and decrease the extinction cross-section depending on the wavelength. The wavelength range of suppressed extinction by porosity is roughly consistent with their results. Voshchinnikov et al. (2006) also showed that the wavelength range where the extinction cross-section is decreased by porosity shifts to longer wavelengths as the grains become larger. This is consistent with our result. At $t\lesssim 100$ Myr, the normalized extinction curves ($A_{\lambda}/A_{V}$) becomes steeper by the porosity because the increase of extinction is more enhanced in the UV than in the $V$ band.

For amC, porosity always increases the extinction; that is, $A_{\lambda}/A_{\lambda,1}$ is always larger than 1. Thus, if we normalize the extinction to $A_{V}$, the porosity effect roughly cancels out, so that the extinction curve shape is insensitive to the porosity (filling factor) for amC.

We also examined the parameter dependence (not shown) and confirmed that the extinction curves are not sensitive to $\xi_{\mathrm{crit}}$ or $\epsilon_{V}$ in the optical and UV with differences less than 5 per cent. However, at $1/\lambda<2~{}\micron^{-1}$, the difference could be as large as 20 per cent, since, at such long wavelengths, the porosity of large grains, which is regulated by $\xi_{\mathrm{crit}}$ and $\epsilon_{V}$, is important. The changes driven by the difference in those parameters are sub-dominant compared with the difference between porous and non-porous grains shown in Fig. 6.

In the pure shattering case, the filling factor simply tends to increase; thus, the resulting extinction curve approaches the one with compact grains. More important is the interplay between shattering and coagulation as shown above. Indeed, coagulation of newly created small grains by shattering produces porous grains, contributing to the increase of porosity in the interstellar dust. Here, we investigate the effect of relative strength between shattering and coagulation; that is, we compare the results for different dense gas fractions, $\eta_{\mathrm{dense}}$.

In Fig. 7, we compare the extinction curves calculated for the models including both coagulation and shattering. The corresponding grain size distributions are shown in Section 3.3 (Fig. 4). We compare the results at the same age $t=100$ Myr. We observe that the steepness of extinction curve is very sensitive to $\eta_{\mathrm{dense}}$. This is a natural consequence of the different grain size distributions for various $\eta_{\mathrm{dense}}$. We also plot the extinction curves with $\phi_{m}=1$ (with the same distribution of $a_{m}$), that is, $A_{\lambda,1}$. Comparing $A_{\lambda}$ with $A_{\lambda,1}$, we observe for silicate that the porosity increases the extinction at short and long wavelengths, but that it rather decreases the extinction at intermediate wavelengths as already noted in the previous subsection. The wavelength range where the porosity decreases the extinction shifts to shorter wavelengths for smaller $\eta_{\mathrm{dense}}$, which is consistent with the above mentioned tendency that porosity of smaller grains reduces the extinction at shorter wavelengths. For amC, the porosity always increases the extinction, but the largest porosity (i.e. $\eta_{\mathrm{dense}}=0.1$) does not necessarily mean the largest opacity enhancement ($A_{\lambda}/A_{\lambda,1}$) in the UV. At long wavelengths, larger porosity indicates larger extinction.

Overall, the effect of porosity on the extinction curve is less than 20 per cent as far as the UV–optical extinction curves are concerned, but this could mean that we ‘save’ up to 20 per cent of metals to realize an observed extinction if we take porosity into account. The difference is not large, though, because the increase of grain volume and the ‘dilution’ of permittivity have opposite effects on the extinction (Li, 2005). For further constraint on the porosity evolution, we need a comprehensive analysis of observed dust properties including the infrared SED (Dwek et al., 1997). Polarization may also help to further constrain the models. Since such a detailed comparison also needs further modelling of the mixture of various dust species and of the interstellar radiation field, we leave it for a future work.

The porosity increases the effective sizes of grains, enhancing the grain–grain collision rate. If only coagulation occurs, however, the porosity, which is large at $a\lesssim 0.1~{}\micron$, does not have a large influence on the grain size distribution, since grains quickly coagulate to achieve $a\gtrsim 0.1~{}\micron$, and are affected by compaction. Note that, if we individually observe clouds denser than $n_{\mathrm{H}}\sim 10^{3}$ cm^-3, we could find very fluffy grains as will be shown by L. Pagani et al. (in preparation) (see also Hirashita & Li, 2013; Wong et al., 2016). The major part of interstellar grains are smaller than 1 $\micron$ (e.g. Weingartner & Draine, 2001b), so that the grains in such dense clouds are not likely to contribute directly to the interstellar dust population.

Coagulation can be balanced by shattering. For $\eta_{\mathrm{cold}}=0.5$, the grain size distribution and the filling factor both converge to equilibrium distributions at $t\sim 100$ Myr, while coagulation continuously increases the grain radii for $\eta_{\mathrm{cold}}=0.9$ (Section 3.3; Fig. 4). Thus, cases with high $\eta_{\mathrm{cold}}$ produce similar results to the pure coagulation case. For the case of $\eta_{\mathrm{cold}}=0.5$, we performed a calculation by forcing $\phi_{m}$ to be always unity (not shown) and found that the grain size distribution changes little by the different treatment of $\phi_{m}$ (i.e. compared with the results shown in Fig. 4a). This is because the evolution of grain size distribution is regulated by the formation of large grains (recall that the grain abundance is dominated by large grains), which have small porosity because of compaction.

Porosity plays a critical role in the case of small $\eta_{\mathrm{dense}}=0.1$. As discussed in Section 3.3 (Fig. 5), the grain size distribution is first dominated by shattering, which decreases the abundance of large grains, while coagulation gradually recovers the large-grain abundance at later stages (typically after $t\sim 100$ Myr). We argued above that coagulation is activated at $t\sim 100$ Myr because the increased porosity enhances the grain cross-sections (the grain–grain collision rate). For demonstration, we compare two calculations in Fig. 8: one is the same as above for $\eta_{\mathrm{dense}}=0.1$ (Fig. 4b), and the other is a calculation with the same setting but with $\phi_{m}=1$ (i.e. without porosity). Since the difference is not prominent on a short time-scale, we focus on the evolution after $t=100$ Myr. If we fix $\phi_{m}=1$, grains at $a\gtrsim 0.1~{}\micron$ are simply depleted by shattering as expected from weak coagulation in $\eta_{\mathrm{dense}}=0.1$. In contrast, if we take the evolution of $\phi_{m}$ into account, grains at $a\gtrsim 0.1~{}\micron$ are re-formed at $t\gtrsim 200$ Myr, reaching roughly an equilibrium at $t\sim 400$ Myr. The filling factor reaches the smallest value ($\phi_{m}\sim 0.1$) at $a\sim 0.1~{}\micron$. There are two effects that promote coagulation at later stages: (i) The increase of grain cross-sections by porosity enhances the grain–grain collision rate because the collision kernel scales with the porosity as $\propto\phi_{m}^{-1/3}$ (Section 2.2). (ii) The decreased grain velocity ($\propto\phi_{m}^{1/3}$; equation 6) makes compaction in coagulation at large grain radii less effective, so that the filling factor is kept small even at $a\sim 0.1~{}\micron$. These two effects are persistent qualitatively even if we consider the compaction of shattered remnants (see Fig. 5).

Recall that our treatment of the two-phase ISM is based on the parameter $\eta_{\mathrm{dense}}$, which sets the fraction of time dust spends in the dense phase. Thus, the above results imply that the dust evolution in a condition where both coagulation and shattering coexist (in a wide area of the ISM and/or on a time-scale longer than the phase exchange time-scale $\sim 10^{7}$ yr; Section 2.5) could be strongly affected by porosity. In other words, dust evolution models which do not include porosity evolution could predict a very different evolution of the grain size distribution from those which correctly take the porosity evolution into account.

For a long-term evolution of dust in the ISM, dust enrichment by stellar ejecta, dust growth by the accretion of gas-phase metals, and dust destruction by supernova shocks are also important (e.g. Dwek, 1998). Thus, it is interesting to additionally model the porosity evolution in these processes. This paper provides a first important step for the understanding of porosity evolution since coagulation and shattering are mechanisms of efficiently modifying the grain size distribution. Indeed, Hirashita & Aoyama (2019) emphasized the importance of these two processes in realizing MRN-like grain size distributions (see also Aoyama et al., 2020).

As shown in Section 4, porosity affects the UV–optical extinction curves by 10–20 per cent. As noted by Voshchinnikov et al. (2006), porosity does not necessarily increase the extinction (see also Jones, 1988; Li, 2005). At long-optical and near infrared wavelengths, the opacity of silicate can decrease owing to the porosity. The wavelength range where this decrease occurs shifts towards shorter wavelengths as the typical grain radius becomes smaller in e.g. a shattering-dominated condition. The extinction of amC is enhanced by 10–20 per cent at all wavelengths. The above results indicate that we could save 10–20 per cent of dust to explain the opacity, especially at long (far-infrared) and short (UV) wavelengths.

The effects of porosity on extinction curves are further pronounced given that the evolution of grain size distribution is affected by the filling factor. As shown above, porosity has the largest effect on the grain size distribution if strong shattering and weak coagulation are both present such as in the case of $\eta_{\mathrm{dense}}=0.1$. Thus, the extinction curve is affected by porosity not only through the optical properties but also through the modification of grain size distribution. To present this effect, we show in Fig. 9 the extinction curves corresponding to the cases shown in Fig. 8. On such a long time-scale as shown in Fig. 8, the loss of dust through the lower grain radius boundary by shattering (recall that we remove the grains which become smaller than $a=3$ Å) leads to a decrease of $A_{\lambda}/N_{\mathrm{H}}$ just because of the boundary condition. Thus, we only show $A_{\lambda}/A_{V}$, which is free from the effect of dust mass loss (i.e. we could purely observe the effect of grain size distribution). For the case without porosity evolution ($\phi_{m}=1$), the extinction curve continues to be steepened because shattering continues to convert large grains to small grains. If the porosity evolution is included, the extinction curve is flattened with age because large grains are ‘recreated’ by enhanced coagulation at later times. Both species show steeper extinction curves than the initial state for the case with porosity evolution although the grain size distributions are similar to the initial condition (Fig. 8). This is because the grain opacity does not increase in the $V$ band, where the extinction curve is normalized, while porosity enhances the UV extinction by 30–40 per cent. Jones (1988) and Voshchinnikov et al. (2006) also showed that porosity does not change the optical extinction much but increase the UV and infrared opacities. This effect makes the UV extinction curve normalized to the $V$-band value steep.

The case shown in Figs. 8 and 9 indicates that porosity evolution can have a dramatic impact on the shape of extinction curves. Even if the grain size distribution becomes the one similar to the initial grain size distribution at later stages of evolution, the small porosity makes the extinction curves significantly steeper than the initial one. In particular, the grains contributing to the optical extinction have radii $a\sim\lambda/(2\pi)\sim 0.1~{}\micron$, where the porosity is the largest. Therefore, if we consider the creation of porosity through the interplay between coagulation and shattering, the evolution of extinction curve could be qualitatively very different. Moreover, the porosity depends on the grain radius; this dependence also creates ‘higher-order’ wavelength dependence of extinction curve. As shown above, the wavelength range where the extinction is reduced by porosity shifts to longer wavelengths if porosity is developed in larger grains.

As shown in Section 3.3 (Fig. 5), if we consider compaction of shattered remnants, the filling factor increases a little at $t\gtrsim 300$ Myr. We also calculated the evolution of extinction curve with remnant compaction (now shown). Since coagulation is less efficient in this case, the extinction curves are steeper than those shown in Fig.  9. If we only see the extinction curve shape in the UV–optical, a larger filling factor and a lower abundance of large grains are degenerate.