Spatial Distribution of O vi Covering Fractions in the Simulated Circumgalactic Medium

The observational sample we adopt comprises the 53 galaxies of the Multiphase Galaxy Halos Survey built from our (Hubble Space Telescope (HST)) program (PID 13398), the HST archive, and the literature. Full details of the sample have been presented in Kacprzak et al. (2015, hereafter KMC15), Kacprzak et al. (2019); Pointon et al. (2019), and Ng et al. (2019); we briefly summarize the sample here.

All galaxy–absorber pairs have spectroscopic redshifts, which range from $0.08\leqslant z\leqslant 0.67$ with a mean redshift of $\langle z\rangle =0.29$. These galaxies all lie within an impact parameter of D = 200 kpc (projected separation of a background quasar). The galaxies are "isolated" in that there are no other similarly bright galaxies within 100 kpc or within an LOS velocity of ±500 $\mathrm{km}\,{{\rm{s}}}^{-1}$. All galaxies are imaged by HST with either the ACS, WFC3, or WFPC2 instruments. Galaxy morphological parameters, inclinations, and sky orientations were obtained using the modeling software GIM2D (Simard et al. 2002).

The ${\rm{O}}\,{\rm\small{VI}}$ absorption-line properties were measured in the HST/COS G130M and G160M spectra of the background quasars. The absorption-line profiles are presented in Nielsen et al. (2017), Pointon et al. (2019), and Ng et al. (2019). Of the 53 galaxies, 29 have detected absorption with rest-frame equivalent widths ${W}_{r}(1031)\geqslant 0.1\,\mathring{\rm A}$ and 24 have upper limits or measured absorption with ${W}_{r}(1031)\lt 0.1\,\mathring{\rm A}$. Hereafter, following the terminology of KMC15 for the computation of the covering fraction, we refer to the former as "absorbers" and the latter as "non-absorbers." We will adopt these definitions for the mock absorption measurements from the simulations.

We adopt a subsample of massive galaxies from the VELA simulations (Ceverino et al. 2014; Zolotov et al. 2015), as shown in Table 1. They use the ART code (Kravtsov et al. 1997; Kravstov 1999; Kravtsov 2003; Ceverino & Klypin 2009). ART combines dark-matter ΛCDM cosmological simulations using an N-body adaptive refinement tree (ART) code and Eulerian methods to treat hydrodynamics while employing the zoom-in technique of Klypin et al. (2001). These simulations can achieve a high-resolution region extending 1–2 Mpc in diameter around simulated galaxies and a cell spatial resolution of approximately 20 pc at z = 1.

Table 1. VELA Galaxy Sample

No.	log(Mvir/M ${}_{\odot }$)	log(M*/M ${}_{\odot }$)	Rvir(kpc)	SFR(M ${}_{\odot }$ yr−1)
22	11.8	10.7	133	1.5
26	11.6	10.4	112	1.0
27	11.6	10.3	110	0.7
28	11.3	9.9	92	0.2

Download table as:  ASCIITypeset image

The VELA simulations are 20 Mpc on a side. We use smaller postproduction boxes that are centered on target galaxies and are roughly four virial radii ($4{R}_{\mathrm{vir}}$) in diameter. The simulations have a maximum cell resolution of 17 pc, a minimum stellar particle mass of 10³ M ${}_{\odot }$, and a dark-matter particle mass of $8\times {10}^{4}$ ${M}_{\odot }$. This high resolution allows us to resolve the regime in which stellar feedback overcomes radiative cooling (Ceverino & Klypin 2009), which results in natural galactic outflows (Ceverino et al. 2010, 2016). This allows for a combination of cold flow accretion, mergers, and galactic outflows that results in galaxy formation and evolution proceeding on physically based principles.

We selected simulated galaxies with halo masses on the order of $11.3\leqslant \mathrm{log}{M}_{\mathrm{vir}}/{M}_{\odot }\leqslant 11.8$, which is in the range of the sample of observed comparison galaxies. While this may introduce a "progenitor bias" due to mass growth of galaxies and the difference in redshift between the observed and simulated sample, our goal was to study galaxies with similar ${M}_{\mathrm{vir}}/{M}_{\odot }$. We further enforced the condition that the galaxies were isolated from other large galaxies in the simulation volume and that none had experienced a major merger (defined by a mass ratio of 0.3 or greater) for at least 1 Gyr prior to $z\,\simeq \,$1. We apply these criteria because the observed comparison galaxies are isolated (see Section 2.1) and their "normal" morphologies suggest none have recently suffered a major merger.

Ideally, we would prefer to study simulated galaxies between $0.1\leqslant z\leqslant 0.7$ to match the redshift range of the observed comparison galaxy sample (with mean $z\simeq 0.3$), but massive VELA galaxies were not evolved below $z\simeq 1$. Dekel et al. (2009a) find that systematic changes in the morphology and mass flux of infalling filaments occurs at $z\lt 1$ such that the subsequent lack of "clumps" carried in these filaments allows enhanced disk stabilization and relaxation at $z\lt 1$. This introduces some concern (which we discuss further in Section 4.3), but, as described above, we have attempted to control for this effect by selecting simulated galaxies for which a major merger has not occurred for at least 1 Gyr prior to $z\simeq 1$. With regard to global CGM properties, there is no observed evolution of the projected radial profiles of absorbing gas strength (Chen 2012), which is a favorable observational result with regard to our being limited to $z\simeq 1$ simulated galaxies.

Finally, because we aim to study the covering fraction as a function of azimuthal angle with respect to the galaxy disk, we selected simulated galaxies for which we could verify a disk morphology. We define the rotation axis of the galaxy disk to be the angular momentum vector of the cold gas, defined to have $T\leqslant {10}^{4}$ K, computed in a volume of $0.1{R}_{\mathrm{vir}}$ around the galaxy center. This angular momentum vector provided the orientation of the galaxy in the simulation box. We then examined sky-projected images of each galaxy's hydrogen number density as a function of galaxy inclination to visually confirm that we could study the galaxies as a function of their "observed" inclination.

Four of the VELA-simulated galaxies fit all our selection criteria. Out of these four galaxies, the last major merger was at z = 2.25, corresponding to ∼3 Gyr prior to z = 1. We list the four simulated galaxies and their general properties in Table 1. Sky-projected images of the galaxies are shown in Figure 1 for both their face-on view, $i=0^\circ$ (LOS parallel to the angular momentum vector), and edge-on view $i=90^\circ$ (angular momentum vector parallel to the plane of the sky).

Zoom In Zoom Out Reset image size

In order to generate the mock quasar spectra and analyze the absorption features, we use the Mockspec pipeline developed by C. Churchill (Churchill et al. 2015) and R. Vander Vliet (Vander Vliet 2017). Mockspec is publicly available on GitHub. ¹¹ A unique capability of Mockspec is the ability to identify (and isolate for further study) the gas cells that give rise to significantly detectable absorption. We can therefore isolate, measure, and analyze the properties of the absorbing gas, such as density, temperature, metallicity, 3D velocities, dynamic relationship to the simulated galaxy, and 3D spatial location.

Each gas cell in ART has a physical size, ${L}_{\mathrm{cell}}$, and a unique 3D spatial coordinate. Also recorded are the cell 3D velocity components, hydrogen number density, ${n}_{{\rm\small{H}}}$, kinetic temperature, T, and metal mass fraction, ${x}_{{\rm\small{M}}}$. To obtain the number densities of all ion stages, we perform postprocessing equilibrium ionization modeling to obtain the ionization fractions. We use the photo+collisional ionization code HartRate (detailed in Churchill et al. 2014). We adopt solar abundance mass fractions for the individual metals up to zinc (Draine 2011; Asplund et al. 2009). This code works best for optically thin, low-density gas, making it well suited for studying the low-density CGM ($\mathrm{log}{n}_{{\rm\small{H}}}/{\mathrm{cm}}^{-3}\lt -1$).

HartRate has been used in previous studies (Churchill et al. 2012; Kacprzak et al. 2012; Churchill et al. 2015; Kacprzak et al. 2019). A comparison to the industry-standard ionization code Cloudy (Ferland et al. 1998, 2013) shows that the ionization fractions are in agreement within ±0.05 across astrophysically applicable densities and temperatures for optically thin gas (Churchill et al. 2014). All gas cells in postproduction boxes (roughly 4R_vir in diameter centered on the galaxy) are illuminated with the ultraviolet background (UVB) spectrum of Haardt & Madau (2005) and the equilibrium solution is obtained. For each cell, HartRate records the electron density, ionization and recombination rate coefficients, ionization fractions, and number densities of all ionic stages from hydrogen through zinc.

Mockspec generates a user-specified number of "quasar" LOSs distributed through each simulated galaxy. For a given galaxy, the LOS path is defined by its position angle on the plane of the sky, its impact parameter, and the sky-projected inclination of the galaxy. The position angle resides in the range $0^\circ \leqslant \phi \leqslant 360^\circ$, and the plane of the sky is defined as the plane through the galaxy center of mass that is perpendicular to the LOS. The impact parameter is defined as the minimum projected distance between the LOS and the galaxy center of mass, which is equivalent to the distance between the point where the LOS intersects the plane of the sky and the center of mass of the galaxy. The inclination of the galaxy is defined by the angle between the angular momentum vector and the plane of the sky.

Observationally, when we speak of the orientation of a given LOS through a galaxy halo, we refer to the galaxy inclination and the LOS azimuthal angle, which is the primary angle ($0^\circ \leqslant {\rm{\Phi }}\leqslant 90^\circ$) between the sky-projected major axis of the galaxy and the LOS. We adopt this definition for the mock LOS through the simulated galaxies, where ${\rm{\Phi }}=0^\circ$ indicates the major axis and ${\rm{\Phi }}=90^\circ$ indicates the minor axis. The conversion of the position angle ϕ to the azimuthal angle Φ is a simple rotation and conversion to the primary angle. The methodology to produce synthetic spectra from quasar sight lines has been outlined in Churchill et al. (2015) and Vander Vliet (2017). Absorption features are detected objectively using the methods described in Churchill et al. (2000), which is derived from methods described in Schneider et al. (1993).

We aim to compare the simulated sample of ${\rm{O}}\,{\rm\small{VI}}$ galaxy–absorber pairs to the observed sample of ${\rm{O}}\,{\rm\small{VI}}$ galaxy–absorber pairs studied by KMC15, focusing on the covering fraction as a function of impact parameter and azimuthal angle. For the impact parameter comparison in this study, we normalized the impact parameter by the virial radius ($D/{R}_{\mathrm{vir}}$). As we are limited to studying only four simulated galaxies, we must pragmatically adopt the philosophy that many LOSs through a few galaxies is equivalent to a single LOS through many galaxies.

This assumption is observationally supported by the facts that (1) the kinematic distribution of ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas is uniform around galaxies (Nielsen et al. 2017) and (2) the spatial bimodality signature is highly significant in a sample of as few as 29 "absorbing" galaxies and 24 "non-absorbing" galaxies, which implies any individual galaxy is likely to exhibit spatial bimodality of its ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas. Further, one of our goals is to learn whether individual simulated galaxies that are isolated and have no recent major merger history have ${\rm{O}}\,{\rm\small{VI}}$-absorbing halos that exhibit covering fraction distributions consistent with the observations. The results will provide insights into the feedback prescriptions, which are physics based.

In order to control our experiment to match observations as closely as is plausible, we create "mock samples" that emulate the number of galaxy–absorber pairs and the distributions of $D/{R}_{\mathrm{vir}}$ and orientation of the LOS in the observed sample. As described below, we account for variations in the realizations of each mock sample and estimate uncertainties in the covering fractions by averaging thousands of mock samples.

We began by creating a grid of LOSs through the four simulated galaxies. For each galaxy, we ran Mockspec at 10 inclination angles ranging from $i=0^\circ$ (face on) to $i=90^\circ$ (edge on) in steps of $10^\circ$. For each galaxy at each inclination angle, we configured Mockspec to generate 1000 LOSs with a random distribution of impact parameters between $0\leqslant D/{R}_{\mathrm{vir}}\leqslant 1.5$ and a random distribution of position angles on the sky between $0^\circ \leqslant \phi \leqslant 360^\circ$. We convert the position angles to the primary azimuthal angles, which range between ${\rm{\Phi }}=0^\circ$ (probing along the projected major axis) and ${\rm{\Phi }}=90^\circ$ (probing along the projected minor axis). This resulted in a library of 40,000 LOSs, with each LOS defined by which galaxy it probes, the galaxy inclination, the LOS impact parameter, and LOS azimuthal angle.

For this study, we chose to normalize the true impact parameter by the virial radius for two reasons. First, we are able to obtain a view of the observed and mock samples normalized to the galaxy mass and mean dark-matter halo overdensity, which is a more commonly adopted view for the theoretical perspective. Second, the observed sample has a mean of $\langle z\rangle \sim 0.3$, whereas the simulated galaxies are at $z\simeq 1$, having not been evolved to lower redshift. As the galaxies span a range of redshifts and masses, normalizing by virial radius places both the observations and simulations on a comparable scale. Given that the virial radii of galaxies grow substantially from z = 1 to $z\sim 0$ (e.g., Bryan & Norman 1998), it is beneficial to do an analysis that normalizes for galaxy mass and mean halo overdensity, as there is evidence, at least for Mg ii absorption, that the extent of the CGM is self-similar with halo mass (e.g., Churchill et al. 2013).

It is not entirely possible to normalize out the evolution in the CGM gas-ionization conditions that would track the evolution in the UVB from z = 1 to z = 0, as the mean intensity of hydrogen ionizing photons decreased by roughly an order of magnitude over this range (Haardt & Madau 2005). This can change the relative fraction of photoionized gas to collisionally ionized gas on the ${\rm{O}}\,{\rm\small{VI}}$ phase. We note that the observed sample spans the approximate redshift range $0.1\leqslant z\leqslant 0.7$, with a mean of $z\simeq 0.3$. For $z\simeq 0.3$, Oppenheimer et al. (2016) showed that the distribution of photoionized to collisionally ionized gas is halo mass dependent. In our halo mass range for the simulations (${10}^{11.3}\leqslant \mathrm{log}{M}_{h}/{M}_{\odot }\leqslant {10}^{11.8}$), they find that the ${\rm{O}}\,{\rm\small{VI}}$ arises primarily in photoionized gas.

Is it possible that this relative balance of photoionized to collisionally ionized ${\rm{O}}\,{\rm\small{VI}}$-bearing CGM gas would be significantly altered at z = 1, the epoch of our simulated galaxies? Further insights are provided by Roca-Fàbrega et al. (2019), who show that the fraction of photoionized to collisionally ionized gas has both a mass and redshift dependence (see their Figure 14); the amplitude in the trends of the ratio are no more than ${\rm{\Delta }}f(\mathrm{phot}/\mathrm{coll})\simeq \pm 0.1$ from z = 1 to z = 0 with substantial overlap between galaxies in our mass range. In summary, though it is difficult to definitively assess the effects of mass growth and UVB evolution in the VELA simulations without running them to lower redshift, the findings of Oppenheimer et al. (2016) and Roca-Fàbrega et al. (2019) are suggestive that the effect is not substantial.

For each LOS, Mockspec generates a synthetic spectrum covering the O vi λλ1031, 1037 absorption. As the observed ${\rm{O}}\,{\rm\small{VI}}$ absorption properties were measured with the HST/COS G130M and G160M gratings, we instructed Mockspec to generate absorption-line spectra with the dispersion and pixel-sampling properties of these gratings and the COS FUV CCDs. We adopted a Poisson fluctuation plus read-noise model (see Churchill et al. 2015) and adopted a signal-to-noise ratio (S/N) of ${\rm{S}}/{\rm{N}}=30$ per pixel in the continuum. This S/N ratio corresponds to a 3σ equivalent-width detection limit of ${W}_{r}\sim 0.05\,\mathring{\rm A}$. Mockspec analyzes the absorption spectra and measures the rest-frame equivalent widths, ${W}_{r}(1031)$, and various other quantities that we discard for this study. Recall that our definition of an ${\rm{O}}\,{\rm\small{VI}}$ "absorber" is a detection of ${W}_{r}(1031)\geqslant 0.1\,\mathring{\rm A}$, whereas a "non-absorber" is defined by ${W}_{r}(1031)\lt 0.1\,\mathring{\rm A}$, which includes systems not detected to the limit of the spectroscopic data. All absorption lines generated by Mockspec are unblended and no additional noise sources are added.

We draw from this library of 40,000 unique LOSs to create our mock samples. Recall that the observed sample comprises 53 galaxies having a distribution of galaxy inclinations and LOS impact parameters and azimuthal angles. Thus, a given mock sample is created by drawing 53 LOSs from the library such that the distribution of inclination angles, impact parameters, and azimuthal angles are consistent with the observed sample.

For each mock sample, we calculated the covering fraction as a function of $D/{R}_{\mathrm{vir}}$ and azimuthal angle. For azimuthal angle, we used the binning employed by KMC15. To account for variations in the covering fraction between realizations of mock samples, we generated and analyzed 15,000 mock samples and compute the mean covering fraction and the 1σ confidence intervals of the mean. This interval was determined from the cumulative distribution function (CDF) of all 15,000 realizations. The upper and lower confidence levels correspond to the 0.8413 and 0.1587 fractional areas under the CDF.

We examined properties of the gas cells that contributed to ${\rm{O}}\,{\rm\small{VI}}$ absorption lines with ${W}_{r}(1031)\geqslant 0.1\,\mathring{\rm A}$ in the mock spectra; we hereafter call these "absorbing cells." The details of how absorbing gas cells are selected are described in Churchill et al. (2015) and Vander Vliet (2017). In brief, these cells contribute more than 5% to the measured equivalent width of the ${\rm{O}}\,{\rm\small{VI}}$ $\lambda 1031$ absorption line. Furthermore, to more broadly characterize the CGM of the VELA galaxies, we compared the absorbing cell properties to the properties of all gas cells within a $\sim 1{R}_{\mathrm{vir}}$ region of the simulated CGM.

The gas mass of a cell is computed from

$$\begin{eqnarray}&&{M}_{{\rm{c}}}=\displaystyle \frac{{n}_{{\rm\small{H}}}{m}_{{\rm\small{H}}}}{{x}_{{\rm\small{H}}}}{L}_{{\rm{c}}}^{3},\end{eqnarray} \tag{ 1 }$$

where ${n}_{{\rm\small{H}}}$ is the number density of hydrogen, ${m}_{{\rm\small{H}}}$ is the mass of hydrogen, ${x}_{{\rm\small{H}}}$ is the hydrogen mass fraction, and ${L}_{{\rm{c}}}$ is the length of the cell face. Employing conservation of mass fractions, ${x}_{{\rm\small{H}}}+{x}_{\mathrm{He}}+{x}_{{\rm\small{M}}}=1$, we estimate ${x}_{{\rm\small{H}}}=(1-{x}_{{\rm\small{M}}})/(1+r)$, where $r={x}_{\mathrm{He}}/{x}_{{\rm\small{H}}}$ is the ratio of the mass fraction of helium to hydrogen. The value of r ranges from 0.3334 for primordial mass fractions to 0.3366 for solar mass fractions (Lodders 2019); this narrow range propagates to a 0.2% difference in in our estimated cell masses. We thus adopt a single value of r = 0.335 for all cells. This value is appropriate for metal mass fractions of roughly ∼0.1–0.01 solar, a value typical of the cells in the simulations.

In Figure 4(a), we show the CGM gas mass distribution on the ${n}_{{\rm\small{H}}}$–T gas-phase diagram for all of the gas enclosed within a 1R_vir volume, where we have coadded all four VELA galaxies. The temperature floor and ceiling of the simulations are $\mathrm{log}T/{\rm{K}}=2.4$ and $\mathrm{log}T/{\rm{K}}=8$, respectively. In Figure 4(b), we show the mass distribution of CGM gas selected by ${\rm{O}}\,{\rm\small{VI}}$ absorption drawn from the catalog of 40,000 sight lines. All four simulated galaxies are coadded, and gas cells in a given galaxy that were probed by multiple LOSs (a rare condition) were counted only once. We emphasize that for the absorbing gas, we do not sample the full volume of the CGM, but only the cells along the LOS.

Zoom In Zoom Out Reset image size

In general terms, for collisionally ionized ${\rm{O}}\,{\rm\small{VI}}$, we expect the gas temperature range $5\leqslant \mathrm{log}T/{\rm{K}}\leqslant 6$ with hydrogen number density in the range $-3\leqslant \mathrm{log}{n}_{{\rm\small{H}}}/{\mathrm{cm}}^{-3}\leqslant -1$, whereas for photoionized ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas, we expect the ranges $4\leqslant \mathrm{log}T/{\rm{K}}\leqslant 4.5$ and $-4.5\leqslant \mathrm{log}{n}_{{\rm\small{H}}}/{\mathrm{cm}}^{-3}\leqslant -3.5$ (e.g., Sutherland & Dopita 1993; Bergeron & Herbert-Fort 2005; Vander Vliet 2017). The majority of the mass of the gas selected by ${\rm{O}}\,{\rm\small{VI}}$ absorption is warm ($\mathrm{log}T/{\rm{K}}\simeq 4.2$), moderately low-density ($\mathrm{log}{n}_{{\rm\small{H}}}/{\mathrm{cm}}^{-3}\approx -4$) gas, which is consistent with photoionized ${\rm{O}}\,{\rm\small{VI}}$.

As computed from Equation (1), the average mass of gas in the CGM of the simulated galaxies within 1R_vir is $\mathrm{log}{M}_{\mathrm{CGM}}/{M}_{\odot }=10.5$, which is comparable to the average stellar mass of $\mathrm{log}{M}_{\star }/{M}_{\odot }\simeq 10.4$. Using Figure 4(b) as a guide, if we crudely assume that ${\rm{O}}\,{\rm\small{VI}}$ absorption arises in all CGM gas in the phase defined by $4\leqslant \mathrm{log}T/{\rm{K}}\leqslant 6.5$ and $-5\leqslant \mathrm{log}{n}_{{\rm\small{H}}}/{\mathrm{cm}}^{-3}\leqslant -1$, as this is the temperature and number density range of hydrogen which produces 98% of the ${\rm{O}}\,{\rm\small{VI}}$ mass, the average mass of gas selected by ${\rm{O}}\,{\rm\small{VI}}$ absorption is $\mathrm{log}{M}_{{\rm\small{OVI}}}/{M}_{\odot }=10.4$. This should be considered an upper limit, but it indicates that a significant amount of gas mass in the CGM is in the phase that would be detected in ${\rm{O}}\,{\rm\small{VI}}$ absorption.

Conversely, comparing to the non-absorbing CGM gas, which we define to be all gas with $\mathrm{log}T/{\rm{K}}\geqslant 4$ that resides outside the phase region defined by $4\leqslant \mathrm{log}T/{\rm{K}}\leqslant 6.5$ and $-5\leqslant \mathrm{log}{n}_{{\rm\small{H}}}/{\mathrm{cm}}^{-3}\leqslant -1$ (which we used to approximate the ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas), we find that the average volume of ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas in the CGM of a simulated galaxy is $\mathrm{log}{V}_{\mathrm{CGM}}/$kpc${}^{3}=6.0$. This is a factor of 5 smaller than the volume filled by non-absorbing gas ($\mathrm{log}{V}_{{\rm\small{OVI}}}/$kpc${}^{3}=6.7$). The statistically equivalent radius of a sphere encompassing the absorbing gas is ≃60 kpc, which corresponds to ≃55% of the average virial radius of the simulated galaxies (${\bar{R}}_{\mathrm{vir}}=112$ kpc). Alternatively, the equivalent radius for the non-absorbing gas is ≃105 kpc, which corresponds to ≃95% of the average virial radius of the simulated galaxies. The relatively small volume of ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas is consistent with an expectation of a relatively low covering fraction in the simulated VELA galaxies.

As with our mass estimates, the volume of the ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas that we estimate here should be viewed as an upper limit because we make these calculations employing a region of the gas-phase diagram defined by $4\leqslant \mathrm{log}T/{\rm{K}}\leqslant 6.5$ and $-5\leqslant \mathrm{log}{n}_{{\rm\small{H}}}/{\mathrm{cm}}^{-3}\leqslant -1$. This is done because the LOSs, being pencil beam probes, do not sample the full volume of the CGM. It is likely that not all gas in the defined phase range gives rise to detectable ${\rm{O}}\,{\rm\small{VI}}$ absorption.

A favorable combination of ionization conditions and chemical enrichment increases the likelihood that CGM gas will give rise to ${\rm{O}}\,{\rm\small{VI}}$ absorption. Again, approximating the gas-phase region and $-5\leqslant \mathrm{log}{n}_{{\rm\small{H}}}/{\mathrm{cm}}^{-3}\leqslant -1$ for ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas and all other gas on the phase diagram with $\mathrm{log}T/{\rm{K}}\gt 4$ as non-absorbing CGM gas, we find that the distribution of the ${\rm{O}}\,{\rm\small{VI}}$ ionization fraction of absorbing gas is highly peaked with a single mode at $\mathrm{log}{f}_{{\rm\small{OVI}}}=-1$ and a flat tail out to $\mathrm{log}{f}_{{\rm\small{OVI}}}=-5$. On the other hand, the distribution for non-absorbing gas is normally distributed with a single mode at $\mathrm{log}{f}_{{\rm\small{OVI}}}=-7$ and an HWHM of ±2.5 dex.

The distribution of metallicity for the absorbing gas is highly peaked, comprising a single mode at $\mathrm{log}Z/{Z}_{\odot }=-0.7$ with a high-metallicity tail up to the solar value and an extended low metallicity tail down to $\mathrm{log}Z/{Z}_{\odot }=-2$. The non-absorbing gas has a bimodal distribution, with a highly peaked and extremely narrow mode at $\mathrm{log}Z/{Z}_{\odot }=+0.1$, and a broader normally distributed peak at $\mathrm{log}Z/{Z}_{\odot }=-0.7$. The high-metallicity mode is primarily very hot and diffuse gas ($\mathrm{log}T/{\rm{K}}\geqslant 6.5$ and $\mathrm{log}{n}_{{\rm\small{H}}}/{\mathrm{cm}}^{-3}\leqslant -6$). This component of the CGM may be highly metal enriched, but it comprises less than 5% of the gas mass within 1R_vir of the galaxy. In these simulations, it is the very hottest gas that carries the markings of recent processing through stellar feedback.

As inflowing and outflowing gas through the CGM are two major components of the baryon cycle, it is of interest to investigate the relative contributions of these components to the covering fraction. In particular, insights into the rapid drop in the covering fractions at $0.375\lt D/{R}_{\mathrm{vir}}\leqslant 0.75$ in the mock sample (see Figure 2) might be connected to the relative contribution of inflowing versus outflowing ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas as a function of $D/{R}_{\mathrm{vir}}$. To quantify the sky-projected contributions of these two kinematic components, we examined the ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas mass surface density of ${\rm{O}}\,{\rm\small{VI}}$ absorption-selected gas, ${\rm{\Sigma }}({\rm{O}}\,{\rm\small{VI}}$-absorbing gas).

For a given sky-projected annular region of the impact parameter, ${D}_{\mathrm{in}}\leqslant D\leqslant {D}_{\mathrm{out}}$, having area A, the ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas mass surface density is

$$\begin{eqnarray}&&{\rm{\Sigma }}({\rm{O}}\,{\rm\small{VI}})=\displaystyle \frac{{M}_{{\rm{O}}{\rm\small{VI}}}}{A}=\displaystyle \frac{{m}_{{\rm\small{O}}}}{2\pi \left({D}_{\mathrm{out}}^{2}-{D}_{\mathrm{in}}^{2}\right)}\displaystyle \sum _{i}{n}_{i,{\rm{O}}{\rm\small{VI}}}{L}_{i,{\rm{c}}}^{3},\end{eqnarray} \tag{ 2 }$$

where ${m}_{{\rm\small{O}}}$ is the mass of oxygen. For each sky-projected annular region, the sum is over absorbing cells that reside within the virial radius, where ${n}_{i,{\rm{O}}{\rm\small{VI}}}$ is the ${\rm{O}}\,{\rm\small{VI}}$ number density in cell i and ${L}_{i,{\rm{c}}}$ is the size of cell i. To obtain the ${\rm{\Sigma }}({\rm{O}}\,{\rm\small{VI}})$ of inflowing absorbing gas, we sum over all cells with a negative radial velocity component, ${v}_{r}\lt 0$, and for outflowing absorbing gas, we sum over all cells having ${v}_{r}\geqslant 0$. The radial velocity is measured with respect to the galaxy center. We also compute ${\rm{\Sigma }}({\rm{O}}\,{\rm\small{VI}})$ for all absorbing cells so we can obtain the total ${\rm{O}}\,{\rm\small{VI}}$ mass surface density.

In Figure 5, we present the average ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas mass surface density, $\bar{{\rm{\Sigma }}}({\rm{O}}\,{\rm\small{VI}}$-absorbing gas), as a function of $D/{R}_{\mathrm{vir}}$, using the same binning shown in Figure 2. The average values are obtained by summing over all LOSs through all simulation galaxies and dividing by the number of galaxies.

Zoom In Zoom Out Reset image size

Overall, the average ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas mass surface density generally exhibits a decline as a function of $D/{R}_{\mathrm{vir}}$ that is qualitatively similar to that of the covering fraction. In fact, as one increases in $D/{R}_{\mathrm{vir}}$, the relative values of $\bar{{\rm{\Sigma }}}({\rm{O}}\,{\rm\small{VI}}$-absorbing gas) are quantitatively consistent (within uncertainties) with the relative values of the simulated covering fractions from the mock samples presented in Figure 2. The total average ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas mass surface density (black points) is roughly a factor of 4 higher at $D/{R}_{\mathrm{vir}}\leqslant 0.375$ as compared to the projected region $0.375\lt D/{R}_{\mathrm{vir}}\leqslant 0.75$, and this matches the ratio of the covering fractions for these regions. That is, in the vicinity of $D/{R}_{\mathrm{vir}}\simeq 0.375$, there is a clear dropoff in the sky-projected ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas mass surface density, and this manifests as a similarly proportioned dropoff in the covering fraction. This is in contrast with the behavior of the observed covering fraction, which remains virtually constant for $0\leqslant D/{R}_{\mathrm{vir}}\leqslant 0.75$. This suggests that, for the simulations, the covering fraction is a good indicator of the average mass surface density of absorbing gas, and vice versa.

For these simulated galaxies, outflowing ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas is the dominant contributor to the total mass surface density at all $D/{R}_{\mathrm{vir}}$. However, the ratio of outflow to inflow decreases rapidly at $D/{R}_{\mathrm{vir}}\simeq 0.375$. Beyond this sky-projected region, the ratio of the outflowing and inflowing mass surface densities remains roughly constant. For $D/{R}_{\mathrm{vir}}\leqslant 0.375$, outflowing gas comprises roughly 70% of the total ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas mass surface density, whereas inflowing gas comprises only about 30%. Taken at face value, we would infer that roughly 70% of the ${\rm{O}}\,{\rm\small{VI}}$ covering fraction in the projected region $D/{R}_{\mathrm{vir}}\leqslant 0.375$ is contributed by outflowing ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas, whereas the remaining 30% is contributed by inflowing ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas.

For $0.375\lt D/{R}_{\mathrm{vir}}\leqslant 0.75$, outflowing and inflowing gas contribute roughly 55% and 45%, respectively, to the total ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas mass surface density. We may infer that these are the approximate proportions that inflowing and outflowing ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas contribute to the covering fraction in the mock sample. Based on the mass surface densities, we find that the contribution to absorption from outflows in the region $0.375\lt D/{R}_{\mathrm{vir}}\leqslant 0.75$ have decreased by a factor of roughly 5 relative to $D/{R}_{\mathrm{vir}}\leqslant 0.375$. That is, there is 20% as much mass surface density of outflowing absorbing gas in the region $0.375\lt D/{R}_{\mathrm{vir}}\leqslant 0.75$ as there is in the region $D/{R}_{\mathrm{vir}}\leqslant 0.375$. In contrast, the contribution to absorption from inflows in the region $0.375\lt D/{R}_{\mathrm{vir}}\leqslant 0.75$ have decreased by a only factor of 2 relative to $D/{R}_{\mathrm{vir}}\leqslant 0.375$.

In summary, though there is a steady and smooth increase in $\bar{{\rm{\Sigma }}}({\rm{O}}\,{\rm\small{VI}}$-absorbing gas) from the outer to the inner halo for inflowing gas, there is a dramatic increase in $\bar{{\rm{\Sigma }}}({\rm{O}}\,{\rm\small{VI}}$-absorbing gas) in the inner halo for outflowing gas. These trends in the average ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas mass surface density with $D/{R}_{\mathrm{vir}}$ inform us that outflows dominate the ${\rm{O}}\,{\rm\small{VI}}$ absorption in the inner regions of the simulated CGM. We infer that the majority of the outflowing gas does not reach the outer regions of the CGM (only $\sim 20$% is reaching $D/{R}_{\mathrm{vir}}\gt 0.375$). We thus attribute the factor of four decrease in the ${\rm{O}}\,{\rm\small{VI}}$ absorption covering fraction at $D/{R}_{\mathrm{vir}}\simeq 0.375$ in the mock sample (see Figure 2) to primarily be a consequence of the majority of the outflowing gas being confined to the inner halo in the simulations.

To further characterize the inflow and outflow kinematics of the ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas, we examined the distribution of the ${\rm{O}}\,{\rm\small{VI}}$ mass for ${\rm{O}}\,{\rm\small{VI}}$ absorption-selected cells as a function of both the radial velocity and the direction cosine of the radial component of the velocities. For a cell with velocity components v_r , v_θ , and v_ϕ , the direction cosine is given by

$$\begin{eqnarray}&&\alpha =\cos (a)=\displaystyle \frac{{\bf{v}}\cdot {\hat{{\bf{v}}}}_{r}}{| {\bf{v}}| }=\displaystyle \frac{{v}_{r}}{\sqrt{{v}_{r}^{2}+{v}_{\phi }^{2}+{v}_{\theta }^{2}}},\end{eqnarray} \tag{ 3 }$$

where ${\hat{{\bf{v}}}}_{r}$ is the radial velocity unit vector. The range of α is $-1\leqslant \alpha \leqslant +1$, where $\alpha =-1$ indicates pure radial inflow and $\alpha =+1$ indicates pure radial outflow. The ${\rm{O}}\,{\rm\small{VI}}$ mass in each cell is computed from ${M}_{{\rm{c}}}({\rm{O}}\,{\rm\small{VI}})={n}_{{\rm\small{O}}}{m}_{{\rm{O}}{\rm\small{VI}}}{L}_{{\rm{c}}}^{3}$.

In Figures 6(a)–(c), we present the ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas mass distribution as a function of radial velocity and, in Figures 6(d)–(f), we present the distribution of the radial velocity direction cosine. The 2D distribution of ${\rm{O}}\,{\rm\small{VI}}$ mass is presented as a function of the radial velocity and the radial velocity direction cosine in Figure 6(g)–(i). The distributions are plotted for the three sky-projected regions $D/{R}_{\mathrm{vir}}\leqslant 0.375$, $0.375\lt D/{R}_{\mathrm{vir}}\leqslant 0.75$, and $0.75\lt D/{R}_{\mathrm{vir}}\leqslant 1.125,$ which correspond to the inner three regions for which the ${\rm{O}}\,{\rm\small{VI}}$ mass surface density was computed (see Figure 5). We did not compute the distributions in the outermost region ($D/{R}_{\mathrm{vir}}\gt 1.125$) as both the ${\rm{O}}\,{\rm\small{VI}}$ covering fraction and mass surface density were quite similar to the third region, and our focus was to interpret the significant drop seen in $0.375\lt D/{R}_{\mathrm{vir}}\leqslant 0.75$.

Zoom In Zoom Out Reset image size

For $D/{R}_{\mathrm{vir}}\leqslant 0.375$, we find two modes in the radial velocity distribution (Figure 6(a)). The first mode is centered on ${v}_{r}\simeq +10$ $\mathrm{km}\,{{\rm{s}}}^{-1}$ and comprises a roughly equal mixture of inflowing and outflowing ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas with a dispersion of $\sigma ({v}_{r})\simeq 100$ $\mathrm{km}\,{{\rm{s}}}^{-1}$. The second is at ${v}_{r}\simeq +180$ $\mathrm{km}\,{{\rm{s}}}^{-1}$ and comprises outflowing ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas with a dispersion of $\sigma ({v}_{r})\simeq 50$ $\mathrm{km}\,{{\rm{s}}}^{-1}$. The distribution of the radial velocity direction cosines shows a very strong peak at $\alpha \simeq +1$ (Figure 6(d)), indicating that a substantial proportion of the outflowing absorbing gas is on a pure radial trajectory. The 2D distribution (Figure 6(g)) clearly reveals a mass peak of high-velocity radially outflowing ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas at $120\leqslant {v}_{r}\,\leqslant 230$ $\mathrm{km}\,{{\rm{s}}}^{-1}$ and $\alpha \simeq +1$.

These distributions indicate that the high covering fraction in the region $D/{R}_{\mathrm{vir}}\leqslant 0.375$ arises due to (1) roughly equal mass of gas both inflowing and outflowing in the velocity range $-100\leqslant {v}_{r}\leqslant 100$ $\mathrm{km}\,{{\rm{s}}}^{-1}$ and (2) a component of radially outflowing gas in the velocity range $120\leqslant {v}_{r}\leqslant 230$ $\mathrm{km}\,{{\rm{s}}}^{-1}$. Per our ${\rm{O}}\,{\rm\small{VI}}$ surface mass density analysis, we infer that, by mass, ∼70% of the absorbing gas is associated with the outflow, and this is consistent with the kinematic analysis. Given that there is roughly equal mass of inflow and outflow associated with the mass distribution of the low-velocity mode (${v}_{r}\sim +10$ $\mathrm{km}\,{{\rm{s}}}^{-1}$; see Figure 6(a)), we can now further infer that roughly ∼40% of the total outflow mass is in the form of $v\sim 200$ $\mathrm{km}\,{{\rm{s}}}^{-1}$ absorbing gas on a highly radial trajectory.

For $0.375\lt D/{R}_{\mathrm{vir}}\leqslant 0.75$, Figures 6(b), (e), and (h), we find that the total mass of ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas has decreased and that the radial outflow bimodality is significantly reduced relative to the bulk of absorbing gas in the velocity range $-100\leqslant {v}_{r}\leqslant 100$ $\mathrm{km}\,{{\rm{s}}}^{-1}$.

It is apparent that the majority of the high-velocity radial outflow detected in ${\rm{O}}\,{\rm\small{VI}}$ absorption within the projected region $D/{R}_{\mathrm{vir}}\leqslant 0.375$ has dynamically stalled prior to reaching the region $0.375\lt D/{R}_{\mathrm{vir}}\leqslant 0.75$. If a majority of the gas had reached this intermediate region and become dynamically mixed with the inflowing and outflowing gas in the velocity range $-100\leqslant {v}_{r}\leqslant 100$ $\mathrm{km}\,{{\rm{s}}}^{-1}$, we would expect the total gas mass to be higher than what is measured. Consistent with our inference from the ${\rm{O}}\,{\rm\small{VI}}$ mass surface density, we see that roughly 45% of the absorbing gas is inflowing and 55% is outflowing. The kinematics indicate that the 10% difference is in the form of $v\sim 200$ $\mathrm{km}\,{{\rm{s}}}^{-1}$ absorbing gas outflowing on a highly radial trajectory.

In the outer projected region $0.75\lt D/{R}_{\mathrm{vir}}\leqslant 1.125$ (Figures 6(c), (f), and (i)), we see a continuation of a drop in total ${\rm{O}}\,{\rm\small{VI}}$ mass for both inflowing and outflowing gas. We also see a remnant of the radial outflow apparent in the extended tail of the radial velocity distribution and the small peak in the direction cosine at $\alpha \sim +1$.

Corroborating the behavior of the ${\rm{O}}\,{\rm\small{VI}}$ mass surface density with increasing $D/{R}_{\mathrm{vir}}$, the kinematics indicate that the majority of the absorbing gas mass in the inner projected region of the halo of the simulated galaxies is outflowing gas. Further, the kinematics suggest that the drop in ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas mass with $D/{R}_{\mathrm{vir}}$ is because the majority of the outflowing gas is not making it out past $D/{R}_{\mathrm{vir}}\simeq 0.375$. There is also a significant decrease in the low-mass, infalling gas at $D/{R}_{\mathrm{vir}}\gt 0.75$ here than was present in the lower $D/{R}_{\mathrm{vir}}$ bins. This suggests a lack of accreting gas in the outer regions of the galaxy, which is potentially a result of the outflowing gas not making it past $D/{R}_{\mathrm{vir}}\simeq 0.375$ and preventing recycled accretion from occurring.

In a study of ${\rm{O}}\,{\rm\small{VI}}$ absorption-line kinematics using a slightly larger sample of simulated galaxies from the same VELA suite as our galaxies, Kacprzak et al. (2019) concluded that the ${\rm{O}}\,{\rm\small{VI}}$-absorbing outflows do not appear to be present beyond ∼50 kpc (approximately $D/{R}_{\mathrm{vir}}\sim 0.4$) and that the outflows eventually decelerate and fall back toward the galaxy. Thus, the characteristics we find for our smaller sample of galaxies (based on the selection criteria adopted for this study) are likely the same characteristics of the VELA simulations in general.

Strawn et al. (2020) also study selected simulated galaxies from the same VELA suite we used. Consistent with our findings, they report ${\rm{O}}\,{\rm\small{VI}}$-bearing gas is predominantly inflowing at $D/{R}_{\mathrm{vir}}\geqslant 0.3$. They further show that this inflowing gas is photoionized, whereas the bulk of the gas halo by volume (not mass) inside $D/{R}_{\mathrm{vir}}\leqslant 0.3$ is collisionally ionized. As we have shown, by mass, the majority of the gas selected by detectable ${\rm{O}}\,{\rm\small{VI}}$ absorption in mock spectra is photoionized (see the phase diagram, Figure 4). For comparison, we further examine the radial distribution of the ionization conditions of our absorption-selected gas by examining the phase diagram for $D/{R}_{\mathrm{vir}}\lt 0.3$ and $D/{R}_{\mathrm{vir}}\geqslant 0.3$. We find results consistent with Strawn et al. (2020). The vast majority (∼90%) of the inflowing ${\rm{O}}\,{\rm\small{VI}}$-absorption-selected gas mass at $D/{R}_{\mathrm{vir}}\geqslant 0.3$ is photoionized. However, for $D/{R}_{\mathrm{vir}}\lt 0.3$, the fraction of ${\rm{O}}\,{\rm\small{VI}}$-absorption-selected gas mass that is photoionized has decreased to ∼60%. Again, the characteristics we find for our smaller sample of galaxies likely reflect the same general characteristics of the VELA simulations.

Considering our findings and those of both Kacprzak et al. (2019) and Strawn et al. (2020), we suggest that it is the inability of the majority of winds to reach the outer halo that explains the dramatic drop in the ${\rm{O}}\,{\rm\small{VI}}$ covering fraction at $D/{R}_{\mathrm{vir}}\gt 0.375$ in the mock sample (see Figure 2). If the radial outflowing absorbing gas were to reach out to $D/{R}_{\mathrm{vir}}\simeq 0.75$, we might expect a flatter profile of the mass surface density across the region $0\leqslant D/{R}_{\mathrm{vir}}\leqslant 0.75$ and therefore a flatter ${\rm{O}}\,{\rm\small{VI}}$ absorption covering fraction out to $D/{R}_{\mathrm{vir}}\simeq 0.75$ (as measured for the observed sample).

If the simulations provide a window on the processes occurring in real galaxies, then we would be in a position to suggest that the galaxies comprising the observed sample have outflows that reach somewhat beyond $D/{R}_{\mathrm{vir}}=0.5$, and perhaps as far out as $D/{R}_{\mathrm{vir}}=0.75$. In a further analysis of the observed galaxy sample of KMC15 and the results of Nielsen et al. (2017), Ng et al. (2019) yielded evidence for the presence of outflows. Interestingly, these outflows might be required by the observations to have a clear spatial structure relative to the projected major axis as measured in the modes of the azimuthal distribution of the covering fraction (see Figure 3).

However, Nielsen et al. (2017) showed that ${\rm{O}}\,{\rm\small{VI}}$ absorption velocity spreads are similar regardless of galaxy inclination or azimuthal angle and suggest that ${\rm{O}}\,{\rm\small{VI}}$ kinematics is apparently not causally connected to the stellar activity of the central galaxy. Nielsen et al. (2017) speculate that the observed bimodal distribution in the ${\rm{O}}\,{\rm\small{VI}}$ covering fraction may be due to the gas-ionization conditions at intermediate azimuthal angles ($20^\circ \leqslant {\rm{\Phi }}\leqslant 60^\circ$) not favoring the ${\rm{O}}\,{\rm\small{VI}}$ phase. If this is true, the VELA simulations have not captured this characteristic of the CGM. Additionally, for a similar sample of observed galaxies, Ng et al. (2019) look at the ${\rm{O}}\,{\rm\small{VI}}$ gas kinematics relative to the galaxy when accounting for the mass of the host galaxy. Their results suggest that there is a strong halo mass dependence on ${\rm{O}}\,{\rm\small{VI}}$ absorber kinematics, and to coax out the signatures of outflows, one must account for the halo mass dependence, which we indirectly do in this work by normalizing both the observed and simulated samples by the virial radius.

As shown in Figure 3 and described in Section 3, the covering fraction of the mock sample was in good agreement with that of the observed sample in the azimuthal range $20^\circ \leqslant {\rm{\Phi }}\leqslant 60^\circ$, but the mock sample did not yield the azimuthal bimodality found for the observed sample. The fact that a significant azimuthal bimodality manifests in the observed sample strongly suggests that, from galaxy to galaxy, the spatial distribution of ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas is similarly organized with respect to the galaxy major axis. It also suggests this spatial structure is not averaged out when galaxies of various inclinations are included in the sample.

In contrast, our lack of an azimuthal bimodality in the mock sample would suggest that the ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas surrounding simulated galaxies either (1) does not exhibit a common spatial distribution from galaxy to galaxy, or (2) does not exhibit any structure to its spatial distribution, i.e., is quasi-isotropically distributed. For the former, each simulated galaxy could exhibit a well-defined spatial distribution, but this distribution would not be oriented similarly from galaxy to galaxy with respect to the galactic disk. In this case, "observing" the absorption through simulated galaxies of various inclinations might average out the covering fraction as a function of azimuthal angle. In the latter case, the amorphous spatial distribution of the absorbing gas would naturally not produce any modes in the covering fraction as a function of azimuthal angle.

In order to improve our understanding of the ${\rm{O}}\,{\rm\small{VI}}$ covering fraction distribution as a function of position angle obtained for the mock sample, we examined the full $360^\circ$ azimuthal distribution of the ${\rm{O}}\,{\rm\small{VI}}$ absorption covering fraction for each of the simulated galaxies as a function of inclination angle. The sky projection of the galaxies was defined so that the axis of rotation for the inclination is in the plane of the sky (perpendicular to the observer).

In Figure 7, we present the full sky galaxy-to-galaxy azimuthal variation in the covering fraction in azimuthal bins of $15^\circ$, where a line connecting the nodes $0^\circ$ and $180^\circ$ defines the projected major axes of the galaxies on the plane of the sky (corresponding to the projection of the galaxy disks), and where the nodes $90^\circ$ and $270^\circ$ define the sky-projected minor axes of the galaxies (corresponding to the projected polar axes of the galaxies.) For each simulated galaxy, we present the covering fraction at three sky-projected inclination angles, $i=0^\circ$, $60^\circ$, and $90^\circ$. For each galaxy projection, all 1000 LOSs with random distributions of $D/{R}_{\mathrm{vir}}\leqslant 1.5$ and azimuthal angle, are included.

Zoom In Zoom Out Reset image size

We see that, around a given galaxy, the ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas is not uniformly distributed and that it has a clear spatial structure as inferred from absorption-line analysis. For each galaxy, the covering fraction exhibits substantial variation with azimuthal angle and galaxy inclination. We also find that, for a fixed inclination (or sky viewing angle), the azimuthal distribution of covering fraction varies substantially from galaxy to galaxy.

It would seem that one outcome of constructing the mock samples is an azimuthal averaging of the covering fraction. Indeed, as presented in Figures 8(a)–(c), we find that when multiple simulated galaxies are observed at a given inclination angle (i.e., many LOSs through several galaxies of identical inclination), the azimuthal distribution of the ${\rm{O}}\,{\rm\small{VI}}$ covering fraction becomes significantly more symmetric for each of the three illustrated inclinations.

Zoom In Zoom Out Reset image size

Observationally, the full $360^\circ$ azimuthal distribution cannot be studied due to the degeneracy of the four polar quadrants; only the primary angle ($0^\circ \leqslant {\rm{\Phi }}\leqslant 90^\circ$) can be studied. For the $360^\circ$ distributions shown in Figures 8(a)–(c), we present the covering fraction as a function of the primary azimuthal angle in Figure 8(d)–(f). These are not comparable to the full mock samples shown in Figure 3, as these distributions are for 1000 LOSs and fixed galaxy inclinations, whereas the mock sample comprises 15,000 realizations of 53 LOSs drawn from the observed distribution of impact parameters and inclinations. The simulated galaxies yield a relatively flat azimuthal distribution of ${\rm{O}}\,{\rm\small{VI}}$ covering fraction for each of the three presented inclination angles ($i=0^\circ$, $60^\circ$, and $90^\circ$), although at both $i=60^\circ$ and $90^\circ$, there is a slightly higher covering fraction along the galaxy major axis (${\rm{\Phi }}=0^\circ ,180^\circ$). For comparison, we also present the observed ${\rm{O}}\,{\rm\small{VI}}$ covering fractions in Figures 8(d)–(f) and note that the ${\rm{O}}\,{\rm\small{VI}}$ covering fraction of the simulated galaxies matches observations well in the azimuthal region $20^\circ \leqslant {\rm{\Phi }}\leqslant 60^\circ$. As shown in Figure 2, there is still a general trend of decreasing ${\rm{O}}\,{\rm\small{VI}}$ covering fraction with increasing azimuthal angle.

In summary, as inferred from Figure 7, we find that individual simulated galaxies have (1) clearly defined and detectable sky-projected spatial structure in ${\rm{O}}\,{\rm\small{VI}}$ absorption and (2) unique spatial structure from galaxy to galaxy with no common orientation relative to the galactic disk. As such, the process of adopting "many LOSs through several galaxies" averages out the sky projection of this spatial structure, even when the inclination of the galaxies is a controlled variable (see Figure 8). Given the complicated nature of the CGM, we might expect that observed galaxies would exhibit a similar, somewhat randomized structure that would wash out any trends in covering fraction with azimuthal angle. However, the observed bimodality tantalizingly suggests that "real-world" galaxies may have similar spatial structure from galaxy to galaxy.

From the work of Kacprzak et al. (2019), we know that the ${\rm{O}}\,{\rm\small{VI}}$-absorbing gas surrounding the simulated galaxies in the VELA simulations resides in two main spatial/kinematic features. The first are outflowing and recycling gas structures within CGM, which reside roughly within 1R_vir of the galaxies. The second are thick inflowing filaments that exist at farther galactocentric distances but penetrate into the CGM. While these structures have been identified in a larger sample of simulated VELA galaxies than studied here, our analysis indicates that such structures are at different spatial locations for different simulated galaxies. Indeed, our analysis suggests that such structures can give rise to ${\rm{O}}\,{\rm\small{VI}}$ absorption and can be "observed" by examining the sky-projected covering fraction of the ${\rm{O}}\,{\rm\small{VI}}$ absorption (as seen in Figure 7).

Using EAGLE simulations and TNG50 of the IllustrisTNG simulations, Péroux et al. (2020) recently showed outflowing gas has preferentially higher metallicity than inflowing gas. They further show higher metallicity gas along the minor axis, suggestive of enrichment of gas from galactic outflows, and lower metallicity accreting gas along the major axis. This bimodality becomes stronger with increasing impact parameter. In the VELA simulations, we do not see an increase of ${\rm{O}}\,{\rm\small{VI}}$ along the minor axis (as is seen in the bimodality in our observed sample); however, we also show (see Figure 6) that ${\rm{O}}\,{\rm\small{VI}}$ outflows do not extend significantly past $D/{R}_{\mathrm{vir}}\gt 0.375$. As the majority of the outflowing gas is not making it out to $D/{R}_{\mathrm{vir}}\approx 0.375$ in the VELA simulations, this may potentially explain why the VELA simulations do not reproduce the bimodal covering fraction distribution with azimuthal angle.

One contributing factor as to why the bimodal spatial covering fraction is not seen in the VELA simulations could be that the simulated galaxies are studied at $z\simeq 1$, whereas the mean redshift of the observed galaxies is $z\simeq 0.3$. As redshift decreases, the angular area subtended by infalling filamentary streams is expected to increase (e.g., Dekel et al. 2009b; Danovich et al. 2012), and this could yield a corresponding evolutionary change in the azimuthal distribution of the ${\rm{O}}\,{\rm\small{VI}}$ covering fraction from $z\simeq 1$ to $z\simeq 0.3$.

A gentler, more spread out accretion could yield a greater influence of the galaxy stellar feedback in governing the spatial distribution of the ${\rm{O}}\,{\rm\small{VI}}$-bearing gas. Potentially, this could cause the spatially stochastic structure at $z\simeq 1$ (see Figure 7) to transform into a systematic structure that gives rise to the observed azimuthal bimodality by $z\simeq 0.3$ (as was seen by Péroux et al. 2020). While we attempted to control for the influence of mergers by selecting simulated galaxies that had no major mergers for up to 1 Gyr prior to the redshift of study, it is currently not possible to directly test the VELA simulations below z = 1 as to whether redshift evolution might resolve the fact that we do not find a bimodal azimuthal covering fraction.

All things considered, we can only speculate how the CGM of the VELA simulations we have studied would evolve from $z\simeq 1$ to $\simeq 0.3$. Observational data and theory are somewhat contradictory in guiding our insights. Though higher ${\rm{O}}\,{\rm\small{VI}}$ column density appears to be associated with star-forming galaxies with specific SFR greater than 10⁻¹¹ yr⁻¹ (Tumlinson et al. 2011), it is hard to reconcile this finding with growing observation evidence (Bielby et al. 2019; Zahedy et al. 2019; Ng et al. 2019) and the theoretical expectation (Oppenheimer et al. 2016; Nelson et al. 2018; Oppenheimer et al. 2018; Qu & Bregman 2018) that the ${\rm{O}}\,{\rm\small{VI}}$ column density peaks for galaxies with halo mass $\mathrm{log}{M}_{h}/{M}_{\odot }\simeq 12$.

It is difficult to interpret the relationship between ${\rm{O}}\,{\rm\small{VI}}$-absorbing CGM gas and stellar-driven outflows and accretion in view of the observed gas kinematics and metallicities. Nielsen et al. (2017) find that the velocity dispersion of the gas has no systematic differences with galaxy orientation or color, Kacprzak et al. (2019) find that the global kinematics of the gas shows no coupling to the host galaxy rotation, and Pointon et al. (2019) find that there is no clear metallicity trend between the minor-axis and major-axis gas.

Despite the difference in redshift of the VELA simulations and the observed galaxy sample, because galaxies mostly undergo passive evolution from redshift $z\simeq 1$ to $z\simeq 0.3$ (e.g., Satoh et al. 2019), and because the cooling time of the hot, diffuse gas is longer than a Hubble time (Sutherland & Dopita 1993; Wiersma et al. 2009), we do not expect that our findings would dramatically change. While there is a dependency of ${\rm{O}}\,{\rm\small{VI}}$ absorption on halo mass (Ng et al. 2019), the mean virial masses of the observed and simulated galaxies are quite similar ($\mathrm{log}({M}_{\mathrm{vir}}/{M}_{\odot })=11.8,11.6$ respectively), so we expect this to not significantly affect our results.

One of the greatest recent challenges is to duplicate the azimuthal spatial dependence, as discussed throughout this work. Interestingly, it is only very recently that a subset of simulations has achieved outflows concentrated along the galaxy minor axis (Nelson et al. 2019; Mitchell et al. 2020; Péroux et al. 2020). However, it remains to be seen exactly how these spatial-kinematic CGM structures translate to covering fractions obtained through mock absorption-line experiments (using many LOS through many simulated galaxies).

It is only recently that pioneering studies by Hummels et al. (2019), Peeples et al. (2019), and van de Voort et al. (2019) have investigated resolution effects on the nature of the simulated CGM gas. Peeples et al. (2019) find that absorption lines arise in gas structures with masses below 10⁴ M ${}_{\odot }$, while Hummels et al. (2019) predict CGM gas structures with sub-M ${}_{\odot }$, both substantially lower than the mass resolution of zoom simulations, such as the VELA simulations.

Improvements of the VELA simulations used for this work, which modeled only thermal-energy supernova feedback and radiative feedback from stars, include adding a nonthermal radiation pressure to the total gas pressure in regions where ionizing photons from massive stars are produced. Modern VELA simulations have been run with enhanced momentum as part of the supernova feedback (see Gentry et al. 2017), and this new feedback mechanism yields improved stellar mass–halo mass relations (D. Ceverino et al. 2020, in preparation). It remains to be seen whether the new feedback mechanisms generate a more realistic CGM relative to the VELA simulations used in this work. We will compare the results of this paper to this new generation of the VELA simulations in future work.