Cloud-by-cloud, multiphase, Bayesian modeling: Application to four weak, low ionization absorbers

The circumgalactic medium (CGM) is the tenuous, multiphase medium in the environment of a galaxy. It plays an important role in steering the evolution of galaxies. The massive reservoir of gas in the CGM seeds the formation of stars which evolve and distribute their energy and content back into the galaxy. CGM gas of galaxies can be probed in intricate detail by using metal absorption lines (e.g., Bergeron 1986; Bergeron & Boissé 1991; Lanzetta & Bowen 1992; Steidel & Sargent 1992; Churchill et al. 1996, 2000; Churchill et al. 2020; Adelberger et al. 2005; Kacprzak et al. 2008, 2015, 2019; Chen & Mulchaey 2009; Steidel et al. 2010; Nielsen et al. 2013a, b; Stocke et al. 2013; Werk et al. 2013; Bordoloi et al. 2014; Turner et al. 2014; Johnson et al. 2015; Richter et al. 2016; Keeney et al. 2017; Muzahid et al. 2018) observed in the spectra of background quasars. Absorption-line spectroscopic studies provide information about the metallicities, densities, temperatures, and kinematic structure of various parcels of gas along the lines of sight. These studies are important because they potentially allow us to learn about the processes that supply inflowing gas to the galaxy, enrich the surroundings with metal-rich outflows, minor and major mergers, and send enriched gas back to the galaxy as recycled accretion.

In this paper, we present a new Bayesian method for extracting the physical properties of the multiphase gaseous medium probed by quasar absorption-line systems. To demonstrate the method, we apply it to four, well-studied, weak, low ionization, systems, also known as weak Mg ii absorbers (Narayanan et al. 2005). These systems were chosen because they are relatively simple in the number of low ionization components, but still have multi-phase structure, and also because the origins of these tiny pockets of high metallicity, far from galaxy centers, are of particular interest.

Weak Mg ii-absorbers are defined by rest frame equivalent width $W_{\_}r^{2796}<$0.3 Å. Unlike strong Mg ii absorbers, which are typically found within $\sim 60$ kpc of a luminous galaxy (Nielsen et al. 2013b, 2018), the weak Mg ii absorbers are found to be associated with lower neutral hydrogen column densities, and are at larger impact parameters within the circumgalactic environment, typically $>40$ kpc (Nielsen et al. 2013b, 2018). Churchill et al. (1999) conducted the first systematic survey of weak Mg ii absorbers in the redshift range $0.4<z<1.4$ using HIRES/Keck spectra, and found that the number of weak absorbers per unit redshift of $d\mathcal{N}/dz\approx 1.74\pm 0.10$ is $\approx$4 times higher than that of the Lyman limit systems (LLSs) which have neutral hydrogen column density $N(\mbox{H\,{\sc i}})>10^{17.2}$ cm^-2, indicating an origin of the weak absorbers in sub-LLSs (Churchill et al. 2000). A tantalizing aspect of the weak Mg ii-absorbers is that, generally, their low ionization phases exhibit near-solar to super-solar metallicities (Rigby et al. 2002; Lynch & Charlton 2007; Misawa et al. 2008; Narayanan et al. 2008) in spite of the fact that luminous galaxies are typically 50$-$200 kpc away. Detailed photoionization models of the intermediate redshift, weak Mg ii-absorbers show that they often have two phases, a high-density region/s that is 1-300 pc thick and produces narrow (a few km s^-1) low-ionization lines, and kiloparsec-scale, lower density region/s that produces broader, high-ionization lines (Charlton et al. 2003; Narayanan et al. 2007). By comparing the relative number of Mg ii and C iv absorbers, Milutinović et al. (2006) concluded that the low ionization absorption must occur in flattened or filamentary structures that occupy more than half the solid angle subtended by structures that produce only C iv absorption.

Studies of weak absorbers at low-$z$ offer the advantage of ease of detection of host galaxies. At low redshift, $0<z<0.3$, Narayanan et al. (2005) measured $d\mathcal{N}/dz=1.00\pm 0.20$ for weak Mg ii absorbers (with $W_{\_}r^{2796}>0.02$ Å) using Si ii and C ii as proxies in $HST$/COS spectra. This is smaller than expected based on the evolution expected for the intermediate redshift population, given the decreasing extragalactic background radiation. This is because the low redshift population should consist both of 1-300 pc absorbers that gave rise to Mg ii absorption at intermediate redshift, and larger clouds that gave rise only to higher ionization (C iv) absorption at intermediate redshift. Understanding which of these types of absorbers dominate the population at low redshift is important to understanding the processes that produce these mysterious, high metallicity objects abundant in the CGM.

More generally, both observational and theoretical studies (e.g., Stinson et al. 2012; Ford et al. 2013; Shen et al. 2013; Suresh et al. 2016) of absorbers reveal the existence of cool (10⁴ K) and warm-hot (10⁵-10⁶ K) phases in the CGM, often all detected along the same sightline to the background quasar. Different ionization states of material trace different phases of gas, which have different densities and temperatures. Cooler, higher-density gas is traced by lower-ionization species (e.g., Mg ii, Si ii) and shows structure on smaller scales compared to the warmer and/or lower-density phase (e.g., C iv or O vi). Length scales of the lower ionization clouds range from $\sim$ 1 pc (McCourt et al. 2012; Liang & Remming 2020) to a few kpc. As mentioned above, strong Mg ii absorption is detected along most sightlines within an impact parameter of 50 kpc from a galaxy (Nielsen et al. 2013a), while weaker Mg ii absorption is found at a median impact parameter of 166 kpc (Muzahid et al. 2018). C iv is often seen out to 100 kpc (Chen et al. 2001; Bordoloi et al. 2014), while O vi is seen out to at least 150 kpc (Prochaska et al. 2011; Tumlinson et al. 2011; Stocke et al. 2013; Liang & Chen 2014; Pointon et al. 2019) of star forming galaxies, but not around quiescent galaxies. Gas traced by H i absorption is seen well beyond (Liang & Chen 2014; Borthakur et al. 2015).

Interpretation of the rich absorption line data is aided by identification of galaxies near the quasar sightlines. It is expected that the orientations of the galaxies will correlate with absorption properties (Bordoloi et al. 2011; Bouché et al. 2012; Kacprzak et al. 2012; Kacprzak et al. 2015; Nielsen et al. 2015). In the simplest view, galactic winds due to episodes of star formation are expected to produce high metallicities in the outflowing gas, and should be seen primarily along the minor axis (Muzahid et al. 2015; Nielsen et al. 2015; Pointon et al. 2019; Schroetter et al. 2019; Peroux et al. 2020). Inflowing pristine gas from the intergalactic medium would be expected to flow in along the major axis. Another important observational result is a bimodality in the metallicity of partial Lyman limit systems, with the high metallicity peak at $\log Z$ = $-0.4$ and the low metallicity peak at $\log Z$ = $-1.7$ (Lehner et al. 2013). It has been hypothesized that this bimodality could be due to the contrast between outflowing and infalling material.

A coherent and logical picture is beginning to emerge from the data (Tumlinson et al. 2017), however, there are some significant problems plaguing interpretations, and the root cause of these problems could be the methods used to derive the properties of the gas from the absorption profiles, particularly the metallicities. Often the total column density (integrating over all components at all velocities) of a given metal line transition is compared to the total column density of hydrogen, in the context of photoionization models, to derive the metallicity. This overlooks the importance of phase structure, with more than one phase contributing to absorption at a single velocity i.e. some transitions could have contributions from two different phases which can lead to an inaccurate estimate of ionization parameter, and thus a systematic error in deriving a single metallicity value. In other words, the average parameters inferred for a system may “wash out” contributions from regions of gas that have metallicities and densities that differ significantly from the means, and the possibility of conditions that significantly change with velocity because of changes in physical conditions along the line of sight. Simulations (Churchill et al. 2015; Peeples et al. 2019) indicate absorption arises in lots of places along the sightline, motivating why multiple metallicities are expected.

In a recent study, Zahedy et al. (2019) have carried out a systematic analysis of 16 intermediate-redshift ($z$=0.21-0.55) luminous red galaxies (LRGs), finding metal-poor absorbing components with $<$ 1/10 solar metallicity in half of the LRG haloes in tandem with solar and super-solar metallicity gas in the same halo, suggestive of poor chemical mixing and complex multiphase structure in these haloes. This highlights the importance of resolving multiple components of an absorption system using high resolution absorption spectra, which otherwise is missed by ionization studies that make use of only the integrated H i and metal column densities along individual sightlines.

This work improves the efficiency of component by component modeling that has been successful in recovering the physical conditions for various individual absorbers (e.g. Churchill & Charlton 1999; Charlton et al. 2000; Ding et al. 2003a; Charlton et al. 2003; Ding et al. 2003b; Zonak et al. 2004; Ding et al. 2005; Masiero et al. 2005; Lynch & Charlton 2007; Misawa et al. 2008; Lacki & Charlton 2010; Jones et al. 2010; Muzahid et al. 2015; Richter et al. 2018; Rosenwasser et al. 2018). Rather than averaging over components and phases, it is possible to determine how much of the H i is associated with these different phases in order to derive separate metallicities for various clouds. Resolving the individual clouds allows us to break the degeneracy for components on the flat part of the Ly$\alpha$ curve of growth, even with coverage of just saturated H i lines, and derive metallicity constraints for different parcels of gas along the line of sight. It is important to do so because different processes e.g., outflows (Bouché et al. 2012; Bordoloi et al. 2014; Rubin et al. 2014; Schroetter et al. 2016), pristine accretion (Rubin et al. 2012; Martin et al. 2012; Danovich et al. 2015), recycled accretion (Ford et al. 2014), minor and major mergers (Martin et al. 2012; Anglés-Alcázar et al. 2017), are surely contributing to the same system, and it is expected that conditions will vary significantly along a line of sight which can span hundreds of kpc spatially (Churchill et al. 2015; Peeples et al. 2019). This will lead to a more meaningful comparison to galaxy properties. For example, Pointon et al. (2019) did not find a difference between the metallicities of absorbers found within an impact parameter of 200 kpc along the major and the minor axes of isolated galaxies. Based on cosmological hydrodynamic simulations, a larger metallicity is expected along the minor axis due to outflows and a lower metallicity along the major axis due to inflows (Peroux et al. 2020). However, an observational trend could exist, for example, for the minor axis to have some, but not all, high metallicity components, or for the minor axis to have one or more low metallicity components. Such results would be “washed out” by deriving an average metallicity for all gas along a line of sight, which clearly often has multiple complex origins. For some datasets/projects the new analysis could be transformative, however to make it feasible to use for large statistical studies it is important that the analysis is semi-automated and robust.

In this paper, we describe a new method to derive the physical properties of quasar absorption line systems using the photoionization code cloudy (ver 17.01; Ferland et al. 2017), Bayesian methods, and parallel processing to efficiently obtain results. Though the code requires some human intervention to choose transitions that define the phase structure, it is mostly automated, and provides uncertainties on parameters in the context of the assumed structure. The paper is organized as follows: in Section 2 we describe the spectral observations that are analysed in this work; in Section 3 we describe the methodology used to determine the physical conditions of an absorber; in Section 4 we present the results and discuss our findings, comparing to previous work, applying tests, and discussing possible limitations. We conclude in Section 5 by discussing the limitations and potential of our cloud-by-cloud, multiphase, Bayesian modeling (CMBM) method. Throughout this work, we assume a cosmology with $H_{\_}{0}\,=\,70$ km s^-1 Mpc^-1, $\Omega_{\_}{\rm M}\,=\,0.3$, and $\Omega_{\_}{\Lambda}\,=\,0.7.$ Abundances of heavy elements are given in the notation $\rm[X/Y]=\log(X/Y)-\log(X/Y)_{\_}{\odot}$ with solar relative abundances taken from Grevesse et al. (2010). All the distances given are in physical units. All the logarithmic values are presented in base-10.

We apply our methods to four low redshift absorption line systems from the COS-Weak sample (Muzahid et al. 2018). The background quasars in this study have UV spectra from the Cosmic Origins Spectrograph (COS) instrument on the Hubble Space Telescope ($HST$) and optical spectra from High Resolution Echelle Spectrometer (HIRES) installed on the Keck telescope or Ultraviolet and Visual Echelle Spectrograph (UVES) installed on the Very Large Telescope (VLT). The four absorption systems that we consider have spectroscopic redshifts $z<0.30$. These weak, low-ionization systems are relatively simple ones compared to some stronger low-ionization absorbers, but they are expected to have multiple gas phases and provide good test cases for our methodology. Table 1 presents the details of the observations of the four systems studied in this work.

We model the observed absorption system using the photoionization code cloudy, and synthesize the expected absorption profiles and compare them to the observed profiles in order to infer the physical conditions of the absorbing gas. What makes this method different from other researchers’ methods is that we model the entire absorption system across all transitions by splitting it into its constituent clouds.

We begin by choosing a “trustworthy” low ionization species (the “optimized ion”) for a Voigt profile (VP) fit, meaning one with lines that are unsaturated, and observed at a high resolution, and with a large $S/N$. The starting point is a fit of Voigt profiles to all transitions within this single ionization species (each parcel of gas described by a single VP will be called a “cloud”). For this optimized transition we determine the parameters of column density ($N$), Doppler parameter ($b$), and redshift ($z$) for each cloud using a Monte Carlo technique employed using PyMultiNest (Buchner et al. 2014). PyMultiNest is a python implementation of MultiNest (Feroz et al. 2009), which is a robust nested sampling method that draws a new, uniformly random point with higher likelihood through an ellipsoidal rejection sampling scheme (Shaw et al. 2007). It is efficient at sampling from a multimodal posterior distribution. We model Voigt profiles using an analytic approximation (Tepper-García 2006) for the Voigt-Hjerting function as implemented in the VoigtFit package (Krogager 2018). We obtain VP fits to the absorption profiles for the optimized ion, allowing for more than one component. We fit the profiles of the optimized ion with a VP model corresponding to the highest log-evidence and the least model complexity. After we determine the observed values of $\log N$, $b$, and $z$ for each cloud, we use the parameters of that VP fit for the optimized ion as a starting point for models in cloudy.

cloudy has many parameters that can be adjusted. The two most important parameters are metallicity relative to solar abundance ($Z$) and ionization parameter ($U$), defined by $U=n_{\_}{\gamma}/n_{\_}{H}$, where $n_{\_}{\gamma}$ is the volume density of ionizing photons. Ionization parameter serves as a direct proxy for the hydrogen number density, but may be more intuitive when thinking of contributions of different clouds to the observed ionization states. We adopt the Haardt & Madau (2012)[HM12] model for the extragalactic background radiation (EBR) to determine $n_{\_}{\gamma}$ at the given redshift. The systematic effects arising from the choice of the EBR are discussed in Keeney et al. (2017). We adopt the HM12 model as it best reproduces the low-redshift H i column density distribution (Danforth et al. 2016) for $\log[N(\mbox{H\,{\sc i}})/\hbox{cm${}^{-2}$}]>$ 14 (Shull 2014). We assume a solar abundance pattern (Grevesse et al. 2010) for our modeling. Detailed constraints on abundance pattern are not possible without higher $S/N$ coverage of a larger number of transitions.

It is a time-consuming endeavor to construct cloudy models component-by-component and phase-by-phase for individual absorbers because of the limitations imposed by serial computing. However, it is clear that doing so will lead to a significant gain in our ability to correctly infer the presence of metallicity structure in the CGM. To increase efficiencies we have developed robust python based codes that employ distributed parallel processing on the CyberLAMP infrastructure ¹¹¹https://wikispaces.psu.edu/display/CyberLAMP to obtain models in a shorter time.

The cloudy photoionization equilibrium (PIE) models are generated for the identified optimized ion tracing the low ionization gas phase for a broad range of values over the $\log Z$ $\times$ $\log U$ grid in the range: $\log Z$ $\in$ [-3.0,1.5] and $\log U$ $\in$ [-6.0, 0.0], obtained with a step-size of 0.1. The cloudy models are then linearly interpolated onto the $\log Z$ $\times$ $\log U$ grid with a spacing of 0.01 dex. The interpolating function is generated by triangulating the input data with Qhull (Barber et al. 1996), and on each triangle performing linear barycentric interpolation. The interpolation error is of the order $\sim O(h^{2})$, $h$ being 0.1, and introduces only small errors $\sim$ 0.01 dex in $\log Z$ and $\log U$. The total column density of hydrogen, $N(H)$, is adjusted for each model until the column density of the optimized ion is matched for each point on the $\log Z$ $\times$ $\log U$ grid. We allow for an uncertainty of 0.05 dex in the column density matching when generating the cloudy models. The uncertainties arising due to the interpolation error and the error in column density are not accounted for in our error estimation, but are usually negligible in comparison to other sources of error. We quote 2$\sigma$ uncertainties in our results for the best model parameter values, which account only for the photon noise in the comparison between the synthesized model profiles and the observed spectra.

We then determine the parameters that best describe the conditions ($\log Z$, $\log U$) present in each one of the absorbers again using PyMultiNest. The posterior distribution of the parameters, which is proportional to the product of priors and likelihood, is determined by using uniform priors on $\log Z$ and $\log U$, spanning the grid range mentioned above. The log-likelihood function, $P$, is formulated as

Typically, all of the observed line transitions cannot be reproduced with one phase of gas. It is often the case that a model that fits the low-ionization transitions (e.g., Mg ii, Fe ii) cannot produce the observed amount of high-ionization absorption (e.g., C iv, O vi), or adequately match the profiles of the full Lyman series. In such cases, we will use more than one optimized ion. The other optimized ions represent the different phases of gas, with distinct regions having different densities along the line-of-sight. For example, we would typically optimize on Mg ii for the low-ionization gas, and on C iv for high-ionization gas. In some intervening absorbers, O vi, because of its comparatively high ionization potential ($\sim$ 114 eV), arises in gas with T $>$ 10⁵ K where the ionization is governed by collisions between energetic free electrons and metal ions (Pradeep et al. 2020). We model the gas phase traced by O vi to be under collisional ionization equilibrium (CIE). For such collisional ionization models, the cloudy grids are run over $\log Z$ $\times$ $\log T$ space: $\log Z$ $\in$ [-3.0,1.5] and $\log[T/K]$ $\in$ [5.0, 5.8]. Similar to the PIE modeling, the CIE models are obtained by varying the total column density of hydrogen, $N(H)$, until the column density of O vi is matched for each point on the $\log Z$ $\times$ $\log T$ grid. In Figure 1, we show for a metallicity of $\log Z=-1.0$ and different gas densities that cloudy can produce accurate results when modeling gas in CIE. The column density ratios of C iv/O vi agree well with the CIE model (Gnat & Sternberg 2007) for temperatures in the range of $10^{5}-10^{5.8}$ K, and are independent of the hydrogen number densities. We contrast the modeling results obtained with cloudy ver6.02d and ver17.01 to show that the apparent discrepancy with the Gnat & Sternberg (2007) model arises because of the evolving atomic database (Chatzikos et al. 2015). The cloudy models are generated for an assumed value of $\log n_{\_}{H}=-3.9$, based on the Milky Way halo density at 50kpc (Kaaret et al. 2020).

To summarize, we start with the assumption that a single, low ionization phase describes an absorption system. If such a simple, one-phase model fails to account for the observed column densities in all the transitions, we incorporate additional phases until the synthesized spectra from cloudy best match the observed spectra for all transitions. Generally, we find that a single phase is insufficient to explain the observed abundances of all the transitions even when the transitions coincide in velocity space. If the synthetic spectra fit the data for all transitions, then it is plausible that the physical conditions in the cloudy model match those in the gas that produces the absorption. This allows us to determine precise allowed ranges of metallicity, density, and temperature for each individual cloud that create the observed transition lines. A summary of the methodology is given in the form of a flowchart in Figure 2.