Mg ii Absorbers in High-resolution Quasar Spectra. I. Voigt Profile Models

We have analyzed 249 high-resolution, high signal-to-noise High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) and Ultraviolet and Visual Echelle Spectrograph (UVES; Dekker et al. 2000) quasar spectra obtained from the Keck and Very Large Telescope (VLT) observatories, respectively. The wavelength coverage of the spectra ranges from approximately 3000–10000 Å, though there is variable coverage in this range from spectrum to spectrum.

The resolving power of both instruments is R = λ/Δλ =45,000, or ∼6.6 km s⁻¹, and the spectra have three pixels per resolution element. The resolution is approximately constant in velocity as a function of observed wavelength. The signal-to-noise ratios are typically 25–80. Such high-quality spectra allow for an in-depth investigation into the distributions and kinematics of the galactic halos and intergalactic structures selected by the presence of metal line absorption.

The Keck/HIRES spectra were obtained from various observing programs prior to the creation of the Keck Observatory Archive, ⁴ and were donated to this work by Charles Steidel, J. Xavier Prochaska, Christopher Churchill, and by Michael Rauch and the late Wallace Sargent. The spectra of Churchill and Steidel were explicitly observed for Mg ii absorption (e.g., Churchill & Vogt 2001; Steidel et al. 2002; Churchill et al. 2003) and were selected based upon previous studies that used low-resolution spectra to detect Mg ii absorption of equivalent widths W_r (2796) ≥ 0.3 Å (Sargent et al. 1988; Steidel & Sargent 1992). Those of Sargent and Rauch were selected for high-resolution analysis of Lyα forest and C iv absorption, and those of Prochaska for damped Lyα (DLA) absorption (Prochaska et al. 2007).

The VLT/UVES spectra were acquired through the efforts of the UVES SQUAD prior to their Data Release 1 (Murphy et al. 2019). These archived UVES spectra were obtained by several researchers for a variety of scientific purposes, including Lyα forest, Lyman-limit and damped Lyα systems, and Mg ii, C iv, O vi, and other metal absorption line studies. The heterogeneous selection bias of our sample is addressed in Section 3.1.

An example UVES spectrum, the z_em = 2.406 quasar J222006−280323, is shown in Figure 1. Strong emission lines are labeled, as are the regions of Lyα forest absorption and A-band and B-band atmospheric absorption. The journal of quasar spectra used in this study is listed in Table 1. For each quasar, the columns list (1) the quasar name, taken from Veron-Cetty & Veron (2001), (2) a common B1950 name (mostly from the catalog of Hewitt & Burbidge 1996), (3) the emission redshift, (4) the lower wavelength limit observed, (5) the upper wavelength limit observed, and (6) the instrument with which the spectrum was obtained.

Zoom In Zoom Out Reset image size

The spectra of Churchill were reduced using the standard Image Reduction and Analysis Facility (IRAF ⁵ ) software, the process of which is detailed in Churchill (1997). Those of Sargent, Rauch, Prochaska, and Steidel were reduced using the Maunakea Echelle Extraction (MAKEE) data reduction package of Barlow (2005), which is optimized for the spectral extraction of single, unresolved point sources. All HIRES spectra were wavelength calibrated, using ThAr lamps, to the vacuum heliocentric standard at rest, and were continuum fit by their respective donors.

The UVES spectra were reduced using the UVES pipeline in the MIDAS environment (Dekker et al. 2000), which is provided by the European Southern Observatory (ESO). The wavelength solution is determined using a standard ThAr lamp exposure. Finally, the quasar flux is extracted and the ThAr wavelength calibration polynomial for each order is attached, having been corrected to vacuum heliocentric velocities. The individual exposures were then combined into one-dimensional spectra using the UVES POst-Pipeline Echelle Reduction (UVES Popler) software (Murphy 2008, 2016; Murphy et al. 2019), which is an extra reduction step that facilitates cosmic ray removal. An initial continuum is derived by iteratively fitting small sections of each spectrum with a low-order Chebyshev polynomial, rejecting points lying many sigma below or above the fit at each iteration. The sections of the continuum are then spliced together by linearly weighting each section from unity at their centers to zero at their edges. Any remaining artifacts, such as from internal spectrograph reflections, are removed and the continuum is refined manually.

For 25 quasars, multiple spectra were available. We optimally combined them in order to exploit all wavelengths covered and to maximize signal-to-noise ratio. Thus, in some cases, the final version of the spectrum we studied may comprise, for example, two unique HIRES observations and one UVES observation. This provided a single spectrum for a quasar. We ensured wavelength alignment by performing cross-correlations on unresolved features in regions of overlapping wavelengths. Possible nonlinear effects were accounted for by fitting a first-order polynomial to the cross-correlation shifts as a function of wavelength. The fluxes in the pixels were optimally averaged, weighted by the inverse of their variances using flux conservation and pixel interpolation.

We limited our search for absorption lines to wavelength regions spanning redward of the quasar Lyα emission line up to 5000 km s⁻¹ blueward of the quasar Mg ii emission line. The lower wavelength limit ensures Mg ii absorption is not confused with Lyα forest lines, and the upper limit is our adopted criterion for detected Mg ii absorption features to not be considered associated with the quasar vicinity (e.g., Weymann et al. 1991). Using our software Search (e.g., Churchill et al. 1999), the spectra are objectively scanned for Mg ii doublet candidates. The initial criteria for a candidate doublet are that the detection of a 5σ feature, which is taken to be the λ2796 line at redshift z = λ/2796.352 − 1, is accompanied by a corresponding feature with a 3σ detection at the projected location λ = 2803.531 × (1 + z). Doublet candidates, preliminary equivalent widths, and detection significance levels follow the formalism of Schneider et al. (1993). The candidates are then checked for doublet ratios of 1.0 ≤ DR ≤ 2.0 consistent within errors.

Once a candidate Mg ii doublet is identified, the Search software is capable of locating and examining numerous associated atomic transitions simultaneously. For this work, examination of the associated absorption features was limited to 13 commonly observed transitions from five abundant chemical elements. Including the Mg ii λλ2796, 2803 doublet, these transitions were Mg i λ2853, the Fe ii λ2344, λ2374, λ2383, λ2587, and λ2600 quintuplet, the Mn ii λ2577, λ2594, and λ2606 triplet, and the Ca ii λλ3935, 3970 doublet. The transitions and their adopted atomic data (Moore 1970; Cashman et al. 2017) are listed in Table 2. The columns are (1) the ion and transition, (2) the transition wavelength in vacuum, (3) the oscillator strength, and (4) the natural broadening, or damping constant.

Candidate Mg ii doublet line profiles were aligned in rest-frame velocity and visually inspected to ensure that they exhibit velocity alignment and flux decrements consistent with the radiative transfer and atomic physics of the two transitions, and to identify and annotate any blends or spurious spectroscopic features. For Mg ii λλ2796, 2803 doublet features passing the criteria for inclusion into the sample, we adopt the definition that two or more Mg ii λ2796 absorption features comprise a single absorption system if they reside within ±800 km s⁻¹ of each other. ⁶ The corresponding transitions listed in Table 2 are also visually checked.

Within the higher wavelength regimes of our sample (above ∼6800 Å), corresponding to Mg ii λ2796 absorption redshifts of z ≥ 1.43, feature identification becomes more difficult due to the presence of telluric lines. The strongest lines occur in the A- and B-bands (7600–7630 and 6860–6890 Å, respectively) and also between 7170–7350 Å (e.g., Barlow 2005). These line complexes are generally distinguishable from the Mg ii λλ2796, 2803 doublet. Except within the highly saturated bands themselves, the A- and B-band lines exhibit distinctive patterns of closely spaced pairs. A chance alignment of two A- or two B-band telluric lines at the precise separation of a candidate Mg ii doublet, along with the required doublet ratio, is not common. Nevertheless, these spurious features may obscure weak Mg ii absorption lines or cause some confusion due to line blending. In addition to our objective criteria for identifying Mg ii doublet candidates using the Search algorithm, we were especially diligent to visually inspect candidates in these wavelength regions. Corroboration by visual inspection of the associated Mg i, Fe ii, Ca ii, and Mn ii features mentioned in Section 2.3 was also also performed, but a corroboration was not required for the candidate to be included in the sample. These procedures were adopted in order to maximize the accuracy of the number of absorbing systems located in the spectra.

Using these criteria, we found a total of 480 Mg ii absorption systems, covering the redshift 0.19 ≤ z ≤ 2.55 and having rest-frame equivalent widths over the range 0.006 ≤W_r (2796) ≤ 6.23 Å. These 480 systems were found in 186 of the 249 quasars that we searched.

Once a Mg ii system was confirmed, we used our code Sysanal (see Churchill et al. 1999; Churchill & Vogt 2001) to analyze the absorption. The code first computed the optical-depth median of Mg ii λ2796 profile, i.e., the wavelength (${\lambda }_{\bar{\tau }}$) at which equal integrated optical depth resides to both sides of the profile. By definition, we adopt the system redshift as ${z}_{\mathrm{abs}}={\lambda }_{\bar{\tau }}/2796.352-1$, which is employed to compute the systemic rest-frame velocity zero point of the absorption system. The code then computed the rest-frame equivalent widths, W_r (2796), the Mg ii doublet ratios, DR, the apparent optical depth column densities, N_aod, the kinematic velocity spreads, ω_v , and errors in these quantities for all transitions. The apparent optical depth column density was obtained by inverting the absorption profile to an optical depth profile using the radiative transfer solution ${\tau }_{\lambda }=\mathrm{log}({I}_{\lambda }^{0}/{I}_{\lambda })$, where ${I}_{\lambda }^{0}$ is the fitted continuum, converting the optical depth to column density per unit velocity, and integrating over the profile (Savage & Sembach 1991),

$$\begin{eqnarray}&&{N}_{\mathrm{aod}}=\displaystyle \frac{{m}_{e}{c}^{2}}{\pi {e}^{2}}\displaystyle \frac{1}{f{\lambda }_{0}^{2}}\int {\tau }_{\lambda }\,d\lambda .\end{eqnarray} \tag{ 1 }$$

Lower limits on τ_λ occur when the pixel at λ and its two adjacent pixels meet the condition ${I}_{\lambda }\lt {\sigma }_{{I}_{\lambda }}$, in which case ${\tau }_{\lambda }\geqslant \mathrm{log}({I}_{\lambda }^{0}/{\sigma }_{{I}_{\lambda }})$. If this condition occurs in three contiguous pixels, corresponding to a resolution element, then the profile is considered to be saturated and N_aod is quoted as a lower limit for the apparent optical depth column density.

The kinematic velocity spread is the proportional to the flux-decrement weighted second-moment of the velocity across the Mg ii λ2796 absorption profile (Sembach & Savage 1992),

$$\begin{eqnarray}&&{\omega }_{v}^{2}=\displaystyle \frac{{V}_{(2)}}{{V}_{(0)}}=\displaystyle \frac{1}{{V}_{(0)}}\int {\left(v-\langle v\rangle \right)}^{2}D(v){dv},\end{eqnarray} \tag{ 2 }$$

where D(v) = 1 − I(v)/I⁰(v) is the flux decrement in velocity coordinates, $\langle v\rangle$ is the mean velocity (the flux-decrement-weighted first moment of the velocity), and V₍₀₎ is the velocity "equivalent width" (zeroth moment),

$$\begin{eqnarray}&&\langle v\rangle =\int {vD}(v){dv},\qquad {V}_{(0)}=\int D(v){dv}.\end{eqnarray} \tag{ 3 }$$

The kinematic velocity spread, ω_v , can be interpreted as the equivalent Gaussian standard deviation of the absorption profile. For each associated transition, if absorption is not formally detected at the 3σ significance level at the expected location in the spectrum, the Sysanal code computes the 3σ upper limits on the equivalent widths and apparent optical depth column densities,

As previously mentioned, we do not examine redshift evolution of the absorption properties in this paper. However, as we do aim to present and discuss the distribution of the kinematic properties over the redshift range we surveyed, it is important we establish that we have a reasonably fair sample of Mg ii-selected absorption systems.

It is well-established that Mg ii λ2796 rest-frame equivalent width, W_r (2796), and various kinematic indicators are correlated. For example, the kinematic velocity spread, ω_v , is positively correlated with W_r (2796) (e.g., Churchill et al. 2000; Churchill & Vogt 2001), as is the number of Voigt profile components (e.g., Petitjean & Bergeron 1990; Churchill et al. 2003). Thus, we adopt the premise that a sample of Mg ii-selected absorbers with a fair distribution of Mg ii λ2796 rest-frame equivalent widths would also represent a fair distribution of Mg ii kinematic properties.

Since none of our quasar spectra were selected based upon knowledge of weak absorption, and because no correlation exists between strong and weak absorbers in a given quasar spectrum (Churchill et al. 1999), our sample is unbiased toward weak systems. However, some of the quasar spectra were observed because the presence of strong Mg ii absorption had already been ascertained from previous low-resolution surveys. Given this partial selection bias and the heterogeneous scientific motivations behind the observations of many of these lines of sight (as discussed in Section 2.1), we cannot a priori expect that the strong Mg ii absorption subset is consistent with an unbiased sample. Some statistically founded reassurance of this would allow us to adopt the view that our sample is a fair sample for studying the kinematic aspects of the strong Mg ii absorption systems.

In order to determine whether our strong sample is fair and unbiased, we quantitatively compare the distribution of equivalent widths for W_r (2796 ≥ 0.3 Å from our sample to those from the large blind SDSS surveys of Nestor et al. (2005), featuring 1331 Mg ii absorbers, and of Zhu & Ménard (2013), covering 40,000 Mg ii absorbers. To allow direct comparisons between all three surveys, we limited our analysis to the redshift range of Nestor et al. (2005), which is in common with our survey and that of Zhu & Ménard (2013). This redshift range is 0.34 ≤ z_abs ≤ 2.27, comprising 469 out of the total of 480 absorbers in our survey. For the Nestor et al. (2005) distribution function, we adopted their exponential fit to the function $f(W)=({N}^{* }/{W}^{* })\exp \{-W/{W}^{* }\}$, where N* = 1.187 ± 0.052 and W* = 0.702 ± 0.017 Å, which applies for W_r (2796 ≥ 0.3 Å. For the Zhu & Ménard (2013) distribution function, we adopted their Equation (5) for dN²/dWdz and best-fit parameters from their Table 2. We integrated this function over the adopted redshift range to obtain f(W) = dN/dW. We then normalized both distributions to the area under the distribution for the observed data from our survey.

In Figure 3(a), we present the binned Mg ii λ2796 rest-frame equivalent width distribution normalized to unity at W_r (2796) =0.3 Å. The normalized Nestor et al. (2005) and Zhu & Ménard (2013) distribution functions are superimposed on the data. The shaded regions account for the uncertainties in the fit parameters. Visual inspection would suggest that our sample of strong Mg ii absorbers is populated by an slight overabundance of absorbers having W_r (2796) ∼ 2 Å and perhaps also having W_r (2796) > 4 Å, but that otherwise, within measurement uncertainties, our distribution has a shape generally consistent with an exponential distribution.

Zoom In Zoom Out Reset image size

In Figure 3(b), we plot the cumulative equivalent width distributions. For the observed data, the relative decrement in the range 1 ≤ W_r (2796) ≤ 2 Å reflects our slight overabundance of W_r (2796) ∼ 2 Å equivalent width systems as compared the unbiased surveys. A Kolmogorov–Smirnov (K-S) statistical test was employed in order to quantitatively measure the consistency (or lack thereof) between the distribution of our sample of strong Mg ii systems and the unbiased distributions. We adopt the criterion that the probability for a K-S statistic indicative of the two distributions being inconsistent with one another is P(KS) ≤ 0.0027, corresponding to a 3σ significance level. Compared to the Nestor et al. (2005) distribution, we obtained P(KS) = 0.115. Compared to the Zhu & Ménard (2013) distribution, we obtained P(KS) = 0.076. Thus, both tests indicate the observed distribution cannot be ruled inconsistent with the unbiased surveys. Even though the observed cumulative distributions exhibit some shape variation compared to the unbiased surveys, this variation is significant only at the 1.6σ and 1.8σ levels for the Nestor et al. (2005) and Zhu & Ménard (2013) distributions, respectively.

The minor discrepancy is likely due to some of our quasar lines of sight being observed because a strong absorber, such as a DLA, was targeted for high-resolution analysis. For example, the HIRES spectra contributed by J. X. Prochaska, as well as some of the UVES spectra from the VLT archive (e.g., Jorgenson et al. 2013), were obtained for DLA studies, which are known to exhibit Mg ii absorption with higher equivalent widths (e.g., Rao & Turnshek 2000). The excess at W_r (2796) ≥ 4 Å may also be indicative of lines of sight targeted for very large W_r (2796) systems, perhaps for galactic wind studies (e.g., Bond et al. 2001a, 2001b; Mas-Ribas et al. 2018). As a result of this analysis, we proceed under the assumption that our sample of Mg ii absorption systems has the characteristics of an unbiased sample for the purpose of examining general kinematic properties.

Studies of the kinematics of the absorption systems must account for variations in the signal-to-noise ratio across the velocity window over which the analysis is performed. A varying ratio can result in variable detection sensitivities for very weak components at various velocities locations across a Mg ii λ2796 profile; this could introduce systematics into the distribution of fitted VP component velocities. It is thus imperative we have a controlled sample for our VP and kinematic analysis; we need a uniform detection sensitivity to ensure weak, blended, and high-velocity components above a fixed equivalent width threshold are detectable in all systems included in the analysis. Here, we describe our selection of this final science subsample (comprising 422 systems), which we use for the VP fitting, and analysis of the column density, b parameter, and kinematic distributions of the VP components.

We created a "kinematic sample" by including only those systems for which a 5σ minimum rest-frame equivalent width detection limit was present over a velocity window of ±600 km s⁻¹ relative to the velocity zero point of the Mg ii λ2796 absorption feature. The criterion for a system to be included in the kinematic sample was that either (1) the average limit is $\langle {W}_{r}^{\mathrm{lim}}\rangle \leqslant 0.02\,\mathring{\rm A} ;$ or (2) if $\langle {W}_{r}^{\mathrm{lim}}\rangle$ was greater than 0.02 Å, then $\langle {W}_{r}^{\mathrm{lim}}\rangle$ minus the standard deviation in $\langle {W}_{r}^{\mathrm{lim}}\rangle$ was less than or equal to 0.02 Å. A total of 422 systems, or 88% our complete sample of 480, met the criterion for inclusion into the kinematic sample.

In Figure 4(a), we show the distribution of the system total velocity width for the kinematic sample. The system total velocity width is defined as ${\rm{\Delta }}{v}_{\mathrm{tot}}=\left|{v}_{\max }-{v}_{\min }\right|$, the absolute difference of the "reddest" absorbing pixel and the "bluest" absorbing pixel. The color scheme for the binned data is the same as for Figure 2.

Zoom In Zoom Out Reset image size

The two populations exhibit markedly different distributions. The distribution of Δv_tot for the weak systems can be modeled as a half-Gaussian centered on Δv_tot = 0 km s⁻¹ with a standard deviation of σ(Δv_tot) ≃ 90 km s⁻¹. Weak systems with Δv_tot ≥ 300 km s⁻¹ are very rare; these systems would comprise several weak components highly separated in velocity (see the z_abs = 2.042606 systems along the J101447+430031 line of sight in Figure 6, which has Δv_tot ≃ 150 km s⁻¹). The distribution for the strong systems can be modeled as an asymmetric Gaussian with a mode at Δv_tot ≃ 200 km s⁻¹ and an average standard deviation of σ(Δv_tot) ≃ 200 km s⁻¹. The long tail extends out to Δv_tot ≃ 800 km s⁻¹ due to a small contribution of kinematically extreme systems. In our sample, there are no strong systems with Δv_tot ≤ 50 km s⁻¹.

The total velocity spread measures the extremes of the absorption velocities; it contains no information about the flux decrement distribution. In Figure 4(b), we show the distribution of the kinematic velocity spread, ω_v , also using the same color scheme as Figure 2. Though the kinematic velocity spread contains information about the velocity distribution of the flux decrements, the weak and strong systems exhibit similar distribution shapes as for Δv_tot.

The distribution of ω_v for the weak systems can be modeled as a half-Gaussian centered on ω_v  = 0 km s⁻¹ with a standard deviation of σ(ω_v ) ≃ 15 km s⁻¹. Weak systems with ω_v  ≥80 km s⁻¹ are very rare. The distribution for the strong systems can be modeled as an asymmetric Gaussian with a mode at ω_v  ≃ 35 km s⁻¹ and a standard deviation of σ (ω_v ) ≃50 km s⁻¹. The long tail extends out to ω_v  ≃ 250 km s⁻¹. Note that there are no strong systems with ω_v  ≤ 10 km s⁻¹.

Central to our fitting approach is that we adopt a minimalist approach to the modeling. We fit the absorption systems using as few "clouds," or VP components, as possible by ensuring all components are statistically significant to the χ² statistic. Details of how we ensure statistical significance of all components are given in Section 4.2.

The ionization potentials of the five ions included in the fit are all within a few to several electron volts (eV) of the H i potential of 13.59 eV. The Mg i and Ca ii ionization potentials are both below that of H i at 7.65 eV and 11.87 eV, respectively, while those of Mg ii, Fe ii, and Mn ii are slightly above at 15.04 eV, 16.19 eV, and 15.63 eV, respectively. This bracketing of the H i ionization potential does mean that the line-of-sight ionization structure of a given parcel of gas could be different for the different ions. We assume that any difference is negligible within the context of the resolution and signal-to-noise ratio of the spectra and the general VP decomposition premise of spatially separated isothermal "clouds." The absorption is therefore assumed, for the purpose of VP fitting, to arise all in the same spatial extent of gas and to reflect the same velocity structure.

Another consideration in VP decomposition is how to fit the Doppler b parameter across ions. For our study, we have chosen to let the user input whether to constrain this parameter as 100% "turbulent" or 100% thermal. In the former case, the b parameter of a given "cloud," or VP component, is enforced to be identical for all ions, even though the b parameter is still fit as a free parameter for each component. In the latter case, the b parameter in a given component is enforced to scale in proportion to the inverse of the square root of the atomic mass for each ion. We adopted the default condition of "turbulent" broadening; departures from this are noted in the descriptions of individual systems in Appendix. The systems were fit under the thermal condition if a satisfactory fit could not be achieved using the "turbulent" condition.

In addition, VP modeling by nature requires the assumptions of isothermal "clouds" that occupy distinct line-of-sight locations in velocity space. These latter two assumptions are probably the least defensible, given our developing understanding of the complexity of gas properties as gleaned from "mock absorption line" studies of hydrodynamic cosmological simulations (e.g., Churchill et al. 2015; Peeples et al. 2019). However, Churchill et al. (2015) did find that, in Eulerian adaptive mesh simulations, for low-ionization ions such as Mg ii, the absorbing gas is cloud-like in that gas cells selected by detected absorption are spatially contiguous and comprise a narrow temperature range.

For each absorption system, the first step is to create an initial model of the VP components. We use our graphical interactive program iVPfit, which is an improved version of Profit (Churchill 1997). The output of iVPfit is a complete model of all components for all transitions for all ions included in an absorption system. The component parameters for all transitions of a given ion are "tied." This means that, for a given "cloud," the N, b, and v for that ion are simultaneously constrained by all transitions of that ion and they all have the same value. The interactive process is streamlined by the use of autoscaling of column densities between ions. We found, in the course of this study, that our final models were not sensitive to our "χ by eye" methodology, in which the user tended to be biased toward introducing more components than are statistically significant. For example, we never experienced a case where adding a greater number of components to an initial model changed the final number of components determined to be significant by Minfit.

Using this initial model, we then employ the least-squares fitter Minfit ⁷ (Churchill 1997), which iteratively eliminates all statistically insignificant components and adjusts the remaining components until the least-squares fit is achieved. This entire process, discussed below and illustrated in Figure 5, is based on a series of objective tests and trials. Minfit utilizes the spectral information from all available transitions to constrain these parameters. The velocity of a given component is constrained to be the same across all transitions of all ions. The column density is constrained to be the same across all transitions of a given ion. Finally, the Doppler b parameter can be set by the user to either vary thermally across ions or to be constant across all ions (the "turbulent" option, as discussed in Section 4.1).

Zoom In Zoom Out Reset image size

In a given system, one or more transitions may be compromised over portions of its velocity extent by either bad/noisy pixels or spurious features that have no corroborating absorption in other transitions. In these cases, Minfit allows for one to mask pixels or pixel regions that then contribute no information to the fit.

Finding physically meaningful VP fits in regions of extreme line saturation can be very challenging. For example, consider the velocity region ≃0–75 km s⁻¹ in the z_abs = 1.405187 system in the quasar spectrum of J035405–272421, illustrated in Figure 6, where both the Mg ii λλ2796, 2803 and the majority of the Fe ii transitions are highly saturated. In extreme cases such as this, the least-squares fitting engine can fixate on a local χ² minimum in which the ratio of N(Fe ii)/N(Mg ii) is unphysical, as it is astrophysically rare for [Mg/Fe] to fall an order of magnitude below the solar value. In their VP decomposition of some two dozen Mg ii absorption systems observed with the HIRES spectrograph, Churchill et al. (2003) found the relation

$$\begin{eqnarray}&&\mathrm{log}N(\mathrm{Fe}\,{\rm\small{II}})=1.37\,\mathrm{log}N(\mathrm{Mg}\,{\rm\small{II}})-0.41\end{eqnarray} \tag{ 5 }$$

for the unsaturated velocity regions of the absorption profiles. For our work here, we constrain the Fe ii column densities to obey this relation in highly saturated velocity regions. The user specifies these velocity ranges when the constraint is deemed necessary. We find that this successfully prevents unphysical Fe ii to Mg ii column density ratios in these specified velocity regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (422 panels)

Tarball of figure set images

As illustrated in Figure 5, once the initial model is constructed, any bad pixels are masked, and any saturated velocity regions are specified, Minfit performs a refinement of the model by minimizing the χ² statistic. Essentially, Minfit is a driver that sets up M nonlinear functions in n parameters and feeds the vector of functions to the netlib.org routine dnls1 (More 1978). The M functions are the individual terms of the χ² statistic, one for each pixel for all transitions, and the n parameters are the VP free parameters, where the value of n depends on the number of components, transitions, and ions. The least-squares fit (the box labeled "LSF" in Figure 5) is performed by the routine dnls1, a modification of the Levenberg–Marquardt algorithm. Two of its main characteristics involve the proper use of implicitly scaled parameters and an optimal choice for the correction terms. The routine approximates the Jacobian by forward differencing. The use of implicitly scaled parameters achieves scale invariance and limits the size of the correction in any direction where the functions are changing rapidly. The optimal choice of the correction guarantees (under reasonable conditions) global convergence from starting points far from the initial guess solution and a fast rate of convergence for problems with small residuals. Further details, including the methods for computing the uncertainties in the fitted parameters, are described in Churchill (1997).

Following the LSF, Minfit computes the errors in the VP parameters, σ_N , σ_b , and σ_z (recall that each component is fitted in redshift space and later converted to rest-frame velocity). The robustness of components is examined in two ways. First, for each ion, each component, i, is checked against its nearest neighbor components, i − 1 and i + 1, for redshift overlap using the conditions $\left|({z}_{i}-{z}_{i-1})/{\sigma }_{{z}_{i}}\right|\gt 1$ or $\left|({z}_{i}-{z}_{i+1})/{\sigma }_{{z}_{i}}\right|\gt 1$. If a component satisfies one of those conditions, it is flagged for significance testing, as described below. For the second robustness check, a "badness" parameter is computed for each component,

$$\begin{eqnarray}&&\mathrm{badness}={\left[{\left(\displaystyle \frac{{\sigma }_{N}}{N}\right)}^{2}+{\left(\displaystyle \frac{{\sigma }_{b}}{b}\right)}^{2}+{\left(\displaystyle \frac{{\sigma }_{z}}{z}\right)}^{2}\right]}^{1/2}.\end{eqnarray} \tag{ 6 }$$

If any components have a badness parameter that exceeds a user-specified value, and/or any components are identified to overlap in redshift, then a series of significance tests are instigated. We have adopted a default value of badness = 1.5, but have relaxed this parameter as needed for various systems (typically those that exhibit extreme saturation). If no VP components exceed the badness threshold and none overlap their nearest neighbor, the current VP model is adopted.

Components flagged for significance checking are rank-ordered, with priority given to redshift overlap followed by a ranking from highest badness to lowest badness. The highest-ranked flagged component is simply deleted from the "current" VP model, and an LSF is obtained for this new "test" model, which has one fewer VP component. Accounting for the different degrees of freedom in the two models, an F-test is performed between the "current" and "test" models to determine whether inclusion of the flagged component provides a statistically significant improvement in the χ² statistic to a user-specified confidence level. We adopt a default confidence level of 97%, though we have relaxed this number for selected systems as needed. No systems have been fitted for which the confidence level of the components is below 90%.

If the tested component was statistically significant, Minfit retains the "current" VP model; it then proceeds to investigate the statistical significance of the next line in the sorted "bad" component array, and so on until the array is exhausted. If the tested component was not statistically significant, then the "test" model is adopted as the "current" model; the number of components required to model the absorption system is now reduced. Minfit then identifies whether any neighboring components exhibit redshift overlap, computes the badness parameters for this newly adopted "current" VP model, and proceeds to test the significance of the components. The entire process is automated and objective. A final VP model is adopted when all components are determined to be statistically significant. Human manipulation of the objective outcome of the process (i.e., the adopted VP model) can occur only via modification of the threshold badness parameter and/or the threshold confidence level. Those systems for which we modified the default values are noted in the descriptions of individual systems in Appendix.

Once the final VP model is adopted and the final uncertainties have been computed, upper limits are computed for the component column densities for associated ions (i.e., Mg i, Fe ii, Mn ii, Ca ii) for which no transitions were formally detected at the velocity positions of Mg ii components. The Doppler b parameter of the specific Mg ii component is adopted and the column density limit is determined from the 3σ equivalent width detection limit across the velocity range of the component

We note that Minfit is a deterministic algorithm, in that, for a given input model and user-specified parameters (masking, saturated regions, badness, and confidence level), the computational path of the least-squares fit and the final solution will always be the same. By exploring the outcomes as a function of variations in the input model and/or the user-specified parameters, we inspected various final models. In the end, the adopted final model for a given system is a human decision (however, see Bainbridge & Webb 2017a).

In Figure 7(a) we show the binned distributions of the number of VP components for all systems, weak systems (W_r (2796) < 0.3 Å), and strong systems (W_r (2796) ≥ 0.3 Å); the color scheme of the histograms is the same as used in Figure 2. The full sample was modeled with an average of $\langle {N}_{\mathrm{VP}}\rangle =7.1$ components, whereas the weak systems were modeled with an average of $\langle {N}_{\mathrm{VP}}\rangle =2.7$ components, and the strong systems with $\langle {N}_{\mathrm{VP}}\rangle =10.3$ components.

Zoom In Zoom Out Reset image size

Churchill (1997) modeled simulated multicomponent Mg ii absorption systems in synthetic spectra having the characteristics of HIRES/Keck spectra and found that, on average, ∼30% of VP components were not recovered using Minfit. Those simulated profiles were generated using the observed distributions of VP component velocity separations, column densities, and Doppler b parameters from the VP decomposition of two dozen systems observed with Keck/HIRES, but the number of components was fixed at N_VP = 10 as a control condition. The signal-to-noise ratio for a given simulation was also held fixed; the quoted results here are for a simulation such that the 5σ equivalent width detection threshold was 0.02 Å. The average number of components recovered in these experiments did vary with signal-to-noise ratio and the presence (or nonpresence) of associated transitions with clear kinematic structure. Component recovery improved to ∼90% with the presence of associated transitions, and always decreased as signal-to-noise ratio decreased. Though the tests are based on the assumption that Mg ii absorption profiles are a complex of VP components (a clearly simplistic scenario), if the outcomes can be applied directly to our sample, it would suggest that the actual mean numbers of components are ∼11%–43% higher than what we report here.

In Figure 7(b), the number of VP components is plotted as a function of the Mg ii λ2796 rest-frame system equivalent width. Blue points represent weak systems and pink points represent strong systems. A linear fit to the full sample resulted in a slope of 8.62 ± 0.23 clouds Å⁻¹, which can be interpreted as the VP component line density. Note the increased scatter for W_r (2796) > 2 Å, where VP fitting can become challenging and problematic for highly saturated or partially saturated absorption profiles. For W_r (2796) ≤ 2 Å, the absolute standard deviation of N_VP about the fitted relation is 3.2 "clouds."

Previous work examined the inverse of the VP component density, i.e., the slope in terms of Å cloud⁻¹. For our sample, our fit corresponds to 0.116 ± 0.003 Å cloud⁻¹. The slope found by Churchill (1997) for a sample of 36 Mg ii systems observed with Keck/HIRES (the same spectral resolution as this study) was 0.076 ± 0.004 Å cloud⁻¹. Their higher number of components per unit equivalent width (13.2 clouds Å⁻¹) is likely due to the modifications we made to the Minfit program. The distribution of the signal-to-noise ratio of the spectra clearly plays a role, as higher-quality data can constrain the VP models to have a larger number of "clouds"; however, the distributions of our survey and that of Churchill (1997) are statistically consistent. In Churchill (1997), the "badness" of only Mg ii components was tested for significance, whereas for this work, all components from all associated ions were also tested; this resulted in a reduction of the number of components required for the final VP models. In a survey of Mg ii absorbers in moderate-resolution spectra (∼30 km s⁻¹), Petitjean & Bergeron (1990) found a linear relationship with slope 0.35 Å cloud⁻¹, corresponding to ≃3 clouds Å⁻¹.

Overall, we see that the measured VP component line density is strongly affected by the spectral resolution, the signal-to-noise ratio, and the VP fitting approach. It is expected that the higher signal-to-noise ratios and resolutions of the future 30 m class telescopes will result in an even higher component line density. As such, if any future works undertake a characterization of the component line density, we recommend adopting the approach of enforcing the minimum number of statistically significant components to a well-defined confidence level.

The VP component column densities are a key input to photoionization models, which are commonly employed to constrain "cloud" ionization conditions and metallicities, and to explore the spectral energy distribution of the local ionizing radiation field (e.g., Werk et al. 2014; Lehner et al. 2019; Pointon et al. 2019). The distribution of column densities also provides key constraints for hydrodynamic cosmological simulations of the circumgalactic and intergalactic medium (e.g., Churchill et al. 2015; Oppenheimer et al. 2018; Peeples et al. 2019). The column density distribution obtained from VP decomposition can more effectively account for unresolved saturation in the absorption profiles than direct profile inversion via the apparent optical depth method (Savage & Sembach 1991; Jenkins 1996). Furthermore, VP decomposition provides an explicit and well-defined segregation of absorbing components.

VP analysis of our full sample of Mg ii systems yields the component column density distribution shown in Figure 8(a). The 422 Mg ii systems comprise a total of 2989 VP components. The distribution has been normalized by this total number of components, so that the quantity f(N) represents the fraction of VP components per unit column density in the sample, a quantity that is reproducible in any survey independent of the number of quasar spectra, absorption line systems, and/or the redshift path sensitivity function of the survey. The pink vertical shaded area indicates the region of partial completeness due to line blending in kinematically complex absorption profiles as determined by the simulations of Churchill et al. (2003). They found that the 90% completeness levels for unblended and blended lines were log N(Mg ii) = 11.6 cm⁻² and 12.4 cm⁻², respectively, for spectra having a mean 5σ equivalent width sensitivity of W_r  = 0.02 Å. "Completeness level" refers to the percentage of simulated components of a given column density recovered during VP analysis. In the shaded region on Figure 8(a), components in complex profiles can be lost due to blending, though the completeness for unblended components is 90%. Below this region, even unblended (single component) absorbers can be lost due to the signal-to-noise ratio of the spectra.

Zoom In Zoom Out Reset image size

The column density distribution can be fit by a power law,

$$\begin{eqnarray}&&f(N)={{CN}}^{-\delta },\end{eqnarray} \tag{ 7 }$$

where f(N) is the fraction of clouds with column density N per unit column density, C is a normalization constant, and δ is the power-law slope. The maximum likelihood method was used to obtain the power-law fit to the unbinned data (see Churchill 1997). The fit was performed only on column densities above the region of partial completeness (log N(Mg ii) > 12.4 cm⁻²) so as to not skew the slope. We obtained δ = 1.45 ± 0.01. In a study of 14 Mg ii systems containing 33 VP components, Petitjean & Bergeron (1990) obtained a significantly shallower slope of δ = 1.0 ± 0.1, although their spectral resolution was lower (≃30 km s⁻¹). In a study with identical resolution, Churchill et al. (2003) obtained δ = 1.59 ± 0.05 by fitting their sample of 175 VP components in 23 Mg ii systems. Our slightly shallower slope is likely due to the modifications we made to Minfit, which resulted in a lower VP component line density (see Section 5.1) skewed slightly toward higher column density components.

In Figure 8(b), we plot the Mg ii Doppler b parameter distribution of the VP components. The median Doppler b parameters and standard deviations are 4.5 ± 3.5 km s⁻¹, 6.0 ± 4.5 km s⁻¹, and 5.7 ± 4.4 km s⁻¹ for the weak (blue), strong (pink), and full (gray) samples, respectively. Churchill (1997) found a median Doppler parameter of ∼3.5 km s⁻¹ in a study of 48 Mg ii systems and Churchill et al. (2003) found 5.4 ± 4.3 in a study of 23 Mg ii systems; in both cases, the data were of comparable quality and resolution (6.6 km s⁻¹). Petitjean & Bergeron (1990) found, in their study of 14 Mg ii systems, that the b distribution peaked between 10 and 15 km s⁻¹; however, their larger value was because their data were of lower spectral resolution (30 km s⁻¹) and they noted that this significantly distorted the observed distribution. Based on simulations designed to test the recovery of the Doppler b distribution in HIRES spectra (Churchill et al. 2003), the observed distribution peak is likely ∼1–2 km s⁻¹ too high relative to the true underlying distribution (assuming Mg ii absorption profiles arise in spatially separated isothermal clouds giving rise to Voigt profiles). In addition, the distribution tail at high b values has been shown in these simulations to be an artifact of component blending and unresolved saturation. As mentioned in Section 5.1, ∼30% of simulated components are not recovered in the VP decomposition for our spectral quality. As a result, some b parameters in the observed distribution are too broad compared to the "true" distribution.

If the Doppler broadening is assumed to be predominantly thermal, then the observed b parameter distribution medians correspond to gas temperatures of ∼30,000 K, ∼53,000 K, and ∼47,000 K for the weak, strong, and full samples, respectively. However, applying the 1–2 km s⁻¹ correction to the mode of the b parameter distribution would produce a median temperature, in the case of the weak sample, of ∼9000–18,000 K, in the strong sample, of ∼23,000–37,000 K, and for the full sample, a median temperature of ∼20,000–32,000 K.

Based on simulations designed to test the recovery of the Doppler b distribution in HIRES spectra (Churchill et al. 2003), the observed distribution peak is likely ∼1–2 km s⁻¹ too high relative to the true underlying distribution (assuming Mg ii absorption profiles arise in spatially separated isothermal clouds giving rise to Voigt profiles). In addition, the distribution tail at high b values has been shown in these simulations to be an artifact of component blending and unresolved saturation. As mentioned in Section 5.1, ∼30% of simulated components are not recovered in the VP decomposition for our spectral quality. As a result, some b parameters in the observed distribution are too broad compared to the "true" distribution.

If the Doppler broadening is assumed to be predominantly thermal, then the observed b parameter distribution medians correspond to gas temperatures of ∼30,000 K, ∼53,000 K, and ∼47,000 K for the weak, strong, and full samples, respectively. However, applying the 1–2 km s⁻¹ correction to the mode of the b parameter distribution would produce a median temperature, in the case of the weak sample, of ∼9000–18,000 K, in the strong sample, of ∼23,000–37,000 K, and for the full sample, a median temperature of ∼20,000–32,000 K.

The velocity clustering of the VP components is quantified using the velocity two-point correlation function (TPCF; Petitjean & Bergeron 1990). The TPCF is the probability, P(Δv), that any randomly selected pair of VP components within a system will have a velocity separation Δv. The velocity TPCF for our full sample is shown in Figure 9 and tabulated in Table 5. This probability distribution is commonly fit with a composite Gaussian function,

$$\begin{eqnarray}&&P({\rm{\Delta }}v)=\displaystyle \sum _{n=1}^{N}{g}_{n}({\rm{\Delta }}v),\end{eqnarray} \tag{ 8 }$$

where N is the number of components,

$$\begin{eqnarray}&&{g}_{n}({\rm{\Delta }}v)=\displaystyle \frac{1}{\sqrt{2\pi }}\displaystyle \frac{{a}_{n}}{{\sigma }_{n}}\exp \left\{-\displaystyle \frac{{\left({\rm{\Delta }}v\right)}^{2}}{2{\sigma }_{n}^{2}}\right\},\end{eqnarray} \tag{ 9 }$$

are the Gaussian function components, and where a_n is a fitted scaling factor and σ_n is a fitted velocity dispersion. The amplitude of each component is then

$$\begin{eqnarray}&&{A}_{n}=\displaystyle \frac{1}{\sqrt{2\pi }}\displaystyle \frac{{a}_{n}}{{\sigma }_{n}}.\end{eqnarray} \tag{ 10 }$$

Petitjean & Bergeron (1990) and Churchill et al. (2003), in their studies of 14 and 23 Mg ii absorption systems, respectively, fit their TPCF distributions of VP components using two-component Gaussian models. The results of Petitjean & Bergeron (1990) were σ₁ = 80 km s⁻¹ and σ₂ = 390 km s⁻¹ for a spectral resolution of 30 km s⁻¹. Those authors attributed the narrower of these two components to motions within galaxy halos, and the broader to the motions of galaxy pairs. However, the Churchill et al. (2003) study, which had 6.6 km s⁻¹ velocity resolution, calculated a best-fit dispersion of σ₁ = 54 km s⁻¹ and σ₂ = 166 km s⁻¹. They suggested that the Mg ii component velocity dispersion might reflect the range observed in face-on galaxy disks and edge-on galaxy rotational motions, as well as infall and outflow in the halos. The larger TPCF dispersions reported by Petitjean & Bergeron were likely due to resolution effects that prevented identification of smaller VP component velocity splittings.

Zoom In Zoom Out Reset image size

We fitted the TPCF from our full sample with a three-component composite Gaussian function; the resulting functions are superimposed on the TPCF in Figure 9. Three components were used because two did not adequately fit the extended tail of our distribution. Our resulting velocity dispersions are σ₁ = 25.4 km s⁻¹, σ₂ = 68.7 km s⁻¹, and σ₃ =207.1 km s⁻¹. Velocity separations of Δv < 10 km s⁻¹ were excluded from the fit because their relative numbers are artificially lowered due to component blending at small Δv.

Small velocity separations are the most probable. The probability drops steeply up to Δv ∼ 150 km s⁻¹; at larger separations, the probability decreases more slowly to our maximum velocity separation of Δv = 572.8 km s⁻¹. Considering the modern view of the kinematics of the low-ionization circumgalactic medium (e.g., Weiner et al. 2009; Kacprzak et al. 2010; Martin et al. 2012; Nielsen et al. 2015, 2016; Ho et al. 2017; Kacprzak 2017; Zabl et al. 2019, 2020), it would be a gross overinterpretation of the TPCF parameterization to identify each Gaussian component with a specific galactic kinematic component or physical phenomenon related to the circumgalactic baryon cycle. Parameterizing the TPCF by a functional fit provides a convenient functional characterization of VP component velocity clustering. One should be cautious to consider the equivalent width detection threshold and the spectral resolution when comparing the TPCF. We remind the reader that we applied a uniform detection threshold criterion for an absorbing system to be included in our analysis of kinematics (see Section 3.2). This criterion ensures that the kinematics analysis has a uniform sensitivity to small equivalent width absorption features at high velocities from system to system.