ACCURATE MODELING OF X-RAY EXTINCTION BY INTERSTELLAR GRAINS

The absorption and scattering of light by interstellar dust and gas has been of interest to astronomers for over a century. Evidence for dust existing between stars in the Milky Way first appeared in the astronomical literature when Herschel (1785) described gaps in the density of stars across the sky. The notion that interstellar dust was responsible for the dimming of starlight appears to have been first proposed by Struve (1847) and independently by Pickering (1897), Kapteyn (1904, 1909a, 1909b), and Barnard (1907, 1910).

Since the early 20th century, our understanding of the dust in the interstellar medium (ISM) has continued to evolve. Extensive studies of the wavelength-dependence of extinction (i.e., reddening) established that submicron grains are present in the ISM. Schalén (1938) estimated a characteristic radius of $\sim 0.05\;\mu {\rm{m}}$, and detailed calculations by Oort & van de Hulst (1946) later put the characteristic radius of ISM dust grains at $\sim 0.3\;\mu {\rm{m}}$. In 1949, Hall (1949) and Hiltner (1949a, 1949b) discovered the polarization of starlight. Their discovery was the first evidence that ISM dust grains are (1) non-spherical, and (2) coherently aligned over large distance scales. In recent decades, the study of interstellar grains has used observations of emission, absorption, and scattering, ranging from microwaves to X-rays.

Interstellar grains absorb and scatter X-rays. Both gas and dust attenuate X-rays propagating through the ISM, but elemental X-ray absorption edges differ between atoms, ions, and solids, so that near-edge X-ray absorption fine structure (NEXAFS) can in principle reveal the composition of interstellar grains (Martin 1970; Martin & Sciama 1970; Evans 1986; Woo 1995; Forrey et al. 1998; Draine 2003; Lee & Ravel 2005; Lee et al. 2009; Costantini et al. 2012; Pinto et al. 2013).

X-ray halos surrounding astrophysical sources provide another valuable tool with which to study the ISM (Hayakawa 1970; Martin 1970). Small angle scattering of X-rays by dust grains along the line of sight produces a halo around the source (Overbeck 1965; Hayakawa 1970; Martin 1970). Measurements of these X-ray halos can be used to test and constrain dust models (see e.g., Smith et al. 2006). Recent work by Seward & Smith (2013) used Chandra observations of Cyg X-1 to look for azimuthal asymmetry in the surrounding X-ray halo, a technique that could potentially be used to constrain dust shape and grain alignment. Observations of X-ray halos around variable sources can constrain the orientation and geometry of dust clouds, allowing astronomers to study the ISM in three-dimensions (Predehl et al. 2000; Vaughan et al. 2004; Tiengo et al. 2010; Heinz et al. 2015), and could even be used for extragalactic distance determination (Draine & Bond 2004).

The focus of this paper will be on the importance of accounting for X-ray scattering when inferring abundances of different grain materials from absorption edge measurements. The Wilms et al. (2000) model for X-ray attenuation, which is employed by XSPEC (Arnaud 1996), Spex (Kaastra et al. 1996), and several other X-ray data analysis codes, does not include scattering; this, together with an approximate treatment of absorption in large grains that assumes each grain is a slab of thickness $4a/3$, can result in incorrect conclusions regarding the composition and abundance of interstellar dust.

The Rayleigh–Gans (RG) approximation (Mauche & Gorenstein 1986), often used to model X-ray halos, is also prone to errant application. As shown in Smith & Dwek (1998), for an MRN (Mathis et al. 1977) size distribution of spherical graphite and silicate grains, the RG approximation substantially overestimates the intensity of the soft X-ray ($\lesssim 1$ keV) scattering halo. Anomalous diffraction theory (ADT), by contrast, is accurate and easy to use.

The paper is organized as follows: in Section 2.1, the problem of modeling dust extinction is discussed, along with popular theoretical and computational techniques for solving the scattering problem. Section 2.2 explains the perils of ignoring dust scattering when modeling observations of X-ray extinction, especially when inferring elemental abundances and mineralogy from X-ray absorption edges. In Section 3, the advantages of using ADT to model dust scattering and absorption are discussed. An open source code suite, "GGADT," which uses ADT to calculate scattering and absorption by grains of arbitrary composition and geometry is presented in Section 3.2. Section 5 summarizes the salient points discussed in this paper. More details about GGADT are given in Appendices A–D.

2.1. Overview of Theoretical and Numerical Techniques

For spheres, Mie theory provides a truncated multipole expansion of the full solution to the scattering problem. As grain radius a increases, the number of multipole terms required also increases. For grain sizes that are much larger than the incident wavelength λ, Mie theory becomes computationally demanding, and computer implementations of Mie theory are limited by roundoff error. Codes are available that can handle size parameters $x\equiv 2\pi a/\lambda$ as large as $\sim {10}^{4}$ (Wiscombe 1980). However, $x=5070(a/\mu {\rm{m}})(E/{\rm{keV}})$ can exceed this limit for large grains at high energies. Mie theory is limited to spheres, but the polarization of starlight indicates that interstellar grains are not spherical (Hall 1949; Hiltner 1949b; Heiles 2000) and thus other methods are needed to model the extinction of light by interstellar dust.

Several approximation schemes have been used to calculate the scattering and absorption of light by non-spherical grains. Among the more popular methods are: the electric dipole approximation (Rayleigh scattering), the Rayleigh–Gans approximation, ADT, and the discrete dipole approximation (DDA). For a given material composition, each of these approximations is valid for a range of grain sizes and electromagnetic wavelengths. The domains of validity for astrosilicate grains are shown in Figure 1.

Zoom In Zoom Out Reset image size

The Rayleigh–Gans approximation is often used in X-ray astronomy to model X-ray halos (see, e.g., Vaughan et al. 2004). The Rayleigh–Gans approximation assumes that each infinitesimal volume element of the grain responds only to the incident electric field. The total scattered field produced by the entire grain is computed by integrating over the dipole scattering contributions from all volume elements.

The Rayleigh–Gans approximation assumes (1) $| m-1| \ll 1$, where m is the complex index of refraction of the grain material, and (2) $2{ka}| m-1| \ll 1$, negligible complex phase shifts in the incident wave as it travels through the grain, where $k\equiv 2\pi /\lambda$. For astrosilicate dust and X-ray wavelengths, these assumptions hold only for small grains ($a\lesssim 0.02\;\mu {\rm{m}}\cdot (E/{\rm{keV}})$).

The DDA (Purcell & Pennypacker 1973; Draine & Flatau 1994) discretizes the grain into a number of finite volume elements, and can be used to calculate scattering and absorption by grains with arbitrary geometries. The finite volume elements must be small enough that they can be treated as dipoles. The DDA, unlike the Rayleigh–Gans approximation, does not assume that the volume elements are non-interacting. The DDA is constrained by computational requirements to problems with $| m| {ka}\lesssim 30$, which limits it to $a\lesssim 0.006\;\mu {\rm{m}}\cdot ({\rm{keV}}/E)$, hence the DDA is not useful at X-ray energies.

ADT (van de Hulst 1957) is applicable to grains of arbitrary geometry that are large compared to the incident wavelength. The approximations used in the derivation of ADT require that $| m-1| \ll 1$, and ${ka}\gg 1$. Draine & Allaf-Akbari (2006) showed that the validity conditions for ADT ($| m-1| \lesssim 0.1$ and ${ka}\gtrsim 10$) are satisfied for silicate grains when $E\gtrsim 60$ eV and $a\gtrsim 0.035\;\mu {\rm{m}}\times \left(\frac{60\;\mathrm{eV}}{E}\right)$. The computational requirements are modest, and it is readily extended to arbitrary geometries.

A number of other approximations exist to efficiently compute scattering and absorption in various limiting cases. For optically "soft" ($| m-1| \ll 1$) particles, Sharma & Somerford (2006) provides a comprehensive comparison of many available approximations, and their accuracy.

2.2. Some Widely Used Models of X-Ray Attenuation by Dust Grains

One popular but flawed technique that has been used in the X-ray astronomy literature is to ignore scattering contributions to the dust extinction. The tbvarabs routine in the XSPEC package¹ and the dabs routine in the Spex package² both use the ISM model of Wilms et al. (2000) (hereafter WAM2000), which approximates the extinction cross section of a single grain as the sum of atomic photoionization cross sections with an approximate correction for grain self-absorption, and integrates over the dust size distribution given by Mathis et al. (1977) (MRN):

$$\begin{eqnarray}{\sigma }_{{\rm{ext}}}^{{\rm{WAM2000}}}(E) & = & \displaystyle \sum _{i}{\displaystyle \int }_{{a}_{{\rm{min}}}}^{{a}_{{\rm{max}}}}{\left(\displaystyle \frac{1}{{n}_{{\rm{H}}}}\displaystyle \frac{{{dn}}_{i}}{{da}}\right)}_{{\rm{MRN}}}\pi {a}^{2}\\ & & \times \left(1-{e}^{-\frac{4}{3}{\alpha }_{i}(E)a}\right){da},\end{eqnarray} \tag{ 1 }$$

where i refers to different grain materials and ${\alpha }_{i}(E)=(4\pi /\lambda ){\rm{Im}}({m}_{i})$ is the attenuation coefficient for grain material i with complex refractive index m_i(E).

The WAM2000 attenuation model ignores scattering contributions to the dust extinction and further approximates a grain of radius a as having a uniform thickness $(4/3)a$. For sufficiently small grains, a grain of volume V has an absorption cross section ${C}_{{\rm{abs}}}(E)\approx \alpha (E)V$ and scattering is unimportant, as assumed by the WAM2000 dust attenuation model. For larger grains, however, scattering contributes significantly to the extinction, and should be taken into account when modeling attenuation.

Figure 2 shows the contributions from scattering and absorption to the extinction cross section of a Weingartner & Draine (2001, hereafter WD01) grain size distribution. Figure 3 shows the ratio of the true extinction cross section and the cross section obtained from the WAM2000 approximation. Figures 2 and 3 show that scattering is an important contributor to grain extinction at X-ray energies for realistic models of ISM grains.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To illustrate the importance of including scattering when computing extinction by interstellar dust grains, we use the WAM2000 model to attempt to recover silicate and carbonaceous masses of interstellar dust from simulated, noise-free observations of extinction.

First, the extinction cross section per H nucleon for ISM dust is calculated for plausible size distributions (either MRN or WD01):

$$\begin{eqnarray}&&{\sigma }_{{\rm{ext}}}(\lambda )=\displaystyle \sum _{i}{\displaystyle \int }_{0}^{\infty }{C}_{{\rm{ext}},i}(\lambda ,a){\left(\displaystyle \frac{1}{{n}_{{\rm{H}}}}\displaystyle \frac{{{dn}}_{i}}{{da}}\right)}_{{mod}}{da},\end{eqnarray} \tag{ 2 }$$

where ${C}_{{\rm{ext}},i}(\lambda ,a)$ is the extinction cross section for a spherical grain composed of material i with radius a, and ${({n}_{{\rm{H}}}^{-1}{{dn}}_{i}/{da})}_{\mathrm{mod}}$ is the size distribution of grain material i for grain model mod.

After calculating ${\sigma }_{{\rm{ext}}}(E)$ for a given grain model near an absorption edge j, we imagine that this has been measured (without noise) and attempt to recover the total volume of material i, ${V}_{i}^{{\rm{true}}}$, by using an attenuation model similar to WAM2000: near absorption edge j we fit the true ${\sigma }_{{\rm{ext}}}$ from Equation (2) with

$$\begin{eqnarray}{\sigma }_{{\rm{ext}}}^{{\rm{fit}}} & = & \displaystyle \sum _{i}\displaystyle \frac{{V}_{{ij}}^{{\rm{fit}}}}{{V}_{i}^{{\rm{true}}}}{\displaystyle \int }_{{a}_{{\rm{min}}}}^{{a}_{{\rm{max}}}}{\left(\displaystyle \frac{1}{{n}_{{\rm{H}}}}\displaystyle \frac{{{dn}}_{i}}{{da}}\right)}_{\mathrm{mod}}\pi {a}^{2}(1-{e}^{-\frac{4}{3}{\alpha }_{i}(\lambda )a}){da}+{C}_{j}.\end{eqnarray} \tag{ 3 }$$

The attenuation coefficients ${\alpha }_{i}$ are presumed to be known, and the shape of the size distribution dn_i/da is assumed to be known, but the multiplier ${V}_{{ij}}^{{\rm{fit}}}$ and the additive offset C_j are free parameters: ${V}_{{ij}}^{{\rm{fit}}}$ is the volume per H of grain material i fit to absorption edge j, and C_j is a constant offset for absorption edge j.³ The domain of the fit contains only wavelengths close to absorption edges, and C_j is fit at each absorption edge.

Fitting is done via the curve_fit function in the scipy Python library (Jones et al. 2001), which implements the Levenberg–Marquardt nonlinear least squares fitting method (Levenberg 1944; Marquardt 1963). For each edge, we fit the extinction over the energy range ${E}_{{\rm{edge}}}\pm {\rm{\Delta }}E$. We tried three values of ${\rm{\Delta }}E$ (10, 20, and 30 eV) to investigate how the results might depend upon the energy range used for the fit.

Two size distributions were considered: the MRN (Mathis et al. 1977) dust grain distribution (${dn}/{da}\propto {a}^{-3.5}$, $5\;{\rm{nm}}\lt a\lt 250\;{\rm{nm}}$) and the WD01 size distribution. Both size distributions use spherical grains. We assume a carbonaceous volume fraction of ${f}_{{\rm{carb}}}=0.488$ for the MRN size distribution and ${f}_{{\rm{carb}}}=0.365$ for the WD01 size distribution. The refractive indices m_i for both materials are taken from Draine (2003).

The fits (using Equation (3)) of the WAM2000-like attenuation model to six X-ray absorption edges for the WD01 size distribution are shown in Figure 4, and for the MRN size distribution in Figure 5. Table 1 shows the ratio of the grain mass estimated by the WAM2000-like fit to the "true" grain mass for various size distribution assumptions and energy ranges.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Material Volumes ${V}_{{ij}}^{{\rm{fit}}}$ Estimated from WAM2000 fit; ${V}_{i}^{{\rm{true}}}$ is the True Volume

ja	Edge	${E}_{{\rm{edge}}}$ b	${\rm{\Delta }}E$	ic	${V}_{{ij}}^{{\rm{fit}}}/{V}_{i}^{{\rm{true}}}$
		(eV)	(eV)		MRNd	WD01 e
1	C K	285.0	10.0	carb.	0.952	0.798
			20.0		1.064	0.847
			30.0		0.949	0.643
2	O K	537.0	10.0	sil.	1.091	1.050
			20.0		0.935	0.911
			30.0		0.843	0.836
3	Fe L	713.5	16.5	sil.	1.084	1.011
			26.5		1.043	0.987
			36.5		0.994	0.954
4	Mg K	1310.0	10.0	sil.	0.985	0.993
			20.0		0.923	0.940
			30.0		0.885	0.909
5	Si K	1845.0	10.0	sil.	0.990	1.021
			20.0		0.962	0.993
			30.0		0.939	0.971
6	Fe K	7123.0	10.0	sil.	1.004	1.066
			20.0		0.995	1.053
			30.0		0.986	1.041

Notes. 

^aj identifies the absorption edge. ^b ${E}_{{\rm{edge}}}=\displaystyle \frac{1}{2}({E}_{{\rm{max}}}+{E}_{{\rm{min}}})$ and ${\rm{\Delta }}E=\displaystyle \frac{1}{2}({E}_{{\rm{max}}}-{E}_{{\rm{min}}})$. ${E}_{{\rm{min}}}$ and ${E}_{{\rm{max}}}$ are the minimum and maximum energies, respectively, over which the fit was performed. ^ci identifies the material: carbonaceous or silicate. ^dMathis et al. (1977). ^eWeingartner & Draine (2001).

Download table as:  ASCIITypeset image

We find that the WAM2000 attenuation model is able to provide moderately accurate estimates of the elemental abundances for Fe L, Mg K, Si K, and Fe K—e.g., errors of $\lesssim 5\%$ for the abundance of Fe based on the Fe L₁, ${{\rm{L}}}_{\mathrm{2,3}}$ edges. However, because the WAM2000 model neglects scattering, the wavelength dependence of the extinction is not well-modeled. Consequently, attempts to identify the chemical state (e.g., Fe metal versus Fe₃O₄) from the details of the edge profile are prone to error, as is evident in Figures 4 and 5 from the poor fits of the WAM2000 model to the "true" extinction profiles calculated for the same material.

As demonstrated in Section 2.2, the WAM2000 attenuation model—which assumes pure absorption—has systematic errors $\gtrsim 10\%$ for energies below the Si K edge. A more accurate attenuation model is necessary to obtain reliable abundances from X-ray attenuation measurements. Accurate modeling of attenuation requires an accurate complex dielectric function and an algorithm for computing absorption and scattering by the dust grains. ADT is a natural choice for the latter: it can handle non-spherical grain geometries, is accurate for large grains, and is more computationally efficient than Mie theory or Rayleigh–Gans.

3.1. Intuition Behind ADT

ADT was invented by van de Hulst (1957) to treat the problem of scattering and absorption by a particle that is optically soft ($| m-1| \ll 1$) but large compared to the wavelength of incident light ($x\equiv {ka}\gg 1$).

Because $x\gg 1$, the concept of independent rays of light passing through the grain is a valid approximation. And, because $| m-1| \ll 1$, refraction and reflection effects are small and may be ignored. The grain can thus be thought of as providing local phase shifts to the incident plane wave. Absorption of the incident plane wave also occurs if ${\rm{Im}}(m)\ne 0$.

Consider a grain of arbitrary geometry. Define a plane S just beyond the extent of the grain and normal to $\hat{z}$, the direction of propagation. The plane, S, is located at z = 0, and the grain is confined to $z\lt 0$. Define $\hat{x}$ and $\hat{y}$ to be orthogonal unit vectors that lie in S with both $\hat{x}$ and $\hat{y}$ orthogonal to $\hat{z}$. Define the "shadow function" on S:

$$\begin{eqnarray}&&f(x,y)=1-{e}^{i{\rm{\Phi }}(x,y)},\end{eqnarray} \tag{ 4 }$$

where ${\rm{\Phi }}(x,y)$ is the complex phase shift:

$$\begin{eqnarray}&&{\rm{\Phi }}(x,y)\equiv k{\displaystyle \int }_{-\infty }^{0}[m(x,y,z)-1]{dz},\end{eqnarray} \tag{ 5 }$$

m is the (complex) index of refraction, and $k\equiv 2\pi /\lambda$.

To obtain far-field solutions to the scattered field, Huygen's principle is applied to the shadow function f on surface S. We refer the reader to Draine & Allaf-Akbari (2006) for a brief derivation of the extinction, absorption, scattering, and differential scattering cross sections in the context of ADT, and to van de Hulst (1957) for a more rigorous and detailed treatment.

The ADT formulae for absorption, scattering, and extinction cross sections are:

$$\begin{eqnarray}&&{C}_{{\rm{abs}}}={\displaystyle \int }_{S}(1-{e}^{-2{{\rm{\Phi }}}_{2}}){dS}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{C}_{{\rm{sca}}}={\displaystyle \int }_{S}| f{| }^{2}{dS}={\displaystyle \int }_{S}(1-2\mathrm{cos}{{\rm{\Phi }}}_{1}{e}^{-{{\rm{\Phi }}}_{2}}+{e}^{-2{{\rm{\Phi }}}_{2}}){dS}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{C}_{{\rm{ext}}}=2{\displaystyle \int }_{S}(1-\mathrm{cos}{{\rm{\Phi }}}_{1}{e}^{-{{\rm{\Phi }}}_{2}}){dS},\end{eqnarray} \tag{ 8 }$$

where ${{\rm{\Phi }}}_{1}\equiv {\rm{Re}}({\rm{\Phi }})$ and ${{\rm{\Phi }}}_{2}\equiv {\rm{Im}}({\rm{\Phi }})$. The differential scattering cross section is given by

$$\begin{eqnarray}&&\displaystyle \frac{{{dC}}_{{\rm{sca}}}}{d{\rm{\Omega }}}=\displaystyle \frac{| S(\hat{{\boldsymbol{n}}}){| }^{2}}{{k}^{2}},\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&S(\hat{{\boldsymbol{n}}})=\displaystyle \frac{{k}^{2}}{2\pi }{\displaystyle \int }_{S}{e}^{{ik}(\hat{{\boldsymbol{n}}}\cdot {\boldsymbol{x}})}f({\boldsymbol{x}})d{\boldsymbol{S}}.\end{eqnarray} \tag{ 10 }$$

3.1.1. ADT for Spheres

For spherical grains, ADT yields analytic expressions for ${C}_{{\rm{ext}}}$, ${C}_{{\rm{abs}}}$, and ${C}_{{\rm{sca}}}$:⁴

$$\begin{eqnarray}{Q}_{{\rm{ext}}} & \equiv & \displaystyle \frac{{C}_{{\rm{ext}}}}{\pi {a}^{2}}=2+\displaystyle \frac{4}{| \rho {| }^{2}}\{\mathrm{cos}2\beta -{e}^{-{\rho }_{2}}[\mathrm{cos}({\rho }_{1}-2\beta )\\ & & +| \rho | \mathrm{sin}({\rho }_{1}-\beta )]\}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{Q}_{{\rm{abs}}}\equiv \displaystyle \frac{{C}_{{\rm{abs}}}}{\pi {a}^{2}}=1+\displaystyle \frac{{e}^{-2{\rho }_{2}}}{{\rho }_{2}^{2}}+\displaystyle \frac{{e}^{-2{\rho }_{2}}-1}{2{\rho }_{2}^{2}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{Q}_{{\rm{sca}}}\equiv \displaystyle \frac{{C}_{{\rm{sca}}}}{\pi {a}^{2}}={Q}_{{\rm{ext}}}-{Q}_{{\rm{abs}}},\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&\rho \equiv 2{ka}(m-1),\end{eqnarray} \tag{ 14 }$$

${\rho }_{1}\equiv {\rm{Re}}(\rho )$, ${\rho }_{2}\equiv {\rm{Im}}(\rho )$, and $\beta \equiv \mathrm{arccos}({\rho }_{1}/| \rho | )$. $S(\theta )$ becomes⁵ ${}^{,}$ ⁶

$$\begin{eqnarray}{S}_{{\rm{ADT}}}(\theta ) & = & {({ka})}^{2}{\displaystyle \int }_{0}^{\pi /2}{{du}{J}}_{0}({ka}\mathrm{sin}\theta \mathrm{cos}u)\\ & & \times (1-{e}^{i\rho \mathrm{sin}u})\mathrm{sin}u\mathrm{cos}u,\end{eqnarray} \tag{ 15 }$$

where J₀ is the Bessel function of the first kind (of order zero); for $\theta =0$,

$$\begin{eqnarray}&&{S}_{{\rm{ADT}}}(\theta =0)=\displaystyle \frac{1}{2}+\displaystyle \frac{1+{{ie}}^{i\rho }(\rho +i)}{{\rho }^{2}}.\end{eqnarray} \tag{ 16 }$$

3.1.2. ADT Accuracy

ADT provides a natural complement to the Rayleigh–Gans approximation in the X-ray regime; small grains are accurately modeled by the Rayleigh–Gans approximation (or Mie theory), while larger grains are accurately modeled by ADT. There is an analytic solution to the ADT equations for spherical scatterers, so the computational cost of using ADT over Rayleigh–Gans is negligible. It is also reasonably straightforward to extend ADT to other geometries.

Even for small grains where the "ray optics" approach of ADT fails, ADT will still produce accurate estimates of extinction, since, in this case, the extinction is dominated by absorption. To illustrate this, Figure 6 (top left) compares ${Q}_{{\rm{abs}}}$, ${Q}_{{\rm{sca}}}$, and ${Q}_{{\rm{ext}}}$ for silicate spheres calculated with Mie theory, ADT, the Rayleigh–Gans approximation, and the WAM2000 absorption-only estimate.

Zoom In Zoom Out Reset image size

3.2. GGADT: General Geometry Anomalous Diffraction Theory

The authors have written a Fortran 95 program GGADT that uses ADT to calculate (1) the energy-dependent scattering and absorption cross sections, and (2) the differential scattering cross section for grains of arbitrary geometry and composition. GGADT is fast, portable, GNU-compliant, and well documented.⁷ GGADT uses the General Prime Factor Algorithm (GPFA) of Temperton (1992) to do fast Fourier transforms. A brief description of GGADT usage can be found in Appendices A and B.

3.2.1. An Application of GGADT: Effect of Grain Geometry on X-Ray Extinction

The geometry of dust grains is not currently well constrained. Polarization of starlight implies that dust grains are not spherically symmetric. The next simplest grain geometry is the spheroid; spheroidal grain models are able to reproduce observations of starlight polarization and extinction (Kim & Martin 1995; Draine & Fraisse 2009).

However, for plausible dust evolution scenarios, dust grains are likely more complicated than single-material spheroids or even ellipsoids. ISM grains could have irregular geometries as well as inhomogeneous composition. Some authors (e.g., Mathis & Whiffen 1989; Henning & Stognienko 1993; Stognienko et al. 1995) have argued for highly porous geometries.

To illustrate the possible effects that grain geometry might have on abundance measurements based on X-ray extinction, we employ GGADT to compute the extinction cross sections for five example grain geometries. The size is specified by the radius of an equal-volume sphere,

$$\begin{eqnarray}&&{a}_{{\rm{eff}}}\equiv {(3V/4\pi )}^{1/3},\end{eqnarray} \tag{ 17 }$$

where V is the volume of the solid material. The five grain geometries used are (1) a sphere, (2) an oblate spheroid, (3) a prolate spheroid, (4) a Ballistic Agglomeration (BA) aggregate, and (5) a BAM2 aggregate each with the same mass (${a}_{{\rm{eff}}}=0.2\;\mu {\rm{m}}$) and silicate composition.

BA aggregates are constructed by single-size spherical monomers arriving on random trajectories and adhering to their initial point of contact; BAM2 aggregates require arriving monomers (after the third) to make contact with a total of three other monomers prior to the arrival of the next monomer. BAM2 aggregates have porosities P significantly less than that of BA aggregates. A detailed description of BA and BAM2 aggregation is given in Shen et al. (2008). Figure 7 shows examples of BA and BAM2 agglomerates.

Zoom In Zoom Out Reset image size

This paper employs the definition of porosity from Shen et al. (2008). For a given grain, define an equivalent ellipsoid (EE) as the uniform density ellipsoid that has the same mass and moment of inertia tensor as the grain. The porosity of the grain P is then defined by ${a}_{{\rm{EE}}}^{3}(1-P)\equiv {a}_{{\rm{eff}}}^{3}$, where $(4\pi /3){a}_{{\rm{EE}}}^{3}$ is the volume of the EE. Thus, ${a}_{{\rm{EE}}}={a}_{{\rm{eff}}}{(1-P)}^{-1/3}$.

For grains with complex geometries, such as BA and BAM2 aggregates, the shadow function, $f(x,y)$ (see Equation (4)) is evaluated on an N_x × N_y grid, to facilitate the computation of discrete Fourier transforms (DFTs). The choice of N_x and N_y determines how accurately $f(x,y)$ will be described on small linear scales. Figure 8 shows how the GGADT result for a sphere depends on the chosen grid size. We see that 64 × 64 yields results that are accurate to a few percent (heights of scattering peaks, and positions of maxima and minima). The reason that high-resolution numerical representations (e.g., 2048 × 2048) of the shadow function are not required to achieve  1% accuracy arises from the nature of Equation (10); for small-angle scattering ($\mathrm{sin}\theta \ll 1$), only the long-wavelength contributions to the Fourier transform of f are relevant to the calculation of $d{\sigma }_{{\rm{sca}}}/d{\rm{\Omega }}$. We will generally use grid resolutions 128 × 128 or higher for calculations in this paper unless otherwise specified.

Zoom In Zoom Out Reset image size

Extinction cross sections calculated for several ${a}_{{\rm{eff}}}=0.2\;\mu {\rm{m}}$ grains are shown in Figures 9 and 10. The dielectric function of the grains was taken to be that of MgFeSiO₄ olivine described in Draine (2003). All cross sections are averaged over 64 random grain orientations. All calculations were done using ADT. For the spherical grain, results calculated with GGADT are indistinguishable from Mie theory. For the spheroidal and agglomerate grains, the extinction was calculated using the GGADT code, with the shadow function $f(x,y)$ (see Equation (4)) sampled on a 512 × 512 grid, providing excellent accuracy (see Figure 8 below). The aggregates used in Figures 9 and 10 are those in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3. Porosity

As stated earlier in this section, higher porosity increases absorption efficiency. This trend, along with the effect of porosity on scattering, is shown in Figure 11. For ${a}_{{\rm{eff}}}=0.2\mu {\rm{m}}$ grains, the scattering efficiency decreases as porosity increases. The extinction efficiency decreases with increased porosity except near the Fe L and O K absorption edges, where the increase in absorption efficiency dominates over the decreasing scattering efficiency.

Zoom In Zoom Out Reset image size

Orientation-averaged X-ray extinction cross sections do not differ appreciably between the spherical and spheroidal dust grains in Figure 9, but larger differences are found for random aggregates. However, the absorption and scattering cross sections of random aggregates, as shown in Figure 11, do deviate significantly from those of spherical grains. The absorption efficiency of random aggregates is increased at all photon energies relative to an equal-mass sphere; volume elements of porous grains are exposed to a larger fraction of the incident light than compact grains, and therefore the grain as a whole absorbs light more efficiently. In compact grains, parts of the grain are "shadowed" (i.e., exposed to less of the incident light) by grain material closer to the light source, and thus these shadowed regions do not absorb as much light as in porous grains.

The effects of porosity on the scattering cross section are more complicated, and we pause to provide a semi-analytical discussion of two limiting cases (the optically thin and optically thick regimes) below.

In the regime where $| {\rm{\Phi }}| \ll 1$ (the optically thin regime), Equation (7) for spherical grains becomes

$$\begin{eqnarray}&&{C}_{{\rm{sca}}}=\int | 1-{e}^{i{\rm{\Phi }}}{| }^{2}{dxdy}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\approx \int | 1-(1+i{\rm{\Phi }}){| }^{2}{dxdy}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&=\ \int | {\rm{\Phi }}{| }^{2}{dxdy}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&=\ {k}^{2}| {m}_{{\rm{eff}}}-1{| }^{2}{\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{0}^{{a}_{{\rm{EE}}}}Z{(r)}^{2}{rdrd}\theta ,\end{eqnarray} \tag{ 21 }$$

where $Z(r)\equiv 2\sqrt{{a}_{{\rm{EE}}}^{2}-{r}^{2}}$ is the length of the grain along the z axis at $r=\sqrt{{x}^{2}+{y}^{2}}$, and ${m}_{{\rm{eff}}}$ is the effective index of refraction which depends upon the porosity of the grain,

$$\begin{eqnarray}&&({m}_{{\rm{eff}}}-1)\approx ({m}_{0}-1)(1-P),\end{eqnarray} \tag{ 22 }$$

where m₀ is the index of refraction for the grain material.

Thus, in the optically thin regime,

$$\begin{eqnarray}&&{C}_{{\rm{sca}}}\approx 2\pi {k}^{2}({a}_{{\rm{eff}}}^{4}{(1-P)}^{-4/3})(| {m}_{0}-1{| }^{2}{(1-P)}^{2})\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&=\ 2\pi {k}^{2}{a}_{{\rm{eff}}}^{4}| {m}_{0}-1{| }^{2}{(1-P)}^{2/3},\end{eqnarray} \tag{ 24 }$$

which implies that an increase in porosity produces a decrease in the scattering cross section. This accounts for the decrease in scattering efficiency at the Fe L edge shown in Figure 11. However, lower energy absorption edges (e.g., the O K edge), exhibit smaller decrements in the scattering cross section close to the absorption edge. This is because, as will be discussed below, porosity has the opposite effect on scattering when not in the optically thin regime.

In the opaque limit $[({ka}){\rm{Im}}(m-1)\gg 1]$, however, we have that ${C}_{{\rm{sca}}}\approx {C}_{{\rm{abs}}}\approx {A}_{{\rm{proj}}}$ so ${C}_{{\rm{ext}}}\approx 2{A}_{{\rm{proj}}}$ (see Bohren & Huffman (1983) for a discussion of the "extinction paradox"), where ${A}_{{\rm{proj}}}$ is the projected area of the grain along the direction of propagation.

For fixed ${a}_{{\rm{eff}}}$, increased porosity P leads to a larger ${A}_{{\rm{proj}}}$. Thus, in the opaque limit, increasing the porosity of a given grain will cause an increase in ${C}_{{\rm{sca}}}$.

As can be seen in Figures 9 and 10 the absorption edge fine structure is significantly affected by grain geometry (especially porosity). Sensitivity to dust geometry means that absorption-edge fine structure is indeed a laboratory for constraining dust geometry as well as composition, but also calls into question the significance of claimed abundance constraints derived from fits to X-ray absorption edge fine structure.

3.4. Effective Medium Theory (EMT) Calculations

EMT have been used since Maxwell Garnett (1904) and Bruggeman (1935) to obtain an effective dielectric function, ${\epsilon }_{{\rm{eff}}}$ for an inhomogeneous particle composed of several materials with dielectric functions ${\epsilon }_{1}$, ${\epsilon }_{2}$, ... . If the particle itself is approximately spherical, Mie theory can then be used to calculate scattering and absorption by the particle.

More sophisticated EMTs (see, e.g., Stognienko et al. 1995) have been constructed to try to take into account the shapes of inclusions and voids in inhomogeneous grains with complex structures. Valencic & Smith (2015) investigated how well observations of X-ray halos could be modeled by a variety of dust models, including those of Zubko et al. (2004) (hereafter ZDA2004), and concluded that the families of porous dust models described in ZDA2004 did not fit observations as well as models without porosity.

At X-ray energies, the wavelength may be comparable or even smaller than the sizes of inclusions and voids, and the validity of EMT is uncertain. Furthermore, EMT aims to reproduce the mean polarization, while scattering depends on the mean square polarization, and hence it is of interest to compare ADT calculations against those of EMT to determine how well EMT approximates the scattering and absorption properties of BA and BAM2 grains.

Figure 12 compares ADT calculations of scattering and absorption by ballistic agglomerates against Mie theory calculations for a sphere with radius $a={a}_{{\rm{EE}}}$ and with a uniform index of refraction $m={m}_{{\rm{eff}}}$ obtained from EMT. At X-ray energies, the different EMTs (e.g., Maxwell Garnett 1904 and Bruggeman 1935) give essentially identical ${m}_{{\rm{eff}}}$. EMT agrees best with the most compact ballistic aggregates (BAM2), but systematically underestimates the scattering cross section and overestimates the absorption cross sections of the ballistic agglomerate grains. These effects partially cancel when the extinction is computed, but the height of the absorption edge—a key diagnostic for ISM abundances—remains overestimated. The accuracy of the EMT approach decreases with increasing porosity.

Zoom In Zoom Out Reset image size

GGADT also calculates the differential scattering cross sections ${{dC}}_{{\rm{sca}}}/d{\rm{\Omega }}\equiv \pi {a}_{{\rm{eff}}}^{2}{{dQ}}_{{\rm{sca}}}/d{\rm{\Omega }}$. Figure 13 shows full two-dimensional (2D) results for ${{dQ}}_{{\rm{sca}}}/d{\rm{\Omega }}$ for BA and BAM2 aggregates with ${a}_{{\rm{eff}}}=0.2\;\mu {\rm{m}}$, each composed of N = 256 astrosilicate monomers (monomer radii $0.0315\;\mu {\rm{m}}$). For each grain, the shadow function for one particular orientation is shown, together with the 2D scattering pattern for that orientation. The BA and BAM2 aggregates have equal masses, but the BAM2 structure is more compact, and its central scattering peak, with ${\rm{\Delta }}\theta \approx \lambda /D$ (where D is a characteristic "diameter"), is noticeably broader than for the BA structure. For the selected orientations, the BA and BAM2 aggregates are each elongated along the $\hat{{\boldsymbol{x}}}$ axis, resulting in scattering patterns that are elongated in the $\hat{{\boldsymbol{y}}}$ direction.

Zoom In Zoom Out Reset image size

The right-hand column of Figure 13 shows the 2D scattering pattern averaged over 300 random (in three-dimensions) orientations of the grain. The BAM2 grain continues to have a significantly broader central peak than the BA grain.

Figure 14 shows the azimuthally averaged ${{dQ}}_{{\rm{sca}}}/d{\rm{\Omega }}$ from the right column of Figure 13, together with ${dQ}/d{\rm{\Omega }}$ for an equal-mass sphere. As expected, the sphere has most of the scattered power at $\theta \lt \lambda /D=320^{\prime\prime}$, but the BAM2 and BA structures have narrower central peaks, because they are more extended (larger D).

Zoom In Zoom Out Reset image size

The similarity of the BA and BAM2 results at $\theta \gt 3000^{\prime\prime}$, with regularly spaced maxima and minima, is striking. At these large scattering angles, ${{dQ}}_{{\rm{sca}}}/d{\rm{\Omega }}$ is determined by the high-spatial frequency portion of the shadow function $f(x,y)$ (see Equation (10)), and this comes from the spherical monomers. Since the BA and BAM2 structures considered here have the same sized monomers, they have very similar ${{dQ}}_{{\rm{sca}}}/d{\rm{\Omega }}$ at large scattering angles. In the ISM, we of course do not expect agglomerate structures to be composed of single-sized monomers, and scattering halos will not show such maxima and minima.

Observations of X-ray absorption and scattering by the ISM provide a valuable opportunity for testing dust models and measuring interstellar abundances of heavy elements, both in the gas and solid phases. By investigating the near-edge fine structure of absorption (extinction) edges, one can deduce information about the (1) geometry (e.g., porosity) and (2) composition of interstellar grains. The dependence of the X-ray scattering halo on both scattering angle and photon energy also constrains the size, structure, and composition of the dust population.

Efforts to infer the dust composition, size distribution, and grain geometry (shape and porosity) rely on detailed comparisons of observations (e.g., energy-dependence of the extinction, especially near absorption edges) with theoretical calculations for various grain models. In this paper we have tested the widely used approach of WAM2000 (Wilms et al. 2000) for calculating dust attenuation.

The WAM2000 method, which ignores scattering, does allow the elemental abundances in dust to be estimated with reasonable accuracy (errors ≲10%–15%). However, because scattering is neglected, the WAM2000 model does not provide accurate profiles of the near-edge fine structure in the extinction, and therefore should not be used when attempting to distinguish between different possible chemical forms in which the element of interest may be present (e.g., using the Fe L edge to try to distinguish metallic Fe, various Fe oxides, or Fe contained in silicates). Inaccurate modeling of the energy dependence of extinction near absorption edges can easily lead to incorrect conclusions regarding the compositions of interstellar grains.

Fortunately, for X-ray energies, the energy-dependent extinction, as well as scattering halos, can readily be calculated accurately using ADT. For spheres, the computational cost of accurate ADT calculations is trivial, effectively equal to that of the WAM2000 approximations. One important advantage of the ADT approach is that it is readily extended to other grain shapes, such as spheroids, clusters of spheres, or indeed whatever shape is of interest to the researcher.

The authors have developed and made available an open source ADT Fortran suite—GGADT—that can perform the necessary ADT calculations for both integrated absorption, scattering, and extinction cross sections, as well as for differential scattering cross sections. GGADT utilizes several numerical techniques to make the ADT calculations quite fast without loss of accuracy (see Appendices C and D).

Above we used GGADT to calculate absorption and scattering by spheres, spheroids, and random aggregates of sphere, to show explicitly that the near-edge absorption fine structure is sensitive to grain geometry, in particular porosity. The fact that grain geometry alone can significantly alter absorption edge profiles adds another complication to future studies of near-edge X-ray extinction fine structure, but in principle also provides another way to obtain information about the shapes and structures of interstellar grains.

Figures 13 and 14 show how X-ray scattering halos are sensitive to the geometric structure of the grains. As shown by Draine & Allaf-Akbari (2006), aligned interstellar grains are expected to produce noncircular X-ray scattering halos, with large enough asymmetries to be measurable by Chandra. Grains with significant porosity produce significantly narrower forward scattering peaks than equal-mass nonporous grains. Heng & Draine (2009) showed that the implied changes to the X-ray scattering halo could be used to test whether interstellar grains are primarily dense, compact structures versus porous "fluffy" aggregates. This will be the subject of future investigations, to see what grain geometries are compatible with observations of extinction and polarization at optical wavelengths, and scattering at X-ray energies.

The principal conclusions of this paper are as follows:

• 1.  
  Neither the WAM2000 approximation, nor the Rayleigh–Gans approximation, should be used to calculate near-edge fine structure, even for spherical grains. Both of these methods introduce systematic errors, which can be large for grain sizes and X-ray energies of astrophysical interest (see Figure 6).
• 2.  
  At X-ray energies, ADT is highly accurate and can be used for spheres (where agreement with exact Mie theory is excellent) as well as for more realistic grain geometries.
• 3.  
  For spheres, ADT extinction, absorption, and scattering cross sections are given by simple analytic formulae (11)–(13). For more general geometries, ADT calculations of absorption, scattering, and extinction cross sections require only evaluation of the integrals in Equations (6)–(8), which is computationally straightforward.
• 4.  
  "Naive" calculation of X-ray scattering halos for general geometries can be computationally demanding (both operations and memory) when high angular resolution is desired. We implement a method that greatly increases the computational efficiency, allowing accurate and high-resolution calculation of scattering halos to be performed with modest requirements of memory and CPU time.
• 5.  
  We make available GGADT, an efficient open-source code to calculate X-ray absorption and scattering by grains with general geometries.
• 6.  
  We use GGADT to calculate X-ray scattering and absorption by spheres, spheroids, and random aggregates. The near-edge extinction versus energy can differ significantly among grains with the same mass and composition, but different shape. Grain geometry must therefore be taken into consideration when seeking to deduce grain composition from observations of near-edge X-ray fine structure.
• 7.  
  In the X-ray regime, for complex grain geometries (e.g., porous grains), estimates of absorption and scattering cross sections made using homogeneous spheres with an "effective" refractive index obtained from an EMT do not, in general, provide accurate results. The scattering tends to be underestimated, and the near-edge fine structure is not accurately reproduced (see Figure 9).
• 8.  
  We use GGADT to calculate X-ray scattering halos for porous grains, extending previous work by Heng & Draine (2009). For fixed grain mass, increasing porosity leads to narrowing of the forward scattering peak and characteristic halo angle (see Figure 14). Observations of X-ray scattering halos can, therefore, provide constraints on the structure of interstellar grains.

We thank Brandon Hensley for stimulating discussion, and the anonymous referee for helpful comments. This work was supported in part by NSF grants AST1008570 and AST1408723.

GGADT is a Fortran 90 code, provided with a GNU autotools configure script⁸ , which should work on most Unix/Linux/BSD-based operating systems (this includes Mac OSX). The full source code and documentation (both pdf and html) are available at http://www.ggadt.org.

There are two fundamental calculations that GGADT can perform:

• 1.  
  Total cross sections of a given grain and grain composition over a range of energies.
• 2.  
  Differential scattering cross section of a given grain and grain composition.

Both calculations can be averaged over a number of orientations; users can control how the orientation averaging is performed; either by choosing random orientations, or by evenly dividing the orientations over euler angles, or by specifying a list of orientations over which to average the calculations.

Input to GGADT can be done either on the command line, in a parameter file, or both (though combining command line arguments and parameter files is not recommended).

To calculate the differential scattering cross section for a cluster of silicate spheres with effective radius $a=0.1\;\mu {\rm{m}}$, at an energy of 500 eV

• ggadt --grain--geometry = 'spheres' --aeff = 0.1 --ephot = 0.5 --norientations = 100 --agglom-file = [..]/BA.256.1.targ --material--file = [..]/index_silD03

Here, [..] should be replaced with filepaths to the locations of these files. Files that describe the geometry of BA/BAM1/BAM2 grains are available online at http://ggadt.org/additional_files.html. The --material--file parameter expects an "index" file; several of these are also provided online at the previous url.

To calculate the total cross section of an ellipsoidal silicate grain with axis ratios x:y:z = 1:2:3 at energies close to the O K edge (520–560 eV), one can run

• ggadt --integrated --grain--geometry='ellipsoid' --aeff=0.1 --grain--axis-x=1 --grain--axis--y=2 --grain-axis-z= 3 --ephot-min=0.52 --ephot--max=0.56 --material--file=[..]/index_silD03

There are of order 30 different parameters and flags that one can set in GGADT; descriptions of all of these are given in the documentation.⁹ An example directory contains two sample cases for running GGADT along with a python script for plotting the results.

The shadow function $f(x,y)$ is represented on a numerical grid. To obtain ${{dC}}_{{\rm{sca}}}/d{\rm{\Omega }}$, we require $S(\hat{n})$, proportional to the Fourier transform of f (see Equation (10)). Naively, one could calculate the orientation-averaged ${{dC}}_{{\rm{sca}}}/d{\rm{\Omega }}$ by obtaining the complex shadow function $f({\boldsymbol{x}})$ on a numerical grid of N_x × N_y points, then using a two-dimensional FFT to calculate $S(\theta ,\phi )$, converting this to ${{dC}}_{{\rm{sca}}}/d{\rm{\Omega }}$, repeating for many orientations, and averaging the results.

For orientation-averaging, however, we can sample a single azimuthal angle, taking $\phi =0$ without loss of generality. In this case, Equation (10) simplifies to

$$\begin{eqnarray}&&S(\theta ,\phi =0)=\displaystyle \frac{{k}^{2}}{2\pi }{\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{e}^{{ikx}\mathrm{sin}\theta }f(x,y){dxdy}\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{{k}^{2}}{2\pi }{\displaystyle \int }_{-\infty }^{\infty }{e}^{{ikx}\mathrm{sin}\theta }\underset{g(x)}{\underbrace{\left[{\displaystyle \int }_{-\infty }^{\infty }f(x,y){dy}\right]}}{dx}\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{{k}^{2}}{2\pi }{\displaystyle \int }_{-\infty }^{\infty }{e}^{{ikx}\mathrm{sin}\theta }g(x){dx}\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{{k}^{2}}{2\pi }\bar{g}(\kappa ),\end{eqnarray} \tag{ 28 }$$

where $\bar{g}$ is the Fourier transform of $g(x)\equiv {\displaystyle \int }_{-\infty }^{\infty }f(x,y){dy}$, and $\kappa =-\displaystyle \frac{1}{2\pi }k\mathrm{sin}\theta$.

By reducing a two-dimensional calculation to a one-dimensional one, we speed up the calculation for a single orientation. The number of operations required to obtain $S(\theta ,\phi =0)$ is $\sim { \mathcal O }({N}_{y})+{ \mathcal O }({N}_{x}\mathrm{log}\;{N}_{x})$, compared to $\sim { \mathcal O }({N}_{x}{N}_{y}\mathrm{log}\;{N}_{x}{N}_{y})$ for the two-dimensional FFT. However, since we are approximating the orientation-averaged differential scattering cross section by averaging over a finite number of orientations, the 1D FFT trick will require more orientations to reach the same level of accuracy as the 2D FFT, since we could sample many ϕ values with a single 2D FFT. However, as shown in Figure 15, only a modest number of orientations (${ \mathcal O }({10}^{2})$) is needed to obtain accurate numerical results. This means that the speedup gained from using a 1D FFT can outweigh the cost of averaging over extra orientations ($\sim {10}^{2}{N}_{x}\mathrm{log}\;{N}_{x}\lt \sim {N}_{x}{N}_{y}\mathrm{log}\;{N}_{x}{N}_{y}$).

Zoom In Zoom Out Reset image size

An additional contribution to the computation time comes from the calculation of the shadow function. For spheres and ellipsoids, the shadow function computation time scales as ${ \mathcal O }({N}_{x}{N}_{y})$. For agglomerations of spheres, this process becomes approximately ${ \mathcal O }({N}_{m}^{1/3}{N}_{x}{N}_{y})$, where N_m is the number of monomers.

The reason that aggregates require a ${ \mathcal O }({N}_{m}^{1/3}{N}_{x}{N}_{y})$ computation time is the following: the contribution to the shadow function from each monomer is computed over a sub-grid of the full shadow function grid, with ${N}_{x}^{\prime }{N}_{y}^{\prime }\approx {({a}_{m}/a)}^{2}{N}_{x}{N}_{y}$, where a_m is the radius of the monomer and a is the effective radius of the grain. Since ${N}_{m}{a}_{m}^{3}={a}^{3}$, we have that ${a}_{m}={{aN}}_{m}^{-1/3}$ and thus ${N}_{x}^{\prime }{N}_{y}^{\prime }={N}_{m}^{-2/3}{N}_{x}{N}_{y}$. The computation time scales as ${N}_{m}{N}_{x}^{\prime }{N}_{y}^{\prime }={N}_{m}^{1-2/3}{N}_{x}{N}_{y}={N}_{m}^{1/3}{N}_{x}{N}_{y}$.

For custom grain geometries, for which no assumptions about the structure of the grain can be made, computation of the shadow function will scale as ${N}_{x}{N}_{y}{N}_{z}$, where N_z is the number of samplings along the $\hat{z}$ direction. In this case, computation of the shadow function dominates the computation time, and thus 2D FFT's are the more efficient solution. In the former cases, (for clusters of spheres, as long as ${N}_{m}\ll {N}_{x}{N}_{y}{N}_{z}$), the FFT dominates the computation time, and thus 1D FFT's are more time efficient.

Recall that the (one-dimensional) discrete Fourier transform (DFT) of an array of values f_n is given by

$$\begin{eqnarray}&&{\hat{f}}_{m}=\displaystyle \sum _{n=0}^{N}{e}^{-2\pi {im}\frac{n}{N}}{f}_{n}\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&=\ \displaystyle \sum _{n=0}^{N}{e}^{-2\pi {im}\frac{{x}_{n}-{x}_{0}}{L}}{f}_{n}\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&=\ {e}^{2\pi {{imx}}_{0}/L}\displaystyle \sum _{n=0}^{N}{e}^{-2\pi i\left(\frac{m}{L}\right){x}_{n}}{f}_{n}\end{eqnarray} \tag{ 31 }$$

Assume (for the sake of simplicity) that ${N}_{x}={N}_{y}={N}_{z}=N$. Here, ${\rm{\Delta }}x$ is the width of a grid element and N is the total number of grid elements along one direction. When we calculate the DFT of g(x), we obtain $S({\theta }_{m},\phi =0)$ for a set of angles ${\theta }_{m}$, where

$$\begin{eqnarray}&&\mathrm{sin}{\theta }_{m}=\displaystyle \frac{2\pi m}{{Nk}{\rm{\Delta }}x}.\end{eqnarray} \tag{ 32 }$$

In astronomical contexts, usually only small angles (${\theta }_{{\rm{max}}}\lesssim {10}^{4}$ arcsec) are relevant, and thus $\mathrm{sin}{\theta }_{m}\approx {\theta }_{m}$. The smallest (non-zero) scattering angle for which we can compute ${{dC}}_{{\rm{sca}}}/d{\rm{\Omega }}$ is

$$\begin{eqnarray}&&\theta \approx \displaystyle \frac{2\pi }{{kL}}=\displaystyle \frac{\lambda }{L},\end{eqnarray} \tag{ 33 }$$

where $L=N{\rm{\Delta }}x$ is the extent of the grid. If L = diameter of the target, then this angle is approximately the size of the central scattering lobe. Generally, we wish to resolve the central lobe with resolution

$$\begin{eqnarray}&&\delta {\theta }_{{\rm{sca}}}=\displaystyle \frac{2\pi /{kL}}{{N}_{{\rm{sca}}}},\end{eqnarray} \tag{ 34 }$$

where ${N}_{{\rm{sca}}}$ is the number of angles sampling the central scattering peak.

For a fixed grid resolution ${\rm{\Delta }}x$ and photon wavenumber k, we can increase the resolution of the differential scattering cross section at small angles by extending the xy grid on which the shadow function $f(x,y)$ is defined, and setting $f(x,y)=0$ beyond the actual shadow. If ${L}_{{\rm{ext}}}$ is the size of this extended grid, then the DFT will have a resolution

$$\begin{eqnarray}&&{\rm{\Delta }}\theta \approx \displaystyle \frac{2\pi }{{{kL}}_{{\rm{ext}}}}.\end{eqnarray} \tag{ 35 }$$

This technique is also known as "padding." However, padding is a memory-intensive (and time-intensive) process. We seek a more efficient algorithm that provides high resolution at small scattering angles and avoids calculating $S(\theta )$ for $\theta \gt {\theta }_{{\rm{max}}}$.

One can throw away calculations of $\theta \gt {\theta }_{{\rm{max}}}$ by performing what is known as a "pruned FFT."¹⁰ Using a pruned FFT changes the asymptotic execution time of a one-dimensional FFT from ${ \mathcal O }(N\mathrm{log}\;N)$ to ${ \mathcal O }(N\mathrm{log}\;K)$, where K is the largest index of the DFT that we care about; in this context, K corresponds to the m value of ${\theta }_{{\rm{max}}}$:

$$\begin{eqnarray}&&K\approx \displaystyle \frac{{Nk}{\rm{\Delta }}x\mathrm{sin}{\theta }_{{\rm{max}}}}{2\pi }.\end{eqnarray} \tag{ 36 }$$

Utilizing a pruned FFT does not speed things up very much on its own unless $K\ll N$, and we are still stuck with needing to pad the grid to increase the angular resolution of ${{dC}}_{{\rm{sca}}}/d{\rm{\Omega }}$. There is, however, another way of improving the angular resolution of ${{dC}}_{{\rm{sca}}}/d{\rm{\Omega }}$ that does not involve extending (padding) the shadow function grid.

Instead of performing a single DFT on an extended (padded) grid, we can perform several offset DFT's on the unpadded grid. A single DFT on the unpadded grid gives a coarse sampling of ${{dC}}_{{\rm{sca}}}/d{\rm{\Omega }}$ over angles ${\theta }_{m}$. We can then perform another DFT to produce another coarse sampling of ${{dC}}_{{\rm{sca}}}/d{\rm{\Omega }}$ over a different set of angles, ${\theta }_{m+\delta }$. We have now doubled the angular resolution of ${{dC}}_{{\rm{sca}}}/d{\rm{\Omega }}$ by performing two unpadded DFTs. This process takes advantage of the following:

$$\begin{eqnarray}&&{\hat{f}}_{m+\delta }=\displaystyle \sum _{n=0}^{N}{e}^{-2\pi i(m+\delta )\frac{n}{N}}{f}_{n}\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&=\;\displaystyle \sum _{n=0}^{N}{e}^{-2\pi {im}\frac{n}{N}}\underset{{F}_{n}}{\underbrace{\left({e}^{-2\pi i\delta \frac{n}{N}}{f}_{n}\right)}}\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&=\;{\widehat{F}}_{m}.\end{eqnarray} \tag{ 39 }$$

By doing offset DFTs of g for as many δ values as necessary to achieve the desired resolution in θ, we avoid the need for any padding. The complexity of this entire algorithm is ${ \mathcal O }({N}_{{\rm{FFT}}}N\mathrm{log}\;K+\beta {N}^{2})$, where the $\beta {N}^{2}$ term comes from the calculation of g(x) (scaled by some constant $\beta ={ \mathcal O }(1)$), and ${N}_{{\rm{FFT}}}=({\theta }_{{\rm{max}}}/\delta {\theta }_{{\rm{sca}}})/K$ is the number of FFTs necessary to achieve the desired resolution in the scattering angle, $\delta \theta$.

The naive calculation of $S(\theta )$, using the full two-dimensional FFT with a padded grid, has a computational time ${t}_{{\rm{naive}}}\propto {N}_{{\rm{pad}}}^{2}\mathrm{log}\;{N}_{{\rm{pad}}}$, where ${N}_{{\rm{pad}}}=2\pi /(k{\rm{\Delta }}x\delta {\theta }_{{\rm{sca}}})$. Compare this to the unpadded grid, where $N\approx 2a/{\rm{\Delta }}x$, where a is the characteristic radius of the grain. We have that

$$\begin{eqnarray}&&{t}_{{\rm{naive}}}\propto {\left(\displaystyle \frac{N\pi }{{ka}\delta {\theta }_{{\rm{sca}}}}\right)}^{2}\mathrm{log}\left(\displaystyle \frac{N\pi }{{ka}\delta {\theta }_{{\rm{sca}}}}\right).\end{eqnarray} \tag{ 40 }$$

Utilizing multiple FFTs $[{N}_{{\rm{FFT}}}=({\theta }_{{\rm{max}}}/\delta {\theta }_{{\rm{sca}}})/K{\rm{]}}$ of length $K={Nk}{\rm{\Delta }}x({\theta }_{{\rm{max}}}/2\pi )=x({\theta }_{{\rm{max}}}/\pi )$, requires a computational time that scales as

$$\begin{eqnarray}&&{t}_{{\rm{GGADT}}}\propto \left(\displaystyle \frac{\pi }{{ka}\delta {\theta }_{{\rm{sca}}}}\right)\mathrm{log}\left({ka}\displaystyle \frac{{\theta }_{{\rm{max}}}}{2\pi }\right)+\beta {N}^{2}.\end{eqnarray} \tag{ 41 }$$

We have turned an ${N}^{2}\mathrm{log}\;N$ problem into an N² problem, and, as shown in Figure 16, GGADT is approximately two orders of magnitude faster than calculations that use a padded two-dimensional FFT.

Zoom In Zoom Out Reset image size

Another key advantage of GGADT over naive padding methods lies in memory requirements. Padded FFT's require storing a two-dimensional array of floats, with dimension ${N}_{\mathrm{pad}}\times {N}_{\mathrm{pad}}$. This corresponds to a memory requirement of

$$\begin{eqnarray}{M}_{{\rm{naive}}} & = & 1.37\;{\rm{GB}}\left(\displaystyle \frac{{M}_{0}}{8\;{\rm{bytes}}}\right){\left(\displaystyle \frac{100^{\prime\prime} }{\delta {\theta }_{{\rm{sca}}}}\right)}^{2}\\ & & \times {\left(\displaystyle \frac{500{\rm{eV}}}{E}\right)}^{2}{\left(\displaystyle \frac{0.1\mu {\rm{m}}}{a}\right)}^{2}{\left(\displaystyle \frac{N}{512}\right)}^{2},\end{eqnarray} \tag{ 42 }$$

where ${M}_{{\rm{naive}}}$ is the total memory required to store the array, and M₀ is the memory requirement for a complex number (we assume 8 bytes, the standard length of a single precision complex number). For a fiducial case of 500 eV, $a=0.1\;\mu {\rm{m}}$, $\delta \theta =100^{\prime\prime}$, and N = 512, we have that ${M}_{{\rm{naive}}}\approx 1.37\mathrm{GB}$, a hefty memory requirement even for computers at the time of this writing. GGADT, on the other hand, requires a mere ${M}_{0}{N}^{2}=2\times {10}^{6}{(N/512)}^{2}$ bytes. It is possible to improve this to $M={M}_{0}N\approx 4\times {10}^{3}(N/512)$ bytes by avoiding storing the entire two-dimensional shadow function.