Giant Lyman-Alpha Nebulae in the Illustris Simulation
AAccepted for publication in ApJ 2016 December 16
Preprint typeset using L A TEX style emulateapj v. 5/2/11
GIANT LYMAN-ALPHA NEBULAE IN THE ILLUSTRIS SIMULATION
Max Gronke
Institute of Theoretical Astrophysics, University of Oslo, Postboks 1029 Blindern, 0315 Oslo, Norway
Simeon Bird
Department of Physics and Astronomy, Johns Hopkins University, 3400 N. Charles St., Baltimore, MD 21218, USA
Accepted for publication in ApJ 2016 December 16
AbstractSeveral ‘giant’ Lyman- α (Ly α ) nebulae with extent (cid:38)
300 kpc and observed Ly α luminosity of (cid:38) erg s − cm − arcsec − have recently been detected, and it has been speculated that their presencehints at a substantial cold gas reservoir in small cool clumps not resolved in modern hydro-dynamicalsimulations. We use the Illustris simulation to predict the Ly α emission emerging from large halos( M > . M (cid:12) ) at z ∼ α spectra to a model where mostgas is clumped below the simulation resolution scale. We find that with Illustris no additionalclumping is necessary to explain the extents, luminosities and surface brightness profiles of the ‘giantLy α nebulae’ observed. Furthermore, the maximal extents of the objects show a wide spread for agiven luminosity and do not correlate significantly with any halo properties. We also show how thedetected size depends strongly on the employed surface brightness cutoff, and predict that furthersuch objects will be found in the near future. Subject headings: galaxies: high-redshift – galaxies: intergalactic medium – line: formation – scatter-ing – radiative transfer – quasars: general INTRODUCTION
The medium immediately surrounding galaxies – oftendubbed the circumgalactic-medium (CGM) – provides agas reservoir for star formation and as such is crucial forthe study of galaxy formation and evolution. Extendedfaint Ly α emission originating from these regions directlyprobes this gas, uniquely so at higher redshifts where theobservation of other emission lines is challenging (for re-views see, e.g., Barnes et al. 2014; Dijkstra 2014; Hayes2015). These ‘Ly α nebulae’ are often dubbed ‘Lyman- α blobs’ (LABs), or when fainter and surrounding a galaxyalso called ‘Lyman- α halos’ (LAHs). While only a fewLABs have been found so far (Fynbo et al. 1999; Steidelet al. 2000; Matsuda et al. 2011; Ao et al. 2015), it hasbeen shown that LAHs surround Lyman-break or ‘drop-out’ galaxies (LBGs), galaxies selected through their Ly α emission (Ly α emitters or LAEs) as well as H α selectedgalaxies (recently shown by Matthee et al. 2016) lead-ing to the conjecture that most (if not all) star-forminggalaxies are associated with a LAH. Due to their very lowsurface-brightness (SB) profiles a stacking technique isoften used in order to study the properties of LAHs. Stei-del et al. (2010) stacked a sample of 92 LBGs at z ∼ . ∼ erg s − cm − arcsec − andfound the LAHs extend out to ∼
80 kpc with a exponen-tial scale length of r α ≈ −
30 kpc which is a factorof a few larger than the scale length of the continuumemission. This result is consistent with the more recentfinding of Matsuda et al. (2012) who used a set of 2128LAEs and 24 LBGs at z ∼ . r α ∼ −
30 kpc and r α ∼
20 kpc, respectively. Otherstudies (Rauch et al. 2008; Momose et al. 2014) support [email protected] this picture but suggest a smaller LAH for LAEs as op-posed to LBGs. This could however be an environmen-tal effect if LBGs reside in higher-density environments(Matsuda et al. 2012). Furthermore, ultradeep MUSEobservations recently revealed the LAHs surrounding 21individual galaxies (Wisotzki et al. 2016). Without theneed of stacking their SB profiles, they confirmed thatthe scale of the Ly α emitting region is several times largerthan the scale of the continuum emission.While the Ly α nebulae surrounding ‘normal’ galax-ies extend ∼ tens of kpc and are relatively faint L α ∼ erg s − , LAHs around a rare z (cid:38) L α ∼ erg s − (Reuland et al. 2003; Villar-Mart´ın et al. 2007; Trainor & Steidel 2012; Martin et al.2015; Hennawi et al. 2015; Borisova et al. 2016). Can-talupo et al. (2014) detected the most massive of thesenebulae (sometimes called the ‘slug nebula’), measuringa maximum extent of ∼
450 kpc and L α ∼ erg s − ( L α ∼ × erg s − excluding the emission directlyfrom the quasar).Although the existence of LABs, LAHs and the giantLAHs is observationally well-established, and some as-pects are theoretically understood (e.g., their existenceimplies a fairly large hydrogen mass around dark-matterhalos, Matsuda et al. 2012), some key questions remainunclear. One which is debated in the literature is theenergy source of these – bright and extended – Ly α neb-ulae. Commonly, two possibilities are discussed: (i) Ly α production in the central region and subsequent scatter-ing in the surrounding gas leading to the observed halo(e.g. Laursen & Sommer-Larsen 2007; Dijkstra & Kramer a r X i v : . [ a s t r o - ph . GA ] D ec TABLE 1Overview of the model parameters
Model Parameter Description Fiducial valueAGN r ion Radius of ionized regionaround AGN 20 kpcStars m ion Cell volume ionized inunits of Str¨omgrenspheres (Eq. (3)) 1 (ii) in situ Ly α production in an extended re-gion. The latter would be possible, via e.g. cooling lumi-nosity of gas falling into the potential well of the galaxy(e.g. Haiman et al. 2000; Dijkstra & Loeb 2009). Al-ternatively, ionizing photons escaping from the centralregion or originating from nearby galaxies could lead torecombination events in the surrounding medium (e.g.Haiman & Rees 2001; Furlanetto et al. 2005; Mas-Ribas& Dijkstra 2016).In particular, the discovery of the giant LAHs posesthe question of how Ly α emission can be located so farfrom an ionizing source. Cantalupo et al. (2014) carriedout radiative transfer simulations on a zoom-in hydro-dynamical simulation and concluded that their modelcould not explain a Ly α object having this SB level andextent. Instead they proposed that hydrodynamical sim-ulations miss a substantial fraction of cold, dense clumpswhich will boost the Ly α luminosity. Specifically, theycalculate the required clumping factor C ≡ (cid:104) n e (cid:105) / (cid:104) n e (cid:105) (where n e is the electron number density) to be between20 and 1000 on scales below a few kpc. In this workwe want to revisit this question using the more mod-ern cosmological hydro-dynamical Illustris simulation(Vogelsberger et al. 2014; Nelson et al. 2015) which fea-tures a fully realised and well-tested model for galaxy for-mation, tuned to produce a realistic galaxy populationat z ∼
0. This model includes efficient supernova feed-back, metal cooling and, importantly for studying thegas around quasars, a recipe for AGN feedback. Usingthis publicly available data allows us to study a statis-tically relevant ensemble of halos instead of focusing onindividual objects.The study is structured as follows. In Sec. 2 we layout our model and numerical methods. We present ourresults in Sec. 3, and conclude in Sec. 4. METHODS
In this section, we first describe the hydrodynamicalsimulation used for our work ( § § α radiative transfer simulation employed ( § The hydrodynamical simulation
The
Illustris simulation (Vogelsberger et al. 2014)is a hydro-dynamical cosmological simulation (simulationbox side length of 106 . AREPO (Springel 2010). It includes gascooling, and photo-ionization from a uniform ultra-violetbackground. Subgrid models are included for black holes& black hole feedback, stochastic star formation (witha density threshold of 0 .
13 cm − above which an ad-hocequation of state is imposed), stellar evolution, and stel-lar feedback, tuned to reproduce the galaxy stellar massfunction at z = 0. Note that molecular cooling is not in-cluded, and so the simulation does not accurately follow the temperature of gas with T < K. In particular, wemade use of the ‘Illustris-1’ simulation which possessesa baryon (dark matter) particle mass of 12 . × M (cid:12) (62 . × M (cid:12) ) and a resulting gravitational softeninglength of ∼
700 pc (1 . Illustris team (Nelson et al. 2015) whichallowed us to post-process individual cutout halos.Inspired by the observations of Cantalupo et al. (2014),we selected a halo with total mass M ∼ . M (cid:12) h − at z ∼
2. This halo possesses an active, central black holewith a mass inflow rate of ˙ M BH ∼ . M (cid:12) yr − , a initialneutral fraction of ∼
26% and a cold (
T < K ) gasfraction of ∼ ∼ . × M (cid:12) h − .The relevant quantities (per particle) for our work are:the cell volume V cell , the neutral hydrogen number den-sity n HI , the ionized hydrogen number density n HII , theelectron density n e , the gas velocity v , the star-formationrate SF R , the metallicity Z , and the temperature T . Wecompute T as T = 2 uµ k B , (1)where u is the internal energy, µ the mean molecularweight of the gas, and k B the Boltzmann constant. Dustplays a crucial role in Ly α radiative transfer and dustreddening is expected to be important in these massivehalos. We follow Laursen et al. (2009) and calculate aneffective number density of dust atoms as n d = ( n HI + f ion n HII ) ZZ (cid:12) (2)which defines a dust optical depth through τ d = n d σ d d ,where σ d is the dust optical depth and d the path lengthconsidered. For the former we use the value of theSmall Magellanic Cloud σ Ly α, SMC ≈ . × − cm Pei (1992) and take furthermore the dust-to-gas ratio inan ionized region to be characterized by f ion = 0 .
01 (see § § § Illustris using the prescriptionof Rahmati et al. (2013a). This includes the neutral frac-tion in the (subgrid) star-forming gas, which is roughlyunity, neglecting the local radiation field of the formingstars.
The model
In addition to the halo (as described in the previoussection), we consider two possibilities for the radiationfield:1.
AGN ionization.
In this scenario, we fully ionize allthe cells falling within a spherical region around thecentral black hole with radius r ion . This clearly is asimplification of the true ionization mechanism byquasars as we neglect radiative transfer effects fromdifferent density structures as well as the likely exis-tence of conical jets. However, in this work we are r ion , AGN (kpc) L L y α , i ( e r g / s ) AGN -2 -1 m ion Stars
Fig. 1.—
Intrinsic Ly α luminosities for the two models. Theblack-dashed line denotes the intrinsic luminosity reached if theentire grid is fully ionised and the red star symbols illustrate thechoice of our fiducial models. From the two leftmost panels, it isclear that in each model reaching the maximal possible luminositywith the gas as given by the hydrodynamical simulation output ispossible with either moderate ionization activity from either thestar forming regions or the black hole. content with a first order approximation as Can-talupo et al. (2014) demonstrated that it was notpossible to resolve the discrepancy between the ob-servations and simulations without extra subgridclumping even in the extreme cases when all thehydrogen is ionized or neutral. Quasars can po-tentially ionize a region up to several Mpc (Cen& Haiman 2000). However, we choose a relativelysmall value for r ion in § (i) show theradiative transfer effect in the outer region (i.e., tobe much smaller than the total extent of the Ly α halo), and (ii) still be large enough so that a smallchange in r ion does not lead to a significant changein the intrinsic luminosity ( § r ion = 20 kpc.2. Ionization from stars.
Here, we ionize all the star-forming cells by subtracting∆ X HI = Q ion m ion n α B ( T ) V cell (3)from the cells’ neutral hydrogen fraction, while en-suring X HI ≥
0. In the above equation, α B isthe case-B recombination coefficient (see below)and Q ion = 2 ×
53 SFR1 M (cid:12) yr for a range of differ-ent stellar models (Rahmati et al. 2013b). Here,the free parameter m ion can be interpreted as thenumber of (fully ionized) Str¨omgren spheres (vol-ume Q ion / ( n α B )) placed in a neutral cell. Thismeans for a cell with X HI = 1 if there is one star(cluster) per cell, the ‘sub-grid’ escape fraction ofionizing photons is unity and if the ionized regiondoes not overlap with the boundary then m ion ∼ −
10 Myr.However, given the uncertainty in calculating theseeffects, and the fact that, as shown in Fig. 1, theLy α luminosity saturates at m ion ≈ . m ion = 1appears to be a reasonable fiducial value and weadopt it for the rest of the paper. These model parameters are summarized in Table 1. Thetotal amount of gas within the simulation is conservedin both the scenarios considered. We assume that theionizing radiation fields are sufficiently intense that thegas will be highly ionized (Hennawi & Prochaska 2013).In fact, we allow the hydrogen neutral fraction in the maximally ionized regions we reach a minimum value of X HI , min = 0. As the emissivity ∝ (1 − X HI ) a slightlylarger value of X HI will not have a strong effect on the SBvalues (however, it can affect the Ly α radiative transfer,see § n HII , theelectron number density n e and a temperature T per cell,we compute the total Ly α luminosity assuming solely‘case-B’ recombination L α, i = (cid:88) cells n HII n e N α ( T ) α B ( T ) V cell , (4)where α B ( T ) is the ‘case-B’ recombination coefficientfrom the fitting formulae by Hui & Gnedin (1997), and N α ( T ) is the average number of Ly α photons producedper ‘case-B’ recombination event. For N α ( T ) we adoptthe fit provided by Cantalupo et al. (2008). As Ly α is produced mainly in high-density environments, thecontribution of ‘case-A’ recombination can be safely ne-glected (Gould & Weinberg 1996). Ly α cooling radia-tion is highly sensitive to the temperature of the coldgas and no consensus has been reached on its contribu-tion (Furlanetto et al. 2005; Goerdt et al. 2010; Faucher-Gigu`ere et al. 2010; Rosdahl & Blaizot 2011; Cen &Zheng 2013). In order to assess its impact a fully cou-pled radiation-hydrodynamical simulation is preferable(as, e.g., in Rosdahl & Blaizot 2011) which can in prin-ciple model the temperature state of the ISM. We do notaddress the impact of cooling radiation in this work, and,thus neglect its contribution. Both choices only reducethe final luminosity. Radiative transfer calculations
Prior to carrying out the full Ly α radiative trans-fer calculation we interpolate the individual halo fromthe hydrodynamical simulation to a Cartesian grid with512 × ×
449 cells (corresponding to a uniform side-length of each cell of ∼ .
63 kpc) while keeping con-stant the number of hydrogen atoms, the number ofdust grains, the number of clouds placed on the gridand the total Ly α luminosity. For the conversion fromthe moving-mesh structure to a Cartesian grid (used by AREPO (Springel 2010) and the radiative transfer calcu-lations, respectively), we used splash (Price 2012) witha smoothing length of l smooth = 4 × V / which cor-responds to ∼
58 neighboring particles. We repeatedour analysis for one model using half the number of cellsper dimension and find the results to be unaffected. Af-ter the conversion the number of clouds in each cell isrounded to the nearest integer, redistributing the cut-off material self-consistently and the clouds are placedrandomly within each cell. ‘case-B’ recombination denotes the recombination in a mediumwhich is optically thick to ionizing photons. This leads to theimmediate re-absorption of an emitted ionizing photon (see, e.g.,Osterbrock 1989; Dijkstra 2014, for details). We used the default M cubic B-spline kernel. l o g S B ( e r g s − c m − a r c s e c − ) kpc k p c L α = 4 . × ergs − Stars, m ion = 1 kpc L α = 3 . × ergs − AGN, r ion , AGN = 20 kpc kpc L α = 2 . × ergs − AGN, r ion , AGN = 20 kpc (no dust)
Fig. 2.—
Surface brightness maps for the fiducial models (from left to right panel): star-driven ionization, AGN-driven ionization, andAGN-driven ionization without dust (see § Distance from BH (kpc) -21 -20 -19 -18 -17 -16 -15 -14 -13 S B ( e r g s − c m − a r c s e c − ) Stars, m ion = 1 AGN, r ion , AGN = 20kpc
AGN, r ion , AGN = 20kpc (no dust)
Fig. 3.—
Radially binned surface brightness profile as a functionof distance from the black hole for the fiducial models as ‘seen’ byan observer directed in the same way as in Fig. 2. The dashed andsolid lines show the (intrinsic) surface brightness profile before andafter the radiative transfer calculations, respectively. The verticalline denotes the virial radius of this halo ( R c ). The Ly α radiative transfer calculations are carried outusing the Monte-Carlo radiative transfer (MCRT) code tlac (Gronke & Dijkstra 2014). General descriptions ofMCRT are given in Dijkstra (2014) or Laursen (2010).The specific settings of tlac employed are identical toGronke & Dijkstra (2016) and we merely summarizethem here. We ran each MCRT simulation using at least2 × photon packages which we placed randomly onthe gas density grid. The probability to position a pho-ton package in a certain cell (i.e., the weight of this cell)was proportional to the cell’s intrinsic luminosity. Wedraw the intrinsic frequency of a photon package from aVoigt profile, the convolution of the natural (Lorentzian)line profile, and the thermal (Gaussian) profile of theemitting atoms of the initial cell. The process of imagemaking with tlac has not been published previously andwe describe the method below.Producing SB maps with MCRT can be challenging asthe number of individual photon packages escaping in aspecific direction is essentially zero. Therefore, we use acommonly used technique sometimes called the ‘peeling’algorithm (Yusef-Zadeh et al. 1984; Zheng & Miralda-´Escude 2002; Laursen 2010). In every N th scatter event(including the emission step), the optical depth, τ obs ( ν v ),in the direction k obs is recorded. Here, ν v is the frequencywhich the photon would have had if it had flown in the direction k obs . This assigns a weight to this scatteringevent which is the escape probability along k obs given by: w = N e − τ obs ( ν v ) , (5)where N is the frequency of recorded scattering eventsdiscussed above. The SB for pixel j can then be calcu-lated using SB j = L α, i S j D L ( z )Ω pix ,j N phot , (6)where N phot is the total number of photon packagesemitted, D L ( z ) is the luminosity distance correspondingto redshift z , Ω pix ,j is the solid angle of pixel j and S j = (cid:80) j w is the sum of weights of the photons fallingwithing this pixel. Naturally, if calculating an observedspectrum, the weights w have to be taken into accountas well, and each frequency bin consists of the sum ofthe photons’ weights. RESULTS
In this section, we present first results from our ‘fidu-cial’ halo ( § § α radiativetransfer simulations for the former part – from whichwe show SB maps and profiles ( § α spectra( § α halos found in the Illustris simulations toobservations.
Intrinsic Ly α luminosities Fig. 1 shows the intrinsic Ly α luminosity for the mod-els presented in § Note, that although we consider the slightly differing redistri-bution functions for the scattering via the 2 P / and the 2 P / states for the full radiative transfer process we assume a uniformangular probability density function for the potential scatteringtowards k obs . However, due to the large amount of photons thatescaped without scattering ( ∼ by the hydrodynamical simulation output. Also notice-able is the proximity of the computed Ly α luminosity( L α, i , max ∼ × erg s − ) and the observed value byCantalupo et al. (2014) ( ∼ × erg s − excluding theemission directly from the quasar). The luminosity satu-rates because the densest regions are already ionised. Forthe ‘stars’ model this occurs at m ion (cid:38) . r ion (cid:38)
15 kpc. For even greater radiithe gas density drops off – and so does the increase inluminosity.
Surface brightness
Fig. 2 shows the surface brightness maps for our twomodels (see § ,
000 recorded events for each ofthe four observers’ directions. In Appendix A we showthe intrinsic SB maps, i.e., without radiative transfer ascomparison.Fig. 2 shows very similar topography in both theAGN and stars ionization models. The stars pro-duce slightly more prominent features in the outskirtsof the halo. Both cases have an maximal extentof ∼
350 kpc for this particular sightline, a centralSB of ∼ − erg s − cm − arcsec − which falls to ∼ − erg s − cm − arcsec − in the outer regions. Over-all, the size of this simulated nebula falls ∼
100 kpc shortof the projected maximum extent measured by Can-talupo et al. (2014). However, if restricted to SB contours (cid:38) − erg s − cm − arcsec − then the measurement isreduced by ∼
20% (Cantalupo et al. 2014) which bringsthe two measures very close to each other. Fig. 2 alsoshows the ‘AGN’ case but without the inclusion of dust(rightmost panel). Here, one can notice that the mor-phology and SB levels are different, that is, in the casewithout dust the halo is extended and overall brighter.Fig. 3 shows a more quantitative comparison of the sur-face brightness between the two models. Here we displaythe SB profiles, as well as the intrinsic SB profiles, as afunction of distance from the black hole, observed fromthe same direction as in Fig. 2. The same features arevisible, with the curves of the ‘stars’ and ‘AGN’ cases fol-lowing each other closely. Note, however, that the SB inthis case is an angular averaged value and thus dependenton the geometry of the halo. This leads to variations inthe radial profile depending on the observer’s direction.The black dashed line in Fig. 3 denotes the ‘virial ra-dius’ R c , i.e., the radius at which the average densityis 200 times the critical density at z ∼
2. In all casesthe central SB values are significantly reduced comparedto their intrinsic values due to dust extinction (also seeLaursen et al. 2009) and re-distribution of photons dueto radiative transfer effects. However, in the outer partsof the halo the intrinsic and post-processed SB valuesapproach each other.In order to be able to distinguish between the effect ofdust attenuation and systematic radial re-distribution ofphotons due to radiative transfer effects we also show inFig. 3 the ‘AGN’ case without any dust (green line). Ra-diative transfer reduces the SB values in the inner partof the halo, but increases them at larger radii. Thus v (km s − ) F l u x ( n o r m a li z e d ) Stars, m ion = 1 AGN, r ion , AGN = 20kpc
Fig. 4.— Ly α spectra taken from both our models in the directionused in Fig. 3. Note that the shapes of the spectra depend stronglyon our fiducial model parameters (see § scattering increases the effective size of Ly α halos (asalso found by Trebitsch et al. 2016, in their dust-freesimulations). Note though that these re-distributed pho-tons will also experience a relatively large dust opticaldepth, leading to lower SB values (as can be seen whencomparing the two ‘AGN’ curves in Fig. 3). When con-sidering, for example, the width of a halo with SB > − erg s − cm − arcsec − , these different effects wouldlead to uncertainties of a factor of ∼
2. In particular,hydrogen scattering tends to increase the apparent sizecompared to the intrinsic emission while dust absorptiondecreases it. This may be of importance when comparingto observations and will be discussed again in § Ly α spectra Fig. 4 shows spectra resulting from the full Ly α radia-tive transfer simulations. These spectra are taken fromthe same direction as the SB values from § α spectra shown in Fig. 4, as the applied resolution isnot sufficient to capture the full Ly α radiative transferdynamics (see, e.g., Verhamme et al. 2012, who showedthat ∼ pc resolution is necessary), and their shape heav-ily relies on the model parameters such as the minimallyallowed neutral fraction for ionized cells X HI , min . Thisvalue – which is zero for our fiducial models – stronglyinfluences the emergent spectral shape. We tested thisby increasing X HI , min to the rather extreme value of 10 − in the ‘AGN’ case and obtained wide (peak separation of ∼ − ) double peaked spectra instead. In addi-tion, we find the spectral shape to be dependent on theobserver’s direction. In particular, the flux at line centeris significantly lower (however, not so low as to form adouble-peaked profile) for other directions in which theoptical depth between the observer and the emitting re-gion is higher. This illustrates the difficulty using ab ini-tio hydro-dynamical simulations to predict the outcomeof Ly α spectra and their comparison with observations(see § L α,i (erg s − ) M a x i m a l e x t e n t ( k p c ) Cantalupo et al. (2014)Hennawi et al. (2015)Borisova et al. (2016) l og M / M fl Fig. 5.—
Maximal extent for halos in mass bins from 10 . M (cid:12) to 10 . M (cid:12) using a SB cutoff of 10 − erg s − cm − arcsec − versus intrinsic Ly α luminosity. The error bars denote thevariance between various directions and the red circle highlightsthe ‘fiducial’ halo shown in the previous figures. Also shownare the observations from Cantalupo et al. (2014), Hennawiet al. (2015) and Borisova et al. (2016) which implicitly assume f α, esc ∼
1. See § Comparison to observations
Cantalupo et al. (2014) observed a giant Ly α halowith luminosity of L α = (2 . ± . × erg s − ( L α ≈ . × erg s − including the central quasar)and maximal projected extent of 460 kpc ( SB > − erg s − cm − arcsec − ) . A similar – slightlysmaller – object was observed by Hennawi et al. (2015)with a total Ly α luminosity of L α ≈ . × erg s − and a maximum extent of 310 kpc. The ∼ . M (cid:12) halowe selected from the Illustris simulation shows verysimilar intrinsic Ly α luminosities ( L i ,α ∼ × erg s − )and extents even without introducing additional clump-ing below the simulation resolution. In our fiducial mod-els, the Ly α escape fraction was always ∼
25% leading toa slightly smaller Ly α luminosity than observed. How-ever, as the escape fraction is given by small-scale ra-diative transfer physics we see the obtained escape frac-tion as a lower limit to the real value. In particular twopoints are rather uncertain. First, the conversion be-tween metallicity and Ly α dust optical depth as givenby Eq. (2) is calibrated to local values and it is unclearhow reliable it is when applied to z ∼ α escape fraction.As already mentioned in § α halo, too. Inparticular, for our ‘fiducial’ halo presented previously wefound that while the intrinsic maximum extent (for which SB > − erg s − cm − arcsec − ) is ∼
250 kpc (forthe AGN and star ionization cases), the post-processedmaximum extent is ∼
318 kpc and ∼
415 kpc with andwithout dust, respectively . That Ly α scattering tends We acknowledge that it is difficult observationally to distin-guish between quasar and nebula emission. Note however the lu-minosity including the quasar can be seen as an upper limit to thenebular emission. These values are for the ‘AGN’ case ( r ion = 20 kpc) and themean of four directions orthogonal to each other. A smoothing to increase the observed size of an object is frequentlyinferred for other objects; galaxies, for instance, showsignificantly larger Ly α halos compared to their UVcounterparts (e.g., Hayes et al. 2014; Wisotzki et al.2016). Other theoretical work seems to confirm this pic-ture (e.g., Trebitsch et al. 2016, found that includingradiative-transfer effects enlarges the LAH in their sim-ulation by ∼ α ra-diative transfer persist as the HI and dust structure onthe smallest scales is unknown – which might alter ourresults. Maximal extents
In order to determine whether this halo is an ex-ceptional case – without running the full Ly α radia-tive transfer simulations on several halos – we used theintrinsic Ly α luminosity as given by Eq. (4) to com-pute ‘intrinsic surface brightness maps’ from which wecan measure the maximum projected extent. This ap-proach is supported by the findings presented in § M/M (cid:12) = 0 . . M (cid:12) and 10 . M (cid:12) which we ionized according to our ‘ionization from stars’model with m ion = 1. To mimic the observations,we then smooth the SB maps with a Gaussian kernel( F W HM = 1 arcsec) and then measure the maximal ex-tent for pixels with SB i > − erg s − cm − arcsec − .Fig. 5 shows the result of this analysis as the maximalextent versus the intrinsic luminosities of the 45 halos.The points and error bars represent the 16th, 50th and84th percentiles of the distribution of extents createdthough 100 randomly drawn observer’s directions, andthe color coding denotes the mass of the halo. Note thatthe size of the error bars thus shows a substantial size-variation with viewing angle. In addition, we displayedthe observations quoted above. These are not, however,the intrinsic values. L α, i for the observed nebula is prob-ably a factor of a few larger, and SB i should also be(slightly) increased. Nevertheless, Fig. 5 illustrates thatLy α halos with similar extent & luminosity can be foundin the Illustris simulation. Another interesting fea-ture is that the extent does not seem to correlate withmass and/or L α, i . We also checked the correlation of theextent versus various other halo properties (gas mass,gas metallicity, black hole mass) and found none of themto correlate significantly with the extent (all the Pear-son as well as the Spearman correlation coefficients were with F W HM = 1 arcsec has been applied. As shown in § m ion above ∼ . L α, i . Becausenot all (sub)halos possess black holes, however, we chose the ‘star’model out of simplicity. L α,i (erg s − ) A r e a ( k p c ) l og M / M fl Fig. 6.—
Projected area for halos in mass bins from 10 . M (cid:12) to 10 . M (cid:12) using a SB cutoff of 10 − erg s − cm − arcsec − versus intrinsic Ly α luminosity. The error bars denote the variancebetween 100 randomly drawn directions and the red circle marksthe ‘fiducial’ halo shown in Figs. 1-3. L α,i (erg s − ) R ( S B > − e r g s − c m − a r c s ec − ) ( k p c ) Cantalupo et al. (2014)Hennawi et al. (2015)Borisova et al. (2016), average l og M / M fl Fig. 7.— Ly α luminosity radius in which the SB (radially av-eraged) is > − erg s − cm − arcsec − versus Ly α luminosityfor halos in mass bins from 10 . M (cid:12) to 10 . M (cid:12) using a SBcutoff of 10 − erg s − cm − arcsec − . As previously, the errorbars denote the variance between 100 randomly drawn directionsand the red circle highlights our ‘fiducial halo’. Also shown arethe observations from Cantalupo et al. (2014) and Hennawi et al.(2015) (estimated from Fig. 12 of Battaia et al. 2016) as wellas the average profile of Borisova et al. (2016) which implicitlyassumes f α, esc ∼ § within [ − . , . Alternative size measurements
As an alternative to the ‘maximal extent’ measure, weplot in Fig. 6 and Fig. 7 the projected area covered bythe LAHs, and the radius for which the radially averagedSB profile falls below SB cut , respectively. Again we useda cutoff of SB cut = 10 − erg s − cm − arcsec − and aFWHM of the convolved Gaussian kernel of 1 arcsec. Itis noticeable that the correlation between area covered aswell as the ‘crossing radius’ and intrinsic Ly α luminosityis slightly better than between the maximal extent and L α, i shown in Fig. 5 (Spearman coefficients of ∼ . ∼ .
72 versus ∼ . -20 -19 -18 -17 SB cut (erg s − cm − arcsec − ) M a x i m a l e x t e n t ( n o r m a li z e d ) M a x i m a l e x t e n t i n k p c a t S B c u t = − e r g s − c m − a r c s e c − Fig. 8.—
Maximal extent versus the surface brightness cutoff SB cut . Each line represents one of 15 randomly selected halosfor which we plot the mean (using 10 random observer directions)normalized to the value at SB cut = 10 − erg s − cm − arcsec − (for 100 random observer directions). The color coding denotes thevalue of normalization. Borisova et al. (2016). The radial extent of the formertwo was estimated from Fig. 12 of Battaia et al. (2016),and thus implicitly assumes a Ly α fraction of unity – asalready discussed above for Fig. 5. Size dependence on SB cutoff
It initially appears puzzling that Ly α halos with thisextent have not been detected previously, as is nicelyillustrated in Fig. 3 of Cantalupo et al. (2014). How-ever, as noted in that work, one has to take into ac-count that previous halos were measured with a sur-face cutoff of SB cut = 3 × − erg s − cm − arcsec − whereas Cantalupo et al. (2014) used SB cut =10 − erg s − cm − arcsec − . That this difference mat-ters is shown in Fig. 8 where we plot how the maxi-mal extent varies with chosen surface brightness cutoff SB cut . In particular, we show the mean of the maximalextent using 10 random observer directions and normal-ized it to the mean of 100 random observer directions for SB cut = 10 − erg s − cm − arcsec − (as used in Fig. 5).Fig. 8 shows that the exact value of SB cut can alter themeasured projected extent dramatically and even a mod-erate change of factor three as discussed above can de-crease the size by ∼ SB cut to 10 − erg s − cm − arcsec − might more than doublethe size but even fainter cutoffs do not increase the extentfurther.While we were finalizing this work, two new studies ap-peared examining extended LAHs. Battaia et al. (2016)examined 15 quasars at z ∼ < z < >
100 kpc in every case examined (seealso Herenz et al. 2015, who found in 4 out of 5 cases noextended nebula). There thus seems at present some dis-agreement as to the abundance of these objects. Whetherthis is due to the differing samples and redshift ranges, or,as discussed in Borisova et al. (2016), differing observa-tional techniques, is beyond the scope of this work. How-ever, we note that our simulations predict fairly ubiqui-tous LAHs in the presence of even moderate ionizing fluxfrom the central black hole; in our simulated sample of For illustration purposes we show only 15 randomly selectedhalos. ¯ b (kpc) -20 -19 -18 -17 -16 -15 -14 S B i e r g s − c m − a r c s e c − Cantalupo+14Borisova+16 (average)Steidel+11 (LBG stack)Wisotzki+16 (LAEs) l og L α , i / ( e r g s − ) Fig. 9.—
Circularly averaged intrinsic surface brightness profilefor our sample of 60 nebulae between 10 . M (cid:12) and 10 . M (cid:12) .Each solid curve represents the average SB i profiles assembled from10 randomly drawn viewing angles and is color coded according tothe intrinsic Ly α luminosity. The dashed lines show the data fromCantalupo et al. (2014), Borisova et al. (2016), Steidel et al. (2011)and Wisotzki et al. (2016) (from top to bottom, extracted fromfigure 5 of Borisova et al. (2016)) as comparison. Note, however,that the data shows naturally the measured SB values and not theintrinsic ones (see § quasars, only ∼
15% had nebular emission below the ob-servable threshold of 10 − erg s − cm − arcsec − < Surface brightness profiles
In Fig. 9 we compare the full circularly averaged sur-face brightness profiles to the results of Cantalupo et al.(2014) and Borisova et al. (2016). Again, we used the‘stars’ model with m ion = 1 and a smoothing kernel withFWHM = 1 arcsec. Overall the agreement is very good,however, note that our profiles (solid lines) show the in-trinsic SB values wheres the observations (dashed lines)measure the SB after radiative transfer effects. Thismeans our curves represent an upper limit on the mea-sured SB levels due to the destruction of Ly α photonsby dust. In particular, comparing the curve of Can-talupo et al. (2014) to our findings suggest an overallhomogeneous escape fraction of (cid:38)
10% as the slopes arecomparable. The (averaged) SB profile of Borisova et al.(2016) is flatter in the inner region which could be dueto a higher dust attenuation the central region and re-distribution to larger radii due to scatterings. This is notunreasonable as the central region possesses a larger op-tical depth which makes photons originating there morevulnerable to destruction by dust and more likely to dif-fuse significantly in space before escape. We do in factfind this in the full radiative transfer simulations in § CONCLUSIONS
Several Ly α halos have recently been detected with un-usually large extent ( (cid:38)
300 kpc) and observed Ly α lumi-nosity ( (cid:38) erg s − ) (Cantalupo et al. 2014; Hennawiet al. 2015). It has been speculated that this suggests the presence of substantial amounts of cold gas in coolclumps at densities which are not resolved by modernhydro-dynamical simulations.Using the publicly available data of the Illustris sim-ulation (Vogelsberger et al. 2014; Nelson et al. 2015) wemodelled the Ly α emission emerging from large halos( M > . M (cid:12) ) at z ∼
2. We considered two simplemodels: an AGN as source of ionizing photons and ion-ization due to stars. We found for a single halo wherewe performed full Ly α radiative transfer that (i) witha moderate strong ionization source the halo showed in-trinsic Ly α luminosities and extents comparable to ob-servational data; (ii) additional clumping does not seemnecessary to explain the first order properties of the ob-served giant LAHs, and, (iii) due to the low optical depthof some routes escaping Ly α photons do not scatter manytimes leading to a single peaked Ly α profile in both mod-els. While the difference in the SB profiles cannot beused to distinguish between the models, the Ly α spec-tra are sensitive to sub-resolution properties such as thekinematics of the ISM and, thus, the question whetheror not the Ly α spectrum contains information about themain ionization source is still outstanding. In particular,changing the ionization state in the maximally ionizedregions slightly leads to a much greater optical depth atline center and, thus, an emergent double peaked as op-posed to the single peaked profiles in our fiducial models.This difficulty might be the cause between observed spec-tra of Ly α nebulae which are often double-peaked withan extended red- or blue-tail (e.g. Matsuda et al. 2006;Vanzella et al. 2016) and cosmological hydro-dynamicalsimulations which often predict a single peaked profile –as in this work (also see, e.g., Trebitsch et al. 2016).In order to put our fiducial halo into a wider context,we computed the intrinsic Ly α SB maps of 45 halos in themass range 11 . < log M/M (cid:12) < . α luminosities. Further-more, we find that that neither the maximal extent northe total projected area of Ly α halos correlates signifi-cantly with other halo properties such as the total (intrin-sic) Ly α luminosity which we attribute to the stochasticnature of the gas-morphology.Varying the surface brightness cutoff used for charac-terizing the extent of the Ly α halos, we find that multi-plying (dividing) this value by factor of three can increase(decrease) the maximal extent by ∼ α halos might not be that different to the smallerones detected previously but also that we expect to findmore of these objects in the near future.This conclusion – i.e., that ‘giant LAHs’ do appear inmodern hydro-dynamical simulations – differs from Can-talupo et al. (2014) who performed a similar analysis onone ∼ . M (cid:12) halo extracted from a hydro-dynamicalsimulation and found that substantial extra small-scalegas clumping was necessary to match their observations.We think this is for two reasons: (i) although their sim-ulation included supernovae feedback it did not includeAGN feedback which leads to a non-negligible distribu-tion of gas for halos in this mass range (see, e.g., vanDaalen et al. 2014; Genel et al. 2014, for an illustrationof this effect); (ii) as we show in Fig. 5 there is signif-icant scatter in Ly α extents – even for fixed luminosityand / or size. Thus, a single halo producing a Ly α SBmorphology as computed by Cantalupo et al. (2014) isentirely possible, even within Illustris.The authors thank the anonymous referee for theconstructive comments that significantly improved themanuscript. We thank M. Dijkstra and Ll. Mas-Ribas for a critical read of the draft. MG thanks the Physics& Astronomy department of Johns Hopkins Universityfor their kind hospitality. SB was supported by NASAthrough Einstein Postdoctoral Fellowship Award Num-ber PF5-160133. This research made use of Astropy, acommunity-developed core Python package for Astron-omy (Astropy Collaboration et al. 2013); matplotlib, aPython library for publication quality graphics (Hunter2007); SciPy (Jones et al. 2001–).
REFERENCESAo, Y., Matsuda, Y., Beelen, A., et al. 2015, Astronomy &Astrophysics, 581, A132Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al.2013, A&A, 558, A33Barnes, L. A., Garel, T., & Kacprzak, G. G. 2014, Publications ofthe Astronomical Society of the Pacific, 126, 969Battaia, F. A., Hennawi, J. F., Cantalupo, S., & Prochaska, J. X.2016, preprint, arXiv:1604.02942Borisova, E., Cantalupo, S., Lilly, S. J., et al. 2016, preprint,arXiv:1605.01422Cantalupo, S., Arrigoni-Battaia, F., Prochaska, J. X., Hennawi,J. F., & Madau, P. 2014, Nature, 506, 63Cantalupo, S., Porciani, C., & Lilly, S. J. 2008, The AstrophysicalJournal, 672, 48Cen, R., & Haiman, Z. 2000, ApJ, 542, L75Cen, R., & Zheng, Z. 2013, ApJ, 775, 112Dijkstra, M. 2014, Publications of the Astronomical Society ofAustralia, 31, 26Dijkstra, M., & Kramer, R. 2012, Monthly Notices of the RoyalAstronomical Society, 424, 1672Dijkstra, M., & Loeb, A. 2009, Monthly Notices of the RoyalAstronomical Society, 400, 1109Faucher-Gigu`ere, C.-A., Kereˇs, D., Dijkstra, M., Hernquist, L., &Zaldarriaga, M. 2010, The Astrophysical Journal, 725, 633Furlanetto, S. R., Schaye, J., Springel, V., & Hernquist, L. 2005,ApJ, 622, 7Fynbo, J. U., Møller, P., & Warren, S. J. 1999, MNRAS, 305, 849Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS,445, 175Goerdt, T., Dekel, A., Sternberg, A., et al. 2010, MNRAS, 407,613Gould, A., & Weinberg, D. H. 1996, ApJ, 468, 462Gronke, M., & Dijkstra, M. 2014, Monthly Notices of the RoyalAstronomical Society, 1103, 10Gronke, M., & Dijkstra, M. 2016, submitted to ApJHaiman, Z., & Rees, M. J. 2001, ApJ, 556, 87Haiman, Z., Spaans, M., & Quataert, E. 2000, ApJ, 537, L5Hayes, M. 2015, Publications of the Astronomical Society ofAustralia, 32, e027Hayes, M., ¨Ostlin, G., Duval, F., et al. 2014, ApJ, 782, 6Hennawi, J. F., & Prochaska, J. X. 2013, ApJ, 766, 58Hennawi, J. F., Prochaska, J. X., Cantalupo, S., &Arrigoni-Battaia, F. 2015, Science, 348, 779Herenz, E. C., Wisotzki, L., Roth, M., & Anders, F. 2015,Astronomy & Astrophysics, 576, A115Hui, L., & Gnedin, N. Y. 1997, Royal Astronomical Society, 292,27Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90Jones, E., Oliphant, T., Peterson, P., et al. 2001–, SciPy: Opensource scientific tools for Python, ,Laursen, P. 2010, PhD ThesisLaursen, P., & Sommer-Larsen, J. 2007, ApJ, 657, L69 Laursen, P., Sommer-Larsen, J., & Andersen, A. C. 2009, TheAstrophysical Journal, 704, 1640Martin, D. C., Matuszewski, M., Morrissey, P., et al. 2015,Nature, 524, 192Mas-Ribas, L., & Dijkstra, M. 2016, preprint, arXiv:1603.04840Matsuda, Y., Yamada, T., Hayashino, T., Yamauchi, R., &Nakamura, Y. 2006, ApJ, 640, L123Matsuda, Y., Yamada, T., Hayashino, T., et al. 2011, MonthlyNotices of the Royal Astronomical Society: Letters, 410, 13—. 2012, Monthly Notices of the Royal Astronomical Society,425, 878Matthee, J., Sobral, D., Oteo, I., et al. 2016, MNRAS, 458, 449Momose, R., Ouchi, M., Nakajima, K., et al. 2014, MonthlyNotices of the Royal Astronomical Society, 442, 110Nelson, D., Pillepich, A., Genel, S., et al. 2015, Astronomy andComputing, 13, 12Osterbrock, D. E. 1989, Astrophysics of gaseous nebulae andactive galactic nucleiPei, Y. C. 1992, Astrophysical Journal, 395, 130Price, D. J. 2012, Journal of Computational Physics, 231, 759Rahmati, A., Pawlik, A. H., Raiˇcevic, M., & Schaye, J. 2013a,Monthly Notices of the Royal Astronomical Society, 430, 2427Rahmati, A., Schaye, J., Pawlik, A. H., & Raiˇcevic, M. 2013b,Monthly Notices of the Royal Astronomical Society, 431, 2261Rauch, M., Haehnelt, M., Bunker, A., et al. 2008, TheAstrophysical Journal, 681, 856Reuland, M., van Breugel, W., R¨ottgering, H., et al. 2003, ApJ,592, 755Rosdahl, J., & Blaizot, J. 2011, preprint, 000, arXiv:1112.4408Springel, V. 2010, MNRAS, 401, 791Steidel, C. C., Adelberger, K. L., Shapley, A. E., et al. 2000, ApJ,532, 170Steidel, C. C., Bogosavljevi´c, M., Shapley, A. E., et al. 2011, TheAstrophysical Journal, 736, 160Steidel, C. C., Erb, D. K., Shapley, A. E., et al. 2010, TheAstrophysical Journal, 717, 289Trainor, R. F., & Steidel, C. C. 2012, ApJ, 752, 39Trebitsch, M., Verhamme, A., Blaizot, J., & Rosdahl, J. 2016,preprint, arXiv:1604.02066van Daalen, M. P., Schaye, J., McCarthy, I. G., Booth, C. M., &Dalla Vecchia, C. 2014, MNRAS, 440, 2997Vanzella, E., Balestra, I., Gronke, M., et al. 2016, preprintVerhamme, A., Dubois, Y., Blaizot, J., et al. 2012, Astronomy &Astrophysics, 546, A111Villar-Mart´ın, M., S´anchez, S. F., Humphrey, A., et al. 2007,MNRAS, 378, 416Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Nature, 509,177Wisotzki, L., Bacon, R., Blaizot, J., et al. 2016, Astronomy &Astrophysics, 587, A98Yusef-Zadeh, F., Morris, M., & White, R. L. 1984, ApJ, 278, 186Zheng, Z., & Miralda-´Escude, J. 2002, The Astrophysical Journal,578, 33APPENDIX
INTRINSIC SB MAPS
Fig. 10 shows the intrinsic (i.e. without radiative transfer effects and (dust) extinction) surface brightness maps.These surface brightness maps have been assembled as described in § SB = 10 − erg s − cm − arcsec − levelsfrom of the full radiative transfer simulation ( § l o g S B ( e r g s − c m − a r c s e c − ) kpc k p c Stars, m ion = 1 kpc AGN, r ion , AGN = 20 kpc kpc
AGN, r ion , AGN = 20 kpc (no dust)
Fig. 10.—
Intrinsic surface brightness maps for the fiducial models (see § SB = 10 − erg s − cm − arcsec − contours from of the full radiative transfer simulation (see § One can note that compared to Fig. 2 the intrinsic SB maps are firstly brighter in the central region due to theeffect of dust extinction. This effect has been discussed in Sec. 3.2 (see Fig. 3). Secondly, the scattering enlarges andwashes out the SB contours. However, the extent for SB (cid:38) − erg s − cm − arcsec −2