Implications of Grain Size Distribution and Composition for the Correlation Between Dust Extinction and Emissivity
DDraft version September 28, 2020
Typeset using L A TEX twocolumn style in AASTeX62
Implications of Grain Size Distribution and Composition for the Correlation Between Dust Extinction and Emissivity
Ioana A. Zelko and Douglas P. Finkbeiner
1, 2 Harvard-Smithsonian Center for Astrophysics60 Garden StreetCambridge, MA 02138, USA Department of Physics, Harvard University17 Oxford StreetCambridge, MA 02138, USA
ABSTRACTWe study the effect of variations in dust size distribution and composition on the correlation betweenthe spectral shape of extinction (parameterized by R V ) and far-infrared dust emissivity (parameterizedby the power-law index β ). Starting from the size distribution models proposed by Weingartner &Draine (2001a), using the dust absorption and emission properties derived by Laor & Draine (1993)for carbonaceous and silicate grains, and by Li & Draine (2001) for PAH grains, we calculate theextinction and compare it with the reddening vector derived by Schlafly et al. (2016). An optimizerand an MCMC are used to explore the space of available parameters for the size distributions. Wefind that larger grains are correlated with high R V . However, this trend is not enough to explain theemission-extinction correlation observed by Schlafly et al. (2016). For the R V − β correlation to arise,we need to impose explicit priors for the carbonaceous and silicate volume priors as functions of R V .The results show that a composition with higher ratio of carbonaceous to silicate grains leads to higher R V and lower β . A relation between E (B − V) /τ and R V is apparent, with possible consequencesfor the recalibration of emission-based dust maps as a function of R V . Keywords: interstellar medium, interstellar dust, interstellar dust extinction, interstellar dust processes INTRODUCTIONDust is an important component of the interstellarmedium, forming structures in the space between thestars in our galaxy. Dust is formed through the deathprocess of stars, through supernovas or stellar winds. Itis composed of elements that have formed in the stars,and it plays an important role in the further formationof complex molecules (Draine 2011).Dust scatters and absorbs ultraviolet and visibile lightcoming from the interstellar radiation field around it.The scattering and absorption together cause extinc-tion of the ultraviolet and visible light. Past studieshave aimed to characterize the wavelength dependen-cies of the extinction. Their work (Savage & Mathis1979; Fitzpatrick 1999; Cardelli et al. 1989) used theparameter R V = A ( V ) / ( A ( B ) − A ( V )) for characteriz-ing extinction functions, based on the observation thatone parameter would capture most of the variation inextinction across the sky.This is a simplifying assumption that holds in certainwavelength regions and breaks down in the UV, where the complexity of the extinction has been shown to betoo great for a single parameter (Peek & Schiminovich2013). In this work, for the extinction we consider thewavelength range 0 . . µ m, for which one parameter issufficient to describe the variation in R V , but see § a r X i v : . [ a s t r o - ph . GA ] S e p Zelko and Finkbeiner Wavelength [ µ m]Extinction Emission Extinction BandsEmission BandsDust Optical Prop. − Frequency [THz]
Figure 1.
The extinction is evaluated at ten bands of the reddening vector from S16. The emission is evaluated at the fourbands used by the Planck satellite team: 353GHz, 545GHz, 857GHz from Planck, and 3000GHz (100 µ m) from IRAS. Theextinction and emission wavelength ranges are far apart, which raises the question of what drives the R V − β relation. and found a correlation between R V and the far in-frared dust emissivity power law, β . The emissivity datawas obtained from the Planck satellite (Planck Collab-oration et al. 2016a). We will use the reddening lawfrom S16 to constrain the dust. It is in good agreementwith the commonly used Fitzpatrick (1999) reddeningcurve, including its variation about the mean. We useS16 because it provides error bars at each of 10 wave-lengths providing an obvious way to compute a likeli-hood, whereas Fitzpatrick (1999) does not.Our goal in this analysis is to see if variations in dustgrain-size distributions and composition can explain theobserved correlation between R V and β . Weingartner &Draine (2001a) (hereafter WD01) fit the size distribu-tions using information for the volume of the grains, ex-tinction A ( λ ), and optical parameters of Laor & Draine(1993), under the assumption of spherical grains. Inaddition to this information, we take into account the R V − β relation found by S16, as well as their reddeninglaw for values ranging between 0.5 µ m and 4.5 µ m (Fig.1). As a result, we can fit a size distribution with thenew R V − β constraints and ask what drives the R V − β relation. We start from the parameters describing thesize distribution of the grains of dust. The scope is to seeif the variation in the 11 parameters of the proposed sizedistribution function can explain the observed correla-tion between R V and β . We can also explore the effectof a having dust grains exposed to different interstellarradiation field intensities.In § § §
4, and the conclusion in § MODELING − − − Radius [ µ m] − n − H a d n g r / d a ( c m ) b C b C a t , a c α, β C PAH Graphite R V = 3.1 Figure 2.
Example of a size distribution for carbonaceousgrains, showing the impact of each of the parameters in themodel: C is an overall factor, related to the abundanceof carbon atoms per hydrogen nucleus; the exponential α in the power law term ( aa t ) α can adjust the slope in d ln n d ln a for a < a t ; β can add a positive (for β >
0) or negative(for β <
0) curvature to the slope. For a > a t , the termexp (cid:8) [ − [( a − a t ) /a c ] ] (cid:9) creates an exponential cutoff whosesharpness can be controlled by a c . In addition, in the func-tion D ( a ),the sum of two log-normal size distributions forsmall radius is controlled by its amplitude b C , which repre-sents the total carbon abundance per hydrogen nucleus inthe log-normal population. Grains with radii smaller than10 − µ m are modeled as PAHs. Our goal is to determine whether variations in thesize distributions of the grains of dust can explain thecorrelation between R V and β found in S16. We usemodels for the size distributions for different types ofgrains. Using these models together with models for theabsorption and scattering cross sections of the grains, weare able to calculate the extinction. Also, together withemission cross section, and an interstellar radiation field(ISRF) we can compute an equilibrium temperature foreach size and type of grain. Using that, we can predict ust extinction-emission correlation. Properties of the Dust Grains
Dust Grain Size Distribution — We use the models forthe dust grain size distributions proposed by WD01 (forwork leading up and related to this, see also Mathis et al.(1977), Greenberg (1978), Cardelli et al. (1989), Desertet al. (1990), Li & Draine (2001), Li & Greenberg (1997)and Jones et al. (2013) for the core-mantle model dustsize distribution, and Wang et al. (2015) or the updatedversion of silicate-graphite model with the addition of apopulation of large, micron-sized dust grains). An alter-native model for the size distribution has been proposed by Zubko et al. (2004), which can be explored in a futurework.In the WD01 model, the dust is modeled using twoseparate grain populations: silicate composition, andgraphite (carbonaceous) composition. For the small car-bonaceous grains (radii smaller than 10 − µ m), differentoptical coefficients are used, corresponding to neutraland ionized polycyclic aromatic hydrocarbons (PAHs).PAHs are structures made of hexagonal rings of carbonatoms with hydrogen atoms attached to the boundary.It is assumed that neutral and ionized PAHs each givehalf of the contribution of the PAHs. Other types ofgrains (such as oxides of silicon, magnesium, and iron,carbides, etc) are not included.The size distributions are modeled by Eq. 1 and 2.1 n H d n gr ( a )d a = D ( a ) + Ca (cid:16) aa t (cid:17) α g × (cid:40) βa/a t , β ≥ − βa/a t ) − , β < (cid:41) × (cid:40) , < a < a t exp (cid:8) [ − [( a − a t ) /a c ] ] (cid:9) , a > a t (1) D ( a ) = (cid:40) , for silicate dust2 . · − bCa e − . a/ . + 9 . · − bCa e − . a/ , for carbonaceous/PAH dust (2)These equations are described by 11 parameters: 5corresponding to the silicate population and 6 to thePAH and graphite (Fig. 2). Optical Parameters of the Dust Grains — To calculate theemission and extinction of dust we need to know theoptical properties of the grains of dust, such as the ab-sorption and scattering coefficients. For silicates andgraphite, we use the values derived by Laor & Draine(1993) and Draine & Lee (1984). For PAH-carbonaceousgrains we use the properties obtained in Li & Draine(2001) .The models contain 81 log-spaced radii between10 − µ m and 10 µ m for silicates, 30 from 3 . × − µ mto 10 − µ m for PAHs, and 61 between 10 − µ m and 10 µ mfor graphite .Laor & Draine (1993) model dust grains as solidspheres of radius a with absorption cross section atwavelength λ of C abs ( λ, a ). They label the scattering ∼ draine/dust/dust.diel.html. The specific files are files Gra 81.gz,PAHion 30.gz, PAHneu 30.gz and Sil 81.gz. The file for graphite, Gra 81.gz, actually has 81 log-spacedradii between 10 − µ m and 10 µ m, but we use only the 61 between10 − µ m and 10 µ m to complement the range of radius for thePAHs. cross section at wavelength λ with C sca ( λ, a ) and theextinction cross section with C ext ( λ, a ) ≡ C abs ( λ, a ) + C sca ( λ, a ). The scattering and absorption efficiencies Q sca and Q abs are defined as: Q sca ( λ, a ) ≡ C sca ( λ, a ) πa ; Q abs ( λ, a ) ≡ C abs ( λ, a ) πa . (3)The wavelength range for the optical parameters forall types of grains is 10 − µ m to 1mm. The graphite andsilicate files have 241 log-spaced wavelength samples.For the PAH files, their wavelength array is five timesmore dense than the wavelength array from the graphiteor silicate files, so we take only every fifth value, corre-sponding to exactly the same values as the sampling ofthe graphite and silicate files.For the wavelength range between 335 µ m and1000 µ m, we model the absorption using a power law, Q sca,abs ( λ, a ) = τ ( a ) · ( λ/λ ) − θ ( a ) (Appendix B). We areinterested in looking at the absorption coefficient behav-ior for different compositions (Fig. 3). What we noticeis that carbonaceous and silicate grains show quite adifferent power law index as a function of the radius.Thus, one can expect to control the resulting emissivitypower law index β for a collection of dust grains bychanging composition or size. We use the power lawfit of the absorption optical properties to extend them Zelko and Finkbeiner − − − Radius [ µ m]1 . . . . . Q a b s P o w e r I nd e x θ GraphiteSilicatePAH neutralPAH ionized
Figure 3.
The absorption optical coefficients for carbona-ceous and silicate grains can be approximated with a powerlaw for wavelength range between 335 µ m and 1000 µ m. Thepower law index has a different dependence on radius for eachtype of grain. In a collection of grains, carbonaceous grainscan contribute higher θ s than the silicate grains, leading tolower β index for the entire collection. to 10 µ m. This range is more in line with future CMBexperiments such as PIXIE (Kogut et al. 2011).2.2. Extinction
Extinction Modeling
For a given collection of dust grains along the line ofsight, we want to calculate the extinction A , defined as: A ( λ ) = m attenuated − m = 2 . F λ F aλ (4)where F aλ , m attenuated are the dust attenuated observedflux and magnitude of the object, and F λ , m are theflux and magnitude that would have been observed ifthere would have been no attenuation from dust. Thus,extinction can be related to the optical depth τ ( λ ) by A ( λ ) = (2 . e ) τ ( λ ). The optical depth is createdfrom the contributions of each grain along the line ofsight. Let i be the index of the grain type, referringto graphite, PAHs, or silicates. For grains of radius a of type i , their impact on the optical depth can be ex-pressed as the product of an effective extinction crosssection C ext ,i ( λ, a ) and the column density N i ( a ). Then,the optical depth given by a distribution of grains of dif-ferent radii a is given by: τ ( λ ) = (cid:88) i (cid:90) d N i d a C ext ,i ( λ, a ) d a (5) Filter g r i z y λ [ µ ]m 0.503 0.6281 0.7572 0.8691 0.9636 ν [THz] 595.8 477.3 395.9 344.9 311.1Filter J H K W1 W2 λ [ µ ]m 1.2377 1.6382 2.1510 3.2950 4.4809 ν [THz] 242.2 183.0 139.4 90.98 66.90 Table 1.
Wavelength and frequency values for the ten pointswhere we compare the modeled extinction with the extinc-tion data coming from S16 reddening vector.
The fraction of dust grains per radius becomes:d N i ( a )d a = dd a (cid:90) n i ( a, s )d s = dd a (cid:90) (cid:18) n i n H (cid:19) ( a, s ) n H ( s )d s, (6)where s is the path length along the direction of integra-tion, n i [grains cm − ] is number of dust particles of type i per volume, n H [atoms cm − ] is the number of hydrogenatoms per volume, and N H [atoms cm − ] = (cid:82) n H ( s )d s isthe hydrogen column density. In this analysis we assumethe dust to gas ratio is constant along the line of sight s . As a result,d N i ( a )d a = (cid:18)(cid:90) n H ( s )d s (cid:19) n H d n i ( a )d a = N H n H d n i ( a )d a . (7)Using the fact that d a = a d log a , the optical depth canthen be calculated as: τ ( λ ) N H = π (cid:88) i (cid:90) n H d n i ( a )d a Q ext ,i ( λ, a ) a d log a . (8)The extinction A λ over the column density is: A ( λ ) N H = (2 . e ) π (cid:88) i (cid:90) n H d n i ( a )d a Q ext ,i ( λ, a ) a d log a (9)2.2.2. Extinction Data
S16 derived the dust extinction curve towards 37,000stars in different directions across the sky. Using pho-tometry from Pan-STARRS1 (Hodapp et al. 2004;Chambers et al. 2016), Two Micron All-Sky Survey(2MASS, Skrutskie et al. 2006) and the Wide-field In-frared Survey Explorer (WISE, Wright et al. 2010; Cutriet al. 2013), and spectra from the APOGEE survey (Ma-jewski et al. 2017; Eisenstein et al. 2011), they performeda principal component analysis and found that the ex-tinction function is well approximated by two principalcomponents, called the vector R (constant across thedirections in the sky) and a perturbation vector d R d x .Both R and d R d x have norm 1. The extinction functioncan be expressed as: A Schlafly = R + x d R d x , (10) ust extinction-emission correlation. x is a parameter that varies across the sky, withvalues between -0.4 and 0.4. Extinction laws are usuallycharacterized by the parameter R V = A (V) A (B) − A (V) . How-ever, since S16 don’t have access to the distances to thestar, the absolute gray component of the extinction isnot known. Instead, they approximate the R V parame-ter with R (cid:48) V = 1 . A ( g ) − A ( W A ( g ) − A ( r ) − .
18. The parameter x is related to R (cid:48) V using equation 11: R (cid:48) V = 3 . . x (11)The intent was that x = 0 ( R (cid:48) V of 3.3) corresponds toa mean reddening vector. However, this results in R (cid:48) V of 3.3 at Fitzpatrick (1999) R V of 3.1. Subsequently, inthis analysis, we use the notation R V to refer to the R (cid:48) V from S16.The reddening vector is specified at the wave-lengths/frequencies showed in Table 1. These wave-lengths/frequencies have been obtained by S16 byweighting over the M-giant star spectrum and over thebandpass of the detectors. ν mean ,b = (cid:82) νS ν F ν,b d ν (cid:82) S ν F ν,b d ν , where S ν is the M-giant spectrum, b represents the index ofthe band (g, r, i, . . . ), and F ν,b represents the filterweight.Li et al. (2014) find that the aliphatic 3.4 µ m C-Hstretch absorption band is seen in diffuse clouds, and ab-sent in dense regions. Therefore, for lines of sight with larger R V , the 3.4 µ m extinction band is weaker or evenabsent. This raises the question of whether R V -basedCardelli et al. (1989) parameterization is valid only at λ < µ m. However, within the range of E (B − V) mea-sured in S16, they do not see evidence for significantvariation in 3.4 µ m W1 and 4.6 µ m W2 Wide-field In-frared Survey Explorer Wright et al. (2010) bands. Sincein this work we employ the extinction laws of S16, the R V parameter is used.2.3. Grain equilibrium temperature as a function ofradius
To calculate the thermal radiation emitted by a collec-tion of dust grains, we need to know the temperatures ofthe grains. The grains are exposed to the ambient radia-tion field. In this calculation, we take into account onlyradiative heating and ignore the collisional heating thatwould be provided to the grains in the situation whenthey are surrounded by gas. In the case of dense clouds,however, this can become a relevant contribution.
Interstellar Radiation Field — We follow § νu ISRF ν = , hν > . . × − erg cm − ( hν/ eV) − . , . < hν < . . × − erg cm − ( hν/ eV) − , . < hν < . . × − erg cm − ( hν/ eV) . , . < hν < . πν/c ) (cid:88) i =1 w i B ν ( T i ) , hν < .
04 eV (12)where w = 1 × − , w = 1 . × − , w = 4 × − ,and T = 7500K, T = 4000K, T = 3000K. In ouranalysis, we would want to modify the radiation field toaccount for inhomogeneities in the interstellar medium,where we can have areas that are hotter than others. Forthat, as an approximation, we will multiply the radiationfield by a factor χ ISRF that varies from 0.5 to 2.
The thermal equilibrium equation — For each grain radius a , we assume thermal equilibrium between the absorbedradiation and emitted radiation ( P in = P out , Fig. 5) Weassume the grain is spherical and emits like a black bodyof unknown temperature T , which we aim to determine.The absorbed radiation is assumed to come from the interstellar radiation field surrounding the grain sphereuniformly.The thermal equilibrium equation for one dust grainof size a is thus: (cid:90) ∞ Q abs ( λ, a ) πa χ ISRF u ISRF ( λ )d λ = (cid:90) ∞ Q abs ( λ, a ) πa πc B λ ( T )d λ (13)The integral in Eq. 13 is taken over the wavelengthrange from 10 − µ m to 10mm, using the extension shownin Appendix B. We obtain the equilibrium temperaturesfor 4 types of grains for different radii (Fig. 6). Zelko and Finkbeiner − Wavelength [ µ m] − u i s r f × − [ J / µ m ] Figure 4.
Mean interstellar radiation field from Mezgeret al. (1982) and Mathis et al. (1983), for χ ISRF = 1. − Wavelength [ µ m]0 . . . λ I λ [ W / µ m ] × − Radius 0.01 µ m gra-Poutgra-Pinsil-Poutsil-Pin Figure 5.
Power input and output for a single grain of ra-dius 0.01 µ m. λI λ is the power per log λ . ie. the area underthe curve is the power. This shows that the power absorbedequals the power emitted, for the case of equilibrium tem-perature. This method of calculation for the equilibrium tem-perature assumes that at each radius for each type ofgrain there is a single temperature. This approximationbreaks down as the radius of the grain becomes smallenough. A grain stays at an equilibrium temperatureif no one photon it absorbs or emits carries enough en-ergy to perturb the temperature much. Big grains havea thermal energy much larger than one photon. Butsmall grains do not: a single photon with several eVcarries more energy than the entire thermal energy ofthe grain, and the emission and absorption of a singlequanta can create temperature spikes. In our calcula-tion, we are considering grains as small as 3 . − µ m (like the PAHs), ina future study, there can be a benefit from replacing the 10 − − − Radius [ µ m]10203040 E q u ili b r i u m T e m p e r a t u r e [ K ] GraphiteSilicatePAH neutralPAH ionized
Figure 6.
The calculated equilibrium temperatures for the4 types of grains for different values of radii, for χ ISRF = 1. approximation with a different method where one con-siders a temperature distribution for each radius size, asdone in the work of Draine & Li (2001). However, inthis work we are mainly interested in long wavelengthemission where (cid:104) B ν ( T ) (cid:105) = B ν ( (cid:104) T (cid:105) ).In addition to the grain size and type, the variation ofthe ISRF throughout the ISM leads to variations in tem-perature. To reproduce this effect, we allow the ISRFmultiplier parameter χ ISRF to vary between 0.5 and 2,and calculate the equilibrium temperature for the range(Fig. 7). 2.4.
Modeling Emission
Calculating the emission intensity from a collectionof grains
The emissivity (power radiated per unit volume perunit frequency per unit solid angle) coming from a col-lection of grains is defined as: j ν = (cid:88) i (cid:90) d a d n i d a C abs ,i ( ν, a ) B ν ( T eq ( a )) (14)where i = index for carbonaceous, silicate and PAHgrains. The spectral intensity I ν is defined as the emis-sivity integrated along the line of sight s : I ν = (cid:90) j ν d s = (cid:90) (cid:18) j ν n H (cid:19) ( s ) n H ( s )d s (15)Since n i n H is assumed to be constant along the line ofsight, j ν n H becomes constant along the line of sight aswell. Using N H = (cid:82) n H ( s )d s , we obtain: I ν = j ν n H (cid:90) n H ( s )d s = j ν n H N H == N H (cid:88) i (cid:90) d a n H d n i d a C abs ,i ( ν, a ) B ν ( T eq ( a )) (16) ust extinction-emission correlation. − . . . log χ ISRF ISRF multiplier . . . . l og T [ K ] E q u ili b r i u m T e m p e r a t u r e Gra. 1.0e-01 µ mSil. 1.0e-01 µ mGra. 1.0e+01 µ m Sil. 1.0e+01 µ mPAH n. 1.0e-03 µ mPAH i. 1.0e-03 µ m − − − Radii [ µ m] . . . . . ( + θ ) - S l o p e GraphiteSilicatePAH n.PAH i.
Figure 7.
Top panel: log of the equilibrium temperature (T)as a function of the log of the multiplier of the interstellarradiation field ( χ ISRF ), for carbonaceous, silicate and PAH(neutral and ionized) grains at select radii. They follow alinear dependence. The temperature curves are calculatedat fixed values for the radius of the grains, as specified inthe legend. Bottom panel: the difference between 1 / (4 + θ )( θ refers to the power law indexes of the absorption opticalcoefficients as seen in Fig. 3) and the slopes of the linear fitsversus grain radii, for each type of grain. The difference issignficant to warrant not using the 1 / (4 + θ ) approximation,and it serves as a good check for our calculations. The Modified Black Body Fit
We aim to compare our analysis with the 2013 Planckrelease (Planck Collaboration et al. 2014b) , as was usedby S16 . The spectral index data can be found in the fileHFI CompMap ThermalDustModel 2048 R1.20.fits at https://irsa.ipac.caltech.edu/data/Planck/release 1/all-sky-maps/previews/HFI CompMap ThermalDustModel 2048 R1.20/index.html. We select the directions in the sky to reproducethe same analysis done by S16, whose data is available athttps://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:
Spectral energy density (SED) of emission from dusthas been modeled in practice by the Planck Collabora-tion et al. (2014b) using a modified black body (MBB)function: I ν = τ (cid:0) ν
353 GHz (cid:1) β B ν ( T ) , (17)where B ν ( T ) is the Planck function for dust of temper-ature T . 353GHz is chosen as a reference frequency.The assumption is used in the optically thin limit. Thepower law parameterized by β models the dependenceof the emission cross-section with frequency. The fit forthe three parameters in Equation 17 is performed usingdata from four photometric bands: 353GHz, 545GHz,857GHz from Planck, and 3000GHz (100 µ m) from IRAS(Schlegel et al. 1998; Beichman et al. 1988). Becausethese are the bandpasses the Planck Collaboration et al.(2014b) used in their analysis, to compare to theirresults, we evaluate the intensity at the same fourbandpasses. We use the weighting given in AppendixB/Table 1 of Planck Collaboration et al. (2014b), andthe corresponding response functions . METHODSThe goal of this work is to explore the space of WD01grain size distributions to find those that are consistentwith our prior knowledge about dust, including:1. the shape of the reddening curve and its variationwith R V ,2. the amount of reddening per H atom, and3. the abundance of metals (C, Si, etc.) per H atomrequired to make dust.For each sample from the WD01 parameter space, wecompute the emission spectrum expected for dust in areference radiation field, and fit the τ , β , and T param-eters of a modified blackbody (MBB) as described in § R V and β parameters. The Planck filter files can be found on the web-site http://pla.esac.esa.int/, in the section called ”Software,Beams, and Instrument Model”. At the time of thispaper, Planck has 3 releases HFI RIMO R1.10.fits (2013),HFI RIMO R2.00.fits (2015), HFI RIMO R3.00.fits (2016). Weuse HFI RIMO R1.10 fits because it was the one used forthe data release from the Planck Collaboration et al. (2014b).For IRAS, the filter files can be found on the websitehttps://irsa.ipac.caltech.edu/IRASdocs/exp.sup/ch2/tabC5.html
Zelko and Finkbeiner
Parameter 10 b C α g β g a t,g [ µ m] a c,g [ µ m] 10 V g [cm H − ] α s β s a t,s [ µ m] a c,s [ µ m] 10 V s [cm H − ]Lower Boundary 3.0 -3.0 -30. 0.000355 0.000355 max b C -3.0 -30. 0.001 0.001 0.1Upper Boundary 5.0-6.5 -0.5 30. 10.00000 10.00000 6 × × Table 2.
The boundaries for the parameter space explored by the MCMC. Since a t controls the position of the exponentialdrop, we allow it to values over the range of the grains. a c controls the smoothness of the exponential factor, so it should beable to get values of similar magnitude to a t . As such, we give it the same range. Gaussian priors are given for the carbonaceous( V g =PAH+graphite) and silicate ( V s ) volumes that are centered within the range. The V g parameter’s low bound is set to belarge enough to account for the maximum possible contribution coming from Eq. 2 with the highest allowed b C value. Themaximum values for V g and V s are set to be 6 times the reference values in WD01. The limits on the b C parameter are explainedin § To perform our analysis, we create a reference extinc-tion function, constrain the gray component of the ex-tinction, normalize to physical values of extinction percolumn density of hydrogen ( A ( I ) /N H ), impose physicalboundaries on the parameters, and keep the total massper H atom of the grain distributions within an expectedrange.An Markov chain Monte Carlo (MCMC)-like algo-rithm is used to explore the available space for the sizedistribution parameters. However, MCMCs are not anoptimization technique, and they can take a long timeto converge on the points of high likelihood. To help thesolver converge, an optimizer is run to obtain an initialguess for the MCMC. A possible downside is that by do-ing this, a region of the parameter space might remainunexplored.Using the posterior points from the MCMC analysiswe create the emissivity corresponding to the Planckbandpass filters and perform a modified black body fitas described in section § R V and β parameters. Creating the reference extinction functions — S16 focusedon the shape of the reddening curve by using the rel-ative extinction in 10 bands. Their work does not in-form us about the gray component (which could not bemeasured without knowing the distance to each star)or the overall amplitude per H atom (which they didnot measure). Therefore our first step is to estab-lish a ”target” reddening curve by setting an additiveand multiplicative term to fix these degrees of free-dom. The condition that A(H)/A(K) = 1.55 (Indebe-touw et al. 2005) constrains the additive term as fol- lows: we define A (cid:48) Schlafly = A Schlafly + C Schlafly with A Schlafly = R + x d R d x A Schlafly ( H ) + C Schlafly A Schlafly ( K ) + C Schlafly = 1 .
55 = rA Schlafly ( H ) + C Schlafly = A Schlafly ( K ) r + C Schlafly rC Schlafly = A Schlafly ( H ) − A Schlafly ( K ) rr − A (cid:48) Schlafly = A Schlafly + C Schlafly (18)Fixing the grey component creates degeneracy between A ( H ) and A ( K ). To maintain the correct number ofdegrees of freedom, we remove the H band from the A vectors and therefore from the covariance matrix andthe ∆ χ calculation. A ( H ) is determined by the otherparameters and is no longer independent, so it can beignored in the calculation and recovered at the end.Having fixed the additive term, we now impose anextinction per N(H) assumption to fix the multiplica-tive term. Cardelli et al. (1989) suggested the conven-tion that A ( I ) /N H = 2 . × − cm . To be consis-tent with the extinction functions presented in WD01,we define A ( I ) /N H = 3 . × − cm . We denotethis quantity by C AIN H , and the normalized extinctionby A (cid:48)(cid:48) Schlafly( λ ) N H = A (cid:48) Schlafly ( λ ) A (cid:48) Schlafly ( I ) × C AIN H . C AIN H is a convention,not a measurement with an error. In reality its valuemost likely varies across the sky. If a future experimentmakes a different measurement of A I /N H , it should betaken into account.Thus, the reference extinction function can be con-structed using: A reference( λ ) N H = A (cid:48)(cid:48) Schlafly( λ ) N H = A (cid:48) Schlafly ( λ ) A (cid:48) Schlafly ( I ) × C AIN H (19) Volume of the dust grains — As we let the MCMC and theoptimizer explore the parameter space, we want to makesure the size distribution doesn’t require more atoms ust extinction-emission correlation. V tot ,g ≈ . × − cm H − , and for silicates V tot ,s ≈ . × − cm H − . The mean of the Gaussian in theprior is fixed for now but later in § b C represents the overall amplitude of the bumps; tocalculate the PAH volumes, one needs to also add partcoming from the non- D ( a ) part of the size distribution,integrated over the range of the PAH radii.Since C is an overall factor, it can be calculated froma proposed combination of volume V and parameters α, β, a t , a c . As a result, we replace the C parameter witha volume parameter V , and calculate the corresponding C when needed for the size distribution calculations.Thus, the 11 parameters explored by the MCMC are b C , α g , β g , a t,g , a c,g , V g , α s , β s , a t,s , a c,s , V s . Table2 lists the boundaries for each parameter. Likelihood — We want to compare to Schlafly’s redden-ing vector. For each proposed set of 11 parameter theMCMC makes, the volume parameters are transformedinto the size distribution parameters (converting from aset of b C , α g , β g , a t,g , a c,g , V g , α s , β s , a t,s , a c,s , V s toa set of b C , α g , β g , a t,g , a c,g , C g , α s , β s , a t,s , a c,s , C s ).Using the size distribution of the grains, the resultingextinction vector A /N H is calculated at the nine wave-lengths (after H was removed when fixing the greycomponent) from S16 using equation 9.Appendix A shows the calculation for the error inextinction that gives the covariance matrix (Σ A (cid:48)(cid:48) , Eq.A11) for each value of x , based on the errors in the red-dening vectors obtained by S16.The extinction vector A /N H is compared to the ref-erence extinction vector obtained with Equation 19: A residual N H = A N H − A reference N H (20)The likelihood function is ln L = − ∆ χ , with the χ given by: ∆ χ = A T residual N H (Σ A (cid:48)(cid:48) ) − A residual N H (21) To the likelihood, we add the Gaussian prior on thevolume:ln prior = − (cid:18) V g − V g, reference . · V g, reference (cid:19) − (cid:18) V s − V s, reference . · V s, reference (cid:19) (22)to obtain the posterior:ln posterior = ln prior + ln L (23)Here V g represents the sum of the volumes for the PAHand carbonaceous grains, and V s the silicate grains.3.1. Exploring the Dust Parameters PosteriorDistribution with an MCMC
We sample from the posterior (Eq. 23) for a target ex-tinction curve at fixed R V ( § R V linearly spaced between 2.94 and 3.67. The posterioris conditional on R V instead of letting R V float, so thatthe uncertainty in the target extinction curve does notdepend on the parameters at each step in the Markovchain. To expedite burn in, we initiate the MCMC at aset of dust grain size distribution parameters determinedby optimization.The MCMC uses the ptemcee Vousden et al. (2016)package that uses Parallel Tempering. This allows for amuch more efficient exploration of the parameter spacethan something like the Metropolis-Hastings algorithm.We experimented with a number of temperatures be-tween 3 and 5, and found that there was no significantdifference in the results, so we settled for 3 temperaturesto reduce the computational time. 300 walkers were runfor each R V , for 100000 steps.3.2. Studying the Correlation Between Dust Emissivityand Absorption
Taking the final posterior distributions from all of thechains, and using the precomputed values of the tem-perature T for each radius of the grain (Fig. 6), weintegrate to calculate the specific intensity, at each ofthe four bandpass frequencies of Planck.The emitted radiation is modeled as a modified blackbody shown using Eq. 17, and fit to find the three pa-rameters ( τ , β, T ) corresponding to each sample fromthe MCMC. We generate the emission at the four wave-length bands corresponding to the Planck satellite, andthen fit the four data points to the modified black bodylaw, using corresponding weighting and bandpass filtersas used by the Planck team. The code can be found at the python repository at https://pypi.org/project/ptemcee/ or at Will Vousden github repositoryat https://github.com/willvousden/ptemcee Zelko and Finkbeiner . . . β . . . . R V -26.90-26.85-26.80-26.75-26.70 l og10 ( V c a r b ) . . . β . . . . R V -26.55-26.50-26.45-26.40 l og10 ( V s il ) . . . β . . . . R V -27.40-27.35-27.30-27.25-27.20-27.15-27.10-27.05-27.00 l og10 ( V P AH ) (a) (b) (c)Figure 8. R V vs. β for 15 MCMC runs, each corresponding to a distinct Rv value. The points are color coded by the logof the volume of (a) the carbonaceous grains, (b) the silicate grains, and (c) the PAH grains. The volume priors are fixed forcarbonaceous grains (graphite+PAH) and silicates to values used in WD01 ( § R V - β correlation. This motivates introducing a dependence of the volume priors on R V . α V sil α V c a r b - . - . - . - . - . - . - . - . - . . . . . . . . . . . . . . . . Smoothed ∆ β and χ for Rv 3.118 α V sil α V c a r b - . - . - . - . - . - . - . - . - . - . . . . . . . . . . . . . . Smoothed ∆ β and χ for Rv 3.209 α V sil α V c a r b - . - . - . - . - . - . . . . . . . . . . . . . . . . Smoothed ∆ β and χ for Rv 3.300 α V sil α V c a r b - . - . - . - . - . - . . . . . . . . . . . . . . . . . . Smoothed ∆ β and χ for Rv 3.391 α V sil α V c a r b - . - . - . . . . . . . . . . . . . . . . . . . Smoothed ∆ β and χ for Rv 3.482 α V sil α V c a r b - . - . - . - . - . . . . . . . . . . . . . . . . . . . . . Smoothed ∆ β and χ for Rv 3.573 Figure 9.
Minimum χ (red-orange-yellow dashed-dotted contours ) and ∆ β (gray and blue continuous/dashed contours) for6 values of R V , as a function of carbonaceous and silicate volume coefficients. The ∆ β = 0 contour (blue) corresponds to theS16 empirical R V - β correlation (Eq. 24). In each panel the minimum χ for ∆ β = 0 is marked by a blue dot. The location ofthese dots shows a monotonic trend as a function of R V . Next, we calculate R V for sample, and see if there isany correlation between β and R V , thus comparing tothe results in S16 (Fig. 12). RESULTS AND DISCUSSION 4.1.
Correlation between dust extinction andfar-infrared emissivity
We consider 2 hypotheses for the origin of the R V − β correlation: ust extinction-emission correlation. R V and β to variations in grainsize distribution, holding the volume (per H) ofeach species fixed .II. the composition hypothesis, which requires varia-tion in the relative volumes of silicates and car-bonaceous grains to vary as a function of R V .Fixing the volume priors in Hypothesis I proved tobe too restrictive: the MCMC does not explore the fullrange of the β parameters (Fig. 8). This fails to yieldthe observed R V − β correlation.If an R V -dependent prior on volumes is needed forhypothesis II, what form should it take?4.1.1. Optimizer Results
To generate a hypothesis about what form of volumepriors might give rise to the R V − β correlation, we beginby mapping out the volume parameter space V s , V g .We use an optimizer to explore the effects of the vol-ume priors on the resulting β . This framework is usedinstead of the MCMC in the beginning of the analysis inorder to take advantage of the great increase in speed,which allows us to make an exploration of the parame-ter space of volume priors on a much finer grid. Thus,we can get an idea fast of what parameter combinationswould provide useful results.The Gaussian volume priors are defined to be centeredat V sil = α V sil · V ,V sil , V carb = α V carb · V ,V carb where α V sil and α V carb are control parameters that we use toscale the total carbonaceous and silicate volumes up anddown. V ,V sil and V ,V carb are the reference fiducial valuesfor the volumes from WD01, as described in section 3.We let α V sil and α V carb take values between 0.35 and4.2, and sample the interval logarithmically at 50 values,creating a 50 by 50 grid for each R V .For each of the points in the grid, the χ returnedby the optimizer is calculated. We smooth over the re-sulting image using a Gaussian filter with σ = 1 .
5, andcalculate the contour plots over the resulting image.For each R V panel, to compare with the expected β ,we use a best fit line obtained from the R V vs. β data However, the Kramers-Kronig relation can be used to deter-mine a lower bound on the volume of the grains for a given extinc-tion function integrated over a finite wavelength interval (Purcell1969). Mishra & Li (2017) applied this relation to approximate thevolumes for the silicate and carbonaceous grains. To perform thiscalculation, we would have to integrate over the UV part of theabsorption, which is a bit uncertain (Peek & Schiminovich 2013),so it was not included in this work. Nevertheless, one should keepin mind that the bounds that we are proposing might be violatingthis relation slightly. . . . . . . . R V × − × − × × α V c a r b o n a c e o u s s ili c a t e α Vsil α Vcarb α Vsil fit α Vcarb fit Figure 10. α V is the ratio of the volume relative to Draineprior. Using the data from the optimizer, we fit ln α V aslinear functions of R V , obtaining a decreasing function forsilicates and an increasing function for carbonaceous grains. set from S16, given by eq 24. β Schlafly ( R V ) = 2 . − . · R V (24)Thus, for each β obtained from the points in the grid,we calculate ∆ β = β − β Schlafly . Using a Gaussian filterwith σ = 1 .
5, we smooth over the resulting image, andcalculate the contour plots. The resulting contour plotsfor both the χ and ∆ β analysis are superimposed (Fig.9). Each grid corresponds to a different R V value.For each grid, we find the point of minimum χ thathas ∆ β = 0. We read the corresponding α V carb and α V sil values, and plot them against R V (Fig. 10). The pointsshow a log-linear dependence. Performing linear fits ofln α V carb vs. R V and ln α V sil vs. R V , the functions shownin equations 25 and 26 are obtained.ln α V sil ( R V ) = − . · R V + 5 .
32 (25)ln α V carb ( R V ) = 1 . · R V − .
85 (26)This relation of the volume priors on R V also de-pends on extinction per N(H), assumed to be A I /N H =3 . × − cm (beginning of § R V - β relationitself does not depend on this assumption, because nei-ther R V nor β depend on the column density per se.However an increase in the assumed A I /N H would re-quire an increase in the volume (per H) of each species,by the same factor. In other words, a different A I /N H convention would simply slide the contours in Fig. 9 upand to the right. If a different A I /N H is chosen, suchas A I /N H = 2 . × − cm (Zhu et al. 2017), the vol-ume results can in turn be scaled by 0.77 = 2.6/3.38 andobtain the same behaviour.2 Zelko and Finkbeiner x R V bC α g β g a t,g [ µm ] a c,g [ µm ] 10 V g [cm H − ] α s β s a t,s [ µm ] a c,s [ µm ] 10 V s [cm H − ]-0.040 2.94 3.00 -0.66 0.07 0.00035 0.68 0.813 -1.83 -5.41 0.16 0.19 8.355-0.034 2.99 3.00 -0.72 0.12 0.06709 0.73 0.871 -1.60 -23.45 0.13 0.21 8.385-0.029 3.04 3.00 -0.83 -0.00 0.00035 1.02 0.975 -1.59 -20.81 0.22 0.11 7.460-0.023 3.09 3.00 -1.47 -0.00 0.00282 1.27 1.047 -1.56 -19.47 0.17 0.16 7.364-0.017 3.14 3.00 -1.05 -0.17 0.00883 1.34 1.111 -1.94 -2.78 0.21 0.14 7.187-0.011 3.20 3.00 -1.49 -0.00 0.00035 1.72 1.204 -1.65 -5.28 0.16 0.18 6.662-0.006 3.25 4.70 -1.30 -0.01 0.00119 1.41 1.315 -1.76 -8.79 0.35 0.00 6.4960.000 3.30 4.66 -1.30 -3.08 0.00036 1.54 1.466 -1.78 -2.72 0.26 0.09 5.6830.006 3.35 4.11 -2.50 0.00 0.00141 1.53 1.526 -1.23 -6.84 0.14 0.20 5.1760.011 3.40 5.40 -1.93 -0.01 0.00241 2.58 1.672 -1.05 -30.00 0.16 0.17 4.9440.017 3.46 5.81 -2.56 0.00 0.00035 1.78 1.868 -0.80 -20.84 0.18 0.14 4.1760.023 3.51 5.94 -2.01 -0.00 0.00065 4.83 1.991 -0.52 -19.04 0.09 0.22 4.1260.029 3.56 6.09 -1.67 -0.53 0.00041 2.48 2.110 -0.61 -4.43 0.15 0.17 3.6840.034 3.61 6.07 -2.85 0.02 0.01434 1.85 2.303 -0.87 -10.25 0.33 0.02 3.6540.040 3.66 6.41 -1.78 -0.25 0.00047 3.68 2.406 -0.50 -3.92 0.19 0.14 3.325 Table 3.
Optimized parameter values of the dust grain size distributions for each x / R (cid:48) V , for hypothesis II. − . − . − . − . − . − . − . ln posterior . . . . . . . . . R v Figure 11.
Log of the posterior for each point in the 15 runsof the MCMC corresponding to a distinct R V value. Whilethere is variation in the posterior values, this variation iscontained in the range between ln posterior of -7.5 and -20.The modest range in ln P is reassuring. MCMC Analysis
The fits described in Section 4.1.1 represent a hypoth-esis for how carbonaceous and silicates volumes might . . . . . Planck β . . . . . . . . . . R ’ ( V ) Schlafly et al.MCMC
Figure 12.
An MCMC analysis if performed for 15 valuesof R V , using dedicated volume priors as a function of R V ( § R V and β are obtained, seen here superimposed on the data fromS16. The anti-correlation relation trend between R V and β is thus reproduced. In this work R V refers to the R (cid:48) V fromS16. change as a function of R V , and now we must validatethat hypothesis using an MCMC analysis. ust extinction-emission correlation. . . . β . . . . R V -27.40-27.30-27.20-27.10-27.00-26.90-26.80-26.70 l og10 ( V c a r b ) . . . β . . . . R V -26.50-26.40-26.30-26.20-26.10-26.00 l og10 ( V s il ) . . . β . . . . R V -27.50-27.40-27.30-27.20-27.10 l og10 ( V P AH ) (a) (b) (c) . . . β . . . . R V V s il / V c a r b . . . β . . . . R V V P AH / V c a r b (d) (e)Figure 13. The points in the MCMC cloud are here colorcoded by (a) the log 10 of the volume of the carbonaceous grains (b) the log 10 of the volume of the silicate grains (c) the log 10 of the volume of the PAH grains (d) the ratio of the volumes of silicategrains and the volume of carbonaceous grains (e) the ratio of the volumes of PAH grains and the volume of carbonaceous grains.We find that a composition with higher ratio of carbonaceous to silicate grains leads to more R V and lower β . Carbonaceousgrains are the sum of PAH and graphite grains. PAH and graphite also increase with R V independently. The resulting values for each parameter can be seenin Table 3. Using the functions found in equations 25and 26, we now set up volume priors correspondingly,and turn to performing an MCMC analysis as describedin the § R V values linearly spaced between 2.936and 3.664 ( x varies between -0.04 and 0.04). For eachof the MCMC runs, we calculated the dust emissivityand the R V for each sample from the posterior, usingthe procedure described in § β is spread between 1.4 and 1.8. The variation inspectral index value is significant, which indicates thathaving different size distributions of dust grains in dif-ferent directions of the sky can motivate the need tomodel these different lines of sight with different spec-tral index values. The values of R V and β obtained arein the same range as the ones obtained by S16. Most im- portantly, our values reproduce the trend of the R V − β anti-correlation. The log posterior of all the end po-sitions of the chains in the runs is contained within amodest range (Fig. 11). The systematic uncertainty inthe Planck β measurements is thought to be of order0.05. Volume and Composition — as expected for hypothesis II,we find that a composition with higher ratio of carbona-ceous to silicate grains leads to higher R V and lower β (Fig. 13).While the functions 25 and 26 seem very precise, theydon’t represent unique solutions. We are making a plau-sibility argument, not a final determination of model pa-rameters for a firmly established model. One might beable to find other solutions that explain the R V − β cor-relation. But it is suggestive that R V -dependent volumepriors give this behavior, and fixed priors do not. Size of grain distribution — One of the intuitive expecta-tions of the analysis was that larger grains would lead to4
Zelko and Finkbeiner − − − Radius [ µ m] . . . . . . C D F a d n g r / d a / n H Carbonaceous
Rv 2.94Rv 3.30Rv 3.66 − − − Radius [ µ m] . . . . . . C D F a d n g r / d a / n H Silicates
Rv 2.94Rv 3.30Rv 3.66
Figure 14.
The CDF corresponding to the volume of the grains. For both silicates and carbonaceous grains, we can see thatas one moves to higher R V , at least 50% of the volume is in grains of larger and larger size. This is in accordance with theexpecatation that larger grains lead to higher R V . For the carbonaceous grains, we represented the CDF separately for PAHgrains at radii smaller than 0.01 µ m and graphite grains at radii larger than 0.01 µ m. For the PAHs at low R V of 2.94, the sizedistribution is constrained tightly through the b C parameter, which results in reduced variation represented by the very thinblack line. − − − Radius [ µ m] − − − n − H a d n g r / d a ( c m ) Carbonaceous
Rv 2.94Rv 3.30Rv 3.66 − − − Radius [ µ m] − n − H a d n g r / d a ( c m ) Silicates
Rv 2.94Rv 3.30Rv 3.66
Figure 15.
The size distributions corresponding the points in the posterior clouds of from 3 MCMCs at different R V s. Thetwo bumps in the carbonaceous grains at radii smaller than 0.01 µ m come from the constraints imposed on the minimum valueof the b C parameter that informs the amount of PAH. We see that larger R V leads to larger grains cutoffs, but also to a largerratio of carbonaceous to silicate grains. higher R V . In order to test this hypothesis, we plot thecumulative density function of the volume of the grainsversus radii. Fig. 14 shows what percentage of the vol-ume of the grains is made up of radii smaller than eachpossible radius value. For example, for the case of sil-icates, 80% of low R V volume is in grains with radiussmaller than 0.1 µ m, but 20% of high R V volume. Forboth silicates and carbonaceous grains, as R V increases,at least 50% of the volume is in grains of larger andlarger size. The size distributions coming from the posterior re-sulting from the MCMC are calculated (Fig. 15). Theyreproduce acceptable size distributions as proposed byWD01 (Fig. 2).There is a broad range of parameters that can produceeach R V , but the distributions are largely distinct fromeach other as R V changes (Fig. 16). Each parameterhas a different impact on the R V − β anti-correlation(Fig.17). Further work may be warranted to isolate the ust extinction-emission correlation. b C − − − α g −
20 0 20 β g − − log 10 a t,g − − log 10 a c,g − . − . − . − . log 10 C g − − − α s −
20 0 20 β s . . . . a t,s . . . . a c,s − − − log 10 C s Figure 16.
Binned histograms of each of the 11 parameters of the size distributions, shown for 3 different R V MCMC runs.The distributions for each of the 11 parameters change substantially as R V changes, and are largely distinct from each other. effect of each parameter on R V independently of theothers. Ultraviolet Extinction — The model is constrained tomatch the S16 extinction curve in 10 bands, but is not6
Zelko and Finkbeiner . . . . . . . R V b C 3 . . . . . . . . . . . . . α g − . − . − . − . . . . . . . . β g − − . . . . . . . R V log a t,g − . − . − . − . − . . . . . . . . . . log a c,g − . − . − . − . − . . . . . . . . . . log C g − − − − − . . . . . . . R V α s − . − . − . − . . . . . . . . β s − − . . . β . . . . a t,s . . . . . . . . . . β . . . . R V a c,s . . . . . .
30 1 . . . β . . . . log C s − − − − − − − Figure 17.
The points from the MCMC clouds colorcoded using the values of the 11 size distribution parameters. a c,s and a t,s are anti-correlated; this anti-correlation can explain the spread of β values at fixed R V . If we hold the ratio of carbonaceous tosilicates, ratio of PAH to graphite and R V constant, we can push β back and forth by about 0.05 (a relatively small amount)by changing the large grains cutoff. ust extinction-emission correlation. − − − /λ [ µ m − ] − − − − − A ( λ ) /3 . × − c m N H R V V V /λ [ µ m − ] A ( λ ) /3 . × − c m N H R V V V Figure 18.
The left panel shows the extinction functions obtained from the MCMC for 3 values of R V (2.993, 3.300,3.664)for wavelengths ranging from far-infrared to x-rays. Looking over the entire wavelength range, we can see that the extinctionfunctions resulting from the MCMC are in agreement with the reference function from WD01. The right panel is a zoomed-inview on the UV feature at 2175 ˚A; the feature varies between the Milky Way, Small Magellanic Galaxy, Large Magellanic GalaxyFitzpatrick & Massa (2007), with the feature being prominent in the Milky way and less prominent outside of it. The 2175˚Abump is often associated with PAHs (Joblin et al. (1992), Li & Draine (2001), Mishra & Li (2015)) and in the WD01 modelits amplitude is explicitly controlled by the b C parameter. Since in our study we were aiming to replicate the conditions in theMilky Way, we restricted the b C parameter to make sure we have a minimum of PAHs involved. sufficient to constrain the 2175˚A feature. We compareextinction functions derived from the MCMC with ref-erence functions from WD01 (calculated for the fourthrow of Table 1 of WD01) and find good agreementacross a wide wavelength range (Fig. 18). In partic-ular the 2175˚A feature varies between the Milky Way,Small Magellanic Galaxy, Large Magellanic Galaxy Fitz-patrick & Massa (2007), with the feature being promi-nent in the Milky way and less prominent outside ofit. The 2175˚A bump is often associated with PAHs(Joblin et al. (1992), Li & Draine (2001), Mishra & Li(2015)) and in the WD01 model its amplitude is explic-itly controlled by the b C parameter. Since in our studywe were aiming to replicate the conditions in the MilkyWay, we restricted b C to values greater than 3 × − tomake sure we have a minimum amount of PAH involved.The maximum value of b C is set between (5.0, 6.5) as R V goes from low to high since the total volume ofthe carbonaceous grains increases with R V as well. Theformula used is max b C = 6 × (1 + 0 . × tanh (cid:0) α V g − (cid:1) ).4.2. Effect of the Interstellar Radiation Field
The effect of the interstellar radiation field on ouranalysis is very significant. As described in Section § . . . . . Planck β . . . . . . . . . R ’ ( V ) χ ISRF
Figure 19.
The emissivity of the collection of dust particlesfor each run is calculated using a modified ISRF multiplica-tion factor χ ISRF . χ ISRF takes 10 log spaced values between0.5 and 2. This results in the R V - β correlation being shiftedleft and right relatively uniformly across R V , with a changeof up to 0.1 in β at low R V and 0.15 at high R V . For each χ ISRF value, the lines on the plot were generated from the15 average R V and β value from the MCMC posteriors. well as the cloud thickness, influence the local ISRF dustgrains are exposed to. To try to get a glimpse of what8 Zelko and Finkbeiner − . − . − .
25 0 .
00 0 .
25 0 .
50 0 . ln χ ISRF . . . . . . . l n T M BB ( K ) R V, slope
Figure 20.
The logarithm of the average T MBB , for eachof the 10 values of χ ISRF , for each of the 15 MCMC runs.This results in a linear relationship between ln T MBB andln χ ISRF . The legend shows the slope of the linear fit foreach R V value. The slope decreases with R V , a result due tothe variable dust composition, as explained in Section § the effects of changing the ISRF look like, we modifiedthe ISRF multiplication factor χ ISRF over a log spacedarray with values between 0.5 and 2 (Figs. 19 and 20).However, the effect is not that drastic, as a factor of 4in χ ISRF leads to a change of up to 0.1 in β at low R V and 0.15 at high R V .In addition, th effect of the variation of the ISRF onthe modified black body temperature fit is explored, av-eraged for each of the 15 MCMC posterior points (Fig.20). For each R V value, a linear relationship is obtainedbetween ln T MBB and ln χ ISRF . The slope of this linearfunction decreases as R V increases. One could expectthe slope to be close to 1 / (4 + β ), and due to the R V - β anti-correlation relation discussed in this paper, the av-erage value of β decreases as R V increases. This wouldresult in the slope increasing at higher R V . However, thediscrepancy is explained by the fact that we are doing amodified black body fit over a collection of grains withvariable dust composition. For each grain individually(Figs. 3 and 7), the relationship 1 / (4 + θ ) is recoveredwhen calculating the equilibrium temperatures at dif-ferent χ ISRF values. However, the value of the opticalpower law index θ varies with the type and size of grain(Fig. 3). At high R V , our runs have a higher ratio ofcarbonaceous to silicate than at low R V . Since carbona-ceous grains have higher θ than silicate grains, it resultsin a lower slope obtained when integrating over the con- tributions from all the particles in the distributions toobtain the modified black body fit.4.3. Effect on A( λ )/ τ Emission-based interstellar dust maps such as Schlegelet al. (1998) have been a very valuable tool for predict-ing extinction across the sky. They make the assumptionthat the ratio of Near Infrared extinction to the emis-sion optical depth does not vary with R V . Using theresults for our 15 MCMC runs, we calculate to ratio ofA( λ )/ τ to see how it varies with R V (Figs. 21 and22).We find that there can be a lot of variation inA( λ )/ τ . For the K band (Fig. 21), the best fitpower law for Hypothesis II is given byln ( A K /τ ) = 4 .
00 ln R V + 4 . . (27)For E (B − V) /τ (Fig. 22) the power law index issmaller, given byln ( E (B − V) /τ ) = 1 .
72 ln R V + 7 . , (28)and there is a tendency for it to be 20% lower at 2.9 and20% higher at 3.7, compared to R V =3.3.We took the E (B − V) of the stars studied by S16,and the corresponding τ data from Planck Collabo-ration et al. (2016a), and obtained the average value of µ E (B − V) /τ = 10501. Figure 22 shows that both hy-pothesis are above this value. Hypothesis II is closer andthus preferred, but this indicates that the overall dustmodel could be improved.This is a potentially impactful result that can moti-vate future research into the effect of R V variation onemission-based interstellar dust map calibrations. Thefact that sign of the trend changes depending on whetherwe assume R V variation is caused only by size variation,or by size and composition variation, suggests that addi-tional research will be required before we can confidentlyderive extinction from emission-based dust maps as Rvvaries.4.4. Correlation between spectral index andtemperature
Researchers have been looking at the T - β correlationand there isn’t quite a consensus on rather it is real (Du-pac et al. 2003; D´esert et al. 2008) or a subtle statisticalartifact (Shetty et al. 2009; Kelly et al. 2012). To ad-dress this question we consider the modified black bodyfit to models at varying R V but with fixed radiation field( χ ISRF = 1). It is important to mention that we do nothave noise in our model, thus we cannot comment on ust extinction-emission correlation. . . . . . . . . . R V A K / τ Hypothesis I . . . . . . . . . R V A K / τ Hypothesis II
Figure 21. R V variation can lead to a significant change in A( λ )/ τ . For the K band shown here, for hypothesis I the powerlaw fit is A K /τ ∝ R . , and for hypthesis II it is A K /τ ∝ R . . In this work R V refers to the R (cid:48) V from S16. studies that presented results on that. In addition, theoptical properties of the dust grains can change with thetemperature of the grains, and this can lead to a changein the spectral index. However, we do not account forthis effect in our models, so we cannot probe this.The intensity emission of the collection of dust parti-cles is fit using a MBB function, as described in Section2.4.2, for each point in the posterior resulting from theMCMC. When plotting the resulting distribution of thespectral index ( β ) as a function of the collective modifiedblack body teperature ( T ) we observe an anti-correlation(Fig. 23).One explanation for this interesting result is that T does not affect the long wavelength end very much, so ifwe boost those points (eg: with more cold dust), to getthe R V correct, the fit must lower β , but then that lower β means a higher T to get the peak in roughly the right place. As a result, there is a natural inverse correlationbetween T and β built into the problem.At high R V one has larger grains, and larger grainstend to have lower equilibrium temperatures (Fig. 6).However, our runs also have different dust composition.Carbonaceous grains have higher temperature than thesilicates, and the ratio of carbonaceous grains to sili-cate increases as R V increases. This explains why inour results, higher R V means higher temperature forthe modified black body fit over the collection of dustgrains. CONCLUSIONWe started from the size distribution of the grainsof dust, and we varied the parameters for it using anMCMC while adding the constraints from the redden-ing vector obtained by S16. We used the MCMC resultsto generate the dust emission spectrum for each sam-0
Zelko and Finkbeiner . . . . . . . . . R V E ( B - V ) / τ Hypothesis I . . . . . . . . . R V E ( B - V ) / τ Hypothesis II
Figure 22.
For hypothesis II, a power law fit gives us a variation of the form E (B − V) /τ ∝ R . . The 1.72 power lawresults in values approximately 20% lower at 2.9 and 20% higher at 3.7, compared to R V =3.3. This can have implications forcalibrations of emission-based interstellar dust maps. In comparison, for hypothesis I, the power law is E (B − V) /τ ∝ R − . .The changing sign of the power index between the case when we keep the volume fixed or varied strengthens the argument forfurther research. ple from size distribution space. We suspected that the R V − β correlation would arise naturally from the sizedistribution alone, but in the case of fixed total volumepriors for each dust species, variation of R V does notproduce an appreciable correlation. This is an interest-ing outcome, and to follow up in the search for a fullexplanation we force the R V − β and see what parame-ters give it. We find again that larger grains are corre-lated with high R V , but in addition we find an explicitfunction of carbonaceous and silicates volume priors asfunctions of R V that gives the R V − β correlation andsatisfy the constraints WD01 used. The properties ofthe optical absorption coefficients for carbonaceous and silicates offer an explanation for the results of the anal-ysis; carbonaceous grains have optical properties thatlower the β for a collection of dust grains, while silicatesraise it (Fig. 3).Widely used dust maps like SFD (Schlegel et al. 1998)and Planck Collaboration et al. (2016b) assume the ratioof E (B − V) /τ is constant. In Section § E (B − V) /τ and A( λ )/ τ on R V . Thisdependence is a testable consequence of our understand-ing of the R V − β relation in the context of the WD01models. Other optical models and size distribution pa-rameterizations are possible, but if this dependence on R V persists in future models, it would have serious con- ust extinction-emission correlation. . . . . . Sp ec t r a l I nd e x ( β ) Figure 23.
The intensity emission of the collection of dust particles is fit using a modified black body function, for each pointin the posterior resulting from the MCMC. Anti-correlation between β and T is observed from the resulting distribution of thespectral index ( β ) as a function of the collective modified black body teperature ( T ). The values of the temperature are spreadbetween 18-22K. Greater R V leads to higher temperature. This is in agreement to the fact that higher R V runs have a higherratio of carbonaceous to silicate grains, since carbonaceous grains have higher individual equilibrium temperature compared tosilicates (Fig. 6). sequences for the recalibration of emission-based dustmaps as a function of R V .Moreover, this result might provide some guidanceon how to improve these dust models in the future. E (B − V) /τ can become an additional constraint usedduring modeling. Reproducing the correct function of E (B − V) /τ versus R V based on real data could be agood target for the next studies.Modeling the size distribution and composition of dustis an area of active research. The parameterizations ofthe size distributions and the optical parameters of thegrains can be revisited. An alternative model for the sizedistribution and optical parameters has been proposedby Zubko et al. (2004), which can be explored in a fu-ture work. In the future, we might have to explore grainsthat are a combination of both carbon and silicate. Themodel we are using here, though it reproduces many em-pirical facts about dust, is necessarily a simplification ofnature. Future work may involve other materials, com-plex grain geometries, composites and coatings, etc. Ourwork is intended as a plausibility argument, not a finaldetermination of parameters for a truly complete model.One might be able to find other solutions that explainthe R V − β correlation. The robust effect we observe is that a composition with higher ratio of carbonaceous tosilicate grains leads to more R V and lower β . It is anopen question whether this tendency is a generic prop-erty of all dust models or if it is a specific feature of theprecise dust models we are using.The fact that larger R V corresponds to smaller silicatevolume can be difficult to understand. Denser regionsthat have larger R V are expected to have depleted Si,Mg and Fe from the gas phase. However, in the denseclouds it may not be possible to know how much hydro-gen there is. Since we perform our calculations per N H ,this could play a significant role. Also, if the carbon iscoming out of the gas faster than the silicate is, theremight be more carbon per N H in the dust cloud. Car-bonaceous grains could also be misidentified with grainscoated with carbons.This work depends critically on the S16 reddening vec-tors. Their pre-Gaia (Gaia Collaboration et al. 2016)analysis of the reddening law was performed in absenceof information regarding distances to the stars whoseextinction they were modeling. As a result, the abso-lute extinction cannot be determined, only the relativedifference of extinction between bands, after the graycomponent had been removed from the analysis. This2 Zelko and Finkbeiner can be improved in future work when the gray compo-nent to the extinction can be fixed using Gaia measure-ments. In addition, this analysis is performed is at fixed A (I) /N (H) , so we still have a free parameter left. If fu-ture data constrains A (I) /N (H) as a function of R V − β ,we can modify the volume relations and get the simi-lar results again with different functions, as A (I) /N (H)scale linearly with the b c , C s , C g parameters combined.The results of this study provide a possible explana-tion of the observed R V − β correlation in the context ofthe WD01, Laor & Draine (1993), Draine & Lee (1984),Li & Draine (2001) family of models. Although this ex-planation may not be unique, it increases our confidencethat the R V − β correlation can be used to our advan-tage. For example, the relation can be use as a crosscheck for CMB experiments: one can start from a sen-sitive map of the sky in R V , like one created from thedatasets from LSST (LSST Science Collaboration et al.2009), and determine the corresponding β . Conversely,one can make predictions of R V given precise measure-ments in β . The R V − β correlation provides valuableinformation about the size distribution and compositionof interstellar dust grains, and may lead us toward amore complete model of the interstellar medium. Acknowledgments — We acknowledge helpful conversa-tions with Ana Bonaca, Blakesley Burkhart, Tansu Day-lan, Bruce Draine, Cora Dvorkin, Daniel Eisenstein,John Kovac, Albert Lee, Karin ¨Oberg, Stephen Por-tillo, Eddie Schlafly, Zachary Slepian, Josh Speagle, JunYin and Catherine Zucker. IZ is supported by the Har-vard College Observatory. DF is partially supportedby NSF grant AST-1614941, Exploring the Galaxy: 3-Dimensional Structure and Stellar Streams. This re-search made use of the NASA Astrophysics Data Sys-tem Bibliographic Services (ADS), the color blindnesspalette by Martin Krzywinski & Jonathan Corum , andthe Color Vision Deficiency PDF Viewer by MarrieChatfield . Facilities:
Odyssey Cluster, Harvard University
Software: ptemcee (Vousden et al. 2016), emcee(Foreman-Mackey et al. 2013), NumPy (van der Waltet al. 2011), Matplotlib (Hunter 2007), pandas McKin-ney (2010), scikit-learn (Pedregosa et al. 2012), IPython(Perez & Granger 2007), Python (Millman & Aivazis 2011;Oliphant 2007)APPENDIX A. ERROR IN EXTINCTIONWe calculate the errors in the extinction function to be used as a reference in the MCMC (Eq. 19). Denote with (cid:15) and (cid:15) the error vectors in R and d R d x , respectively. The values for these error vectors are given in Table 2 of S16.Firstly, we calculate the error propagation in the extinction formula, given by A ( λ ) = R ( λ ) + x d R ( λ )d x .Let us assume we have a function y expressed as a linear combination of the variables x l : y ( x ) = (cid:88) l a l x l (A1)Let Σ x be the covariance matrix for the parameters x l , such that Σ x kl = E [( x l − µ x l )( x k − µ x k )]. The mean (firstmoment) of y is then given by equation A2: E [ y ] = µ y = E (cid:34)(cid:88) l a l x l (cid:35) = (cid:88) l a l E [ x l ] = (cid:88) l a l µ x l , (A2) http://mkweb.bcgsc.ca/biovis2012/color-blindness-palette.png https://mariechatfield.com/simple-pdf-viewer/ ust extinction-emission correlation. σ y = E (cid:2) ( y − µ y ) (cid:3) = E (cid:32)(cid:88) l a l x l − (cid:88) l a l µ x l (cid:33) = E (cid:32)(cid:88) l a l ( x l − µ x l ) (cid:33) = n (cid:88) k n (cid:88) l a k a l E [( x k − µ x k )( x l − µ x l )]= n (cid:88) k n (cid:88) l a k a l Σ x kl (A3)Let Σ A be the covariance matrix of A ( λ ). For our case, the coefficient vector is (1 , x ), and the variable vector x = (cid:16) R ( λ ) , d R ( λ )d x (cid:17) . S16 make the assumption that there is no covariance between the errors in the vectors R and d R d x . As a result, only the diagonal terms of the covariance matrix are non-zero. Then the error σ A ( λ ) is given byequation A4: σ A ( λ ) = (cid:15) ( λ ) + x (cid:15) ( λ ) . (A4)S16 also assume there is no covariance between the errors at different wavelength. As a result, the covariance matrixΣ A is given by equation A5: Σ A kl = δ kl σ A ( l ) (A5)Secondly, we calculate the error introduced by fixing the gray component A ( λ ) → A ( λ )+ C such that A (cid:48) ( H ) /A (cid:48) ( K ) =1 .
55 = r (Indebetouw et al. 2005): A (cid:48) ( λ ) = A ( λ ) + CC = 1 r − A ( H ) − rr − A ( K ) A (cid:48) ( λ ) = A ( λ ) + 1 r − A ( H ) − rr − A ( K ) (A6)We calculate the error in A (cid:48) ( λ ). Since the errors in the different bands from the reddening vector were not correlated,neither are the errors of A ( λ ). The r parameter comes from a measurement given with 10% precision. We keep r fixed in this analysis. The effect of having r values with ±
10% difference can be easily explored by fixing r to differentvalues and re-running the analysis.The covariance matrix of A (cid:48) ( λ ) is given by equation A7:Σ A (cid:48) kl = E [( A (cid:48) ( k ) − µ A (cid:48) ( k ) )( A (cid:48) ( l ) − µ A (cid:48) ( l ) )]= Σ A kl + 1 r − (cid:0) Σ A kH + Σ A lH (cid:1) − rr − (cid:0) Σ A kK + Σ A lK (cid:1) + 1( r − Σ A HH + (cid:18) rr − (cid:19) Σ A KK (A7)Due to fixing the gray component, the covariance matrix of A (cid:48) ( λ ) is not diagonal, but it is still symmetrical.4 Zelko and Finkbeiner
The indexes of the matrix run over the 10 filters, in the order g,r,i,z,y,J,H,K,W1,W2. Labeling with b = σ A ( H ) + r σ A ( K ) ( r − , c = rσ A ( H ) + r σ A ( K ) ( r − , d = σ A ( H ) + rσ A ( K ) ( r − :Σ A (cid:48) = σ A ( g ) + b b b b b b c d b bb σ A ( r ) + b b b b b c d b bb b σ A ( i ) + b b b b c d b bb b b σ A ( z ) + b b b c d b bb b b b σ A ( y ) + b b c d b bb b b b b σ A ( J ) + b c d b bc c c c c c r ( σ A ( H ) + σ A ( K ) ) ( r − r ( σ A ( H ) + σ A ( K ) ) ( r − c cd d d d d d r ( σ A ( H ) + σ A ( K ) ) ( r − σ A ( H ) + σ A ( K ) ( r − d db b b b b b c d σ A ( W + b bb b b b b b c d b σ A ( W + b (A8)Due to fixing of the grey component by requiring A ( H ) /A ( K ) = r , the rows (and columns) of H and K in thecovariance matrix are now related by the constant r , making the covariance matrix singular. As such, we remove therow and column corresponding to the H band, and redefine the matrix in equation A9:Σ A (cid:48) = σ A ( g ) + b b b b b b d b bb σ A ( r ) + b b b b b d b bb b σ A ( i ) + b b b b d b bb b b σ A ( z ) + b b b d b bb b b b σ A ( y ) + b b d b bb b b b b σ A ( J ) + b d b bd d d d d d σ A ( H ) + σ A ( K ) ( r − d db b b b b b d σ A ( W + b bb b b b b b d b σ A ( W + b (A9)The third step is to normalize at the third band-pass value of 7572 ˚A = 0 . µ m, corresponding to the I filter.This is done in order to fix the extinction per hydrogen column density ( N H ) to the chosen prior convention value atthe I band. As a result, we want to divide by the average of A (cid:48) ( I ), µ A (cid:48) ( I ) , and multiply by our chosen prior conventionvalue of C AIN H : A (cid:48)(cid:48) ( λ ) = A (cid:48) ( λ ) µ A (cid:48) ( I ) C AIN H (A10)Since we treat the average of A (cid:48) ( I ) as a fixed quantity without errors, the covariance matrix for the elements of thevector A (cid:48)(cid:48) ( λ ) is given by Equation A11: Σ A (cid:48)(cid:48) kl = Σ A (cid:48) kl µ A (cid:48) ( I ) C AIN H (A11) B. EXTENDING THE DUST OPTICAL PROPERTIESThe optical properties of the dust grains were extended beyond the wavelengths given by Laor & Draine (1993),Draine & Lee (1984), and Li & Draine (2001), to values between 10 µ m and 10 µ m. The extension was done bymodeling Q radii ( λ ) = τ ( λ/λ ) − θ , with λ = 1mm and fitting for the θ and τ parameters for each radii, using the last20 bins. Then, the optical parameters were calculated for the new range. The calculation was performed separatelyfor the scattering and absorption coefficients (Figs. 24, 25, 26 and 27). ust extinction-emission correlation. − − − λ [ µ m]10 − − − − − − − G r a ph i t e O p t i c a l P r o p e r t i e s case Q abs; radii [ µ m] case Q sca; radii [ µ m] Figure 24.
Extending the graphite optical properties for absorption and scattering between 10 µ m and 10 µ m. The graydotted line indicates the boundary of the extension. − − − λ [ µ m]10 − − − − − − S ili c a t e O p t i c a l P r o p e r t i e s case Q abs; radii [ µ m] case Q sca; radii [ µ m] Figure 25.
Extending the silicate optical properties for absorption and scattering between 10 µ m and 10 µ m. The gray dottedline indicates the boundary of the extension. REFERENCES6
Zelko and Finkbeiner − − − λ [ µ m]10 − − − − − − N e u t r a l P AH O p t i c a l P r o p e r t i e s case Q abs; radii [ µ m] case Q sca; radii [ µ m] Figure 26.
Extending the neutral PAH optical properties for absorption and scattering between 10 µ m and 10 µ m. The graydotted line indicates the boundary of the extension. − − − λ [ µ m]10 − − − − − − I o n i ze d P AH O p t i c a l P r o p e r t i e s case Q abs; radii [ µ m] case Q sca; radii [ µ m] Figure 27.
Extending the ionized PAH optical properties for absorption and scattering between 10 µ m and 10 µ m. The graydotted line indicates the boundary of the extension.Beichman, C. A., Neugebauer, G., Habing, H. J., Clegg,P. E., & Chester, T. J. 1988, Infrared astronomicalsatellite (IRAS) catalogs and atlases. Volume 1:Explanatory supplement, 1.https://ui.adsabs.harvard.edu/abs/1988iras....1.....B Bennett, C. L., Hill, R. S., Hinshaw, G., et al. 2003, TheAstrophysical Journal Supplement Series, Volume 148,Issue 1, pp. 97-117., 148, 97, doi: 10.1086/377252 ust extinction-emission correlation. Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989,Astrophysical Journal v.345, p.245, 345, 245,doi: 10.1086/167900Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016,The Pan-STARRS1 Surveys.https://ui.adsabs.harvard.edu/abs/2016arXiv161205560CCutri, R. M., Wright, E. L., Conrow, T., et al. 2013,Explanatory Supplement to the AllWISE Data ReleaseProducts, Tech. rep.https://ui.adsabs.harvard.edu/abs/2013wise.rept....1CDesert, F. X., Boulanger, F., & Puget, J. L. 1990,Astronomy and Astrophysics, Vol. 500, p. 313-334 (2009),500, 313. https://ui.adsabs.harvard.edu/abs/1990A { & } A...237..215DD´esert, F. X., Mac´ıas-P´erez, J. F., Mayet, F., et al. 2008,Astronomy and Astrophysics, 481, 411,doi: 10.1051/0004-6361:20078701Draine, B. T. 2011, Physics of the Interstellar andIntergalactic Medium.https://ui.adsabs.harvard.edu/abs/2011piim.book.....DDraine, B. T., & Lee, H. M. 1984, Astrophysical Journalv.285, p.89, 285, 89, doi: 10.1086/162480Draine, B. T., & Li, A. 2001, The Astrophysical Journal,551, 807, doi: 10.1086/320227Dupac, X., Bernard, J. P., Boudet, N., et al. 2003,Astronomy and Astrophysics, v.404, p.L11-L15 (2003),404, L11, doi: 10.1051/0004-6361:20030575Eisenstein, D. J., Weinberg, D. H., Agol, E., et al. 2011,The Astronomical Journal, Volume 142, Issue 3, articleid. 72, 24 pp. (2011)., 142, 72,doi: 10.1088/0004-6256/142/3/72Finkbeiner, D. P., Davis, M., & Schlegel, D. J. 1999, TheAstrophysical Journal, Volume 524, Issue 2, pp. 867-886.,524, 867, doi: 10.1086/307852Fitzpatrick, E. L. 1999, The Publications of theAstronomical Society of the Pacific, Volume 111, Issue755, pp. 63-75., 111, 63, doi: 10.1086/316293Fitzpatrick, E. L., & Massa, D. 2007, The AstrophysicalJournal, 663, 320, doi: 10.1086/518158Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman,J. 2013, Publications of the Astronomical Society of thePacific, Volume 125, Issue 925, pp. 306 (2013)., 125, 306,doi: 10.1086/670067Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al.2016, Astronomy and Astrophysics, 595, A1,doi: 10.1051/0004-6361/201629272Greenberg, J. M. 1978, in Cosmic dust. (A78-38226 16-90)Chichester, Sussex, England and New York,Wiley-Interscience, 1978, 187–294.https://ui.adsabs.harvard.edu/abs/1978codu.book..187G Hodapp, K. W., Kaiser, N., Aussel, H., et al. 2004,Astronomische Nachrichten, Vol.325, Issue 6, p.636-642,325, 636, doi: 10.1002/asna.200410300Hunter, J. D. 2007, Computing in Science and Engineering,9, 90, doi: 10.1109/MCSE.2007.55Indebetouw, R., Mathis, J. S., Babler, B. L., et al. 2005,The Astrophysical Journal, Volume 619, Issue 2, pp.931-938., 619, 931, doi: 10.1086/426679Joblin, C., Leger, A., & Martin, P. 1992, The AstrophysicalJournal, 393, L79, doi: 10.1086/186456Jones, A. P., Fanciullo, L., K¨ohler, M., et al. 2013,Astronomy and Astrophysics, 558, A62,doi: 10.1051/0004-6361/201321686Kelly, B. C., Shetty, R., Stutz, A. M., et al. 2012, TheAstrophysical Journal, Volume 752, Issue 1, article id. 55,17 pp. (2012)., 752, 55, doi: 10.1088/0004-637X/752/1/55Kogut, A., Fixsen, D. J., Chuss, D. T., et al. 2011, Journalof Cosmology and Astroparticle Physics, Issue 07, id. 025(2011)., 2011, 25, doi: 10.1088/1475-7516/2011/07/025Laor, A., & Draine, B. T. 1993, Astrophysical Journalv.402, p.441, 402, 441, doi: 10.1086/172149Li, A., & Draine, B. T. 2001, The Astrophysical Journal,Volume 554, Issue 2, pp. 778-802., 554, 778,doi: 10.1086/323147Li, A., & Greenberg, J. M. 1997, Astronomy andAstrophysics, 323, 566. https://ui.adsabs.harvard.edu/abs/1997A { & } A...323..566LLi, Q., Liang, S. L., & Li, A. 2014, Monthly Notices of theRoyal Astronomical Society, 440, L56,doi: 10.1093/mnrasl/slu021LSST Science Collaboration, Abell, P. A., Allison, J., et al.2009, LSST Science Book, Version 2.0.https://ui.adsabs.harvard.edu/abs/2009arXiv0912.0201LMajewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al.2017, The Astronomical Journal, Volume 154, Issue 3,article id. 94, 46 pp. (2017)., 154, 94,doi: 10.3847/1538-3881/aa784dMathis, J. S., Mezger, P. G., & Panagia, N. 1983,Astronomy and Astrophysics, Vol. 500, p. 259-276 (2009),500, 259. https://ui.adsabs.harvard.edu/abs/1983A { & } A...128..212MMathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977,Astrophysical Journal, Vol. 217, p. 425-433 (1977), 217,425, doi: 10.1086/155591McKinney, W. 2010, in { P } roceedings of the 9th { P } ythonin { S } cience { C } onference, ed. S. van der Walt &J. Millman, 56–61, doi: 10.25080/Majora-92bf1922-00a Zelko and Finkbeiner
Mezger, P. G., Mathis, J. S., & Panagia, N. 1982,Astronomy and Astrophysics, vol. 105, no. 2, Jan. 1982,p. 372-388., 105, 372. https://ui.adsabs.harvard.edu/abs/1982A { & } A...105..372MMillman, K. J., & Aivazis, M. 2011, Computing in Scienceand Engineering, 13, 9, doi: 10.1109/MCSE.2011.36Mishra, A., & Li, A. 2015, The Astrophysical Journal, 809,120, doi: 10.1088/0004-637X/809/2/120—. 2017, The Astrophysical Journal, 850, 138,doi: 10.3847/1538-4357/aa937aOliphant, T. E. 2007, Computing in Science andEngineering, 9, 10, doi: 10.1109/MCSE.2007.58Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2012,Scikit-learn: Machine Learning in Python.https://ui.adsabs.harvard.edu/abs/2012arXiv1201.0490PPeek, J. E. G., & Schiminovich, D. 2013, The AstrophysicalJournal, 771, 68, doi: 10.1088/0004-637X/771/1/68Perez, F., & Granger, B. E. 2007, Computing in Scienceand Engineering, 9, 21, doi: 10.1109/MCSE.2007.53Planck Collaboration, Ade, P. A. R., Aghanim, N., et al.2014a, Astronomy and Astrophysics, 571, A30,doi: 10.1051/0004-6361/201322093Planck Collaboration, Abergel, A., Ade, P. A. R., et al.2014b, Astronomy and Astrophysics, Volume 571, id.A11,37 pp., 571, A11, doi: 10.1051/0004-6361/201323195Planck Collaboration, Adam, R., Ade, P. A. R., et al.2016a, Astronomy and Astrophysics, 594, A10,doi: 10.1051/0004-6361/201525967Planck Collaboration, Ade, P. A. R., Aghanim, N., et al.2016b, Astronomy and Astrophysics, Volume 586,id.A132, 26 pp., 586, A132,doi: 10.1051/0004-6361/201424945Reach, W. T., Dwek, E., Fixsen, D. J., et al. 1995,Astrophysical Journal v.451, p.188, 451, 188,doi: 10.1086/176210Savage, B. D., & Mathis, J. S. 1979, Annual Review ofAstronomy and Astrophysics, 17, 73,doi: 10.1146/annurev.aa.17.090179.000445 Schlafly, E. F., Meisner, A. M., Stutz, A. M., et al. 2016, \\