Modelling the cosmological Lyman–Werner background radiation field in the early Universe
Abstract
The Lyman–Werner (LW) radiation field is a key ingredient in the chemothermal evolution of gas in the early Universe, as it dissociates H2 molecules, the primary cooling channel in an environment devoid of metals and dust. Despite its important role, it is still not implemented in cosmological simulations on a regular basis, in contrast to the ionizing UV background. This is in part due to uncertainty in the source modelling, their spectra and abundance, as well as the detailed physics involved in the propagation of the photons and their interactions with the molecules. The goal of this work is to produce an accurate model of the LW radiation field at z ≥ 6, by post-processing the physics-rich high-resolution FiBY simulation. Our novelties include updated cross-sections for H2, H− and H+2 chemical species, IGM absorption by neutral Hydrogen and various spectral models for Population III and Population II stars. With our fiducial set of parameters, we show that the mean LW intensity steadily increases by three orders of magnitude from z ∼ 23 to z ∼ 6, while spatial inhomogeneities originate from massive star-forming galaxies that dominate the photon budget up to a distance of ∼100 proper kpc. Our model can be easily applied to other simulations or semi-analytical models as an external radiation field that regulates the formation of stars and massive black hole seeds in high-z low-mass haloes.
Full text
MNRAS 522, 330–349 (2023) https://doi.org/10.1093/mnras/stad1008 Advance Access publication 2023 April 6 Modelling the cosmological Lyman–Werner background radiation field in the early Uni v erse Andrea Incatasciato , 1 ‹Sadegh Khochfar 1 and Jose O ˜ norbe 2 1 Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK 2 Facultad de F ´ ısicas, Universidad de Sevilla, Avda. Reina Mercedes s/n, Campus Reina Mercedes, E-41012 Seville, Spain Accepted 2023 March 31. Received 2023 March 31; in original form 2023 January 23 A B S T R A C T The Lyman–Werner (LW) radiation field is a key ingredient in the chemothermal evolution of gas in the early Universe, as it dissociates H 2 molecules, the primary cooling channel in an environment devoid of metals and dust. Despite its important role, it is still not implemented in cosmological simulations on a regular basis, in contrast to the ionizing UV background. This is in part due to uncertainty in the source modelling, their spectra and abundance, as well as the detailed physics involved in the propagation of the photons and their interactions with the molecules. The goal of this work is to produce an accurate model of the LW radiation field at z ≥6, by post-processing the physics-rich high-resolution FiBY simulation. Our no v elties include updated cross-sections for H 2 , H −and H + 2 chemical species, IGM absorption by neutral Hydrogen and various spectral models for Population III and Population II stars. With our fiducial set of parameters, we show that the mean LW intensity steadily increases by three orders of magnitude from z ∼23 to z ∼6, while spatial inhomogeneities originate from massive star-forming galaxies that dominate the photon budget up to a distance of ∼100 proper kpc. Our model can be easily applied to other simulations or semi-analytical models as an external radiation field that regulates the formation of stars and massive black hole seeds in highz low-mass haloes. Key words: astrochemistry –molecular processes –radiative transfer –methods: numerical –stars: Population III – early Universe. 1 INTRODUCTION Molecular Hydrogen (H 2 ) is a key ingredient of the early-Universe chemistry, as it represents the main cooling channel of pristine gas at T < 10 4 K (Saslaw & Zipoy 1967 ; Peebles & Dicke 1968 ). Light primordial elements such as Hydrogen and Helium are efficient coolants in their atomic form only abo v e that temperature. On the other hand, heavier elements (collectively referred to as metals ) do not form during the Big Bang Nucleosynthesis and are a product of the evolution and explosion of stars (Kobayashi, Karakas & Lugaro 2020 ), either in isolation or in binary systems; hence cooling due to metal-line transitions (Smith, Sigurdsson & Abel 2008 ), C-, F-, and O-based molecules and dust grains (Hirashita & Ferrara 2002 ) starts dominating the energy budget of the interstellar medium (ISM) only after the first chemical enrichment episodes. The abundance of H 2 (and secondarily of other simple molecules, e.g. HD and HeH + ) strongly influences the thermodynamical evolution of the gas that condenses in the first mini-haloes forming at redshift z ≤30 (see e.g. Abel, Bryan & Norman 2000 ; Galli & Palla 2013 , for a re vie w). Molecular cooling allo ws the gas to reach temperatures as low as ∼200 K, condense to high densities and form the first Population III (PopIII) stars (Haiman, Thoul & Loeb 1996 ; Tegmark et al. 1997 ). Analytical models, 1D and 3D simulations all show that the compressional heating that develops E-mail: [email protected] while gas falls into dark matter haloes is efficiently dissipated with a central H 2 fractional abundance of at least 10 −5 –10 −4 (Abel et al. 2000 ; Machacek, Bryan & Abel 2001 ; Yoshida et al. 2006 ; Latif & Khochfar 2019 ). This sets a clear consensus about the initial phase of metal-free PopIII star formation episodes, while models diverge on the final outcome of this process (the multiplicity and the Initial Mass Function –IMF - of PopIII stars), due to differences in the spatial and mass resolution, and in the treatment of accretion, gas chemistry, and turbulence (see Bromm & Larson 2004 , for a re vie w, or e.g. Hirano et al. 2015 ; Chiaki & Yoshida 2022 ; Latif, Whalen & Khochfar 2022b , for more recent discussions). Nevertheless, PopIII stars are generally thought to be massive and hot (Bromm, Coppi & Larson 1999 ; Abel, Bryan & Norman 2002 ) and are predicted to emit a copious amount of energetic photons during their very short lifetime (Schaerer 2002 ). The y e xplode as violent supernovae, leaving black hole remnants with masses ∼10– 100 M (Fryer, Woosley & Heger 2001 ; Madau & Rees 2001 ) and enriching the Universe with metals (Heger & Woosley 2002 ), that pave the way for the formation of the first proto-galaxies made of metal-poor Population II (PopII) stars (Bromm & Loeb 2003 ). Due to their peculiar features, PopIII stars represent also the most important source of Lyman–Werner (LW) photons at the Cosmic Dawn (e.g. Haiman, Abel & Rees 2000 ; Agarwal et al. 2012 ). The LW radiation lies within the soft-UV part of the electromagnetic spectrum (its range is commonly indicated as 11–13.6 eV, or 911– 1150 Å) and is able to efficiently dissociate H 2 through the twostep Solomon process (Solomon 1965 ; Stecher & Williams 1967 ). ©The Author(s) 2023. Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( http://cr eativecommons.or g/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/mnras/article/522/1/330/7110413 by Biblioteca Universidad de Sevilla user on 20 May 2024
Modelling the Lyman–Werner background 331 MNRAS 522, 330–349 (2023) H 2 formation can be also prevented with the detachment of H − and the dissociation of H + 2 , due to NIR-VIS-NUV photons with a few to ∼10 eV (Glo v er 2015a , b ). H −and H + 2 indeed represent the two main H 2 formation channels in the ISM at moderate densities and devoid of dust grains. Radiation abo v e the Lyman limit, that in principle would be able to directly dissociate H 2 molecules, is rapidly absorbed by atomic H in the diffuse intergalactic medium (IGM), that is still completely neutral at this stage. LW photons, on the contrary, have a very long mean free path ( ∼100 cMpc; Haiman et al. 2000 ; Ahn et al. 2009 ), as they can only be absorbed when redshifted to the exact frequencies of the atomic Lyman transitions. H 2 molecules, instead, are not dense enough in the IGM to play any role in this context. This leads to the definition of a spatially nearly homogeneous background, whose intensity at the Lyman limit at z ∼25 −10 is often bracketed by J 21 ∼10 −3 and J 21 ∼10 2 , where J 21 is the Lyman–Werner background (LWB) intensity normalized to 10 −21 erg s −1 sr −1 Hz −1 cm −2 (e.g. Haiman, Rees & Loeb 1997 ; Machacek et al. 2001 ; Ahn et al. 2009 ; Trenti & Stiavelli 2009 ; Johnson, Dalla Vecchia & Khochfar 2013 ). The build-up of a homogeneous LWB during the formation of the first cosmological structures has important implications on the PopIII star formation (Haiman et al. 2000 ), as it makes molecular cooling inefficient in low-mass haloes. Without H 2 molecules, star formation is delayed until dark matter haloes reach virial temperatures of T vir ∼10 4 K, when atomic H cooling becomes efficient and the collapse can start (Haiman et al. 1997 ). Recently, many theoretical efforts have been focused on trying to quantify this effect. The interplay between a time-varying LWB and PopIII star formation has been explored with cosmological hydrodynamical simulations (Wise et al. 2012b ; Johnson et al. 2013 ), that are designed to accurately capture the highly non-linear evolution of cosmic structures, and with semianalytical/seminumerical models (Haiman et al. 2000 ; Ahn et al. 2009 ; Trenti & Stiavelli 2009 ; Agarwal et al. 2012 ; Qin et al. 2020 ; Visbal, Bryan & Haiman 2020 ), that on the other hand require a certain number of approximations and a priori assumptions, but allow a fast parameter exploration. In addition, Latif & Khochfar ( 2019 ), Kulkarni, Visbal & Bryan ( 2021 ), Schauer et al. ( 2021 ) employed high-resolution small-scale cosmological simulations to explore the minimum halo mass required for PopIII star formation in molecular cooling haloes under a range of constant LWB intensities. The modelling of the LW radiation is usually approximated, due to the technical complexity of the calculation from first principles (Abel et al. 1997 ; Wolcott-Green, Haiman & Bryan 2017 ) and the lack of constraints on the spectra of the stellar populations responsible for the LW emission (Bromm & Larson 2004 ). Often, only very young PopIII and PopII stars are considered in the radiative budget, stellar evolution is neglected, and the emissivity is kept constant (Greif & Bromm 2006 ). A fully self-consistent treatment of the closed loop between star formation and the growth of a LWB, that e x erts a ne gativ e feedback on the subsequent star formation episodes, is also made difficult by the computational cost of radiative transfer methods o v er large cosmological volumes (Johnson et al. 2013 ). Another matter of debate is the importance of the LW radiation in the context of the Direct Collapse Black Hole (DCBH) scenario (Begelman, Volonteri & Rees 2006 ; Lodato & Natarajan 2006 ; Dijkstra et al. 2008 ; Agarwal et al. 2012 ), that represents one of the most promising formation channels of the initial seeds of the supermassive black holes observed at z > 6 (Fan et al. 2006 ; Mortlock et al. 2011 ; Ba ˜ nados et al. 2018 ). Haloes illuminated by high LW intensity, such as small star-less satellites of massive highredshift galaxies, where the radiation from the neighbouring galaxies pre v ails by orders of magnitude o v er the large-scale background, have been proposed as birthplaces of black holes with initial masses of 10 4 −6 M (Agarwal et al. 2014 , 2019 ; Wise et al. 2019 ; Lupi, Haiman & Volonteri 2021 , see also Bonoli, Mayer & Callegari 2014 ; Fernandez et al. 2014 ). A critical value J 21, crit of LW intensity is usually assumed to express the minimum level of radiation needed to efficiently prevent H 2 molecular cooling, the first-order requirement of the DCBH scenario together with a pristine chemical composition. In the last few years many studies have explored its feasibility with 1D and 3D hydrodynamical simulations that employ non-equilibrium chemistry, high spatial and temporal resolution, and in some cases a selfconsistent treatment of the radiative feedback from the central object (Omukai, Schneider & Haiman 2008 ; Shang, Bryan & Haiman 2010 ; Regan, Johansson & Haehnelt 2014 ; Luo et al. 2018 ). Ho we ver, a consensus on the value of J 21, crit is still lacking. Recent studies (e.g. Sugimura, Omukai & Inoue 2014 ; Latif et al. 2015 ; WolcottGreen et al. 2017 ) have shown that, if the interstellar radiation field is modelled as a blackbody, J 21, crit can vary by many orders of magnitude (from 10 to 10 5 ), depending on the assumed blackbody temperature, usually 10 5 K (10 4 K) if PopIII (PopII) stars dominate the radiation field. Agarwal & Khochfar ( 2015 ) have highlighted that considering the evolution of the spectral shape across the lifetime of a stellar population has an important impact on J 21, crit , especially when also long-li ved lo w-mass stars are included. Furthermore (Glo v er 2015a , b ; Agarwal et al. 2016 ; Sugimura et al. 2016 ; Luo et al. 2020 ) all proposed that the chemical netw ork emplo yed in the simulations should also include H −detachment and H + 2 dissociation, to provide a more accurate estimate of the H 2 formation rate. Further degrees of freedom include the H 2 self-shielding treatment, that in the optically thick regime can reduce the effect of the LW photons by up to three orders of magnitude and strongly depends on the accuracy of the calculation (Draine & Bertoldi 1996 ; Wolcott-Green, Haiman & Bryan 2011 ; Hartwig et al. 2015b ; Wolcott-Green & Haiman 2019 ), and the impact of additional fields, such as X-rays or cosmic rays, that can increase the fraction of free electrons, thus facilitating the formation of H 2 (Inayoshi & Omukai 2011 ; Inayoshi & Tanaka 2015 ; Glo v er 2016 ; Re gan, Johansson & Wise 2016 ; Park, Ricotti & Sugimura 2021 ). These uncertainties lie on top of other aspects, such as unresolved fragmentation and long-term sustainability of the mass accretion flow, whose role the scientific community still has to have a final say on (e.g. Ge & Wise 2017 ; Bhowmick et al. 2022 ). With this work we tackle some of the current limitations of the studies on the effect of the LW radiation on the formation of stars and black hole seeds in the early Universe. In particular, we aim at sho wing ho w the LWB and the other associated photochemical rates can be accurately modelled given a star formation history, that can be either derived from a simulation or a semi-analytical model. We also study the spatial inhomogeneities of the radiation field (Haiman et al. 2000 ; Ahn et al. 2009 ; Dijkstra, Ferrara & Mesinger 2014 ). To do so we post-process the First Billion Year (FiBY) suite of cosmological simulations. We describe the FiBY project in Section 2.1 , together with the methods employed in the post-processing algorithm. We keep an agnostic approach with regards to the IMF of PopIII and PopII stars (Section 2.2 ), in order to show the intrinsic uncertainty due to the current lack of constraints on the stellar models. We outline our code that accurately calculates the photochemical rates in Section 2.3 , while deferring to a companion paper (Incatasciato et al., in preparation) for a in depth comparison of the specific methods to calculate the H 2 dissociation rate given a stellar spectrum. IGM absorption is described in Section 2.4 . Our homogeneous LWB model is outlined in Section 3.1 , while its spatial inhomogeneities Downloaded from https://academic.oup.com/mnras/article/522/1/330/7110413 by Biblioteca Universidad de Sevilla user on 20 May 2024
332 A. Incatasciato, S. Khochfar and J. O ˜ norbe MNRAS 522, 330–349 (2023) are investigated in Section 3.2 . Finally, our considerations on the impact of the LWB on PopIII star formation are reported in Section 3.3 . We then complete the paper with further discussions and our conclusions in Section 4 . 2 METHODS In this work, we use the simulations of the First Billion Year (FiBY) project to e v aluate the e volution of the LW background 1 at z ≥6. The FiBY suite is described in Section 2.1 . To obtain an estimate of the LW background we sample the radiation field by choosing random points ( observers ) within the simulation box. For each observer we sum the radiation emitted by all the sources, taking into account various stellar models (described in Section 2.2 ), a detailed calculation of the photochemical rates, including recent updates to take into account molecular level populations (Section 2.3 ) and the absorption by the IGM (Section 2.4 ). The number of points used in each snapshot is selected such as the mean and the standard deviation of J 21 converge to per cent level and corresponds to O(10 4 ). We repeat this e x ercise for each snapshot available at z ≥6. 2.1 FiBY The FiBY project (Johnson et al. 2013 ; Paardek ooper, Khochf ar & Dalla 2013 ; Agarwal et al. 2014 ; Paardek ooper, Khochf ar & Dalla Vecchia 2015 ; Cullen et al. 2017 ; Phipps et al. 2020 ) includes a set of high-resolution and physics-rich cosmological simulations of the early Universe. The simulations were run with a modified version of the GADGET-3 code (Springel, Yoshida & White 2001 ; Springel 2005 ), already employed for the OWLS project (Schaye et al. 2010 ). The code was updated further to include the rele v ant physical processes and stellar models required for a better modelling of the formation of the first stars and proto-galaxies at z ∼30 −6. Substructures within the simulations are identified with the SUBFIND algorithm (Springel et al. 2001 ) and merger trees are generated with the method described in Neistein et al. ( 2012 ). We refer the reader to the original FiBY papers and the references therein for a detailed description of the subgrid models and provide here only a brief summary. For completeness and convenience of the reader, we also summarize the key parameters of all the runs (e.g. box size and mass resolution) used in this work in Table 1 . All the simulations were run using the following cosmological parameters, consistent with those reported by the Wilkinson Microwave Anisotr opy Pr obe ( WMAP ) team in Komatsu et al. ( 2009 ): m = 0.265, b = 0.0448, = 0.735, H 0 = 71 km s −1 Mpc −1 , and σ8 = 0.81. The same cosmological parameters are assumed throughout this work, unless otherwise stated. Collisionless dark matter particles and SPH gas particles are the two main constituents of the simulated volumes. The thermodynamical evolution of the gas particles is described with the usual atomic cooling due to H and He, but also with metal line cooling (C, N, O, Ne, Mg, Si, S, Ca, and Fe; Wiersma, Schaye & Smith 2009 ) and H 2 and HD non-equilibrium chemistry (Abel et al. 1997 ; Yoshida et al. 2006 ). The multiphase interstellar medium (ISM) is modelled with an ef fecti ve equation of state (EOS), following Schaye & Dalla Vecchia ( 2008 ), explicitly designed to yield star formation rates consistent with the observed Schmidt–Kennicutt law (Schmidt 1959 ; Kennicutt 1 Here and in the following, when we refer to the LW background, we implicitly consider not only the photons responsible for the H 2 dissociation, but also the ones rele v ant for the H −detachment and the H + 2 dissociation. 1998 ). The density threshold for the ef fecti ve EOS is 10 cm 3 , that represents also the threshold for the star formation. Depending on the metallicity of the star forming gas, collisionless particles representing metal-free Population III or metal-poor Population II stars are spawned. Each stellar particle is treated as a single stellar population with a perfectly sampled IMF. PopIII stars are assigned a Salpeter ( 1955 ) IMF with stellar masses in the range 21–500 M , consistent with the top-heavy IMF predicted e.g. by Bromm & Larson ( 2004 ), while PopII stars have a Chabrier ( 2003 ) IMF that extends down to subsolar masses. The critical metallicity to distinguish PopIII and PopII stellar particles is Z crit = 10 −4 Z (Maio et al. 2011 ) 2 , where Z = 0.02. Massive stars explode as supernovae at the end of their lives. Metal enrichment due to core-collapse (CCSNe, 8 M < M < 100 M ) and pair-instability (PISNe; 140 M < M < 260 M ) supernov ae follo ws He ger & Woosle y ( 2002 , 2010 ). The thermal energy due to the explosions is stochastically injected to the neighbouring particles following the scheme described by Dalla Vecchia & Schaye ( 2012 ). The cosmic reionization is modelled with a time-dependent and spatially uniform UV radiation background (Haardt & Madau 2001 ), while high-density gas is shielded against the UVB as proposed by Nagamine, Choi & Yajima ( 2010 ). One specific run (FiBY LW, see Table 1 ) includes also an on-the-fly LW background, that comprises both a homogeneous component dependent on the cosmic star formation rate (Greif & Bromm 2006 ) and the contribution from the local sources (Johnson et al. 2013 ). In this run, the H 2 self-shielding follows the prescriptions by WolcottGreen et al. ( 2011 ). The main sources of LW photons in the early Universe are PopIII and PopII stars. The simulations predict high-redshift UVluminosity function and a star formation main sequence that are in good agreement with observational constraints (Cullen et al. 2017 ), as well as an o v erall star formation rate density ( ρSFR ) that is in fair agreement with observational bounds (Johnson et al. 2013 ). This gives us confidence that to first order stars form within the simulations at the right rate and in the right objects. For the sake of completeness we show in Fig. 1 ρSFR of the XL (solid purple), L (dashed blue), M (dotted green), and S (dash-dotted red) FiBY simulations, superimposed o v er the one obtained with deep HST observations at z ∼4–10 from two collaborations, in grey (Oesch et al. 2014 ; Bouwens et al. 2016 ; Oesch et al. 2018 ) and in red (McLure et al. 2013 ; Bowler et al. 2015 ; McLeod et al. 2015 ; McLeod, McLure & Dunlop 2016 ), and recent constraints from ground (COSMOS/UltraVISTA) and JWST NIRCam photometry (Donnan et al. 2023 ) at z ∼8–15 in blue. We make use of the shaded regions to highlight the uncertainties due to the underlying assumption on the stellar IMF, where the higher (lower) values are for a Salpeter 1955 (Chabrier 2003 ) IMF, that have different UV luminosity per stellar mass conversion factors (Madau & Dickinson 2014 ). To ensure a consistency between observations and simulations, we include only galaxies with SFR SFR min = 0.3 M yr −1 , that corresponds to the integration limit of the UV luminosity function down to M UV = −17 as in Oesch et al. ( 2018 ) and Donnan et al. ( 2023 ). The simulations employ the same number of particles to evolve the dark matter and baryonic density fields (684 3 each) within cubic volumes of different sizes (from 32 to 4 cMpc). They therefore 2 The specific choice of the critical metallicity does not strongly impact our results, as metal pollution proceeds quickly and increases the metallicity of the interstellar medium (ISM) to large values in the hosting and neighbouring haloes (Maio et al. 2010 , 2011 ; Smith et al. 2015 ). Downloaded from https://academic.oup.com/mnras/article/522/1/330/7110413 by Biblioteca Universidad de Sevilla user on 20 May 2024
Modelling the Lyman–Werner background 333 MNRAS 522, 330–349 (2023) Table 1. Compilation of the different FiBY simulations used in this work. We report the box size in the second column, where the corresponding value in brackets is the size when assuming h = 0.71. Columns 3, 4, and 5 show the number of particles of each component (dark matter and gas) and their (initial) mass. The final redshift reached is listed in Column 6, while the next one indicates whether the LW background is self-consistently calculated starting from the star formation rate and the local contribution of young stars (see Section 2.2 of Johnson et al. 2013 ). The last column reports a rough estimate of the mass of the smallest haloes resolved (with at least 50 dark matter particles). Name Boxlength (cMpc h −1 ) (cMpc) N 1/3 M dm (M h −1 ) M g (M h −1 ) z end LW background M halo,min (M ) FiBY XL 22.72 (32) 684 2.24 ×10 6 4.56 ×10 5 4 N 1.6 ×10 8 FiBY L 11.36 (16) 684 2.80 ×10 5 5.70 ×10 4 4 N 2 ×10 7 FiBY M 5.68 (8) 684 3.50 ×10 4 7.12 ×10 3 6 N 2.5 ×10 6 FiBY S 2.84 (4) 684 4.37 ×10 3 8.90 ×10 2 6 N 3 ×10 5 FiBY LW 2.84 (4) 684 4.37 ×10 3 8.90 ×10 2 6 Y 3 ×10 5 Figure 1. The star formation rate density ( ρSFR ) in the following FiBY simulations: XL (solid purple), L (dashed blue), M (dotted green), and S (dash-dotted red line). Coloured shaded areas indicate the ρSFR derived from deep HST and JWST + COSMOS observ ations: Bo wler et al. ( 2015 , z = 5 −7), McLure et al. ( 2013 , z = 8), McLeod et al. ( 2015 , z = 9), and McLeod et al. ( 2016 , z = 10) are in red, Bouwens et al. ( 2016 , z = 4 −8), Oesch et al. ( 2014 , z = 9) and Oesch et al. ( 2018 , z = 10) in grey and Donnan et al. ( 2023 ) in blue. To ensure a consistency between observations and simulations, we estimate the FiBY ρSFR only from galaxies with SFR 0.3 M yr −1 , corresponding to the usual integration limit of M UV = −17 in the observed UV luminosity function. The horizontal lines indicate the minimum value that can be predicted by FiBY due to the limited volume of each simulation. The shaded regions quantify the uncertainties due to the assumptions on the stellar IMF: the higher (lower) values are for a Salpeter 1955 (Chabrier 2003 ) IMF, that give a slightly different UV luminosity per stellar mass, with a correction factor of 0.63 as suggested by Madau & Dickinson ( 2014 ). investigate different sections of the halo mass function. In particular, only the M and S boxes properly resolve H 2 -cooling haloes with M h < 10 7 −8 M (last column in Table 1 ), but lack the rarer massive galaxies due to the limited volume. On the other hand, the larger L and XL boxes focus on atomic-cooling haloes and include a wide range of cosmic environments, hence they contain massive galaxies above the observational limit of SFR min from earlier times ( z ∼11–14, exactly the redshift range that is currently being studied for the first time with JWST ; Donnan et al. 2023 ; Harikane et al. 2023 ). The minimum ρSFR that can be estimated from each simulation (corresponding to only one galaxy with SFR = SFR min in the entire volume) is shown with the horizontal thin lines. The ρSFR in FiBY shows a reasonable convergence between the XL, L, and M volumes, and an evolution with redshift that is consistent with observations by Oesch et al. ( 2018 ) and collaborators. Ho we ver, the slightly different absolute values suggest that FiBY might o v erproduce stars in massiv e galaxies at z ࣠ 10. Nev ertheless, deviations on this level cannot be too surprising, given the uncertainties both on the observational side (dust correction, incompleteness) and in simulations (resolution, LW radiation, stellar feedback, only to name a few, see e.g. Vogelsberger et al. 2020 ). The shallower evolution found by McLeod et al. ( 2016 ) and Donnan et al. ( 2023 ) indeed demonstrates that observations still do not provide a unique solution to this problem. Finally, the S box shows a peculiar evolution at z = 6–9 that matches McLeod et al. ( 2016 ), despite being very close to the limits set by the small simulated volume. In conclusion, the FiBY simulations produce a fairly realistic highz Universe and we consider it a useful tool to model the evolution of the LWB in the pre-Reionization Era. 2.2 Stellar emission We use nine different models for the spectral energy distribution (SED) of the stars. For PopIII and PopII stars we employ the models described in Table 2 and Table 3 , respectively, calculated with the publicly available stellar population synthesis (SPS) codes Yggdrasil (Zackrisson et al. 2011 ), Slug2 (da Silva, Fumagalli & Krumholz 2012 , 2014 ), and BPASS (Stanway & Eldridge 2018 ). Yggdrasil uses models for PopIII stars from Schaerer ( 2002 ) and Raiter, Schaerer & Fosbury ( 2010 ) and provides pre-computed SEDs for a very topheavy IMF (Salpeter 1955 between 50 and 500 M ) and a more moderate one (lognormal with characteristic mass equal to 10 M ). Slug2, instead, allows the user to calculate stellar SEDs with a wide variety of IMFs and evolutionary tracks (e.g. P ado va Bressan et al. 1993 ; Gene v a Eggenberger et al. 2008 ; and MIST Dotter 2016 ), and atmosphere models resembling the Starburst99 SPS code of Leitherer et al. ( 1999 ). 3 BPASS provides a large set of pre-computed SEDs with an in-depth treatment of stellar binary systems. The minimum available metallicity is 5 ×10 −4 Z , hence we use BPASS models only for metal-poor PopII stars. To include the approximation of the stellar spectra commonly assumed in the literature, we additionally consider two blackbody (BB) spectra with T rad = 10 5 K and 10 4 K for PopIII and PopII stars, respectively. The normalization of these spectra is chosen such that the number of emitted photons in the LW range per stellar baryon ηLW is 2 ×10 4 and 4 ×10 3 , respectively, as adopted in Greif & Bromm ( 2006 ) and Johnson et al. ( 2013 ). When using the BB spectra we do not consider stars older than 5 Myr, in order to match the model used in the FiBY (Johnson et al. 2013 ). PopIII and PopII stars often coexist in simulated highz galaxies. The total LWB is therefore calculated as the sum of the contributions 3 Another important feature of Slug2 is to allow a stochastic sampling of the IMF . W e do not mak e use of it in this w ork. Downloaded from https://academic.oup.com/mnras/article/522/1/330/7110413 by Biblioteca Universidad de Sevilla user on 20 May 2024
334 A. Incatasciato, S. Khochfar and J. O ˜ norbe MNRAS 522, 330–349 (2023) Table 2. SEDs for PopIII stars, where the masses are in M . Nebular emission and extinction are neglected. The Slug2 spectra are calculated with the MISTv1.0 (Choi et al. 2016 ; Dotter 2016 ) non-rotating stellar tracks and are generated with logarithmic time-steps of 0.05 dex from 1 Myr to 1 Gyr, in order to accurately resolve the rapid evolution of young stellar populations. The presented metallicity is the minimum available for these stellar tracks. Name IMF Parameters Metallicity (Z ) SPS code References PopIII Ygg1 Salpeter ( 1955 ) M min = 50, M max = 500 0 Yggdrasil Schaerer ( 2002 ); Zackrisson et al. ( 2011 ) PopIII Ygg2 Lognormal M min = 1, M max = 500 0 Yggdrasil Raiter et al. ( 2010 ); Zackrisson et al. ( 2011 ) M c = 10, σ= 1 PopIII Slug Salpeter ( 1955 ) M min = 21, M max = 300 10 −4 Slug2 da Silva et al. ( 2012 , 2014 ) Table 3. Same as Table 2 , but for for PopII stars. α1 and α2 are the low-mass and the high-mass slopes, respectively, and the masses are in M . The BPASS SEDs include binaries as according to Stanway & Eldridge ( 2018 ) and have M min = 0.1 M . Nebular emission and extinction are neglected. Name IMF Parameters Metallicity (Z ) SPS code References PopII BPASS TH Double power law α1 = −1.3, α2 = −2 5 ×10 −4 BPASS Stanway & Eldridge ( 2018 ) M t = 0.5, M max = 300 PopII BPASS Chab Chabrier ( 2003 ) M t = 1, M max = 100 5 ×10 −4 BPASS Stanway & Eldridge ( 2018 ) PopII BPASS BH Double power law α1 = −1.3, α2 = −2.7 5 ×10 −4 BPASS Stanway & Eldridge ( 2018 ) M t = 0.5, M max = 100 PopII Slug Chabrier ( 2003 ) M t = 1, M max = 120 10 −3 Slug2 da Silva et al. ( 2012 , 2014 ) from these two distinct stellar populations, where we consider 5 out of the 20 possible combinations: (i) FID : PopIII Ygg2 + PopII BPASS Chab , this is our ’fiducial’ choice; see bottom panel of Fig. A1 for an example of the SEDs 1 Myr after the star formation episode; (ii) TH : PopIII Ygg1 + PopII BPASS TH , with top-heavy IMFs; (iii) BH : PopIII Ygg2 + PopII BPASS BH , with bottom-heavy IMFs; (iv) SLUG : PopIII Slug + PopII Slug , where both SEDs are calculated with the Slug2 SPS code; (v) BB : PopIII BB5 + PopII BB4 , with single-temperature blackbody spectra. In particular ‘TH’ and ‘BH’ should bracket the level of uncertainty introduced by the choice of IMF, where the contribution from highmass and low-mass stars, respectively, is enhanced with respect to our fiducial model and to all the other combinations neglected in this work. We show in Fig. 2 the emission rate per stellar baryon of LW photons in the range 11–13.6 eV for each SED used in this work. Black and red lines indicate PopIII and PopII models, respectively. Our fiducial choice for PopIII stars (black thick solid line) is conserv ati ve, as it could have been predicted from the IMF, since the other PopIII models emit ∼5 times more LW photons in the early stages, but die off very quickly after 5–10 Myr. The differences in the IMFs for PopII stars can be noticed in the first 10 Myr, where the number of high-mass stars determines a factor of 4–5 higher (lower) LW emission of BPASS TH (BPASS BH) with respect to the fiducial case (red thin solid, dashed, and thick solid lines, respectively), while they all show pretty similar evolution at later times. PopII BB is hardly distinguishable from the fiducial SED, but it shows a totally different shape when it comes to low energy photons that determine the H + 2 dissociation and H −detachment rates (see bottom panel in Fig. A1 ). Once the radiation is emitted by a star, it takes a not-negligible time to reach an observer at a giv en como ving distance r com . By Figure 2. LW photon emission rate per stellar baryon for each SED used in this work. Models for PopIII and PopII stars are shown in black and red, respectively: with reference to the definitions in the main text and Tables 2 –3 , the ’FID’ (thick solid), the ’TH’ (thin solid), the ’SLUG’ (dot-dashed), the ’BB’ (dashed), and finally the BPASS bottom-heavy IMF (red dashed line). reverting equation (8) of Ahn et al. ( 2009 ), in the high-redshift limit z em = (1 + z obs ) −1 / 2 −H 0 1 / 2 m r com 2 c −2 −1 (1) expresses the redshift of emission of a photon observed at z obs (here the redshift of a given simulation snapshot). We hence account for the light time travel t = t ( z obs ) −t ( z em ) when choosing the age of the emitted spectrum of each stellar particle, where t ( z) is the age of the Universe at a given redshift. The time resolution of the spectra generated with the SPS codes listed in Tables 2 –3 allows us to properly follow the spectral evolution of a stellar population. The physical motion of the stars is instead negligible in this context. 2.3 Optically thin photochemical rates Three photoreactions need to be taken into account to accurately e v aluate the formation and destruction of H 2 in pristine gas (e.g. Glo v er 2015a ): H 2 + γ→ H ∗ 2 → H + H hν 6 . 7 eV (2a) Downloaded from https://academic.oup.com/mnras/article/522/1/330/7110413 by Biblioteca Universidad de Sevilla user on 20 May 2024
Modelling the Lyman–Werner background 335 MNRAS 522, 330–349 (2023) H −+ γ→ H + e −hν 0 . 76 eV (2b) H + 2 + γ→ H + H + hν 0 . 1 eV (2c) Reaction ( 2a ) represents the indirect dissociation of H 2 molecules by LW photons via the two-step Solomon process (Solomon 1965 ), while reactions ( 2b ) and ( 2c ) are the detachment of H −and the dissociation of H + 2 , respectively; these two chemical species are the main catalysts that lead to the formation of H 2 at moderate densities in a gas devoid of metals and dust. The full frequency-dependent computation of these rates requires the following integration over the rele v ant range of photon energies: k [s −1 ] = νmax νmin 4 πJ νσ( ν) hν d ν, (3) where J νis the radiation intensity in erg s −1 sr −1 Hz −1 cm −2 , σνis the frequency-dependent cross section in cm 2 , and νmax is the Lyman limit that corresponds to the ionization energy of Hydrogen atoms at 13.6 eV. As commonly assumed, photons abo v e this threshold are ne glected, as the y are quickly absorbed by the ISM in the proximity of the source. We stress here that for the dissociation of molecules in Reactions ( 2a ) and ( 2c ) the minimum energies required are lower than the threshold energies usually adopted in the literature ( ∼11 and ∼2.65 eV, respectively, see e.g. Abel et al. 1997 ; Glo v er 2015a ). The latter are valid when only the ground state roto-vibrational level of the respective molecule is taken into account; ho we ver, we include in our model also the appropriate population of e xcited lev els, that have lower bounding energies, hence the lower threshold energies. A detailed discussion on the molecular level populations is deferred to a companion paper (Incatasciato et al., in preparation), while here we limit the description of the rates calculation to a more general level. 2.3.1 H 2 photodissociation rate The indirect photodissociation of H 2 molecules takes place through the Solomon process (Solomon 1965 ): a molecule in the rotovibrational level with quantum numbers ( v, J ) of the electronic ground state X 1 + g is excited to a ( v , J ) state of the B 1 + u or C 1 + / − u electronic level, due to the absorption of a Lyman or Werner photon, respectiv ely. A fraction ( ∼15 per cent , on av erage) of the excited molecules then decay into the vibrational continuum of the ground state, resulting in its subsequent dissociation (Abgrall, Roueff & Drira 2000 ). To compute the optically thin dissociation rate, we follow the approach described in Draine & Bertoldi ( 1996 ), Abel et al. ( 1997 ), Wolcott-Green et al. ( 2011 ), for the continuum limit: the ’ef fecti ve cross section’ can be calculated with σ( ν) = C v,J ⎡ ⎣ v ,J i ∈ LW f osc , i V( ν−ν0 , i ) f diss ,v ,J N X ( v, J ) ⎤ ⎦ , (4) where C = πe 2 4 πm e c0 (Corney 1977 ) 4 and the summation runs over all the possible LW transitions, excited, and ground state levels. V( ν− ν0, i ) is the Voigt line profile of the i-th transition between the ground state level ( v, J ) and the excited state level ( v , J ), whose width takes into account both the natural damping coefficient and the thermal 4 e is the electron charge, m e is the electron mass, c is the speed of light in vacuum, and 0 is the electric constant. broadening; f osc,i is the transition oscillator strength, f diss ,v ,J is the fraction of molecules that dissociate after the excitation, and N X ( v, J ) is the fraction of molecules initially in the level ( v, J ). The LW transitions are taken from the data bases of Ubachs et al. ( 2019 ) and Salumbides et al. ( 2015 ), where the transition frequency ν0 , the oscillator strength f osc , and the natural damping coefficient are reported for each transition. These two data sets are complementary and are updated versions of the widely used data base by Abgrall et al. ( 1993a , b , c ). The fraction f diss ,v ,J of excited molecules ( v , J ) that dissociate is instead derived from Abgrall et al. ( 2000 ) as A c / A t , where A c is the probability of decay to the vibrational continuum and A t is the total probability of decay of an electronically excited state B/C. We also include the data from Abgrall et al. ( 1997 ) 5 to derive the mean kinetic energy of the products of the dissociation (two H atoms), that in turn allows us to estimate the average heating rate due to H 2 photodissociation. We find ∼0.4 eV per dissociated molecule, similar to Black & Dalgarno 1977 , but with some variations of the order of 30 per cent depending on the gas temperature and density and the shape of the incident spectrum. We also take into account how the roto-vibrational levels of the electronic ground state X 1 + g are populated, for a given combination of gas temperature and density. These levels can be excited and deexcited both due to collisions or the absorption/emission of photons. When these processes are frequent enough (abo v e a certain density threshold) the local thermodynamical equilibrium (LTE) is reached and the level population follows the Boltzmann distribution. The LTE density threshold for H 2 molecules lies between 10 3 and 10 6 cm −3 depending on the gas temperature. At lower densities only the first rotational levels of the ground vibrational level ( v = 0, J = 0–3) are populated. At intermediate densities, we interpolate between the non-LTE ( k H 2 , 0 ) and the LTE case ( k H 2 , LTE ) as in Glo v er ( 2015a ): k H 2 = k H 2 , LTE k H 2 , 0 k H 2 , LTE α , (5) where α= (1 + n / n crit ( T gas )) −1 . Interested readers will find a detailed description of n crit in Incatasciato et al. (in preparation). In this work, we do not vary the gas temperature (set at 10 3 K) and density (set at 10 2 cm −3 ), that is well within the non-LTE limit, n crit being approximately 3 orders of magnitude higher at 10 3 K. This choice of gas temperature and density ensures that the H 2 and H + 2 dissociation rates are representative of the initial stages of collapse of a gas cloud in a low-metallicity environment (Omukai et al. 2005 ). The effect of the LW radiation during the subsequent evolution at densities n 10 4 cm −3 would instead involve other physical processes, such as the H 2 self-shielding (Wolcott-Green et al. 2011 ; Hartwig et al. 2015b ; Wolcott-Green & Haiman 2019 ), that are beyond the scope of this work. 2.3.2 H + 2 photodissociation Updated state-resolved cross sections for the H + 2 photodissociation and the inverse process (radiative association) are available in the literature (e.g. Babb 2015 ; Zammit et al. 2017 ) for photons with energies as high as 40 eV. We choose to use the data from Zammit et al. ( 2017 , 2018 ) as they are including the cross sections for all the 423 roto-vibrational levels of the electronic ground state of H + 2 . This makes the calculation more reliable in the LTE limit at high temperatures ( T gas ∼10 3 –10 4 K). They also take into account transitions to 23 different electronic excited states, while Babb ( 2015 ) 5 H 2 continuum emission probabilities at https://molat.obspm.fr Downloaded from https://academic.oup.com/mnras/article/522/1/330/7110413 by Biblioteca Universidad de Sevilla user on 20 May 2024
336 A. Incatasciato, S. Khochfar and J. O ˜ norbe MNRAS 522, 330–349 (2023) considers only the first excited state 2p σu : this is less crucial for the purpose of our work, as the total cross section is essentially due to transitions to the continuum of the first excited state, other than at energies 12 eV ( λ< 1000 Å), where the contribution from the other states is noticeable, and we neglect photons above the Hydrogen ionization limit. We again follow the approach of Glo v er ( 2015a ) as in equation ( 5 ) to interpolate between the non-LTE and the LTE rates. For the non-LTE limit, we assume that all the molecules are in the rotovibrational level with the lowest energy ( v = 0, J = 0) (see e.g. Shapiro & Kang 1987 ; Glo v er 2015a ; Latif et al. 2015 ). The LTE limit is assumed for a gas density abo v e the critical value n crit , that is determined as in Glo v er ( 2015a ) Section B1.2, assuming that H atoms and free electrons are the most important collisional partners. 2.3.3 H −photodetachment Several cross sections are available in the literature for this process: Shapiro & Kang ( 1987 ); John ( 1988 ); Chuzhoy, Kuhlen & Shapiro ( 2007 ); McLaughlin et al. ( 2017 ). We choose the latter, that for the first time includes the resonances at 11 eV. This gives an increase of ∼20 per cent on the H −detachment rate for energetic spectra (e.g. a blackbody spectrum with T rad = 10 5 K; Glo v er 2015b ). 2.4 IGM optical depth Haiman et al. ( 2000 ) show that the LW radiation can be efficiently absorbed by the diffuse neutral gas of the IGM. In the pre-reionization Universe, in fact, the Hydrogen optical depth in the Lyman lines (energy range 10.2–13.6 eV) is very high ( τ∼10 6 at z ∼20); hence, LW photons are absorbed by H atoms as soon as they are cosmologically redshifted into the closest atomic Lyman transition. As suggested by Haiman et al. ( 2000 ) and Ahn et al. ( 2009 ), the contribution of H 2 molecules to the IGM optical depth is subdominant and can be neglected due to its low abundance in the diffuse gas. If we consider a LW photon emitted at redshift z em with energy h νem and its closest Lyman line with energy h νline , the maximum distance at which the photon can be observed corresponds to a minimum redshift z obs expressed as 1 + z obs 1 + z em = νline νem . (6) Under the assumption of a homogeneous LW background, this formula leads to the definition of the ’sawtooth modulation’ (see fig. 1 in Haiman et al. 2000 ; Ahn et al. 2009 , for a flat emitted spectrum). Ho we ver, this works aims at studying the LW background beyond the homogeneous-universe approximation. We need to consider the so-called ’picket-fence’ modulation factor by Ahn et al. ( 2009 ), that describes (from the point of view of a single source) the fraction of unabsorbed spectrum in the LW energy range at a comoving distance r com : f mod = Max (0 , A exp [ −( r com /Bα) C ] −D) , (7) where A = 1.7, B = 116.29, C = 0.68, D = 0.7, r com is in (comoving) Mpc and α= h 0 . 7 −1 m 0 . 27 −1 / 2 1 + z em 21 −1 / 2 (8) contains the dependence on the cosmological parameters and the redshift of emission. The parameters in the expression were estimated by Ahn et al. ( 2009 ) considering a flat spectrum in the energy interval Figure 3. H 2 modulation factor assuming a PopIII Ygg2 stellar population at different ages. The black dot-dashed line is the Ahn et al. ( 2009 ) fit, that reproduces fairly well f mod for young stars ( < 10 Myr, in blue), but fails to match it for increasing ages (green, yellow, and red solid lines). 11.5–13.6 eV. From equation ( 7 ) the maximum distance that a LW photon can travel is R LW ≃ 97 αcMpc. This defines a ’LW horizon’, that represents the largest volume that needs to be considered in order to e v aluate a self-consistent LW background. Fialkov et al. ( 2013 ) showed that a more accurate e v aluation of f mod is obtained if the full frequency-dependent calculation is performed without simplifying assumptions, such as the LW transitions being uniformly distributed in frequency or a flat incident spectrum. We build further on this, by recalculating the fit of equation ( 7 ) for each SED mentioned in Section 2.2 and for each of the three photochemical rates described in Section 2.3 . In particular, we calculate here for the first time the modulation factor for k H −and k H + 2 , that have a much larger ‘horizon’, as the corresponding threshold energies are lower than than E Ly α= 10.2 eV. Given a stellar spectrum and the assumption of a Universe at the mean density, we perform the full frequency-dependent rate calculation for the transmitted spectrum at distances from 0.1 αto 500 αcMpc. The result is then fitted using the same functional form as in equation ( 7 ), where we fix D = A −1, but we re-e v aluate the other parameters for each SED at each stellar age. By automatically incorporating the appropriate cross section, in this work (as also in Fialkov et al. 2013 ) f mod represents the true correction factor to the photochemical rates 6 and not just the fraction of unabsorbed spectrum as in Ahn et al. ( 2009 ). Fig. 3 shows how f mod for the H 2 dissociation rate changes if spectra with different shapes are assumed. In particular, we use PopIII Ygg2 SEDs at different ages, from young energetic spectra (in blue, < 10 Myr) to intermediate and old stellar populations (green, yellow, and red solid lines). The dot-dashed black line represents the fit from Ahn et al. ( 2009 ), that matches reasonably well only the modulation factor for young stellar populations at r com < 40 αcMpc. At larger distances the fit o v erestimates it, ev en though less than in Fialko v et al. ( 2013 ), as in our larger data set of LW transitions some at ∼11 eV do not end up in any absorption window until very large 6 We neglect here the resonant photons that are absorbed by the neutral IGM and re-emitted at lower frequencies. The data set of LW transitions adopted in this work is more extended than the one in Fialkov et al. ( 2013 ) and some low-energy photons might still excite H 2 molecules through LW transitions at energies < 10.2 eV. Ho we ver, our choice is still a reasonable assumption, given that those transitions are not important outside the LTE regime of dense gas (defined as abo v e a critical density of n ∼10 4 cm −3 ). Downloaded from https://academic.oup.com/mnras/article/522/1/330/7110413 by Biblioteca Universidad de Sevilla user on 20 May 2024
Modelling the Lyman–Werner background 337 MNRAS 522, 330–349 (2023) Figure 4. Modulation factors (solid lines) and the corresponding fits (dashed lines) for the three reactions and a 1-Myr-old PopIII Ygg2 SED. The optimal parameters valid for this particular stellar model are listed in Table 4 . As in Fig. 3 the black dot-dashed line is the Ahn et al. ( 2009 ) fit for H 2 . It can be noted here that f mod for H 2 is slightly non-monotonic, unlike what is proposed by Ahn et al. ( 2009 ) and Fialkov et al. ( 2013 ), due to the presence of LW transitions at energies lower than the Lymanαtransition, that do not enter any absorption windo w. This, ho we ver, represents only a second-order effect, while the trend first shown by Ahn et al. ( 2009 ) is confirmed and is valid for young stellar populations in general (see the text for the relative discussion). For the same reason f mod for H 2 does not reach exactly zero at the LW horizon, but values as low as 10 −3 –10 −4 that can be approximated with zero. Table 4. Optimal parameters for the fit shown in Fig. 4 . The parameters are defined as in equation ( 7 ) and we have fixed D = A −1. Rate A B C H 2 1.548 79.733 0.719 H −0.357 217.445 0.776 H + 2 1.188 220.907 0.831 H 2 - Ahn et al. ( 2009 ) 1.7 116.29 0.68 distances, hence they contribute to the dissociation rate. For older stellar populations, instead, f mod e volves more and more slo wly with the distance, as hard photons at ∼13 eV, that are absorbed closer to the emitting source due to the high density of Lyman lines, have a minor impact to the total H 2 dissociation rate. We hav e v erified that the same trend is found for all the other SEDs we use in this work, for both PopIII and PopII stars. Assuming Ahn et al. ( 2009 ) fitting function regardless of the stellar age would hence lead to an underestimation of the contribution of old stars to the LW background by a factor of 2–3. In Fig. 4 , we show f mod for the H 2 dissociation, the H −detachment and the H + 2 dissociation (blue, orange, and green solid lines, respectively). The dashed lines show the relative fits, with the black dotdashed line being the Ahn et al. ( 2009 ) fit for H 2 . The SED used in this illustrativ e e xample is the one for a 1-Myr-old PopIII Ygg2 stellar population and the corresponding fitting parameters are reported in Table 4 . As already discussed abo v e, the Ahn et al. ( 2009 ) fit closely describes the H 2 f mod for the radiation emitted by young stars (our fitting parameters are only slightly different) and the LW horizon at ∼100 cMpc is still valid and independent from the spectral shape (see Fig. 3 ). f mod for the other two reactions, instead, decreases much more slowly with the distance and never actually reaches zero, as photons below the Lyman αline are not absorbed by the neutral Hydrogen. This in principle would imply that the volume employed for the calculation of the background is not finite. Ho we ver, we choose to Figure 5. Top panel: Cumulative contribution to the LW background of the sources inside a sphere of radius r (normalized to R LW ), for the XL simulation at z = 8.9 and for the ’FID’ choice of SEDs, with (solid lines) and without (dotted lines) considering the IGM absorption. The three rates, H 2 dissociation, H −detachment, and H + 2 dissociation, are colour coded as in Fig. 4 . As expected, including the IGM optical depth has a stronger impact on the H 2 dissociation rate, while the H −detachment rate is almost unaffected. Bottom panel. Cumulative fractional contribution to the LW background, inside a sphere of radius r as in the top panel. Here we account for the IGM absorption. By definition the H 2 dissociation rate converges within R LW (grey dashed vertical line). We consider a maximum distance of ∼5 R LW (black dashed vertical line) to have a convergence of k H + 2 and k H −too. limit it to a sphere of radius 500 αcMpc. We moti v ate our strategy in the next section. 2.4.1 Cosmological volume needed for the LWB evaluation In the top panel of Fig. 5 we report the cumulative contribution to the LW background of spheres of increasing radius and centred on random points within the simulation volume. We use here the XL box at z = 8.9 and the ‘FID’ choice of SEDs. The solid and dotted lines indicate whether the IGM optical depth is included (solid) or not (dotted). The largest sphere has a radius of ∼500 αcMpc, the same maximum distance considered for the e v aluation of the modulation factor (Fig. 4 ). In order to reach a volume that is larger than the simulated box, we stack several copies of the box until the target sphere is reached. As expected, including the IGM optical depth has the strongest impact on the H 2 dissociation rate, that is reduced by a factor of 5–10. The H + 2 dissociation rate is only moderately reduced, while the H −detachment rate is almost unaffected. In the bottom panel, we demonstrate instead the convergence of the three rates, in terms of cumulative fractional contribution to the LW background. By definition, the H 2 dissociation rate converges within R LW (gre y dashed v ertical line). H −detachment and H + 2 dissociation instead converge somewhere between 1 and 5 R LW , and the exact distance slightly depends on the star formation history and the stellar models. We hence assume that all three photochemical rates converge within a maximum distance of 5 R LW . This sets the volume that needs to be considered around a giv en observ er in order to determine the LW background in that point. Only the total H −detachment rate could be underestimated by not more than 5 per cent at low redshift, due to f mod being significantly greater than zero at any distance. Downloaded from https://academic.oup.com/mnras/article/522/1/330/7110413 by Biblioteca Universidad de Sevilla user on 20 May 2024
338 A. Incatasciato, S. Khochfar and J. O ˜ norbe MNRAS 522, 330–349 (2023) Figure 6. The evolution with redshift of the mean H 2 dissociation rate, given by the LW radiation measured in the XL (solid purple), L (dashed blue), M (dotted green), and S (dash-dotted red) FiBY simulations. The colour and line-style scheme is consistent with Fig. 1 . On the secondary y-axis we express the LW radiation intensity at 13.6 eV in units of J 21 = 10 −21 erg s −1 Hz −1 sr −1 cm −2 , where we use the approximate relation k H 2 = 1 . 38 ×10 −12 J 21 s −1 commonly used in the literature. The simulation boxes described in Table 1 (with sizes ranging from 4 to 32 cMpc) are smaller than 5 R LW . As already done in Ahn et al. ( 2009 ), we account for this by attaching multiple copies of the simulation box next to the central one, until the maximum distance is reached. This ensures that we include all the sources that contribute to the radiation background as measured in the central box. The drawback of our method is that the conclusions we can draw on the inhomogeneities of the radiation field are certainly limited, as the simulated volume is not able to capture the total variance of the cosmic structures that we would expect to find in a sphere with radius ∼500 αcMpc. 3 RESULTS 3.1 Mean LW background In Fig. 6 , we show the LW background obtained with the postprocessing method described in Section 2 , using the ’FID’ choice for the stellar models. The four coloured lines represent the mean H 2 dissociation rate measured in the XL (solid purple), L (dashed blue), M (dotted green), and S (dash-dotted red) FiBY simulations. The secondary y-axis expresses the common conversion between the H 2 dissociation rate and the LW radiation intensity at the Lyman limit (13.6 eV) J 21 , in units of 10 −21 erg s −1 Hz −1 sr −1 cm −2 . With this definition, and under the assumption of low-density gas in the optically thin limit and a flat incident spectrum, the H 2 dissociation rate is k H 2 = 1 . 38 ×10 −12 J 21 s −1 (Abel et al. 1997 ) 7 The LWB intensity generally increases with time, as primordial low-mass haloes and then proto-galaxies grow in mass and trigger the formation of more and more stars. The LWB reaches mean values well abo v e J 21 ∼1 at z ࣠ 10 in all simulations. We attribute the small difference in the XL run to the slightly lo wer ρSFR in massi ve galaxies (Fig. 1 ), that in turn can be explained with a systematic degradation of the star formation efficiency in lower resolution simulations, as already found in IllustrisTNG (Pillepich et al. 2018 ). The S box shows a more irregular evolution at low redshift due to low number 7 Note that this scaling is generally only valid for young stellar populations (Shang et al. 2010 ; Agarwal & Khochfar 2015 ; Glo v er 2015b ). Figure 7. Top panel: the evolution of the three photochemical rates considered in this work, in the M simulation and for our ’FID’ choice of SEDs. Bottom panel: the ratio of H −detachment rate (orange) and H + 2 dissociation rate (green) to H 2 dissociation rate, for the FiBY simulations. of high-mass haloes that dominate the photon budget in the small volume. On the other hand, at z > 12 the different halo mass range that is resolved by each simulation determines when the LWB starts to build up and its intensity. Both M and S simulations resolve star formation in low-mass haloes and hence show a good convergence from early times, even if there is a hint of a missing contribution from M h ࣠ 10 6 M haloes (not resolved in M, see the last column in Table 1 ) at z ≥20. L and XL, instead, have delayed PopIII star formation, as the y resolv e only haloes with M h 10 7 M and M h 10 8 M , respectively. This is reflected into a delayed build-up of the LW intensity, that is 5–10 (100) times lower in L (XL) than in M at z ∼20–15. Overall, we find that the LWB from the FiBY simulations, when the rele v ant halo mass range is resolved, is well fitted by the following polynomial, with 6 < z < 23: log J 21 = A + B(1 + z) + C(1 + z) 2 (9) with A = 2.119, B = −1.117 ×10 −1 , and C = −2.782 ×10 −3 . 3.1.1 Effective LW spectral shape As already stated abo v e, in this work we include also the H − detachment and H + 2 dissociation rate, that are important to determine the rate at which H 2 molecules form during the initial phases of gas collapse. The top panel of Fig. 7 illustrates the evolution of these rates in the M box, for our ‘FID’ choice of SEDs. The blue line is the LWB previously shown in Fig. 6 , while the orange and the green lines are k H −and k H + 2 , respectively. The rates concurrently grow in time, as the UV photons emitted by young massive stars are the major contributors to all of them; nevertheless, an increasing additional contribution of IR photons is present in k H −and k H + 2 . The bottom panel of the same figure presents the ratio between the latter two rates and k H 2 , for the same four FiBY runs as in Fig. 6 . At the zero-th order (see e.g. Latif et al. 2015 ), the H −detachment rate (in orange) and the H + 2 dissociation rate (in green) are approximately two and and one order(s) of magnitude higher than the H 2 dissociation Downloaded from https://academic.oup.com/mnras/article/522/1/330/7110413 by Biblioteca Universidad de Sevilla user on 20 May 2024
Modelling the Lyman–Werner background 345 MNRAS 522, 330–349 (2023) consistent with ours at z ࣠ 10, and Haardt & Madau ( 2012 )-derived LWB only increases by a factor of 2 between z = 15 and z = 6. (iii) For the first time we show the mean H −detachment and H + 2 dissociation rates (Fig. 7 ), necessary to properly model H 2 formation (Glo v er 2015a ; Sugimura et al. 2016 ). Based on how these rates evolve with z with respect to the H 2 dissociation rate (equation 10 ), we find that the resulting spectral shape of the LWB can be approximated with a blackbody spectrum with an ef fecti ve temperature that evolves from 6 ×10 4 K at z ∼23 to 2 ×10 4 K at z ∼6. The same analysis performed on the rates derived from the UV backgrounds of Haardt & Madau ( 2012 ) and Faucher-Gigu ` ere et al. ( 2009 ) gives a consistent evolution for the former (from 3 ×10 4 K at z ∼15 to 1.5 ×10 4 K at z ∼6), and a very hard spectrum with T eff ∼10 5 K for the latter, that would require massive PopIII stars to dominate the UV radiation field even at z ∼10. (iv) The high-resolution FiBY simulations suggest that the contribution from PopIII stars is dominant at z > 12 (Fig. 8 ). Ho we ver, it is worth noting that the exact transition time somewhat depends on the resolution of the simulations and the associated metal enrichment: lower mass and spatial resolution delays the transition from a PopIIIto a PopII-dominated star formation to as late as z ∼10 (Maio et al. 2010 ). (v) Young stellar populations undoubtedly provide the largest contribution to the LWB (Fig. 9 ), thanks to the presence of shortlived hot massive stars (Eldridge & Stanway 2022 , and references therein). None the less, since we follow the evolution of the stellar spectra during their entire lifetime (up to 1 Gyr, more than the age of the Universe at z = 6), we are able to determine also the contribution from old stars. We find that stars older than 20 Myr account for 20 per cent −30 per cent of the H 2 dissociation rate at z ࣠ 10, and up to 60 per cent in the H −detachment rate. This shows that the LW radiation intensity is often underestimated in simulations in the literature. FiBY LW, for example, considers only the LW radiation emitted by stars younger than 5 Myr (Johnson et al. 2013 ): ho we ver, those stellar populations never contribute to more than 60 per cent of the total LW intensity and can reach as low as 20 per cent of the H −detachment rate at z ࣠ 10. The Renaissance suite of simulations (O’Shea et al. 2015 ), as well as its progenitor (Wise et al. 2012b ) and descendant (Phoenix, Wells & Norman 2022 ), adopts a similar approach, where PopII stars contribute to the UV radiation field only if younger than 20 Myr, with a constant luminosity equal to their lifetime-averaged luminosity. In their case, ho we ver, PopIII stellar e volution is followed, under the assumption of a delta function as IMF, centred on 40 M , 100 M , and 20 M , respectively . Similarly , Ahn et al. ( 2009 ) in their semi-analytical model consider a constant stellar emissivity that approximates the spectrum of a young ( < 20 Myr) stellar population. (vi) Our fiducial choice for stellar models and IMFs reflects the fact that PopIII stars have generally a higher characteristic mass and hotter atmospheres compared to solar-metallicity models with standard Salpeter ( 1955 ) or Chabrier ( 2003 ) IMFs. Nevertheless, the ongoing debate on the IMF of metal-free and metal-poor stars (e.g. Abel et al. 2002 ; Frebel, Johnson & Bromm 2007 ; Hirano et al. 2015 ; Stacy, Bromm & Lee 2016 ; Rossi, Salvadori & Sk ´ ulad ´ ottir 2021 ) leads us to relax the initial hypothesis and consider multiple sets of stellar models (see Tables 2 –3 ). The mean LW intensity is increased (decreased) by a factor of 2–3 with more top-heavy (bottom-heavy) IMFs (Fig. B1 ); interestingly, a blackbody with T = 10 4 K, as often assumed in the literature (Shang et al. 2010 ; Johnson et al. 2013 ; Glo v er 2015a ), giv es a factor of 400 and 30 higher H −detachment and H + 2 dissociation rate, respectively, due to the extremely different spectral shape that does not resemble the common SED models of young PopII stellar populations (Fig. A1 ). Present and future observations will be key to complement highresolution simulations of PopIII star formation (Hirano et al. 2015 ; Stac y et al. 2016 ; P ark et al. 2021 ) and specific models for PopIII stars (Schaerer 2002 ; Raiter et al. 2010 ; Gessey-Jones et al. 2022 ; Larkin, Gerasimov & Burgasser 2023 ), to put tighter constraints on their IMF. JWST and ALMA will investigate spectral signatures of PopIIIdominated galaxies (Yajima & Khochfar 2017 ; Woods et al. 2021 ; Latif et al. 2022b ; Nakajima & Maiolino 2022 ), while precise metal abundances measurements in Damped Lymanαsystems (Welsh, Cooke & Fumagalli 2019 ; Welsh et al. 2022 ) model the chemical enrichment from PopIII supernovae and stellar archaeology is already constraining the IMF low-mass end with local observations of extremely metal-poor stars (Frebel et al. 2007 ; Hartwig et al. 2015a ; Rossi et al. 2021 ; Hartwig et al. 2022 ). This will help reducing the uncertainties in the LWB modelling presented here. (vii) equation ( 9 ) models the LWB determined with the methods presented in this work and can be safely used as a realistic homogeneous background in simulations that do not resolve low-mass haloes or do not have enough volume to reach the LW horizon. Ho we ver, it is tightly connected with the modelling of star formation in FiBY. Explicitly taking into account the amount of young and old stars allows us to evaluate on-the-fly a LWB that is more general and can be applied to any simulation. For this reason we include equation ( 11 ) and the fitting parameters in Table 7 , that if applied to the stellar densities in FiBY reconstruct the LWB presented in this work (Fig. 10 ). (viii) We also study the spatial fluctuations of the LW radiation field: the H 2 dissociation rate presents a right-skewed distribution (Fig. 11 ), with a minimum that is very close to the mean value and a long tail extending 1–2 orders of magnitude abo v e the mean. Such tail progressively shrinks at z ࣠ 12, while the H −detachment and H + 2 dissociation rates show lower fluctuations due to the IGM being more transparent in their corresponding energy range (Fig. 12 ). (ix) The highest peaks in the LWB, usually in close proximity to massive star-forming galaxies, have been proposed as the birthplaces of massive black hole seeds at z 10 (see e.g. Dijkstra et al. 2008 ; Agarwal et al. 2014 , 2019 ; Wise et al. 2019 ; Lupi et al. 2021 ). We model the critical distance at which a single galaxy dominates o v er the homogeneous background (Fig. 13 ), that depends on the stellar mass and can be easily rescaled for a LWB intensity that in general can vary with z (Fig. C1 ). Our recommendation to simulators is to go beyond the time-varying homogeneous LWB and to include spatial fluctuations due to local sources. Equations ( 12 )–( 13 ) can be easily incorporated in a simulation, limiting at the same time the computational domain where the UV radiation needs to be accounted for on-the-fly with computationally e xpensiv e radiativ e transfer methods. (x) We use the homogeneous background found in this work to estimate the minimum halo mass for PopIII star formation under the influence of the LW radiation, with the use of a number of analytical and numerical studies available in the literature (Machacek et al. 2001 ; Latif & Khochfar 2019 ; Kulkarni et al. 2021 ; Lupi et al. 2021 ; Schauer et al. 2021 ). Using the LWB in equation ( 9 ), we obtain a minimum mass that approximately evolves from ∼3 ×10 5 M at z = 20 to ∼10 7 M at z ࣠ 10 (Fig. 14 ). The numerical experiments cited here unfortunately suffer from important limitations, such as: (i) H −detachment and H + 2 dissociation are often neglected, and (ii) the LW intensity is assumed as constant throughout simulations that run from z = 30 to z = 15, a wide time window during which the LWB can grow by ∼3.5 orders of magnitude (if equation ( 9 ) is extrapolated to z = 30). We then show in Fig. 15 our preliminary results on the effect of LW radiation in delaying PopIII star formation Downloaded from https://academic.oup.com/mnras/article/522/1/330/7110413 by Biblioteca Universidad de Sevilla user on 20 May 2024
346 A. Incatasciato, S. Khochfar and J. O ˜ norbe MNRAS 522, 330–349 (2023) in molecular-cooling haloes. In particular, when J 21 ∼0.05 at z = 20 the minimum mass of star-forming haloes in FiBY LW starts diverging from the case where the LW radiation is neglected. By z = 13 ( J 21 ∼0.2), only atomic-cooling haloes abo v e T vir = 10 4 K can form stars. Such results should be treated with caution, though, as further analysis is required to distinguish inefficient H 2 -cooling due to LW radiation from other hydrodynamical processes that might have a similar effect, such as dynamical heating (Fernandez et al. 2014 ) or pressure support by turbulence (Latif et al. 2022a ). We highlight here three important caveats of this work. First, the calculation of the LWB is performed in post-processing. Compared to FiBY S, the star formation history and the balance between PopIII and PopII stars change in FiBY LW (Johnson et al. 2013 ), not included in the analysis of Section 3.1 . Similar trends are also observed in recent highz simulations (Sarmento & Scannapieco 2022 ; Wells & Norman 2022 ), where PopIII star formation is significantly reduced at z 10 and does not abruptly decrease afterwards, while PopII form at a slower pace because of the delayed metal enrichment, up until z ∼10 when star formation quickly grows and o v ercomes the one in the case where LW radiation is neglected. Secondly, none of the FiBY simulations encompasses enough volume to reach the LW horizon at ∼100 cMpc. Following Ahn et al. ( 2009 )’s approach, we stack the necessary number of copies of the central box until the LW horizon is reached, but in doing so the cosmic variance is certainly lower than the one expected in a full ∼(100 cMpc) 3 box. Finally, several studies have shown that X-rays with energies of the order of ∼1 keV generally fa v our H 2 formation, driving H and He ionization even at the centre of dense gas clouds and thus increasing the abundance of free electrons, that in turn catalyse the formation of H −and H 2 molecules (Inayoshi & Omukai 2011 ; Inayoshi & Tanaka 2015 ). X-rays ultimately counterbalance the effect of the LW radiation, enabling star formation (Haiman et al. 2000 ) and increasing the J crit needed for the DCBH scenario (Glo v er 2016 ). The FiBY simulations, ho we ver, do not include any highz X-ray background due to sources such as massive X-ray binaries and accreting light black holes (see e.g. emission models in Tanaka, Perna & Haiman 2012 ). In conclusion, with this work we hope to provide a useful contribution to the discussion around the LW radiation in the highz Universe and in particular we aim at assisting theoretical astrophysics to include a realistic LWB into their cosmological simulations. The LW radiation is a key ingredient of our current model of galaxy formation at z 10, but still nowadays most of the numerical efforts in the literature do not include it (see e.g. BlueTides Feng et al. 2016 , SPHINX Rosdahl et al. 2018 , OBELISK Trebitsch et al. 2021 , FLARES Lo v ell et al. 2021 ; Wilkins et al. 2023 , THESAN Garaldi et al. 2022 ; Kannan et al. 2022 ; Smith et al. 2022 , ASTRID Bird et al. 2022 ). At the same time, the first recent results from JWST (Donnan et al. 2023 ; Harikane et al. 2023 ) suggest the existence of massive galaxies even at z 15, when the LW background is rapidly growing ( J 21 ∼0.1–1), is still dominated by PopIII stars and greatly affects star formation in low-mass haloes. As for ourselves, we plan to continue studying the build-up of the LW radiation and its interplay with PopIII and PopII star formation in low-mass haloes in a future set of cosmological simulations, tailored to represent the best trade-off between a high spatial and mass resolution and a large volume and to address all the limitations reported in this work. ACKNOWLEDGEMENTS AI acknowledges support from a STFC-ScotDIST studentship. This work made e xtensiv e use of python3 (Van Rossum & Drake 2009 ) and of the following open-source libraries: IPython (Perez & Granger 2007 ), matplotlib (Hunter 2007 ), numpy (Harris et al. 2020 ), astropy (Astropy Collaboration 2013 , 2018 ), scipy (Virtanen et al. 2020 ), h5py (Collette et al. 2021 ), emcee (ForemanMackey et al. 2013 ); we are grateful to the respective communities of developers. AI is grateful to Dr Eric Tittley for keeping the ROE Cuillin HPC cluster healthy and running. This work made e xtensiv e use of the NASA Astrophysics DataSystem, the astro-ph pre-print archive, the Feedly aggregator, and the Mendeley reference manager. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission. DATA AVAILABILITY The data used in this work is available upon reasonable request. REFERENCES Abel T., Anninos P., Zhang Y., Norman M. L., 1997, New Astron. , 2, 181 Abel T., Bryan G. L., Norman M. L., 2000, ApJ , 540, 39 Abel T., Bryan G. L., Norman M. L., 2002, Science , 295, 93 Abgrall H., Roueff E., Launay F., Roncin J. Y., Subtil J. L., 1993a, A&AS, 101, 273 Abgrall H., Roueff E., Launay F., Roncin J. Y., Subtil J. L., 1993b, A&AS, 101, 323 Abgrall H., Roueff E., Launay F., Roncin J. Y., Subtil J. L., 1993c, J. Mol. Spectrosc. , 157, 512 Abgrall H., Roueff E., Liu X., Shemansky D. E., 1997, ApJ , 481, 557 Abgrall H., Roueff E., Drira I., 2000, A&AS , 141, 297 Adelberger K. L., Steidel C. C., Pettini M., Shapley A. E., Reddy N. A., Erb D. K., 2005, ApJ , 619, 697 Agarwal B., Khochfar S., 2015, MNRAS , 446, 160 Agarw al B., Khochf ar S., Johnson J. L., Neistein E., Dalla Vecchia C., Livio M., 2012, MNRAS , 425, 2854 Agarwal B., Dalla Vecchia C., Johnson J. L., Khochfar S., Paardekooper J.-P., 2014, MNRAS , 443, 648 Agarwal B., Smith B., Glo v er S., Natarajan P., Khochfar S., 2016, MNRAS , 459, 4209 Agarwal B., Cullen F., Khochfar S., Ceverino D., Klessen R. S., 2019, MNRAS , 488, 3268 Ahn K., Shapiro P. R., Iliev I. T., Mellema G., Pen U.-L., 2009, ApJ , 695, 1430 Astropy Collaboration, 2013, A&A , 558, A33 Astropy Collaboration, 2018, AJ , 156, 123 Ba ˜ nados E. et al., 2018, Nature , 553, 473 Babb J. F., 2015, ApJS , 216, 21 Begelman M. C., Volonteri M., Rees M. J., 2006, MNRAS , 370, 289 Bhowmick A. K., Blecha L., Torrey P., Kelley L. Z., Vogelsberger M., Nelson D., Weinberger R., Hernquist L., 2022, MNRAS , 510, 177 Bird S., Ni Y., Di Matteo T., Croft R., Feng Y., Chen N., 2022, MNRAS , 512, 3703 Black J. H., Dalgarno A., 1977, ApJS , 34, 405 Bonoli S., Mayer L., Callegari S., 2014, MNRAS , 437, 1576 Bouwens R. J. et al., 2016, ApJ , 833, 72 Bowler R. A. A. et al., 2015, MNRAS , 452, 1817 Bressan A., Fagotto F., Bertelli G., Chiosi C., 1993, A&AS, 100, 647 Bromm V., Larson R. B., 2004, ARA&A , 42, 79 Bromm V., Loeb A., 2003, Nature , 425, 812 Bromm V., Yoshida N., 2011, ARA&A , 49, 373 Bromm V., Coppi P. S., Larson R. B., 1999, ApJ , 527, L5 Chabrier G., 2003, PASP , 115, 763 Chiaki G., Yoshida N., 2022, MNRAS , 510, 5199 Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ , 823, 102 Chuzhoy L., Kuhlen M., Shapiro P. R., 2007, ApJ , 665, L85 Downloaded from https://academic.oup.com/mnras/article/522/1/330/7110413 by Biblioteca Universidad de Sevilla user on 20 May 2024
Modelling the Lyman–Werner background 347 MNRAS 522, 330–349 (2023) Collette A. et al., 2021, h5p y/h5p y: 3.3.0 Corney A., 1977, Atomic and Laser Spectroscopy. Cambridge Univ. Press, Cambridge Cullen F., McLure R. J., Khochfar S., Dunlop J. S., Dalla Vecchia C., 2017, MNRAS , 470, 3006 Dalla Vecchia C., Schaye J., 2012, MNRAS , 426, 140 da Silva R. L., Fumagalli M., Krumholz M. R., 2012, ApJ , 745, 145 da Silva R. L., Fumagalli M., Krumholz M. R., 2014, MNRAS , 444, 3275 Dijkstra M., Haiman Z., Mesinger A., Wyithe J. S. B., 2008, MNRAS , 391, 1961 Dijkstra M., Ferrara A., Mesinger A., 2014, MNRAS , 442, 2036 Donnan C. T. et al., 2023, MNRAS , 518, 6011 Dotter A., 2016, ApJS , 222, 8 Draine B. T., Bertoldi F., 1996, ApJ , 468, 269 Eggenberger P., Meynet G., Maeder A., Hirschi R., Charbonnel C., Talon S., Ekstr ¨ om S., 2008, Ap&SS , 316, 43 Eldridge J. J., Stanway E. R., 2022, ARA&A , 60, 455 Fan X. et al., 2006, AJ , 132, 117 Faucher-Gigu ` ere C.-A., Lidz A., Zaldarriaga M., Hernquist L., 2009, ApJ , 703, 1416 Feng Y., Di-Matteo T., Croft R. A., Bird S., Battaglia N., Wilkins S., 2016, MNRAS , 455, 2778 Fernandez R., Bryan G. L., Haiman Z., Li M., 2014, MNRAS , 439, 3798 Fialkov A., Barkana R., Visbal E., Tseliakhovich D., Hirata C. M., 2013, MNRAS , 432, 2909 F oreman-Macke y D., Hogg D. W., Lang D., Goodman J., 2013, PASP , 125, 306 Frebel A., Johnson J. L., Bromm V., 2007, MNRAS , 380, L40 Fryer C. L., Woosley S. E., Heger A., 2001, ApJ , 550, 372 Galli D., Palla F., 2013, ARA&A , 51, 163 Garaldi E., Kannan R., Smith A., Springel V., Pakmor R., Vogelsberger M., Hernquist L., 2022, MNRAS , 512, 4909 Ge Q., Wise J. H., 2017, MNRAS , 472, 2773 Gessey-Jones T. et al., 2022, MNRAS , 516, 841 Glo v er S. C. O., 2015a, MNRAS , 451, 2082 Glo v er S. C. O., 2015b, MNRAS , 453, 2901 Glo v er S. C. O., 2016, preprint ( arXiv:1610.05679 ) Greif T. H., Bromm V., 2006, MNRAS , 373, 128 Guzzo L. et al., 2014, A&A , 566, A108 Haardt F., Madau P., 2001, in Neumann D. M., Tran J. T. V.eds, Clusters of Galaxies and the High Redshift Universe Observed in X-rays. CEA, Saclay, France, p. 64 Haardt F., Madau P., 2012, ApJ , 746, 125 Haiman Z., Thoul A. A., Loeb A., 1996, ApJ , 464, 523 Haiman Z., Rees M. J., Loeb A., 1997, ApJ , 476, 458 Haiman Z., Abel T., Rees M. J., 2000, ApJ , 534, 11 Harikane Y. et al., 2023, ApJS , 265, 5 Harris C. R. et al., 2020, Nature , 585, 357 Hartwig T., Bromm V., Klessen R. S., Glo v er S. C. O., 2015a, MNRAS , 447, 3892 Hartwig T., Clark P. C., Glo v er S. C. O., Klessen R. S., Sasaki M., 2015b, ApJ , 799, 114 Hartwig T. et al., 2022, ApJ , 936, 45 Heger A., Woosley S. E., 2002, ApJ , 567, 532 Heger A., Woosley S. E., 2010, ApJ , 724, 341 Hirano S., Hosokawa T., Yoshida N., Omukai K., Yorke H. W., 2015, MNRAS , 448, 568 Hirashita H., Ferrara A., 2002, MNRAS , 337, 921 Hunter J. D., 2007, Comput. Sci. Eng. , 9, 90 Iliev I. T., Scannapieco E., Martel H., Shapiro P. R., 2003, MNRAS , 341, 81 Inayoshi K., Omukai K., 2011, MNRAS , 416, 2748 Inayoshi K., Tanaka T. L., 2015, MNRAS , 450, 4350 John T. L., 1988, A&A, 193, 189 Johnson J. L., Dalla Vecchia C., Khochfar S., 2013, MNRAS , 428, 1857 Kannan R., Garaldi E., Smith A., Pakmor R., Springel V ., V ogelsberger M., Hernquist L., 2022, MNRAS , 511, 4005 Kennicutt Robert C. J., 1998, ARA&A , 36, 189 Kobayashi C., Karakas A. I., Lugaro M., 2020, ApJ , 900, 179 Komatsu E. et al., 2009, ApJS , 180, 330 Kulkarni M., Visbal E., Bryan G. L., 2021, ApJ , 917, 40 Larkin M. M., Gerasimov R., Burgasser A. J., 2023, AJ , 165, 2 Latif M. A., Khochfar S., 2019, MNRAS , 490, 2706 Latif M. A., Bovino S., Grassi T., Schleicher D. R. G., Spaans M., 2015, MNRAS , 446, 3163 Latif M. A., Whalen D. J., Khochfar S., Herrington N. P., Woods T. E., 2022a, Nature , 607, 48 Latif M. A., Whalen D., Khochfar S., 2022b, ApJ , 925, 28 Lee K.-S., Giavalisco M., Conroy C., Wechsler R. H., Ferguson H. C., Somerville R. S., Dickinson M. E., Urry C. M., 2009, ApJ , 695, 368 Leitherer C. et al., 1999, ApJS , 123, 3 Liu B., Bromm V., 2020, MNRAS , 497, 2839 Lodato G., Natarajan P., 2006, MNRAS , 371, 1813 Lo v ell C. C., Vijayan A. P., Thomas P. A., Wilkins S. M., Barnes D. J., Irodotou D., Roper W., 2021, MNRAS , 500, 2127 Luo Y., Ardaneh K., Shlosman I., Nagamine K., Wise J. H., Begelman M. C., 2018, MNRAS , 476, 3523 Luo Y., Shlosman I., Nagamine K., Fang T., 2020, MNRAS , 492, 4917 Lupi A., Haiman Z., Volonteri M., 2021, MNRAS , 503, 5046 Machacek M. E., Bryan G. L., Abel T., 2001, ApJ , 548, 509 Madau P., Dickinson M., 2014, ARA&A , 52, 415 Madau P., Rees M. J., 2001, ApJ , 551, L27 Maio U., Ciardi B., Dolag K., Tornatore L., Khochfar S., 2010, MNRAS , 407, 1003 Maio U., Khochfar S., Johnson J. L., Ciardi B., 2011, MNRAS , 414, 1145 McLaughlin B. M., Stancil P. C., Sadeghpour H. R., F orre y R. C., 2017, J. Phys. B Atom. Mol. Phys. , 50, 114001 McLeod D. J., McLure R. J., Dunlop J. S., Robertson B. E., Ellis R. S., Targett T. A., 2015, MNRAS , 450, 3032 McLeod D. J., McLure R. J., Dunlop J. S., 2016, MNRAS , 459, 3812 McLure R. J. et al., 2013, MNRAS , 432, 2696 Mortlock D. J. et al., 2011, Nature , 474, 616 Nagamine K., Choi J.-H., Yajima H., 2010, ApJ , 725, L219 Nakajima K., Maiolino R., 2022, MNRAS , 513, 5134 Neistein E., Khochfar S., Dalla Vecchia C., Schaye J., 2012, MNRAS , 421, 3579 O’Shea B. W., Norman M. L., 2008, ApJ , 673, 14 O’Shea B. W., Wise J. H., Xu H., Norman M. L., 2015, ApJ , 807, L12 Oesch P. A. et al., 2014, ApJ , 786, 108 Oesch P. A., Bouwens R. J., Illingworth G. D., Labb ´ e I., Stefanon M., 2018, ApJ , 855, 105 Omukai K., Tsuribe T., Schneider R., Ferrara A., 2005, ApJ , 626, 627 Omukai K., Schneider R., Haiman Z., 2008, ApJ , 686, 801 Paardekooper J. P., Khochfar S., Dalla C. V., 2013, MNRAS , 429, L94 Paardek ooper J.-P., Khochf ar S., Dalla Vecchia C., 2015, MNRAS , 451, 2544 Park J., Ricotti M., Sugimura K., 2021, MNRAS , 508, 6176 Peebles P. J. E., Dicke R. H., 1968, ApJ , 154, 891 Perez F., Granger B. E., 2007, Comput. Sci. Eng. , 9, 21 Phipps F., Khochfar S., Varri A. L., Dalla Vecchia C., 2020, A&A , 641, A132 Pillepich A. et al., 2018, MNRAS , 475, 648 Qin Y., Mesinger A., Park J., Greig B., Mu ˜ noz J. B., 2020, MNRAS , 495, 123 Raiter A., Schaerer D., Fosbury R. A. E., 2010, A&A , 523, A64 Regan J., 2022, preprint ( arXiv:2210.04899 ) Regan J. A., Johansson P. H., Haehnelt M. G., 2014, MNRAS , 439, 1160 Regan J. A., Johansson P. H., Wise J. H., 2016, MNRAS , 461, 111 Rosdahl J. et al., 2018, MNRAS , 479, 994 Rossi M., Salvadori S., Sk ´ ulad ´ ottir ´ A., 2021, MNRAS , 503, 6026 Salpeter E. E., 1955, ApJ , 121, 161 Salumbides E. J., Bagdonaite J., Abgrall H., Roueff E., Ubachs W., 2015, MNRAS , 450, 1237 Sarmento R., Scannapieco E., 2022, ApJ , 935, 174 Saslaw W. C., Zipoy D., 1967, Nature , 216, 976 Sassano F., Schneider R., Valiante R., Inayoshi K., Chon S., Omukai K., Mayer L., Capelo P. R., 2021, MNRAS , 506, 613 Schaerer D., 2002, A&A , 382, 28 Schauer A. T. P., Glo v er S. C. O., Klessen R. S., Clark P., 2021, MNRAS , 507, 1775 Schaye J., Dalla Vecchia C., 2008, MNRAS , 383, 1210 Downloaded from https://academic.oup.com/mnras/article/522/1/330/7110413 by Biblioteca Universidad de Sevilla user on 20 May 2024
348 A. Incatasciato, S. Khochfar and J. O ˜ norbe MNRAS 522, 330–349 (2023) Schaye J. et al., 2010, MNRAS , 402, 1536 Schmidt M., 1959, ApJ , 129, 243 Shang C., Bryan G. L., Haiman Z., 2010, MNRAS , 402, 1249 Shapiro P. R., Kang H., 1987, ApJ , 318, 32 Smith B., Sigurdsson S., Abel T., 2008, MNRAS , 385, 1443 Smith B. D., Wise J. H., O’Shea B. W., Norman M. L., Khochfar S., 2015, MNRAS , 452, 2822 Smith A., Kannan R., Garaldi E., Vogelsberger M., Pakmor R., Springel V., Hernquist L., 2022, MNRAS , 512, 3243 Solomon P. M., 1965, PhD thesis, The University of Wisconsin-Madison Spinoso D., Bonoli S., Valiante R., Schneider R., Izquierdo-Villalba D., 2023, MNRAS , 518, 4672 Springel V., 2005, MNRAS , 364, 1105 Springel V., Yoshida N., White S. D. M., 2001, New Astron. , 6, 79 Stacy A., Bromm V., Lee A. T., 2016, MNRAS , 462, 1307 Stanway E. R., Eldridge J. J., 2018, MNRAS , 479, 75 Stecher T. P., Williams D. A., 1967, ApJ , 149, L29 Sugimura K., Omukai K., Inoue A. K., 2014, MNRAS , 445, 544 Sugimura K., Coppola C. M., Omukai K., Galli D., Palla F., 2016, MNRAS , 456, 270 Tanaka T., Perna R., Haiman Z., 2012, MNRAS , 425, 2974 Tegmark M., Silk J., Rees M. J., Blanchard A., Abel T., Palla F., 1997, ApJ , 474, 1 Tornatore L., Ferrara A., Schneider R., 2007, MNRAS , 382, 945 Trebitsch M. et al., 2021, A&A , 653, A154 Trenti M., Stiavelli M., 2009, ApJ , 694, 879 Tseliakhovich D., Hirata C., 2010, Phys. Rev. D , 82, 083520 Ubachs W., Salumbides E. J., Murphy M. T., Abgrall H., Roueff E., 2019, A&A , 622, A127 Van Rossum G., Drake F. L., 2009, Python 3 Reference Manual. CreateSpace, Scotts Valley, CA Virtanen P. et al., 2020, Nat. Methods , 17, 261 Visbal E., Bryan G. L., Haiman Z., 2020, ApJ , 897, 95 Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nat. Rev. Phys. , 2, 42 Wells A. I., Norman M. L., 2022, ApJ , 932, 71 Welsh L., Cooke R., Fumagalli M., 2019, MNRAS , 487, 3363 Welsh L., Cooke R., Fumagalli M., Pettini M., 2022, ApJ , 929, 158 Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS , 393, 99 Wilkins S. M. et al., 2023, MNRAS , 519, 3118 Wise J. H., Abel T., 2005, ApJ , 629, 615 Wise J. H., Abel T., 2007, ApJ , 671, 1559 Wise J. H., Abel T., Turk M. J., Norman M. L., Smith B. D., 2012a, MNRAS , 427, 311 Wise J. H., Turk M. J., Norman M. L., Abel T., 2012b, ApJ , 745, 50 Wise J. H., Regan J. A., O’Shea B. W., Norman M. L., Downes T. P., Xu H., 2019, Nature , 566, 85 Wolcott-Green J., Haiman Z., 2019, MNRAS , 484, 2467 Wolcott-Green J., Haiman Z., Bryan G. L., 2011, MNRAS , 418, 838 Wolcott-Green J., Haiman Z., Bryan G. L., 2017, MNRAS , 469, 3329 Woods T. E., Willott C. J., Regan J. A., Wise J. H., Downes T. P., Norman M. L., O’Shea B. W., 2021, ApJ , 920, L22 Xu H., Wise J. H., Norman M. L., Ahn K., O’Shea B. W., 2016, ApJ , 833, 84 Yajima H., Khochfar S., 2017, MNRAS , 467, L51 Yoshida N., Omukai K., Hernquist L., Abel T., 2006, ApJ , 652, 6 Zackrisson E., Rydberg C.-E., Schaerer D., ¨ Ostlin G., Tuli M., 2011, ApJ , 740, 13 Zammit M. C. et al., 2017, ApJ , 851, 64 Zammit M. C. et al., 2018, in Mendoza C., Turck-Chi ` eze S., Colgan J., eds, ASP Conf. Ser. Vol. 515, Workshop on Astrophysical Opacities. Astron. Soc. Pac., San Francisco, p. 145 APPENDIX A: SEDS EXAMPLES In Fig. A1 we compare the different stellar spectra employed in this work. The top panel reports the 1-Myr-old SEDs included in the fiducial setup: Ygg2 for PopIII stars in black and BPASS Chab Figure A1. Top panel: the UV section of the SEDs, rele v ant for the H 2 dissociation, for a 1-Myr-old PopIII (black line) and a PopII (red line) population. Our fiducial choice of stellar models is shown here, hence Ygg2 for PopIII and BPASS Chab for PopII. The total mass is 10 6 M for each population. Bottom panel: ratio of the other SEDs considered to the fiducial ones, again at an age of 1 Myr. Black lines represent PopIII models: Ygg1 (solid), Slug (dot-dashed), and BB5 (dotted). Red lines represent PopII models: BP ASS TH (solid), BP ASS BH (dashed), Slug (dot-dashed), and BB4 (dotted). Apart from the spectral features, that are influenced by the spectral resolution, the most noticeable feature is the completely different shape the all the PopII SEDs from the 10 4 K blackbody spectrum, that is commonly assumed in the literature to approximate them. for PopII stars in red. The UV spectral range shown in the x-axis is the rele v ant one for the H 2 dissociation rate. The additional SEDs described in Section 2.2 are shown in the bottom panel, as ratios relative to the corresponding fiducial choices. PopIII models are in black: Ygg1 (solid), Slug (dot-dashed), and BB5 (dotted). PopII models are in red: BP ASS TH (solid), BP ASS BH (dashed), Slug (dot-dashed), and BB4 (dotted). As expected, a young stellar population has a higher emission if a more top-heavy IMF is chosen (Ygg1 and BPASS TH) and the opposite is true for a bottom-heavy one (BPASS BH). The spectral shape seems consistent within the PopIII and PopII SEDs separately (ratios are mostly parallel to the horizontal line), apart from the BB4, i.e. the blackbody spectrum with T = 10 4 K that is commonly used in the literature to approximate PopII stellar emission (e.g. Johnson et al. 2013 ; Glo v er 2015a ). The red dotted line in Fig. A1 shows instead that the ratio between BB4 and BPASS Chab varies by ∼3 orders of magnitude in the 6–13.6 eV energy interval, hence implying a much softer spectral shape. This has important consequences in the H −detachment rate, that is determined by photons in a wide energy range, from the UV to the IR, the reaction energy threshold being at 0.75 eV. APPENDIX B: SEDS CHOICE Throughout the paper, we have been showing the results for our fiducial choice of stellar SEDs, although in Section 2.2 (Tables 2 and 3 ) we have listed all the additional SEDs included in the postprocessing algorithm. Here we show the potential, but limited, impact of a different choice of stellar emission models on some of our results. In particular, in Fig. B1 we report the ratio between the rates with the other combinations and the ’FID’ one, for the M simulation. The Downloaded from https://academic.oup.com/mnras/article/522/1/330/7110413 by Biblioteca Universidad de Sevilla user on 20 May 2024
Modelling the Lyman–Werner background 349 MNRAS 522, 330–349 (2023) Figure B1. Ratio of the H 2 dissociation (left-hand panel) and H − detachment rate (right-hand panel) of the alternative SED combinations to the ’FID’ one. left-hand panel shows this for the H 2 dissociation rate: at early times the mean value is al w ays within a factor of 2 from the fiducial case and is approximately constant with time (simply reflecting the small differences between PopIII SEDs), while, when PopII stars dominate at z ࣠ 10, the ratio shows a different behaviour that is more prominent in the ’BB’ case, that by construction accounts only for very young stars (<5 Myr) emitting a constant blackbody spectrum at 10 5 and 10 4 K for PopIII and PopII stars, respectively. The latter case leads to underestimating the LWB by up to a factor of 5 at z ∼8. The same ratios are shown for the H −detachment rate in the right-hand panel of Fig. B1 . ’SLUG’, ’TH’, and ’BH’ cases all lie within a factor of 2–3 abo v e or below the fiducial case., while the ’BB’ spectra give a mean rate (dashed turquoise line) more than two orders of magnitude abo v e the ’FID’ case. We have verified that this is due to the 10 4 K blackbody spectrum for PopII stars, normalized as suggested in Greif & Bromm ( 2006 ): such a soft spectrum, as already shown in Latif et al. ( 2015 ) (their fig. 1) and highlighted in Appendix A , when inte grated o v er the wide wavelength range of the H −cross section, gives a rate that is several orders of magnitude abo v e the same rate for a harder spectrum with the same normalization at the Lyman limit. Consequently, our results on the ef fecti ve LW spectral shape (Section 3.1.1 ) do not change appreciably with a different SEDs choice, with the exception of ’BB’ that converges to T eff = 10 4 K even before PopII stars dominate the LWB (Section 3.1.2 ). In conclusion, the radiation background depends only mildly on the choice of the SEDs, as long as realistic stellar models are employed; on the other hand, approximations such as the ’BB’ case gi ve dif ferent results that in turn can lead to inaccurate e v aluations of the H 2 abundance in the early Universe. The contribution from old stellar populations, normally neglected in the literature (Ahn et al. 2009 ; Wise et al. 2012a ; Johnson et al. 2013 ), also depends on the specific choice of stellar IMF and SED. In Fig. B2 we show the fraction of the rates that is due to stellar populations older than 20 Myr, by combining the two ’oldest’ bins described in Section 3.1.2 . The figure again refers only to M, but the same results are valid for the other simulations, with only subtle variations depending on the specific star formation history. As already shown in Fig. 9 , the contribution from old stars increases with time and is larger for the H −detachment rate (dashed lines). This is true for all the combinations that include old stars (namely, for everyone but the ’BB’ case) and the differences can be Figure B2. Similarly to Fig. 9 , fraction of the H 2 dissociation rate (solid lines) and H −detachment rate (dashed lines) in the M FiBY simulation originated by stars older than 20 Myr. The SED choices are colour-coded in the same way as in Fig. B1 . explained with the different number of low-mass stars in top-heavy (magenta-for PopIII-and gold lines) and bottom-heavy (orange and black lines) IMFs. Given the current uncertainties on the IMF in metal-free and metal-poor environments, old stellar populations can account for up to ∼20 per cent –40 per cent of the H 2 dissociation and ∼40 per cent –80 per cent of the H −detachment rate during the Epoch of Reionization, while their contribution is limited to ∼10 per cent –20 per cent at z ࣠ 14. APPENDIX C: NORMALIZED CRITICAL DISTANCE As shown in Section 3.2.1 , the radiation emitted by a single galaxy can exceed the LW radiation intensity in a volume whose size depends on the galaxy stellar mass and the mean LWB le vel. We sho w in Fig. C1 the normalized critical distance D( M ) as defined in equation ( 12 ), that is then fitted as shown in equation ( 13 ). Figure C1. D( M ), i.e. the critical distance shown in Fig. 13 normalized by the mean LW radiation intensity, for all the FiBY simulations. Each line is colour-coded depending on the redshift. With the exception of the galaxies close to the resolution limit, D( M ) depicts a clear increasing trend with the galaxy stellar mass and hence can be fitted with equation ( 13 ). The small evolution with redshift at a fixed M is explained by the lower UV emission per stellar mass from PopII stars with respect to PopIII-dominated galaxies. This paper has been typeset from a T E X/L A T E X file prepared by the author. Downloaded from https://academic.oup.com/mnras/article/522/1/330/7110413 by Biblioteca Universidad de Sevilla user on 20 May 2024