Detection of diffuse HI emission in the circumgalactic medium of NGC891 and NGC4565 -- II

We probe the neutral circumgalactic medium (CGM) along the major axes of NGC891 and NGC4565 in 21-cm emission out to $\gtrsim 100$kpc using the Green Bank Telescope (GBT), extending our previous minor axes observations. We achieve an unprecedented $5\sigma$ sensitivity of $6.1\times 10^{16}$ cm$^{-2}$ per 20 km s$^{-1}$ velocity channel. We detect HI with diverse spectral shapes, velocity widths, and column densities. We compare our detections to the interferometric maps from the Westerbork Synthesis Radio Telescope (WSRT) obtained as part of the HALOGAS survey. At small impact parameters, $>31-43\%$ of the emission detected by the GBT cannot be explained by emission seen in the WSRT maps, and it increases to $>64-73\%$ at large impact parameters. This implies the presence of diffuse circumgalactic HI. The mass ratio between HI in the CGM and HI in the disk is an order of magnitude larger than previous estimates based on shallow GBT mapping. The diffuse HI along the major axes pointings is corotating with the HI disk. The velocity along the minor axes pointings is consistent with an inflow and/or fountain in NGC891 and an inflow/outflow in NGC4565. Including the circumgalactic HI, the depletion time and the accretion rate of NGC4565 are sufficient to sustain its star formation. In NGC891, most of the required accreting material is still missing.


INTRODUCTION
The circumgalactic medium (CGM), the multiphase gaseous extended halo surrounding the stellar disk and the interstellar medium (ISM), represents a large reservoir of baryonic material and plays a key role in the evolution of galaxies (Putman et al. 2012;Faucher-Giguere & Oh 2023).The inflow of H i, the raw fuel of star formation, from the intergalactic medium (IGM) to the ISM is regulated by the CGM through bimodal processes (Kereš et al. 2005;van de Voort et al. 2011;Nelson et al. 2013).In "hot mode" accretion, the accreting gas is shock-heated to the virial temperature, eventually, a fraction of it cools and fragments into H i clouds and precipitates to form the H i disk.It is predominant in massive (virial mass M 200 > 10 12 M ⊙ ) galaxies and galaxies at lower redshift ( ≲ 3).In lower mass galaxies or galaxies at higher redshift, the gas is predicted to accrete through "cold mode", i.e., it is never heated to the virial temperature and it flows through filaments/sheets from the IGM to the ISM bypassing any major thermodynamic interaction with the surrounding hot CGM.
The depletion time, i.e., the ratio of neutral gas (H i + H 2 ) mass ★ snskriti@stanford.edu in the ISM and the star formation rate (SFR) of galaxies is too small to sustain star formation over the cosmic time (Sancisi et al. 2008).
The stellar mass density surpasses the neutral gas density since  ≲ 2 (Walter et al. 2020).It implies that there has to be an external supply of gas to keep the ISM in a quasi-equilibrium state.The accretion/inflow rate measured from H i in and immediately around the disk of local (D<25 Mpc) galaxies is more than an order-of-magnitude smaller than their star formation rate (Kamphuis et al. 2022).It indicates that the rest of the "missing" accreting material could be at larger galactocentric distances, i.e., in the CGM.In addition to the accretion, circumgalactic H i also traces the tidal interaction, ram-pressure stripping, mergers between galaxies, and galactic outflows (e.g., Hibbard & van Gorkom 1996;Wolfe et al. 2013;Odekon et al. 2016;Veilleux et al. 2020).Ly absorption studies of H i in the CGM of the Milky Way have shown evidence for a large amount of atomic gas at N(H i) ≥ 10 18 cm −2 (e.g., Richter et al. 2017).The circumgalactic H i in external galaxies has been studied using Ly absorption lines at  ⪅ 0.2−0.85(e.g., Borthakur et al. 2016;Péroux et al. 2022).These studies have revealed a strong correlation between the atomic gas mass fraction (M HI /M ★ ) and the impact-parameter-corrected equivalent width of Ly absorbers, suggesting a physical connection between the atomic ISM and the CGM.However, the Ly absorption lines are often found saturated but undamped, i.e., 10 16 <N(H i)<10 18.5 cm −2 , a ≈3 orders of magnitude uncertainty.Absorption studies are limited by the availability of bright background sources like quasars and the covering fraction of the Ly absorbers.The average of many galaxyquasar pairs is challenging to interpret as the stellar mass, SFR, and inclination of the galaxy disks, and azimuthal angles of CGM absorbers to the disks are often entangled with each other.Also, the sensitivity to Ly absorbers at  ≈ 0 with existing instruments is limited, and it is unknown how the accretion history from  ≈ 0.2 − 0.85 to  ≈ 0 evolves, if at all.21-cm H i emission studies can constrain the large-scale spatial distribution and kinematics of neutral gas in the CGM of an individual  ≈ 0 galaxy in an unbiased and complete way.
Over the past decade, there have been several large interferometric programs to map the H i emission in the CGM of nearby spiral galaxies at high angular resolution and high point source sensitivity, including Hydrogen Accretion in LOcal GAlaxies Survey (HALOGAS, Heald et al. 2011), The Local Volume H i Survey (LVHIS, Koribalski et al. 2018), MeerKAT HI Observations of Nearby Galactic Objects; Observing Southern Emitters (MHONGOOSE, Sorgho et al. 2019) survey, and Widefield ASKAP L-band Legacy All-sky Blind surveY (WALLABY, Kleiner et al. 2019).However, due to the finite size of the dishes, there is a minimum possible spacing (i.e., the diameter of each dish) between neighboring telescopes of an interferometric array.This results in a gap in the uv-coverage, which is called the "short-spacing" problem.It limits the sensitivity of the interferometers to low surface brightness and large angular scales.For example, the smallest baseline length of the Westerbork Synthesis Radio Telescope (WSRT) is 36 m, which translates to a maximum recoverable angular scale of 20 ′ .Because extended H i structures are more likely to survive than small H i clouds in the CGM (Armillotta et al. 2017;Gronke et al. 2022), the "diffuse" circumgalactic H i with low column density across large angular scales might be a significant reservoir of gas around galaxies that interferometers would miss.Single dishes have full uv-coverage, allowing detection of structure at all angular scales.Thus, single-dish observations complement interferometric studies of the CGM.
The unblocked aperture design of the Green Bank Telescope (GBT), decent angular resolution (FWHM 9.1 ′ ), its low sidelobes, and high surface brightness sensitivity (T sys ⩽ 20 K) make it ideal to search for low column density structures over large scales.There have been several single-dish observations with the GBT targeting the CGM of nearby galaxies, as discussed below.
The CGM of NGC 2403, NGC 2997 and NGC 6946 have been mapped with the GBT down to a 5 detection limit of ∼10 18 cm −2 over a 20 km s −1 channel (de Blok et al. 2014;Pisano 2014).Observations of NGC 2403 revealed a low column density, 8 kpc extended filament outside the H i disk.It has recently been found to be part of a 20 kpc stellar stream that formed due to a recent (2 Gyr) tidal interaction with a dwarf satellite galaxy (Veronese et al. 2023).The halo of NGC 2997 did not show any filamentary feature, and the H i mass, as measured with the GBT, was only 7% higher than that derived from past interferometric measurements.The H i observations of NGC 6946 revealed a filament of N(H i) = 5 × 10 18 cm −2 connecting the galaxy with its nearest companions.It shows that the H i distribution in the CGM varies from galaxy to galaxy.
As part of the project AMIGA (Absorption Maps In the Gas of Andromeda), Howk et al. (2017) looked for the 21-cm emission with the GBT in the CGM of M 31 along 48 directions at a 5 sensitivity of 4×10 17 cm −2 .Except for the 2 sightlines passing through the Magellanic Stream, they were not able to detect any H i emission across an impact parameter of 25 kpc to 340 kpc.The non-detection indicated that unless M 31 is an exception, the covering fraction of the circumgalactic H i might be too low to be detected in the galaxies as nearby as M 31, where the beam of GBT corresponds to a physical length of 2 kpc.Pingel et al. (2018), hereafter P18, observed the CGM of four galaxies at a distance of 9-18 Mpc: NGC 891, NGC 925, NGC 4414 and NGC 4565 with the GBT, and compared it with the interferometric data from the HALOGAS survey at equal spatial and velocity resolution.They did not detect any considerable amount of excess H i from the GBT data than the WSRT data.It indicated that the column density of diffuse H i, if present, is lower than their 5 GBT detection limit of 0.9-1.4×10 18cm −2 over a 20 km s −1 line width.
The aforementioned single-dish surveys integrate for a few minutes per point and cover a large area, usually out to the virial radius of each galaxy.Das et al. (2020) adopted a complementary approach and integrated for 3-4 hours on-source toward individual pointings along the minor axes of NGC 891 and NGC 4565, achieving an unprecedented 5 sensitivity of 8.2×10 16 cm −2 over a 20 km s −1 line width.It is more than an order of magnitude deeper observation compared to previous single-dish observations.It led to detections of H i at each pointing.By comparing the sensitive GBT spectra to the deep HALOGAS WSRT measurements, they found that ≈ 50% of the emission detected by the GBT cannot be explained by the interferometric maps.This confirmed the presence of low column density diffuse circumgalactic H i along the minor axes of these galaxies.
In the GBT census of 18 local MHONGOOSE galaxies, the azimuthally averaged cumulative H i profiles of 15 galaxies do not truncate at large impact parameters (Sardone et al. 2021).This indirectly suggests the presence of diffuse H i at a large scale.With a 5 H i column density sensitivity of 10 17.9 cm −2 , Wang et al. (2023) mapped the CGM of NGC 4631 using Five-hundred-meter Aperture Spherical Telescope (FAST) and detected H i out to 120 kpc from the galaxy.By comparing it with HALOGAS WSRT measurements, they found that more than one-fourth of the emission detected by the FAST is diffuse.These results are consistent with the findings of Das et al. (2020), hereafter D20.
In this paper, we expand the study of D20 along the major axes of NGC 891 and NGC 4565.We achieve a 5 H i column density sensitivity of 6.1×10 16 cm −2 , equivalent with H i mass sensitivity of 2.1 × 10 5 M ⊙ , calculated over a 20 km s −1 line width.These are among the deepest H i observations around external galaxies obtained to date in 21-cm H i.
We outline our observation and the data reduction in section 2, and the analyses in section 3. We interpret our results in section 4. We summarize our conclusions and comment on future directions in section 5.

OBSERVATIONS AND DATA REDUCTION
We observed the CGM of NGC 4565 and NGC 891 in 21-cm (1.42 GHz) emission as part of GBT projects 20B-360 (PI: Das) and 21B-324 (PI: Das).These observations add a total of eight new, deep spectra around NGC 891 and eight new, deep spectra around NGC 4565 (Figure 1a, triangles and squares).We combine these with observations from GBT project 15B-257 (Figure 1a, pentagons) published in D20.Including all observations, there are three pointings on each side of each target, with pointings stepping along the major axes (20B-360) and minor axes (15B-257, 21B-324), separated from one another by the FWHM beam size.
To minimize the contamination of the observed spectrum by emis-

Figure 1a
. The 21-cm emission spectra for the observed GBT pointings in the CGM of NGC 891, smoothed and re-binned at v (see Table 1).The green contours on the Digitized Sky Survey (DSS) image show the integrated 21-cm intensity from the HALOGAS survey, with the lowest contour showing an H i column density of 10 19 cm −2 .The GBT pointings observed in 15B-257, 20B-360, and 21B-324 are illustrated with pentagons, triangles, and squares on the DSS image, with red/blue symbols implying the mean line-of-sight velocity of the detected emission being positive/negative (arrows directed from pointings toward respective spectra).The direction of disk rotation has also been marked with '+/−' symbols at the edges of the H i disk.The GBT beam (hatched circle) and the scale are shown at the top and the bottom of the DSS image, respectively.The companion galaxy UGC 1807 is visible to the northwest of NGC 891.In each panel of the spectrum, the maximum rotational velocity of the H i disk is marked with vertical dashed gray lines, with the systemic velocity of NGC 891, 528 km s −1 , subtracted from the actual line-of-sight velocity.The horizontal gray patch denotes the ±RMS value of T B .The horizontal gray line has been drawn at T B = 0 to guide the eye.The vertical hatched region shown in some of the panels is the velocity region of the Milky Way that has been masked throughout the data reduction and analysis.We calculate the integrated intensity at each pointing by summing T B over the velocity range shaded in yellow.At the pointings far from the disk, we detect emission signature at a level of N( sion from the disk, we choose the major axes pointings at least 1 FWHM beam size away from the edge of the H i disk measured in interferometric surveys (Figure 1a, contours), and minor axes pointings at least 1 FWHM beam-size away from the center of the galaxy.
The details of the observations, including the sky direction of the pointings, the effective integration time at each pointing, t eff , and the system temperature, T sys are provided in Table 1.
We observed by position switching, using the L-band receiver with the Versatile GBT Astronomical Spectrometer (VEGAS; bandwidth = 23.5 MHz) as the backend.The off-source pointings are chosen to have the same declination as the target pointings and are placed ≈ 1.5 • away from the disk.This places the off-source pointing beyond the virial radii of the target galaxies.We also make sure that there is no known H i source within 1 FWHM beam size of the off-source pointings in the velocity range of concern.We used the quasars 3C 48 (in 15B-257, 20B-360) and 3C 123 (in 21B-324) as the primary flux calibrator for NGC 891 and 3C 286 for NGC 4565.The H i disk of each galaxy was observed before the observations of the CGM to verify the setup.
We use an improved version of the routine developed in GBTIDL by D20 to reduce the position-switched data.We obtain the aperture efficiency and the atmospheric opacity at the zenith, which are used to derive T cal , based on the elevation during each integration and at each frequency.We use that to construct an off-source spectrum from each off-source scan, then subtract this from the on-source spectrum corresponding to the adjacent on-source scan.The calibrated, offsubtracted spectrum of each scan in a session is then weighted by the combined integration time and system temperature of the on-off pair, t eff /T 2 sys , to calculate a combined spectrum.As part of the processing, we also identify channels visibly affected by RFI (radio frequency interference) and replace these with values linearly interpolated from nearby RFI-free parts of the spectrum.
We construct a spectrum at each pointing in each session for each polarization in units of brightness temperature, T B .We shift the velocity axis of each spectrum so that v = 0 km s −1 corresponds to the systemic velocity of the galaxy, v sys .This is 528 km s −1 for NGC 891 and 1230 km s −1 for NGC 4565 (Rupen 1991).We accumulate the spectra from individual sessions at the native velocity resolution of 1.25 km s −1 .
For each pointing, we combine the spectra for each polarization from all sessions, again weighted by the total t eff /T 2 sys for each singlesession spectrum.Thus, we obtain the mean spectrum for each po-larization, in units of brightness temperature.Finally, we average the two polarizations, which increases the final S/N by a factor of √ 2. After constructing the final spectrum, we apply an iterative signal finding and baseline fitting technique, which we describe in Appendix A. The algorithm tests a variety of channel widths and baseline orders.It identifies regions of candidate emission, fits regions away from candidate emission, and subtracts the lowest order robustlyfitted polynomial baseline.Table 1 reports the optimum channel width , velocity range used for the baseline fit Δv base , the velocity range of identified emission v emit , and the order of the polynomial fit  fit .Because this differs somewhat from the approach used by Das et al. (2020), we also re-reduce, re-stack, and then fit new baselines to the spectra from that paper (project 15B-257).This gives us a single homogeneously reduced data set.
As a final check, we compare the RMS noise measured from signal-free regions of the baseline-subtracted spectrum with the noise expected based on the radiometer equation using t eff and T sys (Table 1).The measured RMS noise is on average a factor of 1.18 +0.08 −0.12 and 1.36 +0.60  −0.10 larger than the expected noise for NGC 891 and NGC 4565.This could be due to residual imperfections in the baseline subtraction that require higher-order polynomial or more sophisticated baseline () "UP" and "DN" denote offset to higher and lower declination than the center of the target galaxy.Pointings along the major and minor axes of the target galaxy end with "J" and "N", respectively.Pointing "nJ" is n×9.1 ′ away from the edge of the H i disk measured in interferometric surveys (see Figure 1a), and pointing "nN" is n×9.1 ′ away from the center of the galaxy, where 9.1 ′ is the FWHM size of the GBT beam.The transverse distances of the pointings from the center of the target galaxy are provided in the 7th column, adopting a distance of 9.2 Mpc/10.8Mpc for NGC 891/NGC 4565 from Heald et al. (2012).
() Effective exposure time t eff = t on ×t off /(t on +t off ).Where t on and t off are the integration times toward the on-source and the off-source pointings, respectively.To maximize t eff , we maintain t on and t off to be the same.
() The expected noise at v = 5 km s −1 based on the radiometer equation.
() The velocity resolution at which the spectrum is smoothed.
() The velocity range over which the H i emission is detected.
(  ) The velocity range that we use to estimate the baseline.
() The order of the best-fitted polynomial to describe the baseline.
fitting techniques.Below we use the measured RMS noise rather than the radiometer equation-based estimates, so we expect our S/N estimation to be conservative and capture this additional uncertainty.

ANALYSIS
In the final spectra, we detect H i emission in most of the pointings (Figure 1a and 1b, yellow shaded area).We integrate each spectrum over the range of Δv emit to obtain the intensity, ∫ T B dv in units of K km s −1 .We calculate the statistical uncertainty in the integrated intensity using  T B , v and Δv emit .We add the systematic uncertainty in flux calibration (see §2) to the statistical uncertainty in quadrature.Except for the pointings immediately next to the H i disk, our detections at all other pointings lie below the 5 sensitivity limit of the GBT-HALOGAS survey (P18).
We convert the integrated intensity at each pointing to column density, N(H i) GBT , assuming that the emission is optically thin and the number density of the gas is well above the critical density of H i for collisions with electrons (7×10 −6 cm −3 at 10 4 K; Draine 2011).We do not apply any correction for the inclination angle in computing the column density, because we do not know the geometry and morphology of the detected emission.We show the derived N(H i) as a function of galactocentric radius in Figure 2.
We calculate the intensity-weighted mean velocity over the range of Δv emit : vGBT . We also calculate the H i mass within each beam from the integrated intensities, adopting a gain of 1.86 K Jy −1 .We quote the N(H i) GBT , the associated H i masses, and vGBT in Table 2.

Comparison with interferometric data
Our GBT spectra include emission both from high column density H i clouds and any diffuse, extended CGM structures.They may also include stray light from the disk of the galaxy, reflecting that even the relatively clean GBT beam can still pick up emission from the bright inner region of the galaxy.To help disentangle these components, we .The pointings at higher/lower declination than the galaxy disk are plotted at positive/negative angular separation from the galaxy center, respectively.The error bars along the y-axis include systematic uncertainties and statistical uncertainties.The error bars along the x-axis denote the GBT beam size.In our GBT data, systematic uncertainties come from the difference in T cal across different sessions.In the WSRT data, systematic uncertainties come from the difference in masking threshold (S/N>4-5-6), and circularization of the GBT beam.The shaded region around the beam model is 1 uncertainty in the beam response due to averaging the beam map to one dimension.The positive offset between our GBT data and the convolved WSRT data indicates the presence of diffuse circumgalactic H i. compare our GBT measurements to sensitive interferometric maps from the HALOGAS survey using WSRT (Heald et al. 2011).These interferometric data are sensitive to compact, high column density structures and they capture the main disk of the galaxy well, but they lack sensitivity to extended, low surface brightness structures.
To carry out a rigorous comparison, we convolve the masked interferometric data cubes with the GBT beam so that the interferometric data are at the angular resolution of the GBT.From the convolved interferometric maps, we can estimate the contribution of compact, relatively high column density structures to our observed GBT spectra, i.e., any extraplanar H i emission, and pickup of the galaxy disk from the sidelobes of the GBT beam.The excess in the GBT spectra compared to the spectra from convolved WSRT data cubes will correspond to the extended low column density structures missed by the interferometer.
Same as D20, we follow the method of 1) global noise estimation and 2) primary-beam correction of the low-resolution1 WSRT data cubes, 3) masking of the primary-beam corrected WSRT data cubes, 4) circularization of the GBT beam model, 5) convolution of the masked WSRT data cube with the circularized GBT beam, and 6) estimation of the statistical uncertainty in the convolved WSRT spectra.The WSRT data cubes masked with S/N⩾ 5 are used to calculate the best estimate of the spectra.The additive systematic uncertainty comes from the convolution of the WSRT data cubes that are masked with thresholds of S/N⩾ 4 and S/N⩾ 6.In Appendix A we discuss why S/N< 4 is not a good mask.The circularized beam map accounts for the continuous change in the orientation of the GBT beam throughout the observation, but it averages out the details of the shape of the GBT beam.To account for the azimuthal asymmetry in the GBT beam, we incorporate a multiplicative systematic uncertainty   () The statistical uncertainty is followed by the multiplicative systematic uncertainty due to the variation in flux calibration from session to session.
() The column density obtained from the masked WSRT data convolved with the GBT beam map, estimated over the velocity range where emission is detected in the GBT.Here, the statistical uncertainty is followed by the systematic uncertainty.The latter includes multiplicative systematic uncertainty due to averaging the azimuthally asymmetric GBT beam map to 1-D and additive systematic uncertainty in the masking threshold, added in quadrature.
() The excess H i detected by the GBT compared to that by the WSRT, i.e., the diffuse N(H i) in the CGM. that depends on the sky position of the pointing and thus is different for every pointing.
By integrating the masked and convolved WSRT spectrum over Δv emit at each pointing, we obtain the H i column density from WSRT, N(H i) WSRT .We compare N(H i) GBT to N(H i) WSRT in Figure 2. We find two sets of results:

1) Outer halo (pointings UP/DN 2 and 3):
There is significant excess in our GBT measurements compared to the WSRT estimates across most of the velocity channels of Δv emit (Figure C.1).The excess emission, N(H i) GBT -N(H i) WSRT , is 75 +9 −11 % and 81 +5 −8 % of the GBT emission, N(H i) GBT , of NGC 891 and NGC 4565, respectively.
Because the column densities at these pointings are way below the sensitivity limit of the WSRT data cubes (≈ 10 19 cm −2 ), the masked WSRT data cube is unlikely to have any high S/N emission within the main lobe 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 extraplanar region.
Therefore, the offset between the GBT and the WSRT estimates at these pointings is due to the diffuse H i emission from the CGM.From the diverse shape, velocity width, and amplitude of the emissions, we cannot determine whether it consists of extended H i structure and/or discrete H i clouds/clumps.It would require high-sensitivity interferometric measurements to distinguish between these cases.

2) Inner halo (pointings UP/DN 1):
The emission detected in GBT at some of the velocity channels is explained by the WSRT.Thus the excess in GBT than WSRT is narrower than Δv emit (Figure C.1).The excess emission accounts for 54±11(42 +10 −11 )% of the GBT emission of NGC 891 (NGC 4565).In addition to the stray light from the H i disk of the target galaxy and its companions (e.g., NGC 4562), the masked and convolved WSRT spectra at these pointings represent the small-scale H i in the extraplanar region.Thus, the offset between our GBT and WSRT columns at these pointings represents the diffuse and large-scale H i emission from the CGM.
We test it more rigorously by redefining noise and masking through intermediate convolution in Appendix D and confirm this inference.

Properties of the CGM
We obtain the H i column density in the CGM, N(H i) CGM , by subtracting the H i column density of the masked and convolved WSRT spectrum from that of our GBT spectrum at each pointing: The subscripts 1 and 2 correspond to the GBT and the WSRT, respectively.The column density of the H i disk, N(H i) disk , is not the same across all observing sessions of the GBT due to systematic errors.Also, they are not the same as the N(H i) disk estimated from the masked and convolved WSRT spectrum.Therefore, we apply their ratio as a correcting factor in the latter term in equation 1.
We fit the radial profile of N(H i) CGM as an exponentially decreasing profile: The two parameters N(H i) o and   are the central column density of H i and the scale radius.We also fit different subsets of the data to test the azimuthal symmetry and axial symmetry around the H i disk.
Using the best-fitted N(H i) CGM profile from equation 2, we calculate the H i mass in the CGM by integrating the surface mass density,   N(HI) CGM (), over the area of concern: To calculate the mass within the main lobe of the GBT beam, we use the area with Δ being the FWHM of the GBT beam.To calculate the mass in each quadrant, i.e., within ±/4 of the major and minor axes, we use the area We obtain the average line-of-sight velocity of the H i emission in the CGM, vCGM , using equation 4.
Here, brightness temperature is integrated over the velocity range of Δv emit , and subscript 1 and 2 corresponds to the GBT and the WSRT, respectively.
We calculate the accretion rate using equation 5a.The assumed location of the circumgalactic H i is illustrated in Figure 6.We assume the velocity in the sky plane to be toward the galaxy disk (i.e., radially inward) and the same as our measured line-of-sight velocity.
Here the radial profile of N(HI) CGM comes from equation 2. For the major axes pointings, we estimate the accretion rate also using the angular velocity, () in the disk plane: where () = vCGM ()/ (5b)

RESULTS AND DISCUSSION
Below we discuss our results on the density, mass, velocity, and accretion rate of the detected H i in the CGM of NGC 891 and NGC 4565 and interpret them in the context of star-forming activities in these galaxies.

Density profile
We note the N(H i) CGM at each pointing in Table 2.We show the bestfitted radial profile of N(H i) in the CGM of NGC 891 and NGC 4565 in Figure 3.
In the CGM of NGC 891 we have excluded two pointings, UP 3J and DN 3J, from the fit (highlighted with red circles in Figure 3, left panel).These pointings are the farthest from the center of NGC 891 along the major axes.Because the velocity of the H i emission at UP 3J pointing is significantly different from its adjacent pointings (see Figure 1a), this emission might not be co-spatial with its adjacent pointings.Despite moving at similar velocities, the N(H i) toward DN 3J is larger than that toward DN 2J, which is not suitable for an exponentially decreasing radial profile.The excess N(H i) toward DN 3J than DN 2J might be a result of unresolved H i cloud(s) embedded in diffuse H i. Or, the radial profile of the diffuse H i column density at this radius might become flatter than an exponential.Verifying this scenario would require further deep and high-spatialresolution observation around DN 2J-DN 3J pointings.
Overall, NGC 891 has a 1.6 ± 0.5× steeper profile with 4.7 +13.3 −3.0 × higher central H i density,  o =N(H i) o /2  , than NGC 4565.It results in 1.6 +9.1 −1.0 times more H i in the extraplanar region (i.e., at and around UP/DN 1 pointings) of NGC 891 than NGC 4565.This ratio is similar to the ratio of SFR of these galaxies.This supports the "bathtub" model of star formation where the main factor controlling the SFR of a galaxy is the available H i in and around the galaxy disk (Dekel & Mandelker 2014).Our result is qualitatively consistent with the finding of D20 based on the four pointing (blue triangles in Figure 1a and 1b) along the minor axes of these galaxies.
In Figure 4, we show the best-fits (panel b) and the best-fit parameters (panel a) for the different subsets of our data: pointings only along the major axes ('J') or minor axes ('N'), only at higher declination ('UP') or lower declination ('DN') than the H i disk, etc.In each galaxy, the scatter of  o across ≈2 orders-of-magnitude and scatter of   across a factor ≈ 3 for these different subsets indicate that the circumgalactic H i is not symmetric in two sides of the disk and along the azimuth.

NGC 891
H i profiles are different along the major and minor axes, with the major axis having a 10.2 +1.8 −4.4 × denser and 1.5 +0.5 −0.4 × steeper H i profile (Figure 4(b), top left).It could be due to the active accretion of H i from the CGM to the disk along the major axis.
The H i profiles along the minor axis are different on two sides of the disk (Figure 4(b), top middle), with the higher declination side of the disk ('UPN') having a 6.0 +3.6 −2.3 × denser and 1.6 ± 0.3× steeper H i profile.There is already a galactic fountain along the minor axis (Oosterloo et al. 2007;Fraternali & Binney 2008).Thus the H i asymmetry in two sides of the disk along the minor axis might be due to a diffuse extended part of the fountain which is not mapped in shallow WSRT observations.We also fit the major and minor axes pointings at higher declination, i.e., UP J and UP N together, and compare it with lower declination pointings (Figure 4  See Table 1 for the terminology of the pointings.The markers with darker shades and uncapped error bars correspond to the global best fits (shown in Figure 3).the data into two diagonal sides of the disk, northwest and southeast, respectively.The H i profile is 9.1 +11.7 −5.5 × denser and 1.6 +0.5 −0.4 × steeper at the higher declination ('UP') pointings.There is a H i filament (Oosterloo et al. 2007) toward UGC 1807 in the northwest of NGC 891 (green contours in the DSS image of Figure 1a).Thus the H i asymmetry in two diagonal sides of the disk might be the result of a diffuse extended arm of that filament which is not captured in shallow interferometric maps.

4.1.1.1
What is the origin of H i at UP3J?Given the shape and width of the H i spectrum and the mass of the H i emission at the UP 3J pointing of NGC 891 (top left panel of Figure 1a), we cannot rule out the possibility of a dwarf galaxy.However, we did not find any known object within 10 ′ of this pointing at the similar line-of-sight velocity of the detected H i emission in NED2 .Its velocity offset from the systemic velocity of NGC 891 is similar to the highvelocity clouds found in the CGM of Milky Way (Putman et al. 2012;Richter et al. 2017).Whether or not the H i emission detected at this pointing is a high velocity H i cloud in the CGM of NGC 891 and/or a yet undetected dwarf galaxy needs further investigation with deep interferometric and optical surveys, respectively.

NGC 4565
Along the major and minor axes, the H i profiles have similar slopes (Figure 4(b), bottom left).The best-fit   is higher by a factor of 4.8 along the major axis possibly due to active accretion, but it is consistent with   along the minor axis due to the large errors in density estimate along the major axis (Figure 4(a)).
The lower declination side of the disk along the major axis ('DNJ') has a 7.6 +6.8 −4.6 × denser and 2.0 +0.7 −0.5 × steeper H i profile than the higher declination side (Figure 4(b), bottom right).The asymmetry might have caused/be caused by the strong warp around the disk of NGC 4565 observed in 21-cm (Zschaechner et al. 2012).Also, NGC 4565 hosts ∼5-6 dwarf galaxies (Carlsten et al. 2020) along the higher declination arm ('UPJ').The flatness in the H i profile on this side might be due to ram-pressure/tidally stripped H i from these dwarf satellite galaxies.
Along the minor axis, the H i profiles have similar slopes on both sides of the disk, but   is 2.7 +4.4  −1.6 × larger on the lower declination side of the disk ('DNN'), resulting in more H i mass in the CGM (Figure 4(b), bottom middle).This enhancement could be due to tidal interactions between the H i disks of NGC 4565 and its companion NGC 4562 to the southwest (contours in the DSS image of Figure 1b).
Thus, the detected circumgalactic H i can be primarily explained by a diffuse extended H i, with some H i clouds embedded in it, and the interaction of the target galaxies with their companions.
Thus, both galaxies have more mass along the major axes than the minor axes, indicating active accretion along the major axes.The masses quoted in D20 based on minor axes pointings and assuming azimuthal symmetry were underestimated.
M(H i) CGM is 11.4 ± 4.2% and 5.4 ± 0.8% of the M(H i) disk of NGC 891 and NGC 4565, respectively.Based on the GBT and WSRT measurements within 50 kpc of NGC 891 and NGC 4565, P18 estimated  19 , the fraction of H i mass below N(H i) = 10 19 cm −2 to be 0.4-0.6% and 0.7-0.9%,respectively.Our measurement is larger than the estimation of P18 by an order of magnitude.It suggests the presence of a large amount of diffuse, extended H i with N(H i) below 10 19 cm −2 , extending beyond 50 kpc.This emission might have been missed due to the lower sensitivity of the GBT observations presented in P18.Also,  19 in P18 was based on azimuthal averages in circular annuli around the galaxies, whereas we focus on pointings along principal axes.Therefore, most of the circumgalactic H i could be concentrated along principal axes of the galaxies, leading to a higher estimate in our case.
The depletion timescale,  dep = (M(H i) disk + M(H i) CGM )/SFR is 2.1 ± 0.8 Gyr for NGC 891 and 11.4 ± 1.8 Gyr for NGC 4565.By combining this information with the amount of H i in the extraplanar region estimated in the previous subsection, we infer the following.NGC 891 might have accreted most of its star-forming fuel by now and will soon run out of fuel unless there is a huge H i reservoir in the regions not probed in our observations -regions away from the principal axes, and regions at the larger galactocentric radius.NGC 4565, on the other hand, is forming stars almost quiescently; it has an ample supply of H i to continue forming stars at this rate.

Velocity
We show the line-of-sight mean velocity of the circumgalactic H i in Figure 5 and quote them in Table 2. Except for the UP 3J pointing of NGC 891, the velocities at all other pointings of both galaxies display interesting trends that we discuss below.

Major axis
In both galaxies, the H i velocity in the CGM along the major axes is aligned with the disk rotation (illustrated in panels (a) and (b) of Figure 6).Such co-rotation has been observed previously in a few Mg ii absorbers at  ≈ 0.2 (Ho et al. 2017;Martin et al. 2019) and at  ≈ 1 (Zabl et al. 2019) and in one Ly absorber at  ≈ 0.6057 (Weng et al. 2023), but not in H i emission at  ≈ 0. It could be the neutral CGM that is accreting onto the disk keeping the angular momentum conserved.
Otherwise, the H i disk might not truncate sharply and could be larger than what is found in shallow (≈ 10 19 cm −2 ) interferometric observations (the green contours on DSS images in Figure 1a and 1b).In that case, our detected H i would be the outer part of an extended H i disk.The column density of this co-rotating H i is consistent with the predicted "proto-disk" extended out to ≈80 kpc, modulo the degeneracy between the volume-filling factor, covering fraction and intensity of the cosmic ionizing background (Bland-Hawthorn et al. 2017).
Distinguishing between the scenarios of a co-rotating neutral CGM or an extended H i disk is beyond the scope of this paper.It would require deep co-spatial observations of H emission and H i emission at high angular resolution in the future.

Minor axis
The velocities of circumgalactic H i along the minor axes show a similar gradient in both galaxies (Figure 5), although the velocities are smaller in NGC 891 than NGC 4565.
Here the velocity profile is complicated by the possibility of inflow, fountain, or outflow.The interpretation of a pure inflow or a pure outflow is degenerate to the location of the H i emitter in the CGM along our line-of-sight.In an inflow, the redshifted/blueshifted emission is on the nearer/farther side of the disk than the observer (panel (d) of Figure 6).The fountain, in principle, has both outflowing and inflowing components.If the fountain is not perfectly aligned along the minor axis, the mean line-of-sight velocity would depend on whether it is titled toward/away from the observer.If the redshifted/blueshifted emission is on the farther/nearer side of the disk, it would be an outflow or fountain if its velocity in the sky plane is away from or toward the disk, respectively.
In NGC 891, the presence of a galactic fountain is already known from the increasing lag in rotation with height (Oosterloo et al. 2007;Fraternali & Binney 2008, also discussed in §4.1).Our observed velocity gradient could be due to a cosmic inflow and/or a diffuse extended part of the fountain that is tilted away from us.In the latter scenario, the smaller velocities compared to the maximum rotational velocity of the disk could be due to H i emitters moving at opposite velocities in the fountain suppressing the mean velocity (panel (c) of Figure 6).a weak outflow-dominated phase to an inflow-dominated phase in the past 40 Myr.Thus we cannot conclusively comment on the origin of the H i emitters along the minor axis of NGC 4565 without constrain-ing their velocity in the sky plane.While inflows and outflows should have distinct metallicity, measuring it might not be useful in reality as the CGM is predicted to be inhomogeneously mixed in simulations and there is no observational evidence for a correlation between absorber metallicity and azimuthal angle (Martin et al. 2019;Weng et al. 2023).

Accretion rate
In NGC 891, the total accretion rate, M CGM , at the position of all observed pointings integrated within the GBT beam is 0.30 ± 0.03 M ⊙ yr −1 .For all pointings except UP3J and DN3J (as they were excluded from the fitting of N(H i) CGM in §4.1), we extend the accretion rate calculation to each quadrant in the CGM (i.e., ±/4 region around the major and minor axes) assuming azimuthal symmetry, and obtain M CGM of 0.73 ± 0.08 M ⊙ yr −1 .Adding it to the accretion rate measured within the GBT beam of the excluded pointings, the total M CGM is 0.80 ± 0.08 M ⊙ yr −1 .M CGM along the major and minor axes are 0.65 ± 0.08 M ⊙ yr −1 and 0.15 ± 0.01 M ⊙ yr −1 .
In NGC 4565, the total M CGM at the position of all observed pointings integrated within the GBT beam is 0.20 ± 0.02 M ⊙ yr −1 .Extending the accretion rate calculation to each quadrant in the halo, we obtain a total M CGM of 0.83±0.08M ⊙ yr −1 .M CGM along the major and minor axes are 0.66 ± 0.07 M ⊙ yr −1 and 0.17 ± 0.03 M ⊙ yr −1 .
Using the second definition of accretion rate (equation 5b), M CGM along the major axes of NGC 891 and NGC 4565 are 0.29 ± 0.04 M ⊙ yr −1 and 0.57 ± 0.06 M ⊙ yr −1 , respectively.
Thus, the total accretion rate in both galaxies is similar, and it is larger along the major axes (irrespective of the definition of accretion rate used) than the minor axes, indicating active accretion onto the disk.The specific star formation rate, sSFR = SFR/M ★ , and surface density of SFR, Σ SFR are an order-of-magnitude larger in NGC 891 than NGC 4565 (see Table 2).Despite this difference, the similarity in H i accretion rate suggests that the intermediate stages between the accumulation of H i and star-formation, e.g., the H i -H 2 conversion, are inefficient in NGC 4565 compared to NGC 891.

Tracing accretion: found or missing?
In NGC 4565, M CGM /SFR is 1.24 ± 0.12.Because we cannot distinguish between an outflow or an inflow along the minor axis of NGC 4565, excluding the estimates along the minor axis, we obtain M CGM /SFR of 0.98 ± 0.10 (or 0.85 ± 0.09 using the second definition of accretion rate in equation 5b).This is arguably the first galaxy where the accretion rate is consistent with the SFR.It implies that the CGM of NGC 4565 is continuously supplying H i to the disk for a sustainable star formation.Excluding NGC 4562 and IC 3571 that are already detected in WSRT maps, NGC 4565 hosts 19 dwarf galaxies down to   ∼ −10 within 150 kpc (Carlsten et al. 2020), with H i masses of each being < 10 6.88 M ⊙ (Karunakaran et al. 2022).Determining if the accreting H i is dominated by the cosmic inflow or the ram-pressure/tidally stripped H i from these satellite galaxies is beyond the scope of this paper.If H i along the minor axis pointings trace an outflow, the mass-loading factor, M/SFR, is 0.25 ± 0.04.
In NGC 891, M CGM /SFR is 0.36 ± 0.04 (or 0.13 ± 0.02 using the second definition of accretion rate in equation 5b).Thus, ≈60-90% of the required accreting material is still missing.Below we discuss some of the possibilities of where the missing accretion could be: (i) Beam smearing: If the length scale of the H i in the sky plane is smaller than the GBT beam, its intensity could be underestimated due to the GBT beam smearing that dilutes the intensity of any probed structure.The smaller the size of the H i emitters, the higher the dilution.The WSRT can capture up to the length scale of 20 ′ which is larger than the GBT beam.Thus the beam smearing is not expected to affect the excess H i detected by the GBT than the WSRT.However, our GBT observation is more than 2 orders-of-magnitude deeper than the WSRT maps.Thus, some of the excess H i detected by the GBT could be smaller than the GBT beam whose column density is below the sensitivity of the WSRT maps.While we can rule it out for pointings close to the disk (Appendix D), it remains plausible for the far-away pointings.(ii) Opacity: If the length scale of the H i along the line of sight is smaller than that assumed here, the number density and hence the opacity of the measured H i would be higher.Thus, for a given intensity, the column density would increase after correcting for the self-shielding.This is more likely to be the case for the major axis pointings if the corotating H i we find in §4.3 is genuinely the outermost disk instead of the CGM.(iii) Sky-plane motion: We have calculated the accretion rate based on the line-of-sight motion.If the circumgalactic H i is spiraling in with the velocity in the sky plane being larger than the line-of-sight velocity, the accretion rate could be larger.(iv) Ionization: Given the mass of NGC 891, it is likely a "hot mode" accretion where the H i forms out of the hot CGM.Near the galaxy disk, H i is pressure confined due to the larger density and temperature of the hot CGM.At the interface of the hot ionized CGM and the neutral CGM, partially ionized warm/cool CGM can form because of thermal conduction and enhanced cooling rate due to condensation (Li et al. 2020).If our detected H i is partially ionized, the mass of H i containing gas would be higher than the current estimate.The hydrogen is rapidly ionized at ≈ 10 4 K, thus a 60-90% ionization is not impossible.(v) Off-axis: If most of the accretion is happening away from the principal axes, the accretion rate might have been underestimated.(vi) Larger distance: There might be more H i, diffuse and/or com-pact, at larger impact parameters.The H i is less likely to be destructed at large galactocentric distances due to the lower density of the surrounding hot CGM (Li et al. 2020).Whether or not it immediately contributes to the SFR would depend on its density and velocity.
Testing these possibilities would require highly sensitive mapping of the whole CGM out to the virial radii of these galaxies using single-dish and interferometers in the future.

In-situ fueling
Most of the discussion has assumed that the circumgalactic H i contributes to the star formation in the disk.However, if there is in-situ star formation in the CGM of these galaxies, the circumgalactic H i might be locally supplying fuel to those star formation hubs.
NGC 891 has abundant stellar structures in its CGM, e.g., a thick star-forming cocoon extended out to 40 kpc/15 kpc along the major/minor axis, a giant ∼90 kpc stellar stream, a few arching stellar streams up to ∼50 kpc looping around the stellar disk and a stellar halo extended out to ∼50 kpc (Mouhcine et al. 2010;Monachesi et al. 2016).NGC 4565 has a ∼60 kpc stellar halo, but it is more massive, flatter, and is a larger fraction of the stellar disk mass in NGC 891 than in NGC 4565 (Harmsen et al. 2017).The stellar metallicity, [Fe/H], steeply declines radially in the stellar halo, implying a comparatively young stellar population in the CGM (Monachesi et al. 2016).NGC 4565, showing no evidence of stellar streams/cocoons, has its stellar disk sharply truncated at 3 kpc/26 kpc along minor/major axis (Martínez-Lombilla et al. 2019).The inner CGM of NGC 891 has an extended dust component likely embedded in an atomic gas (Howk & Savage 1997;Bocchio et al. 2016;Yoon et al. 2021) while NGC 4565 does not have any (Howk & Savage 1999).Also, NGC 891 has only 7 (3 confirmed) dwarf satellite galaxies within 200 kpc (Carlsten et al. 2022), while NGC 4565 has 21 satellite galaxies within 150 kpc (Carlsten et al. 2020).Whether some of the stellar structures in the CGM of NGC 891 are disrupted satellites, and whether NGC 891 might have accumulated H i through minor mergers with its satellites in the past are not currently constrained.
Thus, the H i in the CGM of NGC 891 and NGC 4565 might have different origins.The circumgalactic H i might be locally fueling star formation in the CGM of NGC 891.In that scenario, the line-ofsight velocities of H i would follow that of the stellar structures in the CGM.NGC 4565, on the other hand, might be accreting through tidally/ram-pressure stripped H i from its satellites.This might explain the different velocities of circumgalactic H i in these galaxies.

CONCLUSIONS
In this paper, we have extended the study of Das et al. (2020) and presented the GBT observations of NGC 891 and NGC 4565 in 21-cm emission along their major axes out to ∼ 100 − 120 kpc in the CGM, and also revisited the observations in D20.Instead of the traditional mapping which is more expensive due to the required exposure time, we stare at a few selected (24) pointings for ≈ 3 − 4 hours and achieve an unprecedented 5 sensitivity of 6.1 × 10 16 cm −2 for 20 km s −1 velocity width at 5 km s −1 velocity resolution.
We substantially revised our GBT data reduction routine and the baseline subtraction method.We follow the methods outlined in P18; D20 to compare our GBT observation and interferometric WSRT maps from the HALOGAS survey.We include the systematic uncertainty to account for the azimuthal asymmetry in the GBT beam and the variation in position angle throughout our GBT observation.We also rigorously test the angular size of the detected diffuse H i by redefining noise and masking through intermediate convolution of the WSRT data cubes.Below we list our scientific results: (i) We detect H i emission at all GBT pointings at > 4.3 significance including systematic uncertainties, with 3 pointings below 5.
(ii) More than 31 − 43% and 64 − 73% of the H i emission detected by the GBT at small and large impact parameters cannot be explained by the WSRT maps.It implies the presence of a diffuse neutral CGM.The shape, velocity width, and column density of the diffuse H i are quite diverse, indicating a complex morphology and kinematics.(iii) The CGM of NGC 891 has a steeper profile than NGC 4565, with the ratio of H i mass in the extraplanar region being 2 +9 −1 .It is consistent with their SFR ratio, supporting the "bathtub" model of star formation.(iv) The N(H i) profiles have different normalizations and scale radii for different subsets of the data, indicating a lack of azimuthal symmetry and axial symmetry around the H i disk.These asymmetries can be attributed to the diffuse extended arm of the filament and fountain of NGC 891, ram-pressure stripped H i from the satellite galaxies of NGC 4565, and tidal interaction with the companion galaxies of NGC 891 and NGC 4565.(v) The mass ratio of H i in the CGM and the disk is 0.11 ± 0.04 and 0.05 ± 0.01 for NGC 891 and NGC 4565, an order of magnitude larger than previous estimates based on shallow (≈ 10 18 cm −2 ) GBT mapping within 50 kpc of these galaxies.(vi) The velocity of the diffuse H i along the major axes pointings are aligned with the disk rotation in both galaxies, suggesting a corotating CGM or the outer part of an extended H i disk.Along the minor axes pointings, the velocities resemble the signatures of an inflow and/or a galactic fountain in NGC 891.In NGC 4565 we cannot distinguish between the scenarios of an inflow or an outflow.(vii) The depletion timescale and accretion rate in NGC 4565 are > 9.6 Gyr and > 0.51 M ⊙ yr −1 , large enough to sustain its star formation.NGC 891, with 60 − 90% of the required accreting material still missing, would run out of its raw fuel by 2.1 ± 0.8 Gyr.
The Parkes-IMAGINE (Imaging Galaxies Intergalactic and Nearby Environments) survey has mapped the entire CGM of 28 nearby galaxies in H i, reaching sensitivities comparable to our GBT observation (Sardone et al. 2020).In the future, coordinated and deeper interferometric surveys and single-dish observations outfitted with phased array feeds or FLAG (The Focal L-band Array for the GBT; Pingel et al. 2021) architecture would substantially improve our understanding of the large-scale H i accretion in the CGM, the true radial extent of H i disk, and the relation between H i accretion and star-formation in the disk.Measurements of low ionized metal absorbers would constrain the gas-phase metallicity and the ionization fraction of the H i emitters.Along with kinematic information, it would help disentangle the origin of circumgalactic H i unambiguously and better understand the role of H i in galaxy evolution.positive fluctuation is considered a candidate emission in a v = 5 km s −1 spectrum.(v) We return to the smoothed spectrum in step # 2, and fit a new baseline to the spectrum.This time we exclude the Δv candidate region from fitting if the candidate emission is not at rest (i.e., not centered around 0 km s −1 ).If the candidate emission is at rest, the algorithm can identify it anyway.For example, we do not exclude any velocity region for all minor axis pointings of NGC 891 (except UP3N, where the 50-150 km s −1 region is excluded).For the major axis pointings of NGC 891, we exclude the 0-250 km s −1 region at lower declination ('DN'), and the -250 to 0 km s −1 region at higher declination ('UP', except 'UP3J', where we exclude the 100-400 km s −1 region).
If needed, we adjust the velocity range used for the baseline fit Δv base to ensure that it is larger than Δv candidate and includes enough spectral range for a good fit.This best-fit baseline is the "intermediate baseline".
(vi) We calculate the velocity range over which T B in the "intermediate baseline"-subtracted spectrum is continuously positive around the center of the Δv candidate region.That velocity range is set to be the new Δv candidate .Because we look for H i emission bound to the halo, the maximum allowed candidate velocity range is set by Δv disk (except 'UP3J' of NGC 891).
(vii) With the new exclusion region defined in step # 6, we repeat step # 5, this time using the "intermediate baseline"-subtracted spectrum as input.We repeat steps #5 and #6 until Δv candidate converge over successive iterations.Then the converged Δv candidate is taken as the velocity range of detected H i emission, Δv emit , for the working v.
The best-fit baseline of the converged fit is the "final baseline".
(viii) We calculate the RMS noise,  T B from the signal-free part of the "final baseline"-subtracted spectrum.Then, we integrate the spectrum over the range of Δv emit to obtain the intensity, I = ∫ T B dv in units of K km s −1 .We calculate the statistical uncertainty in the integrated intensity,  I =  T B √︁ vΔv emit .(ix) We repeat steps # 2 to 8 for different values of v.(x) We choose the optimum velocity resolution, the final v, to be the value where the H i emission is detected at maximum signal-tonoise ratio (S/N = I/ I ) in the smoothed and final-baseline-subtracted spectrum.
In Table 1, we provide the optimum velocity resolution, the final velocity range of H i emission Δv emit , the final velocity range of baseline Δv base , and the best-fit polynomial orders  fit that we find for each pointing.

APPENDIX B: DECIDING THRESHOLD FOR MASKING WSRT DATA CUBES
(i) Noise Calculation: To estimate the noise in the WSRT data cubes, we consider two "emission-free" regions by excluding data within a radius of 20 ′ and 30 ′ from the center of the target galaxy.A 20 ′ mask is appropriate because signals can be visually identified within this radius, although some of our pointings lie beyond it.A 30 ′ mask provides a truer representation of the "emission-free" region, although it leaves a smaller amount of data to estimate the noise.Thus there is a trade-off between the accuracy and amount of the data used to estimate the noise.Because the positive values in the "emission-free" region might be affected by the faint (and visually unidentifiable) signals, we fit the histogram of only negative values in the "emission-free" region with a Gaussian.We consider the bestfitted standard deviation of that Gaussian as the global noise of the data cube.In the following steps, We consider the minimum of the noise values calculated using the two "emission-free" regions.(ii) Masking: Spectra from the unmasked and convolved WSRT cube would be heavily contaminated by noise, especially at the pointings more than 1 GBT beamsize away from the H i disk.It might lead to the WSRT column densities being significantly larger than the GBT column densities, which is physically impossible.That makes it useless to compare unmasked and convolved WSRT spectra with the GBT spectra.To mask the WSRT data cubes based on noise calculated in the previous step, we consider the S/N of each pixel in each velocity channel.We retain the value of a pixel at a given velocity channel if the S/N of that pixel at that velocity channel and in at least one of the adjacent velocity channels is above the required threshold.This allows us to avoid spurious emissions at a deserted velocity channel.We consider thresholds of 1-7 at a step of 1. Masked pixels in the cube are rendered with a value of zero.(iii) Convolution and Spectra Extraction: We convolve the WSRT cubes masked at different thresholds with the circularized GBT beam and extract the spectrum at each pointing.In Fig. B.1, we show the convolved WSRT spectrum toward the UP 1N pointing of NGC 891 from the WSRT cubed masked at seven different masking thresholds.The spectra corresponding to the cube masked at S/N ⩾ 1 and 2 are discarded since the baseline is far from zero and is systematically positive by construction.In the S/N ⩾ 3 masking, the baseline estimate is significantly improved.But it still retains some residual noise in the cube (see the mask in Therefore, we discard this threshold as well.Additionally, S/N ⩾ 7 does not offer any significant improvement compared to S/N ⩾ 6 and is therefore discarded.Thus, we discard thresholds of S/N⩾ 1-3 and ⩾ 7, while retaining S/N ⩾ 4 − 6.We consider the spectrum from the cube masked at S/N⩾5 as the best estimate.We use the spectra from S/N ⩾ 4 and 6 masked cubes to calculate the additive systematic uncertainty.The GBT and WSRT spectra for our observed GBT pointings of NGC 891.The vertical yellow patch is the region where the H i emission is detected in the GBT data.The maximum rotational velocity of the H i disk is marked with vertical dashed gray lines.The green and red curves correspond to WSRT cubes masked at S/N>5 and S/N>3, respectively, and convolved with the circularized GBT beam.The vertical extent of the green shaded regions around the green curves is due to masking the WSRT data cube at S/N> 4 (upper extent) or S/N> 6 (lower extent), and it includes 1) statistical uncertainty in WSRT data, and 2) systematic uncertainty due to the azimuthal asymmetry of the GBT beam.The masking at S/N> 3 (or any lower threshold) is discarded because it picks up a lot of spurious noise at UP/DN 2 and 3 pointings, and does not pick up any true excess than the S/N> 4 masking at UP/DN 1 pointings.

Figure 1b .
Figure 1b.The 21-cm emission spectra for the observed GBT pointings in the CGM of NGC 4565.The systemic velocity of NGC 4565, 1230 km s −1 , has been subtracted from the actual line-of-sight velocity in each spectrum.The companion galaxy NGC 4562 is visible to the southwest of NGC 4565.The vertical hatched region shown in some panels of the spectrum is the velocity region of NGC 4562 that has been masked throughout the data reduction and analysis.See the caption of Fig.1afor more details.

Figure 2 .
Figure2.The radial N(H i) profile of NGC 891 (top) and NGC 4565 (bottom) along their major axes (left) and minor axes (right).The pointings at higher/lower declination than the galaxy disk are plotted at positive/negative angular separation from the galaxy center, respectively.The error bars along the y-axis include systematic uncertainties and statistical uncertainties.The error bars along the x-axis denote the GBT beam size.In our GBT data, systematic uncertainties come from the difference in T cal across different sessions.In the WSRT data, systematic uncertainties come from the difference in masking threshold (S/N>4-5-6), and circularization of the GBT beam.The shaded region around the beam model is 1 uncertainty in the beam response due to averaging the beam map to one dimension.The positive offset between our GBT data and the convolved WSRT data indicates the presence of diffuse circumgalactic H i.

Figure 3 .Figure 4 .
Figure 3. N(H i) in the CGM of NGC 891 (left) and NGC 4565 (right), with best-fit models.The filled and unfilled symbols correspond to pointings at higher and lower declinations than the galaxy disk.The physical separation corresponding to the angular separation from the disks is noted in the top x-axis, adopting a distance of 9.2 Mpc and 10.8 Mpc for NGC 891 and NGC 4565 from us, respectively(Heald et al. 2012).In the left panel, the points marked within red circles have been excluded from the fit.NGC 891 has a steeper profile than NGC 4565, with more H i closer to the disk of NGC 891.
Figure 4. (a)The H i density at the center,   =N(H i)  /2  vs the scale radius,   , for the models best-fitted to different subsets of the data, labeled accordingly.See Table1for the terminology of the pointings.The markers with darker shades and uncapped error bars correspond to the global best fits (shown in Figure3).(b) N(H i) in the CGM of NGC 891 (top 3 panels) and NGC 4565 (bottom 3 panels), with best-fit models for different subsets of the data, labeled in the legend of each panel.

Figure 5 .Figure 6 .
Figure 5.The line-of-sight mean velocity of H i in the CGM of NGC 891 (left) and NGC 4565 (right) in the rest frame of the galaxy.The horizontal dotted lines denote the maximum rotational velocity of the H i disk measured in the GBT data.The gray circles show the rotation curve of the disk along the major axis, measured in the WSRT data.The spatial extent of the disk along the minor axis is shown with the gray vertical patch.The horizontal dashed line is drawn at 0 km s −1 to guide the eye.
Figure A.1.Illustration of the baseline fitting and subtraction routine.Top 3 rows: 'DN1N' pointing of NGC 891 for v = 7.5 km s −1 , 4th and 5th rows: 'DN2N' pointing of NGC 891 for v = 10 km s −1 , bottom 2 rows: 'UP2N' pointing of NGC 4565 for v = 5 km s −1 .Top 2 rows, 4th and 6th row: Each iteration, noted in parentheses, consists of two panels.The left panels show the original spectrum with the candidate emission region indicated and the best-fit baselines of different polynomial orders plotted as colored lines.The right panels show the best-fit baseline-subtracted spectrum, with the best-fit polynomial order written in the bottom left corner.The vertical gray dotted lines denote the line width of the disk, Δv disk , for the target galaxy, NGC 891.The value of Δv candidate (shown with vertical yellow bands) at the start and end of each iteration are noted as well.The fitting algorithm at this pointing converges in the fourth iteration.3rd, 5th, and bottom row (from left to right): Δv candidate , N(H i), S/N, and number of velocity channels (=Δv candidate / v) of the candidate emission for different v.The gray horizontal band denotes Δv disk of NGC 891 in the leftmost panel.The filled and unfilled circles mark the upper and lower boundary of the candidate emission, respectively.In the rightmost panel, the horizontal line is drawn at 5, the threshold to identify a positive fluctuation as a candidate emission.The vertical dashed line in each panel corresponds to the optimum v.
Figure B.1.The WSRT spectrum toward the UP 1N pointing of NGC 891 convolved with the GBT beam.The legend details the S/N thresholds to mask the WSRT data cube before convolving.The opacity of the red shade is inversely proportional to the S/N threshold.The yellow vertical patch depicts the velocity range of the H i disk.
Fig. D.1, top row, left panel), leading to the residual emission in the convolved spectrum.While it is visible in the baseline region in contrast to the expected constant zero value of the baseline, it can affect any velocity channel, thus leading to an overestimated column density of the WSRT spectrum.It becomes more problematic for pointings farther out from the disk (Fig. C.1).
Figure C.1.The GBT and WSRT spectra for our observed GBT pointings of NGC 891.The vertical yellow patch is the region where the H i emission is detected in the GBT data.The maximum rotational velocity of the H i disk is marked with vertical dashed gray lines.The green and red curves correspond to WSRT cubes masked at S/N>5 and S/N>3, respectively, and convolved with the circularized GBT beam.The vertical extent of the green shaded regions around the green curves is due to masking the WSRT data cube at S/N> 4 (upper extent) or S/N> 6 (lower extent), and it includes 1) statistical uncertainty in WSRT data, and 2) systematic uncertainty due to the azimuthal asymmetry of the GBT beam.The masking at S/N> 3 (or any lower threshold) is discarded because it picks up a lot of spurious noise at UP/DN 2 and 3 pointings, and does not pick up any true excess than the S/N> 4 masking at UP/DN 1 pointings.

Figure C. 2 .
Figure C.2.The difference in brightness temperature weighted velocity (top) and brightness temperature (bottom) between our GBT data and convolved WSRT data of the UP 1N pointing of NGC 891.The hatched region showswhere the GBT has excess detection than WSRT.The area under the curve within the hatched region in the bottom panel is the excess integrated intensity of H i detected by the GBT.The ratio of areas under the curves within the hatched regions in the top and bottom panels is the mean line-of-sight velocity of that excess H i.
Figure D.1.Mask cubes summed over all velocity channels of the WSRT data cube of NGC 891, for S/N>3 (left), S/N>4 (middle), and S/N>5 (right).Black regions with a value of 1 are extracted from the data cube.Top row: Mask for the original data cube, where pixels above the required S/N in two adjacent velocity channels are retained in the masked data cube.In addition to the H i disk of NGC 891, the filament, and UGC 1807, there are a few pixels in the CGM with S/N>4.However, S/N>3 masking allows a lot of pixels across the field which is likely spurious.This was one of the reasons for choosing S/N> 4 − 5 − 6 masking in Appendix A instead of S/N> 3 − 4 − 5 as shown here.Bottom three rows: Mask based on the noise estimated from the data cube that is convolved with 1 ′ (second row), 3 ′ (third row), and 5 ′ (fourth row) Gaussian beams.At every velocity channel, pixels above the required S/N and all pixels within 1 ′ , 3 ′ , or 5 ′ of those pixels are retained in the masked data cube, respectively.

Table 1 .
Details of observations and data reduction

Table 2 .
Details of measurements