Modeling the Dust Properties of z ~ 6 Quasars with ART^2 -- All-wavelength Radiative Transfer with Adaptive Refinement Tree

1.2 MB·English

by  Yuexing Li

DRAFTVERSIONFEBRUARY5,2008
PreprinttypesetusingLATEXstyleemulateapjv.9/08/03
MODELINGTHEDUSTPROPERTIESOFZ 6QUASARSWITHART2—ALL-WAVELENGTHRADIATIVE
∼
TRANSFERWITHADAPTIVEREFINEMENTTREE
YUEXINGLI1,PHILIPF.HOPKINS1,LARSHERNQUIST1,DOUGLASP.FINKBEINER1,THOMASJ.COX1,
VOLKERSPRINGEL2,LINHUAJIANG3,XIAOHUIFAN3,NAOKIYOSHIDA4
1Harvard-SmithsonianCenterforAstrophysics,HarvardUniversity,60GardenStreet,Cambridge,MA02138,USA
2Max-Planck-InstituteforAstrophysics,Karl-Schwarzschild-Str. 1,85740Garching,Germany
3DepartmentofAstronomy,UniversityofArizona,933N.CherryAve,Tucson,AZ85721and
4NagoyaUniversity,Dept.ofPhysics,Nagoya,Aichi464-8602,Japan
DraftversionFebruary5,2008
8 ABSTRACT
0
Thedetectionoflargequantitiesofdustinz 6quasarsbyinfraredandradiosurveyspresentspuzzlesforthe
0 ∼
formationandevolutionofdustintheseearlysystems. Previously(Lietal.2007), weshowedthatluminous
2
quasarsatz&6canformthroughhierarchicalmergersofgas-richgalaxies,andthatthesesystemsareexpected
n toevolvefromstarburstthroughquasarphases.Here,wecalculatethedustpropertiesofsimulatedquasarsand
a theirprogenitorsusingathree-dimensionalMonteCarloradiativetransfercode,ART2–All-wavelengthRadia-
J
tiveTransferwithAdaptiveRefinementTree.ART2incorporatestheradiativeequilibriumalgorithmdeveloped
9
byBjorkman&Wood(2001)whichtreatsdustemissionself-consistently,anadaptivegridmethodwhichcan

efficientlycoveralargedynamicrangeinspatialscalesandcancaptureinhomogeneousdensitydistributions,
]
h amultiphasemodeloftheinterstellarmediumwhichaccountsfortheobservedscalingrelationsofmolecular
p clouds,andasupernova-originmodelfordustwhichcanexplaintheexistenceofdustincosmologicallyyoung,
- starburstingquasars.ByapplyingART2tothehydrodynamicsimulationsofLietal.(2007),wereproducethe
o
observedspectralenergydistribution(SED)andinferreddustpropertiesofSDSSJ1148+5251,themostdistant
r
t quasar detected in the Sloan survey. We find that the dust and infraredemission are closely associated with
s the formationand evolutionofthe quasarhost. As the system evolvesfroma starburstto a quasar, the SED
a
[ changes from being dominated by a cold dust bump (peaking at 50µm) to one that includes a prominent
hotdustcomponent(peakingat 3µm),andthegalaxyevolvesfr∼omacoldtoawarmultraluminousinfrared

3 galaxy (ULIRG) owing to heati∼ng and feedback from star formation and the active galactic nucleus(AGN).
v Furthermore,theAGNactivityhassignificantimplicationsfortheinterpretationofobservableaspectsofthe
6 hosts. Thehottestdust(T &103K)ismostnoticeableonlyduringthepeakquasaractivity,andcorrelateswith
0
thenear-IRflux.However,wefindnocorrelationbetweenthestarformationrateandfar-IRluminosityduring
7
this phase owingto strongAGN contamination. Our results suggestthat vigorousstar formationin merging
3
progenitorsisnecessarytoreproducetheobserveddustpropertiesofz 6quasars,supportingamerger-driven
6. originforluminousquasarsathighredshiftsandthestarburst-to-quasa∼revolutionaryhypothesis.
0
Subject headings: quasar: formation — quasar: evolution — quasar: high redshift — galaxies: starbursts
7
— infrared: galaxies — radiative transfer — interstellar medium — dust, extinction —
0
individual:SDSSJ1148+5251
:
v
i
X 1. INTRODUCTION ture” quasars and the formation of their hosts in rapid star-
r High-redshift quasars are important for understanding the burstsatevenhigherredshifts(z&10).
a Follow-up, multi-wavelength observations have been car-
formation and evolution of galaxies and supermassive black
riedoutforthesez 6quasars,fromX-ray(e.g.,Brandtetal.
holes(SMBHs) in the earlyUniverse. In the pastfewyears,
∼
2002; Stratevaetal. 2005; Vignalietal. 2005; Steffenetal.
nearly two dozen luminous quasars have been discovered
2006; Shemmeretal. 2005, 2006), to optical /infrared (e.g.,
by the Sloan Digital Sky Survey (SDSS; Yorketal. 2000)
Barthetal.2003;Pentericcietal.2003;Freudlingetal.2003;
and the Canada-France High-z Quasar Survey (CFHQS;
Whiteetal.2005;Willottetal. 2005;Hinesetal.2006), and
Willottetal. 2007) at z 6, corresponding to a time when
∼ radio wavelengths (e.g., Carillietal. 2001; Bertoldietal.
the Universe was less than a billion years old (Fanetal.
2003b; Walteretal. 2003; Carillietal. 2004; Wangetal.
2003, 2004, 2006a,b). As summarized by Fan (2006),
these quasars are rare ( 10- 9Mpc- 3 comoving); believed 2007). A noteworthy result of these studies is their detec-
to be powered by SMB∼Hs with masses 109M (e.g., tion of dust in these high-redshift objects. In particular,
∼ ⊙ deep infrared and radio surveys (e.g., Robsonetal. 2004;
Willottetal. 2003; Barthetal. 2003); reside near the end
Bertoldietal. 2003a; Carillietal. 2004; Charmandarisetal.
of the epoch of reionization, as indicated by Gunn-Peterson
2004; Beelenetal. 2006; Hinesetal. 2006) have revealed a
absorption troughs(Gunn&Peterson 1965) in their spectra;
large amount of cold dust in SDSS J1148+5251 (hereafter
and have similar spectral energy distributions (SEDs) and
J1148+5251),the most distant Sloan quasar detected at red-
comparable metallicity to lower-redshift counterparts (e.g.,
shiftz=6.42(Fanetal.2003). Thedustmassisestimatedto
Elvisetal. 1994; Glikmanetal. 2006; Richardsetal. 2006;
be 1- 7 108M . The detection of dust is also reported
Hopkinsetal. 2007d), implying the early presence of “ma- ∼ × ⊙
in the first four CFHQS quasars at z>6, including the new
record holder CFHQS J2329-0301 at z=6.43 (Willottetal.
Electronicaddress:[email protected]
2007). Maiolinoetal. (2004) argue, from the observeddust
2 Lietal.
extinctioncurveofSDSSJ1048+4637atz=6.2,thatthedust Robertsonetal. 2006). In this picture, quasars are, there-
inthesehigh-zsystemsmaybeproducedbysupernovae. fore, descendents of starburst galaxies (Sandersetal. 1988;
Jiangetal. (2006) have performeda comprehensivestudy Norman&Scoville1988;Scoville2003).
ofthirteenz 6quasarsbycombiningnewSpitzerobserva- To make detailed contact with observations of z 6
∼ ∼
tionswithexistingmulti-banddata. Itappearsthatnearlyall quasars, it is necessary to theoretically predict the SEDs of
13ofthesequasarsexhibitprominentinfraredbumpsaround these systems and how they evolve with time. Solving ra-
both λ 3µm and λ 50µm. In a dusty galaxy, most of diation transport in dusty, starbursting quasars is, however,
∼ ∼
the radiation from newly formed stars or an AGN is ab- difficultowingtothenon-localityofthesources,theopacity
sorbed by dust and re-emitted at infrared wavelengths. The andmultiphase characterof the interstellar medium, andthe
SEDs of star-forminggalaxiestypicallypeak at 50 – 80 µm complexmorphologyofthesesystems. Overthelastseveral
(rest frame) and can be approximated by a modified black- decades, a variety of numericaltechniqueshave been devel-
body spectrum with dust temperatures below 100 K (e.g., opedtosolvetheradiativetransferproblemwithdifferentlev-
Dunne&Eales 2001), while in quasars, some dust can be elsof approximationandin multi-dimensions. Based onthe
heated directly by the AGN to temperatures up to 1200 K specific algorithms employed, the approaches can be classi-
(e.g.,Glikmanetal.2006)anddominatethe near-tomid-IR fied into two general categories: finite-differenceand Monte
emission. Therefore, the observations by Jiangetal. (2006) Carlo methods (e.g., Jonsson 2006; see also Pascuccietal.
indicate the presence of large amounts of both hot and cold 2004 who describe the codes as “grid-based” or “particle-
dustinthesehigh-zsystems. based”byanalogytohydrodynamicsolvers).
However,thenatureofthedust,anditsformationandevo- Finite-difference codes solve the equations of radiative
lution in the contextof the quasar host, are unclear. For ex- transfer (RT) iteratively using finite convergence criteria.
ample,thedustdistributionisunknown.Ithasbeensuggested They either solve the moment equations for RT, originally
thatthehotdustliesinthecentralregionsandisheatedbyan formulated by Hummer&Rybicki (1971) for spherical
AGN to producenear-IRemission (Rieke&Lebofsky1981; geometrywith a centralpointsource (e.g., Scoville&Kwan
Pollettaetal. 2000; Haasetal. 2003), while the warm and 1976; Leung 1976; Yorke 1980; Wolfire&Cassinelli 1986;
colddustcanextendtoafewkpcanddominateatmid-IRand Men’shchikov&Henning 1997; Dullemond&Turolla
far-IRwavelengths(Pollettaetal.2000;Nenkovaetal.2002; 2000), or employ ray-tracing methods on a discrete
Siebenmorgenetal.2005). However,therelativeimportance spatial grid for complex density configurations (e.g.,
ofheatingbystarsandAGNactivityisuncertain. Thisises- Rowan-Robinson 1980; Efstathiou&Rowan-Robinson
sentialtoanaccuratedeterminationofthestarformationrate 1990, 1991; Steinackeretal. 2003; Folinietal. 2003;
(SFR), but cannot be inferred simply from observed SEDs. Steinackeretal. 2006). This technique provides full error
Currently, the most common method used to derive a star controlbutcanbetime-consuming.
formation rate is to assume that most of the FIR luminosity Monte Carlo methods sample and propagate photons
comes from young stars. However, if the FIR luminosity is probabilistically (e.g., Witt 1977; Lefevreetal. 1982,
mainlycontributedbyanAGN, thentheSFR wouldbesub- 1983; Whitney&Hartmann 1992; Wittetal. 1992;
stantiallyreduced. Finally,itisnotknownhowtheSEDand Code&Whitney 1995; Lopezetal. 1995; Lucy 1999;
dustcontentofaquasaranditshostevolvewithtime. Wolfetal. 1999; Bianchietal. 2000; Harries 2000;
In order to address these questions, we must combine a Bjorkman&Wood 2001; Whitneyetal. 2003b; Jonsson
modelfortheformationofaquasarwithradiativetransfercal- 2006;Pinteetal.2006). TheMonteCarlotechniqueismore
culations that treat dust emission self-consistently, to follow flexible than the finite-difference one because it tracks the
thephysicalproperties,environment,andevolutionofquasar scattering,absorptionandre-emissionofphotons(orphoton
hostsandtheirdustcontent. packets)indetail,andcanhandlearbitrarygeometries,butat
Earlier (Lietal. 2007), we developed a quasar formation the cost of computational expense to reduce Poisson noise.
modelwhichself-consistentlyaccountsforblackholegrowth, However,advancesin computingtechnologyand algorithms
starformation,quasaractivity,andhostspheroidformationin have made high-accuracy Monte Carlo RT calculations
thecontextofhierarchicalstructureformation. Weemployed feasibleandpopular.
asetofmulti-scalesimulationsthatincludedlarge-scalecos- In particular, Bjorkman&Wood (2001) have developed a
mologicalsimulationsandgalaxymergersongalacticscales, Monte Carlo code to handle radiative equilibrium and tem-
together with a self-regulated model for black hole growth, perature corrections, which calculates dust emission self-
to produce a luminous quasar at z 6.5, which has a black consistently. It conserves the total photon energy, corrects
holemass( 2 109M )andanum∼berofpropertiessimilar thefrequencydistributionofre-emittedphotons,andrequires
⊙
∼ ×
to J1148+5251(Fanetal. 2003). In our scenario, luminous, no iteration as the dust opacity is assumed to be indepen-
high-redshift quasars form in massive halos that originate dent of temperature. This code has been used in a num-
from rare density peaks in the standard ΛCDM cosmology, ber of applications, including protostars that have a disk
andtheygrowthroughhierarchicalmergersofgas-richgalax- and an envelope with a single heating source in the cen-
ies. Gravitational torques excited in these mergers trigger ter (e.g., Whitney&Hartmann 1992, 1993; Whitneyetal.
largeinflowsof gasand producestrongshocksthatresultin 2003a,b, 2004); circumstellar envelopes (e.g., Woodetal.
intensestarbursts(Hernquist1989;Barnes&Hernquist1991, 1996a,b, 1998); protoplanetary systems (e.g., Woodetal.
1992, 1996; Hernquist&Mihos 1995; Mihos&Hernquist 2002; Riceetal. 2003); and galaxies that include a bulge
1996) and fuel rapid black hole accretion (DiMatteoetal. and a disk with multiple heating sources from these two
2005; Springeletal. 2005a). Moreover, feedback from populations(e.g.,Wood&Jones1997;Wood&Loeb2000).
the accreting black hole disperses the obscuring material, This code has also been able to generate optical–far-IR
briefly yielding an optically visible quasar (Hopkinsetal. SEDs that reproduce those of a sample of 21 X-ray se-
2005c,b,a, 2006), and regulating the M - σ correlation be- lected AGNs(Kuraszkiewiczetal. 2003), and the versionof
BH
tween the SMBHs and host galaxies (DiMatteoetal. 2005; Whitneyetal. (2003b) has been applied to simulations of
ModelingDustinz 6Quasars 3
∼
galaxy mergers with black hole feedback to study the local inal code. More important, it captures the inhomogeneous
ULIRGs (Chakrabartietal. 2007a) and submillimetergalax- densitydistributioningalaxymergersandreproducestheob-
ies(Chakrabartietal.2007b). servedSED anddustpropertiesofJ1148+5251basedonthe
InordertoproduceaccurateSEDsofquasarsandtheirhosts simulationsofLietal.(2007). Therefore,ART2 canbeused
formedbygalaxymergers,asinLietal.(2007),theRTcode topredictmulti-wavelengthpropertiesofquasarsystemsand
mustsatisfythefollowingrequirements:(1)tobeabletohan- theirgalaxyprogenitors.
dle arbitrary geometries and distributed heating sources; (2) Thispaperisorganizedasfollows. In§2,wedescribeour
resolvethelargedynamicrangesin spatialscalesanddensi- computationalmethodsandmodels,includingthemulti-scale
ties in the merger simulations; (3) incorporate a multiphase simulations of Lietal. (2007) for z 6 quasar formation,
descriptionoftheinterstellarmedium(ISM);and(4)employ and our implementation of ART2 wh∼ich incorporates radia-
adustmodelappropriateforayoungstarburstsystem. tiveequilibriumasinBjorkman&Wood(2001),anadaptive
In many earlier applications, radial logarithmic or uni- grid, a multiphase ISM, and a supernova-origindust model.
formly spaced meshes were used to generate density grids. In§3,wepresentthemulti-wavelengthSEDfromopticalto
However,in merginggalaxies, the ISM is clumpyandirreg- submillimeterofasimulatedquasaratz 6.5. Thedustdis-
∼
ular owing to shocks and tidal features produced during the tributionandpropertiesofthequasararediscussedin§4and
interactions, and has multiple density centers, making loga- wedescribetheirevolutionin§5.Weconsidertherobustness
rithmic algorithms inefficient. In such a situation, an adap- ofourresultsin§6andsummarizein§7.
tivegridapproachappearsidealforresolvinglocalizedhigh-
density regions, while still covering the large volumes of 2. METHODOLOGY
mergingsystems(Wolfetal.1999;Kurosawa&Hillier2001; In order to capture the physical processes underlying the
Harriesetal.2004;Jonsson2006). formationandevolutionofquasarsintheearlyUniverse,and
In addition, dust models commonly used in previous tocomparetheirmulti-wavelengthpropertiestoobservations,
worksarebasedonobservedextinctioncurvesfortheMilky we combine quasar formation and radiative transfer calcu-
Way (e.g., Weingartner&Draine 2001; Calzettietal. 1994; lations. We first perform a set of novel multi-scale simula-
Kimetal. 1994; Mathisetal. 1977; Calzettietal. 2000), in tionsthatyieldaluminousquasaratz 6.5(Lietal.2007).
which the dust is assumed to be produced mainly by old, We thenapplythe 3-D,Monte Carlo ra∼diativetransfercode,
low-mass stars with ages > 1 Gyr (Mathis 1990; Whittet ART2, to the outputs of the hydrodynamic galaxy mergers
2003; Dwek 2005). However, a large amount of dust simulationsto calculate the SED of the system. The models
( 1- 7 108M⊙) is detected in the z 6.42 quasar andsimulationsofquasarformationaredescribedindetailin
∼ × ≃
host (Bertoldietal. 2003a) at a time when the Universe Lietal.(2007). Here,webrieflysummarizethemodelingof
was only 850 Myr old. It has been suggested by ob- quasarformation,andthespecificationsofART2.
∼
servations (e.g., Maiolinoetal. 2004, 2006; Moseleyetal.
1989;Dunneetal.2003;Morganetal.2003;Sugermanetal. 2.1. FormationModelofz 6Quasars
2006) and theoretical studies (e.g., Todini&Ferrara 2001; ∼
2.1.1. Multi-scaleSimulations
Nozawaetal. 2003; Schneideretal. 2004; Nozawaetal.
2007; Bianchi&Schneider 2007) that supernovae can pro- Thez 6quasarsarerare(spacedensity 1Gpc- 3comov-
∼ ∼
vide fast and efficient dust formation in the early Universe, ing), andappearto bepoweredby supermassiveblack holes
motivatingotherchoicesforthedustmodelintheRT calcu- of mass 109M⊙. Therefore, simulations of high-redshift
∼
lations. quasarformationmustconsideralargecosmologicalvolume
Finally, in a multiphase description of the ISM toaccommodatethelowabundanceofthispopulation,havea
(McKee&Ostriker 1977), as adopted in our simula- largedynamicrangetofollowthehierarchicalbuild-upofthe
tions (Springel&Hernquist 2003a,b; Hopkinsetal. 2006), quasar hosts, and include realistic prescriptions for star for-
the “hot-phase”(diffuse)and“cold-phase”(densemolecular mation, black hole growth, and associated feedback mecha-
andHIcore)componentsco-existunderpressureequilibrium nisms. Themulti-scalesimulationsinLietal.(2007)include
but have different mass fractions and volume filling factors bothN-bodycosmologicalcalculationsinavolumeof3Gpc3
(i.e., the hot-phase gas is 10% in mass but & 99% in toaccountforthelownumberdensityofquasarsatz 6,and
volume),sobothphasescon≤tributetothedustextinctionand hydrodynamicalsimulationsofindividualgalaxyme∼rgerson
shouldthereforebeincludedintheRTcalculations. galactic scales to resolve gas-dynamics, star formation, and
We have refined the Monte Carlo RT code developed blackholegrowth.
by Bjorkman&Wood (2001) and Whitneyetal. (2003b) First, we perform a coarse dark matter-only simulation in
by implementing: an adaptive grid algorithm similar to a volume of 3Gpc3. The largest halo at z=0, within which
that of Jonsson (2006) for the density field and the ar- early, luminousquasarsare thoughtto reside (Springeletal.
bitrarily distributed sources; a multiphase ISM model 2005b), is then selected for resimulation at higher resolu-
(Springel&Hernquist 2003a) which accounts for observed tion. The evolution of this halo and its environment is re-
scalingrelationsofmolecularcloudsforthedustdistribution; simulated using a multi-grid zoom-in technique (Gaoetal.
and a supernova-origindust model for the opacity using the 2005) that provides much higher mass and spatial resolu-
grain size distribution of Todini&Ferrara (2001). We refer tion for the halo of interest, while following the surround-
toournewcodeasART2 (All-wavelengthRadiativeTransfer ing large-scale structure at lower resolution. The merging
withAdaptiveRefinementTree). ART2 iscapableofproduc- history of the largest halo at z 6, which has then reached
ing SEDs and images in a wide range of wavelengths from a mass of 7.7 1012M⊙ thr∼ough seven major (mass ra-
X-ray to millimeter. In the present paper we focus on the tio < 5:1)∼merge×rs between redshifts 14.4 and 6.5, is ex-
dust properties from optical to submillimeter bands. As we tracted. These major mergers are again re-simulated hydro-
show in what follows, ART2 reproduces the spectrum of a dynamicallyusinggalaxymodelswhichincludeaHernquist
single galaxy with bulge and disk calculated with the orig- (1990) haloand anexponentialdiskscaled appropriatelyfor
4 Lietal.
redshift (Robertsonetal. 2006), and adjusted to account for
mass accretion through minor mergers. Each of these eight
galaxy progenitors has a black hole seed assumed to orig-
inate from the remnants of the first stars (Abeletal. 2002;
Bromm&Larson 2004; Tan&McKee 2004; Yoshidaetal.
2003,2006;Gaoetal.2007).
The simulations were performed using the parallel, N-
body/Smoothed Particle Hydrodynamics (SPH) code GAD-
GET2 (Springel 2005), which conserves energy and
entropy using the variational principle formulation of
SPH (Springel&Hernquist 2002), and which incorpo-
rates a sub-resolution model of a multiphase interstellar
medium to describe star formation and supernova feedback
(Springel&Hernquist2003a). Starformationismodeledfol-
lowingtheSchmidt-KennicuttLaw(Schmidt1959;Kennicutt
1998a).Feedbackfromsupernovaeiscapturedbyaneffective
equation of state for star-forming gas (Springel&Hernquist
2003a). A prescription for supermassive black hole growth
and feedback is also included, where black holes are rep-
resented by collisionless “sink” particles that interact gravi-
tationally with other componentsand accrete gas from their
surroundings. The accretionrate is estimated from the local
gas density and sound speed using a spherical Bondi-Hoyle
(Bondi&Hoyle1944; Bondi1952) modelthatis limited by
the Eddington rate. Feedback from black hole accretion is
modeledasthermalenergyinjectedintothesurroundinggas
(Springeletal. 2005a; DiMatteoetal. 2005). We note that
implementationsofourmodelforblackholegrowthandfeed-
backthatdonotexplicitlyaccountforEddington-limitedac-
cretion achieve similar results to the method employed by
us (e.g. compare the works of DiMatteoetal. 2007 and
Sijackietal.2007).
These hydrodynamic simulations adopted the ΛCDM
model with cosmological parameters from the first year
Wilkinson Microwave Anisotropy Probe data (WMAP1,
Spergeletal.2003),(Ωm,Ωb,ΩΛ,h,ns,σ8)=(0.3,0.04,0.7,
0.7,1,0.9).Inthispaper,weusethesameparameters.
2.1.2. HierarchicalAssemblyoftheQuasarSystem FIG. 1.—Evolutionofthestarformationrate, blackholeaccretionrate,
andmassesoftheblackholesandstars,respectively,ofthesimulatedquasar
In the simulation analyzed here, the quasar host galaxy
systematz∼6.5,adoptedfromLietal.(2007).
buildsuphierarchicallythroughsevenmajormergersofgas-
rich progenitorsbetween z=14.4- 6.5. Gravitational inter-
actionsbetweenthemerginggalaxiesformtidaltails, strong
shocksandefficientgasinflowsthattriggerintensestarbursts.
Thehighlyconcentratedgasfuelsrapidblackholeaccretion. sities bring the SMBH accretion and feedback to a climax.
Figure1 showsthe evolutionof some aspectsof the system. The black hole reaches a mass of 2 109M , and has
⊙
∼ ×
Betweenz 14–9,themerginggalaxiesarephysicallysmall a peak bolometric luminosity close to that of J1148+5251.
∼
and the interactions occur on scales of tens of kiloparsecs. Blackholefeedbackthendrivesapowerfulgalacticwindthat
Byz 9-7.5,whenthelastmajormergerstakeplace,thein- clears the obscuring material from the center of the system.
∼
teractions have increased dramatically in strength. Galaxies TheSMBHbecomesvisiblebrieflyasanoptically-luminous
arelargelydisruptedincloseencounters,tidaltailsofgasand quasar similar to J1148+5251. Once the system relaxes, the
starsextendoverhundredsofkiloparsecs,andpowerfulbursts SMBH and the host satisfy the relation, M 0.002M ,
BH star
≈
of star formation are triggered, resulting in an average star similarto thatmeasuredin nearbygalaxies(Magorrianetal.
formationrateof 103M⊙yr- 1thatpeaksatz 8.5.During 1998;Marconi&Hunt2003),asaresultofco-evalevolution
∼ ∼
thisphase,theblackholesareheavilyobscuredbycircumnu- of both components (Lietal. 2007; Robertsonetal. 2006;
clear gas. The luminosity from the starbursts outshines that Hopkinsetal.2007c).
from the accreting black holes. So, we refer to this period After z<6 (the “post-quasar phase”), feedback from star
(z 14- 7.5)asthe“starburstphase.” formation and the quasar quenches star formation and self-
∼
Oncetheprogenitorshavecoalesced,themultipleSMBHs regulatesSMBHaccretion.Consequently,bothstarformation
from the galaxies merge and grow exponentially in mass and quasar activity decay, leaving behind a remnant which
and feedback energy via gas accretion. During this period rapidlyreddens. TheobjectwilleventuallyevolveintoacD-
(z 7.5- 6,hereafterreferredtoas“quasarphase”),theblack likegalaxybythe presentday. (Foran overviewofthissce-
∼
holeluminosityoutshinesthatofthestars. Atredshiftz 6.5, nario,seee.g.,Hopkinsetal.2007a,b.)
≈
whenthegalaxiescoalesce,theinducedhighcentralgasden- The photometric calculations by Robertsonetal. (2007)
ModelingDustinz 6Quasars 5
∼
expensive(Bjorkman&Wood2001;Pascuccietal.2004).
Bjorkman&Wood(2001)proposedasolutiontothisprob-
lembyformulatingan“immediatereemission”algorithm,in
whichthedusttemperatureisimmediatelyupdateduponab-
sorptionofaphotonpacket,andthefrequenciesofre-emitted
photonsaresampledfromaspectrumthattakesintoaccount
themodifiedtemperature. Theadvantageofthisapproachis
that dust radiative equilibrium and the radiative transfer so-
lutions are obtained simultaneously without iteration. This
algorithmisdescribedindetailinBjorkman&Wood(2001);
herewebrieflyoutlinethesteps.
AssumingthateachphotonpacketcarriesanenergyE ,and
γ
after absorptionof N packetsin the i- thgrid cell, the total
i
energyabsorbedinthecellis
Eabs=NE . (1)
i i γ
In radiative equilibrium, this energy must be re-radiated,
witha thermalemissivityof j =κ ρB (T),whereκ isthe
ν ν ν ν
absorptive opacity, ρ is the dust density, and B (T) is the
ν
PlanckfunctionattemperatureT,
FIG. 2.—Testrun ofthe radiative equilibrium algorithm adopted from
Bjorkman&Wood(2001)on3-DCartesiancoordinatesforasimplegalaxy 2h ν3
with a bulge and a disk. The input spectrum of the stars is a black- B (T)= P , (2)
body with a temperature of15000 K (black curve). The red curve is our ν c2 ehpν/(kT)- 1
output, the blue curve is simulation data from Kenny Wood, while the
where h is Planck’s constant, c the speed of light, and k is
crosses are observations of a lenticular galaxy NGC 5866, which has an P
unusual extended dust disk seen exactly edge-on. Both Wood’s data and Boltzmann’sconstant.Theemittedenergyinthetimeinterval
the observations come with the code release package from Kenny Wood ∆t is
(http://www-star.st-and.ac.uk/∼kw25/research/montecarlo/gals/gals.html).
Eem=4π∆t dV ρκ B (T)dν
i i ν ν
Z Z
=4π∆tκ (T)B(T)m , (3)
P i i i
showthatthisquasarsystemsatisfiesavarietyofphotometric where κ = κ B dν/B is the Planck mean opacity, B =
P ν ν
selectioncriteriabasedonLyman-breaktechniques.Themas- σT4/π is theRfrequency integrated Planck function, and mi
sivestellarspheroiddescendedfromthesez 6quasarscould isthedustmassinthecell.
∼
bedetectedatz 4byexistingsurveys,whilethegalaxypro- Equatingtheabsorbed(1)andemitted(3)energies,weob-
∼
genitorsat higher redshiftswill likely require futuresurveys tainthedusttemperatureasfollowsafterabsorbingN pack-
i
oflargeportionsofthesky(&0.5%)atwavelengthsλ&1µm ets:
owingtotheirlownumberdensities. σT4= NiL , (4)
i 4N κ (T)m
2.2. ART2: All-wavelengthRadiativeTransferwithAdaptive γ P i i
RefinementTree whereNγ isthe totalnumberofphotonpacketsin the simu-
lation,andListhetotalsourceluminosity. Notethatbecause
ART2 is based on a unification of the 3-D Monte Carlo
the dust opacity is temperature-independent, the product
radiative equilibrium code developed by Bjorkman&Wood
κ (T)σT4 increases monotonically with temperature. Con-
(2001)andWhitneyetal.(2003b),anadaptivegrid,amulti- P i i
sequently,T alwaysincreaseswhenthecellabsorbsanaddi-
phaseISMmodel,andasupernova-origindustmodel.Below, i
tionalpacket.
wedescribeourimplementationofART2.
Theaddedenergyto be radiatedowingto the temperature
increase∆T isdeterminedbya temperature-correctedemis-
2.2.1. MonteCarloRadiativeTransferforDustinRadiative sivity∆j inthefollowingapproximationwhenthetempera-
Equilibrium ν
tureincrease,∆T,issmall:
The Monte Carlo RT method follows the propagation,
dB (T)
scattering, absorption, and reemission of "photon packets" ∆j κ ρ∆T ν . (5)
ν ν
(groups of equal-energy, monochromaticphotons that travel ≈ dT
in the same direction), by randomly sampling the various There-emittedpackets,whichcomprisethediffuseradiation
probability distribution functions that determine the optical field,thencontinuetobescattered,absorbed,andre-emitted
depth, scattering angle, and absorption rates. For the dusty untilthey finally escape from the system. This methodcon-
environmentsweareconcernedwith,there-emittedspectrum servesthetotalenergyexactly,anddoesnotrequireanyiter-
depends on the temperature of the dust, which is assumed ation as the emergent SED, νL =κ B (T), corresponds to
ν ν ν
to be in thermal equilibrium with the radiation field. In the the equilibriumtemperature distribution (Bjorkman&Wood
traditionalscheme,thefrequenciesofthere-emittedphotons 2001).
aresampledfromlocalreemissionspectrawithfixedtemper- ThisMonteCarloradiativeequilibriumcodeworksasfol-
atures. Thedustradiativeequilibriumisensuredbyperform- lows: first, the photonpacketsare followedto randominter-
ingtheMonteCarlotransferiterativelyuntilthedusttemper- actionlocations,determinedbytheopticaldepth. Thenthey
ature distribution converges. Such methods often require a are either scattered or absorbed with a probability given by
largenumberof iterationsandare thereforecomputationally thealbedo(a=n σ /(n σ +n σ ),wherenandσarenumber
s s s s a a
6 Lietal.
density and cross section for either scattering or absorption, The adaptive-mesh refinement grid serves as an efficient
respectively). If the packet is scattered, a randomscattering treeformassassignmentowingtofastneighborfindingwithin
angle is obtained from the scattering phase function. If in- the grids for SPH smoothing. After the grid is constructed,
steadthepacketisabsorbed,itisreemittedimmediatelyata the gas properties at the center of each grid cell, such as
new frequency determined by the envelope temperature, us- density,temperature,andmetallicity,arecalculatedusingthe
ingthealgorithmdescribedabove. Aftereitherscattering,or SPHsmoothingkernel(Hernquist&Katz1989)oftheorigi-
absorptionplusreemission, thephotonpacketcontinuestoa nalsimulation. Allphysicalquantitiesareassumedtobeuni-
newinteractionlocation.Thisprocessisrepeateduntilallthe formacrossasinglecell.
packets escape the dusty environment. Upon completion of Figure 3 gives an example of the adaptive grid applied to
the Monte Carlo transfer, the codeproducesemergentSEDs a snapshot from the galaxy merger simulations in Lietal.
andimagesatanygiveninclinationforawiderangeofbroad- (2007), and compared with a uniform grid. This particu-
bandfilters,includingthoseofHubbleSpaceTelescope(NIC- larsnapshotrepresentsthetimewhenthesystemreachesthe
MOS bands), Spitzer (IRAC and MIPS bands), and SCUBA peakquasarphase,whenthegalaxiesareinthefinalstagesof
submillimeterbands. coalescence. Thesystemishighlydynamicalasmuchgasis
Bjorkman&Wood(2001)testedthisalgorithmextensively falling into the center, while feedbackfromthe centralmas-
by comparing a set of benchmark calculations to those of siveblackholedrivesanoutflow. Thegasdistributionisthus
Ivezicetal.(1997)forsphericalgeometry. Baesetal.(2005) inhomogeneous.
criticallystudythefrequencydistributionadjustmentusedby AsisapparentinFigure3,auniformgridof503(toppanel)
Bjorkman&Wood (2001), and give a firm theoretical basis barelycapturesthedensitydistributionwitha spatialresolu-
fortheirmethod,althoughitmayfailforsmalldustgrains. tion of 4kpc. The resolution is worse than an adaptive grid
One drawback of this code, however, is its fixed geome- withamoderatemaximumrefinementlevel(i.e.,RL=5,mid-
try andlimited dynamicrange. Itcan handleonly 2-Dor 3- dle panel), which has a finest cell length of 3.1kpc. The
∼
D spherical-polargrids, which are not suitable for capturing gridwithamaximumrefinementlevelof12showninthebot-
thearbitrarygeometryandlargedynamicrangescharacteris- tom panel fully captures the large dynamic range of the gas
ticofmerginggalaxies,whichhavemultipledensitycenters, density distribution in three dimensions. It has a minimum
and where gas and stars extend hundredsof kpc and shocks celllengthof 24.4pcforthefinestcells, whichiscompa-
∼
producehighlycondensedgasonscalesofpcs. Wehavede- rable to the spatial resolution of the original hydrodynamic
velopedanimprovedversionofthiscodebyimplementingan simulationofLietal.(2007). Thisresolutionisequivalentto
adaptiveCartesiangridontopofthecodereleasedbyWood1. a 100003 uniformgrid, which is impracticalwith existing
Toensurethatourimplementationofthe algorithmiscor- co∼mputationalfacilities.
rect, we rerunthe test problemthat comeswith the code re-
lease packageof Wood using a uniformCartesian grid. The
2.2.3. AMultiphaseModelfortheInterstellarMedium
testproblemconsistsofasimplegalaxywithbulgeanddisk
components. As is shown in Figure 2, our code reproduces In determining the dust distribution, we adopt the multi-
theresultofWoodverywell. phase model of Springel&Hernquist (2003a) for the ISM.
The ISM is then comprised of condensed clouds in pres-
2.2.2. Adaptive-meshRefinementGrid sure equilibrium with an ambient hot gas, as in the picture
The SPH simulations using GADGET2 (Springel 2005) of McKee&Ostriker (1977). In the hydrodynamic simula-
outputhydrodynamicinformationasparticledata. However, tions, individual SPH particles represent regions of gas that
theraytracingintheRTcalculationisdoneonagrid.There- containcold,densecoresembeddedinahot,diffusemedium.
fore, it is necessary to interpolate the particle-based density ThehotandcoldphasesoftheISMco-existinpressureequi-
field onto a grid. Jonsson (2006) performed radiative trans- librium but have different mass fractions and volume filling
fercalculationsonSPH simulationsofgalaxymergersusing factors(i.e., the hot-phasegasis 10%in mass but &99%
≤
anadaptivegridasimplementedinhisMonteCarloRTcode, in volume). In our previous studies, Lietal. (2007) and
SUNRISE.Unfortunately,self-consistentcalculationsofdust Hopkinsetal. (2006) used only the hot-phase density to de-
radiative equilibrium and emission are not yet included in terminee.g.theobscurationofthecentralAGN,asthemajor-
SUNRISE. ityofsightlineswillpassthroughonlythiscomponentowing
Our algorithmforconstructingadaptivegridsis similar to toitslargevolumefillingfactor. However,thismethodgives
thatofJonsson(2006). Wetypicallystartwithabasegridof onlyaneffectivelowerlimitonthecolumndensity.
a 43 boxcoveringthe entiresimulationvolume. Eachcellis Here,weconsidertwocomponentsofthedustdistribution,
then adaptively refined by dividing it into 23 sub-cells. The havingdifferentdust-to-gasratios and beingassociated with
refinement is stopped if a predefined maximum refinement the two phases of the ISM. Within each grid cell in the RT
level, RL, is reached, or if the total number of particles in calculation, the hot gas is uniformly distributed, while the
thecellbecomeslessthanacertainthreshold,whichevercri- cold, dense cores are randomly embedded. Because of the
terion is satisfied first. The maximum refinementlevel used muchhigherdensityandsmallervolumeofthecoldphase,it
inthepresentworkis12,andthemaximumparticlenumber isimpracticaltoresolvethesecloudseitherinhydrodynamic
allowedinthecellis32,halfthenumberoftheSPHsmooth- simulations or in our radiative transfer calculations. Con-
ingkernelneighborsusedintheGADGET2simulations.The sequently, we implement an observationally-motivated,sub-
resolution of the finest level is therefore L =L /2(RL+1), resolutionprescriptiontotreatthecoldclouds,constrainedby
min box
where L is the box length, and RL is again the maximum theobservedmassspectrumandsizedistributionofmolecular
box
refinementlevel. Forexample,fortheparametersusedinthis cloudsingalaxies.
simulation,L =200kpc,RL=12,wehaveL =24.4pc. Observationsofgiantmolecularclouds(GMCs)ingalaxies
box min
showthatthe GMCsfollowsimplescaling relations(Larson
1http://www-star.st-and.ac.uk/∼kw25/research/montecarlo/gals/gals.html 1981), namely a power-law mass distribution, dN/dM
∝
ModelingDustinz 6Quasars 7
∼
FIG. 3.—ExampleoftheadaptivegridappliedtoasnapshotofthegalaxymergerfromLietal.(2007). Fromtoptobottomis: uniformgrid,andadaptive
gridswithmaximumrefinementlevelsofRL=5,andRL=12,respectively. ThelabelGNindicatesthetotalgridnumbercoveringtheentirevolume(incaseof
adaptivegrids,itprovidesthenumberofthefinestgridandthatofthetotalgrid). Thisplotdemonstratesthatadaptivegridswithsufficientrefinementcapture
thedensitydistributionwell,whiletheuniformgridfailstodosounlessanunreasonablylargenumberofcellsisemployed. Throughoutthispaper,weusea
maximumrefinementlevelof12,whichhasagridresolutioncomparabletothatofthehydrodynamicsimulationsinLietal.(2007).
8 Lietal.
M- 2(e.g.,Fuller&Myers1992;Ward-Thompsonetal.1994; Integratingoverthecloudradius,oneobtains
Andreetal. 1996; Blitz&Rosolowsky 2006), as well as a
p1o9w85e;r-Dlaawmmeeatssa-lr.a1d9iu8s6;reSlactoivoinlleMet∝alR.219(8.e7.g;,SSoalonmdeornseettaall.. N=πLCR13- 3γ-- γR03- γ. (16)
1987;Rosolowsky2005,2007). Ithasbeensuggestedbythe-
Therefore,theaveragedistancethephotonmusttraveltohita
oreticalmodelingthatthemassfunctionisproducedbyturbu-
coldcloud(themeanfreepath)isgivenby
lenceinself-gravitatingclouds(e.g.,Elmegreen&Falgarone
1996; Elmegreen 2002; Ballesteros-Paredes&MacLow 3- γ
L = . (17)
2002; Lietal. 2003), while the mass-radius relation is at- m πC R3- γ- R3- γ
tributedtovirialequilibrium(Larson1981). 1 0
(cid:16) (cid:17)
Hereweincorporatethesetwoempiricalrelationsintoour
InthedustRTcalculation,weassumethatdustisassociated
ISM model for the RT calculations. For a given cell, the
withboththe coldandhotphase gasesthroughcertaindust-
two-phase break-down in Springel&Hernquist (2003a) de-
to-gasratios. Whenaphotonentersacell,wefirstdetermine
terminesthehotandcoldphasegasdensityaccordingtopres-
thedistanceL ittravelsinthehotphasegasbeforehittinga
sureequilibrium. Thecoldcloudsareassumedtofollowthe h
coldcloud,whichisanexponentialdistributionfunction
Larsonscalingrelations:
1 L
dn p(L)= exp - . (18)
=AM- α, (6) Lm (cid:18) Lm(cid:19)
dM
M=BRβ, (7) Therefore, Lh =- Lmlnξ, where ξ is a random number uni-
formlydistributedbetween0 and 1. Theradiusofthe cloud
where dn/dM is the number density of the clouds differen- thephotonjusthitsisalsodeterminedrandomlyassumingthe
tial in cloud mass M, and R is the cloudradius. From these distributionfunctionofEquation(15);i.e.,
equations,oneobtainsthecloudsizedistribution
dn =βAB1- αR- (αβ+1- β) R= R03- γ+(R13- γ- R03- γ)ξ 1/(3- γ). (19)
dR h i
=CR- γ, (8) ThedistanceLcthephotontravelsinthiscoldcloudisagain
whereC=βAB1- α,andγ=αβ+1- β. a random variable given by Lc =2R√ξ, because clouds are
assumed to be uniformly distributed. These equations com-
Assume that the minimum and maximum values of the
pletely define the statistical procedure for determining the
cloud mass are M and M , then the normalization constant
0 1 dustcolumndensitiesassociatedwithL andL as:
ofthemassspectrumforeachcellcanbedeterminedusing h c
dn N =ρ L
h h h
MdM=x ρ , (9)
Z dM c c 3BRβ- 3L
or Nc= c. (20)
2- α 4π
A=x ρ , (10)
c cM2- α- M2- α Withagivendustopacitycurve,theseequationsallowone
1 0 tocalculatetheopticaldepths,τ andτ forthehotandcold
h c
where ρc and xc are the cold gas density and volume filling dust,respectively. Theyrelatethephotonpathlengthsinthe
factor, respectively. The normalization constant of the M-R multiphase ISM to the total optical depth τ , which is then
tot
relationmaybedeterminedusing
compared with a randomly drawn number τ, to determine
i
4π dn whether the photon should be stopped for scattering or ab-
x = R3 dM, (11)
c Z 3 dM sorption.Indetail,theMonteCarloray-tracingprocedurefor
or theradiativetransferthereforeinvolvesthefollowingsteps:
4πA β/3
B= Mη- Mη , (12) 1. Aphotonpacketisemittedfromeitherablackholeora
(cid:20)3ηx 1 0 (cid:21)
c(cid:0) (cid:1) stellarsourcewithrandomfrequenciesconsistentwith
where the source spectra. The photon is emitted with a uni-
3
η=1+ - α. (13) formlydistributedrandomdirection.Theprobabilityof
β a photon being emitted by any given source is deter-
Theminimumandmaximumcloudradiiinthecellarethere- minedbyitsluminosityrelativetothetotal.
fore
2. A random optical depth over which the photon must
R = M0 1/β travel before an interaction with the dust occurs, τi =
0 (cid:18) B (cid:19) - lnξ, is drawn to determine the interaction location.
The interaction includes scattering and absorption. In
M 1/β
R = 1 . (14) ourmethod,thephotonenergiesarenotweighted,only
1
(cid:18) B (cid:19) one eventis allowed. Thatis, at anygiven interaction
ForaphotontravelingadistanceLinthecell, theaverage site, thephotoniseitherscatteredorabsorbed,butnot
numberofcoldcloudsofradiusRthephotonwillintersectis both.
givenby
3. Starting from the location of the photon emission, the
dN dn
=πR2L cumulative optical depth of the photon, τtot, is calcu-
dR dR latedstochasticallyusingEquation20forbothhotand
=πCR2- γ. (15) cold dusts. If the photon is stopped for interaction
ModelingDustinz 6Quasars 9
∼
within a single cell, then τ is the sum of contribu-
tot
tionsfrompossiblymultiplesegmentsofbothhotand
colddustswithinthiscell. Ifthephotonpassesthrough
multiple cells before an interaction occurs, then τ is
tot
the sum of allcontributionsfromrelevantsegmentsin
thesecells.
4. At each boundarybetween the hotand cold phase gas
clouds, or at the boundary of the grid cell, the next
interaction point is determined by the comparison be-
tween τ and τ . If τ τ , then the photon is ei-
i tot i tot
≤
ther scatteredorabsorbed,witha probabilitygivenby
the scattering albedo. The exact interaction location
is then determinedinside either hotor coldphase gas,
suchthatτ becomesexactlyτ. Ifthephotonisscat-
tot i
tered, its direction and polarization state are altered
using the Henyey-Greenstein phase function, and the
ray-tracingof the new photonis repeatedfromstep 2.
If the photon is absorbed, depending on whether the
absorption site occurs in the hot or cold phase dust,
the temperature of the appropriate dust cell is raised,
FIG. 4.— Comparison of the dust absorption opacity curves from a
andanewphotonisreemittedaccordingtothescheme supernova-origindustmodel(SNmodel,withsilicatefraction fs=10%,20%
of Bjorkman&Wood (2001). The ray-tracing of the in red and black, respectively) and those of Weingartner&Draine (2001)
newlyemittedphotonagainrestartsfromstep2. (WD01 model, with RV =3.1, 4.0, 5.5, in blue, cyan, and green, respec-
tively).Notetherearetwodifferencesbetweenthesetwomodels:thesilicate
5. Ifthephotonescapesfromthesystemwithoutreaching featureat∼9.7µm,andahigheropacityintheopticalband.
the optical depth τ, it is then collected in the output
i
spectrumandimage.Thenextphotonwillbepickedup
fromthesource,andthewholeMonteCarloprocedure
fromstep1willberestarted.
We assume the cold clouds to be in the mass range of
M =103M and M =107M , which is similar to that of
AfterallN photonshavebeentraced,oneobtainsthedust 0 ⊙ 1 ⊙
γ
protoclusters clouds in star forming galaxies, as shown in
temperaturedistributionforbothhotandcolddustsinradia-
Lietal. (2004, 2005a,b, 2006). For the cold clouds, we
tive equilibrium, as well as the output spectra and images.
have enhanced the pressure by a factor of 10 over the ther-
Note that such a stochastic procedureoutlinedhere may not
mal pressure to account for the effects of turbulence (e.g.,
be as accurate as physically tracking the locations and sizes
Blitz&Rosolowsky2006;Chakrabartietal.2007a).Thisen-
ofthecoldphasegasclouds,whichcanberatherdifficult,if
hancementfactor is referredto as “CP” hereafter. The dust-
notimpossible. However,thismethodworksefficientlyfora
to-gas ratio of the cold clouds is chosen to be the same as
largenumberofcoldcloudsuniformlydistributedinthecells,
theMilkyWayvalue(1:124;Weingartner&Draine2001),as
anditgivescorrectaverageextinctionpropertiesforamulti-
found in a large sample of ULIRGs (Dunne&Eales 2001;
phaseISMmodelweareinterestedhere.
Klaasetal.2001),whilethatofthehot,diffusegasischosen
Hereafter,wereferto“HPG-dust”asthedustthatoriginates
to be 1% of the Milky Way value, consistent with the dust
from the hot-phase gas, and “CPG-dust” as that originating
survivalrateofsputteringinahot,diffuseISM(Burke&Silk
from the cold-phase gas. Note that the gas only determines
1974;Reynoldsetal.1997). Observationsofsomeobscurred
thedistributionandmassofthedustthroughagivendust-to-
orredAGNsalsosuggestaratiosignificantlylowerthanthat
gasratio. Thedusttemperatureisnotassociatedwiththegas
ofMilkyWay(e.g.,Maiolinoetal.2001;Kuraszkiewiczetal.
temperature, and is calculated self-consistently according to
2003;Halletal.2006). Wewillperformasystematicparam-
radiativeequilibrium. Wefurtherreferto“colddust”,“warm
eterstudyofthesechoicesin§4. However,wenotethatthese
dust” and “hot dust” as dust with temperatures T .100 K,
valuesare foundtobestreproducetheobservedquasarSED
100.T .1000K,and1000.T .1200K,respectively.
IntheRTcalculations,weadoptamassspectrumwithα= at z 6.5. Because the dust opacity is proportional to the
∼
gasmetallicity, thereforethedustopacityis weightedbythe
1.8,assuggestedbytheobservationsofBlitz&Rosolowsky
metallicityofthegasforbothhotandcoldgasphases.
(2006), and an observed mass-radius relation with β =2.0,
which is also a result of the virial theorem. The resulting
mass-radiusrelationinoursimulationshasanormalizationin 2.2.4. Supernova-originDustModel
therange 10- 104M⊙/pc2. Foracloudsizeof 1pc,the Ourunderstandingoftheformationanddistributionofdust
∼ ∼
normalizationis 300M⊙/pc2,similartoobservationsofthe hasbenefitedfromdustmapsof theMilkyWay (e.g., Gehrz
∼
Milky Way (Scovilleetal. 1987; Rosolowsky 2007). It has 1989;Schlegeletal. 1998). In particular,Gehrz(1989)con-
beenshownthatthenormalizationofthemass-radiusrelation- cludesfromadetailedGalacticsurveyofdust-producingstars
shipdependsongalacticenvironment(e.g.,Elmegreen1989; thatdustintheMilkyWayoriginatesinthreeprincipleways:
Rosolowsky2005;Blitzetal.2007). Forexample,itisabout (1)bycondensationinwindsofevolved,post-main-sequence
50M⊙/pc2intheLargeMagellanicCloudbutcouldbetwo objects, which accountsof 90%of the stellar dust; (2) by
∼ ∼
ordersofmagnitudehigherinULIRGs. Inextremestarburst condensation in ejecta from massive novae, supernovae and
galaxies,thenormalizationcouldgoupto104M⊙/pc2,aswe Wolf-Rayetstars,whichamountsto<10%;and(3)byslow
findinoursimulations. accretioninmolecularclouds(seealsoMarchenko2006fora
10 Lietal.
2003), and SNe 2003gd (Sugermanetal. 2006). Theoreti-
cally, several groups have calculated dust formation in the
ejecta of Type-IISNe (Todini&Ferrara2001; Nozawaetal.
2003;Schneideretal.2004;Hirashitaetal.2005;Dweketal.
2007). Inparticular,Todini&Ferrara(2001)havedeveloped
a dustmodelbasedonstandardnucleationtheoryandtested
it on the well-studied case of SN1987A. Theyfind that SNe
withmassesintherangeof12- 35M produceabout1%of
⊙
themassindustpersupernovaforprimordialmetallicity,and
3%forsolarmetallicity.
∼
In the presentwork, we adoptthe dust size distributionof
Todini&Ferrara(2001)forsolarmetallicityandaM=22M
⊙
SNmodel,asinFigure5intheirpaper.Thissizedistribution
isthencombinedwiththedustabsorptionandscatteringcross
sections of Weingartner&Draine (2001), to calculate dust
absorptionopacitycurves(Finkbeiner1999;Finkbeineretal.
1999). WenotethatBianchi&Schneider(2007)re-analysize
the model of Todini&Ferrara (2001) and follow the evolu-
tionofnewlycondensedgrainsfromthetimeofformationto
theirsurvival. Thisnewfeatureshowsbetteragreementwith
FIG. 5.—DistributionoftheopticaldepthsofbothHPG-andCPG-dust, observationsandfurthersupportsourmotivationtouseaSN
respectively,atλ=0.1µmwherethedustopacitypeaksasshowninFigure4. dust model. Figure 4 shows dust absorption opacity curves,
Thedensitygridcorrespondstothesimulationsnapshotatz=6.5shownin in the range 10- 2- 104µm, from both the supernova-origin
Figure3(bottompanel). The“HPG-dust”and“CPG-dust”refertothedust
dust model (hereafter SN model) and Weingartner&Draine
associatedwithhot-phaseandcold-phase gas,respectively, asdescribedin
§2.2.3. (2001) (hereafter WD model), respectively. The SN models
include silicate fractions f =10% (in red) and f =20% (in
s s
black),whichshowslightdifferencesinthesilicatefeatureat
9.7µm(f =20%modelhasaslightlyhigheropacity).The
s
∼
WD models include all the curves commonly used with ex-
review). tinctionR = A(V) =3.1,4.0,5.5,whichdifferintheopac-
V A(B)- A(V)
Over the past several decades, various dust models have ityintheUV-NIR(0.01µm- 1µm)bands(R =3.1modelhas
V
been developed based on the observed extinction curves thehighestopacity).
of the Large Magellanic Cloud (LMC), the Small Magel- BetweentheSNandWDmodels,therearetwonoticeable
lanic Cloud (SMC) and the Milky Way (e.g., see reviews differences; one is the silicate feature at 9.7µm, and the
bySavage&Mathis1979;Mathis1990;Calzettietal.1994; secondistheUV-NIRband.Ononehand,∼theWDmodelhas
Dorschner&Henning 1995; Calzettietal. 2000; Whittet strongpeaksaround 9.7µm(silicatefeature),whiletheSN
2003; Draine 2003). In these classical models, dust is modelproducesanop∼acitylowerbynearlyoneorderofmag-
assumed to form in the envelopes of old, low-mass stars nitude,becausetheejectaofSNewithmassabove15M pre-
⊙
such as asymptotic giant branch (AGB) stars with ages &1 dominantlyformamorphouscarbongrains, which decreases
Gyrs (Mathis 1990; Morgan&Edmunds 2003; Dwek 2005; thesilicatefractioninthedustgrains. Ontheotherhand,the
Marchenko2006). SNmodelincreasestheopacityintheopticalbandbya fac-
However, this picture may be different for young, high- torofafewowingtothesmallergrainsize. Wenotethatthe
redshift objects. Recent deep millimeter and submillime- smallgrainsizesmaycausequantumfluctuationsinthetem-
ter observations of several z 6 quasars, which trace the perature and high-temperature transients (e.g., Draine&Li
∼
far-infraredthermaldustemissionfromthese systems, show 2001). However,thedustproducedbytheSNmodelisdom-
large masses of dust in these quasar hosts when the Uni- inatedbygraphitegrains,whichhaveasizedistributionfrom
verse was less than 1 Gyr old (e.g., Bertoldietal. 2003a; 100- 2000 Å, comparable to the size range of the dust
Charmandarisetal. 2004; Robsonetal. 2004; Carillietal. ∼
grainsintheWDmodel,sothetemperaturefluctuationisex-
2004; Maiolinoetal. 2004; Beelenetal. 2006; Hinesetal.
pectedtobeinsignificant. Notewedonotincludepolycyclic
2006; Jiangetal. 2006; Willottetal. 2007). In particular,
aromatichydrocarbon(PAH)featuresat 7.7µm,whichwas
Maiolinoetal. (2004) find that the extinction curve of the ∼
modeled in detail by Li&Draine (2002) and Draine&Li
reddenedquasarSDSS J1048+46atz=6.2is differentfrom
(2007). The PAH feature is a good diagnostic of starbursts
those observedat z<4 (similar to that of the Small Magel-
atlow redshifts(Houcketal. 2005), howeveritis verydiffi-
lanicCloud,Hopkinsetal.2004),butmatchestheextinction
culttodetectinz&6systems. Alsodustspinisnotincluded
curveexpectedfordustproducedbysupernovae.
in our modeling, which could be significant around 10 mm
Ithasbeensuggestedthatcore-collapsesupernovae(Type-
(Draine&Lazarian1998;Finkbeineretal.2002).
II SNe) may provide a fast and efficient mechanism for
The differences between these two models, in particular
dust production. The observational evidence for dust for-
the silicate feature, may be diminished by choosing differ-
mation in SNe comes from observations of SN1987A (e.g.,
entgraincompositionsandsizedistributions(Laor&Draine
Gehrz&Ney 1987; Moseleyetal. 1989; Rocheetal. 1993;
1993). For example, the WD model would converge to the
Spyromilioetal. 1993; Colganetal. 1994), Cassiopeia A
SN model if the silicate fraction is reduced to 5%. We
(Dunneetal. 2003; Dwek 2004; see however Krauseetal. ∼
notethatElvisetal.(2002)proposeanalternativemechanism
2004 who argue that the dust emission is not associated
for dust productionin quasars. They suggest that the physi-
withtheremnant),Kepler’ssupernovaremnant(Morganetal.