Evolution of dust porosity through coagulation and shattering in the interstellar medium
Hiroyuki Hirashita, Vladimir B. Il'in, Laurent Pagani, Charlène Lefèvre
aa r X i v : . [ a s t r o - ph . GA ] J a n MNRAS , 1–17 (2020) Preprint 8 January 2021 Compiled using MNRAS L A TEX style file v3.0
Evolution of dust porosity through coagulation and shattering in theinterstellar medium
Hiroyuki Hirashita, ★ Vladimir B. Il’in, , , Laurent Pagani and Charlène Lefèvre , Institute of Astronomy and Astrophysics, Academia Sinica, Astronomy-Mathematics Building, No. 1, Section 4,Roosevelt Road, Taipei 10617, Taiwan St Petersburg State University, Universitetskij Pr. 28, St Petersburg 198504, Russia Pulkovo Observatory, Pulkovskoe Sh. 65/1, St Petersburg 196140, Russia St Petersburg University of Aerospace Instrumentation, Bol. Morskaya 67, St Petersburg 190000, Russia LERMA & UMR8112 du CNRS, Observatoire de Paris, PSL University, Sorbonne Universités, CNRS, F-75014 Paris, France Institut de Radioastronomie Millimétrique (IRAM), 300 rue de la Piscine, 38400 Saint-Martin d’Hères, France
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
The properties of interstellar grains, such as grain size distribution and grain porosity, are affected by interstellar processing, inparticular, coagulation and shattering, which take place in the dense and diffuse interstellar medium (ISM), respectively. In thispaper, we formulate and calculate the evolution of grain size distribution and grain porosity through shattering and coagulation.For coagulation, we treat the grain evolution depending on the collision energy. Shattering is treated as a mechanism of formingsmall compact fragments. The balance between these processes are determined by the dense-gas mass fraction 𝜂 dense , whichdetermines the time fraction of coagulation relative to shattering. We find that the interplay between shattering supplying smallgrains and coagulation forming porous grains from shattered grains is fundamentally important in creating and maintainingporosity. The porosity rises to 0.7–0.9 (or the filling factor 0.3–0.1) around grain radii 𝑎 ∼ . 𝜇 m. We also find that, in thecase of 𝜂 dense = . Key words: methods: numerical – dust, extinction – ISM: evolution – galaxies: ISM – scattering – ultraviolet: ISM
Dust grains are important for various processes in the interstellarmedium (ISM) and in galaxies. Not only the total dust abundancebut also the distribution function of grain radius, referred to as thegrain size distribution, is important, since under a given dust abun-dance, the total surface area and the cross-section of light absorptionand scattering (extinction) are governed by the grain size distribu-tion. Since H forms on dust surfaces (e.g. Gould & Salpeter 1963;Cazaux & Tielens 2004), the total surface area of dust grains affectsthe H formation rate, which could influence the galaxy evolution(e.g. Yamasawa et al. 2011; Chen et al. 2018). The extinction cross-section as a function of wavelength, referred to as the extinctioncurve, gives us a clue to the grain size distribution (e.g. Mathis et al.1977, hereafter MRN). The extinction curve is strongly modifiedby dust evolution (e.g. Asano et al. 2014). The evolution of grainsize distribution also affects the spectral energy distributions (SEDs)of galaxies (e.g. Désert et al. 1990; Takeuchi et al. 2005). Thus, theevolution of grain size distribution is of fundamental importance in ★ E-mail: [email protected] understanding the chemical and radiative processes in the ISM andin galaxies.The evolution of dust in the ISM or in galaxies is governed byvarious processes (e.g. Lisenfeld & Ferrara 1998; Draine 2009). Themajor processes for the dust enrichment are stellar dust productionand dust growth in the ISM (e.g. Dwek 1998). Dust destructionin supernova shocks (e.g. McKee 1989) is considered to dominatethe decrease of dust in the ISM. In the Milky Way environment,since the metallicity is high enough, dust growth by the accretionof gas-phase metals can be efficient enough to balance the dust de-struction (e.g. Draine 1990; Dwek 1998; Hirashita 1999; Inoue 2011;Mattsson et al. 2014). This accretion process, which is also verifiedin some experiments (Rouillé et al. 2014, 2020), occurs in the denseISM. The grain size distribution, on the other hand, is strongly mod-ified by coagulation in the dense ISM and shattering in the diffuseISM (e.g. O’Donnell & Mathis 1997; Hirashita & Yan 2009). Al-though some of the above processes occur either in a certain phaseof the ISM, all dust processes eventually affect the mean propertiesof interstellar dust through a quick ( ∼ yr) exchange between var-ious ISM phases (e.g. McKee 1989). The abundance ratio betweenlarge and small grains is roughly determined by the balance betweencoagulation and shattering (Hirashita & Aoyama 2019). © H. Hirashita
If coagulation and shattering govern the functional shape of grainsize distribution as suggested by previous studies, these processescould put a clear imprint on the grain properties. In particular, coag-ulation could form inhomogeneous grains in terms of the compositionand the shape. Such a possibility is modelled by Mathis & Whiffen(1989), who constructed a model of interstellar grains that are as-sumed to be aggregates of multiple species containing vacuum. Theysucceeded in obtaining a model that is consistent with the observedMilky Way extinction curves. Although composite aggregates (com-posed of multiple grain species) could be disrupted into single speciesin the diffuse ISM (Li & Greenberg 1997; Hoang 2019), it is atleast expected that coagulation develops non-sphericity with vac-uum, which could lead to porous (or fluffy) grains.Indeed, fluffy large grains are observationally indicated by theX-ray scattering halos surrounding point sources, such as X-ray bi-naries (Woo et al. 1994) and Nova Cygni 1992 (Mathis et al. 1995).Mathis et al. (1995) attempted to reproduce the X-ray halo strength ofNova Cygni 1992 under the observed optical extinction towards thatobject. Smith & Dwek (1998) reconsidered the treatment of grainoptical properties using the Mie solution, and concluded, contrary tothe above conclusion, that compact silicate and carbon grains are con-sistent with the observation of the X-ray scattering halo around NovaCygni 1992. Draine & Tan (2003) also indicated that a compact-grainmodel based on Weingartner & Draine (2001b) is consistent with thedust towards Nova Cygni 1992 if a lower optical extinction is adopted.Draine (2003) also discussed the X-ray halo of Cyg X-1 in additionto that of Nova Cygni 1992 and supported the same compact-grainmodel; however, they also mentioned the uncertainty arising fromthe distribution of dust in the line of sight. Within the uncertainty,grains with moderate porosities could still exist in the diffuse ISMConsidering porous grains was also motivated by the severe con-straint on the available metals (especially carbon) historically. Mathis(1996) investigated a possibility of reproducing the Milky Way ex-tinction curve with fluffy grains. As a consequence, they ‘saved’ theheavy elements the most (i.e. used the minimum amount of metals fordust) by including 25–65 per cent of vacuum in the dust grains. Veryfluffy grains are excluded because their extinction (absorption + scat-tering) per unit dust mass is too low. However, Dwek et al. (1997) ar-gued that, if we take the fitting to the infrared dust emission SED intoaccount, Mathis (1996)’s model overproduces the optical–ultraviolet(UV) extinction. Li (2005) showed, based on the argument on theKramers-Kronig relations, that the observed optical–UV extinctionis still underproduced with the Galactic metal abundance derivedfrom B stars. We should also keep in mind that recent dust modelswith updated extinction-to-column density ratios ( 𝐴 𝑉 / 𝑁 H ) and theprotosolar abundance with a correction for Galactic chemical evolu-tion after the Sun formation indicate that the above severe metallicityconstraint may not be an issue any more (Draine & Hensley 2020;Hensley & Draine 2020; Zuo et al. 2020).Clarifying how porous the interstellar dust could be is in itself aproblem independent of the metal abundance constraint. Since theabove mentioned studies did not necessarily give a definite con-straint on the grain porosity, it is still meaningful to investigatethe effect of porosity on the extinction curves in a general context.Voshchinnikov et al. (2005, 2006) studied the dependence of extinc-tion curve on the porosity. These studies give a basis on which weconstrain the allowed ranges of porosity through comparison withobservations. In reality, porosity could depend on the grain size, be-cause the evolution of grain size distribution plays an important rolein determining the grain porosity. However, it is difficult to relategrain size and porosity, since the physical link between these twoquantities is still missing. Because of this ‘missing link’, the poros- ity still remains to be a free parameter. It is definitely necessary tounderstand the interplay between grain size and porosity in the evo-lution of interstellar dust, if we would like to obtain a meaningfulinsight into the physical mechanisms of dust evolution through theobservationally constrained grain porosities.The evolution of porous grains by coagulation has been modelledand discussed in the field of protoplanetary discs. Ormel et al. (2007)developed a Monte Carlo method for the evolution of grain size dis-tribution by coagulation in protoplanetary discs taking the porosityevolution into account. In each collision, depending on the collisionenergy, they treated the difference between the hit-and-stick (stick-ing without modifying the grain shapes) and compaction (decreaseof grain volume by compression) regimes. Okuzumi et al. (2009,hereafter O09) formulated this problem by using the 2-dimensional(2D) Smoluchowski equation that treats distribution function of twovariables, grain radius and porosity. Since directly solving the 2Dequation is computationally expensive, they only concentrated onthe mean porosity at each grain radius while they fully solved thegrain size distribution (see also Okuzumi et al. 2012, hereafter O12).They finally obtained a set of moment equations for the grain vol-ume: the zeroth order describes the grain size distribution and thefirst order the ‘grain volume distribution’. A set of grain radius andgrain volume gives information on the grain porosity. The calculatedporosity is also useful to predict the radiative properties of dust grainsin protoplenetary discs (e.g. Kataoka et al. 2014; Tazaki et al. 2016).The above methods to trace the evolution of porous grains are ap-plied to the grain evolution in dense molecular clouds. Ormel et al.(2009) applied the above mentioned Monte Carlo method (but in-cluding also grain fragmentation) and showed that the filling factorof grain could become as small as ∼ . ∼ . & cm − ) environments, and therelation to the grains actually contributing to the interstellar extinc-tion is not clear. As mentioned above, shattering is also importantas a modifying mechanism of grain size distribution. The functionalshape of the grain size distribution could be determined by the bal-ance between shattering and coagulation. This balance could also beimportant to determine the porosities at various grain sizes.The goal of this paper is to clarify the evolution of porous grainsin the ISM by considering both coagulation and shattering. We havealready formulated the evolution of grain size distribution using theSmoluchowski equation in our previous works (e.g. Hirashita & Yan2009; Asano et al. 2013; Hirashita & Aoyama 2019). For the poros-ity evolution, we utilize the framework developed by O09 and O12 forprotoplanetary discs with some modifications mainly to include shat-tering. Our model developed in this paper is aimed at being used tocalculate the dust evolution on galaxy scales in the future. Althoughthe Smoluchowski equation usually refers to the one describing co-agulation, we broadly use this name also for an extended version forshattering.As a result of the modelling in this paper, we will be able topredict the expected grain porosity as a function of grain radius inthe ISM. We emphasize that the previous models of grain porosityintroduced above mostly treated porosity as a free parameter to befitted to observed dust properties such as extinction curves. Thus, thisstudy will provide the first prediction on the grain-radius-dependentporosity expected theoretically from coagulation and shattering. Wealso calculate extinction curves to clarify the effect of predictedporosity on an observable dust property. This paper is not aimed MNRAS , 1–17 (2020) orosity of interstellar dust at complete modelling but focused on the construction of a basicframework for clarifying how large the porosity could be in thepresence of coagulation and shattering.This paper is organized as follows. In Section 2, we formulate theevolution of grain size distribution and porosity. We show the resultsin Section 3, and present the effects on extinction curves in Section 4.We provide some extended discussions, especially focusing on theeffects of porosity in Section 5. We give the conclusion of this paperin Section 6. First of all, we clarify some terms. The filling factor of a grain isdefined as the volume fraction occupied by the material composingthe grain. The rest (unity minus the filling factor) is referred to asthe porosity , which is the volume fraction of vacuum. A grain is compact if the filling factor is unity (i.e. no porosity). In our previouspapers (e.g. Hirashita & Aoyama 2019) and others (e.g. Jones et al.1994, 1996; Mattsson 2016), the evolution of grain size distributionin the ISM is formulated for compact grains. Here, we extend theformulation to include the evolution of grain porosity.As mentioned in the Introduction, we concentrate on the colli-sional processes for grain evolution, namely coagulation and shatter-ing, since they are the main drivers for the evolution of both grain sizedistribution and grain porosity. We neglect dust evolution processesother than coagulation and shattering. This treatment assumes thatcoagulation and/or shattering occur in a ‘closed box’ without anysupply from stellar dust production and from accretion of gas-phasemetals, and without any loss of dust by supernova shock destruction.For simplicity, we consider a single dust species in solving the evo-lution of grain size distribution and porosity to avoid the complexityarising from compound species. That is, we concentrate on the effectof porosity in a single dust species.Coagulation and shattering occur in the different ISM phases(dense and diffuse ISM phases, respectively). However, since theexchange of these phases occurs on a short time-scale ( ∼ yr;McKee 1989), we consider the mixture of coagulation and shatteringin addition to treating each process separately. For the dense anddiffuse ISM phases, we consider representative phases that couldaffect the major dust populations in the ISM and the physical con-ditions are given in Section 2.5. The grain motions are assumed tobe driven by interstellar turbulence. We do not trace coagulation in astar-forming cloud with densities & cm − , since, as mentioned inthe Introduction, the relation of the grains formed in such extreme en-vironments to the general ISM grains is not clear. Moreover, we needto consider physical processes such as ambipolar diffusion in con-sidering the grain motion in star-forming clouds (e.g. Silsbee et al.2020; Guillet et al. 2020). Therefore, we leave the detailed studies ofdense star-forming clouds for a future work.Coagulation in the dense ISM could be complex sinceaccretion of gas-phase materials could occur at the sametime (Hirashita & Voshchinnikov 2014; Voshchinnikov & Hirashita2014). Accretion could decrease the porosity because part of thevacuum could be filled (Li & Lunine 2003, appendix B). As a re-sult of accretion in an environment where UV radiation is shielded,ice mantles could also develop on the grain surfaces. However, suchmantles are evaporated in the diffuse ISM so that aggregates attachedthrough the ice mantles could dissolve once they are injected into thediffuse ISM. As mentioned later, whether or not ice mantles enhancecoagulation is still debated, and refractory surfaces could be moresticky (Section 2.3.1). Therefore, it is hard to clarify with the cur- rently available knowledge how ice mantles affect the observed grainproperties (e.g. extinction curves) in the ISM through coagulation.Since we are interested in the general population of grains in the ISM(especially, silicate and carbonaceous dust), we neglect the formationof volatile mantles, and assume that the collisions between refractorygrains in the dense ISM lead to formation of refractory coagulatedgrains. We also neglect the possible development of core-mantlestructures as a result of accretion, which could explain the MilkyWay dust properties (Li & Greenberg 1997; Jones et al. 2017). Wesimply examine the effect of porosity on the extinction curves of basicmaterials (silicate and carbonaceous dust in this paper; Section 2.4),but do not aim at reproducing the observed Milky Way extinctioncurve.In this paper, we also neglect the electric charge of grains. The graincollision cross-section could be enhanced or reduced depending onthe grain charge. The charge effect is not important for large grainsbecause they have high kinetic energy (i.e. the Coulomb force doesnot alter the grain motion) as commented in Hirashita & Aoyama(2019). The charge may be important for coagulation of smallgrains, which are positively or negatively charged depending onthe physical condition of the ambient medium and radiation field(Weingartner & Draine 2001a). For example, in a physical conditionappropriate for molecular clouds, small grains tend to be chargednegatively, while large grains positively (Yan et al. 2004). This situ-ation could suppress the collisions between small grains but enhancethose between a small and a large grain. At the same time, the graincharge is very sensitive to the assumed physical conditions in thedense clouds (especially, the ionization degree and the UV radiationfield intensity, which is affected by the dust attenuation). Although itis interesting to investigate the detailed dependence on these physi-cal parameters, in this paper we choose to neglect the effect of graincharge for the simplicity of our treatment. The evolution of grain size distribution and grain porosity by coag-ulation is formulated by O09 and O12. Although their main targetwas dust evolution in protoplanetary discs, the generality of theirformulation allows us to apply it to interstellar dust. We also includeshattering in this paper. The evolution of grain size distribution andgrain porosity by coagulation and shattering is described by the dis-tribution functions of grain mass ( 𝑚 ) and grain volume ( 𝑉 ) at time 𝑡 , and is governed by the 2D Smoluchowski equation. However, asnoted by O09, solving the 2D Smolchowski equation requires a hugecomputational expense. Therefore, we concentrate on the momentequations for 𝑉 , but we still solve the full distribution function for 𝑚 .To treat porous grains, it is convenient to introduce the follow-ing two types of grain radius: characteristic radius ( 𝑎 ch ) and mass-equivalent radius ( 𝑎 𝑚 ). The characteristic radius is related to thevolume as 𝑉 = ( 𝜋 / ) 𝑎 , while the mass-equivalent radius is re-lated to the grain mass as 𝑚 = ( 𝜋 / ) 𝑎 𝑚 𝑠 , where 𝑠 is the bulkmaterial density. For compact grains, 𝑚 = 𝑠𝑉 , so that 𝑎 ch = 𝑎 𝑚 . Forporous grains, 𝑎 ch > 𝑎 𝑚 . Using these two grain radii, the filling fac-tor, 𝜙 𝑚 , is expressed as 𝜙 𝑚 ≡ ( 𝑎 𝑚 / 𝑎 ch ) (≤ ) (see also Ormel et al.2009), while the porosity is 1 − 𝜙 𝑚 . Since it is necessary to specify 𝑠 ,we simply adopt the same value as used in HM20 ( 𝑠 = .
24 g cm − taken from graphite), but adopting 𝑠 = − appropriatefor interstellar grains (e.g. Weingartner & Draine 2001b) does not We use the upper case 𝑉 for volumes and the lower case 𝑣 for velocitiesthroughout this paper. MNRAS , 1–17 (2020) H. Hirashita change our conclusion significantly. We vary some parameters rel-evant for porosity, which are more important in this paper (Section2.5). Thus, we simply use the same parameters as in HM20 for thematerial properties (specifically 𝑠 and 𝑄 ★ D introduced later) that donot directly affect the porosity evolution.Before presenting the basic equations, we need to introduce somequantities. We use the grain mass distribution instead of the grainsize distribution since there is ambiguity in the grain size for porousgrains. We denote the distribution function of 𝑚 and 𝑉 at time 𝑡 as 𝑓 ( 𝑚, 𝑉, 𝑡 ) . The grain mass distribution at time 𝑡 , ˜ 𝑛 ( 𝑚, 𝑡 ) is definedsuch that ˜ 𝑛 ( 𝑚, 𝑡 ) d 𝑚 is the number density of grains whose massis between 𝑚 and 𝑚 + d 𝑚 , and is related to the above distributionfunction as˜ 𝑛 ( 𝑚, 𝑡 ) ≡ ∫ ∞ 𝑓 ( 𝑚, 𝑉, 𝑡 ) d 𝑉 . (1)We also introduce the first moment of 𝑓 for 𝑉 as¯ 𝑉 ( 𝑚, 𝑡 ) ≡ 𝑛 ( 𝑚, 𝑡 ) ∫ ∞ 𝑉 𝑓 ( 𝑚, 𝑉, 𝑡 ) d 𝑉 . (2)This quantity ( ¯ 𝑉 ) is the mean volume of grains with mass 𝑚 .We define the following two quantities, 𝜌 ( 𝑚, 𝑡 ) ≡ 𝑚 ˜ 𝑛 ( 𝑚, 𝑡 ) and 𝜓 ( 𝑚, 𝑡 ) ≡ ¯ 𝑉 ( 𝑚, 𝑡 ) ˜ 𝑛 ( 𝑚, 𝑡 ) , which represent the grain mass distri-bution weighted for the grain mass and volume, respectively (thelatter can also be regarded as the volume distribution function). Theintegration of 𝜌 ( 𝑚, 𝑡 ) for 𝑚 is related to the dust-to-gas ratio, D , as 𝜇 H 𝑚 H 𝑛 H D = ∫ ∞ 𝜌 ( 𝑚, 𝑡 ) d 𝑚, (3)where 𝜇 H = . 𝑚 H is the hydrogenatom mass, and 𝑛 H is the hydrogen number density.The Smolchowski equation contains the collision kernel, which isthe product of the collision cross-section and the relative speed ofcolliding pair. The collision kernel for colliding grains with masses 𝑚 and 𝑚 is denoted as 𝐾 𝑚 ,𝑚 , which is evaluated in Section 2.2.It is straightforward to derive the following general equations ap-plicable for both coagulation and shattering based on equations (6)and (7) of O09: 𝜕𝜌 ( 𝑚, 𝑡 ) 𝜕𝑡 = − 𝑚𝜌 ( 𝑚, 𝑡 ) ∫ ∞ 𝐾 𝑚,𝑚 𝑚𝑚 𝜌 ( 𝑚 , 𝑡 ) d 𝑚 + ∫ ∞ ∫ ∞ 𝐾 𝑚 ,𝑚 𝑚 𝑚 𝜌 ( 𝑚 , 𝑡 ) 𝜌 ( 𝑚 , 𝑡 ) 𝑚 ¯ 𝜃 ( 𝑚 ; 𝑚 , 𝑚 ) d 𝑚 d 𝑚 , (4) 𝜕𝜓 ( 𝑚, 𝑡 ) 𝜕𝑡 = − ¯ 𝑉 ( 𝑚, 𝑡 ) 𝜓 ( 𝑚, 𝑡 ) ∫ ∞ 𝐾 𝑚,𝑚 ¯ 𝑉 ( 𝑚, 𝑡 ) ¯ 𝑉 ( 𝑚 , 𝑡 ) 𝜓 ( 𝑚 , 𝑡 ) d 𝑚 + ∫ ∞ ∫ ∞ 𝐾 𝑚 ,𝑚 ¯ 𝑉 ( 𝑚 ) ¯ 𝑉 ( 𝑚 ) 𝜓 ( 𝑚 , 𝑡 ) 𝜓 ( 𝑚 , 𝑡 )( 𝑉 + ) 𝑚𝑚 ,𝑚 ¯ 𝜃 ( 𝑚 ; 𝑚 , 𝑚 ) d 𝑚 d 𝑚 , (5)where ˜ 𝜃 ( 𝑚 ; 𝑚 , 𝑚 ) describes the distribution function of grainsproduced from 𝑚 in the collision between grains with masses 𝑚 and 𝑚 , and ( 𝑉 + ) 𝑚𝑚 ,𝑚 is the volume of the newly produced grainwith mass 𝑚 in the above collision. Since the expression is not exactlythe same as that used by O09 and O12, we explain the derivation of theabove equations in Appendix A. In particular, in our formulation, wecount the collisional products originating from 𝑚 and 𝑚 separately.O09, in contrast, considered the collisional pair 𝑚 and 𝑚 only once.This is why O09 has a factor 1/2 in front of the second term on theright-hand side in equations (4) and (5).We aim at solving equations (4) and (5) for 𝜌 and 𝜓 . Since 𝜌 / 𝑚 = 𝜓 / ¯ 𝑉 is always satisfied, ¯ 𝑉 is not an independent quantity. Thus, these two equations are closed if we give 𝐾 , 𝑉 + and ¯ 𝜃 . In the followingsubsections, we explain how to determine these three quantities.In computing the grain size distribution, we discretize the entiregrain radius range ( 𝑎 = × − –10 𝜇 m) into 128 grid points withlogarithmically equal spacing. We set 𝜌 d ( 𝑚, 𝑡 ) = The collision kernel is determined by the product of the grain cross-section and the relative velocity. The cross-section ( 𝜎 , ) in the col-lision of grains with masses 𝑚 and 𝑚 , referred to as grain 1 and 2,respectively, is estimated as 𝜎 , = 𝜋 ( 𝑎 ch1 + 𝑎 ch2 ) , where the char-acteristic radii of grains 1 and 2 are 𝑎 ch1 and 𝑎 ch2 , respectively (seethe beginning of Section 2.1 for the definition of the characteristicradius; i.e. 𝑎 ch = 𝑎 𝑚 𝜙 − / 𝑚 ).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 thesame simple model, originally taken from Ormel et al. (2009), forthe grain velocity 𝑣 gr as adopted in our previous papers: 𝑣 gr ( 𝑚 ) = . M / 𝜙 / 𝑚 (cid:18) 𝑎 𝑚 . 𝜇 m (cid:19) / (cid:18) 𝑇 gas K (cid:19) / (cid:18) 𝑛 H − (cid:19) − / × (cid:18) 𝑠 . − (cid:19) / km s − , (6)where M is the Mach number of the largest-eddy velocity, and 𝑇 gas is the gas temperature. This is similar to equation (18) ofHirashita & Aoyama (2019), but modified to include the filling fac-tor 𝜙 𝑚 (see appendix A of Ormel et al. 2009). The derivation of theabove equation does not strictly hold for highly supersonic regime,but we practically use M here as an adjusting parameter for the grainvelocity. Under fixed ISM conditions and grain material properties,the velocity is reduced to a function of 𝑚 . Since the filling factor 𝜙 𝑚 changes as a function of time, 𝑣 gr ( 𝑚 ) also varies with time. In con-sidering the collision rate between grains 1 and 2 with 𝑣 gr ( 𝑚 ) = 𝑣 and 𝑣 gr ( 𝑚 ) = 𝑣 , we estimate the relative velocity 𝑣 , by 𝑣 , = q 𝑣 + 𝑣 − 𝑣 𝑣 𝜇 , , (7)where 𝜇 = cos 𝜃 ( 𝜃 is an angle between the two grain velocities)is randomly chosen between − 𝐾 𝑚 ,𝑚 = 𝜎 , 𝑣 , . (8)Note that the collision kernel depends on the filling factor of the largergrain as ∝ 𝜙 − / 𝑚 . Thus, the grain–grain collision rate is enhanced ifthe grains become more porous (fluffy). The treatment of collisional products is the key for the evolutionof porosity. The porosity is strongly affected by the production ofvacuum volume in a grain after grain–grain sticking (coagulation)and the compaction associated with compression in grain–grain colli-sions. However, these processes are strongly nonlinear and dependenton the actual grain shapes. Moreover, some assumptions in O09 and
MNRAS , 1–17 (2020) orosity of interstellar dust O12 are not applicable to interstellar dust. In particular, their for-mulation is based on single-sized ( ∼ . 𝜇 m) monomers, since theyare interested in protoplanetary discs, where grains grow into muchlarger sizes than interstellar grains. In the ISM, in contrast, we alsoconsider the evolution of grains much smaller than the ‘typical’ grainsize ( ∼ . 𝜇 m). To make the problem more complicated, shatteringalso occurs, producing a large number of small grains. Because ofthe above difficulty and the difference from O09 and O12, we takethe following approach that simplifies the treatment of porosity evo-lution but still preserves the essence of physical processes regardinggrain–grain sticking and compression. For coagulation, two characteristic energies are important: impact(kinetic) energy and the rolling energy. The impact energy ( 𝐸 imp ) inthe collision between grains 1 and 2 is estimated as 𝐸 imp = 𝑚 𝑚 𝑚 + 𝑚 𝑣 , . (9)The rolling energy ( 𝐸 roll ) is defined as the energy necessary for amonomer to roll over 90 degrees on the surface of another monomer(Dominik & Tielens 1997), and is given by (Wada et al. 2007) 𝐸 roll = 𝜋 𝛾𝑅 , 𝜉 crit , (10)where 𝛾 is the surface energy per unit contact area, 𝑅 , is the re-duced particle radius, and 𝜉 crit is the critical displacement of rolling.The surface energy depends on the surface material: 𝛾 =
25, 75,and 100 erg cm − for silicate (originally from quartz), graphite,and water ice, respectively (Dominik & Tielens 1997; Israelachvili1992; Wada et al. 2007). The critical displacement 𝜉 crit lies broadlyin the range of ∼ 𝛾 and 𝜉 crit are degenerate, we fix 𝛾 to an intermediate value (75 erg cm − ) and vary 𝜉 crit . We choosea fiducial value 𝜉 crit =
10 Å unless otherwise stated but we laterexamine a case where 𝜉 crit is varied.The reduced radius, 𝑅 , , is difficult to determine rigidly for inter-stellar dust, since the typical monomer size is not obvious. Thus, wesimply assume that1 𝑅 , = 𝑎 𝑚 + 𝑎 𝑚 . (11)This equation gives a reasonable estimate for 𝑅 , for compact grains.As shown later, small grains tend to be compact. Since the reducedradius is dominated by the smaller grain, the above estimate couldbe justified. Large (submicron/micron) grains are also affected bycompaction, so that the above reduced radius could be applied forcollisions between large grains. However, the above reduced massmay be overestimated in collisions between porous grains. Grainswith intermediate radii ( 𝑎 ∼ . 𝜇 m) tend to be porous. Over-estimation of 𝑅 , leads to less efficient compaction. On the otherhand, the uncertainty in 𝑅 , could be absorbed by the variation of 𝜉 crit and 𝑛 c ( 𝑛 c is introduced below). Thus, we vary 𝜉 crit and 𝑛 c toeffectively investigate the uncertainties caused by collisions betweenintermediate-sized grains.We now estimate the volume of the coagulated grain. We donot directly adopt O12’s formulation (their equation 15; see alsoSuyama et al. 2012) because, as mentioned above, the monomer sizeis not well determined for interstellar dust. Nevertheless, we apply the following essential points (Dominik & Tielens 1997): (i) The phys-ical outcome of the collision is characterized by the ratio between 𝐸 imp and 𝐸 roll . (ii) If 𝐸 imp ≪ 𝐸 roll , the grains stick without modify-ing their structures (hit-and-stick collisions). (iii) If 𝐸 imp ∼ 𝐸 roll , thenewly added void tends to be compressed because the contact pointsof monomers are moved. (iii) If 𝐸 imp & 𝑛 c 𝐸 roll , where 𝑛 c is thenumber of contact points, significant compaction occurs. Since wedo not trace each monomer, 𝑛 c is uncertain. Moreover, not all contactpoints are equally important in compaction. Thus, we treat 𝑛 c as afree parameter, and regard it as the number of major contacts whosemovement contributes to compaction significantly. Since 𝑛 c alwaysappears in the form of product 𝑛 c 𝐸 roll , 𝑛 c is degenerate with 𝐸 roll as mentioned above. Considering the above three points, (i)–(iii), weestimate the volume ( 𝑉 + ) of coagulated product in the collisionbetween grains with volumes 𝑉 and 𝑉 (assuming that 𝑉 ≥ 𝑉 ; if 𝑉 < 𝑉 , we exchange grains 1 and 2, so that the following formula-tion holds) as 𝑉 + = 𝑉 + 𝑉 + 𝑉 void exp (cid:2) − 𝐸 imp /( 𝑏𝐸 roll ) (cid:3) + ( 𝑉 , comp − 𝑉 ) (cid:8) − exp (cid:2) − 𝐸 imp /( 𝑛 c 𝐸 roll ) (cid:3) (cid:9) , (12)where 𝑉 void is the volume of the void newly created in the collision, 𝑏 = .
15 is an adjusting parameter (Wada et al. 2008), and 𝑉 , comp ≡ 𝑚 / 𝑠 (the volume with perfect compaction). The void volume isestimated as (O12) 𝑉 void = min (cid:20) . − .
03 ln (cid:18) 𝑉 / 𝑉 + (cid:19) , . (cid:21) 𝑉 . (13)Equation (12) satisfies the above conditions (i)–(iii): 𝑉 + ≃ 𝑉 + 𝑉 + 𝑉 void for 𝐸 imp ≪ 𝐸 roll . The newly created volume ( 𝑉 void ) iscompressed if 𝐸 imp & 𝐸 roll . If 𝐸 imp becomes comparable to or higherthan 𝑛 c 𝐸 roll , the grain is compressed further (note that 𝑉 , comp − 𝑉 ≤ 𝐸 imp ≫ 𝑛 c 𝐸 roll , 𝑉 + = ( 𝑉 − 𝑉 ) + 𝑉 , comp (notethat 𝑉 ≥ 𝑉 ), which means that 𝑉 becomes compact while 𝑉 iscompressed by the equivalent volume to 𝑉 . According to Wada et al.(2013), in a collision of grains with different sizes, the larger grainis not entirely compressed. We assume that the compressed volumein the larger grain is determined by the volume of the smaller grain.In the above expression, 𝑉 + can be smaller than 𝑉 , comp + 𝑉 , comp ,where 𝑉 , comp = 𝑚 / 𝑠 (the volume with perfect compaction). Thus,we finally adopt the following expression for the coagulated volume, 𝑉 + : 𝑉 + = ( 𝑉 + if 𝑉 + > 𝑉 , comp + 𝑉 , comp , ( + 𝜖 𝑉 )( 𝑉 , comp + 𝑉 , comp ) otherwise, (14)where 𝜖 𝑉 ≥ ( 𝑉 + ) 𝑚𝑚 ,𝑚 in equation (5).The mass distribution of collisional products is described by (seeequation 26 of Hirashita & Aoyama 2019) 𝑚 ¯ 𝜃 ( 𝑚 ; 𝑚 , 𝑚 ) = 𝑚 𝛿 [ 𝑚 − ( 𝑚 + 𝑚 )] , (15)where 𝛿 (·) is Dirac’s delta function. As mentioned above and inAppendix A, we separately treat the collision product originatingfrom 𝑚 and that from 𝑚 (i.e. we count the same collision twice). There has not been any formulation for the evolution of porosityby shattering in the ISM. Shattering produces a lot of fragments,and sometimes leaves a remnant if the original grain is much larger
MNRAS000
MNRAS000 , 1–17 (2020)
H. Hirashita than the colliding partner. We expect that shattered fragments havesmall porosities because weakly bound parts become unbound indisruptive collisions. Thus, we assume that the collisional fragmentsare compact. For the remnant, it is probable compaction occurs;however, in our formulation, the remnant remains only if a graincollides with a much smaller grain: in this case, it is not clear if theentire remnant is efficiently pressed. Thus, we simply assume thatthe shattered remnant has the same porosity as the original grain. Weexamine separately a case where we take compaction of the remnantsinto account.The microscopic processes for shattering of porous grains in high-velocity ( & − ) collisions are not well known. Moreover, in ourformulation, we have to treat collisions of both porous and compactgrains. Thus, we use our previous formulation in Hirashita & Aoyama(2019), which is appropriate for compact grains. We also expect thatthe resulting grain size distribution is not sensitive to detailed as-sumptions on the fragment size distribution (Hirashita & Kobayashi2013). We explain the treatment of shattered fragments in what fol-lows.We consider a collision of two dust grains with masses 𝑚 and 𝑚 (grains 1 and 2) based on Kobayashi & Tanaka (2010)’s model. Theejected mass ( 𝑚 ej ) from grain 1 is estimated as 𝑚 ej = 𝜑 + 𝜑 𝑚 , (16)with 𝜑 ≡ 𝐸 imp 𝑚 𝑄 ★ D , (17)where 𝑄 ★ D is the specific impact energy required to disrupt half of themass ( 𝑚 / 𝑄 ★ D = . × erg g − following HM20.The ejected mass is distributed into fragments, for which we assumea power-law size distribution with an index of 𝛼 f = . 𝜃 as ¯ 𝜃 ∝ 𝑚 −( 𝛼 f + )/ .The maximum and minimum masses of the fragments are assumedto be 𝑚 f , max = . 𝑚 ej and 𝑚 f , min = − 𝑚 f , max , respectively(Guillet et al. 2011). We adopt the following mass distribution func-tion of collisional product from grain 1 in the collision with grain 2as 𝑚 ¯ 𝜃 ( 𝑚, 𝑚 , 𝑚 ) = ( − 𝛼 f ) 𝑚 ej 𝑚 (− 𝛼 f + )/ (cid:20) 𝑚 − 𝛼 f3 f , max − 𝑚 − 𝛼 f3 f , min (cid:21) Φ ( 𝑚 ; 𝑚 f , min , 𝑚 f , max )+ 𝑚𝛿 ( 𝑚 − 𝑚 + 𝑚 ej ) , (18)where Φ ( 𝑚 ; 𝑚 f , min , 𝑚 f , max ) = 𝑚 f , min ≤ 𝑚 ≤ 𝑚 f , max , and 0otherwise. Grains which become smaller than the minimum grainsize ( 𝑎 𝑚 = × − 𝜇 m) are removed. For the volumes of shatteredproducts, we adopt 𝑉 + = ( ¯ 𝑉 ( 𝑚 )/( + 𝜑 ) if 𝑚 = 𝑚 − 𝑚 ej = 𝑚 /( + 𝜑 ) , 𝑚 / 𝑠 otherwise. (19)The first case of this expression indicates that the porosity (or fillingfactor) of the remnant is the same as that of the original grain.The second case means that the fragments are compact (that is, thevolumes are simply estimated by the mass divided by the materialdensity). We use equation (19) for ( 𝑉 + ) 𝑚𝑚 ,𝑚 in equation (5). To investigate the effects on observed dust properties, we also cal-culate extinction curves. We examine two representative grain ma-terials: silicate and carbonaceous species, and investigate how the porosity evolution driven by coagulation and shattering affects the ex-tinction properties of these materials. For silicate, we use the opticalconstants 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 car-bonaceous material, and we confirmed that it produces similar resultsto amC except at wavelengths where the 2175 Å feature is prominent.We also found that the wavelength where the feature peaks shifts withthe porosity, but such a shift is not observed in the Milky Way ex-tinction curves. Thus, special care should be taken of the modellingof the 2175 Å feature, e.g. by treating the 2175 Å carrier such asgraphite and PAHs (e.g. Li & Draine 2001) as a component sepa-rated from porous grain species (Voshchinnikov et al. 2006), whichis out of the scope of this paper. Thus, we adopt amC in this paper,noting that, as mentioned above, amC and graphite produce similarresults at wavelengths not affected by the 2175 Å feature.The optical properties of dust is calculated using the effectivemedium theory (EMT), which averages the dielectric permittivityusing a mixing rule. We adopt the Bruggeman mixing rule. Thismixing rule as well as the Garnett mixing rule gives reasonableextinctions 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 grainssmaller than ∼ . 𝜇 m and we are interested in wavelengths longerthan 0.1 𝜇 m. Thus, the above condition is satisfied in this paper.Even if we model the inhomogeneity using multi-layered spheres, thedifference in the extinction is expected to be less than ∼
10 per cent(Voshchinnikov et al. 2005; Shen et al. 2008). Using the Bruggemanmixing rule, the averaged dielectric permittivity ¯ 𝜀 is obtained from ( − 𝜙 𝑚 ) − ¯ 𝜀 + 𝜀 + 𝜙 𝑚 𝜀 − ¯ 𝜀𝜀 + 𝜀 = , (20)where 𝜀 is the dielectric permittivity of the material (note that thefirst material is assumed to be vacuum so 𝜀 = We assumethat each dust grain is a sphere with radius 𝑎 ch and refractive in-dex ¯ 𝑚 = √ ¯ 𝜀 . The cross-section 𝐶 ext ,𝑚 , which is a function of 𝑚 under a given 𝜙 𝑚 at each epoch, is calculated by the Mie theory(Bohren & Huffman 1983).The extinction at wavelength 𝜆 in units of magnitude ( 𝐴 𝜆 ) can becalculated using the grain size distribution as 𝐴 𝜆 = ( . e ) 𝐿 ∫ ∞ ˜ 𝑛 ( 𝑚 ) 𝐶 ext ,𝑚 d 𝑚, (21)where 𝐿 is the path length. We present the extinction in the followingtwo ways: 𝐴 𝜆 / 𝑁 H and 𝐴 𝜆 / 𝐴 𝑉 (the 𝑉 band wavelength correspondsto 𝜆 − = . 𝜇 m − and 𝑁 H = 𝑛 H 𝐿 is the column density of hydrogennuclei). The first quantity indicates the extinction per hydrogen (thus,it is proportional to the dust abundance), while the second is usefulto show the shape of extinction curve. In both quantities, the pathlength 𝐿 is canceled out. Coagulation and shattering are governed by the same form of equa-tion (equations 4 and 5). These two processes occur in different ISMphases. We assume that coagulation occurs in the dense clouds ofwhich the physical conditions are described by 𝑛 H = cm − and 𝑇 gas =
10 K (typical of molecular clouds); and that shat-tering takes place in the diffuse ISM characterized by 𝑛 H = . We adopt Gaussian-cgs units.MNRAS , 1–17 (2020) orosity of interstellar dust cm − and 𝑇 gas = K. These values are similar to the onesadopted for coagulation and shattering by Hirashita & Yan (2009).The density affects the grain collision time-scale, which scales with ( 𝑛 H 𝑣 gr ) − ∝ 𝑛 − / (see equation 6); that is, the adopted duration forshattering/coagulation is degenerate with the gas density. Therefore,we fix the above physical conditions for the gas and examine the grainsize distributions at various times (i.e. for various durations of theprocess), keeping in mind that the same value of 𝑛 / 𝑡 gives similarresults.For the normalization of the grain velocities adjusted by M inequation (6), we apply M = M = M are adopted to obtain a similar levelof grain velocities to those calculated for magnetized turbulence byYan et al. (2004).We not only treat shattering and coagulation separately, but alsoexamine some cases where these two processes occur at the sametime. The coexistence of shattering and coagulation is a reasonableapproximation for a volume of the ISM wide enough to samplea statistically significant amount of both ISM phases and/or for atime much longer than the mixing time-scale of the two ISM phases( ∼ yr; e.g. McKee 1989). Since coagulation and shattering occurin the media with different densities, it is convenient to define thegrain mass and volume distributions per hydrogen number density as˜ 𝜌 ≡ 𝜌 / 𝑛 H and ˜ 𝜓 ≡ 𝜓 / 𝑛 H , respectively. In this case, we calculate theevolution of ˜ 𝜌 and ˜ 𝜓 by 𝜕 ˜ 𝜌 ( 𝑚, 𝑡 ) 𝜕𝑡 (cid:12)(cid:12)(cid:12)(cid:12) tot = 𝜂 dense 𝜕 ˜ 𝜌 ( 𝑚, 𝑡 ) 𝜕𝑡 (cid:12)(cid:12)(cid:12)(cid:12) coag + ( − 𝜂 dense ) 𝜕 ˜ 𝜌 ( 𝑚, 𝑡 ) 𝜕𝑡 (cid:12)(cid:12)(cid:12)(cid:12) shat , (22) 𝜕 ˜ 𝜓 ( 𝑚, 𝑡 ) 𝜕𝑡 (cid:12)(cid:12)(cid:12)(cid:12) tot = 𝜂 dense 𝜕 ˜ 𝜓 ( 𝑚, 𝑡 ) 𝜕𝑡 (cid:12)(cid:12)(cid:12)(cid:12) coag + ( − 𝜂 dense ) 𝜕 ˜ 𝜓 ( 𝑚, 𝑡 ) 𝜕𝑡 (cid:12)(cid:12)(cid:12)(cid:12) shat , (23)where 𝜂 dense , which is treated as a free parameter, is the dense gasfraction determining the weight of each process, the subscript ‘tot’indicates the total changing rates of ˜ 𝜌 and ˜ 𝜓 , 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 aboveequations, we apply equations (4) and (5). The change by coagulationand that by shattering are mixed with a ratio of 𝜂 dense : ( − 𝜂 dense ) ;in other words, 𝜂 dense determines the fraction of time dust spendsin the dense ISM. In calculating shattering and coagulation, we use 𝜌 and 𝜓 , but when we sum up the contributions from coagulationand shattering, we divide them by 𝑛 H and obtain ˜ 𝜌 and ˜ 𝜓 . Here weimplicitly neglect gas which hosts neither coagulation nor shattering.Existence of such a gas component effectively lowers the efficienciesof both shattering and coagulation (or makes the time-scales of thesetwo processes longer) equally and does not affect the relative rolesof these two processes. As mentioned in Section 2.3.1, because of the degeneracy between 𝛾 and 𝜁 crit (equation 10), we fix 𝛾 =
75 erg cm − and vary 𝜁 crit = 𝜁 crit =
10 Å as a fiducial value. Since thegrain–grain collision velocities are typically less than 50 m s − inthe dense ISM, we assume for coagulation that the grains alwaysstick when they collide (Wada et al. 2013).The parameter that regulates the maximum compaction, 𝜖 𝑉 (equa-tion 14) is also an unknown parameter. As shown later, this parameterbasically determines the filling factor of grains larger than submicron.Although this parameter is uncertain, we argue later that 𝜖 𝑉 .
1. If 𝜖 𝑉 is larger than 1, it imposes an artificial fluffiness for submi-cron grains as we discuss in Section 3.1. We adopt 𝜖 𝑉 = . 𝜖 𝑉 = 𝑛 c , is also unknown since we donot trace each monomer. We interpret this parameter as the numberof contacting points whose motion significantly reduces the porosity.We assume 𝑛 c =
30 unless otherwise stated. We also examine theeffect of 𝑛 c by changing its value. In the above, we assumed for shattering that the fragments are com-pact while the remnant has the same filling factor as the original grain(equation 19). However, compaction could occur for the remnants.Here, for an experimental purpose, we consider a model in which thecomparable volume to the ejecta suffers compaction. Noting that theejecta have a fraction 𝜑 /( + 𝜑 ) of the original grain, the compactionof the volume corresponding to that fraction can be written as 𝑉 + = max (cid:20) − 𝜑 + 𝜑 𝑉 ( 𝑚 ) + 𝜑 + 𝜑 𝑚 𝑠 , 𝑚 ( + 𝜑 ) 𝑠 (cid:21) if 𝑚 = 𝑚 /( + 𝜑 ) . (24)The first case in the max function describes the replacement of thevolume [ 𝜑 /( + 𝜑 )] 𝑉 ( 𝑚 ) (the original total volume of the ejecta)with [ 𝜑 /( + 𝜑 )] 𝑚 / 𝑠 (the corresponding compact volume). How-ever, this breaks down for large 𝜑 , and 𝑉 + could even becomenegative. Thus, we set the second case in the max function, which de-scribes the fully compact remnant. By default, we use equation (19),but when we include compaction of remnants, we use equation (24)instead of the first condition in equation (19). Coagulation and shattering basically conserve the total grain mass.Strictly speaking, because we set the upper and lower boundariesfor the grain radii ( 𝑎 min = 𝑎 max = 𝜇 m, respectively;Section 2.1) and remove all grains coagulated or shattered beyondthe boundaries, the total grain mass is not strictly conserved. Exceptfor this effect, our algorithm guarantees the conservation of the totalgrain mass (see appendix B of Hirashita & Aoyama 2019). We adopta power-law grain size distribution similar to the MRN distribution: 𝑛 init ( 𝑎 𝑚 ) ∝ 𝑎 − 𝑝𝑚 ( 𝑝 = .
5) with the lower and upper bounds ofgrain radii being 𝑎 min , ini = . 𝜇 m and 𝑎 max , ini = . 𝜇 m,respectively. This initial grain size distribution is related to theabove grain mass distribution as 𝑛 init ( 𝑎 𝑚 ) d 𝑎 𝑚 = ˜ 𝑛 ( 𝑚, 𝑡 = ) d 𝑚 .With the above power-law grain size distribution, we obtain, recallingthat 𝜌 ( 𝑚, 𝑡 ) = 𝑚 ˜ 𝑛 ( 𝑚, 𝑡 ) , 𝜌 ( 𝑚, 𝑡 = ) = ( − 𝑝 ) 𝜇 H 𝑚 H 𝑛 H D h 𝑚 ( − 𝑝 )/ , ini − 𝑚 ( − 𝑝 )/ , ini i 𝑚 (− 𝑝 + )/ (25)for 𝑚 min , ini < 𝑚 < 𝑚 max , ini , where 𝑚 max , ini = 𝜋𝑎 , ini 𝑠 / 𝑚 min , ini = 𝜋𝑎 , ini 𝑠 /
3. We used equation (3) for normalization.Thus, if we give D , we set the initial condition. We adopt D = . Note that 𝑎 min / max , ini and 𝑎 min / max are different with the latter chosen tobe similar to the constraint from MRN. We also confirmed that the resultsbelow are insensitive to 𝑎 min , ini as long as 𝑎 min , ini . . 𝜇 m. By the natureof coagulation, grains with 𝑎 < 𝑎 min , ini do not appear in the pure coagulationcases below. MNRAS000
3. We used equation (3) for normalization.Thus, if we give D , we set the initial condition. We adopt D = . Note that 𝑎 min / max , ini and 𝑎 min / max are different with the latter chosen tobe similar to the constraint from MRN. We also confirmed that the resultsbelow are insensitive to 𝑎 min , ini as long as 𝑎 min , ini . . 𝜇 m. By the natureof coagulation, grains with 𝑎 < 𝑎 min , ini do not appear in the pure coagulationcases below. MNRAS000 , 1–17 (2020) H. Hirashita as a typical dust-to-gas ratio in the Milly Way, but remind the readerthat the time-scales of coagulation and shattering scale with D − .The initial filling factor is assumed to be the same (unity unlessotherwise stated) for all grains because we do not know the fillingfactor in advance.Note that the initial condition here is not meant to represent theinitial grain size distribution in a galaxy. It is aimed at setting a ‘stan-dard’ grain size distribution, based on which we investigate the effectof coagulation and shattering. Therefore, the MRN grain size distri-bution, which roughly reproduces the dust properties (e.g. extinctioncurves) in the Milky Way is simply used as a starting point. Thisinitial condition also serves to understand how the extinction curvesshould be modified by coagulation and shattering. This approachis similar to the one taken by Hirashita & Yan (2009) for compactgrains. We show the resulting grain size distributions, filling factors andextinction curves in this section. We first examine coagulation andshattering separately to clarify the effects of each process (these casesare referred to pure coagulation and pure shattering). After that, weclarify the combined effects arising from the ISM phase exchange byincluding both processes simultaneously (Section 2.5).In the following, the grain size distribution is always shown inthe form of 𝑎 𝑚 𝑛 ( 𝑎 𝑚 )/ 𝑛 H (the variable 𝑡 in 𝑛 is omitted), where thegrain size distribution 𝑛 ( 𝑎 𝑚 ) is defined as 𝑛 ( 𝑎 𝑚 ) d 𝑎 𝑚 = ˜ 𝑛 ( 𝑚 ) d 𝑚 [recall that 𝑚 = ( 𝜋 / ) 𝑎 𝑚 𝑠 ]. Since 𝜌 ( 𝑚 ) d 𝑚 = 𝑚 ˜ 𝑛 ( 𝑚 ) d 𝑚 ∝ 𝑎 𝑚 𝑛 ( 𝑎 𝑚 ) d 𝑎 𝑚 ∝ 𝑎 𝑚 𝑛 ( 𝑎 𝑚 ) d log 𝑎 𝑚 , 𝑎 𝑚 𝑛 ( 𝑎 ) is proportional to themass distribution function per log 𝑎 𝑚 . We further divide this quan-tity by 𝑛 H to cancel out the overall density difference between thedense and diffuse ISM. We refer to 𝑎 𝑚 𝑛 ( 𝑎 𝑚 )/ 𝑛 H also as the grainsize distribution whenever there is no risk of confusion. We examine the evolution of grain size distribution and filling factorby coagulation in the dense ISM (the pure coagulation case). We showthe results in Fig. 1. Coagulation continuously forms large grains, de-pleting small grains. The bump in the grain size distribution at largegrain radii is prominent after coagulation takes place significantly.The height of the bump is almost constant, reflecting the mass conser-vation in coagulation. Micron-sized grains are formed on a time-scaleof 100 Myr, which is consistent with Hirashita & Voshchinnikov(2014, see their figure 4a). Note that this time-scale scales with theassumed density roughly as ∝ 𝑛 − / as mentioned in Section 2.5(recall that we adopt 𝑛 H = cm − for the dense ISM).We also show the filling factor 𝜙 𝑚 in Fig. 1. We observe that thefilling factor steadily decreases up to 𝑡 ∼
100 Myr because coagu-lation creates porosity. The decrease of filling factor is governed bycoagulation of small grains, which are, however, effectively depletedby coagulation. As the abundance of small grains becomes lower,the decrease of the filling factor is saturated. For large ( 𝑎 & . 𝜇 m)grains, the porosity is determined by the assumption of maximumcompaction regulated by 𝜖 𝑉 in equation (14). The filling factor atlarge grain radii roughly approaches 𝜙 𝑚 ∼ /( + 𝜖 𝑉 ) ≃ .
67 (re-call that 𝜖 𝑉 = . 𝑉 , comp + 𝑉 , comp < 𝑉 + < ( + 𝜖 𝑉 )( 𝑉 , comp + 𝑉 , comp ) . In thiscase, grains whose filling factors are between 1 /( + 𝜖 𝑉 ) and 1 form, Figure 1.
Evolution of grain size distribution (upper window) and filling factor(lower window) for the pure coagulation case. The grain size distributionis multiplied by 𝑎 𝑚 and divided by 𝑛 H , so that the resulting quantity isproportional to the grain mass abundance per log 𝑎 𝑚 relative to the gasmass. The solid, dotted, short-dashed, dot–dashed, triple-dot–dashed, andlong-dashed lines show the results at 𝑡 = 𝜙 𝑚 (filling factor) is alsoshown in the same line species as in the upper window. Note that we usethe mass-equivalent grain radius 𝑎 𝑚 , while the characteristic grain radius isobtained by 𝑎 ch = 𝑎 𝑚 𝜙 − / 𝑚 . contributing to raising the averaged 𝜙 𝑚 above 1 /( + 𝜖 𝑉 ) at largegrain radii.We also examine the dependence on the parameters relevant forcoagulation. First, we present the effect of 𝐸 roll , which regulatescompaction. As indicated in Section 2.3.1, 𝐸 roll is determined bythe critical displacement 𝜉 crit in our model. In Fig. 2a, we show theresults for 𝜉 crit =
2, 10, and 30 Å (note that 𝐸 roll is proportional to 𝜉 crit ). We only show the results at 𝑡 =
100 Myr since the effect of 𝜉 crit is qualitatively similar at all ages. We observe that the filling factorat 𝑎 ∼ . 𝜇 m is affected by the change of 𝜉 crit (or 𝐸 roll ), with larger 𝜉 crit showing smaller 𝜙 𝑚 . This is because compaction is less efficientin the case of larger 𝜉 crit . In contrast, the filling factor at 𝑎 . . 𝜇 mis not affected by 𝜉 crit because compaction is not important in thatgrain radius range. The filling factor also converges to the samevalue determined roughly by 1 /( + 𝜖 𝑉 ) at the largest grain radii asmentioned above. The grain size distribution, on the other hand, isinsensitive to 𝜉 crit . Since the grain abundance is dominated by largegrains, small grains collide predominantly with large grains. In thissituation, the collisional cross-section is governed by large grainswhich are almost compact. Thus, the difference in the porosity atsubmicron radii does not affect the grain size distribution.Grain compaction is also affected by the number of contacts 𝑛 c (equation 12). In Fig. 2b, we show the results for 𝑛 c =
10, 30, and 100at 𝑡 =
100 Myr. Naturally, the effect of 𝑛 c appears at grain radii wherecompaction is important. As 𝑛 c becomes larger, the filling factor iskept lower against compaction, so that the increase of 𝜙 𝑚 becomesshallower at large 𝑎 𝑚 . The grain size distribution is insensitive alsoto the change of 𝑛 c .We also examine the dependence on 𝜖 𝑉 , which regulates the max-imum compaction. As we argued above, the filling factor roughlyconverges to 1 /( + 𝜖 𝑉 ) at large grain radii, where grain velocitiesare high enough for significant compaction. As shown above, the MNRAS , 1–17 (2020) orosity of interstellar dust Figure 2.
Parameter dependence for the pure coagulation cases. We presentthe grain size distributions (upper window) and the filling factors (lowerwindow) at 𝑡 =
100 Myr. (a) Dependence on 𝜉 crit , which regulates 𝐸 roll . Thesolid, dotted, and dashed lines show the results for 𝜉 crit =
2, 10, and 30 Å,respectively. (b) Dependence on 𝑛 c , which regulates compaction. The solid,dotted, and dashed lines show the results for 𝑛 c =
10, 30, and 100, respectively.(c) Dependence on 𝜖 𝑉 , which regulates the maximum compaction in high-velocity collisions (relevant for large grains). The solid, dotted, and dashedlines show the results for 𝜖 𝑉 =
0, 0.5, and 1, respectively. Other than thevaried parameter in each panel, we adopt the fiducial values ( 𝜉 crit =
10 Å, 𝑛 c =
30, and 𝜖 𝑉 = . Figure 3.
Same as Fig. 1 but for the pure shattering model. We show the caseswhere we (a) do not consider and (b) consider compaction for the shatteredremnants. minimum value of 𝜙 𝑚 is around 0.3–0.5. Thus, 𝜖 𝑉 should be smallerthan ∼
1; otherwise, the filling factor at large grain radii decreasesbelow the minimum value, which means that we add artificial (orcontradictory) porosity to the grains in compaction. Thus, we ex-amine 𝜖 𝑉 =
0, 0.5, and 1 at 𝑡 =
100 Myr in Fig. 2c. We confirmthat the effect of 𝜖 𝑉 appears at large grain radii ( 𝑎 & . 𝜇 m). Asmentioned above, 𝜙 𝑚 at large grain radii is roughly 1 /( + 𝜖 𝑉 ) . Thegrain size distribution is affected by 𝜖 𝑉 since larger porosity makesthe collision kernel larger by a factor of 𝜙 − / 𝑚 (Section 2.2), whichincreases the grain–grain collision rate. We show the evolution of grain size distribution and filling factorby shattering in the diffuse ISM. Since shattering does not createnew porosity in grains, we start with 𝜙 𝑚 = . 𝜙 𝑚 = MNRAS000
100 Myr in Fig. 2c. We confirmthat the effect of 𝜖 𝑉 appears at large grain radii ( 𝑎 & . 𝜇 m). Asmentioned above, 𝜙 𝑚 at large grain radii is roughly 1 /( + 𝜖 𝑉 ) . Thegrain size distribution is affected by 𝜖 𝑉 since larger porosity makesthe collision kernel larger by a factor of 𝜙 − / 𝑚 (Section 2.2), whichincreases the grain–grain collision rate. We show the evolution of grain size distribution and filling factorby shattering in the diffuse ISM. Since shattering does not createnew porosity in grains, we start with 𝜙 𝑚 = . 𝜙 𝑚 = MNRAS000 , 1–17 (2020) H. Hirashita grain radius 𝑎 = 𝑡 ∼
100 Myr, the grain abundance isdominated by small grains. The filling factor increases, especiallyat small grain radii, because shattered fragments are compact. Sincewe assume that the remnants have the same filling factor as theoriginal grains, the filling factors of the largest grains, which aredominated by shattered remnants, do not change. Recall that thelargest fragment size is ( . ) / ≃ .
27 times the original grainsize (Section 2.3.2). Since the maximum grain radius in the initialcondition is 0.25 𝜇 m, the largest grain radius where 𝜙 𝑚 is raised byfragments is 𝑎 ≃ . × . 𝜇 m ≃ . 𝜇 m. Thus, the porosity isonly changed from its initial value at 𝑎 < . 𝜇 m. All the grainswith 𝑎 < 𝑎 min , ini = . 𝜇 m are newly formed by shattering; thus,they always have 𝜙 𝑚 = If we consider the grain size distribution in a region large enoughto include both the dense and diffuse ISM (or if we consider time-scales much longer than the phase-exchange time; Section 2.5), theeffect of coagulation and that of shattering coexist. In our one-zonemodel, we could simulate this coexisting effect by simultaneouslytreating coagulation and shattering at each time-step with a weightof 𝜂 dense : ( − 𝜂 dense ) as formulated in equations (22) and (23).We adopt 𝜂 dense = . 𝑡 ∼
100 Myr. The filling factor changes in a way similar to the purecoagulation case (Fig. 1) but the decrease of the filling factor pro-ceeds more if both coagulation and shattering are present becausesmall grains which coagulate to form porosity are continuously sup-plied. The filling factor at large grain radii converges to ∼ /( + 𝜖 𝑉 ) as explained in Section 3.1.To examine the balance between coagulation and shattering, wealso show the results for different values of 𝜂 dense (0.1 and 0.9) inFig. 4. Naturally, a higher abundance of large grains is obtained forlarger 𝜂 dense . For 𝜂 dense = .
9, the result is similar to that in the purecoagulation case (Fig. 1), except that small grains continue to be pro-duced by shattering. The case of 𝜂 dense = . 𝑡 ∼
100 Myr but large grains increase (‘re-form’)
Figure 4.
Same as in Fig. 1, but including both shattering and coagulation.Panels (a), (b) and (c) show the results for 𝜂 dense = . 𝜙 𝑚 between this figure and Fig.1.MNRAS , 1–17 (2020) orosity of interstellar dust Figure 5.
Same as Fig. 4b but with compaction of shattered remnants. after that. (ii) The filling factor becomes very small at 𝑡 >
100 Myr.This is because the production of small grains by shattering contin-uously activates coagulation, which creates porosity. The increasedporosity makes the geometrical cross-section larger, and enhancescoagulation. Moreover, since the velocities of porous grains are re-duced, compaction does not occur as efficiently as in the pure co-agulation case. Thus, efficient shattering together with coagulationmakes the interstellar dust highly porous and further activates co-agulation. This interesting interplay is further discussed in Section5. The dependences on the parameters related to coagulation ( 𝜉 crit , 𝑛 c and 𝜖 𝑉 ) are similar to those shown in Section 3.1 (Fig. 2). On theother hand, compaction of shattered remnants (Section 2.5.2) coulddecrease the porosity if 𝜂 dense is low. In Fig. 5, we show the resultswith compaction of shattered remnants for 𝜂 dense = .
1. This figureshould be compared with Fig. 4b. We observe that the filling factorbecomes higher at 𝑡 ∼
300 Myr at large grain radii, although theresults are almost unchanged at younger ages. Compaction of rem-nants also makes coagulation at later epochs less efficient, producingless large grains at 𝑡 ∼
300 Myr. However, the filling factor is still assmall as 0.1–0.2 at 𝑎 ∼ . 𝜇 m at 𝑡 ∼
300 Myr; thus, the conclusionthat efficient shattering together with coagulation helps to increaseporosity is robust.
Based on the grain size distributions and the filling factors presentedabove, we calculate the extinction curves by the method in Section 2.4for silicate and amC. As mentioned above, the extinction 𝐴 𝜆 is pre-sented in two ways: 𝐴 𝜆 / 𝑁 H and 𝐴 𝜆 / 𝐴 𝑉 . To clarify the effect ofporosity, we also calculate extinction curve by forcing the filling fac-tor of all grains to be unity with 𝑎 𝑚 fixed. The extinction calculatedin this way (i.e. 𝜙 𝑚 =
1) is denoted as 𝐴 𝜆, . The effect of porosity ispresented by 𝐴 𝜆 / 𝐴 𝜆, . In Fig. 6, we present the extinction curves corresponding to the grainsize distributions and the filling factors for the pure coagulation case
Figure 6.
Extinction curves for silicate and amC in the upper and lower panels,respectively. The solid, dotted, short-dashed, dot–dashed, triple-dot–dashed,and long-dashed lines (thick lines) show the extinction curves at 𝑡 = 𝐴 𝜆, (extinction with 𝜙 𝑚 = 𝐴 𝜆 = 𝐴 𝜆, at 𝑡 =
0. In each panel, the upper, middle, and lower windows present theextinction per hydrogen, the extinction normalized to the 𝑉 -band value, andthe ratio of 𝐴 𝜆 to 𝐴 𝜆, (an indicator of the porosity effect), respectively.MNRAS000
0. In each panel, the upper, middle, and lower windows present theextinction per hydrogen, the extinction normalized to the 𝑉 -band value, andthe ratio of 𝐴 𝜆 to 𝐴 𝜆, (an indicator of the porosity effect), respectively.MNRAS000 , 1–17 (2020) H. Hirashita (Fig. 1). As expected from the evolution of grain size distribution,the extinction curve becomes flatter as coagulation proceeds. Sincea large fraction of the grains become bigger than ∼ . 𝜇 m aftercoagulation, the extinction per hydrogen declines at wavelengthsshorter than ∼ 𝜋 × . 𝜇 m for both materials. The effect of porosity,which appears in the difference between 𝐴 𝜆 (thick lines) and 𝐴 𝜆, (thin lines) or in the ratio 𝐴 𝜆 / 𝐴 𝜆, shown in the bottom window,is clear. The extinction is 10–20 per cent higher in the porous casethan in the non-porous case in a large part of wavelengths for bothmaterials.There are some differences between the two materials. For sili-cate, porous grains have smaller extinction than compact ones in acertain wavelength range (e.g. 1 / 𝜆 ∼ . 𝜇 m − at 𝑡 =
10 Myr) andthis range shifts to longer wavelengths with age. Voshchinnikov et al.(2006) and Shen et al. (2008) showed that the porosity could bothincrease and decrease the extinction cross-section depending onthe wavelength. The wavelength range of suppressed extinction byporosity is roughly consistent with their results. Voshchinnikov et al.(2006) also showed that the wavelength range where the extinctioncross-section is decreased by porosity shifts to longer wavelengthsas the grains become larger. This is consistent with our result. At 𝑡 .
100 Myr, the normalized extinction curves ( 𝐴 𝜆 / 𝐴 𝑉 ) becomessteeper by the porosity because the increase of extinction is moreenhanced in the UV than in the 𝑉 band.For amC, porosity always increases the extinction; that is, 𝐴 𝜆 / 𝐴 𝜆, is always larger than 1. Thus, if we normalize the extinction to 𝐴 𝑉 ,the porosity effect roughly cancels out, so that the extinction curveshape is insensitive to the porosity (filling factor) for amC.We also examined the parameter dependence (not shown) andconfirmed that the extinction curves are not sensitive to 𝜉 crit or 𝜖 𝑉 in the optical and UV with differences less than 5 per cent. However,at 1 / 𝜆 < 𝜇 m − , the difference could be as large as 20 per cent,since, at such long wavelengths, the porosity of large grains, whichis regulated by 𝜉 crit and 𝜖 𝑉 , is important. The changes driven by thedifference in those parameters are sub-dominant compared with thedifference 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 com-pact grains. More important is the interplay between shattering andcoagulation as shown above. Indeed, coagulation of newly createdsmall grains by shattering produces porous grains, contributing to theincrease of porosity in the interstellar dust. Here, we investigate theeffect of relative strength between shattering and coagulation; that is,we compare the results for different dense gas fractions, 𝜂 dense .In Fig. 7, we compare the extinction curves calculated for themodels including both coagulation and shattering. The correspond-ing grain size distributions are shown in Section 3.3 (Fig. 4). Wecompare the results at the same age 𝑡 =
100 Myr. We observe thatthe steepness of extinction curve is very sensitive to 𝜂 dense . This is anatural consequence of the different grain size distributions for vari-ous 𝜂 dense . We also plot the extinction curves with 𝜙 𝑚 = 𝑎 𝑚 ), that is, 𝐴 𝜆, . Comparing 𝐴 𝜆 with 𝐴 𝜆, ,we observe for silicate that the porosity increases the extinction atshort and long wavelengths, but that it rather decreases the extinc-tion at intermediate wavelengths as already noted in the previoussubsection. The wavelength range where the porosity decreases theextinction shifts to shorter wavelengths for smaller 𝜂 dense , which isconsistent with the above mentioned tendency that porosity of smallergrains reduces the extinction at shorter wavelengths. For amC, the Figure 7.
Comparison among the extinction curves in the models includingboth coagulation and shattering. The relative strengths of these two processesare regulated by the dense gas fraction 𝜂 dense . We show the same quantities asin Fig. 6 at 𝑡 =
100 Myr. The solid, dotted, and dashed lines show the resultsfor 𝜂 dense = .
1, 0.5, and 0.9, respectively. The thick and thin lines presentthe cases using the calculated filling factors (Fig. 4) and those with 𝜙 𝑚 = 𝐴 𝜆, ). The upper and lower panelsare for silicate and amC, respectively.MNRAS , 1–17 (2020) orosity of interstellar dust porosity always increases the extinction, but the largest porosity (i.e. 𝜂 dense = .
1) does not necessarily mean the largest opacity enhance-ment ( 𝐴 𝜆 / 𝐴 𝜆, ) in the UV. At long wavelengths, larger porosityindicates larger extinction.Overall, the effect of porosity on the extinction curve is less than20 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 torealize an observed extinction if we take porosity into account. Thedifference is not large, though, because the increase of grain volumeand the ‘dilution’ of permittivity have opposite effects on the extinc-tion (Li 2005). For further constraint on the porosity evolution, weneed a comprehensive analysis of observed dust properties includingthe infrared SED (Dwek et al. 1997). Polarization may also help tofurther constrain the models. Since such a detailed comparison alsoneeds further modelling of the mixture of various dust species andof the interstellar radiation field, we leave it for a future work. The porosity increases the effective sizes of grains, enhancing thegrain–grain collision rate. If only coagulation occurs, however, theporosity, which is large at 𝑎 . . 𝜇 m, does not have a large in-fluence on the grain size distribution, since grains quickly coagulateto achieve 𝑎 & . 𝜇 m, and are affected by compaction. Note that,if we individually observe clouds denser than 𝑛 H ∼ cm − , wecould 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 𝜇 m (e.g.Weingartner & Draine 2001b), so that the grains in such dense cloudsare not likely to contribute directly to the interstellar dust population.Coagulation can be balanced by shattering. For 𝜂 cold = .
5, thegrain size distribution and the filling factor both converge to equilib-rium distributions at 𝑡 ∼
100 Myr, while coagulation continuouslyincreases the grain radii for 𝜂 cold = . 𝜂 cold produce similar results to the pure coagulationcase. For the case of 𝜂 cold = .
5, we performed a calculation byforcing 𝜙 𝑚 to be always unity (not shown) and found that the grainsize distribution changes little by the different treatment of 𝜙 𝑚 (i.e.compared with the results shown in Fig. 4a). This is because theevolution of grain size distribution is regulated by the formation oflarge grains (recall that the grain abundance is dominated by largegrains), which have small porosity because of compaction.Porosity plays a critical role in the case of small 𝜂 dense = . 𝑡 ∼
100 Myr). We arguedabove that coagulation is activated at 𝑡 ∼
100 Myr because theincreased porosity enhances the grain cross-sections (the grain–graincollision rate). For demonstration, we compare two calculations inFig. 8: one is the same as above for 𝜂 dense = . 𝜙 𝑚 = 𝑡 =
100 Myr. If we fix 𝜙 𝑚 =
1, grains at 𝑎 & . 𝜇 m are simply depleted by shatteringas expected from weak coagulation in 𝜂 dense = .
1. In contrast, ifwe take the evolution of 𝜙 𝑚 into account, grains at 𝑎 & . 𝜇 mare re-formed at 𝑡 &
200 Myr, reaching roughly an equilibrium at 𝑡 ∼
400 Myr. The filling factor reaches the smallest value ( 𝜙 𝑚 ∼ . Figure 8.
Long-term evolution of grain size distribution and filling factorwith both coagulation and shattering under 𝜂 dense = .
1. Panels (a) and (b)show the results with including the evolution of filling factor (i.e. the samemodel as in Fig. 4b) and with fixing 𝜙 𝑚 =
1, respectively. We present theresults at 𝑡 =
0, 100, 200, 300, 400, and 500 Myr (at 𝑡 <
100 Myr, the tworesults have little difference) by the solid, dotted, short-dashed, dot–dashed,triple-dot–dashed, and long-dashed lines, respectively. at 𝑎 ∼ . 𝜇 m. There are two effects that promote coagulation at laterstages: (i) The increase of grain cross-sections by porosity enhancesthe grain–grain collision rate because the collision kernel scaleswith the porosity as ∝ 𝜙 − / 𝑚 (Section 2.2). (ii) The decreased grainvelocity ( ∝ 𝜙 / 𝑚 ; equation 6) makes compaction in coagulation atlarge grain radii less effective, so that the filling factor is kept smalleven at 𝑎 ∼ . 𝜇 m. These two effects are persistent qualitativelyeven if we consider the compaction of shattered remnants (see Fig.5).Recall that our treatment of the two-phase ISM is based on theparameter 𝜂 dense , which sets the fraction of time dust spends in thedense phase. Thus, the above results imply that the dust evolutionin a condition where both coagulation and shattering coexist (in awide area of the ISM and/or on a time-scale longer than the phaseexchange time-scale ∼ yr; Section 2.5) could be strongly affectedby porosity. In other words, dust evolution models which do notinclude porosity evolution could predict a very different evolution MNRAS000
100 Myr, the tworesults have little difference) by the solid, dotted, short-dashed, dot–dashed,triple-dot–dashed, and long-dashed lines, respectively. at 𝑎 ∼ . 𝜇 m. There are two effects that promote coagulation at laterstages: (i) The increase of grain cross-sections by porosity enhancesthe grain–grain collision rate because the collision kernel scaleswith the porosity as ∝ 𝜙 − / 𝑚 (Section 2.2). (ii) The decreased grainvelocity ( ∝ 𝜙 / 𝑚 ; equation 6) makes compaction in coagulation atlarge grain radii less effective, so that the filling factor is kept smalleven at 𝑎 ∼ . 𝜇 m. These two effects are persistent qualitativelyeven if we consider the compaction of shattered remnants (see Fig.5).Recall that our treatment of the two-phase ISM is based on theparameter 𝜂 dense , which sets the fraction of time dust spends in thedense phase. Thus, the above results imply that the dust evolutionin a condition where both coagulation and shattering coexist (in awide area of the ISM and/or on a time-scale longer than the phaseexchange time-scale ∼ yr; Section 2.5) could be strongly affectedby porosity. In other words, dust evolution models which do notinclude porosity evolution could predict a very different evolution MNRAS000 , 1–17 (2020) H. Hirashita of the grain size distribution from those which correctly take theporosity evolution into account.For a long-term evolution of dust in the ISM, dust enrichmentby 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 porosityevolution in these processes. This paper provides a first importantstep for the understanding of porosity evolution since coagulationand shattering are mechanisms of efficiently modifying the grainsize distribution. Indeed, Hirashita & Aoyama (2019) emphasizedthe importance of these two processes in realizing MRN-like grainsize distributions (see also Aoyama et al. 2020).
As shown in Section 4, porosity affects the UV–optical extinctioncurves by 10–20 per cent. As noted by Voshchinnikov et al. (2006),porosity does not necessarily increase the extinction (see also Jones1988; Li 2005). At long-optical and near infrared wavelengths, theopacity of silicate can decrease owing to the porosity. The wavelengthrange where this decrease occurs shifts towards shorter wavelengthsas the typical grain radius becomes smaller in e.g. a shattering-dominated condition. The extinction of amC is enhanced by 10–20per cent at all wavelengths. The above results indicate that we couldsave 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 pronouncedgiven that the evolution of grain size distribution is affected by thefilling factor. As shown above, porosity has the largest effect on thegrain size distribution if strong shattering and weak coagulation areboth present such as in the case of 𝜂 dense = .
1. Thus, the extinctioncurve is affected by porosity not only through the optical propertiesbut also through the modification of grain size distribution. To presentthis effect, we show in Fig. 9 the extinction curves corresponding tothe 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 shat-tering (recall that we remove the grains which become smaller than 𝑎 = 𝐴 𝜆 / 𝑁 H just because of the boundarycondition. Thus, we only show 𝐴 𝜆 / 𝐴 𝑉 , which is free from the effectof dust mass loss (i.e. we could purely observe the effect of grainsize distribution). For the case without porosity evolution ( 𝜙 𝑚 = 𝑉 band, where the extinction curve isnormalized, while porosity enhances the UV extinction by 30–40 percent. Jones (1988) and Voshchinnikov et al. (2006) also showed thatporosity does not change the optical extinction much but increasethe UV and infrared opacities. This effect makes the UV extinctioncurve normalized to the 𝑉 -band value steep.The case shown in Figs. 8 and 9 indicates that porosity evolutioncan have a dramatic impact on the shape of extinction curves. Even ifthe grain size distribution becomes the one similar to the initial grainsize distribution at later stages of evolution, the small porosity makesthe extinction curves significantly steeper than the initial one. Inparticular, the grains contributing to the optical extinction have radii 𝑎 ∼ 𝜆 /( 𝜋 ) ∼ . 𝜇 m, where the porosity is the largest. Therefore,if we consider the creation of porosity through the interplay between Figure 9.
Evolution of extinction curve for the model with 𝜂 dense = . 𝜙 𝑚 to be always unity (Fig.8b). Note that 𝜙 𝑚 = 𝑡 = coagulation and shattering, the evolution of extinction curve couldbe qualitatively very different. Moreover, the porosity depends on thegrain radius; this dependence also creates ‘higher-order’ wavelengthdependence of extinction curve. As shown above, the wavelengthrange where the extinction is reduced by porosity shifts to longerwavelengths if porosity is developed in larger grains.As shown in Section 3.3 (Fig. 5), if we consider compaction ofshattered remnants, the filling factor increases a little at 𝑡 &
300 Myr.We also calculated the evolution of extinction curve with remnantcompaction (now shown). Since coagulation is less efficient in thiscase, the extinction curves are steeper than those shown in Fig. 9. Ifwe only see the extinction curve shape in the UV–optical, a largerfilling factor and a lower abundance of large grains are degenerate.
We formulate and calculate the evolution of grain size distributionand filling factor (porosity) through coagulation and shattering in theISM. We adopt the 2D Smoluchowski equation to solve the distribu-
MNRAS , 1–17 (2020) orosity of interstellar dust tion functions of grain size and filling factor. To save the computa-tional time, we only treat the mean filling factor for each grain radiusbased on O09 and O12. For coagulation, the transition from the hit-and-stick to compaction regime are characterized and modelled bycomparing the impact energy with the rolling energy. Shattering istreated as a formation mechanism of small compact fragments. Forshattered remnants, we basically assume the same porosity as theoriginal grain but we also examine the case where the volume equalto the colliding particle suffers perfect compaction. We assume thatcoagulation and shattering are hosted by the dense and diffuse ISM,respectively.For the pure coagulation case (without shattering), the porositydevelops around 𝑎 ∼ . 𝜇 m, where a low filling factor of 𝜙 𝑚 ∼ . 𝑎 > . 𝜇 m, where compaction occurs. Therefore, the porosity littleaffects the evolution of grain size distribution if only coagulationis present. For the pure shattering case, we confirm that shatteringtends to make the filling factors asymptotically approach 𝜙 𝑚 = 𝑎 . . 𝜇 m, although those at 𝑎 & . 𝜇 m depend on thetreatment of compaction for shattered remnants.Next, we examine the case where coagulation and shattering areboth present. This corresponds to a situation where we consider awide enough area in the ISM which contains both the dense anddiffuse ISM, or where the time-scale of interest is much longer thanthe exchange time of the ISM phases. We find that the filling factordrops even below 0.1 around 𝑎 ∼ . 𝜇 m. The porosity evolutionis sensitive to the relative efficiency of coagulation to shattering,which is regulated by the dense gas fraction, 𝜂 dense . The filling factortends to be small if shattering is stronger (e.g. 𝜂 dense = . 𝜂 dense = .
1, coagulation is activated in later stages(after porosity develops) because the high porosity enhances the graincross-sections. Compaction is not efficient in this case since the grainvelocity is diminished by the increased porosity. Thus, large grainsare kept porous in this case.We also calculate the evolution of extinction curve using the EMT.Porosity formed as a result of coagulation enhances the UV and in-frared extinction by ∼ 𝜂 dense = .
1) shows a recreation of large grains at later stages asmentioned above. In this case, although the grain size distributionitself is similar to the MRN distribution, the extinction curve shapestays steep if we normalize the extinction to the 𝑉 band value. This isbecause porosity makes the grains relatively ‘transparent’ in the op-tical, while the extinction is enhanced in the UV. Thus, the steepnessof extinction curve is also affected by the porosity evolution.Although the predicted features in the extinction curves could becompared with observations, our model developed in this paper isstill premature for detailed comparison. There are two necessary ex-tensions of our modelling. First, we could include other processeswhich also play an important role in the evolution of grain size dis-tribution; that is, dust production by stellar ejecta, dust destructionby supernova shocks, and dust growth by the accretion of gas-phasemetals. Secondly, since the filling factor and the grain size distri- bution could be degenerate in the resulting extinction curve shape,it is desirable to predict other independent properties such as in-frared emission SED and polarization. We emphasize that the basicframework developed in this paper provides a basis on which we willextend our predictions. ACKNOWLEDGEMENTS
We are grateful to S. Okuzumi, Y. Matsumoto, and the anonymousreferee for useful comments. HH thanks the Ministry of Science andTechnology for support through grant MOST 107-2923-M-001-003-MY3 and MOST 108-2112-M-001-007-MY3, and the AcademiaSinica for Investigator Award AS-IA-109-M02. VBI acknowledgesthe support from the RFBR grant 18-52-52006 and the SUAI grantFSRF-2020-0004. LP acknowledges support from the ProgrammeNational ‘Physique et Chimie du Milieu Interstellaire’ (PCMI) ofCNRS/INSU with INC/INP co-funded by CEA and CNES and fromAction Fédératrice Astrochimie de l’Observatoire de Paris. HH ac-knowledges the financial support and the hospitality of l’Observatoirede Paris during his stay.
DATA AVAILABILITY
The data underlying this article are available in Figshare at https://doi.org/10.6084/m9.figshare.12917375.v1 . REFERENCES
Aoyama S., Hirashita H., Nagamine K., 2020, MNRAS, 491, 3844Asano R. S., Takeuchi T. T., Hirashita H., Nozawa T., 2013, MNRAS, 432, 637Asano R. S., Takeuchi T. T., Hirashita H., Nozawa T., 2014, MNRAS, 440, 134Bohren C. F., Huffman D. R., 1983, Absorption and Scattering of Light bySmall Particles. Wiley, New YorkCazaux S., Tielens A. G. G. M., 2004, ApJ, 604, 222Chen L.-H., Hirashita H., Hou K.-C., Aoyama S., Shimizu I., Nagamine K.,2018, MNRAS, 474, 1545Désert F. X., Boulanger F., Puget J. L., 1990, A&A, 500, 313Dohnanyi J. S., 1969, J. Geophys. Res., 74, 2531Dominik C., Tielens A. G. G. M., 1997, ApJ, 480, 647Draine B. T., 1990, in Blitz L., ed., Astronomical Society of the PacificConference Series Vol. 12, The Evolution of the Interstellar Medium.Astron. Soc. Pac., San Francisco, pp 193–205Draine B. T., 2003, ApJ, 598, 1026Draine B. T., 2009, in Henning T., Grün E., Steinacker J., eds, AstronomicalSociety of the Pacific Conference Series Vol. 414, Cosmic Dust - Nearand Far. Astron. Soc. Pac., San Francisco, p. 453Draine B. T., Hensley B. S., 2020, ApJ, submitted, arXiv:2009.11314Draine B. T., Tan J. C., 2003, ApJ, 594, 347Dwek E., 1998, ApJ, 501, 643Dwek E., et al., 1997, ApJ, 475, 565Gould R. J., Salpeter E. E., 1963, ApJ, 138, 393Guillet V., Pineau Des Forêts G., Jones A. P., 2011, A&A, 527, A123Guillet V., Hennebelle P., Pineau des Forêts G., Marcowith A., CommerçonB., Marchand P., 2020, A&A, 643, A17Heim L.-O., Blum J., Preuss M., Butt H.-J., 1999, Phys. Rev. Lett., 83, 3328Hensley B. S., Draine B. T., 2020, ApJ, submitted, arXiv:2009.00018Hirashita H., 1999, ApJ, 510, L99Hirashita H., Aoyama S., 2019, MNRAS, 482, 2555Hirashita H., Kobayashi H., 2013, Earth, Planets, and Space, 65, 1083Hirashita H., Li Z.-Y., 2013, MNRAS, 434, L70Hirashita H., Voshchinnikov N. V., 2014, MNRAS, 437, 1636Hirashita H., Yan H., 2009, MNRAS, 394, 1061Hoang T., 2019, ApJ, 876, 13 MNRAS000
Aoyama S., Hirashita H., Nagamine K., 2020, MNRAS, 491, 3844Asano R. S., Takeuchi T. T., Hirashita H., Nozawa T., 2013, MNRAS, 432, 637Asano R. S., Takeuchi T. T., Hirashita H., Nozawa T., 2014, MNRAS, 440, 134Bohren C. F., Huffman D. R., 1983, Absorption and Scattering of Light bySmall Particles. Wiley, New YorkCazaux S., Tielens A. G. G. M., 2004, ApJ, 604, 222Chen L.-H., Hirashita H., Hou K.-C., Aoyama S., Shimizu I., Nagamine K.,2018, MNRAS, 474, 1545Désert F. X., Boulanger F., Puget J. L., 1990, A&A, 500, 313Dohnanyi J. S., 1969, J. Geophys. Res., 74, 2531Dominik C., Tielens A. G. G. M., 1997, ApJ, 480, 647Draine B. T., 1990, in Blitz L., ed., Astronomical Society of the PacificConference Series Vol. 12, The Evolution of the Interstellar Medium.Astron. Soc. Pac., San Francisco, pp 193–205Draine B. T., 2003, ApJ, 598, 1026Draine B. T., 2009, in Henning T., Grün E., Steinacker J., eds, AstronomicalSociety of the Pacific Conference Series Vol. 414, Cosmic Dust - Nearand Far. Astron. Soc. Pac., San Francisco, p. 453Draine B. T., Hensley B. S., 2020, ApJ, submitted, arXiv:2009.11314Draine B. T., Tan J. C., 2003, ApJ, 594, 347Dwek E., 1998, ApJ, 501, 643Dwek E., et al., 1997, ApJ, 475, 565Gould R. J., Salpeter E. E., 1963, ApJ, 138, 393Guillet V., Pineau Des Forêts G., Jones A. P., 2011, A&A, 527, A123Guillet V., Hennebelle P., Pineau des Forêts G., Marcowith A., CommerçonB., Marchand P., 2020, A&A, 643, A17Heim L.-O., Blum J., Preuss M., Butt H.-J., 1999, Phys. Rev. Lett., 83, 3328Hensley B. S., Draine B. T., 2020, ApJ, submitted, arXiv:2009.00018Hirashita H., 1999, ApJ, 510, L99Hirashita H., Aoyama S., 2019, MNRAS, 482, 2555Hirashita H., Kobayashi H., 2013, Earth, Planets, and Space, 65, 1083Hirashita H., Li Z.-Y., 2013, MNRAS, 434, L70Hirashita H., Voshchinnikov N. V., 2014, MNRAS, 437, 1636Hirashita H., Yan H., 2009, MNRAS, 394, 1061Hoang T., 2019, ApJ, 876, 13 MNRAS000 , 1–17 (2020) H. Hirashita
Inoue A. K., 2011, Earth, Planets, and Space, 63, 1027Israelachvili J., 1992, Intermolecular and Surface Forces, 2nd edn. AcademicPress, LondonJones A. P., 1988, MNRAS, 234, 209Jones A. P., Tielens A. G. G. M., Hollenbach D. J., McKee C. F., 1994, ApJ,433, 797Jones A. P., Tielens A. G. G. M., Hollenbach D. J., 1996, ApJ, 469, 740Jones A. P., Köhler M., Ysard N., Bocchio M., Verstraete L., 2017, A&A,602, A46Kataoka A., Okuzumi S., Tanaka H., Nomura H., 2014, A&A, 568, A42Kimura H., Wada K., Senshu H., Kobayashi H., 2015, ApJ, 812, 67Kobayashi H., Tanaka H., 2010, Icarus, 206, 735Kusaka T., Nakano T., Hayashi C., 1970, Progress of Theoretical Physics,44, 1580Lefèvre C., Pagani L., Ladjelate B., Min M., Hirashita H., Zylka R., 2020,in Mayet F., Catalano A., Macías-Pérez J., Perotto L., eds, EuropeanPhysical Journal Web of Conferences Vol. 228, European Physical JournalWeb of Conferences. Grenoble, p. 00013Li A., 2005, ApJ, 622, 965Li A., Draine B. T., 2001, ApJ, 554, 778Li A., Greenberg J. M., 1997, A&A, 323, 566Li A., Lunine J. I., 2003, ApJ, 590, 368Lisenfeld U., Ferrara A., 1998, ApJ, 496, 145Mathis J. S., 1996, ApJ, 472, 643Mathis J. S., Whiffen G., 1989, ApJ, 341, 808Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ, 217, 425Mathis J. S., Cohen D., Finley J. P., Krautter J., 1995, ApJ, 449, 320Mattsson L., 2016, Planet. Space Sci., 133, 107Mattsson L., De Cia A., Andersen A. C., Zafar T., 2014, MNRAS, 440, 1562McKee C., 1989, in Allamandola L. J., Tielens A. G. G. M., eds, IAU Sympo-sium Vol. 135, Interstellar Dust. Kluwer Academic Publishers, Dordrecht,p. 431O’Donnell J. E., Mathis J. S., 1997, ApJ, 479, 806Okuzumi S., Tanaka H., Sakagami M.-a., 2009, ApJ, 707, 1247Okuzumi S., Tanaka H., Kobayashi H., Wada K., 2012, ApJ, 752, 106Ormel C. W., Spaans M., Tielens A. G. G. M., 2007, A&A, 461, 215Ormel C. W., Paszun D., Dominik C., Tielens A. G. G. M., 2009, A&A,502, 845Ormel C. W., Min M., Tielens A. G. G. M., Dominik C., Paszun D., 2011,A&A, 532, A43Rouillé G., Jäger C., Krasnokutski S. A., Krebsz M., Henning T., 2014,Faraday Discuss., 168, 449Rouillé G., Jäger C., Henning T., 2020, ApJ, 892, 96Shen Y., Draine B. T., Johnson E. T., 2008, ApJ, 689, 260Silsbee K., Ivlev A. V., Sipilä O., Caselli P., Zhao B., 2020, A&A, 641, A39Smith R. K., Dwek E., 1998, ApJ, 503, 831Steinpilz T., Teiser J., Wurm G., 2019, ApJ, 874, 60Suyama T., Wada K., Tanaka H., Okuzumi S., 2012, ApJ, 753, 115Takeuchi T. T., Ishii T. T., Nozawa T., Kozasa T., Hirashita H., 2005, MNRAS,362, 592Tanaka H., Inaba S., Nakazawa K., 1996, Icarus, 123, 450Tazaki R., Tanaka H., Okuzumi S., Kataoka A., Nomura H., 2016, ApJ,823, 70Voelk H. J., Jones F. C., Morfill G. E., Roeser S., 1980, A&A, 85, 316Voshchinnikov N. V., Hirashita H., 2014, MNRAS, 445, 301Voshchinnikov N. V., Il’in V. B., Henning T., 2005, A&A, 429, 371Voshchinnikov N. V., Il’in V. B., Henning T., Dubkova D. N., 2006, A&A,445, 167Wada K., Tanaka H., Suyama T., Kimura H., Yamamoto T., 2007, ApJ,661, 320Wada K., Tanaka H., Suyama T., Kimura H., Yamamoto T., 2008, ApJ,677, 1296Wada K., Tanaka H., Okuzumi S., Kobayashi H., Suyama T., Kimura H.,Yamamoto T., 2013, A&A, 559, A62Weingartner J. C., Draine B. T., 2001a, ApJS, 134, 263Weingartner J. C., Draine B. T., 2001b, ApJ, 548, 296Williams D. R., Wetherill G. W., 1994, Icarus, 107, 117Wong Y. H. V., Hirashita H., Li Z.-Y., 2016, PASJ, 68, 67 Woo J. W., Clark G. W., Day C. S. R., Nagase F., Takeshima T., 1994, ApJ,436, L5Yamasawa D., Habe A., Kozasa T., Nozawa T., Hirashita H., Umeda H.,Nomoto K., 2011, ApJ, 735, 44Yan H., Lazarian A., Draine B. T., 2004, ApJ, 616, 895Zubko V. G., Mennella V., Colangeli L., Bussoletti E., 1996, MNRAS,282, 1321Zuo W. B., Li A., Zhao G., 2020, ApJS, in press, arXiv:2011.09440
APPENDIX A: DERIVATION OF THE BASIC EQUATIONS
We derive equations (4) and (5) based on O09. We start from the2D Smoluchowski equation generalized to treat shattering as well ascoagulation: 𝜕 𝑓 ( 𝑰 , 𝑡 ) 𝜕𝑡 = − 𝑓 ( 𝑰 , 𝑡 ) ∫ 𝐾 ( 𝑰 ; 𝑰 ′ ) 𝑓 ( 𝑰 ′ , 𝑡 ) d 𝑰 ′ + ∬ 𝐾 ( 𝑰 ; 𝑰 ) 𝑓 ( 𝑰 , 𝑡 ) 𝑓 ( 𝑰 , 𝑡 ) 𝜃 ( 𝑰 ; 𝑰 , 𝑰 ) d 𝑰 d 𝑰 , (A1)where 𝑓 is the distribution function of 𝑰 ≡ ( 𝑚, 𝑉 ) (grain mass andvolume, respectively) at time 𝑡 , 𝐾 ( 𝑰 , 𝑰 ) is the collisional kernel(product of the collisional cross-section and the relative velocity ofthe two colliding grains), 𝜃 ( 𝑰 ; 𝑰 , 𝑰 ) is the distribution function ofthe produced grains by the collision, and the integration is executedfor all the relevant range of 𝑰 , which is usually [ , ∞] × [ , ∞] . Wedistinguish the two colliding grains with subscipts 1 and 2, referred toas grain 1 and 2, respectively [i.e. 𝑰 = ( 𝑚 , 𝑉 ) and 𝑰 = ( 𝑚 , 𝑉 ) ].The normalization of 𝜃 is determined by ∫ 𝑚𝜃 ( 𝑰 ; 𝑰 , 𝑰 ) d 𝑰 = 𝑚 . (A2)We note that when we consider the collision of grain 1 with grain2, we only consider the redistribution of 𝑚 (i.e. we separately treatthe collision of grain 2 with grain 1). This is why only 𝑚 enters thenormalization. Because of this, we do not have a factor of 1 / 𝑉 , that is, weintegrate it for 𝑉 . We adopt the following form for 𝜃 for simplicityand for analytical convenience: 𝜃 ( 𝑰 ; 𝑰 , 𝑰 ) = ¯ 𝜃 ( 𝑚 ; 𝑚 , 𝑚 ) 𝛿 [ 𝑉 − 𝑉 + ( 𝑚 ; 𝑰 , 𝑰 )] , (A3)where ¯ 𝜃 describes the mass distribution function of the grains formedafter the collision between grains with 𝑰 and 𝑰 , 𝛿 is Dirac’s deltafunction, and 𝑉 + describes the volume of the grain formed fromthe collision between grains 1 and 2 (and the produced grain has amass of 𝑚 ). This expression assumes that ¯ 𝜃 is independent of thevolume (or porosity) of the original grains and that the volume ofthe collisional product is determined by 𝑰 and 𝑰 and is a functionof 𝑚 . By integrating both sides of equation (A1) for 𝑉 , we obtainthe equation for the zeroth moment, ˜ 𝑛 ( 𝑚, 𝑡 ) , defined in equation (1).We also take the first moment of equation (A1); that is, we multiplyboth sides of equation (A1) by 𝑉 and integrate them for 𝑉 to obtainthe equation for ¯ 𝑉 defined by equation (2). The resulting momentequations are written as (see also O09) 𝜕 ˜ 𝑛 ( 𝑚, 𝑡 ) 𝜕𝑡 = − ˜ 𝑛 ( 𝑚, 𝑡 ) ∫ ¯ 𝐾 ( 𝑚 ; 𝑚 ) ˜ 𝑛 ( 𝑚 , 𝑡 ) d 𝑚 + ∬ ¯ 𝐾 ( 𝑚 ; 𝑚 ) ˜ 𝑛 ( 𝑚 , 𝑡 ) ˜ 𝑛 ( 𝑚 , 𝑡 ) ¯ 𝜃 ( 𝑚 ; 𝑚 , 𝑚 ) d 𝑚 d 𝑚 , (A4) MNRAS , 1–17 (2020) orosity of interstellar dust 𝜕 ¯ 𝑉 ( 𝑚, 𝑡 ) ˜ 𝑛 ( 𝑚, 𝑡 ) 𝜕𝑡 = − ˜ 𝑛 ( 𝑚, 𝑡 ) ∫ 𝑉 𝐾 ( 𝑚 ; 𝑚 ) ˜ 𝑛 ( 𝑚 , 𝑡 ) d 𝑚 + ∬ 𝑉 + 𝐾 ( 𝑚 ; 𝑚 ) ˜ 𝑛 ( 𝑚 , 𝑡 ) ˜ 𝑛 ( 𝑚 , 𝑡 ) ¯ 𝜃 ( 𝑚 ; 𝑚 , 𝑚 ) d 𝑚 d 𝑚 , (A5)where¯ 𝐾 ( 𝑚 ; 𝑚 ) ≡ ∬ 𝐾 ( 𝑚 , 𝑉 ; 𝑚 , 𝑉 ) 𝑓 ( 𝑚 , 𝑉 ) ˜ 𝑛 ( 𝑚 ) 𝑓 ( 𝑚 , 𝑉 ) ˜ 𝑛 ( 𝑚 ) d 𝑉 d 𝑉 , (A6) 𝑉 𝐾 ( 𝑚 ; 𝑚 ) ≡ ∬ 𝑉 𝐾 ( 𝑚, 𝑉 ; 𝑚 , 𝑉 ) 𝑓 ( 𝑚, 𝑉 ) ˜ 𝑛 ( 𝑚 ) 𝑓 ( 𝑚 , 𝑉 ) ˜ 𝑛 ( 𝑚 ) d 𝑉 d 𝑉 , (A7) 𝑉 + 𝐾 ( 𝑚 ; 𝑚 , 𝑚 ) ≡ ∬ 𝑉 + ( 𝑚 ; 𝑚 , 𝑉 , 𝑚 , 𝑉 ) 𝐾 ( 𝑚 , 𝑉 ; 𝑚 , 𝑉 )× 𝑓 ( 𝑚 , 𝑉 ) ˜ 𝑛 ( 𝑚 ) 𝑓 ( 𝑚 , 𝑉 ) ˜ 𝑛 ( 𝑚 ) d 𝑉 d 𝑉 . (A8)The integration is performed for [ , ∞] × [ , ∞] . The moment equa-tions, in general, are not closed since a higher-order moment alwaysappears. To close this hierarchy, we adopt the same assumption asO09; that is, the volume is replaced with the mean value at each 𝑚 . Under this assumption, the distribution function is written as 𝑓 ( 𝑚, 𝑉 ) = ˜ 𝑛 ( 𝑚 ) 𝛿 [ 𝑉 − ¯ 𝑉 ( 𝑚 )] . O09 refer to this approximation asthe volume-averaging approximation, and confirmed that it gives aconsistent result with the full solution of the 2D Smoluchowski equa-tion for their coagulation problem. Although there is no guaranteethat this approximation is valid for shattering, the increasing fillingfactor at small grain radii (Section 3.2) is at least qualitatively con-sistent with the evolution expected from the production of compactfragments by shattering.Adopting the above volume-averaging approximation, we obtain¯ 𝐾 ( 𝑚 ; 𝑚 ) = 𝐾 [ 𝑚 , ¯ 𝑉 ( 𝑚 ) ; 𝑚 , ¯ 𝑉 ( 𝑚 )] , (A9) 𝑉 𝐾 ( 𝑚 ; 𝑚 ) = ¯ 𝑉 ( 𝑚 ) ¯ 𝐾 ( 𝑚 ; 𝑚 ) , (A10) 𝑉 + 𝐾 ( 𝑚 ; 𝑚 , 𝑚 ) = 𝑉 + [ 𝑚 ; 𝑚 , ¯ 𝑉 ( 𝑚 ) , 𝑚 , ¯ 𝑉 ( 𝑚 )] ¯ 𝐾 ( 𝑚 ; 𝑚 ) . (A11)We apply these relations to equations (A4) and (A5),reorganize the notations as ¯ 𝐾 ( 𝑚 ; 𝑚 ) = 𝐾 𝑚 ,𝑚 , 𝑉 + [ 𝑚 ; 𝑚 , ¯ 𝑉 ( 𝑚 ) , 𝑚 , ¯ 𝑉 ( 𝑚 )] = ( 𝑉 + ) 𝑚𝑚 ,𝑚 . and mul-tiply both sides in equation (A4) by 𝑚 . Finally, introducing 𝜌 ( 𝑚, 𝑡 ) ≡ 𝑚 ˜ 𝑛 ( 𝑚, 𝑡 ) and 𝜓 ( 𝑚, 𝑡 ) ≡ ¯ 𝑉 ( 𝑚, 𝑡 ) ˜ 𝑛 ( 𝑚, 𝑡 ) , we obtainequations (4) and (5). This paper has been typeset from a TEX/L A TEX file prepared by the author. MNRAS000