Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.mrao.cam.ac.uk/~bn204/mk2/publications/2008/dustgal1.pdf Äàòà èçìåíåíèÿ: Wed Nov 25 23:11:21 2009 Äàòà èíäåêñèðîâàíèÿ: Tue Oct 2 19:09:24 2012 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: spitzer space telescope

Mon. Not. R. Astron. Soc. 000, 1­16; (2008)

Printed 12 May 2008

A (MN L TEX style file v2.2)

A Semi-Empirical Model of the Infrared Emission from Galaxies
D. C. Ford, B. Nikolic and P. Alexander
Astrophysics Group, Cavendish Laboratory, J J Thomson Avenue, Cambridge CB3 0HE, UK

arXiv:0805.1732v1 [astro-ph] 12 May 2008

Accepted 2008 April 17. Received 2008 March 19; in original form 2007 July 14

ABSTRACT

We present a semi-empirical model for the infrared emission of dust around star-forming sites in galaxies. Our approach combines a simple model of radiative transfer in dust clouds with a state-of-the-art model of the microscopic optical properties of dust grains pioneered by Draine & Li. In combination with the S TA R B U R S T 9 9 stellar spectral synthesis package, this framework is able to produce synthetic spectra for galaxies which extend from the Lyman limit through to the far-infrared. We use it to probe how model galaxy spectra depend upon the physical characteristics of their dust grain populations, and upon the energy sources which heat that dust. We compare the predictions of our model with the 8- and 24-µm luminosities of sources in the Spitzer First Look Survey, and conclude by using the models to analyse the relative merits of various colour diagnostics in distinguishing systems out to a redshift of 2 with ongoing star formation from those with only old stellar populations. Key words: infrared: galaxies ­ infrared: ISM ­ dust, extinction ­ radiative transfer ­ galaxies: starburst ­ galaxies: high-redshift

1 INTRODUCTION It is increasingly clear that the infrared spectra of galaxies hold vital clues concerning galaxy energetics. Observations of our nearest neighbours tell us that around 60 per cent of their star formation is obscured by dust at visible wavelengths (Takeuchi et al. 2006), and the cosmic infrared background (CIB) indicates by its brightness that this is also true of sources at cosmological redshifts (Hauser & Dwek 2001). Furthermore, deep sub-millimetre surveys (Smail et al. 1997; Hughes et al. 1998; Smail et al. 1998; Blain et al. 1999) have revealed a large population of ultra-luminous infrared galaxies (ULIRGs) at z 2 (Blain et al. 2002) ­ clearly very dusty systems, and, if not harbouring active nuclei, very actively star-forming also (Smail et al. 1997). Given the degree of optical extinction in these systems (Swinbank et al. 2004), it is apparent that UV/optical studies of the cosmic star formation history are subject to substantial incompleteness. The advent of the Spitzer Space Telescope has allowed a much greater understanding of the sources which comprise the CIB. Its resolution and sensitivity has allowed more than two million infrared galaxies to be resolved in 49 deg2 of sky by the Spitzer WideArea Infrared Extragalactic (SWIRE) Survey (Lonsdale et al. 2003, 2004). Moreover, Spitzer's wavelength coverage, 3.6­160 µm, encompasses three emission regimes in the spectra of normal galaxies, each yielding information complementary to the others. In the rest-frame far-infrared (FIR) ­ here taken to extend from around 30 to 300 µm ­ thermal emission from large dust grains dominates. In the rest-frame mid-infrared (MIR) ­ here taken to extend from around 4 to 30 µm ­ emission from transiently-heated dust grains


dominates, marked by a series of broad emission features (see, e.g., Draine 2003, and references therein), the catalogue of which has recently been greatly expanded by Spitzer (Draine & Li 2007; Smith et al. 2007). These are attributed to polycyclic aromatic hydrocarbon (PAH) molecules, and so we shall refer to them as `PAH features'1 . Finally, in the rest-frame near-infrared (NIR) ­ here taken to refer to 4 µm ­ stellar emission dominates. This wealth of available information has motivated many studies which have sought to provide a framework in which this emission can be interpreted. Some of these take an empirical approach, matching unresolved sources to template spectra derived from a variety of local galaxies (Rowan-Robinson & Crawford 1989; Xu et al. 2001; Rowan-Robinson et al. 2004, 2005). These yield fast diagnostics which are readily applicable to large numbers of sources. But they provide little information about the physical processes which fundamentally shape the spectra. Others have sought to develop semi-empirical models of infrared spectra, considering the propagation of radiation through dusty media (e.g. Silva et al. 1998; Efstathiou et al. 2000; Takagi et al. 2003; Efstathiou & Rowan-Robinson 1995). This task can be split into two parts: modelling the optical properties of individual dust grains, and modelling the large-scale transport of radiation through some realistic dust geometry. Some authors (e.g. Fritz et al. 2006; Piovan et al. 2006a,b) have incorporated a detailed consideration of the radiative transport, including a treatment not only of absorption, but also of elastic photon scattering, which becomes the principal source of complexity. Others (e.g. Li & Draine

Email: dcf21@mrao.cam.ac.uk

In the literature, they are also commonly referred to as `Aromatic Infrared Bands' (AIBs), or, historically, as `Unidentified Infrared Bands' (UIBs)

1


2
2001; Draine & Li 2007) have developed highly sophisticated models of the optical properties of individual dust grains, but simplify radiative transfer. Essentially a computational trade-off exists: studies which include a detailed treatment of radiative transfer involve many evaluations of dust emissivities in differing environments. This is prohibitively time consuming with state-of-the-art models of the microscopic optical properties of individual dust grains and so simpler alternatives must be sought. In this paper, we use a state-of-the-art model for the individual dust grains ­ that of Draine & Li (2001) and Li & Draine (2001) ­ and then seek to make a minimal set of simplifications in our treatment of the radiative transfer problem such that it becomes computationally viable. In Section 2 we outline the dust geometries which we consider. In Section 3 we go on to present our modelling of the radiative transfer, and in Section 4 we construct UV­visible radiation fields appropriate for dust-heating by star-forming galaxies. In Sections 5 and 6 we draw upon models of dust grain populations and their microscopic optical properties from the literature. In Section 7 we conclude the development of our model with a simple framework for modelling the evolution of the metallicity and mass of gas in passively-evolving star-forming gas clouds. In Sections 8 and 9 we present the basic predictions of our model. In Section 10 we use the model to predict 8- and 24-µm luminosities for a sample of galaxies and compare to observations. Finally, in Section 11, we develop a simple model for the evolving colours of early- and latetype galaxies. Wherever required, we assume a flat CDM cosmology with H0 = 72 km s-1 Mpc-1 and = 0.7.
Dust Shell H number density nH (r) Grain number density ngr (r)

Ray Y

s Heat Source

r

0

X r

r

1

Figure 1. Dust in a spherical shell of dust around a single isotropic heating source.

region, essentially devoid of dust due to the sublimation of grains by energetic photons.

3 RADIATIVE TRANSFER 2 THE MODEL GEOMETRY In this paper we consider two basic geometries for the spatial distribution of dust and the source of illumination. The first geometry is a shell of dust grains ­ which we shall term a `circumnuclear' grain population ­ surrounding a point-like heating source, representing dust around a star-forming region. The second is a uniform distribution of dust grains ­ which we shall term a `diffuse' grain population ­ spread throughout a diffuse inter-stellar medium (ISM) within which the radiation field is assumed spatially uniform. In a future paper, we shall use them as components of a composite model for star-forming galaxies. The geometry adopted for our diffuse populations is the simpler. We assume the dust-bearing ISM to be optically thin at infrared wavelengths, such that the re-absorption of dust emission can be neglected. We further assume the UV­visible heating radiation field within it to be spatially uniform, hence all grains of any given size and composition have the same emissivity. Under these conditions, total dust emission is directly proportional to the number of grains present, and we therefore scale all quantities per unit volume of ISM. The circumnuclear geometry is illustrated in Figure 1. A point-like heating source lies at the centre of a spherical shell of dust, of inner radius r0 and outer radius r1 . Within the shell, we trace the density of the medium via the number density of hydrogen nuclei, assumed to be spherically symmetric and denoted nH (r), such that the column density Nc of hydrogen nuclei along a line of sight passing through a dust shell to its nucleus is given by:
r

We adopt a highly-simplified treatment of the radiative transfer which we argue is sufficiently accurate for a wide range of problems.

3.1 The circumnuclear geometry The evolution of the surface brightness I along ray Y in Figure 1 is governed by the time-independent equation of radiative transfer: dI ds (r)I + (r)nH (r) + nH (r) 4
I C

(2) where the first term on the right-hand side describes the absorption of radiation by dust, the second dust emission, and the third the scattering of photons into the ray. The integral in the third term is over solid angle at X ; is the scattering angle between d and the direction of Y . As mentioned above, we use nH (r) to parameterise the varying spatial density of material. C ,abs and C ,sca are the cross sections to absorption and scattering respectively, expressed per hydrogen atom, and averaged as described below over the compositions and sizes of particles within the grain population. (r) is the emissivity of the grain population ­ the power emitted per unit frequency into unit solid angle, normalised in the same way as for the cross sections. s measures distance along the ray. The extinction cross section, C ,ext , is the sum of the absorption and scattering cross sections: C
,ext

X

= -C

,ext nH

,sca

( ) d ,

Nc =

1

r

nH (r) d r.

(1)
II

0

= C

,abs

+ C

,sca

.

(3)

Astrophysically, the cavity at r < r0 might correspond to an H

The averaging of the quantities above over grains of varying com-


Infrared Spectra of Galaxies
positions and sizes is performed as follows: C
,j

3

= =


i

a

i C , j (a) i (a, r)

1 dni r (a) g nH da 1 nH dni r g (a) da

through the sphere of constant radius passing through X is related to H (r) via: (4) F (r) = 4 H (r). Equation (6) can thus be re-written: 1 r2 F (r) = nH (r) 4 (r) - C r2 r
,abs F

da, da,

(10)

(r)


i

(5)

a

where j {ext, abs, sca}, i denotes a population of grains of given composition, and ni r (a) denotes the spatial number density of g grains of each composition with radius smaller than a. The first simplification that we make is to neglect scattering, such that C ,sca = 0. We discuss the validity of this assumption in the following section. Integrating the remaining terms of Equation (2) over all rays passing through the point X in Figure 1, we obtain (Chandresekhar 1960; Rowan-Robinson 1980): 1 r2 H (r) = nH (r) - C r2 r where: J (r) = and: H (r) = 1 2
1 -1 ,abs J

(r) .

(11)

This equation is integrated numerically from the inner to the outer radius.

3.2 Assumptions made in the circumnuclear geometry The assumptions made in the previous section ­ i.e. the neglect of scattering and the radial beaming of dust emission ­ will have negligible effect upon the predictions of our model for column densities of dust which are optically thin at all wavelengths, i.e. for Nc 1023 H m-2 . For column densities of dust which are optically thick in the UV, but not in the infrared, our neglect of scattering will lead us to under-estimate the path lengths of UV/optical photons through the dust shell by a factor of 1­2, and so to under-estimate the absorption of UV/optical radiation by a similar factor. Since this effect is essentially the same as that of reducing the column density of dust, the effect when using these models to fit the spectral shape of real sources will be that we will over-estimate the dust masses of these objects. Our assumption that dust emission is beamed radially outwards only begins to fail for dust shells with higher column densities still, when they become optically thick even at infrared wavelengths. As the dust emission is assumed to take the shortest path out of the dust shell, we will under-estimate its re-absorption in these optically thick cases. In practice, this effect becomes significant for dust column densities 1026 H m-2 , as will be shown in Figure 11. Rowan-Robinson (1980) studied the effects of a similar set of assumptions in their calculation of radiative transport in hot-centred star-forming clouds, and for the range of systems they consider, they report errors of around 10 per cent. In addition to the two assumptions just discussed, it is also apparent that the adopted geometry is simplistic, having only one single heat source. We note, however, that this geometry is observationally indistinguishable from an ensemble of N smaller circumnuclear geometries, each heated by its own central heat source with luminosity scaled by a factor 1/N with respect to the single shell, and each co aining a dust shell with inner and outer radii scaled nt by factor 1/ N and dust density distribution n (r) scaled according to: n (r) = N n(r N ) (12) where n(r) is the density distribution of the single shell. The column density of dust around each heat source in this latter ensemble is the same as that in the former single circumnuclear shell; the radiation field incident upon grains on the inner edge of each dust shell is the same; the total mass of dust in the two cases is the same; and to a remote observer, the total solid angle subtended by the dust in the two cases is the same. In summary, although our circumnuclear model is nominally of dust around a single heat source, it is also a good model of systems where that luminosity production is distributed between several discrete sources.

(r) ,

(6)

1 2

1 -1

I dµ ,

µ = cos ,

(7)

I µ dµ ,

(8)

with denoting the angle made between each ray and the radial direction, as shown in Figure 1. Geometrically, J (r) may be visualised as the average surface brightness along all rays passing through X , averaged over 4 steradians. H (r) may similarly be visualised as the average projected onto the radial direction. Rowan-Robinson (1980) introduced what has become a widely-used decomposition of the surface brightness I along each ray into three components, depending upon where photons last interacted with matter: I = I + I + I ,
(1) I (1) (2) (3)

(9)

is radiation from the central heat source which has not where (2) (3) been absorbed by dust, I is radiation emitted by dust and I is scattered radiation ­ which we have already neglected. While we do not use this in our mathematical treatment of Equation (2), it is useful in our discussion presently. The relationship between J (r) and H (r) encodes the angular distribution of the radiation flux passing through X . In the limiting case of a radiation field propagating exclusively in the radial direction, H (r) = J (r). In the opposite limit of an isotropic radiation field, H (r) = 1/2J (r). Given a point heat source, as in Figure 1, (1) the former limit is applicable to the component I ; the radiation (2) field emanating from the heat source is purely radial. For I , however, H (r) < J (r). Our second, and final, assumption, is that J (r) = H (r) in Equation (6). For UV­visible wavelengths, this assumption holds because the heating radiation field is expected to dominate over (1) dust emission at these wavelengths, and so I I . In the infrared, (2) where I is expected to dominate I , J (r) is under-predicted, but the assumption continues to hold if (r) C ,abs J (r), that is to say, if the dust emission from the shell is not appreciably reabsorbed. Geometrically, this assumption is equivalent to assuming that dust emission is beamed along the outward radial direction. Finally, we note that the net outward flux F (r) of radiation


4
3.3 The diffuse geometry Our treatment of radiative transfer in diffuse dust grain populations is simpler than the above. Instead of having a heat source of luminosity L , we have an interstellar radiation field (ISRF), whose (1) energy density E we normalise to 0 , where 0 is that of the solar neighbourhood interstellar radiation field less the cosmic microwave background (CMB)2 , and is a free parameter. We adopt 0 = 7.46 â 10-14 J m-3 , derived from integration of the ISRF of Mathis et al. (1983) and Mezger et al. (1982). The total luminosity emerging from the model can then be written: L = 4 nHV + E
(1)

stellar population of age t with our star formation history:
t

L (t ) =

0

(t ) 106 M

L

SSP

(t - t ) dt .

(14)

5 THE DUST MODEL Whilst there exist fairly tight observational constraints on the composition and size distribution of dust grains in the Milky Way, relatively little is known about those in other galaxies (Draine 2003). In this paper, we therefore base our dust grain population upon that inferred from observation of our own galaxy. Following Li & Draine (2001, hereafter, LD01), we consider binary populations of `carbonaceous' and silicate grains. The former sub-population includes both PAH molecules and larger graphitic grains; the optical properties of these grains exhibit a smooth transition with grain radius, centred around a radius of a . The absorption cross sections of carbonaceous grains of radius a is taken to be:
cr Caabs (a) = , PAH P (a)CAHs (a) + [1 - ,ab PAH PAH

cA

4

(13)

where V is the volume of the dust-bearing ISM, and A its surface area through which the interstellar radiation field leaves the galaxy.

4 THE HEATING RADIATION FIELD In this paper, we consider models for the heating radiation field due to star formation. We generate these using Version 5.1 of the S TA R B U R S T 9 9 stellar spectral synthesis package (Leitherer et al. 1999; Vazquez & Leitherer 2005). This package offers the facil´ ity to model stellar populations with two classes of star formation histories (SFHs). The first is ongoing star formation, proceeding at constant star-formation rate (SFR) , which began at some time t previously; we hereafter term these models `continuous' SFHs. The second is a delta-function SFH, representing an instantaneous burst of star formation at some time t previously, of total mass m; we term this the `instantaneous' SFH. The stellar populations modelled for the instantaneous SFHs may be referred to as single stellar populations (SSPs), which is to say that all of the stars within them are coeval. S TA R B U R S T 9 9 models both stellar emission and also nebular continuum. For the former it uses the stellar evolution tracks of the Padova group (Fagotto et al. 1994), with the addition of tracks for thermally pulsing asymptotic giant branch (TP-AGB) stars to improve the accuracy of the modelling of low and intermediate mass stars.3 For the nebular continuum, the emission coefficients of Ferland (1980) are used. This is the source of a problem in S TA R B U R S T 9 94 : the data of Ferland (1980) do not extend beyond 4.5 µm, however they are extrapolated to 160 µm, yielding large, unphysical, infrared luminosities. For 4 µm, we use a more physical extrapolation, taking L -0.1 . In addition to the continuous and instantaneous SFHs modelled by S TA R B U R S T 9 9, we have also considered arbitrary SFHs, modelled by convolving the luminosity, LSSP (t ), of a 106 M single

(a)]C

gra ,abs

(a),

(15)

where the weighting parameter

(a) is given by: (16)



PAH

(a) = 1 - qgra â min 1, (a /a)3 ,

and qgra = 0.01 sets even the smallest PAH molecules to exhibit 1 per cent of the continuum absorption of graphitic grains. The tran° sition radius, a , is set by default to 50 A. For the size-distribution of grains in each of these populations, we follow the parametric forms used by Weingartner & Draine (2001); for the carbonaceous grains, we use: 1 nH Z Z 1, exp - a - at,g /ac
,g 3

dnca gr da

r

= D(a) +

Cg a

a at,g ° 3.5 A



g

F (a, g , at,g ) (17) a a < at,g

â

< <

,

at,g

where D(a) represents two log-normal peaks:

D(a) =

i=1



2

Bi 1 ln a/a0,i exp - a 2

,

(18)

which were introduced by LD01 to reproduce the mid-infrared luminosities observed by ISO, the Infrared Telescope in Space (IRTS) and by IRAS at 60 µm. The term F (a, g , at,g ) provides curvature: F (a, g , at ) = 1 + a/at (1 - a/at )-1

0 < 0,

(19)

We neglect the CMB in this normalisation because, in contrast to the starlight component of the ISRF, it would make no sense to enhance it by a factor . It should be noted that the CMB is also absent from all models presented in this paper. 3 It should be noted that this is a significant departure from previous versions of S TA R B U R S T 9 9, which used the stellar evolution tracks of the Geneva group, and did not model low mass stars, introducing serious errors in the modelling of old stellar populations. 4 See notice by Hunt, October 30, 2006, in the Knowledge Base of the S TA R B U R S T 9 9 website.

2

Z is the mass ratio of metals, Z = 0.02 is the solar mass ratio of metals, and all other symbols are as defined in Weingartner & Draine (2001). The normalisation constants Bi are given by: 3 exp -4.5 a3,i (2 )3/2 0
2

Bi =

bC,i m 1 + erf
3 2

C

+

° ln a0,i /3.5 A 2

(20)

where mC = 1.99 â 10-26 kg is the mass of a carbon atom, = ° 2.24 â 103 kg m-3 is the density of graphite, a0,1 = 3.5 A and


Infrared Spectra of Galaxies
Parameter Value -1.54 -0.165 0.0107 µm 0.428 µm 9.99 â 10- -2.21 0.300 0.164 µm 1.00 â 10- 0.75bC 0.25bC 6.0 â 10-5

5

100 Carbonaceous grains 10
12
3

1029

g g at,g ac,g Cg s s at,s Cs bC,1 bC,2 bC

13

-1 H

a4 / cm n

1

100 Silicate grains 10

Table 1. The parameters of the default size distribution which we adopt ­ the Weingartner & Draine (2001) preferred distribution for RV = 3.1 Milky Way sight lines.

dngr da

1

° a0,2 = 30 A are the wavelengths of the centres of the two log-normal peaks, and = 0.4. For the silicate grains, we use: 1 nH 1, exp - [(a - at,s ) /ac,s ]
3

10-9

10-8 a/m

10

-7

10

-6

Z Z

dnsirl g da

Cs = a

a at,s



s

F ( ; s , at,s )

(21)

â

° 3.5 A , at,s

Figure 2. The adopted size distributions for silicate (top) and carbonaceous (bottom) grains, assuming a solar metallicity environment. Unit areas under each distribution represent unit masses of grain material.

< <

a a

<

at,s of neutral and ionised PAH molecules with a weighting parameter f describing the ionisation fraction. We take this to have a default value of 80 per cent, matching that which Draine & Li (2001) find in their model fits to Galactic photo-dissociation regions. The recr sulting absorption cross section Caabs (a) is shown for a range of , s grain radii in Figure 3(c), and the absorption cross section Ci,labs (a) of the silicate grain population Figure 3(d).

By default, we set the parameters of these size distributions to those preferred by Weingartner & Draine (2001) for RV = 3.1 Milky Way sight lines, as given in Table 1. The effect of using instead the size distributions preferred by those authors for RV = 4.0 and RV = 5.5 Milky Way sight lines will be discussed in Section 8. The appearance of Z in Equation (17) is a departure from Weingartner & Draine (2001), who only considered Galactic environments. Thus, their size distributions are normalised to a dust-togas-mass ratio appropriate for solar metallicity environments. As the variation in this ratio with Z is quite poorly understood, we make the assumption that it is linearly proportional to Z , which is implicit in our renormalisation above. These distributions are plotted in Figure 2. To calculate the absorption cross sections of silicate and graphitic grains, we follow LD01 and use dielectric functions for these species (Draine & Lee 1984) and Mie theory (see, e.g., Bohren & Huffman 1998) to estimate the absorption cross sections of spherical particles of radius a. The treatment of graphitic grains is slightly complicated by the anisotropy of graphite's dielectric function. We follow LD01 in calculating an averaged absorption cross section using the `1/3-2/3 approximation' (Draine & Malhotra 1993). For the PAH molecules, LD01 give algebraic fits to terrestrial P laboratory measurements of CAHs (a) for neutral and ionised sam,ab ples. Draine & Li (2007, hereafter, DL07) revise these in the light of new near-infrared data (Mattioda et al. 2005b), and in order to fit the high-fidelity spectra of nearby star-forming galaxies observed by the Spitzer Infrared Nearby Galaxies Survey (SINGS; Kennicutt et al. 2003) project. We implement the cross sections given by both LD01 and DL07, which are shown in Figures 3(a) and 3(b) for ° neutral and ionised grains respectively, both of radius 5 A. In the remainder of this paper, we use the DL07 cross sections throughout, except in Section 10. In both cases, we average the cross sections

6 MODELLING THE EMISSIVITY OF THE DUST
i In this section, we outline how we model the emissivities (a, r) of dust grains as a function of the energy density E (r) of radiation to which they are subjected.

6.1 Transiently-heated grains Modelling emission from transiently-heated grains requires the calculation of the time-averaged probability distributions P(E ) for their internal energies. Exact treatment of this problem would require knowledge of all of their vibrational energy levels and transition probabilities. Our simplified analysis follows Draine & Li (2001). We use Debye models for the normal modes of the C/Si skeletons of PAH and silicate particles, with Debye temperatures as used by Draine & Li (2001). We use Einstein models for the stretching, in-plane bending, and out-of-plane bending modes of the peripheral C­H bonds of PAH molecules: we assume the modes of each bond to be quantum harmonic oscillators with the same fundamental frequencies, as given in Draine & Li (2001). To reduce the resulting mode spectra to a computationally tractable number of energy states, we follow Guhathakurta & Draine (1989) and Draine & Li (2001) in dividing them into Nbin bins (for our choice of bins, see Appendix A), with mean energies


6
(a) 0.01 (b)

Cabs / a2

10-4

10-6 Neutral PAHs; LD01 model Neutral PAHs; DL07 model C ; a = 5° A Ionised PAHs; LD01 model Ionised PAHs; DL07 model C ; a = 5° A

1

(c)

(d)

0.01 Cabs / a2

10-4

10-6

a = 5 â 10-10 m a = 1 â 10-8 m C -6 a = 1 â 10 m 0.1 1 /µm 10 100 0.1

a = 5 â 10 a = 1 â 10 a = 1 â 10

-10 -8 -6

m m Si m 1 /µm 10 100

Figure 3. The adopted absorption cross sections. Panels (a) and (b) show those for neutral and ionised PAH grains respectively, as given by LD01 and DL07 ° for grains of radius 5 A. A discussion of the near-infrared feature introduced into the cross sections of ionised PAHs by DL07 at 1.05 µm , and the negative feature at 1.905 µm , can be found in Mattioda et al. (2005a). Panels (c) and (d) show those for carbonaceous and silicate grains respectively for a variety of grain radii a, assuming a PAH ionisation fraction f = 0.8. Each trace is normalised with respect to the classical grain cross section of a2 .

U i , widths U i , and time-averaged occupation probabilities Pi . We denote as T ji the transition rate between bins i and j. The time evolution of Pi is then given by: dPi = dt

constraint:


i

Pi = 1,

(25)



j=i

Ti j P j -



T j i Pi .

(22)

j=i

The time-averaged steady-state probability distribution which we seek is that to which the above converges over time, and for which dPi /dt = 0. The elements of T ji with j > i describe the upward transitions of grains that result from photon absorption; we model these using equations (15­25) of Draine & Li (2001); in Appendix B we reproduce these relations and describe a numerical optimisation that we use in their calculation. The elements with j < i describe the radiative cooling of grains; here we use the `thermal continuous' approximation (equation 41 of Draine & Li 2001, reproduced here as Equation B4), which models the cooling of grains as a continuous process, where each state i only makes downward transitions to the adjacent state i - 1. This allows much faster solution of Equation 22 to find Pi . The diagonal terms are chosen (Draine & Li 2001) so that: Tii = -

using the method of Guhathakurta & Draine (1989). Given the vector Pi , we calculate the time-averaged emissivity of each grain usi ing the thermal approximation, under which (a, r) can be calculated using equation (56) of Draine & Li (2001):

=

2h c2

3


i

Pi , exp(h /k i ) - 1

(26)

where i is the characteristic temperature of bin i, as defined in Draine & Li (2001), the sum is over all bins i whose central energies are greater than h , and we have neglected the factor (1 + 3 uE /8 ) shown by those authors; this represents stimulated emission and may straightforwardly be shown to be negligible in all of the models presented in this paper.

6.2 Large grains For sufficiently large grains, the approach outlined above becomes inefficient. Their internal energies become much larger than the energies of the photons they absorb, and so their temperature fluctuations are not significant. Their internal energy probability distributions P(E ) tend towards delta functions (Li & Draine 2001). In this limit, we can model the energetics and emission of these grains by numerically solving the equation of radiative balance to find their equilibrium temperatures:
0



T ji ,

(23)

j=i

hence Equation (22) can be re-written in the form: dPi = dt
Nb

j=0



in

Ti j P j = 0.

(24)

These equations are solved, subject to the additional normalisation

C

,abs

(a)cu d =

0

4C

,abs

(a)B (T ) d ,

(27)


Infrared Spectra of Galaxies
where c is the speed of light, and B (T ) the Planck function at temperature T . Li & Draine (2001) show that when grains are bathed in the local ISRF of the Solar Neighbourhood, this continuum ap° proximation is valid for grain radii greater than 250 A. We adopt this same transition radius for the range of models considered here. The emissivities of these grains are calculated by assuming them to radiate as modified blackbodies with a single characteristic temperature:
Mass of gas / 109 M Lookback redshift, z 10 6 (a) Constant Z Model 4 2 1.5 1 2 0.5 6 (b) Perfect Mixing Model 4 m Z 2 0.5 0 0.4
-1 c

7

7

5

4

3

2

1

= C

,abs B

(T ).

(28)

Metallicity / Z


0 2 1.5 1

7 EVOLUTION OF GAS MASS AND METALLICITY To follow the evolution of a star-forming system we use a simple closed box model. We denote the total stellar mass-loss rate (due to stellar winds and supernovae) as (t ), with metallicity Zoutflow . The total mass of gas and dust in a cloud is mc (t ), of which metals comprise a mass mz (t ), such that Z (t ) = mz (t )/mc (t ). Thus we may write: dmc = -(t ) + (t ). dt (29)

(c) Star Formation Rate

We consider two models for the mixing of metal-enriched material in the cloud. In the simpler, we assume a time-invariant metallicity, such that mz (t ) = Z mcloud (t ). We call this model the `constant Z ' mixing model. In the second, we assume perfect mixing, such that the metallicity of the material which goes into forming new stars is representative of that of the whole cloud, in which case: dmz = -Z (t )(t ) + Z dt
outflow

0.3 0.2 0.1 0.0

(t )(t ).

(30)

(t)/M yr

0.5

1

2 Time / Gyr

5

10 13

We call this the `perfect mixing' model. To calculate (t ) and Zoutflow (t ) we used the same S TA R B U R S T 9 9 models as described in Section 4. The evolving metallicity of the gas forming into new stars could not be treated smoothly, as the Padova stellar-evolution tracks used by S TA R B U R S T 9 9 are only available for five stellar metallicities: Z = 0.0004, 0.004, 0.008, 0.02 and 0.05. Instead, we modelled each SSP within our calculation of the SFH using the track whose metallicity was logarithmically closest to that desired. As a simple starting point, we model two dwarf-galaxysized clouds. Firstly, we model a system composed initially of 5 â 109 M of gas and dust of metallicity 0.05, but with no stars. This cloud begins to undergo star formation at a rate of 0.3 M yr-1 at a redshift5 of z = 10, and maintains this rate of star formation until the current epoch. We show the evolution of the mass and metallicity of this system in Figure 4. Secondly, we form a model of an early-type galaxy. We assume a system of total mass 5 â 109 M , composed initially of zero-metallicity gas. We assume it to undergo star formation at a constant rate in the redshift range 10­5, such that 90 per cent of this material is converted into stars over this period. The evolution of the mass and metallicity of this system is shown in Figure 5. We note that our perfect mixing model predicts such galaxies to have decreasing metallicities with time (Figure 5b). After z 5, mass loss from low-mass stars dominate the return of material to the ISM. This material has lower metallicity than supernova ejecta.

Figure 4. A model of the ISM of a galaxy of total mass 5 â 109 M , initially of metallicity 0.05 solar, but containing no stars. It begins to form stars at a rate of 0.3 M yr-1 at z = 10, and maintains this rate of star formation until the current epoch. mc is the total mass of gas and dust in the galaxy's ISM. Panels (a) and (b) represent differing models of the mixing of material in the ISM, as described in the text.

8 RESULTS: VARIATIONS IN DUST EMISSION WITH PHYSICAL CONDITIONS In this section, we present model results using the diffuse geometry described in Section 3.3. The simple optically-thin approximation allows us to investigate how the emissivity of dust depends upon the physical properties of the grain population, and the radiation field which heats the dust. Throughout this section, we take for our interstellar heating radiation field a S TA R B U R S T 9 9 model of a continuous SFH of age 1 Gyr, and assume the ISM to be contained within a homogeneous spherical volume of radius r0 . This leaves three free parameters: nH , and r0 . We choose, however, to re-parameterise the problem in terms of more physically interesting quantities ­ the column density of hydrogen nuclei along a diameter of the ISM, Nc , the total star-formation rate of the heating radiation field, , and ­ by means of the relations: r0 = L ,0 , 0 c (31)

5

See Section 1 for a definition of our adopted cosmology.


8
Lookback redshift, z 10 7 5 4 3 2 1
1022

(a) Constant Z Model 1

1.5 1 0.1 Mass of gas / 109 M 0.5 0.01 (b) Perfect Mixing Model 1 m m Z 0.1 0.5 0.01 8
-1

L /W Hz

-1

m 1018 Total emission PAH emission Graphite emission Silicate emission 2 5 10 20 /µm 50 100 200 1016

0 2 1.5
c Z

1

Figure 6. Model spectrum of the emission from dust in our diffuse geometry, neglecting stellar emission. We assume an ISRF appropriate for = 1, Nc = 1024 H m-2 , = 1 M yr-1 , for a system of age 1 Gyr. We decompose the total emission into the contributions from grains composed of graphite, PAHs and silicates.

-3

2

1020

Metallicity / Z


0 (c) Star Formation Rate
-1

2 â 1022 1022

= 0.2 = 1.0 = 10.0

(t)/M yr

4 2 0 0.5 1 2 Time / Gyr 5 10 13

L /W Hz

6

5 â 1021

2 â 1021 1021

5 â 1020 2 5 10 20 /µm 50 100 200

Figure 5. A model of the ISM of an early-type galaxy of total mass 5 â 109 M , which converts 90 per cent of this material into stars in the redshift range 10­5. mc is the total mass of gas and dust in the galaxy's ISM and mz the mass of metals. Panels (a) and (b) represent differing models of the mixing of material in the ISM, as described in the text.

Figure 7. Spectra modelled using our diffuse geometry with Nc = 1024 H m-2 ; we plot the sum of stellar and dust emission. We illustrate the effect of varying the intensity of the ISRF with respect to that of the solar neighbourhood.

and nH = Nc /2r0 , (32)

where L ,0 is the specific luminosity of a S TA R B U R S T 9 9 model of unit SFR. The former relation may be derived by setting the second term on the right-hand side of Equation (13) to equal L ,0 . We truncate our adopted ISRF at an upper energy bound Emax , which we set to the ionisation energy of hydrogen, 13.6 eV; astrophysically, we would expect inter-stellar neutral hydrogen to severely attenuate the radiation field at shorter wavelengths. To avoid altering the bolometric luminosity of the ISRF as a result of this truncation, we assume that the truncated radiation is reprocessed to longer wavelengths with a spectrum matching that of the nebula continuum component of the S TA R B U R S T 9 9 model. In Figure 6 we illustrate the contributions made by grains composed of silicates, graphite, and PAH molecules to the total dust emission. The division between the contributions made by graphitic and PAH grains is somewhat blurred because of their being incorporated into a single carbonaceous grain population with smoothlyvarying optical properties. For the purposes of the decomposition shown in Figure 6, we take `graphitic' grains to be those with radii

° ° > 50 A, and `PAH molecules' to be those grains with radii < 50 A. We see that the emission shortwards of 20 µm is dominated by grains in the carbonaceous log-normal peaks, whilst the silicate grain population contributes around 70 per cent of the FIR thermal emission. The graphitic grain population is seen to play a relatively minor role in shaping the SED, in agreement with the finding of Li & Draine (2001) that the upper limit on the sizes of the carbonaceous grains in their model was poorly constrained. In Figure 7, we show the effect of varying the intensity of the ISRF with respect to that of the Solar Neighbourhood, whilst holding Nc and the star-formation rate of the galaxy constant. This means effectively altering the volume occupied by the dust and hence the energy density in the ISRF at fixed total luminosity. As decreases, the peak of the FIR thermal emission is seen to move to longer wavelengths. The temperature of the grains contributing to this emission can be estimated from peak of the grey body emission, showing that the temperature of the large grains falls from about 24 K to 13 K as varies from 10 to 0.2. In the MIR, the luminosities of the PAH features are essentially independent of . This behaviour is explained by considering the transient nature of the heating of the small grains. For the range of conditions considered, these small grains cool efficiently between photon absorption events. Thus, the MIR spectra of galax-


Infrared Spectra of Galaxies
3 â 1021

9

2 â 1021

1021

6 â 1020

f = 0% f = 50% f = 100%

4 â 1020

2

3

4

6 /µm

10

20

30

Figure 8. Spectra modelled using our diffuse geometry with Nc = 1024 H m-2 , including both stellar and dust emission. The effect of varying the ionisation fraction f from its default of 0.8 is illustrated.

ies are independent of whether a small number of grains are each undergoing frequent excitations, or a larger number of grains are each undergoing less frequent excitations. If we consider a single photon of stellar origin, the probability of its being absorbed depends only upon the integrated column density of dust it traverses in its passage through the ISM. Thus, the emission of transientlyheated grains depends only upon Nc and the form of the spectral energy distribution and hence , independent of . In Figure 8, we show the effect of changing the ionisation fraction f of the PAH grains, as defined in Section 5, from its default of 80 per cent. For a discussion of observational evidence that PAH ionisation states vary between differing astrophysical environments, see, e.g., Draine & Li (2001). We see that the MIR features in the 5 < /µm < 9 region are enhanced in ionised PAHs; those in the 10 < /µm < 20 region are enhanced when PAHs are neutral. The effect of PAH ionisation can be understood in terms of the adopted absorption cross sections shown in Figure 3. The enhancement of the luminosity of the 5 < /µm < 9 PAH features when PAHs are ionised can be explained simply from our adopP tion of a CAHs (a) which is larger at these wavelengths for ionised ,ab PAHs; the enhancement of the 3.3 µm feature when PAHs are neutral may similarly be explained. The enhancement of the features in the 10 < /µm < 20 region when PAHs are neutral, however, canP not be so explained; CAHs (a) is essentially independent of PAH ,ab ionisation state at these wavelengths. The explanation instead lies in the energetics of these grains. The UV­visible absorption cross sections of PAH molecules are essentially independent of their ionisation states; neutral and ionised grains re-process UV­visible radiation into the MIR at essentially identical rates. If neutral PAHs are less luminous than ionised PAHs in the 5 < /µm < 9 region, they must be more luminous at other wavelengths to have the same bolometric luminosities. We now examine the variation in the MIR feature strengths as a function of the physical properties of the population of PAH molecules. The size distribution of PAH grains may depend critically upon environment; in H I I regions, for example, strong fluxes of ionising photons may destroy many of the smallest grains. Such speculation is supported by the observation (see, e.g., Weingartner & Draine 2001) that a diverse set of grain size distributions are required to fit Milky Way sight lines. In Figure 9(a1) we decrease the number bC of interstellar carbon nuclei per hydrogen nucleus locked up in PAH molecules from its default value of bC = 6 â 10-5 . We meanwhile hold the powerlaw component of the grain size distribution unchanged. We see

that the result of removing some of the smallest grains from the model in this way is that the MIR luminosity decreases ­ eventually to leave only the Jeans tail of the stellar emission in its place. This ° result is as might be expected as carbonaceous grains with a < 15 A contribute in excess of 90 per cent of emission at < 10 µm. In Figure 9(b1), the effect of altering the relative proportions of the numbers of carbon nuclei composing PAH molecules in the two peaks is illustrated; as defined in Section 5, bC,1 is the propor° ° tion in the 3.5 A peak and bC,2 = 1 - bC,1 that in the 30 A peak. Our default model has bC,1 = 0.75. We see that the effect this change on dust emission is to enhance emission at 17 µm when a greater ° fraction of the carbon nuclei are placed in the 3.5 A peak, while reducing emission at 17 µm. A pivoting motion is seen around 17 µm, where the emissivities of grains in the two log-normal peaks, normalised per unit carbon atom, approximately equal one another. In Figure 9(c1), we consider the effect of varying a ­ the grain radius, as defined in Equation (15), at which the optical properties of the carbonaceous grain population shift from those of graphite to those of PAH molecules. We see that the features at ° ° < 11 µm are little altered when a is changed from 50 A to 10 A, ° but that further reducing a to 4 A does suppress those PAH features with respect to the MIR continuum. As shown in Figure 9(a1) ° above, emission at 17 µm is dominated by dust in the 3.5 A log-normal distribution, and so is only affected once a is reduced to a comparable size. We conclude that our model predictions are relatively insensitive to a , a result in agreement with Li & Draine (2001), who found it to be relatively poorly constrained by observation. Thus far we have used as our default grain size distribution that preferred by Weingartner & Draine (2001) for Milky Way sight lines with RV = 3.1. These are typically sight lines which pass through infrared cirrus, away from dense molecular clouds. However, those authors also fit grain size distributions to the dust along sight lines which pass through more dense parts of the Milky Way, characterised by larger values of RV . For each value of RV , they present a range of fits for different assumed values of bC . In Figure 10, we show the spectra produced by our model when we use the grain size distributions given by those authors for RV = 4.0 and 5.5, in each case using the value of bC which they find to best fit the sight lines studied. For these enhanced values of RV , we find that the infrared emission is reduced in line with what would be expected from the reduced number of dust grains per unit column density in these size distributions.

L /W Hz

-1

9 RESULTS: CIRCUMNUCLEAR SHELLS In this section, we present spectra modelled using our circumnuclear geometry. As in the previous section, our default heating radiation field is a stellar population with a continuous SFH, starformation rate = 1 M yr-1 and age 1 Gyr. We again truncate this at an upper energy bound of Emax = 13.6 eV and scale up the nebula continuum component of the heating radiation field accordingly. We find that the predictions of our model depend only weakly upon the geometrical thickness, r1 - r0 , of the dust shell. For example, the emission from a dust shell of inner radius r0 = 1 kpc and outer radius r1 = 3 kpc is little different from that of an infinitesimally thin shell at a radius of 2 kpc. Within the accuracy of our numerical integration, these two models are indistinguishable at 14 µm. At longer wavelengths, they do differ, the thicker


10
30 bC = 6 â 10-5

2 â 10

21

bC = 1 â 10-5

10

bC = 2 â 10-6 3

10

21

1 6 â 10 4 â 10 3 â 10 2 â 10
20

0.3
20

(a1)
20

(a2)

0.1 30

b
21

C ,1 C ,1 C ,1

= 0.90b = 0.50b = 0.10b

C C C

b b

10 1029 n
-1 4 Ha

L /W Hz

-1

3 10
21

dngr /da/cm

1 6 â 10 4 â 10 3 â 10 2 â 10
20

3

0.3
20

(b1)
20

(b2)

0.1 30

a = 4 ° A
21

A a = 10 ° a = 50 ° A

10

3 10
21

1 6 â 10 4 â 10 3 â 10
20

0.3
20

(c1)
20

(c2) 2 3 4 6 /µm 10 20 30 10
-9

0.1 10
-8

10 a/m

-7

10

-6

Figure 9. In each row, the left plot shows spectra modelled using our diffuse geometry with Nc = 1024 H m-2 , = 1 and = 1 M yr-1 . These include both stellar and dust emission, and illustrate the effect of varying the the size distribution of carbonaceous grains upon the overall dust emission. The right plot shows the size distribution of carbonaceous grains used in the plot to its left. Panel (a1) illustrates the effect of varying the strengths of the log-normal peaks in the size distribution. Panel (b1) illustrates the effect of varying the relative strengths of the two log-normal peaks. Panel (c1) shows the effect of varying the grain radius a at which the carbonaceous grain population makes its smooth transition from having the optical properties of PAH molecules to those of macroscopic graphite particles.

2 â 10 10

22

22

RV = 3.1 RV = 4.0 RV = 5.5

5 â 10 L /W Hz
-1

21

2 â 10 10

21

21

5 â 10

20

2

5

10

20 /µm

50

100

200

shell having an enhanced luminosity at 15­70 µm, enhanced by 23, 72 and 29 per cent at 20, 40 and 60 µm respectively. This can be understood in terms of two effects: an enhanced probability of multiple-photon heating of small grains close to the inner edge of the shell, and radial variations in the temperatures of the large dust grains. For the remainder of this section, we adopt a geometry with r0 = 1 kpc and r1 = 3 kpc. Figure 11 shows model spectra for dust shells with column densities Nc = 1024 and 1026 H m-2 ; we assume solar metallicity for both the stellar populations and the dust shells. For each, we show the contribution made by attenuated stellar radiation ­ i.e. (1) I in Equation (9) ­ given by: L
stellar

Figure 10. Spectra modelled using our diffuse geometry with Nc = 1024 H m-2 ; we plot the sum of stellar and dust emission. We use the grain size distributions preferred by Weingartner & Draine (2001) for Milky Way sight lines with RV = 3.1, RV = 4.0 and RV = 5.5. The values of bC in these models are, respectively, 6 â 10-5 , 4 â 10-5 and 3 â 10-5 .

= L e-C

,abs

Nc

,

(33)

where L is the luminosity of the central heating source. When Nc = 1024 H m-2 , stellar emission contributes 95 per cent of emission at 4 µm. Convolution of this spectrum with the passbands of Spitzer's InfraRed Array Camera (IRAC) reveals that


Infrared Spectra of Galaxies
10
24

11

1025 (a) Nc = 1024 H m r0 = 1 kpc r1 = 3 kpc
-2

(a) Early-type system 1024 Age = 0.89 Gyr (z = 6) Age = 2.5 Gyr (z = 2.5) Age = 5.6 Gyr (z = 1) 1023 Age = 13 Gyr (z = 0)

10

23

10

22

1022 10
-1 21

Stellar emission Combined emission

1021

L /W Hz

10 10

20 24

1024 (b) Nc = 1026 H m r0 = 1 kpc r1 = 3 kpc
-2

(b) Late-type system 1023 Age = 0.89 Gyr (z = 6) Age = 2.5 Gyr (z = 2.5) Age = 5.6 Gyr (z = 1)
-1

10

23

Age = 13 Gyr (z = 0) 1022

10

22

L /W Hz Stellar emission Combined emission 0.1 1 /µm 10 100

1021 10
21

1020 10
20

1024 (c) Late-type system; Nc = 1025 H m 1023 Age = 0.89 Gyr (z = 6) Age = 2.5 Gyr (z = 2.5) Age = 5.6 Gyr (z = 1) Age = 13 Gyr (z = 0) 10
22 -3

Figure 11. Spectra of dusty star-forming regions modelled using our circumnuclear geometry. For the heating radiation field, a stellar population with continuous SFH, star-formation rate = 1 M yr-1 and age 1 Gyr is used. Panels (a) and (b) show spectra of dust shells of column densities Nc = 1024 and 1026 m-2 respectively.

1021

stellar emission contributes 92, 91, 34 and 11 per cent in the channels at 3.5, 4.5 5.8 and 8.0 µm respectively. At Nc = 1026 H m-2 , these contributions reduce to 60 per cent of emission at 4 µm, and 48, 45, 3.9 and 0.88 per cent in the IRAC channels at 3.5, 4.5, 5.8 and 8.0 µm respectively. An absorption feature in the attenuated stellar emission component is seen around 10 µm, indicating that the silicate 9.7 µm feature has an optical depth of 0.4 in this model. In these model spectra, two PAH features appear well-isolated from their neighbours: that at 3.3 µm and the new 1.05 µm feature added by DL07. In the Nc = 1026 H m-2 model, these have equivalent widths 0.3 and 0.04 µm respectively. In the Nc = 1024 H m-2 model, the 3.3 µm feature has equivalent width 0.07 µm; the 1.05 µm feature is not apparent. In these plots we have, as previously, assumed the grain size distribution which Weingartner & Draine (2001) fit to Milky Way sight lines with RV = 3.1, assuming bC = 6 â 10-5 . From Figure 10, we can see that the effect of using a size distribution appropriate for a larger value of RV would be a small reduction in dust emission at all wavelengths. In Figure 12 we trace the evolution of our model spectra as a function of age for a range of SFHs. Figure 12(a) traces the evolution of the early-type galaxy modelled in Figure 5; Figure 12(b) traces that of the system with ongoing star formation modelled in Figure 4. Figure 12(c) shows a system with the same SFH as that modelled in Figure 12(b), but assuming a constant dust shell column density Nc = 1025 H m-2 . The time-invariance of the emission at 300 nm can be understood in terms of an equilibrium being reached between star formation and stellar death among the mas-

1020

0.1

1 /µm

10

100

Figure 12. The time-evolution of the spectra of dust-enshrouded starforming regions modelled using our circumnuclear geometry. In all cases, r0 = 1 kpc and r1 = 3 kpc. Panel (a) shows the evolution of the SED of the early-type galaxy modelled in Figure 5; the star-formation history and dust column density assumed are as shown in that figure for the perfect mixing model. Panel (b) shows the evolution of the system with ongoing star formation modelled in Figure 4. Panel (c) shows a system with the same SFH as Panel (b), but assuming a constant column density Nc = 1025 H m-2 .

sive stars which produce it. This contrasts with the behaviour seen in Figure 12(b), where the system becomes more metal-rich and dusty as it evolves, increasing extinction at ultraviolet wavelengths.

10 THE DEPENDENCE OF MIR LUMINOSITY UPON S FR To probe the dependence of the MIR luminosities of our models upon SFR, we constructed a set of circumnuclear models, each having a continuous SFH of age 1 Gyr but different star formation rates, ranging from 0.02 to 100 M yr-1 . We adopt otherwise constant model parameters: r0 = 1 kpc, r1 = 3 kpc and Nc = 1025 H m-2 . We constructed two sets of models, using the LD01 and the DL07 PAH cross sections, to compare their predictions. We compare these to the sources in the Spitzer extragalac-


12
tic First Look Survey (FLS) field. We correlated this source catalogue with the Fourth Data Release of the Sloan Digital Sky Survey (SDSS; York et al. 2000; Adelman-McCarthy et al. 2006) using a 2 matching radius to obtain redshifts for these sources, and used emission line measurements from the MPA/JHU catalogue6 . To remove AGN and spurious/marginal detections from the catalogue, and to correct for aperture effects and extinction, we followed Nikolic et al. (2004): AGN were identified using the spectral classification scheme of Veilleux & Osterbrock (1987); aperture and extinction corrections were applied to the emission line data using the prescription of Hopkins et al. (2003). We estimated the SFRs of the sources by comparing their H luminosities with that predicted by S TA R B U R S T 9 9 for a 1 Gyr-old continuous SFHs. K -corrections for the Spitzer 8 µm and 24 µm fluxes were derived using our models described above, selecting for each object that model which most closely matched the H -derived SFR. We convolved our model spectra with the 8-µm and 24-µm passbands of Spitzer to produce matching data from our model and the observed source population. The results are shown in Figure 13; also shown are power-law fits to the observed source population, the fit coefficients for which are listed in Table 2 (c.f. a similar fit performed by Wu et al. 2005). In the upper panels, we use the LD01 PAH cross sections (see Section 5) in our models; in the lower panels we use the DL07 cross sections. The non-linearity of our model L (24 µm) - relations at > 10 M yr-1 results from the increasing emission of grains excited by multiple photons at the highest star formation rates. These results are insensitive to the chosen geometry, as the emission from our models at both 8 and 24 µm is dominated by grains in the transiently heated regime. The luminosities are, however, dependent upon the column density of dust; this is likely to be one cause of the scatter of the FLS sources from our model relations. Increasing Nc by an order of magnitude, to 1026 H m-2 , increases our predictions of L (8 µm) by a factor of around 1.75, and L (24 µm) by a factor of around 2. Increasing Nc further causes our predictions for L (8 µm) to decrease as self-absorption becomes significant. Further scatter in these plots is likely to result from the uncertainly in the conversion of H luminosities into star formation rates. It is interesting to note that a better fit to the 24-µm luminosities of the FLS sources is obtained using the LD01 PAH cross sections compared to the DL07 cross sections (see Figures 13b and 13d). This difference can be explained by the enhanced PAH cross sections at 24 µm in the LD01 model, as seen in Figures 3(a) and 3(b). The fit to the 24 µm luminosities in Figure 13(d) can be improved by increasing Nc , but this leads to an over-prediction of the 8-µm luminosities in Figure 13(c).

8 µm 24 µm a 9.45 ± 0.05 9.06 ± 0.04 b 0.87 ± 0.08 0.97 ± 0.07 N 75 81

Table 2. Coefficients of the power-law fits shown in Figure 13, where log10 ( L ( )/L ) = a + b â log10 / M yr-1 . The number of sources used in each fit is N . Although the K -corrections derived for the FLS galaxies differ depending upon whether the LD01 or DL07 PAH cross sections are used, these differences are modest; the changes to the power-law coefficients above are much smaller than the quoted errors.
4 (IRAC4.5 - IRAC8.0) / magnitudes Late-type systems Early-type systems Nc = 1026 H m Nc = 10
25 -2

3

Hm

-2

2 Nc = 1026 H m 1
-2

Nc = 1024 H m

-2

0

Nc = 1025 H m

-2

-1

Nc = 1024 H m -1 -0.5

-2

0

0.5

1

1.5

2

(IRAC3.6 - IRAC5.8) / magnitudes Nc = 1026 H m
-2

4

(J - IRAC8.0) / magnitudes

Nc = 1025 H m 2

-2

Nc = 1024 H m 0

-2

Nc = 1026 H m

-2

Nc = 1025 H m -2 Nc = 1024 H m -4
-2

-2

-2

0

2

(J - MIPS24) / magnitudes

Figure 14. Plot of the positions of galaxies modelled by our circumnuclear geometry in colour-colour space. Each trace represents the evolution of a model galaxy from an age of 1.6 Gyr (top-right end) to 13.8 Gyr (lowerleft end). Late-type galaxies, with continuous SFHs, are shown with solid traces; early-type galaxies, with instantaneous SFHs, are shown with dashed traces. Three traces are shown for each galaxy type, representing galaxies enshrouded with three different column densities of dust. All colours are shown in magnitudes, and are calculated in the galaxy rest frame.

11 PREDICTED GALAXY COLOURS Using the models presented in Section 9, we can make predictions of the relative merits of colour diagnostics in the discrimination of systems with ongoing star formation from those with only old stellar populations. In Figure 14, we show the colours of lateand early-type models, with continuous and instantaneous SFHs respectively, enshrouded by three different column densities of dust. In the lower panel we consider colours formed between the J -band of the UK InfraRed Telescope (UKIRT), as used by the UKIRT

6

Available from: http://www.mpa- garching.mpg.de/SDSS/

Infrared Deep Sky Survey (UKIDSS), in conjunction with Spitzer MIR photometry. In the upper panel, we consider colours which use only IRAC photometry. In all cases, the late-type models appear redder than systems with only old stellar populations, a result of the enhanced UV-pumping of PAH grains. Models which are enshrouded by higher column densities of dust appear systematically redder, as a result of their enhanced dust emission. All of these colour predictions are robust to changes in starformation rate . In our model, both stellar and PAH emission are linearly proportional to except at the highest star formation rates (see Figure 13); only the far-infrared greybody varies in colour with . These predictions are also robust to changes in the assumed


Infrared Spectra of Galaxies
(a) LD01 PAH model (b) LD01 PAH model

13

10 L (8 µm)/L

10

10

10

L (24 µm)/L

109

10

9

108

FLS galaxies Circumnuclear models Power-law fit

10

8

107 (c) DL07 PAH model (d) DL07 PAH model

10

7

10 L (8 µm)/L

10

10

10

L (24 µm)/L

109

10

9

108

10

8

107 0.1 1 /M yr-1 10 100 0.1 1 /M yr-1 10 100

10

7

Figure 13. Solid lines: The Spitzer -convolved 8-µm (Left panels) and 24-µm (Right panels) luminosities of circumnuclear models with continuous SFHs, age 1 Gyr, r0 = 1 kpc, r1 = 3 kpc and Nc = 1025 H m-2 . Dashed lines: Power-law fits to the Spitzer extragalactic FLS source population (points). In the upper panels we use the LD01 PAH cross sections; in the lower panels we use the DL07 cross sections.

masses of the galaxies if is assumed to scale linearly with mass whilst Nc remains constant. In Figure 15 we investigate the redshift and column density dependence of the colours shown in Figure 14 by tracing the redshift evolution of systems of constant Nc . The diagnostic power of colours formed between IRAC channels is limited beyond z 0.5; at higher redshifts these bands are increasingly dominated by stellar emission, especially at lower dust column densities. The power of J - 8 µm and J - 24 µm colours extends to higher redshifts, though a degeneracy exists between highly obscured old stellar populations and recent star formation. In all cases, the two SFHs considered in Figure 15 converge at high redshift; both contain relatively young stellar populations shortly after their formation at z = 10.

sections of Li & Draine (2001) than those which use the cross sections of Draine & Li (2007). In Section 11, we provide estimates of the near-to-midinfrared colours of early- and late-type galaxies, demonstrating the power of our model as a tool for interpreting colour-colour diagrams.

ACKNOWLEDGMENTS DCF acknowledges the receipt of a PPARC studentship. This work has made use of the distributed computation grid of the University of Cambridge (C A M G R I D). We are grateful to the anonymous referee for helpful comments.

12 SUMMARY We have developed a semi-empirical model of the emission of dust in star-forming galaxies, capable of producing model spectra which extend from the Lyman limit through three orders of magnitude in wavelength to the far-infrared. A particular strength is its state-ofthe-art model of the mid-infrared PAH features: we treat the transient heating of small dust grains using the sophisticated method of Draine & Li (2001) and use the grain interaction cross sections of Draine & Li (2007). In Section 8, we showed how the predictions of our model depend upon the microscopic properties of the dust grain population ­ its grain size spectrum and ionisation state, for example. In Sections 9 and 10, we went on to demonstrate the spectra produced by our model depend upon the large scale properties of galaxies: their star formation histories and dust column densities. In Figure 13, we found that the Spitzer extragalactic First Look Survey (FLS) source population is better matched by models which use the PAH cross

APPENDIX A: CHOICE OF ENERGY BINS The energy bins used in Section 6.1 are chosen so that P(E ) is well-sampled in regions where it is changing rapidly. For the smallest grains, this means that the ground state, and those close to it, must be well-sampled. To ensure that this is the case, we set the N1 lowest-energy bins in our scheme to contain the N1 lowest-energy vibrational states of dust grains, one in each bin, thus treating them exactly. At the opposite extreme, we require the occupation probability of our highest-energy bin to be very small, such that the high-energy tail of P(E ) is not truncated. To meet these requirements, we place the next N2 bins at linear intervals of U N1 - U N1 -1 above U N1 . We then place the remaining


14
6 (a1) Continuous SFH (J - IRAC8.0) / magnitudes 4 (a2) Instantaneous SFH

2

0

-2 6 (J - MIPS24) / magnitudes (b1) Continuous SFH (b2) Instantaneous SFH

4

2

0

-2 -4

(IRAC3.6 - IRAC5.8) / magnitudes

0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 3

(c1) Continuous SFH

(c2) Instantaneous SFH

(IRAC4.5 - IRAC8.0) / magnitudes

(d1) Continuous SFH

Nc = 1024 H m Nc = 1025 H m Nc = 1026 H m

-2 -2 -2

(d2) Instantaneous SFH

2

1

0

-1

0

0.5

1 z

1.5

2

0

0.5

1 z

1.5

2

Figure 15. The colours of galaxies modelled by our circumnuclear geometry with continuous (left panels) and instantaneous (right panels) SFHs; in each case, star formation starts at z = 10. Each trace represents the redshift evolution of a system of constant Nc . All colours are shown in magnitudes, and are calculated in the frame of the observer.

Nbins - N1 - N2 bins at logarithmic intervals above U multiplicative spacing given by: 1.1 1.5 a a a a

N1 +N2

wi t h

In this paper, we use N1 = 20, N2 = 50 and Nbins = 600.

=

1.5 + 7.0

° a-50 A ° ° 350 A-50 A

(7.0 - 1.5)

° 10 A ° 50 A ° 350 A

< < <

° 10 A ° 50 A ° 350 A (A1)

APPENDIX B: CALCULATING THE TRANSITION MATRIX T Several approximations for calculating the emission from stochastically heated dust grains are discussed in detail by Draine & Li (2001). They all rely on dividing the possible enthalpy content of


Infrared Spectra of Galaxies
the grain into of the order of 500 bins and calculating the probability of a grain transitioning from an enthalpy that belongs to one bin to an enthalpy belonging to a different bin. The upward transition rates T ji are given by (Draine & Li 2001, equation 15): T ji = cU j U j -U
W4 i W1

15
U
max j

U

j

U

m in j

G ji (E )Cabs (E )uE (E ) dE
Wc

for j < Nbins , (B1)
W1 W2 W3 W4 Uim Ui
ax

T

ji

=

c U j -U +
Wc

i

W1

E - W1 Cabs (E )uE (E ) dE Wc - W1 for j = Nbins (B2)

Cabs (E )uE (E )dE

Uim Figure B1. Illustration of the energy boundaries used in function G ji (E ), where j and i are the indices of the upper and lower bins respectively and U min and U max are the boundaries of the ith bin. Note that W2 is the smaller i i and W3 the larger of the pair of values: (U max - U max ) and (U min - U min ). j i j i

in

where Cabs is the grain density of the radiation finite width of bins and mi G ji =

absorption cross section, uE is the energy field and G ji is a correction factor for the is given by: E - W1 U j U n U j , U i U j U i W4 - E U j U
i

W1 < E < W2 , W2 < E < W3 , W3 < E < W4 , otherwise.

(B3)

RE F E RE NCE S Adelman-McCarthy J. K., Agueros M. A., Allam S. S., Anderson ¨ K. S. J., Anderson S. F., Annis J., Bahcall N. A., Baldry I. K., Barentine J. C. a., 2006, ApJS, 162, 38 Blain A. W., Smail I., Ivison R. J., Kneib J.-P., 1999, MNRAS, 302, 632 Blain A. W., Smail I., Ivison R. J., Kneib J.-P., Frayer D. T., 2002, Phys. Rep., 369, 111 Bohren C. F., Huffman D. R., 1998, Absorption and Scattering of Light by Small Particles. Wiley Science Paperback Series Chandresekhar S., 1960, Radiative Transfer. Dover, New York Draine B. T., 2003, ARA&A, 41, 241 Draine B. T., Lee H. M., 1984, ApJ, 285, 89 Draine B. T., Li A., 2001, ApJ, 551, 807 Draine B. T., Li A., 2007, ApJ, 657, 810 Draine B. T., Malhotra S., 1993, ApJ, 414, 632 Efstathiou A., Rowan-Robinson M., 1995, MNRAS, 273, 649 Efstathiou A., Rowan-Robinson M., Siebenmorgen R., 2000, MNRAS, 313, 734 Fagotto F., Bressan A., Bertelli G., Chiosi C., 1994, A&AS, 105, 29 Ferland G. J., 1980, PASP, 92, 596 Fritz J., Franceschini A., Hatziminaoglou E., 2006, MNRAS, 366, 767 Guhathakurta P., Draine B. T., 1989, ApJ, 345, 230 Hauser M. G., Dwek E., 2001, ARA&A, 39, 249 Hopkins A. M., Miller C. J., Nichol R. C., Connolly A. J., Bernardi M., Gomez P. L., Goto T., Tremonti C. A., Brinkmann ´ J., Ivezic Z., Lamb D. Q., 2003, ApJ, 599, 971 ´ Hughes D. H., Serjeant S., Dunlop J., Rowan-Robinson M., Blain A., Mann R. G., Ivison R., Peacock J., Efstathiou A., Gear W., Oliver S., Lawrence A., Longair M., Goldschmidt P., Jenness T., 1998, Nature, 394, 241 Kennicutt Jr. R. C., Armus L., Bendo G., Calzetti D., Dale D. A., Draine B. T., Engelbracht C. W., Gordon K. D., Grauer A. D., Helou G., Hollenbach D. J., Walter F., 2003, PASP, 115, 928 Leitherer C., Schaerer D., Goldader J. D., Delgado R. M. G., Robert C., Kune D. F., de Mello D. F., Devost D., Heckman T. M., 1999, ApJS, 123, 3 Li A., Draine B. T., 2001, ApJ, 554, 778 Lonsdale C., Polletta M. d. C., Surace J., Shupe D., Fang F., Xu

i

0

The integration limit quantities W1 , W2 , W3 and W4 are defined in m Figure B1. Wc is equal to UNbiins - Uimin . n Throughout this paper, we use the continuous cooling approximation for the downward transitions (Draine & Li 2001, equation 41): T
ji

=

1 8 U i - U j h3 c2 â E
0 u

E 3Cabs (E ) dE exp(E /k u ) - 1

for i > 1; j = i - 1 (B4) (B5)

T ji = 0

where u is the characteristic temperature of bin u, as defined by Draine & Li (2001). The computational cost of evaluating the transition matrix T ji is dominated by evaluating the upward transitions, in particular the integrand in Equation (B1), separately for each {i, j} pair. It may be reduced by noting that, although G ji is a function of i and j, it always remains a linear function of E. Hence the integrand in Equation (B1) is, for any {i, j} pair, a linear combination of:
b

for all other {i, j}

I1 (a, b) = and
b

a

d ECabs (E )uE (E )

(B6)

I2 (a, b) =

a

d E · E · Cabs (E )uE (E ).

(B7)

where I1 and I2 are independent of i and j. We can then write T ji as a linear combination of I1 (W1 , W2 ), I1 (W2 , W3 ), I1 (W3 , W4 ), I2 (W1 , W2 ), I2 (W2 , W3 ) and I2 (W3 , W4 ), with coefficients which are independent of E . We evaluate I1 and I2 by non-adaptive integration using a pre-computed grid of the two integrands appearing in I1 and I2 . By using this grid, we reduce the number of necessary evaluations of Cabs (E ) and uE (E ) from O (N 2 ) to O (103 ). Non-adaptive integration is appropriate since uE , the output of stellar population models, is computed on a fixed grid.


16
C. K., Smith H. E., Siana B., Rowan-Robinson M., Babbedge T., Puetter R., 2004, ApJS, 154, 54 Lonsdale C. J., Smith H. E., Rowan-Robinson M., Surace J., Shupe D., Xu C., Oliver S., Padgett D., Fang F., Conrow T., Franceschini A., Condon J. J., Dole H., Serjeant S., 2003, PASP, 115, 897 Mathis J. S., Mezger P. G., Panagia N., 1983, A&A, 128, 212 Mattioda A. L., Allamandola L. J., Hudgins D. M., 2005a, ApJ, 629, 1183 Mattioda A. L., Hudgins D. M., Allamandola L. J., 2005b, ApJ, 629, 1188 Mezger P. G., Mathis J. S., Panagia N., 1982, A&A, 105, 372 Nikolic B., Cullen H., Alexander P., 2004, MNRAS, 355, 874 Piovan L., Tantalo R., Chiosi C., 2006a, MNRAS, 366, 923 Piovan L., Tantalo R., Chiosi C., 2006b, MNRAS, 370, 1454 Rowan-Robinson M., 1980, ApJS, 44, 403 Rowan-Robinson M., Babbedge T., Surace J., Shupe D., Fang F., Lonsdale C., Smith G., Polletta M., Siana B., Gonzalez-Solares E., Stacey G., Vaccari M., 2005, AJ, 129, 1183 Rowan-Robinson M., Crawford J., 1989, MNRAS, 238, 523 Rowan-Robinson M., Lari C., Perez-Fournon I., Gonzalez-Solares E. A., La Franca F., Vaccari M., Oliver S., Gruppioni C., Ciliegi P., Heraudeau P., Carraminana A., Mujica R., 2004, MNRAS, ´ ~ 351, 1290 Silva L., Granato G. L., Bressan A., Danese L., 1998, ApJ, 509, 103 Smail I., Ivison R. J., Blain A. W., 1997, ApJ, 490, L5+ Smail I., Ivison R. J., Blain A. W., Kneib J.-P., 1998, ApJ, 507, L21 Smith J. D. T., Draine B. T., Dale D. A., Moustakas J., Kennicutt Jr. R. C., Helou G., Armus L., Roussel H., Sheth K., Bendo G. J., Murphy E. J., Walter F., 2007, ApJ, 656, 770 Swinbank A. M., Smail I., Chapman S. C., Blain A. W., Ivison R. J., Keel W. C., 2004, ApJ, 617, 64 Takagi T., Vansevicius V., Arimoto N., 2003, PASJ, 55, 385 Takeuchi T. T., Ishii T. T., Dole H., Dennefeld M., Lagache G., Puget J.-L., 2006, A&A, 448, 525 Vazquez G. A., Leitherer C., 2005, ApJ, 621, 695 ´ Veilleux S., Osterbrock D. E., 1987, ApJS, 63, 295 Weingartner J. C., Draine B. T., 2001, ApJ, 548, 296 Wu H., Cao C., Hao C.-N., Liu F.-S., Wang J.-L., Xia X.-Y., Deng Z.-G., Young C. K.-S., 2005, ApJ, 632, L79 Xu C., Lonsdale C. J., Shupe D. L., O'Linger J., Masci F., 2001, ApJ, 562, 179 York D. G., Adelman J., Anderson Jr. J. E., Anderson S. F., Annis J., Bahcall N. A., Bakken J. A., Barkhouser R., Bastian S., Berman E., Boroski W. N., Yanny B., Yasuda N., 2000, AJ, 120, 1579