Detection of the Diffuse H i Emission in the Circumgalactic Medium of NGC 891 and NGC 4565

We compare our GBT measurements with sensitive interferometric maps from the HALOGAS survey using WSRT. These interferometric data are very sensitive to any compact emission, but they lack our single-dish sensitivity to extended structure. Our GBT observation at any pointing includes both the high-column H i  clumps and the diffuse large-scale structures in the CGM, if any, as well as the stray light from the disk of the galaxy.

To carry out a rigorous comparison, we convolve the masked interferometric data cubes with the GBT beam. This moves the interferometric data to the angular resolution of the GBT. From the convolved interferometric maps, we can estimate the contribution of compact, relatively high-column structures to our observed GBT spectra. The convolved interferometer map will reflect any emission detected in the interferometer cube, i.e., any bright, compact CGM structure, extraplanar emission, and pickup of the disk of the galaxy from the extended wings of the GBT beam. The excess in the GBT spectra compared to the spectra from convolved WSRT cubes will correspond to the extended structures missed by the interferometer.

For this exercise, we perform the following steps:

• 1.  
  Noise estimation. In the interferometric maps, we do not see any emission in the region beyond r = 20' from the center of the target galaxy. Thus, we extract the pixel values beyond this radius in all velocity channels and use these to determine the noise. We estimate a single global noise value, σ_g, by fitting the negative half of the histogram of those pixel values with a Gaussian.
• 2.  
  Primary-beam correction. We use the MIRIAD task pbcor to correct for the primary-beam response of WSRT in the interferometric data cubes.¹⁷
• 3.  
  Masking. We are interested in emission securely detected by the interferometer. To identify such emission, we set a masking threshold and remove all emission in pixels with S/N below that. The interferometric cubes masked with S/N $\geqslant \,5$ are used to calculate the best estimate of the spectra. We vary the masking threshold from S/N $\geqslant \,3.5$ to S/N $\geqslant \,6.5$ to estimate the systematic uncertainty associated with this step.
• 4.  
  Circularization of the GBT beam model. The GBT beam is not a sharp tophat with width 91. It is not azimuthally symmetric either (e.g., Pingel et al. 2018). As the sky position of a given pointing changes during the course of an observation, the orientation of the GBT dish also changes. To account for that, we take the azimuthal average of the GBT beam response. The circularized GBT beam model along with its 1σ uncertainty from Pingel et al. (2018) is shown in Figure 3. Even separated by more than an FWHM, the GBT beam picks up emission from nearby bright sources, including the disk of the target galaxy.
• 5.  
  Convolution. For each pointing, we convolve the primary-beam corrected, masked WSRT cube with the circularized GBT beam model.¹⁸ We set the center of the circularized GBT beam response at the sky location of each GBT pointing (Table 1), and treat any missing or masked data in the WSRT cube as having zero intensity. The resulting spectrum represents what we would expect in the GBT spectrum if GBT observed only the securely detected emission in the WSRT cube. This includes stray light from the disk of the galaxy, extraplanar emission and any circumgalactic cloud with high enough column density to be detected by the interferometer (Figure 3, green points).
• 6.  
  Uncertainties. The statistical uncertainty in the convolved WSRT cube, σ_con is calculated from the global noise of the unmasked unconvolved WSRT cube, σ_g (see step 1) using the response distribution of the GBT beam across the sky: ${\sigma }_{\mathrm{con}}={\sigma }_{g}\times \tfrac{\sqrt{{\rm{\Sigma }}{w}^{2}}}{{\rm{\Sigma }}w}$, where w is the GBT beam response. The systematic uncertainty, as mentioned in step 3, comes by convolving the primary-beam corrected WSRT cubes masked with thresholds from S/N ≥ 3.5 to S/N ≥ 6.5.

We compare the shape and the integrated intensity of our GBT spectrum to the convolved WSRT spectrum at each pointing. We find three sets of results.

• 1.  
  Inner halo (${r}_{\perp }\leqslant 20\mbox{--}30$ kpc): Here the H i  emission is dominated by the disk in both the GBT and WSRT spectrum. The WSRT measurements are comparable with those from the GBT.
• 2.  
  Large impact parameters (pointings up/down 1.5, 2, and 3; ${r}_{\perp }\geqslant 40\mbox{--}50$ kpc): There is a significant excess in our GBT measurements compared to those from the WSRT (Figure 3). The WSRT columns account for 42 ± 20%(50 ± 11%) of the GBT columns of NGC 891 (NGC 4565).Note that the column densities at these pointings are below the sensitivity limit of WSRT ($\approx {10}^{19}\,{\mathrm{cm}}^{-2}$). Therefore, the masked WSRT cube does not have any high significance emission within ≈91 (the FWHM of the GBT beam) of these pointings. That means the masked and convolved WSRT spectra at these pointings represent the stray light from the H i  disk and the immediate extraplanar region. Therefore, we infer that the offset between our GBT and WSRT columns at these pointings represents the diffuse H i  emission from the CGM.This diffuse emission could consist of a mixture of small H i  clumps/clouds, and the diffuse, extended H i, as is evident from the spectral shape of GBT at these pointings (see Figures 1 and 2). The large line width of the spectrum implies that whatever the form, this emission must be widespread along the line of sight. Deeper interferometric data are necessary to distinguish between the small scale and extended H i  structures at these pointings.
• 3.  
  Intermediate impact parameters (pointings up/down 1; $20\mbox{--}30\leqslant {r}_{\perp }\leqslant 40\mbox{--}50$ kpc): Because the column density of H i  sharply drops off near pointings "up1/down1" in both galaxies, the convolved WSRT spectra at these pointings are extremely sensitive to the adopted GBT beam model. Here the assumption of a circular beam might not be adequate. Despite this uncertainty, the shape of the spectrum gives some clue as to the origin of the observed H i. The contribution of the H i  disk at any of our extended pointings shows a spectrum that resembles the disk spectrum. Mainly the amplitude, not the shape, depends on the details of the beam response. According to this criterion, at pointing "up1" around both galaxies and at pointing "down1" around NGC 891, we find that the emission cannot be solely explained by the disk contamination.In Figure 4, we plot the GBT spectrum at pointing "up1" around NGC 891, and compare it with the WSRT spectra convolved with two different GBT beam models, one with FWHM 86 and one with FWHM 96. Both WSRT spectra show 500 km s⁻¹ broad emission (shaded in gray) from the H i  disk. The GBT spectrum shows a similar component, with an amplitude between that from the 86 and 96 GBT beam models. This indicates that the true GBT beamwidth must be somewhere in between 86 and 96. The shape and the column densities of the GBT and the convolved WSRT spectra at the center of the galaxies are similar for the 96 beam, so the discrepancy is unlikely due to calibration offset between the single dish and the interferometer. Also, as the beam response is not azimuthally symmetric, circularizing the beammaps might affect the spectral shape and amplitude at these sensitive pointings.Our GBT spectrum shows a 200 km s⁻¹ wide emission (shaded in yellow) in addition to the 500 km s⁻¹ wide emission produced by the convolved WSRT. Because the interferometric data have limited sensitivity to any smooth, extended H i  component, we interpret this excess in the GBT compared to the WSRT as the emission from the diffuse, extended H i  structures.

Zoom In Zoom Out Reset image size

We compare our emission-based measurements with absorption measurements from the literature. From the Lyman series analysis of 44 COS-Halos systems, Tumlinson et al. (2013) set an upper limit of 10¹⁹ cm⁻² and a lower limit of 10¹⁶ cm⁻² on the H i  column densities. The large range reflected the difficulty of extracting column densities from Lyα absorption in this regime, when the damping wings are not yet strong but the lines are optically thick. Our column densities at $r\gtrsim 40\,\mathrm{kpc}$ are in the range of $5\times {10}^{17}\mbox{--}3\times {10}^{18}$ cm⁻², in good agreement with Tumlinson et al. (2013).

Based on more detailed modeling, Prochaska et al. (2017) estimated the N(H i) of 14 COS-Halo Lyman limit systems. The median and range of their estimated typical N(H i) profile are shown in Figure 5 (blue curve and hatched area). A direct comparison between our emission-based results and the absorption work by Prochaska et al. (2017) is complicated by the different targets and distinct nature of the observations. The N(H i) profile of galaxies is not necessarily independent of their galactic properties, and may not be azimuthally symmetric. The COS-Halo sample does not include our targets, and spans two orders of magnitude in stellar mass and specific star formation rate (sSFR). Tumlinson et al. (2013) and Prochaska et al. (2017) do not take the relative orientation of galaxy-quasar pairs into account in their analysis, so we cannot restrict our comparison to the minor axes of COS-Halos galaxies. All of this means that we can only compare average trends in the COS-Halo sample with our targets. Also, the absorption lines sample a "pencil beam" through the CGM, while the emission measurements convolve a large area, including potentially stray light pickup from the disk (see above). A direct comparison of absorption with our GBT measurement must take into account this different spatial sampling, which is not trivial.

Zoom In Zoom Out Reset image size

As a zeroth-order comparison, we subtract the convolved WSRT values from our GBT values at each pointing. When doing so, we propagate the statistical and systematic errors of each. If this excess does come from the CGM including the diffuse, extended H i  and small clouds/clumps, we might expect this to also be the column density sampled by absorption measurements. Unfortunately, imperfect knowledge of the beam model and azimuthal structure leave uncertainties in our subtraction of the WSRT. For the pointings "up1/down1," the calculation is more sensitive. Therefore, we make a rough estimate, considering only the 200 km s⁻¹ wide line sitting above the flat component as discussed in the previous section. By approximating the shape of the line as a triangle, we estimate the N(H i)_diff ∼ 5.0 (3.4) × 10¹⁸ cm⁻² for pointing "up1 (down1)" around NGC 891 and N(H i)_diff ∼ 1.4 $\times \,{10}^{18}$ cm⁻² for pointing "up1" around NGC 4565.

Figure 5 shows that the column density of the H i  that we detect in emission is generally larger than the median absorption seen by COS-Halos at a similar impact parameter. However, our measurements do lie well within the maxima of absorption. Therefore, our GBT measurements are roughly consistent with those from Lyman limit H i  absorption. By fitting the profile with an exponential: N(H i) = N(H i)_oexp(−r/r_s), we obtain the scale heights r_s to be 22 ± 7 kpc (NGC 891) and 36 ± 8 kpc (NGC 4565). This shows that the profile of diffuse H i  in the CGM is steeper in NGC 891 compared to that in NGC 4565. Also, we fit the galaxies together, and get an average r_s = 28 ± 6 kpc.

The H i  that we detect in the CGM could originate from the IGM, tidal interactions with other galaxies, dwarf satellite/companion galaxies, cold outflows, or condensation of cold gas from the hot CGM. Below, we consider each possibility.

Emission from companion galaxies: The H i  mass within the GBT beam is often a few times ${10}^{6}\mbox{--}{10}^{7}$ ${M}_{\odot }$ (Table 1). This is comparable to the H i  mass of a dwarf galaxy. However, the linewidths that we observe at those pointings (Figures 1 and 2, gray areas) are comparable with the velocity range of the target galaxies, which is also close to the velocity range expected for random motions in the halos of the galaxies. The spectrum of a dwarf companion (e.g., in Figure 2 the yellow area shows the spectrum of a dwarf companion to NGC 4565) would be noticeably narrower. We would also expect the deep interferometric observations to detect any reasonably compact emission like that expected for a dwarf galaxy. This disfavors the possibility of confusion with a satellite galaxy.

Tidal interactions: NGC 891 has a filament toward its companion UGC 1807 in the northwest (Oosterloo et al. 2007). NGC 4565 has two very nearby companions: NGC 4562 to the southwest, IC 3571 directly to the north (Zschaechner et al. 2012). The filament and the companion galaxies are captured in our pointings (Figure 1). If the diffuse H i  traces the interaction between the galaxies, the column density might be biased in the direction of the companions. However, the column density profiles are quite symmetric in two sides of the disks along the minor axes, lacking any noticeable excess toward the companions (Figure 3). Also, tidal features are often compact enough to be detectable in interferometer maps (e.g., the low column density (N(H i) $\approx {10}^{17}\,{\mathrm{cm}}^{-2}$) H i  filament between M31 and M33 (Braun & Thilker 2004; Lockman et al. 2012) was found to be clumpy and made up of small ≈kiloparsec-scale H i  clouds in high resolution observations by Wolfe et al. 2013). Given the morphology of the detected emission, symmetry of the profiles and the low likelihood of lining tidal features exactly along the minor axis in both sides of the disk, we consider it unlikely that these observations are driven by tidal interactions.

Inflow/outflow/ galactic fountain: We chose the minor axis to avoid contamination from the extended H i  disk. However, outflows and galactic fountains tend to be aligned with the minor axis, which could lead to preferential enhancement of column density along the minor axis. In principle, this can be checked by looking for velocity gradients or a correlation with the launching mechanism (SFR or AGN).

There is no appreciable difference between the column density profiles of the galaxies (Table 1, Figure 3). Given that NGC 891 appears to have a much higher sSFR (specific star formation rate) and ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ (surface density of SFR) than NGC 4565 (Table 1, Pingel et al. 2018), we might expect NGC 891 to launch more powerful outflow and/or show a more prominent signature of condensation. However, NGC 4565 hosts a LINER type AGN (Veron-Cetty & Veron 2000), which is a potential contributor to the outflow. Because the combined effect from star formation and AGN is quantitatively uncertain, the relation between the detected H i  and outflow/precipitation based on feedback argument remains inconclusive.

The metal content of the gas may help constrain its origin. From the UV absorption analysis of the warm halo along two sightlines 5 and 108 kpc above the disk near the minor axis of NGC 891, Qu et al. (2019) infer metallicities of 0.5 and 0.1Z_⊙, respectively. We cannot eliminate the possibility of a cold metal-enriched outflow from this. However, the difference in the metallicity at 5 and 108 kpc shows that the metal-enriched gas, even if ejected from the disk, has not reached a large distance.

Also, we try to infer the origin of the gas from kinematic information. In Figures 1 and 2, we have plotted the brightness temperature-weighted average velocity with respect to the systemic velocity of the galaxy along our line of sight. Because the galaxies are edge-on, the observed velocity would represent a small component of the true velocity for any motion perpendicular to the disk. If the dominant component of the emission is a quasi-static medium spread evenly across the halo, the average velocity should be consistent with the systemic velocity of the galaxy. If the observed emission is from the far side of the galaxy, positive (or negative) velocity would imply outflow (or accretion/precipitation). Similarly, positive (or negative) velocity of the emission from the near side of the galaxy would imply accretion/precipitation (or outflow). If a bipolar outflow or inflow is present, we should expect to see a velocity gradient with positive values in one side, and negative values in the other side of the disk.

We observe such gradient in the halo of NGC 4565 (Table 1). However, as we do not know the geometry and the position of the observed gas, we cannot determine whether it is an outflow or an inflow.

In the halo of NGC 891, the velocity is positive at all pointings and usually increases in magnitude with impact parameter (Table 1). Outflow or inflow alone cannot explain such a trend. Therefore, it could be produced by a combination of outflow, inflow, and galactic fountain activity. This is consistent with the observation and modeling of NGC 891 (Oosterloo et al. 2007; Fraternali & Binney 2008), where a galactic fountain is inferred from the increasing lag in rotation with height. It is also possible that the mean velocity of the disk of NGC 891 is slightly lopsided, so that it is not perfectly representative of the halo mean velocity. However, such warped H i  profiles are actually visible in both of our targets, with the stronger in NGC 4565 (Rupen 1991; Oosterloo et al. 2007; Zschaechner et al. 2012). Though not conclusive, it would be unexpected if the less prominent warp in NGC 891 affected the kinematic alignment significantly.

Condensation: Note that the absolute values of the average velocities are within 10% of the velocity width, i.e., $| {v}_{\mathrm{avg}}| \lt 0.1\times {v}_{\mathrm{width}}$. That is indicative of a quasi-static neutral medium permeating the halo in addition to any inflow, outflow, and galactic fountain. In the "hot mode" accretion common for galaxies as massive as our targets, the infalling gas from the IGM becomes shock-heated and virialized and forms a quasi-spherical hot halo (Rees & Ostriker 1977; Kereš et al. 2005). The diffuse H i  we detect could be a neutral phase mixed with this predominantly hot medium, which would eventually cool and accrete onto the galaxy.

Previous studies have explored extraplanar emission from both of our targets at other wavelengths. The studies find bright radio continuum, dust filaments, Hα recombination, and X-ray emission in the extraplanar region of NGC 891 but no corresponding features in NGC 4565 (Rand et al. 1992; Howk & Savage 1997, 1999; Rand 1998; Strickland et al. 2004; Tüllmann et al. 2006; Hodges-Kluck et al. 2018; Schmidt et al. 2019). These differences have been attributed to the stark difference in the star formation activity of these galaxies. This also suggests that the H i  halos of NGC 4565 and NGC 891 may have different origins.

All of these studies focused within 10 kpc of the disk. It is not clear if the effect of disk-halo interaction remains similar out to the large impact parameters as probed by our observations. The H i  profile is flatter in NGC 4565 (Section 3.2), with smaller H i  column within 40 kpc (Figure 5), as compared to NGC 891. This implies that despite having similar H i  masses within the volumes probed (Table 1), in the inner CGM NGC 4565 has less H i  mass than NGC 891 by a factor of a few. It is worth noting that the SFR of NGC 4565 is smaller than the SFR of NGC 891 by about the same factor as the H i  mass difference (Pingel et al. 2018). If our detected H i  is from the condensing hot CGM accreting to the galaxy as discussed in the previous paragraph, the mass difference of H i  in the inner CGM between the galaxies might imply a corresponding difference in H i  supply to the disk. This could be consistent with the "bathtub" model of star formation, in which the main factor controlling the SFR of a galaxy is the rate of accretion of H i  onto the galaxy disk (Dekel & Mandelker 2014).

We interpret the excess H i  observed by the GBT compared to the convolved WSRT as the H i  in the CGM (see Section 3.1). In this section, we discuss the contribution of this H i  to the total mass budget of neutral CGM, along the minor axes and to the whole CGM of our targets.

We define the overall halo mass fraction of the H i, f_h as follows:

$\begin{eqnarray}&&{f}_{h}=\displaystyle \frac{{\text{}}M{({\rm{H}}{\rm{I}})}_{\mathrm{CGM}}}{{\text{}}M{({\rm{H}}{\rm{I}})}_{\mathrm{disk}}+{\text{}}M{({\rm{H}}{\rm{I}})}_{\mathrm{CGM}}}.\end{eqnarray} \tag{ 1 }$

Here, M(H i)_CGM is the total H i  mass in the CGM, traced by our GBT observation. M(H i)_disk is the H i  mass of the disk. Therefore, M(H i)_disk + M(H i)_CGM is the total H i  mass. As we have probed out to the impact parameter of 90–120 kpc from the galaxies, we consider the CGM within this region only.

First, we consider the CGM along the minor axis, probed directly by our pointings. To calculate M(H i)_CGM, we subtract the convolved WSRT intensities from the GBT intensities. Then, we convert the diffuse column density estimates to the mass. For this exercise, we use only the pointings spaced by full beam size (i.e., up/down 1, 2, 3). Then we sum the mass from each of these pointings (see Table 1). For the pointings "up1/down1," we use the rough estimates obtained in Section 3.1. By comparing M(H i)_CGM with M(H i)_disk of NGC 891 (Oosterloo et al. 2007) and NGC 4565 (Zschaechner et al. 2012), we obtain f_h = 0.014 ± 0.001 ± 0.002 and 0.005 ± 0.001 ± 0.001 for NGC 891 and NGC 4565, respectively. This f_h is a lower limit, as this is measured only along the minor axis, not in the whole CGM. It accounts for about 1% of the total H i  mass, which illustrates why this emission has been so hard to detect.

We estimate the total mass of diffuse, extended H i  in the CGM assuming an azimuthally symmetric diffuse H i  distribution. At some impact parameters, measurements are available in both sides of the disk (e.g., up/down 1). We consider the average of the intensities at these pointings. Finally, we obtain f_h = $0.052\pm {0.003}_{-0.010}^{+0.008}$ and $0.020\pm 0.005\pm 0.006$ for NGC 891 and NGC 4565, respectively. As a sanity check, we compare these estimates to those implied by integrating the N(H i) profiles that we fit in Section 3.2. Integrating those profiles yields f_h = 0.047 ± 0.014 and 0.024 ± 0.007 for NGC 891 and NGC 4565, in excellent agreement with the values obtained directly from the measured column densities. Thus, adopting simple assumptions, the neutral CGM could, in total, account for roughly 2–5% of the total H i  mass.

From the area within 50 kpc around NGC 891 and NGC 4565, Pingel et al. (2018) estimated that f₁₉, the fraction of H i  mass below N(H i) = 10¹⁹ cm⁻² is 0.004–0.006 and 0.007–0.009, respectively. Our measurement is larger than the estimation of Pingel et al. (2018) by almost an order of magnitude. We suggest several potential causes for this discrepancy: (1) The simplest explanation would be that there is a large mass of diffuse, extended H i  that has column density $\lt {10}^{19}$ cm⁻² and that extends beyond 50 kpc. This emission might have been missed due to the lower sensitivity of the Pingel et al. (2018) observations. (2) We assumed an azimuthally symmetric diffuse H i  distribution. If our observed H i  is from a feedback-driven outflow or a large-scale galactic fountain, the mass could be preferentially concentrated along the minor axis. In that case, assumption of azimuthal symmetry would lead to an overestimated H i  mass. The lack of offset between the GBT and WSRT data in Pingel et al. (2018) might be due to the extended H i  disk contamination in the azimuthal averages. The diffuse extended H i  is likely abundant away from the major axis. As we focus only on the minor axis, we cannot confirm this in the whole CGM. Deep observations toward other points in the halo, including off-axis directions, will help clarify the situation. In any case, f_h differs from f₁₉ by definition, so it is not totally clear whether the results are in conflict or just defined differently.