Signatures of Gas Flows-I: Connecting the kinematics of the Hi circumgalactic medium to galaxy rotation

In order to connect the CGM absorption to the host galaxies, we need the galaxy morphologies and the alignment of the quasar sightline relative to the galaxy disks (i.e., inclination and azimuthal angles). We obtained previously published inclination angles and azimuthal angles for 23 galaxies (French & Wakker, 2020), which are listed in Table 1. We further obtained previously published morphologies/geometries for 28 galaxies (Kacprzak et al., 2015; Kacprzak et al., 2019a; Kacprzak et al., 2019b), which were computed using GIM2D (Simard et al., 2002) models of HST images with ACS, WFC3, or WFPC2 in the F702W, F814W, or F625W filters as listed in Table 2. Here we add new models for 5 galaxies having HST images and 5 galaxies with Pan-STARRS (hereafter PS; Chambers et al., 2016) images, which are listed in Table 2. The orientation of the galaxies, such as their inclination ($i$) and azimuthal angle ($\Phi$), were modelled following the methods adopted from Kacprzak et al. (2011a) and Kacprzak et al. (2015). We fit two-component disc+bulge models using GIM2D (Simard et al., 2002) to the HST and PS images using modelled point spread functions (see Kacprzak et al., 2015). The galaxy disk component is modelled with an exponential profile and the bulge component has a Sersic profile with $0.2<n<4.0$. The modelled inclination and azimuthal angles, and their errors, for all galaxies are listed in Table 1. We adopt the convention of the azimuthal angle ${{\Phi}}=0^{\circ}$ to be along the galaxy major axis and ${{\Phi}}=90^{\circ}$ to be along the galaxy minor axis.

Our galaxies have a full range of azimuthal and inclination angles (see Table 1). We test for potential effects of a biased distribution in azimuthal and inclination angle on our results by conducting a one-dimensional Kolmogorov-Smirnov (KS) test on $\Phi$ and $i$ distributions. Our analysis indicates that the azimuthal angle is consistent with the expected flat distribution for a random sample of galaxies at a significance level of $1.36\sigma$. We also find that the inclination angles are consistent with the expected sin($i$) distribution for a random sample of galaxies (Law et al., 2009), at the significance level of $1.84\sigma$. Although a preference towards edge-on galaxies is ideal for examining outflows and inflows, we have determined that our results remain unchanged within the errors reported here when we exclude galaxies with $i<30$ deg (14 galaxies in total). Thus, we include all galaxies in our analysis.

The behaviour of the CGM is dependent on galaxy mass and its properties vary with location within the virial radius of the halo (Chen et al., 2010; Churchill et al., 2013a, b; Tumlinson et al., 2013; Oppenheimer et al., 2016; Ng et al., 2019). Therefore we have computed the stellar masses for all galaxies using the rest-frame $g-r$ colour and $g$-band mass-to-light ratio ($M/L$) relation from Bell et al. (2003). Galactic extinction corrections are not applied as they are on the order of the uncertainties in our method. Given the range of redshifts and galaxy angular extents in our sample, we have used a range of catalogues, such as PS, SDSS (York et al., 2000), and DESI Legacy surveys (Dey et al., 2019), to obtain the galaxy photometry and colours.

For galaxies with $z>0.05$, we obtained $g$ and $r$ Kron magnitudes from PS. In cases where PS photometry was not available, we used either SDSS model magnitudes or HST photometry. To apply the $M/L$ ratio relation, the rest-frame colours of galaxies are required. We applied $K$-corrections following the methods described by Nielsen et al. (2013) to obtain rest-frame absolute $g$- and $r$-band magnitudes. For the galaxies with no $g-r$ observed colour, we assumed an Sbc type, which represents the typical galaxy colour found for galaxies associated with CGM absorption. (Steidel et al., 1994; Zibetti et al., 2007; Nielsen et al., 2013; Kacprzak et al., 2015).

Galaxies with $z<0.05$ have larger angular extents, which typically results in underestimated $g$-band magnitudes from SDSS and PS. To obtain their $g$-band absolute magnitudes we adopted the colour transformation from Blanton & Roweis (2007), ${g=B-0.2354+0.3915((g-r)-0.6102)}$, to convert $B$-band magnitudes to $g$-band magnitudes. The $B$-band magnitudes are computed using the $B$-band galaxy luminosity function of Marzke et al. (1994) with the galaxies’ luminosities adopted from French & Wakker (2020). The galaxy colours are measured using DESI Legacy Survey imaging in $g$ and $r$ bands. Here, $K$-corrections are not applied for these low redshift galaxies since they are negligible. The measured $g-r$ colour for all the galaxies is presented in Table 1.

Using our uniformly computed photometry, we show the g-band absolute magnitude distribution of our galaxies as a function of redshift in Fig. 1(c). We find that the vast majority of our galaxies reside above $M_{g}=-18$ with a few less luminous galaxies. Our sample appears to be mass complete near $M_{g}=-19$, however this cannot be concluded since this plot is limited to absorbers only and does not account for the Hi absorption–mass dependence (Bordoloi et al., 2018). Since this work attempts to eliminate objects with major companions that could have a larger influence on the CGM kinematics, we are less concerned about lower mass companions like LMCs, dwarfs, etc., which can be considered as part of the more massive halo and their kinematics. So having a complete survey to the same depths in each field (i.e., similar central-to-satellite mass ratio limit) is more important than equal mass sensitivity (i.e., $10^{9}$ M_⊙) across all fields, which permits us to examine lower masses at lower redshifts. We have also verified that our results remain unchanged within the errors reported here when we exclude low luminosity/mass galaxies ($M_{g}>-19$, 9 galaxies in total). Thus, we include all galaxies in our analysis.

To test the validity of our computed masses, we compared our values with the stellar masses of the 11 galaxies that overlap with the COS-Halos sample (Werk et al., 2013). We found a mean difference of 0.065 dex between the two samples., which provides confidence in our mass estimates.

We also converted the galaxy stellar masses ($M_{\ast}$) to halo masses ($M_{\rm h}$) using the stellar-to-halo mass relation (SHMR) from Girelli et al. (2020). We adopted the parameterised SHMR in two redshift bins of $0.0\leq z<0.2$ and $0.2\leq z<0.5$. The best-fit parameters with a relative scatter of 0.2 dex from their Table 2 are used for our conversions and uncertainty calculations. The distribution of galaxy masses is shown in Fig. 1. The halo masses have a median value of $\log(M_{\rm h}/M_{\odot})=11.5$ and span a full range of $10.5<\log(M_{\rm h}/M_{\odot})<12.7$. The virial radius of all the galaxies in our sample is calculated following the formalism of Bryan & Norman (1998). The virial radii span a range of $61<R_{\rm vir}<324$ kpc with a median value of $R_{\rm vir}=131.5$ kpc. The virial radius normalised impact parameters have a range of $0.1\leq D/R_{\rm vir}\leq 5.0$, with a median value of $D/R_{\rm vir}=0.91$. The galaxy masses, $R_{\rm vir}$ and $D/R_{\rm vir}$ can be found in Table 1.

To compare the CGM kinematics to the kinematics of the galaxies, we require galaxy redshift zeropoints and their rotation curves. The galaxy kinematics for 23 galaxies were obtained from French & Wakker (2020). We further obtained spectra for 39 galaxies using the Keck/ESI over the course of 10 observing nights across 2010, 2014, 2015, and 2016. The wavelength coverage of ESI is 4000 – 11000 Å, which covers a range of emission lines like the [Oii] doublet, $\rm{H}\beta$, the [Oiii] doublet, $\rm{H}\alpha$, and Nii doublet. The width of the ESI slit was set to $1^{\prime\prime}$ and it is $20^{\prime\prime}$ long. The slit position angle was selected to be aligned with the optical major axis of each galaxy to acquire their full range of rotation velocities (see Fig. 2). The echellete spectra obtained over $2014-2016$ were binned on-chip by two in the spatial and spectral directions resulting in a pixel size of $0\aas@@fstack{\prime\prime}27-0\aas@@fstack{\prime\prime}34$ and spectral resolution of $R\sim 4600$ with a sampling rate of 22  km s^-1 pixel^-1 (${\rm FWHM}\sim 65$  km s^-1). The spectra obtained in 2010 were binned only spatially on-chip by two.

The standard echelle package in IRAF was used to combine, to perform flat-field correction, and to extract the ESI spectra. The wavelength solutions were derived using a list of known sky-lines having vacuum wavelengths, where our wavelength solutions have a rms scatter of $\sim 0.03$ Å or about 2  km s^-1. The spectra were also heliocentric velocity corrected.

The galaxy rotation curves were extracted following the method described in Kacprzak et al. (2010) with a similar approach used by Vogt et al. (1996) and Steidel et al. (2002). In summary, we adopted a three-pixel-wide aperture size and shifted the aperture by one pixel intervals along the spatial direction and extracted a series of spectra along the major axis of each galaxy. We performed Gaussian fits to galaxy emission lines (mainly $\rm{H}\alpha$ and $\rm{H}\beta$), which provided the wavelength centroids used to derive the galaxy systemic redshifts and rotation curves. The galaxy redshifts are listed in Table 1. Fig. 2 shows the extracted rotation curve of a galaxy at $z_{\rm gal}=0.20419$ associated with Hi absorption in J$113910-135043$, where the $\rm{H}\alpha$ emission line was used to extract this rotation curve. In this particular geometry, the quasar sightline is positioned in the negative direction along the slit along the galaxy major axis.

Our sample contains 58 quasars, with some quasars probing multiple galaxies and some galaxies having multiple quasar sightlines. The background quasars in each field were observed with the HST/COS and Table 2 provides details of the quasar spectroscopy with the coordinates, redshifts, HST program IDs, and the COS gratings used. The far-ultraviolet gratings G130M and/or G160M have a moderate resolving power of $R\sim 20,000$, giving a full width at half maximum of $\sim 18$  km s^-1 and wavelength coverage of $1410-1780$ Å. We used the STScI CALCOS V2.21 pipeline (Massa & et al., 2013) to reduce and flux calibrate all spectral data acquired from the HST archive. All spectra are heliocentric velocity corrected and in vacuum wavelengths. We co-added multiple integrations with the IDL code coadd_x1d¹¹1http://casa.colorado.edu/~danforth/science/cos/costools.html (Danforth et al., 2010) and binned spectrally by three pixels to enable an increased signal-to-noise ratio. Continuum normalisation was performed by fitting low-order polynomials to the spectra while excluding regions with strong absorption lines.

We implemented the interactive SYSANAL code (Churchill, 1997; Churchill & Vogt, 2001) to define the velocity bounds of Ly$\alpha$ absorption profiles, compute the optical depth-weighted mean systemic redshifts of absorption ($z_{\rm abs}$), and to compute the rest-frame equivalent widths ($W_{r}({1215})$). Column densities were adopted from Sameer et al. (2024) and are listed in Table 3. They use a cloud-by-cloud, multi-phase, Bayesian ionisation modelling approach to determine the physical properties of the absorption systems. It has been demonstrated to produce reliable column densities even when a saturated Ly$\alpha$ is the only Hi line available (Sameer et al., 2021). We do note however, that our results do not depend on the accuracy of the Hi column densities as we have selected our highest data bins to account for saturation of Ly$\alpha$. Where column densities were not available in the literature, we used VPFIT²²2http://www.ast.cam.ac.uk/~rfc/vpfit.html (Carswell & Webb, 2014) to measure the Hi column densities by fitting Voigt profile models to the absorption lines. We used the appropriate line spread function³³3https://www.stsci.edu/hst/instrumentation/cos/performance/spectral-resolution at the corresponding lifetime position when fitting the data. The Hi column densities for all absorption systems are listed in Table 3.

Fig. 1 shows the distribution of rest-frame equivalent widths and column densities as a function of the impact parameter. Both show a strong anti-correlation between the absorption strength and $D$. While high column density systems tend to exist only within the inner halos of galaxies, lower column density systems tend to reside at low and high impact parameters.

We developed a new method for measuring the co-rotation fraction of the Hi halo. For each Ly$\alpha$ absorption system, we compute the fraction of the total equivalent width that is consistent with our co-rotation model. The only dependent choice required is the velocity window we consider when determining whether the gas is consistent with co-rotation. We choose a velocity window that includes all gas from the systemic velocity onward in the direction of galaxy rotation towards the quasar sightline, defined as $f_{\rm EWcorot}$. A value of 1 indicates all of the gas is consistent with a co-rotation scenario, while 0 suggests none of the gas is consistent with a co-rotation scenario. Errors on the co-rotation fraction were calculated by bootstrapping the errors associated with galaxy redshift and the absorption profile and range from $0.001-0.008$.

Fig. 2 (middle) shows the galaxy rotation curve. In this galaxy-quasar pair, the quasar resides on the negative side of the slit position, where the galaxy’s rotation is redshifted with respect to its systemic velocity. For our velocity window criterion, we include all the absorption between the galaxy systemic velocity and the most positive velocity of the absorption boundary defined by SYSANAL as highlighted in pink in Fig. 2 (bottom). In this case, 69% of the absorption equivalent width is consistent with a co-rotation model. This value is comparable to other works where they state that the bulk of the absorption is consistent with co-rotation models (e.g., Ho et al., 2017).

Our new method still allows for a comparison to other works even though different variations of kinematic methods are used (e.g., Steidel et al., 2002; Kacprzak et al., 2010; Ho et al., 2017; Kacprzak et al., 2019a; French & Wakker, 2020) since they discuss co-rotation in a binary form and only when the majority/bulk of the gas is consistent with the model is it co-rotating. Here, we can state that the bulk of the absorption is co-rotating when $f_{\rm EWcorot}\geq 0.5$. Our result now provides a quantification of the amount of gas that is consistent with a co-rotation model. However, we do note that the $f_{\rm EWcorot}$ should be considered an upper limit, and it is plausible that the true co-rotation fraction could be lower since there is the possibility of selecting gas at higher velocities than the galaxy maximum rotation velocity.

Simulations show that the CGM has a vast range of Hi column densities and that different column density regimes may probe different components of the CGM (Fumagalli et al., 2011; Suresh et al., 2019). For example, simulations have shown that Lyman Limit Systems (LLSs) may be the best probe of CGM gas flows, including accretion and outflows (van de Voort et al., 2012; Faucher-Giguère & Kereš, 2011; Faucher-Giguère et al., 2015; Hafen et al., 2017). Here we explore how the gas co-rotation fraction of Ly$\alpha$ behaves as a function of Hi column density.

Figure 3 presents $f_{\rm EWcorot}$ as a function of Hi column density. The grey data points are individual absorption system column densities. We find that the low column density systems span the full range of $f_{\rm EWcorot}$ whereas higher column density systems, particularly those above ${\log(N(\hbox{{\rm H}\kern 1.00006pt{\sc i}})/{\rm cm}^{-2})}\sim 16$ tend toward higher $f_{\rm EWcorot}$. To better investigate this trend, we divided the data into three bins of column density: ${\log(N(\hbox{{\rm H}\kern 1.00006pt{\sc i}})/{\rm cm}^{-2})}<14.5$, $14.5\leq{\log(N(\hbox{{\rm H}\kern 1.00006pt{\sc i}})/{\rm cm}^{-2})}<16.2$, and ${\log(N(\hbox{{\rm H}\kern 1.00006pt{\sc i}})/{\rm cm}^{-2})}\geq 16.2$. The large pink squares represent the mean $f_{\rm EWcorot}$ in each column density bin where the vertical error bars are calculated using 10,000 bootstrapped realisations of the data and their errors to measure the mean and its 1$\sigma$ error. We find that the $f_{\rm EWcorot}$ is correlated with ${\log(N(\hbox{{\rm H}\kern 1.00006pt{\sc i}})/{\rm cm}^{-2})}$, where the co-rotation fraction increases from 0.42$\pm$0.06 for the lowest column density bin to 0.59$\pm$0.05 in the highest column density bin. Therefore stronger absorbers are more likely to have kinematics that are consistent with having gas with line of sight velocities aligned with the rotation curve of the galaxy, increasing from 40% to 60% co-rotation.

Quasar absorption line studies have shown that Hi column densities decrease with increasing impact parameter (e.g., Tumlinson et al., 2013; Borthakur et al., 2015; Kacprzak et al., 2021). Combined with our previous results showing the co-rotation fraction may also be dependent on column density, we explore how it behaves as a function of impact parameter as well as column density. To do this, we bifurcated our sample into two sub-samples with splits at $\log(N({\hbox{{\rm H}\kern 1.00006pt{\sc i}}})/{\rm cm}^{-2})=14.5$ and $\log(N({\hbox{{\rm H}\kern 1.00006pt{\sc i}}})/{\rm cm}^{-2})=16.2$. A column density of 14.5 appears to be where the transition occurs between the CGM and IGM (e.g., Rudie et al., 2012; Wakker et al., 2015; Bouma et al., 2021) and so significant differences in the co-rotation fraction above and below this value can indicate whether the galaxy influences its surroundings beyond $R_{\rm vir}$. The higher column density cut at 16.2 was selected since this is the lower limit for partial Lyman limit systems (pLLs) following the classification by Lehner et al. (2018), which is where a bimodality in CGM metallicity is the most apparent (Lehner et al., 2013; Wotta et al., 2016) and where simulations suggest that signatures of gas flows may become dominant.

Figure 4(a) shows the first-order polynomial fits to $f_{\rm EWcorot}$ as a function of impact parameter for absorption systems bifurcated by $\log(N({\hbox{{\rm H}\kern 1.00006pt{\sc i}}})/{\rm cm}^{-2})=14.5$. The dark and light purple fits present the low (dashed line) and high (dotted line) column density system subsamples, respectively. The $1\sigma$ error on the fits is measured by bootstrapping the fit and calculating the average and standard deviation of 10,000 realisations of the data and their errors. We find that the $f_{\rm EWcorot}$ for lower column density systems, which reside at larger impact parameters and are more IGM-like, appear to be decreasing but the error bars in the slope are consistent with a flat distribution ($f_{\rm EWcorot}=(-0.17\pm 0.26)\log(D/{\rm kpc})+(0.80\pm 0.58)$). Higher column density systems, which span a large range of impact parameters, also have a distribution that is consistent with being flat ($f_{\rm EWcorot}=(-0.04\pm 0.13)\log(D/{\rm kpc})+(0.65\pm 0.22)$).

Fig. 4(b) shows the first-order polynomial fits to $f_{\rm EWcorot}$ for absorption systems bifurcated by $\log(N({\hbox{{\rm H}\kern 1.00006pt{\sc i}}})/{\rm cm}^{-2})=16.2$. The light and dark pink fits represent low (dashed line) and high (dotted lines) column densities, respectively. The 1 $\sigma$ errors are also measured with the bootstrapping method. We see that the larger column density systems reside within 100 kpc, while the lower column density systems extend to larger impact parameters. We find that the statistical behaviour of low column density ($f_{\rm EWcorot}=(-0.11\pm 0.13)\log(D/{\rm kpc})+(0.72\pm 0.27)$), and high column density ($f_{\rm EWcorot}=(-0.01\pm 0.16)\log(D/{\rm kpc})+(0.64\pm 0.27)$) systems are not significantly different, aside from their extent. The trend in low column density systems (light pink dashed line) shows a slightly decreasing $f_{\rm EWcorot}$ with increasing impact parameter, yet it is consistent with a flat distribution. The high column density systems remain roughly constant with impact parameter (dark pink dotted line).

These trends, or lack thereof, should be taken with caution given that our sample covers a wide range of galaxy masses across 2.5 dex (Fig. 1). In fact, the CGM seems to be self-similar over a large mass range (Churchill et al., 2013a, b), where more massive galaxies host CGM gas out to larger distances, but similar absorption strengths are found at similar fractions of the virial radius across the mass range. Therefore, we normalise the impact parameter by the galaxy virial radius and present the computed values of $D/R_{\rm vir}$ in Table 1.

We investigate this mass dependence in Fig. 5(a), which presents the $f_{\rm EWcorot}$ (purple) versus $D/R_{\rm vir}$ for absorption systems split at a column density of $\log(N({\hbox{{\rm H}\kern 1.00006pt{\sc i}}})/{\rm cm}^{-2})=14.5$. We find that higher column density systems (dark purple dotted line) have a flat distribution ($f_{\rm EWcorot}=(-0.02\pm 0.13)\log(D/R_{\rm vir})+(0.57\pm 0.07)$) that extends beyond the virial radius of the galaxies. For the lower column density systems (light purple dashed line), we find a slightly decreasing trend ($f_{\rm EWcorot}=(-0.20\pm 0.26)\log(D/R_{\rm vir})+(0.45\pm 0.07)$) where the $f_{\rm EWcorot}$ could decrease beyond the virial radius. Overall, both low and high column density systems are roughly consistent with each other over the overlapping $D/R_{\rm vir}$ range.

In Fig. 5(b) high column density systems with ${\log(N(\hbox{{\rm H}\kern 1.00006pt{\sc i}})/{\rm cm}^{-2})}\geq 16.2$ are fitted with a first-order polynomial (dotted dark pink line) that has a positive slope ($f_{\rm EWcorot}=(0.10\pm 0.16)\log(D/R_{\rm vir})+(0.66\pm 0.11)$), yet is still consistent with a flat distribution. We find a slightly decreasing trend for the low column density systems with ${\log(N(\hbox{{\rm H}\kern 1.00006pt{\sc i}})/{\rm cm}^{-2})}<16.2$ plotted as a dashed light pink line ($f_{\rm EWcorot}=(-0.16\pm 0.13)\log(D/R_{\rm vir})+(0.47\pm 0.04)$), showing that the $f_{\rm EWcorot}$ decreases for absorption detected mostly beyond the virial radius of galaxies. This indicates that $f_{\rm EWcorot}$ anti-correlates with impact parameter for low column density absorbers and is flat for high column density systems.

The low significance of these (anti-)correlations could be improved by increasing the sample size. Nevertheless, it seems likely that normalising by the virial radius is an important step in understanding how the gas behaves within galaxy haloes. Regardless of the column density cut selected, a large fraction of the gas ($\sim 60\%$) has kinematics consistent with co-rotation within the virial radius. Outside of the virial radius, there is a decline in the $f_{\rm EWcorot}$. Thus, a key factor that determines how much of the CGM is co-rotating is its location within the halo and not its column density.

Our current picture of the CGM is one in which cool gas enters the galaxy halo preferentially along the galaxy’s major axis and likely accretes onto the galaxy while co-rotating with the disk (Stewart et al., 2011a, b, 2013; Nelson et al., 2016; Stewart et al., 2017; Suresh et al., 2019; Péroux et al., 2020). On the other hand, stellar winds and galaxy feedback will be ejected biconically along the minor axis of the galaxy with higher velocities than the disk (Bouché et al., 2012; Schroetter et al., 2016; Lan & Mo, 2018; Schroetter et al., 2019; Reichardt Chu et al., 2022). Some observations have also shown that the spatial distribution of CGM around galaxies appears to be bimodal (Bordoloi et al., 2011; Bouché et al., 2012; Kacprzak et al., 2012b; Kacprzak et al., 2015; Lan et al., 2014; Dutta et al., 2017; Zabl et al., 2019). In order to further examine this picture and the relationships between gas flows with respect to their host galaxies, we explore how the co-rotation fraction relates to galaxy orientation and its behaviour as a function of column density and distance away from the galaxies.

Fig. 6 presents $f_{\rm EWcorot}$ as a function of $D/R_{\rm vir}$ and the azimuthal angle, $\Phi$. Panels (a) and (b) are similar; however, while panel (a) shows the full sample, panel (b) plots only Hi absorption systems within $R_{\rm vir}$ to focus on gas within the “halos” of these galaxies. The host galaxies are located at $D/R_{\rm vir}=0$ with their projected major axis aligned with $\Phi=0$ deg. Each point represents Hi absorption in a background quasar sightline and their sizes represent the absorption column densities to emphasise where the higher and lower column densities tend to reside. The anti-correlation between the Hi column density and $D/R_{\rm vir}$ is clearly visible (also see Fig. 1). The points are also colour-coded based on the measured $f_{\rm EWcorot}$ in each system. The dark pink points in Fig. 6 have the highest consistency with co-rotation kinematics, while dark green is least consistent with the host galaxy rotation velocity.

From the figure, it is clear that the majority (56%) of our absorption systems reside within $R_{\rm vir}$. On average, the absorption systems beyond $R_{\rm vir}$ tend to have a much lower $f_{\rm EWcorot}$ than within $R_{\rm vir}$, with average $f_{\rm EWcorot}$ values of $0.39$ and $0.58$ outside and within $R_{\rm vir}$, respectively. This is consistent with the idea that gas flows may be more organised nearer to galaxies and even that outflows mostly tend not to escape the halos of galaxies (Oppenheimer & Davé, 2008; Tumlinson et al., 2011, 2013; Stocke et al., 2013). Within $R_{\rm vir}$, absorption systems are consistent with higher $f_{\rm EWcorot}$ and tend to be found within 20 degree of galaxy major and minor axes and cover a large range of Hi column densities. Furthermore, these high corotation fractions extend out to $R_{\rm vir}$ along both the major and minor axes. In the remaining part of this subsection, we further explore how $f_{\rm EWcorot}$ behaves with azimuthal and inclination angles as a function of column density and $D/R_{\rm vir}$.

We next examine how $f_{\rm EWcorot}$ and azimuthal angle behaves with column density cuts. We apply the same bifurcation in column density of $\log(N({\hbox{{\rm H}\kern 1.00006pt{\sc i}}})/{\rm cm}^{-2})=14.5$ and $\log(N({\hbox{{\rm H}\kern 1.00006pt{\sc i}}})/{\rm cm}^{-2})=16.2$. Fig. 7 shows the $f_{\rm EWcorot}$ as a function of azimuthal angle with column densities bifurcated at $\log(N({\hbox{{\rm H}\kern 1.00006pt{\sc i}}})/{\rm cm}^{-2})=14.5$. The dark purple circles and light purple diamonds present the averaged $f_{\rm EWcorot}$ in bins of azimuthal angles (horizontal bars) for absorbers with ${\log(N({\hbox{{\rm H}\kern 1.00006pt{\sc i}}})/{\rm cm}^{-2})>14.5}$ and ${\log(N({\hbox{{\rm H}\kern 1.00006pt{\sc i}}})/{\rm cm}^{-2})\leq 14.5}$, respectively, with $1\sigma$ bootstrap errors (vertical bars). Here the fainter data points in the background present individual systems where the smaller purple circles show the systems with higher column densities (${\log(N({\hbox{{\rm H}\kern 1.00006pt{\sc i}}})/{\rm cm}^{-2})>14.5}$) and the smaller light purple diamonds present absorbers with lower column densities (${\log(N({\hbox{{\rm H}\kern 1.00006pt{\sc i}}})/{\rm cm}^{-2})\leq 14.5}$). The high column density systems only have three bins due to the number of systems and sampling of the azimuthal angles. In general, we find that $f_{\rm EWcorot}$ is higher for high column density systems and lower for low column density systems, except within $\Phi=20-30$ deg of the galaxy major axis where low column density systems have a higher $f_{\rm EWcorot}$. Similarly, we find that the co-rotation fraction ($f_{\rm EWcorot}$) for the $\log(N({\hbox{{\rm H}\kern 1.00006pt{\sc i}}})/{\rm cm}^{-2})=16.2$ cut has a higher average value for the high column density systems ($f_{\rm EWcorot}\sim 0.62$) than for the low column density systems ($f_{\rm EWcorot}\sim 0.46$). This trend remains flat across all azimuthal angles for the 16.2 column density cut (not shown here). We will discuss the implications of these results in the next section.

Given that we found $f_{\rm EWcorot}$ is more dependant on $D/R_{\rm vir}$ than $D$ (see Fig. 5), we explore how $f_{\rm EWcorot}$ and azimuthal angle varies with $D/R_{\rm vir}$. In Fig. 8 we present the $f_{\rm EWcorot}$ as a function of $\Phi$ where the sample is bifurcated at ${D/R_{\rm vir}=1}$. The orange circles and blue diamonds present the averaged $f_{\rm EWcorot}$ in bins of azimuthal angles (horizontal bars) for absorbers detected at ${D/R_{\rm vir}\leq 1}$ and ${D/R_{\rm vir}>1}$, respectively, with $1\sigma$ bootstrap errors (vertical bars). The smaller data points in the background show individual systems where the blue diamonds are absorption detected inside the virial radius and the orange circles are the absorption detected outside the virial. The data show that along the projected galaxy major axis ($\Phi<30$ deg), $f_{\rm EWcorot}$ has the same value inside and outside the virial radius with just over half of the gas consistent with co-rotation. This may be expected if gas accretes along filaments, which are co-rotating/accreting from large-scale structure in the IGM down to the galaxy. As the azimuthal angle increases, $f_{\rm EWcorot}$ diverges. We find a high co-rotation fraction for Hi gas in the ${R_{\rm vir}}$ CGM (orange) that slightly increases to a peak of $0.6$ along the projected galaxy minor axis ($\Phi>75$ deg). In contrast, Hi detected at ${D/R_{\rm vir}>1}$ has a decreasing $f_{\rm EWcorot}$ with increasing $\Phi$, where the value drops to $0.27$ along the projected minor axis. This difference in $f_{\rm EWcorot}$ along the minor axis within and outside the virial radius is a factor of $\sim 2$.

Thus, we find overall that only minor axis absorption (${\Phi\geq 60}$ deg) yields significant variations of $f_{\rm EWcorot}$ with $D/R_{\rm vir}$, whereas the major axis Hi gas appears to have a flat distribution of $f_{\rm EWcorot}$ at all radii. If outflows are the dominant source of CGM gas along the minor axis, then outflowing gas co-rotates within the virial radius and either loses angular momentum with increasing height or fails to get to distances beyond the virial radius.

As our galaxies have a range of inclination angles, we investigate how $f_{\rm EWcorot}$ varies with galaxy inclination. In Fig. 9 we present $f_{\rm EWcorot}$ as a function of galaxy inclination angle, $i$, where the sample is bifurcated at $D/R_{\rm vir}=1$. We have verified that our sample’s distribution of inclination angles is consistent with a random distribution of galaxy inclination angles, and because of this, we have fewer galaxies at low inclination angles. The sample is split into three bins based on the inclination angle of the host galaxies: $i\leq 30$ deg, $30<i\leq 60$ deg, and $i>60$ deg. The orange circles show the averaged $f_{\rm EWcorot}$ in each inclination angle bin for systems detected within the $R_{\rm vir}$ of the host galaxies (smaller pale orange data points in the background), while the blue diamonds present the averaged $f_{\rm EWcorot}$ in inclination angle bins for systems detected beyond the virial radius (smaller pale blue diamonds in the background). The vertical bars present the $1\sigma$ bootstrap errors. We find that the $f_{\rm EWcorot}$ of Hi absorption within ${R_{\rm vir}}$ remains almost constant across all inclination angles within uncertainties. However, there is a possible trend that beyond $R_{\rm vir}$, $f_{\rm EWcorot}$ drops at high inclination angles when compared to low and intermediate inclination angles. In highly inclined galaxies, $f_{\rm EWcorot}$ within the virial radius is a factor of $\sim 2$ higher when compared to Hi gas beyond the virial radius. This result could be due to outflows not being able to travel beyond the virial radius and is consistent with what we found in Fig. 8 for the azimuthal angle trends.

Some simulations predict that cold-mode accretion halts for halo masses of log($M_{\rm h}/M_{\odot})>12$ (Dekel & Birnboim, 2006; Kereš et al., 2009; Stewart et al., 2011a), which could present as lower Hi co-rotation fractions for higher mass galaxies. Therefore, it is important to test how $f_{\rm EWcorot}$ varies with the halo mass. Our sample spans over two decades of halo mass from $10.5<\log(M_{\rm h}/M_{\odot})<12.7$, which allows us to test the simulation prediction. In Fig. 10 (a) and (b) we show the $f_{\rm EWcorot}$ as a function of stellar and halo mass, respectively. The grey data points in the background represent individual galaxies and the grey bars represent the stellar and halo mass errors, respectively. The pink squares are the averaged $f_{\rm EWcorot}$ in mass bins (horizontal bars) with $1\sigma$ bootstrap errors (vertical bars). We do not find a significant dependence of $f_{\rm EWcorot}$ on galaxy stellar or halo mass. The data are consistent with being drawn from a flat distribution. Around $\sim$45% of the Hi gas is consistent with a co-rotation model at masses $M_{\rm h}\geq 10^{12}M_{\odot}$, where simulations predict a truncation or halting of cold-mode co-rotating gas accretion. This could imply that the Hi is consistent with being coupled to the kinematics of the galaxy at all masses and/or that accretion is present in galaxies of all masses for our sample. It is also possible that the kinematic connection at higher masses is due to the motions of the larger scale environment that those galaxies live in.

The Hi column density is an important measure of the CGM as it provides insight into where the gas is located and where it originated from. Rudie et al. (2012) determined that ${\log(N(\hbox{{\rm H}\kern 1.00006pt{\sc i}})/{\rm cm}^{-2})}<14.5$ most probably traces distant IGM gas while cosmological simulations have shown that ${\log(N(\hbox{{\rm H}\kern 1.00006pt{\sc i}})/{\rm cm}^{-2})}>16$ is primarily found in the CGM gas flows within the halos of galaxies (e.g., Fumagalli et al., 2011; Suresh et al., 2019). We found that there is a correlation between $f_{\rm EWcorot}$ and N(Hi) in Fig. 3 where there is an increase in $f_{\rm EWcorot}$ by a factor of $\sim 1.5$ from the lowest column density systems to the highest column density systems. This correlation likely implies that there is some dependence on co-rotation and physical processes in a halo, such as outflows from the galaxy, accretion from the IGM, recycled accretion and diffuse components of the CGM. If CGM gas flows exhibit the bulk of the co-rotating gas, then our results are in line with the simulation predictions that higher column density systems are better tracers of gas flows. However, the Hi column density is strongly correlated with the impact parameter and is a hidden additional parameter not accounted for in Fig. 3.

To examine dependencies with absorption location with respect to the host galaxy, we studied the $f_{\rm EWcorot}$ as a function of $D$ and $D/R_{\rm vir}$. We found no significant difference between the statistical behaviour of low and high column density systems over a large range of impact parameters (see Fig. 4). However, the scatter in these trends could be significant given the large range of galaxy masses in the sample. We removed the impact of the host galaxies’ mass by probing the co-rotation fraction as a function of the virial radius normalised impact parameter (see Fig. 5). We showed that the co-rotation fraction is almost constant within the CGM and drops by a factor of two outside of the virial radii of galaxies. This implies that the Hi is most kinematically connected to galaxies within the halo, where most of the physical processes are expected to occur (e.g., outflows, tidal streams, recycling, accretion, etc.). Thus, both column density and distance from the galaxy play critical roles in where gas is kinematically connected to their galaxies. Nonetheless, there is still 35% of Hi that is consistent with co-rotation out to $3R_{\rm vir}$. This gas could be probing filamentary accretion or larger-scale movements of the local environment.

Our results are consistent with the findings of previous works. French & Wakker (2020) reported a co-rotation fraction of $59\pm 5\%$ for low column density Ly$\alpha$ absorbers. Although we use a new method, and our samples overlap at low column densities, the results do not show any significant differences. The authors also reported a decrease in the co-rotation fraction as a function of the impact parameter. However, what we additionally noted here is that both column density and $D/R_{\rm vir}$ are significant factors in the behaviour of $f_{\rm EWcorot}$, with an important transition occurring at $R_{\rm vir}$, where the co-rotation fraction changes from a roughly constant value to rapidly decreasing.

Our co-rotation fractions are also consistent with those found for Mgii absorption systems (Steidel et al., 2002; Kacprzak et al., 2010, 2011b; Ho et al., 2017; Martin et al., 2019), which is not so surprising since Mgii and Hi likely trace similar structures and densities. The main difference between the findings for the two gas tracers is that Mgii absorption tends to have the majority of the absorption aligned with the rotation of the galaxy (Steidel et al., 2002; Kacprzak et al., 2010; Ho et al., 2017), while Hi exhibits a wide range of co-rotation values (see scatter in Fig. 1). This could be due to the fact that Mgii directly traces higher column density Hi that has some metal enrichment where, as we have shown, this higher column density Hi has a higher co-rotation fraction. On the other hand, the Hi absorption is more spatially widespread than Mgii, and traces a larger range of column densities, which could be tracing CGM/IGM or a diffuse component of the CGM, etc., along the same sightline, which has been shown to occur within cosmological simulations (Churchill et al., 2015; Peeples et al., 2019; Marra et al., 2021; Marra et al., 2022). The lower column density Hi can also be traced by Ovi, which has lower co-rotation values of $\sim 50\%$ (Kacprzak et al., 2019a). We will explore the co-rotation fraction of the metals for our sample in an upcoming paper. Overall, our observations are consistent with those from the literature, which provide a picture of the CGM where we find a stronger kinematic connection to the CGM with higher column densities close to the galaxies and weaker kinematic connection to the CGM with lower column densities further from the galaxies.

Our kinematic analysis of Hi absorption also supports a non-uniform and likely bi-modal picture of CGM around galaxies. By accounting for the distance away from galaxies, we found that the co-rotation fraction increases with increasing azimuthal angle within $R_{\rm vir}$ and decreases with increasing azimuthal angle beyond $R_{\rm vir}$ (Fig. 8). This is further supported by a similar trend for galaxy inclination angle, where the co-rotation fraction decreases the most at large inclination angles for gas outside $R_{\rm vir}$ (Fig. 9). There appears to be a geometric preference for co-rotating gas around galaxies, especially along the minor axis and within $R_{\rm vir}$. Compared to previous results, French & Wakker (2020) also reported a sharp decrease in their measure of co-rotation fraction above $i>70$ deg, which is consistent with our results for gas outside $R_{\rm vir}$. They likely only saw a decrease since the majority of their sample is low column density Hi gas that resides near-to-outside of the virial radius.

Tying together all of our results and motivations from previous works, we further examine the geometric distribution of co-rotating gas. We computed the frequency of absorption systems as a function of azimuthal angle in Fig. 11, focusing on systems dominated by co-rotating gas (e.g., the "bulk" of the gas where $f_{\rm EWcorot}\geq 0.5$, which would be similar to other works) within the virial radius of galaxies. Following the methods of Kacprzak et al. (2012b) and Kacprzak et al. (2015), we model the measured azimuthal angles and their uncertainties for each of the galaxies as asymmetric univariate Gaussian PDFs (see Kato et al., 2002). We then compute the mean PDF of all galaxies as a function of $\Phi$. The mean PDF represents the absorption frequency of co-rotating gas at a given $\Phi$. The resulting PDF plotted in Fig. 11 may be bimodal, where the frequency of highly co-rotating gas within $R_{\rm vir}$ is elevated along the major axis and is highly elevated along the minor axis. This distribution mimics the Mgii and Ovi covering fraction bi-modalities, which were assumed to be caused by accretion along the major axis and outflows along the minor axis (Kacprzak et al., 2012b; Kacprzak et al., 2015).

The most surprising aspect of our results with azimuthal angle is that we see significant co-rotation along the minor axis, where gas is often assumed to be outflowing from the galaxy. Fig. 8 further supports the inference that the minor axis co-rotation is dominated by outflows since we find that $f_{\rm EWcorot}$ diverges with increasing azimuthal angle inside and outside the virial radius. If outflows dominate along the minor axes of galaxies, then they appear to have a higher level of co-rotation within the virial radius, suggesting that outflows travel up to $R_{\rm vir}$, but the drop in co-rotation beyond the virial radius suggests that most outflows do not escape. This is consistent with the relative galaxy–absorption velocities which tend to be below the escape velocity of galaxy halos (Stocke et al., 2013; Mathes et al., 2014).

Our co-rotating outflow signatures are in contrast to previous emission line maps of outflowing gas, which tend to show a rapid decrease in the rotation velocities along the minor axis to within a few $10-20$ kpc. Outflows from local starbursting galaxies suggest that signatures of co-rotation diminish beyond 1 kpc from M82 in CO (Leroy et al., 2015). However, rotating gas was found along the outflow axis in recent Enzo (FOGGIE) simulations (Lochhaas et al., 2023). Mapping the gas around a more normal star-forming galaxy in Mgii, [Oii], and [Oiii] emission, Zabl et al. (2021) did not find any co-rotation signatures beyond 10 kpc. However, the ionisation mechanism for the Mgii and nebular emission is still unknown, so our CGM observations do not necessarily trace the same gas, and the starbursting galaxies have higher star formation rates than our sample, which could result in different outflow kinematic properties. In more comparable work, Martin et al. (2019) reported no correlation between the sign of the Doppler shift of Mgii absorption and the rotation of galaxies along the minor axis. This difference could be due to a difference in Mgii/Hi gas tracers or most likely, how co-rotation fractions are defined and measured. Kacprzak et al. (2012a) reported low metallicity ($\sim-2$dex) multi-phase accretion along the minor axis. They found that the low ionisation ions, like Mgii and Siii, were consistent with co-rotation, while the higher ionisation features, like Ciii, Civ, and Ovi, contained both co-rotating gas as well as a fraction of gas that is inconsistent with co-rotation. This again points to a picture where co-rotation likely depends on Hi column density.

If gas accretion via filaments is instead driving the co-rotation seen along the minor axis, then we would expect high co-rotation fractions both inside and outside the virial radius since the accretion is originating from the IGM. However, Fig. 8 clearly shows this to not be the case. Accretion more likely explains the major axis kinematic and co-rotation fraction trends where we do not see any significant transition in kinematics as a function of $R_{\rm vir}$, which is expected from simulations where co-rotation is seen beyond the virial radius (Danovich et al., 2015; Stewart et al., 2013; Stewart et al., 2017).

From another perspective, the significant drop in co-rotation at larger distances along the minor axis, is interesting in itself. It is intriguing that there is a significant amount of gas that have kinematics opposite to the rotation of galaxies. One possibility is that this gas could arise from ancient debris from previous galaxies interactions, which can produce retrograde orbits or reversed kinematics. For example, N-body simulations of NGC 7252 and NGC 3921 are able to reproduce the reversal kinematics observed in their low Hi column density tidal tails (Hibbard & Mihos, 1995; Hibbard & van Gorkom, 1996). Numerical studies of Milky-way also find counter-rotating material around the Galaxy as a result of interactions between disk galaxies. Modelling of interactions implies that merger or fly-by can produce material with retrograde orbits (Pawlowski et al., 2011). This retrograde motion of satellites seems to be quite common in a sample of $z<0.04$ SDSS galaxies, which was shown that 40% of the galaxies have retrograde motions, which is the same fraction seen in cosmological simulations (Azzaro et al., 2006). Therefore, it is plausible to find low column density gas in retrograde motion at large distances, but it is strange that we see this occurring around the minor axes of our more isolated sample of galaxies. This may only occur if the alignment of the large scale structure is in the direction of the minor axes of these galaxies (Bailin et al., 2008). A larger statistical sample would help to explore this discovery.

Another interesting finding shown in Fig. 7 is that while there is the expected lower $f_{\rm EWcorot}$ for lower column density systems, it is not true along the major axes of galaxies. $f_{\rm EWcorot}$ jumps by nearly a factor of two for lower column density systems, while high column density systems have a lower $f_{\rm EWcorot}$. This could possibly be due to the fragmentation of clouds as they approach the disk or within an accretion stream. One would expect this to occur along the minor axis as well since simulations have shown clouds can fragment as they move through outflows (McCourt et al., 2018; Sparre et al., 2018; Nelson et al., 2020). So it is unclear why we see a transition along the major axis and not the minor axis.

Churchill et al. (2013b) studied the relation between Mgii equivalent width and host galaxies’ virial mass and found no anti-correlation. Despite the prediction of simulations (van de Voort et al., 2011; Stewart et al., 2011a), the strength of low-ionisation CGM absorption is not dependent on halo mass and there is no sudden truncation of cold-mode accretion at a mass of $\log(M_{\rm h}/M_{\odot})=12$. It is interesting that our results for $f_{\rm EWcorot}$ as a function of $\log(M_{\rm h}/M_{\odot})$ are roughly constant. If we assume that the co-rotating gas is associated with cold-mode accretion, then gas accretion is not dependent on halo mass and cold-mode accretion is still occurring in the most massive galaxy halos.

It is plausible that the existence of cool gas for high mass galaxies could be arising from their environment, i.e., arising from the galaxies surrounding them. However, here we find the majority of the gas within $R_{\rm vir}$ is kinematically consistent with galaxy co-rotation. Given that the vast majority of Mgii absorption is found within $R_{\rm vir}$, and that we find most of the gas with $R_{\rm vir}$ is kinematically coupled to the galaxy, then it is less likely that the CGM detected in high mass galaxies arises from larger-scale environments like groups and clusters since they would have larger velocity offsets and could be inconsistent with a co-rotation model. Environmental effects may be more significant for gas outside $R_{\rm vir}$, where there is 35% out at $3~{}R_{\rm vir}$, but this does not represent the bulk of the total absorption seen in the CGM.