Evidence for Late-time Feedback from the Discovery of Multiphase Gas in a Massive Elliptical at z = 0.4

A visual inspection of our HST/COS spectrum of HE 0047−1756A reveals the presence of numerous absorption features consistent with the Lyman and Werner bands at redshift , or approximately −430 from the systemic redshift of the lens galaxy (see Figure 2).⁵ The absorption features coincide in velocity with the strongest low-ionization absorption component identified in Z16 (component 1 in their Table 6). We are able to identify more than 140 absorption transitions originating from the ground state of at different rotational levels from J = 0 to J = 5. Each of these transitions has a vibrational quantum number of for the lower state and () for the upper state in the Lyman (Werner) band.

Zoom In Zoom Out Reset image size

To characterize the molecular gas properties, we perform a Voigt profile analysis using custom software that models the observed transitions in each J level simultaneously. We adopt the line list from Ubachs et al. (2019), which was made available to us by P. Petitjean (2020, private communication). Although we detect more than 140 transitions, a significant fraction of these lines are blended with each other or other absorption lines. To ensure robust fitting results, we perform our absorption analysis on a subset of available lines (between five and 14 transitions) for each J value, which are selected to contain minimal blending and have unambiguous local continuum level. While only a fraction of observed transitions are used to find the best-fit model, we find that the resulting full absorption model reproduces the absorption profiles of most of the excluded transitions reasonably well.

For each rotational J level, we first generate a model spectrum for a single-component line profile, motivated by both the narrow linewidths and lack of kinematic substructures in the observed absorption profiles (see Figure 2). The Voigt profile is uniquely defined by three free parameters: the line centroid redshift z_abs, the absorption column density , and the Doppler parameter b. To reduce the number of free parameters, all transitions from a given J level are tied to have the same and b. We further require different J levels to share the same line centroid redshift. Once a theoretical absorption spectrum has been generated, it is convolved with the relevant COS LSF and subsequently binned to match the pixel resolution of the data. Finally, this model spectrum is compared to the data, and the best-fit model parameters for each J level are found by minimizing values at the selected transitions. We estimate the model uncertainties by constructing a marginalized posterior probability distribution for each model parameter, based on a Markov Chain Monte Carlo (MCMC) analysis done with the Emcee package (Foreman-Mackey et al. 2013). Each MCMC run consists of 500 steps performed by an ensemble of 250 walkers, which are seeded in a small region of the parameter space around the minimum solution to speed up convergence.

We present the best-fit model absorption profiles and compare them to the data in Figure 2. In addition, we summarize the results of the Voigt profile analysis in Table 1, where we report the model parameters and estimated 68% confidence intervals for the J = 0 to J = 5 levels. For J = 6, which does not exhibit any detectable absorption, we report in Table 1 the 95% upper limit on the absorption column density for a b = 10 line profile (matching the linewidth of the J = 5 level), estimated using the error array at the strongest available J = 6 transition in the COS data. The best-fit model yields a total column density of and a best-fit redshift of . The centroid of the line profile is consistent within uncertainties ( ) with the strongest low-ionization metal component identified in ground-based optical echelle spectra (Z16), which indicates their association.

Table 1.  Properties at d = 4.6 kpc from the Lens Galaxy

zabs	J	log N/	b
			()
0.405985	0	17.34+0.09−0.26	2.7+0.8−0.5
	1	17.58+0.07−0.21	3.9+0.5−0.3
	2	15.71+0.44−0.18	6.3+0.8−1.1
	3	15.86+0.28−0.23	6.2+1.0−0.8
	4	14.82+0.06−0.04	8.7+1.6−1.7
	5	14.61+0.06−0.05	10.5+5.0−2.5
	6	<14.4a	10
	Total	17.8+0.1−0.3

Note.

^a95% upper limit (see Section 3.1).

Download table as:  ASCIITypeset image

As shown in Table 1, our analysis also identifies a trend of increasing Doppler parameter with increasing J value, from at J = 0 to at J = 5. The trend of rising velocity dispersion with rotational level has been reported in a number of -bearing damped Lyα absorbers (DLAs) at low and high redshifts (e.g., Ledoux et al. 2003; Albornoz Vásquez et al. 2014; Boettcher et al. 2020). The kinetic temperature needed to thermally broaden the line profiles to a linewidth of is K, significantly higher than the temperatures at which a significant amount of molecular gas is expected to be present. Therefore, our measurements indicate that nonthermal line broadening is dominant for both low- and high-J states in the gas, with increasing turbulence toward higher rotational states. We discuss the possible origins of this trend in Section 4.1.

With both the neutral and molecular hydrogen contents of the absorber known, we can estimate the molecular gas fraction according to the following expression, . The extremely strong low-ionization absorber observed along HE 0047−1756A is resolved into 15 kinematic components (Figure 3;Z16). While the total N(H i) can be measured robustly from the strong Lyα damping wings (Figure 1), it is impossible to constrain the H i column densities of these individual components because all available H i Lyman series lines are heavily saturated and the different components blended. Thus, we first estimate by attributing all the observed N(H i) to the -bearing component. Although this assumption is unrealistic because it would result in highly asymmetric Lyα damping wings, owing to the -bearing component occurring at the blue extremum of the profile, it yields a conservative lower limit on of log(. To estimate an upper bound on the molecular gas fraction, we note that the -bearing component contains % of the total column densities of the low-ionization species probed by Mg i, Mg ii, and Fe ii absorptions (see Table 6 of Z16). If we assume that all 15 components have similar metallicities and dust content, which is justified by the relatively uniform elemental abundance ratio observed across all components (Z16), then the inferred N(H i) of the -bearing component is log N(H i). Consequently, the implied molecular gas fraction is log(.

Zoom In Zoom Out Reset image size

The observed at from the lens galaxy is comparable to nearby ellipticals with CO detections (Welch et al. 2010; Young et al. 2014) but is among the highest known for DLAs, where of absorbers with log N(H i) have log (e.g., Crighton et al. 2013; Muzahid et al. 2015b, 2016; but see Boettcher et al. 2020). Considering that forms on the surface of dust grains, the high can be explained by the high gas metallicity, (Z17), which results in an elevated dust-to-gas ratio relative to the general DLA population.

Z17 previously noted possible absorption features from high-ionization metal lines associated with the lens galaxy along both sightlines of HE 0047−1756. However, the low spectral resolution of the Z17 data prevented a detailed investigation into this tentative detection of high-ionization gas. The HST/COS spectrum of HE 0047−1756A clearly resolves different metal absorption profiles, enabling precise measurements of gas kinematics and column densities of the highly ionized species.

As shown in the top three panels of Figure 3, the new COS spectrum confirms that O vi absorption is indeed detected and resolved into multiple components in the lens galaxy. In addition, N v absorption is also detected with a kinematic structure that is consistent with O vi. To constrain their absorption properties, we perform a joint Voigt profile analysis of the O vi line and the N v doublet, following the method of Zahedy et al. (2019).⁶ To ensure the robustness of the fit, we tie both the component structure and Doppler linewidths of the two ions. We summarize the results from our Voigt profile analysis of these high-ionization species in Table 2. The continuum-normalized absorption profiles and best-fit models for O vi and N v are presented in the top three panels of Figure 3. To compare the kinematics between low- and high-ionization species, we also show the observed and modeled absorption profiles of Mg i, Mg ii, and Fe ii in the bottom five panels of Figure 3, from the previous absorption analysis reported in Z16.

Table 2.  High-ionization Absorption Properties at d = 4.6 kpc

Component	Species	dvca	b	log N/
		(km s−1)	(km s−1)
1	O vi	+32.0 ± 3.4	51.4 ± 4.9	14.31 ± 0.03
	N v			13.75 ± 0.11
2	O vi	+220.1 ± 5.6	47.9 ± 8.9	14.02 ± 0.06
	N v			13.58 ± 0.16
3	O vi	+372.3 ± 6.6	42.8 ± 3.5	14.71 ± 0.06
	N v			13.76 ± 0.15
4	O vi	+423.4 ± 14.7	64.2 ± 21.4	14.60 ± 0.13
	N v			14.30 ± 0.06
5	O vi	+507.3 ± 12.7	43.4 ± 8.4	14.13 ± 0.18
	N v			13.59 ± 0.18
6	O vi	+672.3 ± 3.0	25.5 ± 4.6	14.10 ± 0.06
	N v			13.33 ± 0.22

Note.

^aRelative velocity shift from the H₂ absorption redshift, z_abs = 0.405985.

Download table as:  ASCIITypeset image

It is clear from Figure 3 that the high ions exhibit a distinct kinematic structure from that of the low ions, which indicates that the high ions arise in a different gas phase (Zahedy et al. 2019). Specifically, our analysis reveals a highly ionized gas phase in the lens ISM that is kinematically complex, comprising six broad kinematic components () that span in line-of-sight velocity. The observed total column densities of these highly ionized species are log N(O vi) and log N(N v). These O vi and N v absorbers are among the strongest known to be in the vicinity of galaxies (see Johnson et al. 2015; Muzahid et al. 2015a; Werk et al. 2016; Rosenwasser et al. 2018; Zahedy et al. 2019), where high-ionization absorbers with log N(O vi) and log N(N v) are rare.

It is also interesting to note the observed N v to O vi column density ratios among the six high-ionization components, which have an arithmetic mean and dispersion of log (N v)(O vi) . These ionic ratios are considerably higher than typical values seen in the Galactic corona (e.g., Wakker et al. 2012), the CGM of external galaxies (e.g., Werk et al. 2016; Zahedy et al. 2019), and the high-redshift intergalactic medium (e.g., Lehner et al. 2014), where a large majority of absorbers in these diverse environments exhibit log N(N v)(O vi) . We argue that a supersolar in the highly ionized gas phase is the most likely explanation, considering that high nitrogen-to-alpha ratios of have been reported in the evolved stellar populations and cool ISM of nearby ellipticals (e.g., Greene et al. 2013; Yan 2018). Similar ratios in both high- and low-ionization gases would also suggest a causal link between different phases of the ISM of the elliptical lens (we discuss this connection in Section 4.3).

The distribution of molecules among different rotational levels reflects the excitation state of the gas and offers insight into the physical mechanisms that are responsible. The relative populations of different rotational levels can be described by a Boltzmann distribution, , where N_J is the column density for the rotational level J, is the excitation temperature from J = 0 to rotational level J, and B_v = 85.36 K. The statistical weight g_J is for even-numbered J or for odd-numbered J. In Figure 4, we show the excitation diagram of the absorber for rotational states between J = 0 and J = 5. The column density ratio between J = 0 and 1 states, which contain % of the total , implies an excitation temperature of K.

Zoom In Zoom Out Reset image size

The observed for the bulk of molecular gas along the lensed sightline is comparable to typical values reported in -bearing DLAs (see Muzahid et al. 2015b). However, while the observed population for J = 2 can also be reproduced well by the same excitation temperature, this single-temperature model fails to explain the observed column densities at (Figure 4, left panel). The predicted column densities for these higher rotational states are orders of magnitude lower than the observations, which indicates a higher excitation temperature for .

It is well-known from Galactic studies that a one-temperature fit typically works only for optically thin absorbers with log (e.g., Spitzer et al. 1974; Spitzer & Jenkins 1975; Jenkins & Peimbert 1997). In contrast, stronger absorbers in the Galaxy and beyond have been found to exhibit elevated populations at higher rotational levels (e.g., Jenkins & Peimbert 1997; Reimers et al. 2003; Noterdaeme et al. 2007; Rawlins et al. 2018; Balashev et al. 2019; Boettcher et al. 2020), which indicates that the gas is bifurcated into two excitation temperatures. Motivated by these prior observations, we perform a simultaneous fit of a two-temperature model to our data and find that the observed level populations are reproduced well by two excitation temperatures of K and K (Figure 4, middle panel).

The observed temperature bifurcation and trend of rising velocity dispersion with J level (see Section 3.1) can be understood to be a consequence of the absorption originating in a gas cloud with an internal density and/or temperature stratifications (e.g., Noterdaeme et al. 2007). In this scenario, most of the column densities at low-J levels originate from the inner layer of the cloud, where the gas is sufficiently dense and shielded from radiation that collisions are the dominant excitation mechanism. Consequently, the low-level populations are essentially thermalized and the lower excitation temperature is highly coupled to the kinetic temperature of the gas.

In contrast, the elevated column densities and broader line profiles of high-J levels indicate that they arise primarily from warmer and more turbulent outer layers of the cloud (e.g., Lacour et al. 2005). At these locations, molecules can be highly excited through collisions triggered by shocks and turbulent dissipation (e.g., Jenkins & Peimbert 1997; Gredel et al. 2002; Gry et al. 2002; Ingalls et al. 2011), as well as through radiation pumping by an external UV radiation field (e.g., Jura 1975; Klimenko & Balashev 2020). A unique prediction of the shock scenario is a systematic shift of up to a few in line centroids with increasing rotational state, which is caused by the different J levels originating from distinct locations moving at slightly different speeds relative to the shock front (e.g., Jenkins & Peimbert 1997; Gredel et al. 2002). Although we do not detect any systematic shift in line centroids with J levels to within the precision of our COS wavelength calibration ( ), we cannot rule out a more modest shift of or less, which may be the result of weaker shocks (e.g., Gredel 1997).

As an alternative, we now explore radiation pumping as an excitation mechanism. We perform a series of calculations using the Cloudy v.13.03 code (Ferland et al. 2013) to simulate a plane-parallel slab of gas with uniform density n_H that is irradiated by two UV radiation fields: the updated Haardt & Madau (2001) extragalactic UV background at z = 0.4, known as HM05 in Cloudy, and the built-in unextinguished Milky Way ISM radiation field from Black (1987). To constrain the strength of UV radiation that is required to reproduce the observations, we vary the overall intensity of the ISM radiation field by a scale factor of between 0.1 and 100. We incorporate dust grains in the calculations following the observed grain abundance and size distribution in the local ISM. For each input radiation field, we construct a grid of Cloudy models spanning a wide range of gas densities () at the observed gas metallicity (Z17). For each grid point, Cloudy calculates the expected column density for each J level assuming thermal and ionization equilibrium. To simulate two-sided illumination of the cloud, we use half the observed as the stopping condition for the calculations and subsequently double the output level populations for comparison with the data.

We summarize the results of our Cloudy calculations in the right panel of Figure 4, where the set of models that best reproduce the observed excitation diagram are shown in thin blue curves. These models have UV radiation fields that are 15–25 times stronger than the local ISM radiation field. The gas densities have a range of , with mean kinetic temperatures (90–130 K) and H i column densities (log N(H i)) that are broadly consistent with the observations. While it is clear that these simple models are only able to roughly reproduce the general trend seen at , this exercise demonstrates that if the elevated populations at higher rotational levels are primarily due to radiation pumping, the required UV radiation field is significantly higher than what is observed in the Galactic ISM (see also, e.g., Boettcher et al. 2020; Klimenko & Balashev 2020).

A benefit of using multiply lensed images of a QSO system as gas probes is the ability to investigate spatial variations in the gas properties of a foreground galaxy. As described in Z16 (see their Figure 1), the doubly lensed images of HE 0047−1756 probe opposite sides of the massive elliptical lens galaxy, with sightline A at (1.8 half-light radii, r_e) and sightline B at (1.3 r_e). The observed H i and low-ionization metal column densities differ by less than 0.1–0.2 dex between the two sightlines (Z16; Z17), despite a separation of in projection. These similarities suggest that the cool (T ∼ 10⁴ K) ISM phase is spatially extended, with a high gas covering fraction at .

While we are unable to perform a detailed analysis on the O vi absorption detected in the low-resolution STIS FUV spectrum of HE 0047−1756B (Z17), we can compare the general absorption properties of the O vi absorbers detected along the two sightlines. Specifically, the total O vi rest-frame equivalent width is Å along sightline B, which is very similar to what is observed in the COS spectrum along sightline A, . Furthermore, the observed FWHM of the O vi profile along sightline B is , which is comparable to the observed kinematic spread of for the O vi absorption profile along sightline A (Section 3.2). The coherent O vi absorption properties between the two sightlines imply that, similar to the low-ionization gas phase, the highly ionized ISM is spatially extended and has a high covering fraction on a scale of ∼5 kpc in the massive elliptical.

The lack of a high-resolution FUV spectrum of HE 0047−1756B prevents a direct search for along this sightline. To assess whether we can constrain spatial variations in molecular gas properties using the available low-resolution STIS FUV spectrum of sightline B, we perform the following experiment. First, we divide the best-fit model for the full Lyman and Werner bands from the high-resolution COS FUV spectrum of sightline A, to remove all absorption from the spectrum. Next, we convolve the resulting "-free" spectrum of the QSO with the STIS LSF and compare the result with our STIS spectrum of sightline A. We find that, while individual lines are unresolved in the STIS spectrum, the combined absorptions from the Lyman and Werner bands result in an overall flux decrement that is detectable across the QSO spectrum.

Motivated by the result of the experiment, we generate a series of model Lyman and Werner bands spanning a wide range of , apply them to the "-free" QSO spectrum, and convolve the results with the STIS LSF. Each of the resulting spectra is then rescaled to the level of sightline B using the mean observed flux ratio of the two lensed images in two absorption-free regions: 1575–1585 and 1640–1650 Å in the observed frame. Finally, we compare the products to the STIS spectrum of sightline B in the spectral region between 1415 and 1435 Å in the observed wavelength, which has a large concentration of strong transitions, and infer the allowed using a analysis. The observed spectrum of HE 0047−1756B is consistent with the presence of with at the 95% confidence level.

The inferred molecular gas fraction of along sightline B is a factor of at least times lower than that observed along sightline A on the opposite side of the galaxy, (Section 3.1). This exercise suggests that, in contrast to the neutral and highly ionized gas phases, the molecular gas distribution in the lens ISM is clumpier. Furthermore, the observed along the two lensed sightlines are consistent with nearby quiescent galaxies found to harbor molecular gas (e.g., Young et al. 2014) but low compared to typical values in star-forming disks (Chen 2017b and references therein). If these constraints are representative of the rest of the galaxy, they imply a low mass fraction of dense, cold molecular gas in the multiphase ISM of the lens.

How the ISM is partitioned by mass into its different gas phases depends sensitively on the gas cooling rate, the heating available to offset this cooling, and the relevant timescales of these processes. The simultaneous detections of multiple gas phases in the lens ISM enable such an investigation for the first time in a distant elliptical, which can offer valuable insight into late-time feedback in massive elliptical galaxies.

Specifically, now that we have robustly detected highly ionized gas in the lens galaxy and constrained its properties, we can calculate the mass budget in the warm (T ∼ 10⁵ K) gas phase and compare it to the previously estimated mass budget in the cool ISM (Z17). For the cool phase, Z17 estimated a total Fe mass of at (, matching the region probed by the doubly lensed QSO), where is the cool gas covering fraction. The corresponding total mass in the cool phase is

where Z_cool is the cool gas metallicity. For and a solar-metallicity gas, which Z17 inferred for the cool phase, the inferred mass in the cool (T ∼ 10⁴ K) ISM is .

Assuming that the observed N(O vi) along sightline A is representative at , the estimated mass in the O vi-bearing phase of the ISM is

where is the covering fraction of the warm phase, Z_warm is its metallicity, and is the ionization fraction of ions. If we further assume a unity covering fraction and solar gas metallicity for the warm (T ∼ 10⁵ K) ISM phase, and adopt a reasonable , which is predicted for a wide range of physical conditions at T ∼ 10⁵ K (e.g., Oppenheimer & Schaye 2013), we find a total gas mass of in the warm ISM phase that is likely traced by O vi absorption. Despite the uncertainties inherent in our simple calculations, the estimated mass budgets in the cool and warm ISM phases are comparable to within a factor of a few if the two phases have similar gas covering fractions and metallicities.

In the physical picture where the observed O vi absorber traces transitional-temperature (T ∼ 10⁵ K) gas that is radiatively cooling from a virialized hot phase (T ∼ 10⁶ K), is proportional to the mass flow rate into the cool ISM following , where t_cool is the cooling timescale of the O vi-bearing gas. The cooling timescale depends on the gas temperature, metallicity, and density. For K and a solar-metallicity gas with a density of , typical of the hot halo of massive ellipticals at (e.g., Singh et al. 2018), the expected cooling time is (Gnat & Sternberg 2007; Oppenheimer & Schaye 2013) with a total cooling rate of . Thus, in a radiative cooling scenario, the estimated translates to a mass flow rate of at . If the bulk of this flow cooled to T ≲ 10⁴ K and remained in this phase, we should expect over a timescale of ∼100 Myr in the absence of significant star formation activity. Considering Z16 found a minimum stellar population age of and no detectable star formation () in the lens galaxy, this calculation suggests that most of the cooling gas is reheated to the coronal phase.⁷ To heat the gas back to virial temperature, the required heating rate is , assuming K, given the estimated mass of the dark-matter host halo of the lens (Z16).

Observations of nearby massive ellipticals show that mechanical feedback (often dubbed "radio-mode feedback") from active galactic nuclei (AGNs) can output as much power as (e.g., Werner et al. 2019). If the lens galaxy of HE 0047−1756 hosts an active nucleus at present, then in principle it has more power than what is required to reheat the O vi-traced cooling gas. There are two caveats to this statement, however. Because the estimated cooling time of a T ∼ 10⁵ K gas is short (∼10⁷ yr), owing to the expected high gas densities at , the actual amount of available heating depends sensitively on the radio-mode duty cycle (i.e., the fraction of time that an AGN is in radio mode). The radio-mode duty cycle in ellipticals has been estimated to be no more than outside of rich cluster environments (O'Sullivan et al. 2017). Furthermore, even if an AGN is currently on, its energy output is likely to be distributed over a large volume in the gaseous halo. Indeed, observations of large X-ray cavities/bubbles and extended radio lobes around nearby giant elliptical galaxies indicate that AGNs deposit kinetic energy on scales of ∼50 kpc or larger in the hot halo (e.g., McNamara & Nulsen 2007; Fabian 2012). In conclusion, whether AGNs are viable as a continuous heating source requires not only a high duty cycle but also that its mechanical energy can be effectively coupled with the ISM on ∼1 kpc scales in the galaxy.

Alternatively, we consider heating sources associated with the old stellar populations themselves. Previous analytic and simulation studies suggest that feedback from SNe Ia and stellar winds from asymptotic giant branch (AGB) stars may offset radiative cooling from diffuse gas in massive elliptical galaxies with (e.g., Conroy et al. 2015; Li et al. 2018). Empirically, the observed high abundance ratios at from quiescent galaxies (Z16; Zahedy et al. 2017a) also support the idea that their ISM has been subjected to significant influence from recent SNe Ia. Using the mean SN Ia rate in nearby ellipticals (e.g., Mannucci et al. 2005), Z17 estimated an integrated SN Ia rate of ∼0.3 per century within from the massive elliptical lens galaxy. Multiplying this rate by a mean energy of per SN Ia, we estimate that the heating rate available from SNe Ia is , which is comparable to the required heating.

In addition to SNe Ia heating, Conroy et al. (2015) also considered how materials ejected from AGB stars can interact with and heat the ambient ISM in elliptical galaxies. Using their analytic formula for the AGB heating rate, we estimate a heating rate of in the lens galaxy of HE 0047−1756. This exercise suggests that heating from SNe Ia and AGB stars may suffice to match the cooling rate inferred from the observed O vi absorption and prevent a large accumulation of cold gas in the ISM, even in the absence of strong feedback from an active nucleus.