The CGM–GRB Study. I. Uncovering the Circumgalactic Medium around GRB Hosts at Redshifts 2–6

Various methods have been employed in past and ongoing observations to extract the diagnostic features of the multiphased CGM. The most popular technique involves using a bright background quasar (QSO) to trace the CGM around an intervening galaxy. Various local CGM surveys (z < 0.5) including COS-Halos (Tumlinson et al. 2013; Werk et al. 2016; Prochaska et al. 2017), COS-Dwarfs (Bordoloi et al. 2014), COS-GASS (Borthakur et al. 2015), and others (Stocke et al. 2013; Zhu & Ménard 2013) utilize UV/optical absorption spectra to study CGM kinematics and physical properties through high-ionization potential species (high-ion) such as O vi, N v, C iv, and Si iv, and low-ionization potential species (low-ion) such as Fe ii, Si ii, C ii, and Ca ii in the CGM. These observations allow the absorbers to be matched to their respective impact parameters from the galaxy with a precision of tens of kiloparsecs. The higher redshift surveys of QSO–galaxy pairings (e.g., Fox et al. 2007; Lehner et al. 2014; Turner et al. 2014; Rudie et al. 2019), QSO–QSO pairings (QPQ; Hennawi et al. 2006; Prochaska et al. 2014; Lan & Mo 2018), and galaxy–galaxy pairings (Steidel et al. 2010; Lopez et al. 2018) use rest-frame UV spectra for similar analysis with limited information about the associated galaxies and/or impact parameters.

In "down-the-barrel" spectroscopy, a star-forming galaxy's own starlight is used as a background illumination to detect absorption from the intervening ISM and CGM (Martin 2005; Steidel et al. 2010; Bordoloi et al. 2011; Kornei et al. 2012; Rubin et al. 2012, 2014; Heckman et al. 2015; Rigby et al. 2018). This technique has been successful in tracing the galactic inflows and outflows in the CGM of star-forming galaxies with the caveat that the radial coordinate of the absorbing component remains unconstrained.

Our approach utilizes the spectra of bright afterglows of long gamma-ray bursts (GRBs) to derive the kinematic properties of the CGM around their host galaxies and use a simple toy model to further constrain the physical properties of the CGM gas. This technique will further enable the investigation of possible relations between the CGM and the galaxy properties that may govern feedback processes and their evolution with redshift. Thanks to the nature of GRBs and the extensive follow-up effort over the past 10 years, we were able to collect enough data to study the CGM of GRB hosts in the redshift range 2 ≲ z ≲ 6.

Long-duration GRBs are the most powerful explosions in the universe, being several orders of magnitude brighter than typical supernovae. The prompt emission is followed by a rapidly fading (∼1–2 days) X-ray, UV, and optical afterglow. GRB afterglows have been detected from low redshift ∼0.01 out to a redshift of 8.2 (Salvaterra et al. 2009; Tanvir et al. 2009; Cucchiara et al. 2011), thus probing the first few billion years of the universe, an era characterized by the formation and early evolution of galaxies that may have had a critical role in enriching the universe with metals.

Because GRB afterglows are bright background sources, their spectra harbor absorption features produced by the material along the line of sight (LOS) including the host galaxy ISM as well as the CGM and intergalactic medium (Prochaska et al. 2007; Fox et al. 2008; Prochaska et al. 2008a; Cucchiara et al. 2015). Thus, it provides an excellent opportunity to probe the chemical composition and physical conditions of the entire galaxy ecosystem. The GRBs also fade rapidly, making it possible to study the host galaxy component separately, in the absence of the bright GRB.

Compared to QSO sight lines, GRB sight lines have key advantages: (1) GRBs happen within their host galaxy, thus probing the host galaxy ISM as well as the associated CGM, the main components of the galaxy ecosystem; (2) GRB discovery is based on gamma- and X-ray detection and hence, is largely independent with respect to host galaxy properties (further discussed in Section 2.2); (3) the optical afterglow fades in one to two days, clearing the path for future deep observations of the host galaxy to study its global properties and surroundings.

However, at the same time, the number of GRBs suitable for CGM investigation is small, despite the detection rate of 100 yr⁻¹ by dedicated space-based missions like the Neil Gehrels Swift Observatory (Swift). Similarly, the fast-decaying nature of their afterglows requires opportune follow-up strategies in order to obtain high signal-to-noise ratio (S/N) spectra. Finally, it requires separating the CGM and ISM contributions in the absorption spectra, which can be challenging. Our approach tackles this problem by using the kinematic information derived by absorption-line spectroscopic data.

In this paper, we present the analysis of a data set of 27 high-z (z ≳ 2) GRBs, out of which six are at z ≳ 4 (Section 2, Section 3). We use column density line profiles to study the kinematics and line ratios of the absorbing gas to distinguish between the CGM and ISM (Section 4). We designed a simple toy model (Section 5) to obtain an estimate of the outflow properties in the CGM. We estimate the CGM mass for a typical GRB host in this sample and summarize its possible evolution with redshift in Section 6. In Section 7, we discuss various implications our findings may have on our current understanding of the CGM kinematics, outflow rates, metal enrichment, and CGM–galaxy coevolution. The key conclusions of this study are summarized in Section 8.

Obtaining a good estimate of ionic column density as a function of velocity requires high-S/N (≳ 5) spectra of medium resolution (R ≳ 8000, ∼50 km s⁻¹). The cutoff in resolution is selected to clearly distinguish the absorption systems at different velocities and minimize the errors in column density estimation due to saturation effects and blending (see Prochaska et al. 2006; Cucchiara et al. 2013). Our CGM–GRB sample consists of 27 long GRBs that satisfy these stringent criteria. The sample properties are summarized in Figure 1. The median redshift of the sample is 2.71 and the median S/N is 10. Notably, six of the GRBs are at z > 4. Table 1 lists the GRBs in this sample along with the observational details.

Zoom In Zoom Out Reset image size

As mentioned earlier, the transient nature of GRBs requires rapid-response facilities capable of observing their afterglow within a few minutes of their discovery. Our data set comprises primarily of archival data acquired by the X-Shooter and UVES spectrographs on the Very Large Telescope (VLT). These spectra provide a wide wavelength coverage (from optical to near-infrared (NIR)) and sufficiently high spectral resolution (R ∼ 8000–55,000). In addition, we retrieve spectra from the archival data set of the Keck telescope's HIRES and ESI spectrographs.

The majority of these spectra were reanalyzed and normalized using the data analysis pipelines in Cucchiara et al. (2013, 2015). More recent data (from 2014) were acquired from the PHASE 3 VLT archive,^{5 ,6} which provides fully reduced, research-ready one- and two-dimensional spectra. We utilized the flux-calibrated one-dimensional spectra and normalized the GRB-afterglow continuum using a spline function. Every spectrum is manually inspected, and the overall continuum is determined using the Python-based linetool package.⁷ The error from the continuum fit is propagated into the flux error spectrum.

Due to the high-S/N and resolution requirement, this sample is biased toward the brighter end of the GRB-afterglow distribution. The afterglow magnitudes at the time of taking the spectra are listed in Table 1, clearly indicating a limit at m_AB ∼ 21.0 (with the exception of GRB 100219). This selection effect is complex because the magnitude at the time of observation depends on the intrinsic brightness and distance of the afterglow as well as the time elapsed since the prompt gamma-ray emission. Regardless, it can be said that this sample selectively avoids intrinsically faint afterglows. However, considering that the apparent brightness of the afterglow also depends on the host galaxy dust extinction, it can be said that this sample selectively avoids heavily dust-obscured (A_V > 0.5) sight lines (Perley et al. 2009; Krühler et al. 2011; Zafar et al. 2018).

In general, long GRBs trace cosmic star formation (Greiner et al. 2015; Schady 2017). At z ∼ 2.5, the typical star formation rate of GRB hosts is ∼10 M_⊙ yr⁻¹ (Krühler et al. 2015), the typical GRB host stellar mass is ∼10^9.3 M_⊙ (Perley et al. 2016), and the typical gas-phase metallicity is ∼0.05–0.5 solar (Trenti et al. 2015; Arabsalmani et al. 2018). Thus, from a CGM perspective, this sample traces star-forming, low-mass galaxies at z > 2.

In order to infer the kinematics of different chemical species, the redshifts of the GRB host galaxies in the sample need to be determined in a precise and uniform manner. Commonly, the galaxy redshifts are obtained using nebular emission lines (e.g., Hα), but this method is not viable for faint GRB hosts at high redshifts. Therefore, we use the local GRB environment (within a few × 100 pc of the GRB) as a proxy for the systemic redshift of the GRB host galaxy.

The fine-structure transitions of the species such as such as Ni ii*, Fe ii*, Si ii*, and C ii* in the rest-frame UV (Bahcall & Wolf 1968) trace the ISM clouds in the vicinity of the GRB (∼few × 100 pc–1 kpc) due to UV pumping. This is further corroborated by temporal variations found in the strength of these lines in multi-epoch spectra of a few GRBs (Vreeswijk et al. 2007; Hartoog et al. 2013; D'Elia et al. 2014). The strongest absorption components (i.e., velocity components) of these fine-structure transitions are therefore good proxies for the rest-frame velocity of the burst environment within a few hundred parsecs (Chen et al. 2005; Dessauges-Zavadsky et al. 2006; Prochaska et al. 2006). Therefore, we choose the redshift by visual inspection such that the strongest absorption components in the fine-structure transitions occur at ∼0 km s⁻¹, i.e., rest frame. We primarily use Si ii* and C ii* transitions for estimating the redshift. In the case of saturation or confusion between Si ii* and C ii*, we use Ni ii* lines due to their lower oscillator strength.

In addition, we visually check for the presence of other low-ionization lines such as Fe ii 1608, Si ii 1526, and Al ii 1670 (especially the weak transitions of Si ii and Zn ii), which are reliable tracers of the host galaxy ISM (Prochaska et al. 2006; Cucchiara et al. 2013), to visually confirm the redshift determination. With the use of the strongest components of fine-structure transitions to estimate the zero point, the GRB redshift is accurate to within 50 km s⁻¹.

GRB spectra show a plethora of signatures, ranging from the circumburst and interstellar media to the galactic winds and CGM. Prior studies have extracted the intervening systems along the GRB sight lines (at z < z_GRB) to trace cosmic chemical evolution. These studies have used the doublets from Mg ii, C iv, and Si iv to determine the chemical enrichment of the universe similar to the analysis of quasars intervening systems (Fox et al. 2008; Prochaska et al. 2008a; Fynbo et al. 2009; Simcoe et al. 2011; Thöne et al. 2012; Sparre et al. 2014; Cucchiara et al. 2015; Vergani et al. 2017). In this paper, we focus on the contrast between the high- and low-ion kinematics(see Figure 2). The GRB spectra were normalized, binned, and fitted to extract the column density, Doppler parameter, and the line center (in velocity space) for each absorbing component within a velocity window of ±400 km s⁻¹ for key high- and low-ion species. The parameters of these species were then used to study the kinematics of the absorbing gas and estimate the likely origin of the absorbing component(s).

Filtering: Various intervening systems previously reported in the literature (references in Table 1) were identified, and the regions where the absorption lines from the intervening system blend with rest-frame GRB absorption were flagged. In addition, the regions containing telluric absorption from the atmosphere were flagged. For lines with strong neighboring transitions (e.g., Si ii 1260 and S ii 1259, etc.), the velocity windows considered for fitting were carefully adjusted to minimize the confusion.

Voigt profile fitting: Thanks to the medium-resolution spectra in this sample, it is possible to resolve the kinematics of absorbing systems residing in the ISM and CGM into individual components in the velocity space. The GRB spectra were analyzed by fitting individual components, where the optical depth of each component is modeled as a Voigt profile. Given the complex nature of GRB sight lines, typically more than five absorption components are observed in these spectra. While there are other χ² grid-search-based codes (e.g., VPFIT⁸ and MPFIT Markwardt 2009) for fitting the Voigt profiles, these methods have fundamental limitations to solve this particular problem. Obtaining useful results with a grid search in such a large parameter space of a nonlinear model is computationally expensive. Further, it is difficult to capture the degeneracy between various parameters in a quantitative manner, for instance, the degeneracy between Doppler parameter and column density for saturated components. The Bayesian approach provides a rigorous way to visualize the degeneracy and estimate the errors around the fit parameters in a systematic way for multicomponent Voigt profiles.

For these reasons, we developed a Bayesian-inference-based, multicomponent Voigt-profile-fitting code in Python to determine the best-fit values of the parameters for each component (i), i.e., the column density (N_i), Doppler parameter (b_i), and line center (v_i). To sample the posterior probability distributions, the Markov Chain Monte Carlo (MCMC) method was used. The MCMC sampling was implemented using emcee library in Python (Foreman-Mackey et al. 2013). A detailed summary of the Bayesian Voigt-profile fitting method and error estimation is given in Appendix A. This code provides a robust and novel approach to fit the complex multicomponent absorption systems such as GRB or QSO sight lines and obtain reliable estimates of the optimal parameters and the associated errors.

To fit multicomponent Voigt profiles to a given transition, a velocity window spanning ±400 km s⁻¹ around the GRB rest frame is extracted. This velocity window enables fair comparisons with the simulations of the CGM as well as previous observations of high- and low-redshift CGM. The spectra are binned by a factor of 2 to 4 depending on the noise level for easier visual inspection. The number of components to be fitted is determined through manual inspection of doublets (e.g., C iv 1548 and 1550) and lines with similar ionization potential (e.g., C iv, Si iv). The line-spread function of the spectrometer is modeled as a Gaussian function and convolved with the synthesized multicomponent Voigt profile to obtain the comparison spectrum for evaluating the residuals. The initial guesses of the parameters are manually provided to the MCMC routine (as priors) to find the optimal line parameters and associated uncertainties corresponding to all the components (N_i, b_i, v_i). Also, doublets are fitted simultaneously. Our optimal parameters are consistent with the results from other references in Table 1 as shown in Figure 3. The Voigt profile fits for all the spectra are shown in Figures 17–43 of Appendix B, and the line profile parameters are listed in Table 2.

Caveats: The kinematic resolution of the analysis is linked to the spectral resolution of each observation, leading to a variance in the precision of the parameter estimation throughout the sample, due to different spectral resolutions. For warm CGM (T ∼ 10^5–6 K), the thermal Doppler parameter is expected to be in the range 5–20 km s⁻¹ for C and Si, with further broadening expected due to turbulence (Lehner et al. 2014). For spectra with R ∼ 10,000, we can marginally resolve b ∼ 15 km s⁻¹. In saturated regions, the determination of optimal parameters has a higher uncertainty, due to degeneracy between the Doppler parameter and column density. In such cases, lines with similar ionization potential and weaker oscillator strength help provide an estimate without breaking the degeneracy. While the Voigt profile fitting works well in the case of mildly saturated lines, it does not fully alleviate the uncertainty for strongly saturated lines (such as saturation spanning ∼100 km s⁻¹). In Table 2, a quality flag of 0 indicates unsaturated or mildly saturated components while a quality flag of 1 indicates strongly saturated components. The tabulated profile parameters for strongly saturated components denote a likely but nonunique solution.

In order to understand the overall kinematics of the sample, the normalized rest-frame spectra were plotted in velocity space. In Figure 2, various high- and low-ion transitions of the individual GRBs are shown in blue, while the median kinematic profiles are shown in green. The key qualitative results from the median plots are:

• (a)  
  There is a significant blueward absorption excess at velocities v ≲ −100 km s⁻¹, which is a clear signature of outflowing gas.
• (b)  
  This blueward asymmetry is stronger in the high-ion transitions than in the low-ion transitions. The median profile for the fine-structure transition of Si ii* is, not surprisingly, more symmetric.
• (c)  
  The median absorption for the low-ion transitions is fairly limited to within ±100 km s⁻¹, unlike the high-ion lines, which extend much farther (especially blueward).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These qualitative observations may indicate outflowing gas predominantly traced by the high-ion transitions. Thus, we expect at least two different phases that are kinematically distinct. In order to test this hypothesis in greater detail, we use a toy model to reproduce the observed kinematic behavior (Section 5).

The absorption lines in the afterglow spectra from various high- and low-ion species were fitted with multicomponent Voigt profiles as described in Section 3.2. In order to quantitatively measure and compare the kinematics of high- and low-ion species as well as compare the observations with models (as described in Section 5), the fitted continuous line profiles were converted into integrated column density profiles. The fitted Voigt profiles profiles (with N_i, b_i, and v_i) were converted to apparent optical depth profiles as a function of velocity (τ_a(v)) such that τ_a(v) = ln [F_c(v)/F(v)], where F_c(v) is the normalized continuum flux level (i.e., 1) and F(v) is the normalized flux from the fitted profile at velocity v (Savage & Sembach 1991). The apparent column density is then evaluated as N_a(v) = 3.768 × 10¹⁴τ_a(v)/fλ, where f is the oscillator strength and λ is the rest-frame wavelength of the line in angstrom. The integrated column density is the integral of N_a(v) over bins of 100 km s⁻¹.

The integrated column density profiles are shown in blue in Figure 4. The strongly saturated components are treated as lower limits evaluated by imposing a maximum cap on τ_a(v) of 4.5 (equivalent to a lower cap of ∼0.01 on the normalized flux), and they are marked as open circles. The saturation issue does not significantly affect the velocity bins beyond ±100 km s⁻¹. The downward triangles show integrated column density for a bin with no detected absorption and is evaluated using τ(v) = 0.05, which denotes the typical detection limit of the sample. The error bars on the integrated column density are evaluated by calculating the difference between the Voigt profiles generated using the optimal parameters and the profiles generated using the 1σ deviation parameters (as shown in Appendix A). This captures both the fitting uncertainty as well as the noise spectrum. The median integrated column density profile for each line is shown in blue. The inner and outer purple-shaded regions show the central 50th and 80th percentiles of the integrated column density distribution for each bin, respectively.

The integrated column density profiles help visualize the fraction of sight lines where various species are detected as a function of velocity. Broadly speaking, the detection fractions can be categorized in three kinematic regions: central ($| v|$ < 100 km s⁻¹), blue wing (v < −100 km s⁻¹), and red wing (v > +100 km s⁻¹). The detection fractions of various high- and low-ion species are calculated as the number of detected sight lines divided by the number of sight lines with the spectral coverage for that ion and are reported in Table 3. The O vi and N v absorption lines are redshifted to the low-sensitivity (blue) regions of the spectrographs for z ∼ 2, which leads to a poor S/N. In addition, the vicinity to the Lyα forest and possible doublets contamination can prevent reliable detection of these species. Therefore, we only consider the detections where both doublets show absorption lines that can be reasonably fitted. The uncertainty on the detection fractions for these ions reflects the small number of cases where these lines are within the spectral coverage and are not affected by low S/N, saturation, or mismatch between the doublets.

Table 3.  Detection Fractions in Various Kinematic Regions

Species	Blue wing	Central	Red wing
	(v < −100 km s−1)	($| v|$ < 100 km s−1)	(v > 100 km s−1)
O vi	${0.65}_{-0.0}^{+0.2}$	${0.8}_{-0.0}^{+0.2}$	${0.6}_{-0.0}^{+0.2}$
N v	${0.52}_{-0.0}^{+0.16}$	${0.76}_{-0.0}^{+0.16}$	${0.56}_{-0.0}^{+0.16}$
Si iv	${0.92}_{-0.0}^{+0.04}$	${0.96}_{-0.0}^{+0.04}$	${0.69}_{-0.0}^{+0.04}$
C iv	0.96	1.00	0.63
Al iii	0.67	1.00	0.62
Si ii	0.73	0.96	0.61
Fe ii	0.67	0.92	0.55

Download table as:  ASCIITypeset image

The central region within ±100 km s⁻¹ shows the highest detection fraction for all the ions. The blue wings have a higher detection fraction compared to the red wings for all ions except N v, as expected from the stacked spectra (Figure 2). We compare the detection fractions in this sample with the KODIAQ sample—a QSO-based survey of CGM absorbers at z ∼ 2–3.5 (Lehner et al. 2014). The overall detection fractions in the KODIAQ survey for O vi (75%), N v (55 ± 20%), Si iv (95%), C iv (90%), and Fe ii (78%) are consistent with the detection fractions in the blue wing region of the CGM–GRB sample within ∼10%.

It is evident from the stacked spectra in Figure 2 that there is a clear excess of absorption in the blue wing (v < −100 km s⁻¹) relative to the red wing (v > 100 km s⁻¹). The excess also appears to be stronger for the high-ionization species such as Si iv and C iv relative to the low-ionization species such as Si ii and Fe ii. We quantify this kinematic asymmetry using the median integrated column density profiles shown in Figure 4. The median of the total integrated column densities are tabulated in three distinct kinematic regions in Table 4: central ($| v|$ < 100 km s⁻¹), blue wing (v < −100 km s⁻¹), and red wing (v > 100 km s⁻¹). The reported uncertainties are evaluated as 2σ intervals of the distribution of the medians of the total integrated column densities derived via a simple resample-with-replacement bootstrap technique. The degree of asymmetry is calculated as the ratio of blue- to red-wing column densities. The median asymmetry in the high ions C iv and Si iv are 0.36 and 0.52 dex, respectively, in contrast to ∼0.17 and 0.19 dex for the low ions Si ii and Fe ii, respectively.

Zoom In Zoom Out Reset image size

The blue absorption excess (or asymmetry) has been interpreted as a galactic outflow signature in previous GRB-afterglow sight-line studies (e.g., with seven high-z GRBs in Fox et al. 2008 and GRB 080810 in Wiseman et al. 2017a). From the kinematic distribution in the blue wing, the typical outflow velocity is ∼150–250 km s⁻¹. Another important aspect is the higher degree of asymmetry in Si iv and C iv relative to Si ii and Fe ii. These high-ionization species are more asymmetric by roughly a factor of ∼1.5–2 than the low-ionization species. This key observation implies that not only is there an excess of outflowing gas compared to infalling gas, but also the outflowing gas is more highly ionized than the infalling gas. Unlike all other ions presented here, N v has no significant asymmetry and therefore could be tracing a phase (or a combination of phases) that is kinematically distinct from the gas phase (or a combination of phases) traced by the other high- and low-ion lines. We return to this point in Section 7.

The ratio of column densities of various high- and low-ion species provides another perspective to learn about the physical conditions of the intervening gas in different kinematic regions. We select Si iv, C iv, Si ii, and Fe ii for this analysis because they have excellent spectral coverage in the sample and are fit reasonably well due to the high S/N in the corresponding observed wavebands. Si ii 1526 and Fe ii 1608 are taken as representative low-ion lines because (a) they have a moderate line strength, thus preventing saturation, and (b) there are no adjacent strong lines in the rest frame that could potentially blend/contaminate the ±400 km s⁻¹ region.

Figure 5 shows the high- to low-ion line ratios as a function of velocity. The line ratios are evaluated as the ratio of the integrated column densities for the two lines as a function of velocity. The solid red symbols indicate that at least one of the two absorption lines is detected in that velocity bin. The empty symbols indicate that one of the two lines is saturated, while the green circles show the cases where both lines are saturated in that velocity bin (double saturation points). The double saturation ratios are evaluated by taking the ratios of the integrated column densities by putting a lower limit on the observed flux (as described in Section 4.2). A circle denotes a detection of both lines, an upward (downward) triangle denotes lower (upper) limits. A cross symbol indicates points where both the lines have nondetections (double nondetection points). In such cases, the ratio is taken as the ratio of their typical detection thresholds evaluated using τ(v) = 0.05 (see Section 4.2). The blue line traces the median of the line ratios including the double nondetection points with the purple shades spanning the 50th and 75th percentile zones. The cyan line traces the median excluding the double nondetection points (i.e., requiring the detection of at least one of the two lines in the given velocity bin).

Zoom In Zoom Out Reset image size

To avoid a large number of double nondetection points, we focus on velocity bins from −250 to 150 km s⁻¹, where the double nondetection cases are limited to less than 40%. In this region, the high- to low-ion ratio is higher in the blue wing relative to the red wing. This is more clearly noticeable in the Si iv/Fe ii and C iv/Fe ii ratios. The actual ratios in the blue wing are likely to be higher, due to a large number of lower limits in the blue wings (i.e., detection of high ions and nondetection of low ions). This can also be seen by comparing the high- and low-ion lines in Table 4. The line ratios in the central region are more uncertain, due to the high occurrence of double saturation cases, but they appear to be commensurate with the blue wing ratios. Qualitatively, the line ratios highlight the distinct physical characteristics of the three kinematic zones and hint toward a general presence of high-ion rich outflowing gas at a projected speed of ∼150–250 km s⁻¹.

To simulate how the GRB sight lines sample the kinematics of the CGM–galaxy system, we constructed a simple geometrical model of the system as shown in Figure 6. The galaxy ISM is modeled as a disk with an exponential gas density distribution in the radial direction. The ISM parameters of the galaxy are chosen to reflect a typical, representative galaxy from the known population of GRB hosts in the redshift range z ∼ 2–3 (Wainwright et al. 2007; Blanchard et al. 2016; Perley et al. 2016; Arabsalmani et al. 2018). The CGM is modeled as an isothermal sphere populated by clouds where the density of the cloud is inversely proportional to the square of the radial coordinate of the cloud (see Equation (17)). This distribution assumes that the clouds originated from a mass-conserving outflow (Steidel et al. 2010; Chisholm et al. 2017). In addition to the isothermal velocity distribution, a certain fraction of clouds (f_out) are randomly selected to have an additional radial outward velocity component to simulate the outflows. The value of the additional radial velocity depends on the radial coordinate of the cloud and follows a ballistic velocity profile with a radial launch velocity of v_out at an outflow launching radius R_launch of 2 kpc. The entire simulation setup is described in detail in Appendix B.

Zoom In Zoom Out Reset image size

To create a complete sample of sight lines for each of these models, 200 GRB sight lines were synthesized by randomly sampling uniform distributions of (a) GRB positions within the galactic disk and (b) the 3D direction vectors of the sight line. The projected velocity of the intervening gas (from the ISM and CGM) along the GRB sight line was calculated with respect to the projected velocity of the gas in the immediate vicinity of the GRB. Setting this velocity reference is important to maintain consistency with the observations, where v = 0 is assigned to the strongest fine-structure line absorption, a tracer of gas in the vicinity of the GRB as evident from the UV-pumping argument (Prochaska et al. 2006 and Vreeswijk et al. 2007).

The most important model parameters affecting the observed CGM–ISM kinematics are listed in Table 5. We approximate the typical stellar mass of GRB hosts at z > 2 as ∼2 × 10⁹ M_⊙ (Perley et al. 2016), thereby constraining the halo mass (Hopkins et al. 2014; Wechsler & Tinker 2018) and thus the typical rotation and dispersion velocities. The volume-filling fraction is approximated as 0.1 from prior CGM studies at lower redshifts (Stocke et al. 2013; Werk et al. 2016). With these assumptions, the free parameters of the model are the CGM mass (traced by C iv), M_CGM; the launching velocity of the outflow, v_out; and the outflow fraction, f_out. Hence, we synthesize a matrix of 27 models with three distinct values of each of these model parameters, as stated in Table 5.

The C iv kinematics are evaluated by making certain assumptions about the CGM and ISM metallicities and the ionization fraction for C iv, f_{C iv}. The C iv column density is evaluated as

$$\begin{eqnarray}&&{N}_{{\rm{C}}{\rm{IV}}}={N}_{\mathrm{total}}\times \displaystyle \frac{{f}_{{\rm{C}}{\rm{IV}}}}{0.3}\times \displaystyle \frac{Z}{{Z}_{\odot }}\times {\left({n}_{{\rm{C}}}/n\right)}_{\odot }.\end{eqnarray} \tag{ 1 }$$

We conservatively assume f_{C iv} = 0.3 as a maximal ionization fraction similar to Bordoloi et al. (2014) and as derived in Oppenheimer & Schaye (2013) for z > 2. The CGM and ISM metallicities are assumed to be 0.25 and 0.15 solar, respectively (Shen et al. 2013; Trenti et al. 2015; Arabsalmani et al. 2018). The solar ratio of carbon number density to total number density is 3.26 × 10⁻⁴.

The kinematic separation of the CGM and ISM has been used in various previous studies (e.g., to study the CGM of Andromeda galaxy) to distinguish the CGM and ISM contributions in the velocity profiles (Fox et al. 2014, 2015; Lehner et al. 2015). In this paper, we performed kinematic simulations to infer the relative contributions of the CGM and ISM as a function of velocity using the characteristic parameters of GRB host galaxies of our sample.

The breakdown of the total column density into ISM and CGM components, with and without an outflow, are shown in Figure 16 in Appendix B. Two key inferences can be drawn from this illustration: (a) an outflow component is necessary to explain the observed kinematic asymmetry, and (b) the ISM component kinematically dominates the central region ($| v| \lt 100$ km s⁻¹), while the CGM component is the main contributor of column density in the red and blue wings (v > 100 km s⁻¹ and v < −100 km s⁻¹, respectively). With this insight, we use the kinematic separation to further estimate the physical properties of the CGM.

Figures 7–9 show the integrated CGM + ISM column density for C iv as a function of velocity bins of width 100 km s⁻¹ in a fashion similar to the observations. The results are shown for 27 synthesized models with different combinations of CGM mass (as traced by C iv), outflow launch velocity (v_out), and outflow fraction (f_out). The red dots indicate the integrated column density for each simulated GRB sight line. The blue line shows the median profile, and the purple shades indicate the central 50th and 80th percentiles. The vertical dotted lines separate the central region dominated by the ISM from the blue and red wings dominated by the CGM.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We compare the blue and red wings with the observations because the central region is often saturated in the observations. As clearly seen in these figures, a higher outflow fraction increases the kinematic asymmetry, while a higher outflow velocity shifts the asymmetry blueward, as expected. A higher CGM mass proportionally increases the median column densities in the red and blue wings, while the increment in the central region is slower as it is dominated by the ISM.

Despite its simplicity, this kinematic model does an excellent job at reproducing the median kinematics profile as well as typical detection fractions and column density distribution as a function of velocity. These simulations rule out an outflow fraction of ≳0.5 and a CGM mass of ≳10^10.3 M_⊙ by simple comparisons of the median profiles and the percentile distributions. However, it should be noted that a lower metallicity for the CGM would favor a higher CGM mass to explain the observations, while a higher ionization fraction would favor a lower CGM mass. Nonetheless, the current assumptions are reasonable within a factor of 2. The model that best explains the observed kinematics is M_CGM = 10^9.8 M_⊙, v_out = 300 km s⁻¹, and f_out = 0.25. Note that v_out is the launch velocity at 2 kpc for a ballistic outflow. In Section 6, we use these kinematic simulations to estimate the contribution of CGM in the central ±100 km s⁻¹ velocity region, which is crucial in estimating the typical CGM mass in our sample.

No extra inflow component has been added in these kinematic simulations. Various feedback simulations suggest that at z ∼ 2.5, cold inflows are hardly detectable in metal lines, due to their low metallicity and low covering factors (Fumagalli et al. 2011; Goerdt et al. 2012). Given the observational challenges in detecting the inflows, we do not include an inflow component in our toy model. However, a plausible inflow scenario is explored in Section 7.3.

Our kinematic toy models provide a reasonable handle on the CGM mass (M_CGM ∼ 10^9.8 M_⊙). In this section, we provide further insights into the CGM mass from the observed column densities, adopting realistic geometrical assumptions, and compare the results with the toy models. For this, we assume the CGM–GRB sample as a set of random GRB sight lines probing the typical ISM and CGM of a GRB host at z ∼ 2–5.

6.1.1. Method

A key uncertainty in this formulation is the contribution of the CGM to the column density in the central region ($| v| \,\lt 100$ km s⁻¹). To resolve the CGM–ISM degeneracy in this region, we make certain assumptions based on the insights we gained from the toy models. Because the central region also suffers from a saturation issue leading to the measurement of only the lower limit of the column density, we avoid making direct use of the CGM-to-ISM ratio in this region. Instead, we use the column density in the −150 km s⁻¹ bin (i.e., −200 km s⁻¹ < v < −100 km s⁻¹) to extrapolate the central column density because the column density in the −150 km s⁻¹ velocity bin is measured more accurately. This ratio is in the range of 2–4 for the models that were found viable in the previous section. Hence, we use a factor of 3 to estimate the CGM-contributed column density in the observed spectra. To calculate the CGM mass, we will make use of C iv column density because it is ubiquitous and a good tracer of the outflow component.

To convert the integrated CGM column density of C iv, as approximated above, into the mass of C iv in the CGM, we assume that the CGM density profile is a power law given by

$$\begin{eqnarray}&&{n}_{\mathrm{CGM}}(r)={n}_{0}{\left(\displaystyle \frac{r}{{R}_{0}}\right)}^{-\alpha }.\end{eqnarray} \tag{ 2 }$$

For convenience, we select the reference radius R₀ as the starting point of the CGM, equivalent to twice the half-light radius (R_e) of the galaxy. We assume a typical R_e of 2 kpc for GRB hosts at z > 2, following previous observations of GRB hosts (Wainwright et al. 2007; Blanchard et al. 2016), and thus R₀ = 4 kpc. Both the detailed simulations and our toy model indicate that the gas that leads to the observed absorption in C iv is spread out to about 2 × R_vir. Beyond that, the gas density is too low to give rise to a detectable absorption component (assuming n(r) ∝ r⁻²). Given a constant volume-filling fraction of f_vol, the mass of the CGM (for α $\ne$ 3) can be stated as

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{CGM}} & = & 4\pi {m}_{{\rm{H}}}{f}_{\mathrm{vol}}{\int }_{{R}_{0}}^{{R}_{\mathrm{CGM}}}n(r)\times {r}^{2}{dr}\\ & = & 4\pi {m}_{{\rm{H}}}{f}_{\mathrm{vol}}{n}_{0}{R}_{0}^{\alpha }\left(\displaystyle \frac{1}{3-\alpha }\right)\left[{R}_{\mathrm{CGM}}^{3-\alpha }-{R}_{0}^{3-\alpha }\right].\end{array}\end{eqnarray} \tag{ 3 }$$

We further define the line-covering fraction f_line as the typical fraction of a sight line that passes through a CGM cloud. The typical line fraction does not strongly vary with radius if the volume-filling fraction is constant. Given a constant f_line, the typical column density (for α $\ne$ 1) can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\mathrm{obs}} & = & {f}_{\mathrm{line}}{\int }_{{R}_{0}}^{{R}_{\mathrm{CGM}}}n(r){dr}\\ & = & {f}_{\mathrm{line}}{n}_{0}{R}_{0}^{\alpha }\left(\displaystyle \frac{1}{1-\alpha }\right)\left[{R}_{0}^{1-\alpha }-{R}_{\mathrm{CGM}}^{1-\alpha }\right].\end{array}\end{eqnarray} \tag{ 4 }$$

Combining the two equations, we get

$$\begin{eqnarray*}\begin{array}{rcl}{M}_{\mathrm{CGM}} & = & 4\pi {m}_{{\rm{H}}}{N}_{\mathrm{obs}}\left(\displaystyle \frac{{f}_{\mathrm{vol}}}{{f}_{\mathrm{line}}}\right)\\ & & \times \ \left(\displaystyle \frac{1-\alpha }{3-\alpha }\right)\left[\displaystyle \frac{{R}_{\mathrm{CGM}}^{3-\alpha }-{R}_{0}^{3-\alpha }}{{R}_{0}^{1-\alpha }-{R}_{\mathrm{CGM}}^{1-\alpha }}\right].\end{array}\end{eqnarray*}$$

In the following, we take α = 2, i.e., we consider the case where a significant fraction of the metal-rich CGM takes part in a mass-conserving outflow (Steidel et al. 2010; Pallottini et al. 2014; Chisholm et al. 2017). Because R₀ ≪ R_CGM, the expressions can be simplified. Further, we evaluate the typical value of f_vol/f_line by simulating various volume-filling fractions in a spherical shell and measuring the distribution of the line-covering fraction for random sight lines. For small filling fractions (∼0.05–0.25; e.g., as derived in Werk et al. 2014), the typical value of the f_vol/f_line ratio is found to be 1.2. Using this value, we estimate the mass of C iv in the CGM (in solar mass units) as

$$\begin{eqnarray}&&{M}_{\mathrm{CGM},{\rm{C}}{\rm{IV}}}=1.2\times {m}_{{\rm{C}}}{N}_{{\rm{C}}{\rm{IV}}}{R}_{0}{R}_{\mathrm{CGM}}\times {10}^{-13},\end{eqnarray} \tag{ 5 }$$

where m_C is the atomic mass number of carbon, R₀ and R_CGM are in kiloparsecs, and M_CGM is in solar masses.

6.1.2. Results

As stated earlier, we estimate R₀ to be 4 kpc and R_CGM to be 2 × R_vir ∼ 100 kpc. We choose the geometric mean of all the sight lines in the sample as a representative value of N_{C iv} because it best captures the large range in column density as a function of velocity whereas an arithmetic mean tends to be skewed to high values by the small number of large column densities. Similar mass estimates were performed for other species. The CGM mass estimates for various species as a function of velocity are shown in Figure 4. To evaluate carbon mass in the CGM, the fraction of carbon in the C iv phase (f_{C iv}), which depends on the temperature, density, and ionization process, needs to be constrained. We assume a conservative maximum f_{C iv} of 0.3 (see Figure 7 in Bordoloi et al. 2014), which gives

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{CGM},{\rm{C}}} & \gtrsim & 4.8\times {10}^{4}{M}_{\odot }\left(\displaystyle \frac{{N}_{{\rm{C}}{\rm{IV}}}}{{10}^{14}\,{\mathrm{cm}}^{-2}}\right)\\ & & \times \ \left(\displaystyle \frac{{R}_{0}{R}_{\mathrm{CGM}}}{100\,{\mathrm{kpc}}^{2}}\right)\left(\displaystyle \frac{0.3}{{f}_{{\rm{C}}{\rm{IV}}}}\right).\end{array}\end{eqnarray} \tag{ 6 }$$

Based on this formulation, the conservative lower limit on the carbon mass in the CGM of GRB hosts in our sample is ∼1.5 × 10⁶ M_⊙. The carbon mass in the CGM can be further extrapolated to derive the total mass of the CGM in the phase traced by C iv (T ∼ 4.5–5.5 × 10⁵ K). For this, we assume that the metallicity in the CGM is roughly 0.25 solar. This assumption is informed by the detailed simulations of CGM at z > 2 (Shen et al. 2013) and the low-z CGM metallicity estimates (e.g., Prochaska et al. 2017). This is further supported by the observed gas-phase metallicity of ∼0.1–0.2 solar for the ISM of the GRB hosts (Krühler et al. 2015; Arabsalmani et al. 2018). Due to the metal-enriched outflows, the typical CGM metallicity tends to be ∼1.5 times the ISM metallicity (Muratov et al. 2017). Thus, the total mass of the CGM gas traced by C iv can be expressed as

$$\begin{eqnarray}&&{M}_{\mathrm{CGM},\mathrm{total}}\sim {10}^{9}\,{M}_{\odot }\times \left(\displaystyle \frac{{M}_{\mathrm{CGM},{\rm{C}}}}{{10}^{6}\,{M}_{\odot }}\right)\left(\displaystyle \frac{0.25}{{Z}_{\mathrm{CGM}}}\right).\end{eqnarray} \tag{ 7 }$$

Thus, the typical mass of the CGM gas traced by C iv in this sample is ∼10^9.2 M_⊙, which is comparable to the typical stellar mass of the GRB hosts at these redshifts (z ∼ 2–4). This implies that the CGM is a very significant reservoir of baryons and metals in the galactic ecosystem at high redshifts. Thus, from a galaxy evolution standpoint, the CGM appears to be already in place at z ∼ 2–4. Despite various uncertainties in the assumed parameters, we can say with high significance that the mass of the CGM in GRB hosts is at least as much as the mass that resides in the stars, and it can be higher by as much as ∼0.3 dex if the conservative assumptions are relaxed.

The CGM mass estimates from the toy models (Section 5.2) and column density profiles are complementary in nature, strengthening our CGM mass estimate of ∼10^9.2–10^9.8 M_⊙. It should be noted that the difference between the CGM mass of the optimal toy model (10^9.8 M_⊙) and the CGM mass estimated here (∼10^9.2 M_⊙) arises for two key reasons: (a) the conservative estimate of the CGM contribution to the central ±100 km s⁻¹ and (b) the use of the geometrical mean of the column densities in the CGM mass measurement instead of the arithmetic mean.

6.1.3. Caveats

To study the evolution of the CGM mass with redshift, we divide the CGM–GRB sample into two roughly equal time bins of 1 Gyr: group 1 (z1 ∼ 2–2.7, midpoint: z = 2.3) and group 2 (z2 ∼ 2.7–5, midpoint: z = 3.6). The number of GRBs in these two bins is also nearly equal. The integrated column density as a function of velocity for high- and low-ion species are plotted in Figures 10 and 11, in the same way as Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It can be seen in these figures that the high-ion kinematics in both redshift groups are quite similar with respect to the blue asymmetry. In the case of Si ii, the blue asymmetry is weaker in the lower redshift bin (Figure 10). There are two possible reasons for this lack of asymmetry in Si ii: (a) the ionization level of the outflows is different at low z leading to a higher [Si iv/Si ii] ratio relative to the high-z bin, or (b) there is more inflowing gas in the low-z bin compared to the high-z bin, which produces a relatively stronger Si ii red wing. From the high-z simulations in Shen et al. (2013), it is clear that Si ii is a much stronger tracer of the inflows compared to outflows by almost an order of magnitude. However, a more rigorous treatment using ionization modeling is required to distinguish between the two scenarios and constrain the physical state of the outflows and inflows.

We follow the same procedure as the one described in Section 6.1 for estimating the mass of the CGM in these two redshift bins. The key changes in the assumed parameters are (a) the typical value of R_vir and (b) the typical value of R_e in these two redshift bins. The other parameters are not expected to change significantly. The virial radius is calculated as the radius within which the normalized density ρ/ρ_cosmic > 200, using the Navarro–Frenk–White (NFW) profile for dark matter distribution and standard cosmological parameters (Ω_m = 0.3, Ω_rad = 0, Ω_Λ = 0.7). Thus, the typical virial radii for the z1 and z2 redshift bins are 53 and 39 kpc, respectively. Therefore, the R_CGM for the two bins is 106 and 78 kpc. The typical half-light radius for the two bins are assumed to be 2 and 1.5 kpc, respectively, following the previous population study of the GRB hosts at these redshifts (Wainwright et al. 2007; Blanchard et al. 2016). Hence, the typical R₀ for the two bins is taken to be 4 and 3 kpc, respectively. With this setup, the masses of various species were calculated in the two redshift bins and are shown in Figures 10 and 11. The C iv in the CGM was estimated to be $\mathrm{log}({M}_{{\rm{C}}{\rm{IV}},{\rm{z}}1}/{M}_{\odot })={5.6}_{-0.2}^{+0.1}$ and $\mathrm{log}({M}_{{\rm{C}}{\rm{IV}},{\rm{z}}2}/{M}_{\odot })={5.1}_{-0.1}^{+0.2}$. The total CGM masses in the two redshift bins (following the same procedure as described in Section 6.1) are estimated to be M_CGM,z1 = 10^9.2 M_⊙ and M_CGM,z2 = 10^8.7 M_⊙.

The C iv mass in the CGM, as shown in Figures 10 and 11, is higher in the lower redshift bin by almost 0.5 dex. However, it is quite likely due to the combined effect of redshift evolution and the difference in typical stellar mass and star formation rate between the two redshifts. This will be explored in more depth by associating the CGM properties with the host properties in an upcoming paper on this sample.

To convert the C iv mass to carbon mass, we make a conservative assumption for f_{C iv}. Despite the extragalactic ionizing UV background flux evolution (within a factor of 2) between the two redshift bins, f_{C iv} does not exceed 0.3 for both photo- and collisional ionization models (Oppenheimer & Schaye 2013). Hence, we consider the maximum f_{C iv} to be 0.3 for both redshift bins (same as in Section 6.1). The M_CGM,C is thus estimated using Equation (6) for both redshift bins and plotted in Figure 12 to compare with the M_CGM,C estimates from the COS-Dwarfs study (at z ∼ 0). The typical GRB host stellar masses in the redshift bins shown in the figure are median values from the SHOALS (Swift GRB Host Galaxy Legacy Survey) sample, which is the largest systematic survey of long GRB hosts (Perley et al. 2016). The stellar mass uncertainties are the bootstrapping 1σ intervals from the SHOALS survey. As described in Section 6.1, our estimates for M_CGM,C are evaluated within 2 virial radii. Hence, for a better picture of redshift evolution, the COS-Dwarfs estimates within ∼1 virial radius are shown in green and the estimates within ∼2 virial radii in light green.

Zoom In Zoom Out Reset image size

From Figure 12, it can be seen that the CGM carbon mass (or by extension, metal mass) increased by only a factor of ∼2 (comparing with the 2R_vir estimates) from z > 2 to z ∼ 0 for dwarf galaxies. Despite the fact that f_{C iv} is assumed to be a conservative upper limit for all of these calculations, it can be clearly observed that most of the metal mass in the CGM of the low-mass galaxies represented by GRB hosts is already in place by z ∼ 2.5. While the COS-Dwarf galaxies are not the descendants of galaxies represented by GRB hosts at z > 2, this comparison has significant implications on the distribution of metals throughout the galaxy ecosystem as a function of redshift, as discussed in Section 7.2.

By comparing the carbon mass in the ISM with that in the CGM, we can infer the level of CGM enrichment for GRB hosts at high redshifts. We estimate the gas-phase carbon mass in the ISM of the GRB hosts by assuming a modest gas-phase metallicity of 0.15 solar (Davé et al. 2017; Arabsalmani et al. 2018) with a solar-like relative abundance pattern. Here, we assume a metallicity slightly lower than the median value reported in Arabsalmani et al. (2018) because the UV/optical absorption- and emission-based metallicity estimates do not account for the potentially lower metallicity gas in the ISM. The total ISM gas mass is estimated as 2 × M_* based on various molecular and neutral hydrogen measurements and simulations of high-redshift main-sequence galaxies, where f_gas = M_gas/M_* is reported to be ∼2 for 10⁹ < M_*/M_⊙ < 10¹⁰ (Daddi et al. 2010; Carilli & Walter 2013; Lagos et al. 2014; Davé et al. 2017; Tacconi et al. 2018). We estimate an uncertainty of roughly ±0.3 dex in the product of metallicity and the gas fraction (f_gas) to account for the uncertainty associated with the mass–metallicity relation for GRB hosts at z > 2 and the variation of f_gas with redshift. Note that here we are considering total gas mass of the ISM to estimate the carbon mass in the ISM (and not just the warm phase of the ISM). With these reasonable assumptions, we plot the ISM mass of carbon in Figure 12.

The carbon mass in the ISM is higher than the minimum carbon mass in the CGM by a factor of ∼3–5. It should be noted again that the carbon mass in the CGM is a lower limit and considers only the warm phase traced by C iv. In principle, the actual value of total carbon in the CGM can be higher by a factor of as much as ∼2–3 (see Section 6.1). This observation indicates that the carbon (metal) content of the CGM is ∼20%–50% of the carbon (metal) content of the ISM, indicating that a significant fraction of metals synthesized in the galaxy are able to escape into the CGM, due to the galactic outflows. In Bordoloi et al. (2014), this fraction is 50%–80% for z ∼ 0 galaxies of similar stellar mass. This finding hints toward a modest evolution in the carbon (metal) content, or in other words, a modest enrichment of the CGM over ∼11 Gyr span.

In Section 6.2, the CGM–GRB sample is divided into two redshift bins each spanning 1 Gyr (bin 1: z ∼ 2–2.7, bin 2: z ∼ 2.7–5). Outflows are clearly detected in both redshift bins as shown in Figures 10 and 11. As described in Section 4.4, the outflows are inferred from the blue asymmetry that is observed in both the high- and low-ion lines, but they are primarily traced by the high ions (C iv and Si iv). Using the kinematic information from the observations and further clarity from the toy models, we can estimate various aspects of the outflows. From the toy models, the optimal model for explaining the observed kinematics has a typical outflow launch speed of ∼300 km s⁻¹ (at 2 kpc), and ∼25% of the CGM clouds in the model contribute to outflows in excess of the virial motion in the CGM of a typical dark matter halo mass M_h ∼ 10^11.2 M_⊙. With this insight, we estimate the key properties of the outflowing material in this section. We assume an NFW profile for the dark matter distribution with a concentration parameter of 4.5 (derived by forcing R₂₀₀ = R_vir and defining the mass enclosed within a virial radius as the halo mass).

7.1.1. Outflow Mass

7.1.2. Mass Outflow Rate

We can now evaluate the mass outflow rate in these two redshift bins. The outflow rate out of a spherical shell is given as

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{\mathrm{out}} & = & \bar{m}n({R}_{\mathrm{shell}})\times 4\pi {R}_{\mathrm{shell}}^{2}\times {v}_{\mathrm{out}}({R}_{\mathrm{shell}})\\ & & \times \ {f}_{\mathrm{vol}}\times {f}_{\mathrm{out}},\end{array}\end{eqnarray} \tag{ 8 }$$

where the average mass per atom or ion is $\bar{m}$ = 1.15 m_H. We consider R_shell = 4 kpc (3 kpc) for the redshift bin z1 (z2) to estimate the outflows coming out of the galaxy disk. Note that we assume 4 kpc (3 kpc) as the boundary of the galaxy disk, R₀, in Section 6.2. Considering M_CGM = 10^9.2 M_⊙ (10^8.7 M_⊙), we get n(R_shell) = 0.028 (0.02) cm⁻³ (following Equation (3)). The outflow launch velocity in the 2–4 kpc range is approximated as ∼300 km s⁻¹ from the optimal toy model and the observed kinematics for both redshift bins. Following the optimal toy model, we assume a volume-filling fraction f_vol = 0.1 and an outflow fraction f_out = 0.25. Within this framework, we estimate the mass outflow rates (in the gas phase traced by C iv) as

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{\mathrm{out}} & = & 0.27{M}_{\odot }{\mathrm{yr}}^{-1}\left(\displaystyle \frac{n(R)}{{10}^{-2}\ {\mathrm{cm}}^{-3}}\right)\\ & & \times \ \left(\displaystyle \frac{{R}^{2}}{10\,{\mathrm{kpc}}^{2}}\right)\left(\displaystyle \frac{{v}_{\mathrm{out}}(R)}{300\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right),\end{array}\end{eqnarray} \tag{ 9 }$$

where R is the radius of the shell, R_shell. Inserting the appropriate values for the two redshift bins, we get ${\dot{M}}_{\mathrm{out},z1}$ = 1.2 M_⊙ yr⁻¹ and ${\dot{M}}_{\mathrm{out},z2}$ = 0.5 M_⊙ yr⁻¹. This warm-phase mass outflow rate is comparable to the latest TNG simulation for galaxies with M_* ∼ 10⁹–10¹⁰ M_⊙ at z = 2 (area under the curve for warm phase in Figure 10 of Nelson et al. 2019). Note that the gas outflow rate in Nelson et al. (2019) is measured at 20 kpc, however. In our density profile assumption, n(R) × R² is independent of R and the cloud velocities at 20 kpc reduce to ∼70% of their launch velocities for v_launch = 200–300 km s⁻¹. Hence, for comparison, ${\dot{M}}_{\mathrm{out},20\mathrm{kpc}}$ ∼ $0.7{\dot{M}}_{\mathrm{out},\mathrm{launch}}$ in our framework.

As a corollary of the escaping mass (M_esc) calculation, we can now estimate the mass escape rate by modifying Equation (8). For this, we consider the boundary of the CGM as R_CGM = 2R_vir. At this radius, a ballistic outflow launched at 2 kpc with a velocity of 250 km s⁻¹ decelerates to 50 km s⁻¹ (from the NFW dark matter profile). Due to our assumption of mass-conserving outflow (i.e., n(R) ∝ R⁻²), the product n(R) × R² is invariant. The modified equation for the mass escape rate becomes

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{\mathrm{esc}} & = & \bar{m}n({R}_{\mathrm{esc}})\times 4\pi {R}_{\mathrm{esc}}^{2}\times {v}_{\mathrm{out}}({R}_{\mathrm{esc}})\\ & & \times \ {f}_{\mathrm{vol}}{f}_{\mathrm{out}}{f}_{\mathrm{esc}}.\end{array}\end{eqnarray} \tag{ 10 }$$

Using 50 km s⁻¹ for v_out(R_esc) and 80% for f_esc (the fraction of outflowing gas with a launch velocity >250 km s⁻¹ from Section 7.1.1), we get ${\dot{M}}_{\mathrm{esc}}\sim 0.14\times {\dot{M}}_{\mathrm{out}}$. Thus, the estimates for the mass escape rate into the IGM for the two redshift bins are ${\dot{M}}_{\mathrm{esc},z1}$ = 0.17 M_⊙ yr⁻¹ and ${\dot{M}}_{\mathrm{esc},z2}$ = 0.07 M_⊙ yr⁻¹. These estimates are robust within a factor of 2 over ±30 km s⁻¹ variations in the launch velocity (v_out) in the optimal toy model.

The calculations described in this section suggest that ∼80% (f_esc) of the mass ejected from a GRB host galaxy of median halo mass (M_h ∼ 10^11.2 M_⊙) escapes to enrich the IGM while only 20% (f_retain) stays within the CGM of a typical GRB host in this sample at z > 2. Combining this result with the 0.5 Gyr timescale to traverse 2 × R_vir, one can say that ∼20% (f_out × 80%) of the total CGM mass of a typical GRB host at z ∼ 2–5 escapes out to the IGM over 0.5 Gyr. Conversely, the CGM mass in the warm phase grows by ∼5% (f_out × 20%) every 0.5 Gyr if we assume that the nonescaping gas becomes a part of the virialized CGM. Such a growth rate, if steady, is sufficient to significantly grow the warm-phase CGM over a 5 Gyr period.

7.1.3. Mass-loading Factor and Energetics

The mass-loading factor is defined as $\eta ={\dot{M}}_{\mathrm{out}}/\mathrm{SFR}$ and serves as an important indicator of the mechanism driving the outflow. We consider a median SFR of 10 M_⊙ yr⁻¹ for GRB hosts from Krühler et al. (2015) at z > 2 (typical range 5–50 M_⊙ yr⁻¹). The mass-loading factors at the launch radii can then be estimated as η_z1 = 1.2/10 = 0.12 and η_z2 = 0.5/10 = 0.05. Several simulations calculate the mass-loading factors at some intermediate radius such as 20 kpc. The mass-loading factors at 20 kpc can be evaluated by using decelerated velocities as (see Section 7.1.2) η_z1,20kpc = 0.7 × 0.12 = 0.084 and η_z2,20kpc = 0.7 × 0.05 = 0.035.

While the outflow velocities in mass ranges comparable to our sample are always in the range of 200–400 km s⁻¹, there is at least an order of magnitude variation in the mass-loading factors reported in various observational studies at high redshifts, due to the diversity in probes, underlying assumptions, and the phases traced in the outflow. Therefore, only an order-of-magnitude comparison can be done.

Crighton et al. (2014) used a QSO sight-line probe at z = 2.5 for an M_* ∼ 10^9.1 M_⊙ galaxy and inferred a mass-loading factor of order 1. Weiner et al. (2009) also inferred an η of order unity at launch for cold-phase outflow (tracer: Mg ii) in star-forming galaxies at z ∼ 1.4 using rest-frame UV spectra of the galaxies. Along similar lines, Martin et al. (2012) reported an η of order 1 for Fe ii-traced outflows in a redshift range z ∼ 0.4–1.4 for star-forming galaxies over a wide range of stellar mass (10^9.5–10^11.5 M_⊙). Using a similar method, Rubin et al. (2014) conservatively estimated a cold gas mass-loading factor η_5kpc ≳ 0.02–0.6 for galaxies with M_* ≳ 10^9.6 M_⊙ and SFR ≳ 2 M_⊙ yr⁻¹ in the redshift range z ∼ 0.3–1.4. Davies et al. (2019) used IFU spectra of star-forming galaxies with z ∼ 2–2.6 and M_* ∼ 10^9.5–10^11.5 M_⊙ to estimate the mass-loading factor η ∼ 0.05–0.5 at launch. This sampling of the literature shows that overall, for galaxies of comparable mass to our sample, the mass-loading factor estimates range from 0.05 to 5. It is likely that the higher mass loading in the down-the-barrel observations is due to the LOS effects (down the barrel versus random sight lines of GRBs). This highlights the need for a multiprobe approach to trace the outflow process in various phases and for various orientations of a galaxy to capture the full picture of the CGM outflows.

The mass-loading factors estimated here are smaller than the estimates from various galaxy evolution and zoom-in simulations by one to two orders of magnitude. For example, cosmological zoom simulations such as those by Anglés-Alcázar et al. (2014) and Shen et al. (2012) suggest η ∼ 1 at z > 2 to reproduce the morphological and dynamical properties of galaxies with M_h ∼ 10¹¹–10¹² M_⊙ at z ∼ 2. However, it should be noted that these mass-loading factors encompass all of the phases, and not just the warm phase traced by C iv or Si iv. The latest TNG simulations appear to resolve this issue by separating the phases in the outflows (Nelson et al. 2019). The total (all-phase) mass-loading factor η_tot,20kpc ∼ 4 (Figure 5 in Nelson et al. 2019) for a main-sequence galaxy of M_* ∼ 10^9.3 M_⊙ at z = 2 whereas the loading factor for the warm phase is η_warm,20kpc ∼ 0.15 (Figure 10 in Nelson et al. 2019). This warm-phase mass-loading factor is within a factor of 2 of the mass-loading factor evaluated here. This also implies that a significant fraction of the outflowing mass is in other phases. It will be of further interest to carry out such comparisons at higher redshifts to synthesize a complete picture of the impact of outflows on galaxy evolution and vice versa.

It should be noted that here we have assumed a ballistic outflow driven by star formation and no halo drag or outflow acceleration (e.g., by cosmic rays, ram pressure, or radiation pressure) is considered (Murray et al. 2011; Girichidis et al. 2016; Hayward & Hopkins 2017). The comparative effects of the two can be nontrivial and need to be explored in the future. At the same time, the observational mass-loading factors have several uncertainties that need to be considered while comparing them to numerical simulations. The current analysis has been conducted for an ensemble, so only typical values are considered here. We will do this for individual hosts in an upcoming paper by connecting the CGM and galaxy properties. While comparing the outflow characteristics with simulations and other surveys, the ionization fraction, outflow fraction in various phases, and the dynamics of an outflow (various accelerations and drags) need to be treated more carefully. In addition, the relative fraction of the reaccretion of the enriched gas versus its virialization in the CGM needs to be explored in the simulations.

The modest mass outflow rates estimated here can be entirely supernova driven. The kinetic energy of the outflow can be estimated as $\tfrac{1}{2}\dot{{M}_{\mathrm{out}}}{{v}^{2}}_{\mathrm{out}}\sim 5.7\times {10}^{40}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. Following the formalism described in Murray et al. (2005; see Equations (34) and 35), we assume a Salpeter IMF and that ξ ∼ 10% of the supernova energy (E_SN ∼ 10⁵¹ ergs) is efficiently thermalized in the ISM, thus driving the displacement of the gas. With a supernova rate (f_SN) of 10⁻² per M_⊙ of star formation, the energy deposition from supernovae becomes

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{SN}}\sim \xi {f}_{\mathrm{SN}}{E}_{\mathrm{SN}}\times \mathrm{SFR}\,\sim \,3\times {10}^{40}\displaystyle \frac{\xi }{0.1}\displaystyle \frac{\mathrm{SFR}}{1{M}_{\odot }\,{\mathrm{yr}}^{-1}}\,\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 11 }$$

Thus, SFR ∼ 2 M_⊙ yr⁻¹ is energetically sufficient to drive the observed warm outflows. While a comparable mass outflow rate may be hidden in other gas phases (not traced by C iv), given the typical SFR ∼ 10 M_⊙ yr⁻¹ for z > 2 GRB hosts, it is very likely that the outflows are predominantly driven by supernova energy injection.

From the standpoint of galaxy evolution, it is interesting to compare the evolution of the CGM with the evolution of the host galaxy. We combined published data from the literature (Steidel et al. 2010; Borthakur et al. 2013; COS-Dwarfs: Bordoloi et al. 2014; KODIAQ: Lehner et al. 2015; Burchett et al. 2016; COS-burst: Heckman et al. 2017; Rudie et al. 2019) with our new observations to investigate the evolution of C iv mass in the CGM relative to the stellar mass $\left({M}_{{\rm{C}}{\rm{IV}}}/{M}_{\ast }\right)$, as a function of redshift. For uniformity, we only considered C iv column density within respective R_vir (estimated using M_*–M_halo relation; Wechsler & Tinker 2018) and converted that into a C iv mass estimate using Equation (4) in Bordoloi et al. (2014). The results are summarized in Figure 13.

Zoom In Zoom Out Reset image size

It can be seen that log (M_{C iv}/M_*) ∼ −4.5 throughout the evolution of current-day galaxies with M_* > 10¹⁰ M_⊙. This relatively uniform ratio within ±0.5 dex indicates that typically the CGM coevolves with the stellar mass of the galaxy all the way from z ∼ 4 to z ∼ 0. It is possible to explain this coevolution with the cosmic star formation history. The star formation builds the stellar mass of the galaxy as well as the warm, ionized CGM traced by C iv. This important aspect of the CGM–galaxy coevolution needs to be explored further in large-scale galaxy evolution simulations.

Another possible explanation for the coevolution of warm ionized CGM and the galaxy stellar mass can be the scaling of the virial radius with redshift and galaxy mass. This is especially important because the virial radius is an arbitrary choice to define the CGM. Chen (2012) studied the radial distribution of L_* galaxies from z ∼ 2 to z ∼ 0 with a similar halo mass (M_h ∼ 10¹² M_⊙) in that redshift range. An alternative definition of the extent of the CGM was used in this case to reveal that the radial distribution of column density (normalized to the "CGM radius") did not evolve substantially from z ∼ 2 to z ∼ 0. So, it is very likely that the radial distribution of column density (normalized to the corresponding virial radius) remains fairly constant, and the evolution of the CGM mass is observed, due to the scaling of virial radius with redshift. From Figure 13, it could also be inferred that any major mergers leading to an increased halo and stellar mass grow the CGM metal mass in the same proportion as M_*. Further theoretical investigation of the CGM–galaxy coevolution is necessary to corroborate these growth pathways.

For current-day low-mass galaxies (M_* < 10¹⁰ M_⊙), there is a seemingly larger scatter in log (M_{C iv}/M_*) at z = 0, going from −3.7 to −5.0. The COS-Dwarfs sample (Bordoloi et al. 2014; primarily made of star-forming dwarf galaxies) has a systematically higher ratio, suggesting that the CGM of low-mass galaxies could be relatively more enriched in terms of either the metallicity or baryons. On the other hand, the low-mass sample (predominantly star-forming dwarfs) in a blind survey for C iv absorbers in QSO spectra (Burchett et al. 2016) shows a systematically lower ratio. However, one should note that the C iv column densities reported in Burchett et al. (2016) refer only to high impact parameters (>0.3 R_vir), which may substantially underestimate the total C iv mass. Therefore, further observations at a range of impact parameters would be required to assess whether the M_Civ/M_* ratio is indeed higher for these low-mass galaxies.

The cold (T < 10⁴ K) phase is almost exclusively traced by low ions (Tumlinson et al. 2017). If the relative proportions of warm (traced by C iv, Si iv) and cold (traced by Fe ii, Si ii) phases were the same for the red and blue wings, we would expect to see a blue asymmetry in the low ions that resembles that of the high ions (∼0.45 dex, see Table 4), but instead the asymmetry is only ∼0.18 dex. We interpreted this relatively stronger blue asymmetry in the high ions to be due to the outflows being more enriched in high ions. While this is our favorite explanation, another plausible explanation is the existence of cold inflows, which make the red wings of low ions stronger compared to the red wings of high ions, thus leading to a lower "blue asymmetry" in low ions. In this section, we briefly explore this possibility.

Under the aforementioned scenario, the inflow contribution of low ions can be evaluated as the fraction of red-wing column density of low ions that is leading to their relatively stronger red wings compared to the high ions. From Table 4, we can estimate that N_{Si ii} ∼ 10^13.2 cm⁻² and N_{Fe ii} ∼ 10^13.4 cm⁻² are contributed by the inflows moving at v ∼ 150 km s⁻¹. This column density is consistent with the feedback simulations in Shen et al. (2013) at z ∼ 3.

Because the characteristics of inflows at high-z are not well understood observationally, we rely on existing inflow models to infer a rough estimate of the inflow mass. Assuming a pristine inflow of 0.03 solar metallicity (Fumagalli et al. 2011; Glidden et al. 2016) and an ionization fraction for Si ii of ∼50% for the cold phase (see Figure 5 in Shen et al. 2013), we derive N_H,inflow ∼ 10^19.8 cm⁻². We get a similar estimate with Fe ii. However, we note a caveat that, although we have considered a pristine metallicity of 0.03 solar, it is not clear how the circulation of metals enrich the inflows by the time they reach the galactic disk. In addition, the metal-enriched recycling flows also manifest as inflowing gas and could dominate the observational signature, due to their high metal content.

Assuming a constant average density in the accretion stream and a typical area-covering fraction of the inflow of ∼5% (mostly concentrated along the direction of rotation of the galaxy; Fumagalli et al. 2011; Goerdt et al. 2012), we can estimate the inflowing mass as

$$\begin{eqnarray}&&{M}_{\mathrm{in}}\,=\,{10}^{9}{M}_{\odot }{\left(\displaystyle \frac{{R}_{\mathrm{vir}}}{50\mathrm{kpc}}\right)}^{2}\left(\displaystyle \frac{{N}_{{\rm{H}},\mathrm{in}}}{{10}^{19.8}\,{\mathrm{cm}}^{-2}}\right)\left(\displaystyle \frac{{f}_{\mathrm{area}}}{0.05}\right)\left(\displaystyle \frac{0.03}{Z/{Z}_{\odot }}\right).\end{eqnarray} \tag{ 12 }$$

Note that we have considered R_vir as the extent of the inflow instead of 2 × R_vir because we are assuming a constant average density for the inflow stream. Thus, we get an inflowing mass of ∼10⁹ M_⊙ within R_vir, which amounts to ∼60% of the diffuse warm CGM traced by C iv that we calculated in Section 6.2. Note that this is a lower limit, because we have only considered the high-velocity inflows (v_los > 100 km s⁻¹). This is an order-of-magnitude estimate owing to uncertainties in metallicity and covering fraction.

The mass accretion flow rate can be approximated by considering the inflow motion to be free fall in the halo. Various theoretical models predict large radial streams for inflows until the inner few percent of the virial radius (e.g., Fumagalli et al. 2011). Thus, the inflow timescale can be well approximated as the free-fall timescale of the halo. From the typical stellar masses of GRBs at z > 2, we estimated a halo mass of 10^11.2 M_⊙ (Section 7.1 and Appendix B). The free-fall timescale for such a halo is 500 Myr. This gives an overall inflow rate estimate of ${\dot{M}}_{\mathrm{in}}$ = M_in/t_ff = 2 M_⊙ yr⁻¹.

How can this modest gas inflow rate drive a typical high-z GRB host SFR of the order of 10 M_⊙ yr⁻¹? It should be noted that the rate evaluated above is over a 500 Myr (free fall) timescale, while the UV-derived star formation rates are fairly instantaneous in comparison (<100 Myr). GRBs typically take place in a transient (age ∼ 10–100 Myr; see Erb et al. 2006 and Levesque et al. 2010) high-SFR phase of its host. Therefore, a lower and steady cold gas accretion may be sufficient to support the typical long-term SFR history of a GRB host galaxy. In addition, large reservoirs of cold gas (from prior accretion) seem to be already present in these low-mass galaxies, due to an order-of-magnitude lower star formation efficiency at z > 4 compared to z < 3 owing to the low metallicity (Reddy et al. 2012).

While this is not direct observational evidence of cold inflows, these calculations provide an order-of-magnitude insight into the inflow rates of the low-mass galaxies at z > 2. Constraining the weak inflows will require high-quality spectra over a much larger sample.

The cases of O vi and N v in the CGM–GRB sample are interesting because of two key reasons: they trace similar regions in the temperature–density phase space (T ∼ 10^5.5 K, n ∼ 10^−4.5 cm⁻³; Tumlinson et al. 2017), but their kinematics appear to be distinct from each other. While O vi, like C iv and Si iv, shows a strong blue asymmetry of ∼0.6 dex, there is no such asymmetry in N v. The overall N v kinematics more closely resemble the low-ion kinematics. Despite the small number of sight lines (13) with O vi in the bandpass and often lower S/N in O vi, due to lower sensitivity and the Lyα forest, the kinematic differences are significant. It is important to understand the origin of N v and its distinct nature compared to O vi in order to study the warm-hot CGM phase (T ∼ 10^5.5 K), as well the nature of photo- or collisional-ionized gas harboring N v.

In various low-z CGM surveys using QSO sight lines such as COS-Halos (Werk et al. 2016), N v is rare. In a similar survey at high z (KODIAQ; Lehner et al. 2014), N v is more prevalent, but still with a detection rate of only ∼50% and with a typical column density log (N(N v)/cm⁻²) ∼ 13.6 (it is ∼14.0 in the CGM–GRB sample). On the other hand, O vi is ubiquitous in both high- and low-z surveys. Several sight lines through the Milky Way galactic halo and disk show N v absorption (Savage et al. 1997; Fox et al. 2015; Karim et al. 2018), but only ∼10% show column densities N(N v) ≳ 10¹⁴ cm⁻².

The nature of the excess N v absorption in the GRB spectra is not well understood. Possible explanations for the N v absorption in the GRB spectra include (1) photoionization of the circumburst medium within r ∼ 10 pc (Fox et al. 2008; Prochaska et al. 2008b), (2) recombination of the promptly ionized nitrogen (all electrons stripped) to N v within 10 pc (Heintz et al. 2018), and (3) N v in the CGM (Heintz et al. 2018; Fox et al. 2008). The kinematic similarity of N v absorption with the UV-pumped fine-structure lines (such as C ii*, Si ii*), which are associated with absorbers within a few hundreds of parsecs from the GRB (Vreeswijk et al. 2013), is considered an indicator for the circumburst origin of N v.

For the GRB spectra presented here, the N v absorption within ±400 km s⁻¹ typically comes from an ensemble of kinematically distinct absorbing components. Therefore, it is difficult to imagine that the absorption is produced solely by the circumburst medium. While most of the absorption components of N v (especially the strongest components) have a counterpart that is kinematically associated with fine-structure transitions (C ii*, Si ii*; as also seen in Prochaska et al. 2008a), there are also weaker components that are not associated with an excited transition within ±30 km s⁻¹ (Figures 17–43). Therefore, we suggest that the low-$| v|$ N v absorption seen in the GRB spectra comes primarily from the highly ionized/recombining gas associated with the circumburst medium while the weaker absorption at higher $| v|$ comes primarily from the warm gas associated with the CGM. This would explain the relatively higher N v detection rate and column density along GRB sight lines.

The typical O vi column densities derived in this survey are comparable to the typical O vi column densities in the high-z QSO sight-line surveys (Lehner et al. 2014) and the CGM of local star-forming galaxies (Tumlinson et al. 2011). While O vi absorption could also have a circumburst component because the recombination rates of N v and O vi are similar (Heintz et al. 2018), the kinematics of N v are distinct. In fact, the O vi kinematics resemble the kinematics of C iv and Si iv more than those of N v. This strongly suggests that the ISM and CGM are the dominant contributors to the O vi column density, although a minor contribution from the circumburst medium cannot be ruled out. Detailed photoionization models addressing the circumburst O vi absorption are required to quantitatively ascertain this.

As shown in Prochaska et al. (2006) and Chen et al. (2007), C iv and Si iv absorption have no association with the circumburst medium up to several tens of parsecs. The ionizing photon flux from the GRB strips electrons from these species (with ionization potentials of 64.5 and 45.1 eV, respectively) and their recombination timescales are of the order of 1 yr for a typical H ii region (Chen et al. 2007). Therefore, the C iv and Si iv lines primarily trace the typical ISM and CGM of the GRB host and not the circumburst medium.

High-resolution rest-frame UV spectra of the afterglows at multiple epochs are required to probe their explosion sites and environments, constrain the ionization models, and thus understand the origin of peculiar N v absorption in the GRB-afterglow spectra.