The Role of Quasar Radiative Feedback on Galaxy Formation during Cosmic Reionization

The simulations use the Adaptive Refinement Tree code (Kravtsov et al. 1997, 2002; Rudd et al. 2008), which uses adaptive mesh to achieve high spatial resolution. The simulations presented here are in a volume of $40\,{h}^{-1}\,\mathrm{cMpc}$, which is comparable to the typical size of a quasar proximity zone. The number of dark matter particles is 1024³, with mass resolution 6.2 × 10⁶M_⊙. The grid size of the minimum refinement level is 39 h⁻¹ ckpc, smaller than the filtering length of the IGM at z > 6 (Gnedin & Hui 1998; Gnedin 2000), therefore capable of capturing the gas motion of large-scale structure formation. There are additionally seven levels of adaptive refinement, down to a peak resolution of 100 physical parsecs at z = 6 (Gnedin 2016), capable to resolve galaxy disks in the radial direction. Radiative transfer in the simulation is coupled with gas dynamics using the Optical Thin Variable Eddington Tensor (OTVET; Gnedin & Abel 2001), an approximate momentum tensor method that lowers the computation cost.

All the simulations in this study use the same parameters of star formation and stellar feedback as the fiducial CROC simulations, which have shown good agreement with observed galaxy UV LFs at z = 6 ∼ 10 and the Gunn–Peterson optical depth at z < 6. I briefly describe the star formation and stellar feedback prescription in CROC and refer readers to the methodology paper for more details (Gnedin 2014). The star formation model is motivated by the empirical Kennicutt–Schmidt law, where the star formation rate (SFR) is proportional to the surface density of molecular hydrogen (Schmidt 1959; Kennicutt 1998; Colombo et al. 2018, and references therein). The molecular hydrogen is calculated using the fitting formula from Gnedin & Draine (2014), which is a function of UV radiation intensity, dust-to-gas ratio, local gas density, and density gradient. The stellar feedback is implemented with the standard "delayed cooling" model (Stinson et al. 2006). Briefly speaking, the gas is prevented to cool for ∼10 Myr near a new star particle. With these prescriptions, in CROC simulations, halos below 10⁸M_⊙ do not form any stars because they lack molecular gas.

I study the quasar radiative feedback by running a suite of simulations with and without quasars. I start these simulations from the snapshot of CROC simulation B40E at redshift z = 8.9. In this snapshot, there are five halos with dark matter mass greater than 1 × 10¹¹M_⊙. Under the abundance matching ansatz, it is likely that they host quasars of modest magnitude ${M}_{1500}$ ∼ −24 (Chen & Gnedin 2018), comparable to the quasar observed by Ota et al. (2018).

I choose the most massive halo (dark matter halo mass M_h = 1.6 × 10¹¹M_⊙) at this snapshot as the quasar host galaxy. The quasar is modeled as a massless ultraluminous point source. To follow the motion of the quasar, I track 50 old stars in the gas center of the halo. The quasar is then placed at the median position of these stars. By doing so, the quasar remains in the gas center of the galaxy, which is a widely accepted nature of real quasars.

The current RT solver in my simulation requires all the radiation sources to use the same type of spectrum to speed up RT calculations to couple it to hydrodynamic calculations. Therefore the spectrum of the quasar is the same as stars in the simulations. More specifically, the spectrum has a spectral index of 1.8 above 1 Rydberg and an exponential cutoff above 4 Rydberg (see the solid line in Figure 4 of Ricotti et al. 2002). The light curve of the quasar is described by a simple light-bulb model, in which they have a fixed luminosity when they are turned on. The quasar's luminosity is controlled by the total luminosity of ionizing photons L_Q. At a given a magnitude ${M}_{1500}$, L_Q can vary widely depending on the spectra index, which is poorly understood for z ∼ 6 quasars. According to Pâris et al. (2011), the spectra index for z < 3 quasars can vary from 1.8 to 0.5. For a ${M}_{1500}$ = −24 quasar, this corresponds to L_Q = 3 × 10⁴⁵ to 3 × 10⁴⁷ erg s⁻¹. In this study, I use L_Q = 1 × 10⁴⁶ erg s⁻¹ and 1 × 10⁴⁷ erg s⁻¹ as the luminosities for low and high-luminosity quasars, respectively.

I run several simulations to explore the dependency of quasar radiative feedback on different variables, as shown in Table 1. The NoQ simulation is the reference run without a quasar. Q_z9_L runs with a low luminosity quasar turned on for 309 Myr from z = 8.9 to z = 6.4, a relatively long quasar lifetime. Q_z9_H is the same as Q_z9_L except that the quasar is 10 times brighter. This simulation is expected to show the maximum effect of quasar radiative feedback on galaxy formation. In simulation Q_z7_H I do not turn on the bright quasar until z = 6.8. By comparing this simulation with Q_z9_L and Q_z9_H, we can study the relative importance of quasar lifetime and quasar luminosity on radiative feedback. In the last simulations Q_z9_L_z8off and Q_z9_L_z7off I turn off the quasar at z = 8.2 and z = 6.8, respectively. These simulations allow us to study the behavior of galaxy formation in the "post-quasar" phase.

In Figure 1 I show the neutral hydrogen map of three simulations NoQ, Q_z9_L, and Q_z9_H at z = 8.6, centered on the quasar host halo. In the left panel where the quasar is not turned on, ionized bubbles are relatively small and clustered around the massive halo. While in the middle and right panels, the quasar radiation dominates near the quasar host. At this time the quasar has been shining for t_Q = 30 Myr, and the ionized bubble in the Q_z9_H run (right panel) has grown to a radius ≳1 pMpc.

Zoom In Zoom Out Reset image size

I use HOP (Eisenstein & Hut 1998), a built-in halo finder in yt³ (Turk et al. 2011), to generate halo catalogs. I use the default density threshold of 160 to obtain dark matter halo mass M_h and halo radius r_h. At above 5 × 10⁸M_⊙, the halo catalog is complete.

To quantify the quasar radiation feedback to a specific halo, it is useful to identify the "same" halos across simulations. This is done by matching the positions and mass of halos in two catalogs. If the distance between two halos (in two halo catalogs respectively) are less than 50 comoving kiloparsecs and their halo mass difference is less than 10%, they are identified as the same halo.

I then calculate halo properties like gas mass and gas temperature within r_h for each halo, as well as the unattenuated galaxy magnitude ${M}_{1500}$. ${M}_{1500}$ is calculated using the Starburst99 model (Leitherer et al. 1999). Specifically, with information of age, metallicity, and mass of a star particle, I search a Starburst99 look-up table for its luminosity at 1500 Å, and sum the luminosity for all star particles within r_h of each halo. Finally, I convert the luminosity to AB magnitude. With the information of mass and age of star particles in each halo in the last snapshot, I also obtain the full SFH for each halo.

If we look at how diverse SFHs are in the simulation without quasars, we will later better appreciate how quasars impact the behaviors of star formation. In this section, I use no quasar simulation NoQ to study these "intrinsic" SFHs. Specifically, I examine the SFHs of halos grouped at z = 6.4 into six halo mass bins (5 × 10⁸M_⊙–1 × 10⁹M_⊙, 1 × 10⁹M_⊙–3 × 10⁹M_⊙, 3 × 10⁹M_⊙–1 × 10¹⁰M_⊙, 1 × 10¹⁰M_⊙–3 × 10¹⁰M_⊙, 3 × 10¹⁰M_⊙–1 × 10¹¹M_⊙, > 1 × 10¹¹M_⊙).

Halos have very different SFHs, even for those in a narrow halo mass bins. This is especially the case for low-mass halos, as shown in the upper panel of Figure 2. The halos shown here are from a sphere of 1 pMpc centered on the quasar host at z = 6.4. In the upper panel, I plot the SFHs of halos in the lowest halo mass bin. Black thin lines are SFHs for individual halos. I randomly select a tenth of the halos to plot, due to the large number of small mass halos. The blue lines mark certain percentiles of the stellar mass at each time. The scatter (10% to 90%) in stellar mass is as large as 2 dexes in these low-mass halos. As halo mass increases, this scatter decreases. In the high halo mass bin 3 × 10¹⁰M_⊙–1 × 10¹¹M_⊙, the scatter reduced to (∼0.5 dex), as shown in the lower panel of Figure 2. The larger scatter of stellar mass in lower mass halos could be due to environmental effects. For example, there are more subhalos in lower halo mass bins that accrete mass early but halt later, thus have a higher concentration and more stars formed in the center.

Zoom In Zoom Out Reset image size

I also calculate the mean stellar mass averaging over all halos in each halo mass bin (dashed line). We can see the stellar mass—halo mass ratio ${M}_{\star }/{M}_{h}$ at z = 6.4 increases from 0.0007 for the lowest halo mass bin 5 × 10⁸M_⊙–1 × 10⁹M_⊙ to 0.004 for higher halo mass bin 3 × 10¹⁰–1 × 10¹¹. This is consistent with Behroozi et al. (2019) and is a natural consequence that higher mass halo has higher gas fraction (Gnedin & Kaurov 2014). I briefly note that because of the order of magnitude difference in stellar mass, the mean stellar mass does not trace the median of halo mass, especially in the low halo mass bins (upper three panels of Figure 2). Halos with larger stellar mass drive the mean stellar mass, and because they have smoother evolution, the mean stellar mass evolution also shows a smoother behavior than the percentile stellar mass evolution. When we study EoR, we are usually more interested in the total stars formed. Therefore, the mean stellar mass is more appropriate to use when calculating the SFR in each halo mass bin. In the following sections SFHs are calculated by differentiating the mean stellar mass evolution.

We have seen that the low-mass halos have a large variation in SFHs, which may be a result of environmental effects. This environmental dependence should also lead to a systematic difference in SFHs between halos close to the quasar and far away from the quasar. In Figure 3, I plot the SFHs of low-mass halos with 5 × 10⁸M_⊙–1 × 10⁹M_⊙ into different distance bins. The upper panel shows the halos within 0.5 pMpc, while the lower panel shows the halos in a shell 1.0–1.5 pMpc from the quasar host. Closer to the quasar host, halos tend to form stars earlier and eventually have more stars. This is partially because compared to normal places, in the denser environment subhalo fraction is higher. These subhalos have early accretion, which leads to a higher concentration, and thus a denser gas core where stars are efficiently formed. The SFHs for halos beyond 1.5 pMpc are very similar to the distance bin 1–1.5 pMpc. This distance dependence is weaker for more massive halos. In halo mass bins larger than 3 × 10⁹M_⊙ the distance dependency diminishes until it is nearly unnoticeable.

Zoom In Zoom Out Reset image size

Quasar radiation can impact star formation in surrounding galaxies, but the level of impact depends on many factors, the most prominent two of which are halo mass and radiation intensity. In this subsection I focus on how quasar radiative feedback impacts SFHs in halos with different halo masses, using simulations NoQ and Q_z9_L. In the next section, I study the dependence on quasar luminosity.

In Figure 4, I group halos within 1 pMpc of the quasar host into six mass bins and plot the mean SFHs in simulations with and without the quasar. Halos in the lowest halo mass bin 5 × 10⁸M_⊙–1 × 10⁹M_⊙ show significant suppression of star formation. As the halo mass increases, the differences between different simulations decrease. For halos larger than 3 × 10⁹, there is no significant suppression.

Zoom In Zoom Out Reset image size

It is possible that the quasar only suppresses star formation in halos with very low stellar mass (lighter blue lines in Figure 2) that causes the suppression we see in Figure 4. To examine this, I further plot the percentile SFHs of simulations NoQ and Q_z9_L, which is shown in the upper panels of Figure 5. For the lowest mass bin 5 × 10⁸M_⊙–1 × 10⁹M_⊙ (upper left panel), we can see a systematic delay of SFHs for halos that form stars only after the quasar turns on, as pointed out by the magenta arrow. Note that the 90% line of the quasar simulation disappears in the left panel because these halos have not formed any stars. In halos that have already formed stars before the quasar turns on (e.g., the 10% line), the quasar suddenly suppresses the ongoing star formation. For halos that form stars only after the quasar turn-on time, the quasar delays the star formation.

Zoom In Zoom Out Reset image size

Interestingly, we can find some similarities across different halo mass bins. For example, the delay of the 30% line in the upper left panel and the 70% line in the upper right panel are similar, and so are the 10% line in the upper left panel and the 30% line in the upper right panel. This indicates the stellar mass at t_Q = 0 is a better proxy to the degree of suppression. In the lower panels of Figure 5, I regroup the halos according to their stellar mass at t_Q = 0. In the lower left panel I show halos with stellar mass between 1 × 10⁴M_⊙–1 × 10⁵M_⊙ and the lower right 1 × 10⁵M_⊙–1 × 10⁶M_⊙. I omitted 30% and 70% lines for clarity. We can see that for halos with M_⋆ < 10⁵M_⊙, there is a visible delay in star formation. However, for a higher stellar mass halo this delay diminishes and totally disappears when M_⋆ > 10⁶M_⊙.

We expect quasar luminosity impacts the degree of change in star formation of halos in quasar fields. In this section, I analyze three simulations NoQ, Q_z9_L, and Q_z9_H to study this dependency on quasar luminosity. In the upper panels of Figure 6, I plot the mean SFHs of halos in three halo mass bins 5 × 10⁸M_⊙–1 × 10⁹M_⊙ (left), 1 × 10⁹M_⊙–3 × 10⁹M_⊙ (middle), and 3 × 10⁹M_⊙–1 × 10¹⁰M_⊙ (right) in these three simulations. The red line shows the SFHs of the high-luminosity run Q_z9_H where the quasar is 10 times brighter than that in Q_z9_L. It shows that the higher luminosity quasar has a greater impact on SFHs than the low luminosity quasar, especially in low-mass halos.

Zoom In Zoom Out Reset image size

In observations, SFRs are more accessible, so in the middle panels, I plot the SFRs. We notice that the suppression happens right after the quasar turns on. This short timescale indicates the main suppression mechanism is photodissociation rather than photoevaporation, because the latter operates in the timescale of ∼100 Myr (Shapiro et al. 2004). In our simulation, SFR is directly proportional to the mass of molecular hydrogen, which depends on the LW radiation. When the quasar turns on, the LW radiation immediately increases, causing the suppression of SFR.

In the lower panels, I show the ratio of SFRs of the two quasar runs to NoQ. Within 100 Myr after the quasars turn on, the SFR ratio only differs by less than a factor of two, although the luminosity of the quasars has an order of magnitude difference. This weak dependence is largely due to the complex shielding of dust and the self-shielding of molecular hydrogen, which in our simulation is done using the subgrid model of Gnedin & Draine (2014). In Section 5 I will discuss it in more detail.

Although photoheating does not impact SFR immediately, it alters the property of gas in halos and has a long-term effect on SFH, as we will better appreciate later in Section 3.5, when I show the quasar turn-off simulations. In this section, I show how the quasar radiation alters the gas properties of low-mass halos in simulations NoQ, Q_z9_L, and Q_z9_H. In the upper panel of Figure 7 I select halos within 1 pMpc from the quasar host and show the averaged gas temperature evolution for three low halo mass bins. When close to the quasar, the ionization front moves very fast and heats the gas very efficiently. Therefore, soon after the quasar turns on, it raises gas temperature significantly. Also, the Q_z9_H run shows a higher gas temperature than Q_z9_L. This may be surprising because Abel & Haehnelt (1999) have shown that the gas heating does not depend on radiation intensity explicitly, but rather on radiation spectrum—harder spectrum heats the gas more efficiently than softer spectrum. The dependency on quasar luminosity we see here is mostly due to the higher luminosity quasar creating faster ionization fronts, which travel deeper inside the halos and heat more gas. Because of the higher optical depth inside the halo, the spectrum also becomes harder and heats the gas more efficiently (Abel & Haehnelt 1999). In the meantime, the recombination rate decreases as temperature increases, which further reinforces heating.

Zoom In Zoom Out Reset image size

In the middle panels, I show the evolution of total gas mass. As a consequence of the increased gas temperature, halo mass in the quasar simulations is significantly lower than NoQ run after ∼100 Myr, roughly the sound-crossing time of these halos. This decrease in total gas mass depends on the halo mass. The smaller the halo, the faster and more significant the decrease is. Especially for the lowest mass bin in Q_z9_H, the halos suffer from net gas loss, or "photoevaporation". In contrast, for halos more massive than 3 × 10⁹M_⊙, the decrease is not significant.

Similar to the middle panel, in the lower panels of Figure 7 I plot the mass of dense gas, defined as n_g > 3 × 10⁻²⁴ g cm⁻³, which is slightly above the density threshold (0.01 cm⁻³) that stars can form in our simulation. For the lowest mass bin, the average dense gas mass in Q_z9_H is only 75% of that in NoQ after ∼200 Myr. This shows that photoheating is also an important mechanism to suppress star formation at long timescales.

If a quasar is on for a shorter period of time, we expect it to impact star formation less. In this section I compare simulations NoQ, Q_z9_H, and Q_z7_H to study the star formation dependency on different quasar lifetimes.

In the left column of Figure 8, I plot the SFHs of low-mass halos (5 × 10⁸M_⊙–1 × 10⁹M_⊙) in these three simulations in the top panel. Similar to Figure 6, I group the halos within 1 pMpc from the quasar according to their halo mass at z = 6.4, or cosmic time = 864 Myr. At this time the quasar in Q_z9_H has been on for 309 Myr and that in Q_z7_H for 64 Myr. It is clear that in Q_z7_H, the total suppression at z = 6.4 is significantly less than that in Q_z9_H. Not only the total suppression, but also the instantaneous suppression when the quasar just turns on, is significantly less, as manifested by the milder drop in SFRs shown in the middle panels. The SFR ratios of the two quasar runs relative to the NoQ run are shown in the bottom left panel. For these low-mass halos, the SFR ratio is only ∼0.3 for the early turn on simulation, compared to 0.7 for the late turn on simulation.

Zoom In Zoom Out Reset image size

This large difference in the suppression of SFR is partially because by the time the quasar turns on at z = 6.8, halos have grown significantly. Most of the halos have doubled their halo mass during the 309 Myr from z = 8.9 to z = 6.8. To eliminate the halo growth and study other impact factors, I regroup the halos by their halo mass at the same quasar lifetime, rather than the same cosmic time. In the right column of Figure 8, I plot the same quantities as the left column, but for halos with 5 × 10⁸M_⊙–1 × 10⁹M_⊙ at the same quasar lifetime t_Q = 64 Myr for both simulations Q_z9_H (thin red) and Q_z7_H (thick purple). The thin (thick) blue lines, as references, are from simulation NoQ at the same time as the thin red (thick purple) lines. Interestingly, comparing the thin and thick lines of this figure, we can find that even for halos with the same halo mass, at lower z = 6.4, the average stellar mass is twice as much as those at z = 8.7. In Section 3.2, we have seen that quasar radiation has a very limited impact on high stellar mass halos. The result here is similar: since the stellar mass of z = 6.4 halos is larger, the suppression from quasar radiation is smaller, as shown in the middle and lower panels.

In Figure 9, I examine what the suppression looks like if I group halos by their stellar mass at t_Q = 0. Specifically, in a 1 pMpc sphere centered on the quasar host in simulation NoQ, I select all the halos with stellar mass between 3 × 10⁴M_⊙ and 4 × 10⁴M_⊙ at z = 8.9 (marked by the vertical purple line). Their mean SFH is shown as the thin blue line. Then I select the same halos in the Q_z9_H run, and their mean SFH is shown as the thin red line. Similarly, at a later time z = 6.8 (also marked by the vertical red line), I select halos with the same stellar mass range (3 × 10⁴M_⊙–4 × 10⁴M_⊙) in simulation NoQ. Their mean SFH is shown as the thick blue line. Then I select the same halos in the Q_z7_H run, and their SFHs are shown as the thick purple line. We can find that the thick purple line overlaps very well with the thin red line. This demonstrates that stellar mass at t_Q = 0 is a good indicator of the degree of suppression: Even the quasar turns on at two different redshifts 8.9 and 6.4 separated by 245 Myr, at a fixed quasar lifetime, halos with the same stellar mass experience a similar degree of suppression.

Zoom In Zoom Out Reset image size

It is possible that quasars change the halo properties so significantly that even after the quasar turns off, they still leave a significant impact. In this section, by comparing simulations NoQ, Q_z9_L, Q_z9_L_z8off, and Q_z9_L_z7off, I will show that if the quasar turns off early (t_Q ∼ 60 Myr), the star formation will recover quickly, while if the quasar shines over 200 Myr the star formation will recover slowly.

Similar to the left columns of Figure 8, I analyze the low-mass halos (5 × 10⁸M_⊙–1 × 10⁹M_⊙) within 1 pMpc of the quasar host in simulations NoQ, Q_z9_L, Q_z9_L_z8off, and Q_z9_L_z7off. In the left column of Figure 10, I show their SFHs. In the middle panels, I plot the SFRs and in the lower panels, I plot the ratio of SFR to the reference Q_z9_L run. When the quasar turns off after 64 Myr, the SFR quickly rises, and after ∼200 Myr, it finally catches up to the NoQ simulation where there is no early quasar phase. However, in simulation Q_z9_L_z7off, where the quasar has shone for more than 200 Myr, the SFR only slightly increases, as shown by the brown line in the lower panel.

Zoom In Zoom Out Reset image size

This dependency of recovery speed on the quasar lifetime is most likely due to the degree of photoheating. In the right column of Figure 10, I show the gas temperature, total gas, and dense gas of low-mass halos in these four simulations. When the quasar turns off at z = 8.2 (black lines), the gas temperature soon decreases to be comparable to the NoQ simulation (blue lines). As a consequence, the total gas mass and the dense gas mass do not differ much from the NoQ simulation. However, if the quasar turns off after ∼200 Myr (brown lines), although the temperature drops significantly within ∼60 Myr, the total gas and dense gas are still almost indistinguishable from the Q_z9_L case where the quasar keeps shining. This explains the slow recovery of star formation in the Q_z9_L_z7off simulation.

Because cosmological simulations typically have spacial resolution ≳100 pc, to model star formation, we need to adopt some subgrid recipes. Different subgrid recipes may result in different properties of galaxies. Therefore it is important to discuss some details of the subgrid model in this suite of simulations.

In my simulations, the stars form from molecular hydrogen, and the SFR is proportional to the H₂ density in each cell:

$$\begin{eqnarray*}&&{\rho }_{s}=\displaystyle \frac{{\rho }_{{{\rm{H}}}_{2}}}{\tau },\end{eqnarray*}$$

where τ = 1.5 Gyr is the depletion time. However, H₂ chemistry is not explicitly modeled in the simulations, to a large degree because the spacial resolution 100 pc is too coarse to allow detailed H₂ formation calculation.

One way to calculate the H₂ mass is to use fitting formulas from higher resolution simulations. Here I use the subgrid model developed by Gnedin & Kravtsov (2011). They have run a suite of zoom-in simulations to a Milky Way sized halo. These simulations have a spatial resolution of 65 pc. In resolved scales, the radiation transfer, including the radiation transfer in Lyman–Werner bands, is done using the OTVET RT solver. In subgrid scales, Lyman–Werner radiation I_LW inside each simulation cell is modeled as

$$\begin{eqnarray*}&&{I}_{\mathrm{LW}}={S}_{{{\rm{H}}}_{2}}\times {S}_{\mathrm{dust}}\times {I}_{\mathrm{LW}}^{0},\end{eqnarray*}$$

where ${S}_{{{\rm{H}}}_{2}}$ is the self-shielding factor of molecular hydrogen and ${I}_{\mathrm{LW}}^{0}$ is the "free-space" radiation field returned directly from the OTVET solver. S_dust is the shielding factor from dust,

$$\begin{eqnarray*}&&{S}_{\mathrm{dust}}={e}^{-D{\sigma }_{0}({n}_{\mathrm{HI}}+2{n}_{{{\rm{H}}}_{2}}){L}_{\mathrm{Sob}}},\end{eqnarray*}$$

where D is the dust-to-gas ratio, σ₀ = 2 × 10⁻²¹ cm is the dust cross-section, and ${L}_{\mathrm{Sob}}=\tfrac{{\rho }_{g}}{2| {\rm{\nabla }}{\rho }_{g}| }$. ${S}_{{{\rm{H}}}_{2}}$ is the self-shielding factor from molecular hydrogen. The calculation of ${S}_{{{\rm{H}}}_{2}}$ is complex since it is line absorption instead of continua absorption. This factor depends on gas density, internal velocity dispersion, molecular cloud size, radiation field, etc. Gnedin & Draine (2014), considering the supersonic turbulence and line overlap, parameterize the model including two free parameters: sonic length and the clumping factor. These two free parameters are then fixed by calibrating against observations of Milky Way and Magellanic Clouds. They further ran a parameter grid of UV radiation intensity and dust-to-gas ratio and provided a fitting formula to calculate the H₂ fraction in each cell. This is the model used in my code.

To give some intuitive examples of how this model shapes the galaxies in my simulations, in Figure 13 I show the neutral hydrogen fraction (left panels), molecular hydrogen (middle panels), and temperature (right panels) of a small mass halo in simulations NoQ (upper panels) and Q_z9_L (lower panels) at z = 8.2. At this moment the quasar has been on for 64 Myr. The overlaid arrows are gas velocities, and the red dots are young stars formed within 20 Myr. This halo is 0.43 pMpc away from the quasar host, and at this moment its halo mass is 5.6 × 10⁸M_⊙. Comparing the upper and lower panels, we can find when the quasar turns on, although in the outskirts the hydrogen becomes highly ionized and temperature significantly increases; in the central region the gas remains cold and neutral. There is still molecular hydrogen in the center. Therefore, the star formation does not stop completely despite the halo being small and close to the quasar. Another thing worth pointing out is that stars do not only form in situ. Stellar mass of a halo also grows when it accretes stars from surrounding subhalos. This can be another reason that halos do not stop increasing their stellar mass under quasar radiation.

Zoom In Zoom Out Reset image size

Kashikawa et al. (2007) used 1D spherical simulations (Kitayama et al. 2000, 2001) to calculate the delay time of star formation of halos with masses ranging from 10⁹M_⊙ to 10¹¹M_⊙ under different UV intensities. They modeled the halo formation as a spherical cloud collapse, and defined the star formation time as 10% of the baryonic mass has cooled below 1000 K or density exceeds 10⁶ cm⁻³. In their Figure 8 they showed the delay time between runs without UV radiation and UV radiation turned on at z = 8.3. They concluded that under UV intensity⁴ ${J}_{21}=1$, star formation is completely suppressed for halos smaller than 10⁹M_⊙. At higher J₂₁ = 72, star formation delayed for 100 Myr for halos as massive as 10¹⁰M_⊙.

It is hard to directly compare their results to this study, since in real scenarios halos are not isolated spheres and radiation intensity is not simply proportional to the inverse square of distance, due to the absorption of IGM. Nevertheless, here I choose halos in my simulation that are close to their input halo mass and UV intensity to compare. First, their halos collapse at z = 4.87 while my simulations are run until z = 6.4, so I use the average halo accretion histories from Fakhouri et al. (2010) to translate the halo mass at different redshifts. Halo masses 10⁹M_⊙ and 10¹⁰M_⊙ at z = 4.87 correspond to 6.5 × 10⁸M_⊙ and 5.9 × 10⁹M_⊙ at z = 6.4. As for UV intensity, on average, J₂₁ = 1 corresponds to 1 pMpc away from the quasar in run Q_z9_L and J₂₁ = 72 corresponds to 0.37 pMpc away from the quasar in run Q_z9_H. I check the SFHs of all halos of mass 6 × 10⁸M_⊙ < M_h < 7 × 10⁸M_⊙ within a shell of distance 0.9–1.1 pMpc from the quasar in run Q_z9_L, which results in 147 halos. These SFHs vary widely similar to the upper left panel of Figure 4. Unlike predicted in Kashikawa et al. (2007), they are not completely suppressed, although they all experience an apparent delay. The delay time is not a single value but varies from ∼10 Myr to infinity, with a typical value of ∼100 Myr. To mimic the input value of J₂₁ = 72 and M_h = 10¹⁰M_⊙ at z = 4.87, I check the SFHs of halos of 5 × 10⁹M_⊙ < M_h < 7 × 10⁹M_⊙ within 0.35 pMpc < d < 0.4 pMpc from the quasar in run Q_z9_H. There are only four halos that meet this requirement. Their SFHs flattens immediately after the quasar turns on, but star formation resumes again within less than 50 Myr.

Solely comparing the star formation delay time, qualitatively, the results agree. However, one important thing we learn from this cosmological simulation suite is that even though the halo mass and UV intensity is the same, the delay time in star formation varies because of the diversity in halo SFHs. To show this point more clearly, I select halos in the same mass bin 5 × 10⁸M_⊙ < M_h < 1 × 10⁹M_⊙ and under similar UV intensity in two quasar simulations. In Figure 14, plotted in the left panel in orange are SFHs of halos within the shell of 0.3–0.4 pMpc from the quasar host in run Q_z9_L, while in the right panel in red are SFHs within the shell of 1–1.32 pMpc in run Q_z9_H. I also plot the same halos in the NoQ run in thin blue lines for comparison. We can see a difference in the "delay time" between these two. This is consistent with Figures 3 and 5: the halos closer to the quasar (left) form stars earlier, and thus become more "resilient" to quasar radiation.

Zoom In Zoom Out Reset image size