# FEASTS Combined with Interferometry. III. The Low Column Density H i Around M51 and Possibility of Turbulent-mixing Gas Accretion

Xuchen Lin, Jing Wang, Lister Staveley-Smith, Suoqing Ji, Dong Yang, Xinkai Chen, Fabian Walter, Hsiao-Wen Chen, Luis C. Ho, Peng Jiang, et al.

► To cite this version:

Xuchen Lin, Jing Wang, Lister Staveley-Smith, Suoqing Ji, Dong Yang, et al. FEASTS Combined with Interferometry. III. The Low Column Density H i Around M51 and Possibility of Turbulent-mixing Gas Accretion. The Astrophysical Journal, 2025, 982 (2), pp.151. 10.3847/1538-4357/adb718. hal-05486745

HAL Id: hal-05486745

https://hal.science/hal-05486745v1

Submitted on 23 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# FEASTS Combined with Interferometry. III. The Low Column Density H I Around M51 and Possibility of Turbulent-mixing Gas Accretion

Xuchen Lin (林旭辰)<fnref n="1" /> ID, Jing Wang (王菁)<fnref n="2" /> ID, Lister Staveley-Smith<fnref n="3" />,<fnref n="4" /> ID, Suoqing Ji (季索清)<fnref n="5" />,<fnref n="6" /> ID, Dong Yang (杨冬)<fnref n="1" /> ID, Xinkai Chen (陈新凯)<fnref n="2" /> ID, Fabian Walter<fnref n="7" /> ID, Hsiao-Wen Chen<fnref n="8" /> ID, Luis C. Ho<fnref n="2" /> ID, Peng Jiang (姜鹏)<fnref n="9" /> ID, Nir Mandelker<fnref n="10" /> ID, Se-Heon Oh<fnref n="11" /> ID, Bo Peng<fnref n="12" /> ID, Céline Péroux<fnref n="13" />,<fnref n="14" /> ID, Zhijie Qu (屈稚杰)<fnref n="8" />,<fnref n="15" /> ID, and Q. Daniel Wang<fnref n="16" /> ID

Received 2024 November 8; revised 2025 February 13; accepted 2025 February 14; published 2025 March 26

## Abstract

With a new joint-deconvolution pipeline, we combine the single-dish and interferometric atomic hydrogen (H I) data of M51 observed by the Five-hundred-meter Aperture Spherical radio Telescope (FAST) (FEASTS program) and the Very Large Array (VLA) (THINGS). The product data cube has a typical line width of 13 km s⁻<fnref n="1" /> and a 2σ line-of-sight (LOS) sensitivity of H I column density $N_{H I} \sim 3.2 \times 18 \text{ cm}^{-2}$ at a spatial resolution of $\sim 18''$ ($\sim 0.7$ kpc). Among the H I detected LOSs extending to $\sim 50$ kpc, $\sim 89\%$ consist of diffuse H I only, which is missed by previous VLA observations. The distribution of dense H I is reproduced by previous hydrodynamical simulations of this system, but the diffuse component is not, likely due to unresolved physics related to the interaction between the circumgalactic and interstellar media. With simple models, we find that these low $N_{H I}$ structures could survive the background ultraviolet photoionization, but are susceptible to the thermal evaporation. We find a positive correlation between LOS velocity dispersion ($\sigma_v$) and $N_{H I}$ with a logarithmic index of $\sim 0.5$. Based on existing turbulent mixing layer (TML) theories and simulations, we propose a scenario of hot gas cooling and accreting onto the disk through a TML, which could reproduce the observed power index of $\sim 0.5$. We estimate the related cooling and accretion rates to be roughly one-third to two-thirds of the star formation rate. A typical column density of diffuse H I ($\sim 10^{19} \text{ cm}^{-2}$) can be accreted within 300 Myr, the interaction timescale previously estimated for the system. Such a gas accretion channel has been overlooked before, and may be important for gas-rich interacting systems and for high-redshift galaxy evolution.

*Unified Astronomy Thesaurus concepts: Interstellar atomic gas (833); Circumgalactic medium (1879); Galaxy accretion (575); Galaxy mergers (608); Deconvolution (1910); Atomic interferometry (119); Single-dish antennas (1460)*

## 1. Introduction

The evolution of galaxies is largely determined by the cosmic baryon cycle (C. Péroux & J. C. Howk 2020), which involves complex hydrodynamic processes (R. A. Crain & F. van de Voort 2023). In particular, the multiphase nature of gas poses challenges to both theoretical and observational understanding of galaxy evolution (T. Naab & J. P. Ostriker 2017). The accretion and depletion of neutral gas regulate the supply of star-forming material, thereby influencing the growth of stellar mass and the formation of disk structures on scales of 1–10 kpc (H. J. Mo et al. 1998). Continuous gas accretion into the interstellar

medium (ISM) (the gas inside galaxy disks) is necessary because the gas within star-forming regions would otherwise be exhausted within a few gigayears (A. Saintonge et al. 2017; C. Péroux & J. C. Howk 2020; L. J. Tacconi et al. 2020).

In theory, gas accretion can occur through various modes, including the cosmological cold- and hot-mode accretion, the redistribution of satellite gas, and the galactic fountain (M. E. Putman 2017; C.-A. Faucher-Giguère & S. P. Oh 2023). Recent multiwavelength observations have revealed that gas is accreted primarily from the intergalactic medium (IGM) (the gas outside dark matter halos) in an ionized phase, while the cooling and condensation into a neutral phase happen later within the circumgalactic medium (CGM) (the gas between IGM and ISM) (M. E. Putman 2017). Possible mechanisms happening during gas accretion and cooling include the precipitation from thermal instability (e.g., G. M. Voit et al.

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Department of Astronomy, School of Physics, Peking University, Beijing 100871, People's Republic of China</footnote>

<footnote n="2">² Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, People's Republic of China; jwang_astro@pku.edu.cn</footnote>

<footnote n="3">³ International Centre for Radio Astronomy Research, University of Western Australia, 35 Stirling Highway, Crawley, WA 6009, Australia</footnote>

<footnote n="4">⁴ ARC Centre of Excellence for All-Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia</footnote>

<footnote n="5">⁵ Center for Astronomy and Astrophysics, Department of Physics, Fudan University, Shanghai 200438, People's Republic of China</footnote>

<footnote n="6">⁶ Key Laboratory of Nuclear Physics and Ion-beam Application (MOE), Fudan University, Shanghai 200433, People's Republic of China</footnote>

<footnote n="7">⁷ Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany</footnote>

<footnote n="8">⁸ Department of Astronomy and Astrophysics, The University of Chicago, 5640 S. Ellis Avenue, Chicago, IL 60637, USA</footnote>

<footnote n="9">⁹ National Astronomical Observatories, Chinese Academy of Sciences, 20A Datun Road, Chaoyang District, Beijing 100101, People's Republic of China</footnote>

<footnote n="10">¹⁰ Centre for Astrophysics and Planetary Science, Racah Institute of Physics, The Hebrew University, Jerusalem 91904, Israel</footnote>

<footnote n="11">¹¹ Department of Physics and Astronomy, Sejong University, 209 Neungdong-ro, Gwangjin-gu, Seoul, Republic of Korea</footnote>

<footnote n="12">¹² Max-Planck Institut für Astrophysik, Karl-Schwarzschild-Straße 1, D-85741 Garching, Germany</footnote>

<footnote n="13">¹³ European Southern Observatory, Karl-Schwarzschild-Straße 2, D-85748 Garching, Germany</footnote>

<footnote n="14">¹⁴ Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388 Marseille, France</footnote>

<footnote n="15">¹⁵ Department of Astronomy, Tsinghua University, Beijing 100084, People's Republic of China</footnote>

<footnote n="16">¹⁶ Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA</footnote>
2017; M. Gaspari et al. 2018), the condensation by mixing (M. Gronke & S. P. Oh 2020a), and the rotating cooling flow (Z. Hafen et al. 2022; J. Stern et al. 2024). However, the details and relative importance of these physical processes remain uncertain (C.-A. Faucher-Giguère & S. P. Oh 2023). Observationally, though not directly tracing the material accreted from the IGM, the neutral gas in CGM serves as an indicator of high-density, low-temperature regions, providing insights into the overall efficiency of multiphase gas accretion (P. Richter 2017; D. Nelson et al. 2020).

Due to their low column density, cold (<10⁴ K) gas components within the CGM have been primarily studied using absorption lines. Statistically, these absorption line detections sample a large area around the associated galaxies, even extending beyond the virial radius (e.g., N. Lehner et al. 2020; M. C. Wilde et al. 2021; M. A. Berg et al. 2023), in line with the large-scale structures traced by neutral gas in simulations (e.g., F. van de Voort et al. 2019; R. Ramesh et al. 2023). However, the current sample size is not large enough to fully explore the dependence of absorption line detections on diverse galaxy properties and other factors, such as the group environment and tidal interactions (H.-W. Chen et al. 2020; C. Péroux et al. 2022; M. A. Berg et al. 2023). One example is the dependence on the azimuthal angle of absorbers around galaxies, which is predicted by theories (C. Péroux et al. 2020) but still inconclusive in observations (R. Bordoloi et al. 2011; G. G. Kacprzak et al. 2015; S. Weng et al. 2023). The geometric degeneracies and the incomplete galaxy optical information complicate the association of absorption systems with galaxies, and further make it difficult to interpret their distances and relative velocities (C. L. Martin et al. 2019). Near the disk–CGM interface, the saturation of Lyα absorption lines at intermediate neutral atomic hydrogen (HI) column densities ($N_{HI} \sim 10^{17}-10^{19} \text{ cm}^{-2}$) and the limited availability of high-quality Lyman limit system observations pose challenges to understanding the ISM–CGM connection (M. A. Berg et al. 2023; V. V. Klimenko et al. 2023; S. Weng et al. 2023).

Directly mapping HI around galaxies across a wide range of column densities is essential to solving many of the aforementioned challenges and uncertainties. The 21 cm emission HI mapping offers a wide field of view around galaxies, covering the various conditions and environments where multiphase gas resides. If the sensitivity is sufficient, it can directly trace the final stages of transition into HI and accretion from the CGM onto the ISM disk. The CGM–ISM interface is where matter, energy, and momentum are actively exchanged and recycled, driven by processes like thermal conduction, radiative cooling, and turbulent mixing (C.-A. Faucher-Giguère & S. P. Oh 2023), which could be powerfully constrained with the direct mapping of HI column densities and velocities. In addition, direct mapping is not limited by the incidence of quasars, making the process of sample selection much simpler than absorption studies. The direct mapping would also help to refine simulations that have difficulties in linking the cosmological accretion history with small-scale gas physics (M. Vogelsberger et al. 2020; C.-A. Faucher-Giguère & S. P. Oh 2023).

With the high sensitivity of the Five-hundred-meter Aperture Spherical radio Telescope (FAST), the FAST Extended Atlas of Selected Targets Survey (FEASTS; J. Wang et al. 2023, hereafter W23)<fnref n="17" /> is able to map the HI in the local Universe

down to a 3σ column density limit of $N_{HI} \sim 5 \times 10^{17} \text{ cm}^{-2}$, nearly reaching the $N_{HI}$ range of Lyman limit systems. W23 detected large-scale low $N_{HI}$ gas associated with the NGC 4631 group, showing the potential of FAST to bridge the gap between ISM and CGM studies. With a larger sample of 11 systems, a follow-up study by J. Wang et al. (2024b, hereafter Paper I) estimated the HI number density based on the assumption of pressure balance with the CGM. They found that the low-density, diffuse component of HI missed by interferometers is ubiquitous, and that the diffuse fraction increases with the strength of tidal interactions. This low-density faint HI is likely a direct result of CGM–ISM material exchange, potentially enhanced by tidal interactions.

Tidally interacting galaxies offer unique laboratories for studying the physical processes connecting the CGM and ISM in the context of both galaxy evolution and multiphase gas dynamics. Observational and theoretical studies have found that intense tidal interactions during mergers could enhance star formation (H.-A. Pan et al. 2019; S. L. Ellison et al. 2022), and at the same time, the gas mass remains constant (S. Quai et al. 2021; Q. Yu et al. 2022; W. Li et al. 2023) or even increases (S. L. Ellison et al. 2015, 2018; J. Moreno et al. 2019; J. Shangguan et al. 2019) with possibly diluted metallicity (M. D. Thorp et al. 2019); simulated and observed major merger remnants also exhibit rejuvenation (M. Sparre & V. Springel 2017). These findings suggest rapid evolution in galaxy mass assembly, driven by active gas cooling, inflow, and consumption at various merger stages. Mergers are known to drag large and long-lived tidal tail structures (e.g., A. Chung et al. 2007; M. Zhu et al. 2021), enlarging the interface between cold HI structures and the hot CGM and possibly catalyzing gas exchange between them. The widespread morphology, perturbed kinematics, and rapid evolution of merging systems offer a wide range of gas conditions, making them ideal for testing gas flow-related physical scenarios. These systems may also serve as analogs for galaxies at higher redshifts, where gas-rich mergers are frequent on average (C. J. Conselice 2014; C. Péroux & J. C. Howk 2020; Q. Duan et al. 2024), and are likely a significant channel of gas accretion.

M51a (NGC 5194, aka the Whirlpool Galaxy) is in the late stage of merging with (M51b (NGC 5195), as supported by simulations reproducing the interferometrically observed HI distribution and kinematics (S. Howard & G. G. Byrd 1990; C. Theis & C. Spinneker 2003; C. L. Dobbs et al. 2010). Strong tidal effects induced the formation of evident spiral arms and an extended tail, observed in both stellar light (W. Parsons 1850) and gas emission (A. H. Rots et al. 1990; J. Pety et al. 2013). Previous studies have mostly focused on the inner disk of M51a (e.g., D. Colombo et al. 2014; S. Oikawa & Y. Sofue 2014; S. Manna & S. Roy 2023), where the galaxy’s spiral structure is rather regular and of grand design, and multi-wavelength data are abundant. Recent FAST observations have revealed extended low $N_{HI}$ gas surrounding the pair (Paper I; see also H. Yu et al. 2023), especially at the unexpected west and south sides. Compared to other galaxies that FEASTS has already observed (Paper I), the faint HI in M51 is not closely attached to the bright and dense HI detected in interferometric observations (A. H. Rots et al. 1990; F. Walter et al. 2008), but extends much further (see Figure 15 in Appendix F). Without nearby satellites, this large extension of relatively diffuse HI gas likely originates from the strong tidal interaction. The diffuse HI may represent an intermediate phase between the

<footnote n="17">¹⁷ Data are publicly available at https://github.com/FEASTS/LVgal/wiki.</footnote>
colder gas and hot CGM, and in this study, we explore its role
in the evolution of this system.

Despite the high sensitivity, FAST is limited by its ~3' spatial resolution due to its single-dish nature. One long-standing goal in radio astronomy is to combine low-resolution single-dish H I data with high-resolution interferometric data, which usually miss total flux because of the short-spacing problem (S. Stanimirovic 2002; A. Plunkett et al. 2023). Among the many possible methods, jointly deconvolving the interferometric dirty cube and single-dish H I data is the one most likely to coherently recover both the spatial and kinematic information, but has rarely been applied beyond the Local Group (e.g., T. J. Cornwell 1988) due to many difficulties. Especially, the apparent size of H I structures decreases significantly at large distances, making them difficult to resolve with single-dish telescopes that have either large beams (e.g., the Green Bank Telescope) or significant side lobes (e.g., Arecibo, R. Giovanelli et al. 2005).

This series of studies focuses on combining FEASTS with interferometric data, leveraging FAST's relatively small beam and low side-lobe level. In the first two papers (Paper I; J. Wang et al. 2024a, hereafter Paper II), the cleaned interferometric models of 11 galaxies are compared or linearly combined with their FAST cubes, which contains the diffuse H I missed by interferometers. The diffuse H I was found to rotate more slowly and have a higher velocity dispersion than the dense one, possibly forming a thick layer extending vertically into the CGM. In this paper (Paper III), we conduct a detailed study of M51, tackling its unexpectedly extended H I structure. We develop a new pipeline using the maximum entropy method (MEM, e.g., S. F. Gull & G. J. Daniell 1978) to jointly deconvolve the FAST cube and the Very Large Array (VLA) data obtained by The H I Nearby Galaxy Survey (THINGS, F. Walter et al. 2008), simultaneously achieving the deep sensitivity and high resolution of the two data sets. Paper IV (D. Yang et al. 2025) presents a statistical study of the extraplanar H I around edge-on galaxies.

This paper is organized as follows. The data sets and analysis methods are introduced in Sections 2 and 3, including the deconvolution pipeline and decomposition of diffuse gas components. We present the results in Section 4. The velocity dispersion is compared pixelwise with column density, and the H I number density and thickness are calculated assuming pressure balance with the CGM. In Section 5, we discuss the possible mechanisms of the H I cooling and accretion driven by turbulent mixing, and summarize in Section 6. Throughout this paper, the variables with a subscript of “H I” include only the contribution from hydrogen, and a factor of 1.4 is applied when helium needs to be considered. Meanwhile in Section 5, the density of CGM, and thus the cooling rate, takes the helium into account so that we can compare them with the star formation rate (SFR).

## 2. Data

We take the coordinates and distance of M51a as those adopted by THINGS, $\alpha_{2000} = 202.4696$, $\delta_{2000} = 47.1952$ (J. L. Turner & P. T. P. Ho 1994) and $D = 8.0$ Mpc (I. D. Karachentsev et al. 2004). Its stellar mass ($10^{10.73} M_{\odot}$) and SFR ($4.5^{+2.6}_{-1.6} M_{\odot} \text{ yr}^{-1}$) are taken from the $z=0$ Multi-wavelength Galaxy Synthesis (zOMGS; A. K. Leroy et al. 2019) as in Paper I, and the B-band optical size $R_{25} = 11.6$ kpc ($5.0$) (G. Paturel et al. 2003). We calculate its virial mass

$M_{\text{vir}} = 10^{12.37} M_{\odot}$ from the stellar mass, using the relation of P. S. Behroozi et al. (2010). Here, we do not use the dynamical mass ($10^{12.190} M_{\odot}$) from E. Kourkchi & R. B. Tully (2017) as in W23 because for M51, they detected a small number of group members (only six). They estimated a mass of $10^{12.685} M_{\odot}$ from the K-band luminosity, and the value we use serves as a middle value. The coordinates of M51b are taken from the NASA/IPAC Extragalactic Database (NED), $\alpha_{2000} = 202.4983$ and $\delta_{2000} = 47.2661$ (M. F. Skrutskie et al. 2006).

### 2.1. FAST Data

The FAST observations of the FEASTS sample were carried out using basket-weaving on-the-fly scans. The data reduction pipeline implements the standard procedures of single-dish H I data processing, and is optimized for the FEASTS data (W23). The H I data cubes have a beam of 3.24 FWHM, and a channel width of $1.610 \text{ km s}^{-1}$. The FEASTS observations of M51 reach an rms level of $0.80 \text{ mJy beam}^{-1}_{\text{FAST}}$, corresponding to a $3\sigma$ column density detection limit of $4.0 \times 10^{17} \text{ cm}^{-2}$, assuming a line width of $20 \text{ km s}^{-1}$. More details of the observations of M51 can be found in Paper I.

### 2.2. VLA Data

The main interferometric data we use are the calibrated and continuum-removed visibility obtained by the THINGS (F. Walter et al. 2008), which covers a heliocentric radial velocity range between 214.9 and $704.5 \text{ km s}^{-1}$ with a channel width of $5.153 \text{ km s}^{-1}$ (96 channels). The continuum subtraction of the THINGS visibility is improved as described in Appendix A. We verify that the contribution of the w-components (A. R. Thompson et al. 2017) is insignificant in this visibility data set.

The naturally weighted dirty cube is generated using MIRIAD task invert, during which we apply an FWHM $= 3''$ Gaussian taper to increase the signal-to-noise ratio (S/N). The pixel size is the same as that of THINGS products, $1.5$, and the image size is $1800 \times 1800$ pixels ($45' \times 45'$). The dirty cube has a noise level of $\sigma_{\text{int}} = 0.379 \text{ mJy beam}^{-1}_{\text{VLA}}$, corresponding to a $3\sigma$ column density limit of $4.6 \times 10^{19} \text{ cm}^{-2}$, assuming a line width of $20 \text{ km s}^{-1}$. The corresponding clean beam has a shape of $17.84 \times 15.40$ with position angle of P.A. $= 67.52$, or a physical resolution of $\sim 0.7$ kpc.

In Figure 1(d), we plot the cleaned spectrum of this corrected THINGS visibility as red dashed-dotted line. Using MIRIAD CLEAN, we largely follow the procedures of F. Walter et al. (2008), including the rescaling of residual, with a cleaning threshold of $1.5\sigma$.

## 3. Analyses

### 3.1. Jointly Deconvolving the FAST and VLA Cubes

We developed a pipeline to jointly deconvolve the VLA dirty cube and FAST cube per channel using MEM. MEM seeks the positive solution model image $m_j$ ($j$ denotes the serial number of pixels) that simultaneously

1. fits the input single-dish or dirty image, i.e., the rms value of residual map $r_j$ is smaller than a threshold $\theta$, and
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The moment maps (a)–(c) and the spectrum (d) of the jointly deconvolved data. The centers of M51a and M51b are marked as black crosses in panels (a)–(c). The stellar disk size ($R_{25}$ = 11.6 kpc) is indicated by the black ellipse, with P.A. = 172° and i = 42° (D. Tamburro et al. 2008, as adopted by THINGS). The gray contour encloses the region with a moment 2 value higher than $10^{19}$ km s⁻¹ ≈ 80 km s⁻¹, and the white contour indicates the 2σ confidence LOS column density limit, ~3.2 × 18 cm⁻². In the lower left corner of panel (a), the open white circle displays the diameter of the instrumental resolution (~57′′) at $N_{H_1}$ ~ 3.2 × 18 cm⁻², and the white ellipse within is the restoring beam (17′′84 × 15′′40). The actual resolution transits between these two limits, reaching the best one when $N_{H_1}$ > 1.6 × $10^{19}$ cm⁻². In panel (d), the spectrum of the jointly deconvolved cube (blue solid) is compared with those of the single-dish FEASTS (orange dashed) and VLA THINGS (red dashed-dotted; recleaned). As reference, the spectrum of the FEASTS cube after regridded into the interferometric coordinates is shown as the green dotted line. The gray shades indicate the channel range of continuum subtraction for THINGS data. The central velocity of $v_{sys}$ = 460.0 km s⁻¹ is taken from the velocity modeling (see Section 3.2).</figcaption></figure>

2. has the maximum information entropy

$$H = -\sum_{j} m_{j} \ln \frac{m_{j}}{e\hat{m}_{j}} = -\sum_{j} m_{j} \left( \ln \frac{m_{j}}{\hat{m}_{j}} - 1 \right) \quad (1)$$

relative to an input reference image $\hat{m}_j$, which is referred to as “a priori image” by T. J. Cornwell & K. F. Evans (1985) and as “default image” in the MIRIAD task mosmem (R. J. Sault 1990). Traditionally, this reference image is instead set to be a flat image with the total flux of the single-dish image.

Condition 1 is the basic chi-square requirement for a possible solution, and condition 2 forces the model to trace the reference image $\hat{m}_j$ when the constraint from condition 1 is weak. In this pipeline, we use the MIRIAD task mosmem (T. J. Cornwell & K. F. Evans 1985; T. J. Cornwell 1988; R. J. Sault 1990;

R. J. Sault et al. 1996), a common implementation of MEM deconvolution for radio astronomy.

Unlike traditional MEM applications (see the summary by S. Stanimirovic 2002), this new pipeline has these two characteristics:

1. A non-flat reference cube $\hat{m}_j$ is passed to the MEM deconvolution, along with the interferometric dirty cube $d_j^{int}$ and the single-dish cube $d_j^{sd}$.

2. We iteratively invoke the MEM deconvolution. The reference $\hat{m}_j$ is initially set to be $d_j^{sd}$, and is adaptively adjusted between each iteration according to the output model $m_j$ and residuals.

By this pipeline, we more exhaustively exploit the available information in the data, and reduce and flatten the residuals more thoroughly than before. The noise level of the jointly deconvolved result is expected to be between the
interferometric ($\sigma_{\text{int}}$) and single-dish ($\sigma_{\text{sd}}$) noise levels, and we estimate it as $\sqrt{\sigma_{\text{int}}\sigma_{\text{sd}}}$. The step-by-step technical details of the pipeline are given in Appendix B for readers who are interested.

The output model $m_j$ is convolved with a restoring beam chosen to be the VLA clean beam (17.84 × 15.40), without adding back either of the two residuals. Figure 1 gives the moment maps and spectrum of the jointly deconvolved M51 H I data. No cube mask is applied except when we calculate the moment 2 map, which therefore appears smaller. The reason is that the line-of-sight (LOS) velocity dispersion measurement is more sensitive to noise than the column density is. The cube mask is generated with SoFiA 2 (P. Serra et al. 2015; T. Westmeier et al. 2021) from the model cube plus a Gaussian noise of $\sqrt{\sigma_{\text{int}}\sigma_{\text{sd}}}$ at the resolution of restoring beam.

Our pipeline is justified by the mock tests in Appendix C.2, which show that the accuracy of moment 2 map increases significantly with a SoFiA mask, and that at the resolution of restoring beam (~18", or ~0.7 kpc), the 2$\sigma$ column density limit is ~3.2 × 18 cm⁻² for a single LOS (white contour in Figure 1(a)). This $N_{\text{HI}}$ limit is roughly the geometric mean of the sensitivities of THINGS and FEASTS (Section 2), justifying our choice of $\sqrt{\sigma_{\text{int}}\sigma_{\text{sd}}}$ as the estimated noise level. In Appendix C.3, we also model the instrumental resolution of the jointly deconvolved cube at a given $N_{\text{HI}}$ level using the noise Fourier spectra of THINGS and FEASTS data. We find that the instrumental resolution at the above $N_{\text{HI}}$ limit (~3.2 × 18 cm⁻²) is ~57" (2.2 kpc), and that a resolution of restoring beam (~18") can be reached at $\gtrsim 1.6 \times 10^{19}$ cm⁻². The mock-verified effective resolution (~18") is much better than the instrumental one (~57"), likely because the intrinsic spatial scale for the H I gas at ~3.2 × 18 cm⁻² is much larger than the instrumental resolution. In Figure 1(a), we plot at the lower left corner both the restoring beam and a circular ~57" beam as the resolution reference.

We define our LOS sample as the LOSs with a finite moment 2 measurement ($\sigma_v$, the LOS velocity dispersion) smaller than $10^{1.9}$ km s⁻¹ ≈ 80 km s⁻¹, i.e., the colored regions out of the gray contour in Figure 1(c). We set the upper limit as 80 km s⁻¹ because these high $\sigma_v$ regions may be related to shocks or the projection of multiple components in the northern region, where M51 encounters M51b.

The total H I mass we jointly deconvolved is $10^{9.59} M_\odot$, the same as the FEASTS value (Paper I) and 1.55 times the THINGS-only value. We largely recovered the extra diffuse H I flux detected by FAST relative to the THINGS data (Figure 1(d)). In Appendix C.4, we compare the moment 0 map from the joint deconvolution with those of FEASTS, recleaned THINGS cube, and the result of linear combination. The channel maps of the model and residual cubes can be found in Appendix G. They consistently show that the joint deconvolution pipeline produces a data cube with smooth transition from high- to low-density regions.

## 3.2. Removing Projection Effects from the Velocity Field

In later analyses of gas entrainment, one key parameter to be used is the velocity of H I moving in the CGM, which is reflected in the moment 1 map (Figure 1(b)) but suffers from projection effects. The velocity fields both within and beyond $R_{25}$ show dipolar and spider-like morphology, indicating the possibility of partially removing the projection effects by modeling the velocity field with tilted rings, as done previously

<figure><img src="image_5.png" /><figcaption>**Figure 2.** The profile of fitted rotational ($v_{\text{rot}}$), radial ($v_{\text{rad}}$), and total ($v_{\text{tot}}$) velocities from the tilted ring modeling using BBarolo. The value of $v_{\text{tot}}$ (thick red solid line) is calculated as $\sqrt{v_{\text{rot}}^2 + v_{\text{rad}}^2}$. The crosses are the fitted values of $v_{\text{rot}}$ (blue) and $v_{\text{rad}}$ (green), and the dotted lines are the median-filtered profiles; the filter sizes are 50" and 210" for $v_{\text{rot}}$ and $v_{\text{rad}}$, respectively. By convention, a positive $v_{\text{rad}}$ in BBarolo indicates inflowing for a galaxy rotating counterclockwise as M51a. The value of the sound speed in CGM ($c_s$, CGM = 156 km s⁻¹) is indicated by the horizontal gray dashed line.</figcaption></figure>

in other large-scale interacting systems (e.g., A. Sorgho et al. 2019). The kinematics of a tidally perturbed system cannot be well described by rings, and the outer H I structure often does not consist of full rings (see Section 4.2), but to the first order, the modeling provides a best-effort correction for projection.

With BBarolo (or 3DBarolo) (E. M. Di Teodoro & F. Fraternali 2015) 3DFIT task, we conduct the fitting in four rounds before obtaining the final results. Each model consists of 130 rings with a radial increment of 10". In the first round, these rings are divided into three groups with a radial range of <300", 300"–570", and >570", ensuring relatively constant P. A.'s within each group. These three groups roughly correspond to three visually identified regions with distinct velocity patterns: the inner disk, the transitional region, and outer tidal structures. We fit the rings in each group separately, and in later rounds, we cancel the grouping. Results of each round are fixed or used as the initial guess for the next round. In the first round, the initial values of central position, systemic velocity ($v_{\text{sys}}$), and inclination (i) take the kinematic value used by THINGS (D. Tamburro et al. 2008), while those of rotational velocity ($v_{\text{rot}}$) and P.A. are visually guessed. The radial velocity $v_{\text{rad}}$ is set as zero initially. The variation ranges of i and P.A. are 30° and 60° in the first two rounds, and decrease to 5° and 15° in the third round. The ring centers and $v_{\text{sys}}$ are fixed after the first round. In the final round, $v_{\text{rot}}$ and $v_{\text{rad}}$ are fitted with all other parameters fixed. By doing so, we gradually tighten the constraints on the model toward the final results.

Two parameters crucial to our later analyses are the total velocity $v_{\text{tot}} = \sqrt{v_{\text{rot}}^2 + v_{\text{rad}}^2}$ and the inclination i. In Figure 2 are the fitted velocity profiles, and the full fitting results are presented in Appendix D. The fluctuation reflects the complexity of this system and the degeneracy between velocities and geometric parameters. To suppress the fluctuation, we smooth the profiles of $v_{\text{rot}}$ and $v_{\text{rad}}$ with median filters of 50" and 210" wide, respectively. The $v_{\text{tot}}$ profile (red solid line in Figure 2) is calculated from the smoothed profiles of $v_{\text{rot}}$ and $v_{\text{rad}}$ (dotted lines). We construct spatial maps of $v_{\text{tot}}$ and i based on the tilted
ring model for later pixelwise analyses. The details of the map construction are given in Appendix D.

### 3.3. Separating Diffuse and Dense H I

One advantage of the M51 system is that a large fraction of H I detected LOSs have a low H I column density. These LOSs are likely to be mainly composed of diffuse H I. Physically, we expect the diffuse H I to have relatively low number densities, high velocity dispersions, high thermal temperatures, and large variance scales. Its properties should lie between those of the CGM hot gas and the “dense” H I typically captured by interferometry (W23). Its temperature may be of ~5 × 10³ K, the typical value of the warm neutral medium detected in the Milky Way (MW; N. M. McClure-Griffiths et al. 2023), so as to roughly remain neutral and thermally stable (G. B. Field et al. 1969).

Because the data have a limited spatial resolution, the observed diffuse H I possibly includes unresolved small-scale dense H I cloudlets (M. Gronke & S. P. Oh 2020b). The observation thus provides an ensemble measure of gas properties (M. Gaspari et al. 2018), meaning that instead of dealing with the behavior of individual cloudlets, we focus on the collective properties and behavior of the gas. The unresolved dense H I cloudlets either are insignificant in mass compared to the diffuse H I, or else behave as an ensemble like the diffuse H I (like cloudlets in a fog/mist). While the large-scale H I distribution and motion in the M51 system should be effectively shaped by the tidal interaction (but see Section 4.2), the localized survival of H I depends on the interaction with the immediate radiative and thermodynamic environment of the CGM (in the off-plane, less self-shielded direction). This environment leaves signatures in the ensemble H I properties, and can be compared with the prediction of simple models in the assumption of quasi-equilibrium.

With these considerations, we develop a simple procedure to separate the LOS H I spectrum into the broad and narrow components, which corresponds to the more diffuse (warmer) and denser (cooler) gas parts. We select bright channels along each LOS, and fit them with one or two Gaussian component(s) plus a flat positive baseline. These Gaussian models are relatively arbitrarily required to have $\sigma_v \lesssim 35$ km s⁻¹ and peak intensity above $2\sigma_{int}$. Two Gaussian components are used only when two distinct peaks above $1.5\sigma_{int}$ are identified; otherwise, only one component is used. The detailed procedure is given in Appendix E.

We separate the H I into the following three components, and refer to them as

1. *Diffuse-only H I.* No fitting is conducted in these LOSs. The spectral peak is smaller than $2\sigma_{int}$, corresponding to a maximum $N_{H\text{I}}$ of ~9.5 × 10¹⁹ cm⁻² × ($\sigma_v$/18 km s⁻¹).

2. *Dense H I* refers to the fitted Gaussian components. For the LOS spectra with two dense components (only 7% of the dense H I bearing LOSs), the one with a higher column density is denoted as the main component, and the other is the secondary one.

3. *Dense-associated diffuse H I* is the residual of the fitted dense H I component. This component and the diffuse-only H I are referred to as *diffuse H I* together.

Among the LOS sample, 29% of the mass is diffuse-only, 20% is dense-associated diffuse, and 51% is dense. In the same sample, ~89% of the LOSs are associated with diffuse-only

H I. Because this study focuses on LOS properties, the spectral decomposition's goal is equivalently reduced to simply separating the diffuse-only LOS from other LOSs. This simple goal mitigates the uncertainties coming from the procedure.

In the following, the LOSs of diffuse-only H I are used as the fiducial sample for discussion, while the dense-associated diffuse H I is only shortly discussed in Sections 4.3 and 4.4 for the benefit of completeness and consistency check. In the future, when extending this study to a sample of galaxies, we will conduct more-careful LOS decomposition using tools like BAYGAUD (S.-H. Oh et al. 2022).

## 4. Results

We present the results of this paper in this section. We first present the global kinematics and morphologies of the H I components that we separated, especially the H I in the diffuse-only region. We then focus on the statistical distributions of the LOS parameters, particularly those related to the diffuse H I, including the relation between the velocity dispersion and H I column density, and the radial profiles of the number density and H I thickness. These relations reflect the physical processes that regulate these parameters, such as the heating or cooling of diffuse H I.

### 4.1. Inspecting the Moment Maps of H I in the Inner Region

We report in Figure 3 the results of our double-Gaussian decomposition. The dense and diffuse H I (panels (a) and (b)) divides roughly at a column density of $10^{19.7}$ cm⁻². As expected from the decomposition criteria, the spatial distribution of dense H I is similar to the interferometric observations (Appendix C.4), and its velocity field (panels (c) and (g)) resembles that of the total H I (Figure 1(b)). Both moment 1 maps have an inflow pattern that the zero-relative-velocity line and the minor axis are misaligned,<fnref n="18" /> while the dense-associated diffuse H I does not have a similar misaligned pattern (panel (h), within the gray contour).

Paper I found that the diffuse H I in the inner region rotates at a smaller velocity than the dense one does, and the velocity difference shown in Figure 3(f) exhibits this “lagging” pattern at a better resolution. There are some fine structures on the velocity difference map at the east or west sides of the dense H I disk, possibly relating to tidal stripping and shocks (also see Section 5.5). Nonetheless, the overall behavior of the dense-associated diffuse H I is corotating with the dense H I, while lagging behind it with a median velocity offset of 12 km s⁻¹ (Figure 3(f)), or a relative extent of 23% (Figure 3(e)).

As discussed in Paper I, the lag in velocity often indicates extraplanar gas accretion through the galactic fountain mechanism (F. Fraternali 2017). It was also found that the denser gas tends to be more responsive to spiral arms and thus to form radial inflows (D. Colombo et al. 2014). These patterns confirm the complex but typical physical processes happening close to the star-forming disk. In the following, we more focus on the low column density H I extending radially far into the CGM beyond the star-forming disks. This outlying H I extends the radial range for the disk radial inflow to happen, and serves as a useful test of H I–CGM interaction.

<footnote n="18">¹⁸ The misalignment is due to the $\sin\phi\sin\phi$ projection of inflow velocities, where $\phi$ denotes the azimuthal angle from the major axis.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The moment 0 and 1 map properties of the dense and diffuse H I. The locations of M51a and b are marked as crosses. (a) The moment 0 map of total dense H I. The LOS that consists of two dense components is enclosed by orange contours. (b) The moment 0 map of diffuse H I. The negative values from noise are blanked. (c) and (d) The moment 1 maps of the main dense H I component and the diffuse H I. The zoomed-in version is in panels (g) and (h). The dense-associated diffuse H I is enclosed by the gray contour in panel (f). (e) The ratio of LOS velocity, relative to the central systemic velocity of $v_{\text{sys}} = 460.0 \text{ km s}^{-1}$, in the diffuse H I to the dense H I. (f) The LOS velocity difference between the dense and diffuse H I.</figcaption></figure>

## 4.2. The Position–Velocity Distribution of Tidal Structures and the Possible Merger History

To infer the possible origins of the large-scale H I structure, we plot the position–velocity (PV) diagrams in Figures 4(a)–(i) along the 1' wide section lines shown in Figure 4(j). The labeled white circles in panel (j) mark the left-hand edge of the corresponding PV diagrams in panels (a)–(i).

We identify three ridges along the tidal tail ((a)–(c)). Their PV diagrams indicate that the tail is a single continuous structure stemming from the southwestern edge of the dense H I disk. The asymmetric wings in panels (a) and (b) in the velocity direction imply the lagging or leading of diffuse H I by tidal and/or ram pressure forces. The structure (d) has a similar PV shape and may have been broken from the major tidal tail.

The H I features traced by lines (e)–(i) are more complicated in their morphology. Along lines (e) and (f), the velocity is almost constant at ~475 km s⁻¹, slightly redshifted relative to the M51a center ($v_{\text{sys}} = 460.0 \text{ km s}^{-1}$). The blueshifted feature on the right side of panel (f) is from the projection of a different structure that is also traced by line (h). Panels (g) and (i) show a positive velocity gradient toward the north, where their lines recede and connect with the inner disk.

In summary, from the inner H I disk protrude out two relatively coherent groups of tail-like structures ((a)–(d) and (e)–(i)). Those in the former are elongated and linked, while those in the latter are intermittent and segmented. These two groups become kinematically similar when reaching the dense disk in the west. It is possible that at the beginning, these two groups were produced together during one major encounter between the two galaxies, but then the second has been further disrupted by another major encounter, as indicated by previous simulations (see Section 4.2.1).

### 4.2.1. The Comparison with the Hydrodynamic Simulation of C. L. Dobbs et al. (2010)

C. L. Dobbs et al. (2010) made a hydrodynamical simulation of M51 based on the orbital initial conditions from C. Theis & C. Spinneker (2003). They largely recovered the H I structure detected by A. H. Rots et al. (1990) with VLA, which was similar to the THINGS image (shown in Figure 13(b) of Appendix C.4), but they missed the radially extended structure in our jointly deconvolved image of Figure 1(a). The simulation did not include a CGM halo and therefore any hydrodynamic CGM–ISM interaction. Thus, it can help to distinguish the roles of tidal interaction with M51b (when the simulation successfully reproduced structures) and hydrodynamic effects of H I–CGM interactions (when the simulation failed to reproduce structures) in forming the observed large-scale H I structure.

In their simulation, M51a formed two large tidal tails as M51b approached from the north ~300 Myr ago (see the top middle panel of Figure 4 in C. L. Dobbs et al. 2010). One of them later became the tail currently observed to the east of M51a (structures (a)–(d)). However, there are two significant differences from the observation, which possibly come from the absence of the interaction with CGM: (1) their simulated tidal tail opens more widely, and (2) the simulated system lacks a diffuse envelope that should have extended to ~50 kpc.

In the simulation, ~120 Myr ago, M51b pierced through and thus disrupted the other tail when both M51b and the tail orbited to the south (the top right panel of their Figure 4). Since then, a great portion of the shattered tail was dragged to the north, but a narrow stripe extending to the west also formed in the simulation. This predicted structure roughly corresponds to but is very different from the broad, intermittent, and segmented structures labeled as (e)–(i) in Figure 4(j), which
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Identified structures in M51 (j) and their corresponding PV diagrams (a)–(i). The structures are divided into three groups and identified with different colors: (a)–(d) (blue), (e)–(f) (orange), and (g)–(i) (violet). The redshift of M51a center ($v_{sys} = 460.0$ km s⁻¹) is plotted on the PV diagrams as gray dotted lines. In panel (j), white circles mark the line end corresponding to the left side of PV diagrams.</figcaption></figure>

extend to the south and then bend eastward. Especially, structures (e) and (f) have no obvious velocity gradient, implying that the radial velocities have previously decreased, possibly due to the drag of ram pressure or newly accreted CGM gas, both of which were not considered by C. L. Dobbs et al. (2010).

In summary, the encounter simulated in C. L. Dobbs et al. (2010) was largely consistent with the dense H I structures, but missed most of the new diffuse H I features shown in FEASTS data, possibly because they did not include the CGM-ISM interaction, which is the major theme of this study. In the simulation, the interaction started ~300 Myr ago, and the second encounter happened ~120 Myr ago. These results set a rough time range for the diffuse H I to grow and survive.

### 4.3. The Relation between H I Gas Column Densities and Velocity Dispersions

From now on, we focus on the statistical behavior of localized properties, such as the column density $N_{H\ i}$ and the velocity dispersion $\sigma_v$. We use the previously defined LOS sample (with $\sigma_v < 10^{1.9}$ km s⁻¹ ≈ 80 km s⁻¹) in the following analyses, and exclude the stellar disk (i.e., $r < R_{25}$).

A change in $N_{H\ i} = n_{H\ i}h/\cos i$ might come from the variation in either H I number density ($n_{H\ i}$) or thickness ($h$), both of which are related to the local turbulence level or velocity dispersion. In Figure 5(a), we plot the pixelwise distribution of $N_{H\ i}$ and $\sigma_v$ of the jointly deconvolved cube with contours enclosing 90%, 75%, 50%, and 25% of the LOS sample. Here, $\sigma_v$ is represented by the value of moment 2 map. The distribution for diffuse-only H I is given in Figure 5(b).

Both distributions in panels (a) and (b) show a positive relation between $\sigma_v$ and $N_{H\ i}$ when $N_{H\ i} \lesssim 10^{20}$ cm⁻². The upper envelopes of the distributions in Figures 5(a) and (b) are less sharp than the lower envelopes, possibly because of the projection of structures along the same LOS, leading to the overestimation of $\sigma_v$ at a given $N_{H\ i}$. This trend is driven by the behavior of diffuse-only H I because the dense H I is detected only when the total $N_{H\ i} > 10^{19.7}$ cm⁻².

We use bisector least-squares linear fitting to evaluate the ridge lines of the distributions in Figures 5(a) and (b), using the relatively confident LOSs within the 75% contour. We plot the result as the black dashed line and orange dashed-dotted line in panels (a) and (b), respectively, with a slope of $0.548^{+0.078}_{-0.037}$ and $0.585^{+0.090}_{-0.019}$. This final fitting result is the bisector line of the following two fittings: (1) we fit $\sigma_v$ against $N_{H\ i}$ with the weighting of $1/n$, where $n$ is the number of LOSs in each log $N_{H\ i}$ bin, and (2) we fit $N_{H\ i}$ against $\sigma_v$, weighted by the inverse of the LOS number in log $\sigma_v$ bins. The systematic uncertainties of slopes are estimated as the range of variations with the threshold contour level changing from 90% to 50% (gray and blue solid lines in two panels). We confirm with mock tests that the random uncertainties from the individual values of $N_{H\ i}$ and $\sigma_v$ do not contribute significantly (<0.020) to the slope uncertainty. For reference, the bisector-fitted line in panel (a) is repeated in all four panels.

This slope of ~0.5 coincides with that of self-gravitating gas like the dense molecular clouds (e.g., M. Heyer et al. 2009; A. K. Leroy et al. 2015; J. Sun et al. 2018, 2020b), but the $N_{H\ i}$ expected for self-gravitating at a given $\sigma_v$ should be much higher (the gray dashed line in panels (a) and (c)), inconsistent with the actual distribution. The H I in this regime is mostly
<figure><img src="image_8.png" /></figure>

**Figure 5.** The pixelwise distribution of H I column density $N_{H\ i}$ and the moment 2 map value of $\sigma_v$ of (a) the jointly deconvolved cube, (b) diffuse-only H I, (c) dense and dense-associated diffuse H I, and (d) the FAST cube, given as the filled colored contours. In panel (c), the dense and dense-associated diffuse H I corresponds to the lower and upper branches, respectively, and the distribution before decomposition is given as black contours. In all panels, the contours enclose 90%, 75%, 50%, and 25% of the corresponding data from the outside in. We conduct bisector least-squares linear fitting between the $\sigma_v$ and $N_{H\ i}$ for pixels within the 75% contour in panels (a) and (b). The fitted line of panel (a) and (b) are plotted as the black dashed line (repeated in all panels) and orange dashed-dotted line, respectively. The slopes are roughly 0.5, reported at the corner of each panel. The uncertainties of slopes are the variation range when the threshold contour changes from 90% to 50% (the gray and blue solid lines in two panels). In panels (a) and (c), we plot as the gray dashed line the anticipated $\sigma_v-N_{H\ i}$ relation if the turbulence pressure is balanced by the self-gravity of gas. The gas thickness is assumed to be 1 kpc as an upper limit (see the thickness estimation in Section 4.4); the line would move down with a smaller thickness.

suggesting the influence of beam-smoothing effects. Therefore, the joint deconvolution is essential for recovering the H I distribution and revealing its correct properties.

located outside the optical disk, and should be strongly regulated by CGM–H I interactions.

The distribution in panel (a) flattens at $N_{H\ i} \gtrsim 10^{20} \text{cm}^{-2}$, where the dense H I is more dominant in mass. In this regime, processes originating in the stellar disks, e.g., the gravity and stellar (and possibly active galactic nuclei) feedback, should play a more important role than the CGM–H I interaction in regulating the relation (e.g., M. R. Krumholz et al. 2018). A deeper investigation of this regime is beyond the scope of this paper.

After we decompose the dense and diffuse components (Figure 5(c)) in the dense H I detected region, the $\sigma_v-N_{H\ i}$ distribution clearly separates into two branches. The diffuse branch (the upper one) roughly follows the same trend of the positive correlation between $\sigma_v$ and $N_{H\ i}$. However, the diffuse branch effectively consists of the residuals of the Gaussian fitting for dense H I, therefore having a low S/N for a confident $\sigma_v$ measurement. We thus refrain from making further interpretations.

In Figure 5(d), the trend of FAST data is shallower compared with the jointly deconvolved one (black dashed line),

### 4.4. The Radial Distribution of Number Density and Thickness of H I

The radial distribution of column density $N_{H\ i}$ is shown in Figure 6(a) for different H I components of the LOS sample, after we exclude the stellar disk ($r < R_{25}$). The $N_{H\ i}$ profiles are largely flat at the ~$10^{19} \text{cm}^{-2}$ level beyond $2 R_{25}$, echoing our earlier point that this system is an interesting laboratory to study the H I–CGM interaction. In comparison, the slope of the M51 $N_{H\ i}$ profile is much shallower than the median radial profile (gray dotted line in panel (a)) of H I in galaxies based on interferometry data (J. Wang et al. 2020), the profiles of nine other galaxies that have both THINGS and FEASTS data (see Appendix F, data from Paper II), and most of the galaxies observed by FEASTS (J. Wang et al. 2025).

In Figure 6(b), the velocity dispersion $\sigma_v$ of diffuse H I decreases with radius. The trend seems like the extension of
<figure><img src="image_9.png" /><figcaption>**Figure 6.** The radial profiles of the following H I properties: (a) column density $N_{H\ i}$, (b) velocity dispersion $\sigma_v$, (c) volumetric density $n_{H\ i}$, and (d) thickness $h$. The radius $r$ from the center is normalized by optical $R_{25} = 11.6$ kpc, after considering the inclination. The diffuse-only H I is given as the violet shading, and the “Y”-shaped symbols and error bars are the median and 16th and 84th percentiles in each $0.2\ R_{25}$ wide bin. Similarly, the orange dots and filled green stars correspond to dense H I and dense-associated diffuse H I, respectively. The orange contour encloses 90% of the dense H I pixels. In panel (a), the median profile of total H I is the white dashed line. The statistical median H I profile from J. Wang et al. (2020) is the gray dotted line, linearly extrapolated to $2\ R_{25}$; the H I radius $R_{H\ i}$ is defined at an H I surface density of $1\ M_{\odot}\ pc^{-2}$. In panel (b), the horizontal gray dotted line is $10\ km\ s^{-1}$ for reference. In panel (d), the minimal thicknesses to survive thermal evaporation ($h_{evap,c}$) and photoionization ($h_{pi,c}$) are given as the red dashed line and the filled red pentagons with error bars, respectively.</figcaption></figure>

what has been observed in the $N_{H\ i} > 10^{20}\ cm^{-2}$ regime (e.g., C. Bacchini et al. 2019). A non-negligible fraction (20%) of LOSs at $r > 3\ R_{25}$ have $\sigma_v < 10\ km\ s^{-1}$ (gray dotted line in panel (b)), implying that the diffuse H I with so low a column density is not very efficiently heated by the CGM and UV background. We will discuss this further in Section 4.4.1.

Assuming that the H I turbulence pressure $1.4m_p n_{H\ i} \sigma_v^2$ is in equilibrium with the sum of CGM pressure and midplane H I gas weight,<fnref n="19" /> we calculated per pixel the volume density $n_{H\ i}$ and the thickness $h = \cos(i)N_{H\ i}/n_{H\ i}$ for different H I components at $r > R_{25}$. Although the whole galaxy is a dynamically unstable interacting system, the localized gas tends to be in thermodynamic equilibrium. Any localized inequilibrium should dissipate within a sound-crossing time of ~10 Myr across the thickness, much shorter than the merging timescale. Here, $m_p$ is the mass of proton, and the local inclination angle $i$ is taken from the tilted ring modeling (Section 3.2 and Appendix D). The CGM pressure, calculated as $p_{CGM} = \rho_{CGM}k_B T_{CGM}/\mu_{CGM}$, is three times higher than the

H I gas weight within $1.5\ R_{25}$ and has a larger contribution beyond. We take the average particle weight $\mu_{CGM} = 0.6m_p$. The virial temperature ($T_{CGM} = 1.06 \times 10^6$ K) and gas density profile ($\rho_{CGM}$) are calculated from the $M_{vir}$, assuming a typical $\beta$ model following Appendices D and E of W23.<fnref n="20" /> We ignore the insignificant contribution (<10%) of H I thermal pressure. We assume that the cool gas is mostly neutral, which is validated a posteriori in Section 4.4.1. The magnetic pressure is not considered because it is not well constrained observationally (e.g., C. L. Van Eck et al. 2015; J. L. Han 2017; A. S. Borlaff et al. 2021; S. Manna & S. Roy 2023), implying an unknown level of over- and underestimation of $n_{H\ i}$ and $h$, respectively.

The results are shown in Figures 6(c) and (d). Both dense and diffuse H I have an almost constant volumetric density, especially out of $2.5\ R_{25}$. The diffuse $n_{H\ i}$ fluctuates around $2 \times 10^{-2}\ cm^{-3}$, and the dense $n_{H\ i}$ is ~1 dex higher, consistent with the theoretical values typical of their respective column densities (A. Rahmati et al. 2013). The constant $n_{H\ i}$ implies a constant cooling rate of diffuse H I, possibly serving as the cooling interface between hot CGM and dense H I (also see Paper I, based on coarser measurements).

<footnote n="19">¹⁹ We calculate the H I gas weight felt by the diffuse H I as $\pi G \Sigma_{\text{diffuse}} \times (\Sigma_{\text{diffuse}}/2 + \Sigma_{\text{dense}})$, i.e., its self-gravity plus the gravity from embedded dense H I (see Equation (7) of E. C. Ostriker et al. 2010). The weight felt by dense H I is calculated as $\pi G(\Sigma_{\text{diffuse}} + \Sigma_{\text{dense}})^2/2$, i.e., the total H I gas weight, because the dense H I is much thinner than the diffuse one; a demonstrative diagram can be found in J. Sun et al. (2020a, Figure A.1). Here, G is the gravitational constant, and $\Sigma_{\text{diffuse}}$ (or $\Sigma_{\text{dense}}$) denotes the deprojected diffuse (dense) H I surface density including the contribution of helium.</footnote>

<footnote n="20">²⁰ The virial mass $M_{vir}$ is assumed to have a critical overdensity of 101. We estimate the total CGM gas mass using Equation (5) of S. Ettori (2015). The density profile has $\beta = 0.64$ and a core radius of $0.15\ R_{500}$, where $R_{500}$ is the radius within which the average density is 500 times the critical density.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The distribution of diffuse-only LOSs in the H I thickness $h$-volumetric density $n_{H\ i}$ space, after we exclude the stellar disk and high $\sigma_v$ regions. The critical thickness $h_{pi,c}$ for a neutral ratio of 50% due to photoionization is plotted as the red dashed line. In total, 78% of LOSs and 91% of diffuse-only H I mass are self-shielded from photoionization.</figcaption></figure>

The H I thickness $h$ decreases with radius, contrary to the H I flaring phenomenon previously found at smaller radii based on a static hydrodynamic modeling of the balance between the gas turbulence pressure and the off-plane gravity (e.g., C. Bacchini et al. 2019). For a flaring disk, the $h$ increases with $r$ due to the increasingly weaker gravitational confinement. From our calculation, the underlying reasons for this contrary trend are the relatively high CGM pressure and the low H I velocity dispersion. The CGM takes over the confining role, pressing the H I to rather high volume densities, and the low velocity dispersion of H I implies the existence of efficient kinetic energy dissipation or cooling (Section 5.1).

### 4.4.1. The Critical Thicknesses of Photoionization and Thermal Evaporation

There is the possibility that the diffuse H I surrounding M51 could be destroyed by the CGM through thermal conduction or photoionization. Meanwhile, the existence of these H I structures indicates that the net heating is not that efficient in causing a change in the phase distribution. In this section, we estimate criteria for these two possible mechanisms to take effect.

**Photoionization.** We estimate the condition for significant photoionization following W23, using their Cloudy simulation results (G. J. Ferland et al. 1998). The simulation assumes that the gas is smooth; therefore, any derived critical condition for surviving photoionization should be loosened if the actual gas has clumpy structures. Given that the leakage of UV photons along the spiral disk plane is weak, we assume that the ionizing photons in the outer disk (>R<sub>25</sub>) mainly come from the cosmic background radiation, described by the ionization parameter U as 10<sup>-6</sup> cm<sup>-3</sup>/n<sub>H I</sub> (J. Tumlinson et al. 2011) with the default Cloudy background setting.

In Figure 7, we plot the critical thickness for a neutral ratio of 50% ($h_{pi,c}$, red dashed line) over the $h-n_{H\ i}$ distribution. With the stellar disk and high $\sigma_v$ regions excluded, 78% of LOSs and 91% of diffuse-only H I mass are self-shielded from photoionization. The diffuse H I is reasonable to be treated as neutral when we calculate the pressure balance in Section 4.4. In Figure 6(d), we also provide the median $h_{pi,c}$ within bins of

width 0.2 R<sub>25</sub>, which roughly follows the lower envelope of the majority of $h$ values.

**Thermal evaporation.** We take the L. L. Cowie & C. F. McKee (1977) model to calculate the minimal survival thickness ($h_{evap,c}$) that makes the saturation parameter $\sigma_0$ = ($T_{CGM}/1.5 \times 10^7$ K)<sup>2</sup>/(n<sub>H,CGM</sub>h<sub>evap,c</sub>/cm<sup>-3</sup>pc) reach the critical value of 0.027 (C. F. McKee & L. L. Cowie 1977), above which there would be evaporation flow.

The $h_{evap,c}$ profile is plotted in Figure 6(d) as the red dashed line. Beyond 2.5 R<sub>25</sub>, most of the H I gas is possibly susceptible to evaporation. Following L. L. Cowie & C. F. McKee (1977), we estimate that at 38% of the diffuse-only LOSs, the H I should evaporate in less than ~120 Myr, the time since the last encounter between M51a and M51b (Section 4.2.1).

In this evaporation model, the H I is assumed to exist as individual spherical cloudlets with a size of ~$h$; if the H I is a slab, the evaporation could be much less effective (L. L. Cowie & A. Songaila 1977). However, the turbulent nature of cool gas in the CGM tends to rule out the slab morphology (H.-W. Chen et al. 2023). Another possible caveat is that the magnetic field could significantly suppress the conduction (L. Spitzer 1962; S. I. Braginskii 1965; R. Kooij et al. 2021) if its orientation mainly follows the surface of H I gas. However, the latest observations suggest that magnetic lines loop perpendicular to the H I disk plane (Y. Stein et al. 2019). If we use calibrated equations of criteria from the latest magnetohydrodynamics simulations (e.g., Figure 5 of Z. Li et al. 2020), the obtained critical thickness is very similar.

One possibility is that strong cooling processes are at work, and are able to compensate for the thermal conduction loss and sustain the extended diffuse H I. We propose a scenario of turbulent-mixing gas cooling and accretion in Section 5.

## 5. Turbulent-mixing Gas Accretion

In this section, we propose a scenario whereby the ISM–CGM interaction is governed by a turbulent mixing layer (TML; M. C. Begelman & A. C. Fabian 1990). The TML theory highlights the efficient radiative cooling of the intermediate-temperature (~10<sup>5</sup> K) gas in TML, which is generated from the mixing at the shearing interface between the cool- and hot-phase gas. In our proposed scenario, the H I not only efficiently condenses from the CGM, but also flows inward to fuel the star formation. Globally, these two processes correspond to the gas accretion onto H I disk and into the star-forming optical disk. This process may help the diffuse H I to survive thermal conduction, and to replenish the galactic H I reservoir.

We emphasize that this is just one possible scenario, but it has largely been overlooked in the literature and is thus interesting to explore. It is worth pointing out that the TML is a physical component generally unresolved in the existing cosmological simulations, and is difficult to directly observe in observations (C.-A. Faucher-Giguère & S. P. Oh 2023). Our study assumes a simple model, and obtains a possible ensemble-averaged view of the phenomenon. Specifically, tidal interaction in the M51 system dragged out its large-scale H I structure, which serves as a seed for future cooling due to shearing in TMLs. Keeping in mind that other physical processes may be at play to influence the results, we nonetheless obtain interesting constraints, which may hopefully inspire future efforts with multiwavelength observations and theoretical recipes.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The schematic plots of the turbulent-mixing gas cooling and accretion scenario. (a) The face-on view of the cool H I disk rotating at a velocity of $v_s$. (b) The edge-on view. The shearing between the H I disk and the CGM forms TML, where the mixed gas at $T_{\text{mix}} \sim 10^5$ K efficiently cools, and is entrained into the H I disk at a velocity of $v_{\text{ent}}$ perpendicular to the disk. The gas condensed from the CGM has lower angular momentum, and would flow in. The inflow rate is denoted as $\dot{M}_{\text{in}}$. (c) A slice of the temperature field in the strong-cooling regime (see the main text in Section 5.1). The direction of $v_s$ and $v_{\text{ent}}$ are shown by the white arrows. The TML thickness $h_{\text{turb}}$ is bracketed by black arrows, and $\sigma_{\text{turb}}$ is the velocity dispersion within the layer. (d) The same as panel (c), but for the weak-cooling regime. The well-mixed region is labeled with thickness $h_{\text{mix}}$ and dispersion $\sigma_{\text{mix}}$. The temperature slices in panels (c) and (d) are taken from Figure 7 of T21.</figcaption></figure>

## 5.1. The Scenario

This paper aims to provide a feasibility test for the turbulent-mixing gas cooling and accretion on galaxy scale, rather than an observational verification of the scaling relations from the numerous 3D idealized simulations of TML (e.g., L. Armillotta et al. 2016; D. Padnos et al. 2018; S. Ji et al. 2019, hereafter J19; N. Mandelker et al. 2019, 2020; D. B. Fielding et al. 2020; M. Gronke & S. P. Oh 2020a; B. Tan et al. 2021, hereafter T21; Y. Yang & S. Ji 2023; H. Aung et al. 2024). These simulations have explored in different geometric configurations the relations between the kinematic and geometric properties of TML, and the rate of hot gas cooling and entrainment. A conclusive census on the exact scaling relations has still not been reached, but these studies share similar physical insights, and give similar numerical predictions for the TML properties.

We summarize the key points of our scenario as the cartoons in Figure 8. We consider the shearing (at speed $v_s$) between the cool H I disk ($T_{\text{cool}} \lesssim 10^4$ K) and the hot and rather static CGM ($T_{\text{hot}} \gtrsim 10^6$ K) from the rotation of the former (Figure 8(a)). Due to the Kelvin-Helmholtz instability, TML forms between them with a temperature of $T_{\text{mix}} \sim \sqrt{T_{\text{cool}}T_{\text{hot}}} \sim 10^5$ K (Figures 8(c) and (d)). From the equation of momentum conservation, the thickness $h_{\text{turb}}$ and LOS velocity dispersion

(in the direction of $h_{\text{turb}}$) $\sigma_{\text{turb}}$ of TML are related as (Z. Kuncic & G. V. Bicknell 2004; M. C. White et al. 2016; J19)

$$ \frac{h_{\text{turb}}}{L} \sim \frac{\sigma_{\text{turb}}^2}{v_s^2}. \qquad (2) $$

Here, L is the typical scale of TML along the shearing interface.

When $v_s$ is subsonic within the CGM$^{21}$ (Y. Yang & S. Ji 2023), the strong cooling at $T_{\text{mix}}$ causes the well-mixed gas in TML to condense into gas of $T_{\text{cool}}$ in the H I phase, which we directly observe in this work. After merging with the cold layer, the condensed H I would flow in because of their lower angular momentum inherited from the CGM. In a dynamical equilibrium, the cooling and condensation rates are expected to be comparable to the radial inflow rate. We expect that the condensed H I replenishes and makes up a significant part of the diffuse H I, and their properties, such as thickness and dispersion, may correlate with those of TML.

T21 categorized the cooling process within TMLs into the weak and strong regimes, according to the relative importance of radiative cooling and turbulence, which is characterized as the ratio of their timescales, i.e., the Damköhler number

<footnote>$^{21}$ The speed of sound in M51 CGM is $156 \text{ km s}^{-1}$, larger than most of the shearing velocity (see Section 3.2).</footnote>
Da = t<sub>turb</sub>/t<sub>cool</sub>.<sup>22</sup> The TMLs in these two regimes have distinct physical structures, and their entrainment rates of hot gas show different scaling relations. In the following, we briefly summarize features of each regime, and discuss in each case properties of possible H I observables, such as the column density N<sub>H1</sub> and LOS velocity dispersion σ<sub>v</sub>.

The strong cooling (Da > 1). In this regime, the cooling is efficient compared with the turbulent diffusion, resulting in a thin but highly corrugated mixing and cooling interface, which separates the TML into distinct cool- and hot-phase regions. Most of the cooling processes take place in the thin wrinkled interface (Figure 8(c)), of which the increased interface area helps the global cooling rate to increase. The cool gas occupies roughly half of the strong-cooling TML volume (T21), and shares the turbulence level at the TML scale.

In this case, the H I that we observe traces the cool-phase region of the TML, approximately satisfying

$$N_{\text{H1}} = 2n_{\text{H I,cool}} \times \frac{1}{2}h_{\text{turb}}, \quad (3)$$

$$\sigma_v = \sigma_{\text{turb}}, \quad (4)$$

when the TML accounts for a significant part of the H I. The factor 2 accounts for both sides of the shearing, and 1/2 is the approximate filling factor of cool gas in the TML (also see Figure 6 of T21). Here, the measured N<sub>H1</sub> and σ<sub>v</sub> are the averaged value within the ~0.7 kpc beam.

The H I in the TML is in direct pressure equilibrium with the hot gas in and outside the TML. Their pressures are mainly contributed by turbulence and thermal pressure, respectively, i.e.,

$$p_{\text{CGM}} = \rho_{\text{H I,cool}} \sigma_{\text{turb}}^2 = 1.4 m_{\text{p}} n_{\text{H I,cool}} \sigma_{\text{turb}}^2. \quad (5)$$

The weak cooling (Da < 1). In this regime, the cooling is inefficient, making the TML thicker than the cooling interface, which is on the cold gas side of the TML (Figure 8(d)). It is in this thin interface of thickness h<sub>mix</sub> that the hot and cold gas fully mixes into gas of temperature T<sub>mix</sub> and significantly cools. Across the whole TML of h<sub>turb</sub> > h<sub>mix</sub>, the temperature smoothly transits from T<sub>hot</sub> to T<sub>cool</sub>, similar to the case within a thermal conduction front. The velocity dispersion within the cooling interface σ<sub>mix</sub> is smaller than the global value σ<sub>turb</sub> within the whole TML (Figure 6 of T21).

The well-mixed gas in the cooling interface condenses into numerous H I cloudlets, of which each has a mass m<sub>cloud</sub> = 1.4m<sub>p</sub>C<sub>H1</sub>, where C<sub>H1</sub> denotes the number of neutral hydrogen atoms in one cloudlet. We treat the cloudlets as individual particles, each having its own internal thermal temperature and velocity turbulence, while they together form a mist with a cloud number density n<sub>cloud</sub> and rms velocity dispersion σ<sub>cloud</sub>. As the ensemble measure within a beam-sized aperture, σ<sub>cloud</sub> is a robust proxy for the turbulence of the mixed gas before condensation, σ<sub>mix</sub> (M. Gaspari et al. 2018).

In this regime, the observed H I is the condensed cloudlets, and the measurements can be written as

$$N_{\text{H1}} = 2n_{\text{H1}}h_{\text{mix}} = 2n_{\text{cloud}}C_{\text{H1}}h_{\text{mix}} < 2n_{\text{cloud}}C_{\text{H1}}h_{\text{turb}}, \quad (6)$$

$$\sigma_v = \sigma_{\text{mix}} < \sigma_{\text{turb}}. \quad (7)$$

The pressure in the cooling layer is dominated by the cloudlets, which fill roughly half the volume and account for most of the

mass. Similar to the Equation (5), the pressure equilibrium can be written as

$$p_{\text{CGM}} = \bar{\rho}_{\text{cloud}} \sigma_{\text{mix}}^2, \quad (8)$$

where the cloudlet density averaged over the whole cooling interface is $\bar{\rho}_{\text{cloud}} = 1.4m_{\text{p}}n_{\text{cloud}}C_{\text{H1}}$.

## 5.2. Expectations for the Observations

Considering the scenario introduced above, we expect the following observational results in a statistical sense.

1. A relation of N<sub>H1</sub> ~ σ<sub>v</sub><sup>2</sup> with scatter. If most of the TMLs are in the strong-cooling regime, we expect that N<sub>H1</sub> ~ h<sub>turb</sub> ~ σ<sub>turb</sub><sup>2</sup> = σ<sub>v</sub><sup>2</sup> (Equations (2)-(4)). Such a relation has been found in Section 4.3 (Figure 5). If some of the TMLs are in the weak-cooling regime, data points would deviate from the relation, as h<sub>mix</sub> < h<sub>turb</sub> and σ<sub>mix</sub> < σ<sub>turb</sub>.

2. The thickness of the H I layer h<sub>H1</sub>. The pressure equilibrium relations of Equations (5) and (8) indicate that an effective H I thickness h<sub>H1</sub> can be derived as the h in Section 4.4, but the interpretation is different.

In the strong-cooling regime, h<sub>H1</sub> is half the TML thickness h<sub>turb</sub> (Equation (3)). Considering that the H I structures have two sides of shearing interface, we could estimate that h<sub>turb</sub> ~ 2h<sub>H1</sub> = h, which is roughly 47–290 pc (25th and 75th percentiles out of 2 R<sub>25</sub>), consistent with the typical values of ~100 pc in simulations (e.g., J19; T21; Y. Yang & S. Ji 2023).

In the weak-cooling regime, h<sub>H1</sub> is roughly the cooling interface thickness h<sub>mix</sub> ~ h/2 (Equation (6)). In this case, the H I is highly clumpy, and is thus more capable of surviving photoionization.

3. The possibility of estimating the CGM gas entrainment velocity (v<sub>ent</sub>). The corresponding cooling rate per surface area can be written as dM<sub>cool</sub>/dA = ρ<sub>CGM</sub>v<sub>ent</sub> (more calculations in Section 5.3). The value of v<sub>ent</sub> is regulated by the turbulence diffusion and cooling rate within the TML, related to the two regimes discussed above. By the conservation of angular momentum, the gas that cooled from the rather static CGM would flow in (Figure 8(b)) at the rate of M<sub>in</sub>, which we estimate in Section 5.4.

## 5.3. The Global Cooling Rate

Under the turbulent-mixing gas cooling scenario, the hot gas would mix and cool in the TML, and then get entrained into the cool phase. Below, we estimate this global cooling/entrainment rate M<sub>cool</sub> from observations, using the expression of dM<sub>cool</sub>/dA = ρ<sub>CGM</sub>v<sub>ent</sub> in Section 5.2.

We calculate M<sub>cool</sub> within the region that has moment 2 values (Figure 1(c)), and exclude the inner disk with r < R<sub>25</sub> and the highly dispersed region with σ<sub>v</sub> > 10<sup>1.9</sup> km s<sup>-1</sup> ~ 80 km s<sup>-1</sup>, where feedback and fountain possibly play a more important role. We take the velocity dispersion of diffuse H I (σ<sub>v</sub>) as σ<sub>turb</sub>, and v<sub>tot</sub> (Section 3.2) as the shearing velocity v<sub>s</sub>. At each pixel, an additional factor of 2/cos i is multiplied to account for the projection of area and the two sides of H I structure. We estimate the cooling time t<sub>cool</sub> as (5/2)(k<sub>B</sub>T<sub>mix</sub>)<sup>2</sup>/pΛ, where p and Λ are the partial pressure of the hydrogen and cooling function at T<sub>mix</sub> = 2 × 10<sup>5</sup> K, with a metallicity of 0.3 Z<sub>☉</sub>, which is typical for the CGM of an

<footnote n="22"><sup>22</sup> T21 estimated t<sub>turb</sub> as the eddy turnover time (L<sub>eddy</sub>/u'), where L<sub>eddy</sub> is the outer scale of eddies, and u' ~ √3 σ<sub>v</sub> is the 3D turbulent velocity.</footnote>
<table><caption>**Table 1**
Global Turbulent-mixing Gas Cooling Rates</caption><thead><tr><th>Global M<sub>cool</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>v<sub>ent</sub> ~ σ<sup>2</sup><sub>turb</sub>/v<sub>s</sub></th><th>J19 Calibration</th><th>T21 Calibration</th></tr></thead><tbody><tr><td>Modeled v<sub>s</sub></td><td>2.1</td><td>1.8</td><td>3.9</td></tr><tr><td>v<sub>s</sub> = 120 km s<sup>-1</sup></td><td>1.8</td><td>1.9</td><td>3.8</td></tr><tr><td>v<sub>s</sub> = 150 km s<sup>-1</sup></td><td>1.4</td><td>1.6</td><td>4.3</td></tr></tbody></table>

**Note.** The cooling rates are summed over all diffuse-only region, excluding the stellar disk ($r < R_{25}$) and the highly dispersed region. Modeled inclination $i$ and shearing velocity $v_s$ are reconstructed using tilted ring modeling (Section 3.2 and Appendix D). Three models are used for the estimation of hot gas entrainment velocity $v_{ent}$. Additionally, the values assuming a constant $v_s = 120$ or $150$ km s$^{-1}$ are reported in the second row. For more information, refer to the text in Section 5.3.

MW-type galaxy such as M51 (e.g., M. J. Miller & J. N. Bregman 2015; N. Lehner 2017; J. X. Prochaska et al. 2017).

We conduct three sets of estimations of $v_{ent}$ and $\dot{M}_{cool}$, using both an idealized model and the latest calibrations from hydrodynamic simulations. Their order-of-magnitude consistency lends support for the estimation, while their differences represent the model uncertainties.

1. *Small-perturbation approximation.* From the equation of continuity, J19 obtained the approximate relation $v_{ent}/v_s \sim h_{turb}/L$. Combining Equation (2), $v_{ent}$ can be estimated as $v_{ent} \sim \sigma_{turb}^2/v_s$. The resultant global cooling rate $\dot{M}_{cool}$ of $2.1 M_{☉} \text{yr}^{-1}$ could cumulatively build a column density of $8.3 \times 10^{18} \text{cm}^{-2}$ in $\sim 300$ Myr averaged over the region included in $\dot{M}_{cool}$ calculation.

2. *The J19 hydrodynamic calibration.* J19 calibrated a scaling relation of $v_{ent}(v_s, t_{cool})$ (their Equation (24)) using 3D hydrodynamical simulations. T21 pointed out that the J19 calibration was carried out for the weak-cooling regime, and should be scaled down by $Da^{-1/4}$ when describing the strong-cooling regime. We calculate $Da$ per pixel as $(h_{turb}/\sigma_{turb}t_{cool})^{4/3}$ (T21), and obtain a global cooling rate of $\sim 1.8 M_{☉} \text{yr}^{-1}$.

3. *The T21 hydrodynamic calibration.* With Equations (53) and (54) of T21, we derived a global cooling rate of $3.9 M_{☉} \text{yr}^{-1}$. In this case, the eddy scale $L_{eddy}$ and the turbulence velocity $u'$ at this scale are required, but are difficult to constrain from observation. Therefore, we estimate them using Equation (55) in T21, with an assumed $L_{eddy}$ of 100 pc and cool gas sound speed of $15 \text{ km s}^{-1}$, the typical values that T21 used.

We summarize the estimations of global cooling rate in Table 1, with additional values using a constant $v_s = 120 \text{ km s}^{-1}$ (roughly the median $v_{tot}$ beyond $R_{25}$) or $150 \text{ km s}^{-1}$ (similar to the CGM sound speed $c_{s,CGM}$) to circumvent the uncertainty of tilted ring model. All of these estimations are comparable to the SFR of $4.5_{-2.6}^{+1.6} M_{☉} \text{yr}^{-1}$ in M51 (A. K. Leroy et al. 2019), indicating that the turbulent-mixing gas cooling may contribute significantly to the fueling of star formation, if they can be efficiently transported inward to the star-forming disk. We discuss the possible inflow in the next section.

### 5.4. The Global Inflow Rate After Entrainment

The cooling of rotationally lagged CGM onto the H I disk reduces the specific angular momentum and induces inflow.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** The estimated inflow rates due to turbulent-mixing gas cooling at each modeled ring of radius $r_{ring}$. The curves with different line styles represent the results from three different models: the continuity model of $v_{ent} \sim \sigma_{turb}^2/v_s$ (blue solid), J19 (orange dashed), and T21 (green dashed-dotted) calibrations. The thin and thick curves correspond to the estimation using the tilted ring-fitted velocity or a constant shearing velocity of $v_s = 120 \text{ km s}^{-1}$. The value of the SFR is the horizontal gray dashed line, and the averaged value of the curves across the radial range in display is the horizontal red dashed-dotted-dotted line; the uncertainties are indicated as the shaded region.</figcaption></figure>

We assume that before condensation, the hot CGM has zero velocity on average, given that it is primarily supported by dispersion and that currently, the CGM velocity field cannot be conclusively modeled from absorption (e.g., C. L. Martin et al. 2019; D. M. French & B. P. Wakker 2020; S. Weng et al. 2023). One study reported that the hot gas within 50 kpc around the MW has a rotational velocity as high as 75% of that of the stellar disk (E. J. Hodges-Kluck et al. 2016). However, they used O VII line measurements, thus biased toward metal-rich gas and likely related to galactic fountains, rather than to the bulk motion of less-enriched low-spin CGM (H. J. Mo et al. 1998), the latter of which is relevant to the gas accretion we discuss here.

For a circularly rotating ring with the radius $r_{ring}$, the inflow rate induced by cooling at this radius would be

$$ \dot{M}_{\text{in,TML}}(r_{\text{ring}}) = 2\pi r_{\text{ring}}^2 \left\langle \frac{d\dot{M}_{\text{cool}}}{dA} \right\rangle, \quad (9) $$

with an assumed flat rotation curve. Here, $\langle d\dot{M}_{cool}/dA \rangle$ is the entrainment rate per area $\rho_{CGM}v_{ent}$ averaged over the ring.

We calculated this value between $R_{25}$ and $4R_{25}$ using the rings from tilted ring modeling, shown in Figure 9. Averaged over this radial range and over six models, the $\dot{M}_{in,TML}$ is $1.6 M_{☉} \text{yr}^{-1}$ (red dashed-dotted-dotted line, with a scatter of $1.0 M_{☉} \text{yr}^{-1}$), roughly one-third of the SFR. Therefore, the inflow rate, with some scatter, could account for a considerable part of the SFR fueling. Interestingly, based on interferometric data, the dense H I inflow rate is measured to be $\sim 0.1-0.2 \times \text{SFR for the } 10^{20} \text{cm}^{-2} \text{ gas at a typical radius of } 1.5-3 R_{25}$ (E. M. Di Teodoro & J. E. G. Peek 2021). Turbulent-mixing gas cooling could serve as an important channel of gas accretion from CGM that replenishes both the outer disk of dense H I and the inner star-forming disk.
We note that by plotting $\dot{M}_{\text{in,TML}}$ as a function of $r_{\text{ring}}$, we do not intend to claim a full description of the continuous radial inflow across the H I structures driven by turbulent-mixing gas cooling, but to infer the feasibility of such a picture from an order-of-magnitude and mass-conservation perspective. If such a gas inflow exists in a less dynamic and quasi-equilibrium system, $\dot{M}_{\text{in,TML}}$ should remain roughly constant throughout the radial range, and fluctuate around the global cooling rate and SFR. It is interesting that this prediction is roughly what we find in Figure 9, which urges the next step extension to more isolated galaxies.

## 5.5. Caveats and Other Possible Effects

During the calculation of cooling and inflow rates, we make many simplifications and approximations. We adopt the TML scaling relations calibrated from idealized TML simulations by J19 and T21, who used typical configurations of radiation, temperature, metallicity, etc. for the CGM of MW-type galaxies. These two studies did not explicitly consider the role of dust, which, however, is poorly constrained in observations and should not have a significant role beyond the stellar disk, the main region of our analyses. While the heating and cooling mechanisms are indeed complex, with possibly varying conditions in real galaxies, their overall effects are practically encapsulated into parameters such as cooling coefficients or cooling timescales. Because the dependence on these cooling parameters has been included in the calibrated scaling relations, which are carried out with CGM-like parameters as mentioned, these scaling relations remain broadly applicable for our case, without being tied to the detailed physical modeling of specific heating and cooling mechanisms.

The entrainment velocity $v_{\text{ent}}$ could contribute to the observed values of the moment 2 map, making $\sigma_{\text{turb}}$ overestimated. Meanwhile, by using the value of $\sigma_v$ as $\sigma_{\text{turb}}$, $\sigma_{\text{turb}}$ is underestimated in the weak-cooling regime (Equation (7)). These effects lead to both the over- and underestimation of the cooling rate. The eddy scale $L_{\text{eddy}}$ is also a source of uncertainty to the estimation based on the T21 model.

The tidal forces may cause the H I global structure to expand, and slow down the radial inflow rate. However, for the observed snapshot of M51, the tidal effects on the diffuse H I is likely weak, as the current radial extension of H I is much larger than the prediction in the simulation of C. L. Dobbs et al. (2010, see Section 4.2.1).

If the H I gas is in the form of small, separate clouds, the ram pressure may slow them down, causing them to flow in (e.g., Y.-H. Huang et al. 2016) until they become interlinked. After that, ram pressure could both strip and compress the diffuse H I (J. E. Gunn & J. R. Gott 1972), causing it to lag behind dense H I as shown in the PV diagram (Figures 4(a)–(c)). Therefore, ram pressure can both assist and harm the growth and inflow of diffuse H I.

Finally, taking into account the feedback-driven outflows in the star-forming disk, the total inflow rates onto the disk should be systematically higher than the SFR under a quasi-steady assumption, so as to sustain the star formation. We thus expect the turbulent-mixing gas accretion to work together with other channels, particularly the galactic fountains within an optical disk, where we exclude for most of the time in this study. The necessity of having different gas accretion channels, as well as the separate role of each, has been the topics of recent cosmological simulations (R. J. J. Grand et al. 2019).

## 5.6. The Turbulent-mixing Cooling in the Context of Gas Accretion and Galaxy Evolution

MW-type galaxies must accrete gas; otherwise, they cannot show their current neutral gas richness after the cosmic time consumption from star formation (C. Péroux & J. C. Howk 2020; L. J. Tacconi et al. 2020). While merger and stripping of gas-rich satellites cannot make a significant contribution (G. Kauffmann et al. 2013), the cooling and accretion from the virialized CGM (the so-called hot-mode gas accretion) is the most likely origin.

There are several theoretically possible ways for the accretion to take place. The most extensively studied mechanisms include the precipitation, the cooling flow, and the fountain. In the precipitation scenario, $t_{\text{cool}} \ll t_{\text{dyn}}$ (the dynamical time), cool clouds condense out of the hot CGM due to thermal instabilities, and drop onto the galaxy due to the loss of buoyancy (e.g., G. M. Voit et al. 2017; M. Gronke & S. P. Oh 2020b). In the cooling flow scenario, $t_{\text{cool}} \sim t_{\text{dyn}}$, the hot gas flows radially inward due to the loss of pressure support, with radiative cooling balancing the compression heating; the cooling into the cold phase only happens when it gets close to the disk (e.g., Z. Hafen et al. 2022). In the fountain scenario, metal-enriched cool gas is launched by supernovae perpendicular to the disk and then rains back, while CGM cooling is boosted through the mixing with it (e.g., F. Fraternali 2017).

In the current study, we advocate an alternative scenario of disk–CGM interface cooling, where the tidal structures are counted as the extension of the H I disk, and the turbulent mixing serves as the primary driver. The key requirement for this scenario is an extended H I structure with subsonic shearing motion relative to the CGM. Since most star-forming galaxies have a large H I disk rotating at a velocity close to the virial one, this scenario can be applied to general star-forming galaxies.

The disk–CGM interface cooling has clear links to the other three scenarios. The precipitation and cooling flow are like the preparing process, bringing gas close to the disk–CGM interface for further turbulent-mixing gas cooling. The cooling at the interface between the launched clouds and hot gas is an essential part of the galactic fountain model, whose localized thermodynamic physics is very close to the disk–CGM interface cooling. The former increases interface with the CGM initially through ejecting outflows, while the latter through tidal interaction or simple mechanisms sets the H I size–mass relation.

Meanwhile, its difference from the other three scenarios is also clear. The proposed cooling process happens at the interface of an existing H I disk or structure with the CGM, without the necessity of requiring the final success of the other three CGM cooling scenarios (e.g., not all precipitated clouds survive the journey from large distances onto the disk, Z. Li et al. 2020). It is a promising way to maintain the gas richness of a galaxy that is already gas-rich, and sustain the existing active star formation.

At cosmic noon, turbulent mixing potentially enhances the cold-mode accretion efficiency (H. Aung et al. 2024). In this scenario, the shearing happens between the cold stream and hot CGM, and could be the high-redshift counterpart of our case.

# 6. Summary

Using a new MEM pipeline, we succeed in jointly deconvolving the single-dish and interferometric H I cubes
from FAST FEASTS and VLA THINGS of the extragalactic M51 system. At a spatial resolution of ~18" (~0.7 kpc at 8.0 Mpc), the LOS column density sensitivity is ~3.2 × 18 cm⁻² at 2σ confidence. With the deconvolved cube and moment maps, we obtain the following results:

1. The high $N_{H1}$ structures are consistent with the prediction of a previous hydrodynamical simulation (C. L. Dobbs et al. 2010), but the newly detected low $N_{H1}$ structures are much more extended. It is thus necessary to take the CGM-ISM interaction into consideration.

2. Roughly half (48%) of the H I mass in M51 is diffuse, and 59% of this diffuse H I is not associated with any dense H I (or is “diffuse-only”), covering 89% of the LOS sample.

3. For diffuse-only H I, the values of $\sigma_v$ and $N_{H1}$ show a power-law relation of $\sigma_v \propto N_{H1}^{0.5}$.

4. The H I gas outside of the galaxy optical radius ($R_{25}$) is found to probably survive photoionization from the UV background, but may be significantly influenced by thermal evaporation under the L. L. Cowie & C. F. McKee (1977) model.

With these observational results, we propose that in this system, turbulent-mixing gas cooling could be an important and previously overlooked hot gas cooling and accretion channel. In this scenario, the vast CGM-ISM interface is generated by tidal interaction, where the shearing between them produces TMLs and induces cooling. We make the following developments in linking the TML theories and simulations with H I observations.

1. We design ensemble measurements and estimators for the attributes of TML-related H I gas, including the thickness, velocity dispersion, and pressure. The properties of M51 diffuse H I, particularly the $\sigma_v-N_{H1}$ relation, are consistent with the scenario that they are primarily the H I that has condensed from the TML.

2. We estimate the rates of turbulent-mixing gas cooling and the subsequent gas inflow using three different recipes. These rates are around one-third of the M51 SFR. Such a cooling could accumulate an averaged H I column density of ~$10^{19}$ cm⁻² in ~300 Myr, potentially accounting for a considerable part of cold gas replenishment.

This work illustrates the possibility of probing the multiscale structures, kinematics, and processes at ISM-CGM interface through H I emission observations. The proposed turbulent-mixing gas accretion is worth investigation for more-general H I-rich galaxies. For a FEASTS sample of 35 moderately inclined H I-rich galaxies, the H I disks down to the $N_{H1} \sim 10^{18}$ cm⁻² level on average extend to 70 kpc or 0.3 times virial radius (J. Wang et al. 2025), providing a large area of shearing interface with the hot CGM. Previously, these $N_{H1} \sim 10^{18}$ cm⁻² regions could only be detected from low-ion metal absorption along selected LOSs.

Our observation and analyses on M51 favor the picture that the H I gas around this H I-rich merging system distributes as a relatively contiguous and large structure extending into the CGM. At the ISM-CGM interface, the proposed turbulent-mixing gas accretion offers a possible explanation for the cold gas enhancement observed statistically in post-merger galaxies (e.g., S. L. Ellison et al. 2015, 2018; J. Shangguan et al. 2019).

We expect systems similar to M51 can be more common in the early Universe, where mergers are more frequent and gas-rich. If so, a large fraction of ~$10^{18}$ cm⁻² gas detected by absorptions could be related to tidally deformed H I structures, particularly at the inner part of the halo, rather than being the projection of a 3D sea of distinct clouds, or forming a number of major filaments. Also, the possible turbulent-mixing gas cooling may contribute more significantly to the mass assembly of these high-redshift galaxies. More future combinations of FEASTS and archival or new interferometric data could help to further constrain this key component of gas accretion.

## Acknowledgments

We thank the anonymous referee for providing constructive and helpful comments. We thank Bi-Qing For, Gregory J. Herczeg, Jin Koda, Thijs van der Hulst, Junfeng Wang, Simon Weng, and Tobias Westmeier for useful suggestions and discussions. We asked Google Gemini (1.5 Flash, accessed on 2024 October 29) for suggestions on the grammar and expression of the Introduction.

J.W. acknowledges research grants and support from the Ministry of Science and Technology of the People's Republic of China (No. 2022YFA1602902), the National Natural Science Foundation of China (No. 12073002), and the China Manned Space Program (No. CMS-CSST-2021-B02). S.J. acknowledges support from the National Natural Science Foundation of China (Nos. 12133008, 12192220, and 12192223), the China Manned Space Program through its Space Application System, and the National Key Research and Development Program of China (No. 2023YFB3002502). N.M. acknowledges support from Israel Science Foundation (ISF) grant 3061/21 and US-Israel Binational Sciences Foundation (BSF) grant BSF 2022281. S.H. O. acknowledges support from the National Research Foundation of Korea (NRF) grant funded by the Korea government (Ministry of Science and ICT: MSIT) (No. RS-2022-00197685).

This work has used the data from FAST (https://cstr.cn/resolver?identifier=31116.02.FAST). FAST is a Chinese national mega-science facility, operated by the National Astronomical Observatories, Chinese Academy of Sciences (NAOC).

Facilities: FAST:500m, NED, VLA.

Software: Astropy 5.2.2 (Astropy Collaboration et al. 2022), BBarolo 1.7 (E. M. Di Teodoro & F. Fraternali 2015), CASA 6.5.3 (CASA Team et al. 2022), Karma 1.7.25 (R. Gooch 1996), mosmem 1.10 and other MIRIAD utilities (R. J. Sault et al. 1995), NumPy 1.24.3 (C. R. Harris et al. 2020), Python 3.10.12, pingouin 0.5.4 (R. Vallat 2018), pvextractor 0.3 (A. Ginsburg et al. 2016), regions 0.8 (L. Bradley et al. 2023), SciPy 1.10.1 (K. W. Vugrin et al. 2007; P. Virtanen et al. 2020), SoFiA 2.5.0 (T. Westmeier et al. 2021).

## Appendix A

### Improving the Continuum Subtraction of the VLA Data in the Visibility Domain

The limited velocity range of THINGS data (from 214.9 to 704.5 km s⁻¹) leads to an improper continuum subtraction at $v \ge 658.1$ km s⁻¹ (see the VLA observation of A. H. Rots et al. 1990 with a broader spectral window). It mainly influences a small region of ~13 arcmin² to the north of M51b (see the contour on the channel maps in Appendix G).
We thus added back the missing interferometric flux in that region using the data from a shallower VLA observation by A. H. Rots et al. (1990). With a beam of $34'' \times 27''$5, the A. H. Rots et al. (1990) data have an rms level of $0.69 \text{ mJy beam}^{-1}$. We clean the data of A. H. Rots et al. (1990) and THINGS<fnref n="23" /> with MIRIAD CLEAN down to their respective $1\sigma$ levels. We obtain the difference of their cleaned models within this small region, and convert the difference into visibility, which is added onto the original THINGS visibility.

We note that this relatively narrow velocity range does not introduce a significant loss of H I flux. Only 0.16% of the M51 flux detected by FEASTS is located beyond the velocity of 704.5 km s$^{-1}$, the maximum velocity observed by THINGS.

## Appendix B
### The Joint Deconvolution Pipeline

In this appendix, we provide the details of the joint deconvolution pipeline introduced in Section 3.1.

*Data preparation.* The interferometric data are the dirty cube $d_j^{\text{int}}$ described in Section 2.2. The dirty beam image $b^{\text{int}}$ is twice as large as the cube, of $3600 \times 3600$ pixels.

For the single-dish cube $d_j^{\text{sd}}$, we use the FAST data after flux cross-calibration (scaled down by a factor of 1.15) with the interferometric data as described in Paper I, and regrid them into the coordinate system of $d_j^{\text{int}}$. The World Coordinate System calibration is similar to Paper I, but is done within the central $600'' \times 600''$ box instead, which fully covers the brightest inner H I disk. The World Coordinate System reference pixel values, CRPIX1/2, of FAST cube are increased by -0.15 and -0.2. The single-dish beam $b^{\text{sd}}$ is the average of the 19-beam FAST receiver beams described in Appendix B.2 of Paper I, after necessary coordinate projections.

The mask for the MEM deconvolution is taken as the mask of the FAST cube generated using SoFiA as described in Section 2.2 of Paper I. The regridded single-dish cube $d_j^{\text{sd}}$ has an rms of $\sigma_{\text{sd}} = 0.489 \text{ mJy beam}^{-1}$ beyond the mask, lower than the one ($0.80 \text{ mJy beam}^{-1}$) reported in Paper I due to the binning during regridding. Only within this deconvolution mask will the MEM deconvolution subroutine generate the nonzero pixels of the solution $m_j$, calculate the entropy $H$, and measure the residual rms.

The cubes of $d_j^{\text{int}}$ and $d_j^{\text{sd}}$ and the deconvolution mask are then separated into individual channels. The following deconvolution procedure is then run per channel as described below. For simplicity, these single-channel images are also referred to as $d_j^{\text{int}}$ and $d_j^{\text{sd}}$.

*The procedures.* The major steps of the single-channel joint deconvolution iterations are as follows. In each of these iterations (labeled as *n*th), the MEM routine *mosmem* is invoked once with the same $d_j^{\text{int}}$ and $d_j^{\text{sd}}$, and different reference images $\hat{m}_j^{(n)}$. *Mosmem* derives the best model $m_j^{(n)}$ for the given $\hat{m}_j^{(n)}$ and threshold $\theta^{(n)}$, as described in Section 3.1.

We always set the rms convergence threshold for single-dish map as $\theta_{\text{sd}}^{(n)} = \theta_{\text{int}}^{(n)} \text{beam}_{\text{FAST}}/\text{beam}_{\text{V LA}}$,<fnref n="24" /> and both $\theta^{(n)}$'s have

a tolerance of 5%. In each iteration, the initial guess for $m_j^{(n)}$ is always a flat estimate filling the deconvolution mask with the total flux of $d_j^{\text{sd}}$.

1. *First iteration.* In the first call of MEM deconvolution, $\hat{m}_j^{(1)}$ is set to be $d_j^{\text{sd}}$, and the interferometric rms threshold $\theta_{\text{int}}^{(1)} = 1.5\sigma_{\text{int}}$.

2. *Generate the new reference image.* The output model $m_j^{(n)}$ (after primary-beam (PB) correction) and residuals of the *n*th iteration are combined to produce the new reference $\hat{m}_j^{(n+1)}$ for the next iteration. We search for residual peaks that are potential signals, and add them onto $m_j^{(n)}$ as $\hat{m}_j^{(n+1)}$, thus increasing the deconvolution quality and speeding up the convergence, which shares the same concept as the peak-searching step between iterations in traditional CLEAN (J. A. Högbom 1974; B. G. Clark 1980; D. G. Steer et al. 1984).

For the interferometric residual $r_{\text{int},j}^{(n)} = d_j^{\text{int}} - b^{\text{int}} * (\text{PB} \times m_j^{(n)})$, we select out the pixels where both $r_{\text{int},j}^{(n)}$ and PB $\times m_j^{(n)}$ are high enough relative to the interferometric noise level $\sigma_{\text{int}}$.<fnref n="25" /> We require the product of these two S/N values to be higher than $3^2$, so as to include pixels where at least one S/N is very significant. A hyperbolic tangent (tanh) function is used to generate a smooth weight map and thus to reduce the clumpiness of the new reference map.

In a similar way, the pixels of single-dish residual $r_{\text{sd},j}^{(n)} = d_j^{\text{sd}} - b^{\text{sd}} * m_j^{(n)}$ are selected out with a threshold of $3^2$. We use two sets of criteria, (a) comparing $r_{\text{sd},j}^{(n)}$ and $m_j^{(n)}$ with $\sigma_{\text{sd}}$ and $\sqrt{\sigma_{\text{int}}\sigma_{\text{sd}}}$ (converted to the same beam area) and (b) comparing $r_{\text{sd},j}^{(n)}$ and the smoothed $r_{\text{int},j}^{(n)}$ with $\sigma_{\text{sd}}$ and $\sigma_{\text{int}}/\sqrt{1+3^2}$, where the Gaussian smooth kernel has three times longer major and minor axes than the clean beam. The motivation of criterion b is to find the region where both data sets have significant positive residuals. We compare the weights from these two criteria, and only use the larger one.

When added onto $m_j^{(n)}$, the weighted $r_{\text{int},j}^{(n)}$ is PB corrected, and is scaled down to account for the difference in the dirty and clean beam areas (S. Jörsäter & G. A. van Moorsel 1995). To avoid repetitively adding fluxes, we subtract the weighted $r_{\text{int},j}^{(n)}$ from the weighted $r_{\text{sd},j}^{(n)}$ values after convolving it with the single-dish beam $b^{\text{sd}}$, and replace the negative pixels in the resultant $r_{\text{sd},j}^{(n)}$ with 0.

3. *The (n + 1)th iteration.* We calculate the interferometric residual rms $\sigma_{\text{int,out}}^{(n)}$ outside the joint deconvolution mask of the previous residual $r_{\text{int},j}^{(n)}$. If $\sigma_{\text{int,out}}^{(n)} < \theta_{\text{int}}^{(n)}$, we update the $\theta_{\text{int}}^{(n+1)}$ to be $\sigma_{\text{int,out}}^{(n)}$; otherwise, $\theta_{\text{int}}^{(n+1)}$ is kept as $\theta_{\text{int}}^{(n)}$. Then, the MEM routine is called with the new interferometric rms threshold $\theta_{\text{int}}^{(n+1)}$ and the new reference image $\hat{m}_j^{(n+1)}$ from step 2.

4. *Stop criterion.* The iteration is repeated for at least eight times. A channel is considered fully deconvolved, and thus the iteration stops, if the interferometric residual rms within the deconvolution mask is lower than the one out

<footnote n="23">²³ When doing the correction, we taper the THINGS cube so that it has the same clean beam as the A. H. Rots et al. (1990) data. The rms level becomes 0.57 mJy beam$^{-1}$ after this tapering.</footnote>

<footnote n="24">²⁴ Thus, both $\theta^{(n)}$ have the same value when they are converted to Jy pixel$^{-1}$. In other words, the interferometric and single-dish images share the same physical rms criterion.</footnote>

<footnote n="25">²⁵ Here, the $m_j^{(n)}$ is compared with $\sigma_{\text{int}}$ instead of $\sqrt{\sigma_{\text{int}}\sigma_{\text{sd}}}$, because the goal is to select out the high S/N region relative to the interferometric residual that will be added back.</footnote>
of the mask, $\sigma_{\text{int,out}}^{(n)}$. We choose this criterion because the major bottleneck of joint deconvolution is the interferometric noise. Otherwise, steps 2 and 3 would be repeated.

5. *The solution for this channel.* The stop criterion of step 4 only considers the interferometric residual. To quantify the overall quality of the output models, we use the average relative residual level $\rho$, calculated as $\sqrt{\rho_{\text{int}} \rho_{\text{sd}}}$, where $\rho_{\text{int}}$ (or $\rho_{\text{sd}}$) is the interferometric (single-dish) residual rms within the deconvolution mask and has been normalized by $\sigma_{\text{int}}$ ($\sigma_{\text{sd}}$). For each channel, the output of the iteration with the lowest $\rho$ is automatically selected as the final solution. We note that two out of 88 channels have a very small deconvolution mask, and thus do not meet the stop criteria within 20 iterations. For them, the final iteration model is used. We also manually select the best iteration for two other channels for a better overall residual map.

*Final data products.* When all channels have finished the above iterations, we convolve each of the output models ($m_j$) with the restoring beam that is chosen as interferometric clean beam ($b^{\text{clean}}$), and combine all of them into one single cube. The medians and 90% distribution intervals of the relative residual levels $\rho$, $\rho_{\text{int}}$, and $\rho_{\text{sd}}$ are $1.54_{-0.43}^{+0.97}$, $1.11_{-0.12}^{+0.04}$ and $2.1_{-1.0}^{+3.4}$, respectively.

# Appendix C

## A Justification of the Joint Deconvolution Procedure

One typical concern about the deconvolution is over-deconvolution. A reasonable expectation of deconvolution is that, similar to the traditional clean used in deconvolving dirty images, it can only sharpen the features well-resolved at the existing resolution. Indeed, one important role of the information entropy $H$ is to suppress overdeconvolution, particularly when the FAST image, which has a large but relatively simple beam, is used as the initial reference image $\hat{m}_j^{(1)}$ (i.e., the zero-point of $H$) to calculate the $H$. Hence, the model $m_j^{(1)}$ by construction should be relatively smooth, at the expense that the localized peaks may be overlooked (i.e., insufficient deconvolution). At early iterations, such a disadvantage of underdeconvolution is manifested as the clear and sharp patterns remaining in the residual images $r_{\text{int,sd},j}^{(n)}$.

To compensate for the insensitivity at small scales, one key feature of our pipeline is to iteratively decide whether the connected positive pixels in $r_{\text{int},j}^{(n)}$ contain statistically significant signals, and thus whether they should be added back to make a new reference image $\hat{m}_j^{(n+1)}$. This step is very similar to the step between multiscale clean iterations (T. J. Cornwell 2008). The difference (advantage) here is that the FAST signals locally help to verify the THINGS residual patterns with only a moderately S/N, which would not be considered statistically significant just by themselves. By construction, the THINGS image contributes the most information where it has sufficient S/N, the FAST image takes over where the THINGS image only contains noise, and these two images determine the model together in the intermediate region. Putting together, the whole procedure is designed to achieve a model that is incompletely but sufficiently deconvolved at the given resolution of raw data.

In this appendix, we discuss the accuracy of our joint deconvolution pipeline from several different points of view. Appendix C.1 gives a qualitative discussion on the PB attenuation. A mock test quantifies the accuracy of $N_{\text{H1}}$ and

$\sigma_v$ at the effective resolution of restoring beam ($\sim 18''$) in Appendix C.2. We estimate the instrumental resolution directly from the raw data noise in Appendix C.3. Finally, we compare the jointly deconvolved moment 0 map to the data before joint deconvolution and to the result of linear combination in Appendix C.4.

### C.1. Influence of the Radially Varying Uncertainties

One complicating factor is the PB attenuation effect on the THINGS image, which increases the noise at large radius and thus obscures the THINGS signal. However, by construction, the THINGS noise only propagates to the regions where THINGS signal is significant, i.e., where PB attenuation is not strong. There, the THINGS data actually lower the systematic uncertainty from the smearing of the FAST beam. Furthermore, for M51, all of the pixels with $N_{\text{H1}} > \sim 3.2 \times 18 \text{ cm}^{-2}$ (the $2\sigma$ detection limit from Appendix C.2) have a PB attenuation factor $>0.5$, indicating that the role of PB is very limited.

However, a following major worry is that the output image (model convolved with the THINGS clean beam, $b^{\text{clean}} * m_j$) has increasingly lower resolution at larger radius. Within the scope of this study, the joint deconvolution result is acceptable, as long as the statistical trends and relations of LOS $N_{\text{H1}}$ and $\sigma_v$ are relatively robust against these uncertainties. Here, we support this argument with metal experiment, and mock tests on the accuracy of the recovered $N_{\text{H1}}$ and $\sigma_v$ is given in Appendix C.2.

For the value of $N_{\text{H1}}$, Paper II demonstrated that the simple technique of linear combination (S. Stanimirovic 2002) already works surprisingly well in statistically recovering its true value at a large radius, or the low $N_{\text{H1}}$ values far below the interferometry detection limit. The stochastic and systematic uncertainties are $<0.1$ dex. The key is that by incorporating the interferometric data, we remove the additional light in the FAST-only region that is scattered from the high $N_{\text{H1}}$ regions. The joint deconvolution procedure effectively improves the removal of scattered light, and should achieve an even more reliable measurement of $N_{\text{H1}}$.

As for the kinematical measurements, at a large radius of $r$, the velocity smearing is insignificant when the velocity field $v$ is relatively smooth. To the first order, the smearing effects can be estimated as $\delta\sigma_v = v \sin i \sin\phi \delta\phi$, where $\phi$ is the azimuthal angle starting from the major axis of a tilted ring, and $\delta\phi = 0.5\text{FWHM}/r\sqrt{\sin^2\phi + \cos^2 i \cos^2\phi}$ is the azimuthal angle spanned relative to the center by a resolution element of half the FAST beam. Assuming that $v = 100 \text{ km s}^{-1}$ (roughly the median velocity of our fitted $v_{\text{tot}}$ profile), $i = 55^\circ$ at $r = 15'$ (35 kpc or $3R_{25}$ for M51) and the azimuthally averaged velocity smearing is $\langle\delta\sigma_v\rangle = 0.5\text{FWHM}/r \times 2\text{vi}/\pi = 6.6 \text{ km s}^{-1}$, and the maximum value is $\delta\sigma_{v,\text{max}} = v \sin i \times 0.5\text{FWHM}/r = 8.8 \text{ km s}^{-1}$ at $\phi = 90^\circ$, both of which are lower than most of the measured $\sigma_v$ at that radius (see Figure 6). These deductions support our statistical analyses of $N_{\text{H1}}$ and $\sigma_v$ based on the joint-devolution output.

### C.2. The Mock Test on the Recovery of Column Densities and Dispersions

We largely follow the procedure in Paper I and Paper II to generate the mock cube of an H1 disk with the size of $R_{\text{HI}} = 6.5$ and the inclination of $45^\circ$. The angular size of this mock H1 disk is larger than M51 ($R_{\text{HI}} = 6.0$, Paper I), and thus
the PB attenuation more significantly impairs the recovery of H I structures in our mock tests than in the real case.

The mock disk has a radial distribution and power spectral slope typical for real galaxies. However, those two works did not consider the velocity axis, and here, we improve them by generating the smooth part of the disk with BBarolo. The rotation curve is set to gradually increase from 200 km s⁻¹ within r = 120'' and then flattens, and the velocity dispersion drops from 35 to 5 km s⁻¹ across the radial range. No warping, P.A. changing, nor radial velocity is included. Following Paper I, we conduct the mock observations with the settings close to those of FEASTS and THINGS. The mock VLA visibility is generated using casa.simobserve. We jointly deconvolve the mock FAST cube and the VLA mock dirty cube as for M51 in this paper.

We compare the measurement of the mock FAST cube and jointly deconvolved results with that of the input model in Figure 10(a). Statistically, both the N_H1 and σ_v are recovered by our joint deconvolution pipeline remarkably well, with the 3σ-clipped median deviations of 0.01 dex. There is no clear correlation between the deviations of two parameters. The clipped scatters are both of ~0.1 dex. In contrast, the values of N_H1 and σ_v measured from the mock FAST cube are both overestimated with larger deviations and scatters.

Our mock test shows that when we calculate the moment 2 map of the joint deconvolution result, a cube mask is necessary, otherwise the σ_v would be systematically overestimated by 0.09 dex due to the noise in the spectral wings. The cube mask for the moment 2 map is generated using SoFiA 2. We add onto the cube Gaussian noise of √σ_vσsd at the resolution of the restoring beam, before inputting the noise cube to SoFiA. We use a SoFiA smoothing kernel that is 30'' wide spatially and three channels wide (15.5 km s⁻¹) spectrally. The detection and reliability thresholds are 2σ and 0.99, respectively. Only the LOSs with moment 2 values after we apply the SoFiA mask are counted in Figure 10(a).

For the statistical accuracy of N_H1, we find that it is better not to apply the cube mask when generating the moment 1 map of the joint deconvolution result. The cube mask would remove the fluxes in the spectral wing, introducing an underestimation of N_H1 of 0.02 dex (not shown in Figure 10(a)). Figure 10(b) shows that when the jointly deconvolved N_H1,JD is higher than ~3.2 × 18 cm⁻², the logarithmic deviation has a scatter smaller than 1/(2 ln 10) = 0.22 dex, i.e., a confidence level of 2σ for LOS N_H1 measurements.

Therefore, in this study, the M51 moment 2 map of the jointly deconvolved cube is generated with a SoFiA cube mask, while no mask is applied when we generate the moment 0 and moment 1 maps.

### C.3. The Resolution at Different Column Density Detection Levels

In this appendix, we estimate the best instrumental resolution that our joint deconvolution pipeline could recover at a given N_H1 detection level. While for a single-dish-only or interferometry-only data set, its beam shape and residual rms level determine the resolution and detection level, our joint deconvolution pipeline takes two different beams and produces two different residual maps. Thus, its relation between resolution and final noise level is not self-evident.

To simplify the question, we first omit the details of the pipeline, thereby deciding the limit set by the data themselves.

<figure><img src="image_13.png" /><figcaption>Figure 10. The accuracy of the column density (N_H1) and LOS velocity dispersion (σ_v) reconstructed by our joint deconvolution pipeline, tested on the mock H I cube (see text of Appendix C.2). (a) The logarithmic deviations from the mock input of the mock FAST cube (blue dashed lines) and jointly deconvolved results. The σ_v values of orange solid and green dotted contours are measured with and without SoFiA mask, respectively. The 2D-distribution contours enclose 90%, 75%, 50%, and 25% of the pixels. The 3σ-clipped medians of the 1D distributions are indicated as solid lines (and the crosses on the 2D map). The corresponding 3σ-clipped standard deviation is also reported in the legend. (b) The logarithmic deviation of jointly deconvolved N_H1 plotted against the jointly deconvolved value (N_H1,JD). The solid line is the 3σ-clipped median value, and shaded area indicates the clipped standard deviation. The dotted vertical line indicates the LOS N_H1,JD threshold of 2σ confidence, ~3.2 × 18 cm⁻².</figcaption></figure>

An ideal combination of interferometric and single-dish data is expected to maximize the S/N at each spatial frequency k in the Fourier u-v space. We model such a combination as weighted average in the u-v space, using the same weighting scheme as in J. Koda et al. (2019, Section 4.1). The averaged dirty beam and noise spectra are

$$ \tilde{b}_{\text{avg}} = \alpha \times (|\tilde{b}_{\text{int}}/\tilde{\sigma}_{\text{int}}|^2 + |\tilde{b}_{\text{sd}}/\tilde{\sigma}_{\text{sd}}|^2), \quad (C1) $$

$$ |\tilde{\sigma}_{\text{avg}}| = \alpha \times \sqrt{|\tilde{b}_{\text{int}}/\tilde{\sigma}_{\text{int}}|^2 + |\tilde{b}_{\text{sd}}/\tilde{\sigma}_{\text{sd}}|^2}, \quad (C2) $$

where the constant factor α is chosen to normalize the peak of b_{avg}, the averaged dirty beam, as 1. Here, we use the symbols b̃ and ˜σ to denote the u-v spectrum of beam and noise, respectively, and add subscripts to differentiate the interferometric (int), single-dish (sd), and averaged (avg) ones.

We transform b̃_{avg} and ˜σ_{avg} back to the image domain, and taper them with Gaussian kernels of varying sizes. From
<figure><img src="image_14.png" /></figure>

**Figure 11.** The $N_{H1}$ detection levels allowed by different resolutions for the combination of FEASTS and THINGS observations, modeled as weighted average in the Fourier $u-v$ space. The thick green line indicates the $2\sigma$ detection limit if the S/N at each spatial frequency is optimized, while the green shaded region gives the whole possible range. The two green dashed-dotted lines correspond to the cases where THINGS data are attenuated by PB. For the $2\sigma N_{H1}$ limit determined from mock test ($\sim 3.2 \times 18 \text{ cm}^{-2}$, horizontal gray dotted line, Appendix C.2), the vertical gray dotted line gives the best combined instrumental resolution, $\sim 39$ ($\sim 1.5$ kpc). For THINGS and FEASTS, the solid and dashed lines represent $2\sigma$ and $3\sigma$ limits, respectively. For reference, the resolution and detection limits of original FEASTS (blue) and THINGS (orange) data are the two crosses, and the solid (dashed) lines in corresponding color are their relations between the detection level and resolution.

FWHMs of tapered $b_{avg}$'s and rms levels of tapered $\sigma_{avg}$'s, we obtain the lowest $2\sigma N_{H1}$ detection levels $N_{H1,lim,avg}$ allowed by each resolution, shown in Figure 11 as the green thick line. The $2\sigma$ detection limits of FEASTS (THINGS) at tapered resolutions are the blue (orange) solid lines, which we denote as $N_{H1,lim,sd}$ and $N_{H1,lim,int}$, and the $3\sigma$ limits are the dashed lines. The original resolution and limits are plotted as two crosses.

In the above analyses, we have used a local-optimization model assuming that the S/N at each spatial frequency *k* is maximized, representing the best case of data combination. On the other hand, if only a global optimization is applied across all *k*, the expected detection level would be $1/\sqrt{N_{H1,lim,int}^2 + N_{H1,lim,sd}^2}$, representing the loosest case. These two cases provide a rather narrow constraint for the possible output $N_{H1,lim}$ varying as a function of resolution, plotted as the green shaded region in Figure 11. For the $2\sigma N_{H1}$ limit of $\sim 3.2 \times 18 \text{ cm}^{-2}$ (gray dotted line) obtained in Appendix C.2, the instrumental resolution range is $\sim 39''-74''$ ($\sim 1.5-2.9$ kpc). We take their average value ($\sim 57''$) as the final estimation for the instrumental resolution.

Although our pipeline conducts joint deconvolution in the image domain, the actual resolution should be closer to the best case ($\sim 39$) than to the loosest case due to the following reasons: (1) Between iterations, we correct the reference image with not only the single-dish residual map $r_{sd,j}$ and interferometric residual map $r_{int,j}$, but also a smoothed $r_{int,j}$. This procedure uses information at three spatial frequencies, effectively improving the sampling of the $u-v$ space. (2) We visually inspect the interferometric residual maps to check the deconvolution convergence, thereby confirming that no

<figure><img src="image_15.png" /></figure>

**Figure 12.** The cumulative distributions of PB attenuation factors. The red solid line is the distribution of pixels within 0.2 dex above the $2\sigma$ detection level of $\sim 3.2 \times 18 \text{ cm}^{-2}$, and the purple dashed line is that of all pixels above the detection level.

obvious moderate- to large-scale features are left in the residual.

We have so far ignored the PB attenuation, which decreases the S/N of interferometric data with a large distance to the image center, but is actually not significant for M51. In Figure 12, we plot the cumulative distributions of LOS with $\log N_{H1}$ between 18.5 and 18.7 (or above 18.5) as a function of the PB attenuation factor, and more than 84% of the pixels have a PB attenuation factor above 0.6 (or 0.8). The corresponding relations between the resolution and detection level with the interferometric noise scaled up by $0.6^{-1}$ (or $0.8^{-1}$) are the green dashed-dotted lines in Figure 11, which do not deviate significantly from the line without PB attenuation.

Meanwhile, some characteristics of the real galaxy data may assist in increasing the actual resolution. The ISM structures usually have a smooth transition between different scales and different densities, as indicated by the power-law power spectra (B. G. Elmegreen & J. Scalo 2004). MEM requires a maximum entropy relative to the reference image, thus helping to give the smoothest solution that at the same time has the lowest residual rms levels. The mock test in Appendix C.2 reveals that effectively, a $2\sigma N_{H1,lim}$ of $\sim 3.2 \times 18 \text{ cm}^{-2}$ is achieved at a high resolution of $\sim 18''$, better than the S/N-based instrumental resolution. A likely reason is that the intrinsic spatial scale for the $\sim 3.2 \times 18 \text{ cm}^{-2}$ HI gas is much larger than the instrumental resolution. This is the same underlying principle for the method of linearly combining interferometric and single-dish data in the Fourier space (S. Stanimirovic 2002; Paper II): because the large-scale HI missed by interferometers is so well-resolved by the single-dish telescope, an image at the single-dish resolution has no essential difference from the one at interferometric resolution. Therefore, adding the missing HI back to the interferometric-only data results in a combined cube with the effective resolution of the interferometric clean beam. Compared to the linear combination method, the joint deconvolution ensures a smoother transition in phases and scales between two types of data.

In summary, we estimate that at an $N_{H1}$ level of $\sim 3.2 \times 18 \text{ cm}^{-2}$, to reach a $2\sigma$ confidence for a single LOS, the individual resolution element should be $\sim 57''$ (2.2 kpc). The intrinsically large scale of low $N_{H1}$ gas helps to reach a
<figure><img src="image_16.png" /><figcaption>**Figure 13.** A comparison of the moment 0 maps of different data and/or methods under the same color coding. (a) From the new joint deconvolution method. (b) The recleaned THINGS moment 0 map. (c) FEASTS single-dish moment 0 map. (d) The moment 0 map linearly combining FEASTS and THINGS data, using the procedure of Paper II. In the lower left corner of panel (a), the open white circle displays the diameter of the instrumental resolution (~57") at $N_{H1} \sim 3.2 \times 18 \text{ cm}^{-2}$, and the white ellipse within is the restoring beam (17"84 × 15"40). The actual resolution transits between these two limits, reaching the best one when $N_{H1} > 1.6 \times 10^{19} \text{ cm}^{-2}$. The corresponding beam of panels (b) and (c) is given at the lower left corner as a white ellipse.</figcaption></figure>

higher resolution, e.g., ~18" as used in the main part of this paper, as supported by our mock tests in Appendix C.2 where the median values show good accuracy. Other factors may also influence the resolution, such as the pointing and calibration uncertainties, which we do not consider here.

### C.4. Comparison to the Data Before Joint Deconvolution and to the Result of Linear Combination

The jointly deconvolved moment 0 map (Figure 13(a)) contains most of the large-scale structure detected by FAST (Figure 13(c)) with the resolution comparable to the interferometric THINGS data (Figure 13(b)). Compared with the result of linearly combining the FEASTS data and the recleaned THINGS model using the procedure of Paper II (Figure 13(d)), the new joint deconvolution method provides a more natural transition from high to low $N_{H1}$ regions and reveals more details at the outskirts.

## Appendix D: Deprojected Velocity and Inclination Maps

The moment 1 map and the residual of the final model are shown in Figures 14(a) and (b), and the profiles of P.A. and i are given in panels (c) and (d). Except for the northern region where M51b encounters M51a, the velocity field is rather well recovered.

By tilted ring modeling, we aim to obtain the values of H I velocity without projection effects. Such reconstructed velocity maps are shown in Figures 14(f)–(h). For $v_{rot}$ and $v_{rad}$, we first smoothed the fitted radial profile with a median filter, the kernel size being 50" (five annuli). Then, we sampled on each annulus with a separation of the annulus width (~10"), linearly interpolated the samples onto a grid, and finally smoothed the grid with a kernel of FWHM = 30". By doing this, the values at the overlapping annuli are averaged out. The total velocity is calculated as $v_{tot} = \sqrt{v_{rot}^2 + v_{rad}^2}$. A similar map of inclination i is in panel (e).
<figure><img src="image_17.png" /><figcaption>**Figure 14.** The full results of tilted ring modeling using BBarolo. (a) The moment 1 map of the best-fitting model. The gray dashed line traces the major axis ends of the rings. (b) The residual between the data and model moment 1 map. (c) The profile of the position angle of the receding side, starting from the north and increasing counterclockwise. (d) The profile of inclination angle. (e) The reconstructed inclination map. ((f)–(h)) The reconstructed maps of radial ($v_{\text{rad}}$), rotational ($v_{\text{rot}}$), and total ($v_{\text{tot}}$) velocities.</figcaption></figure>

In the literature, velocity field fitting for M51a was typically limited within the relatively unperturbed part within ~300" (e.g., D. Colombo et al. 2014; S. Oikawa & Y. Sofue 2014; F. Sylos Labini et al. 2019, 2023). They generally reported a P. A. of ~170°, consistent with ours, but tended to estimate a lower *i* of 20°–35° in this radial range, possibly limited by the boundary at ~300" where *i* changes abruptly. The maximum rotation velocity in their fittings typically reaches 150 km s⁻¹, while in our fitting the value is ~125 km s⁻¹ within 300". Our fitting is not contradictory to these previous results considering the different *i* values.

## Appendix E
### Spectral Component Fitting Procedures

In the following, we provide the implementation of decomposing the diffuse- and dense H I components mentioned in Section 3.3.

1. The spectrum $I_v$ should have at least four channels brighter than 0.5 times the interferometric noise level of $\sigma_{\text{int}}$, and the peak intensity $I_{\text{max}}$ needs to be higher than $2\sigma_{\text{int}}$. Otherwise, this LOS is considered to have only diffuse H I.

2. A channel mask (mask 1) selecting those brighter than $1.5\sigma_{\text{int}}$ are generated.<fnref n="26" /> If all the selected channels in mask 1 are connected, we decide there is only one dense Gaussian component in this LOS, or else double-Gaussian fitting will be conducted.

(a) For single-Gaussian fitting, the moment 0, 1, and 2 values are used to generate the initial guess.

(b) For double-Gaussian fitting, the two regions with highest intensities are used to determine the initial

<footnote n="26">²⁶ This channel mask and the $0.5\sigma_{\text{int}}$ one below, i.e., masks 1 and 2, have both had one-channel binary closing and dilation successively applied to them, so as to link and extend adjacent features.</footnote>
guess for each of the two Gaussian models, respectively.

3. We use SciPy `optimize.curve_fit` (K. W. Vugrin et al. 2007) to get the least-squares optimized models for the region within mask 2 (brighter than $0.5\sigma_{\text{int}}$). The model function is one or two positive Gaussian function(s) of dense H I, plus a positive constant baseline. The data cubes for the dense and diffuse components are generated accordingly.

All of the fittings converge within 170 iterations. On the whole, the LOS line structures are relatively simple. Among the LOS sample, only 11% of LOSs require a decomposition, 7% of which (0.8% of the LOS sample) contain two Gaussian components.

## Appendix F

### Averaged Column Density Profile of the FEASTS Sample

We compare the H I surface density ($\Sigma_{H\text{I}}$) profiles of M51 and nine other galaxies in Figure 15 (data from Paper II). M51 has a distinct almost flat profile between 2 and $3R_{H\text{I}}$ at $\Sigma_{H\text{I}} \approx 0.2 M_{\odot} pc^{-2}$, and the slope at larger radii is also shallower than other profiles. Here, the H I radius $R_{H\text{I}}$ is defined at an H I surface density of $1 M_{\odot} pc^{-2}$.

<figure><img src="image_18.png" /></figure>

**Figure 15.** Comparing the H I surface density ($\Sigma_{H\text{I}}$) profiles of M51 (orange thick line) and nine other galaxies (gray lines), normalized by H I radius $R_{H\text{I}}$ of $1 M_{\odot} pc^{-2}$. These galaxies have all been observed by both THINGS and FEASTS, and the profiles are taken from Paper II, where THINGS and FEASTS data are linearly combined. The profile of M51 is almost flat between 1.5 and 2.5 $R_{H\text{I}}$ with $\Sigma_{H\text{I}} \approx 0.3 M_{\odot} pc^{-2}$, different from most of the other galaxies. The median H I profile from compiled interferometric observations are given as the blue dashed line for reference (J. Wang et al. 2020).
Appendix G
Channel Maps of Deconvolved Model and Residuals

In Figures 16–18, we provide the every fifth channel map of
the jointly deconvolved model, interferometric residual, and
single-dish residual cubes.

<figure><img src="image_19.png" /></figure>

**Figure 16.** The channel maps of the jointly deconvolved cube of M51, combining the FEASTS and THINGS data. Every fifth channel is shown (channel width: 5.153 km s⁻¹). The color coding is in asinh scale, with the luminosity below √σintσsd = 3.665 × 10⁻² mJy beam<sup>1</sup>LA (the horizontal black dotted line on the color bar) being approximately linear. In the first channel map, the black contour indicates the region where we correct the THINGS visibility using A. H. Rots et al. (1990) data (see Section 2.2).
<figure><img src="image_20.png" /><figcaption>**Figure 17.** The same as Figure 16, but for the interferometric residual (with PB attenuation). The color coding is linear, with $3\sigma_{\text{int}} = 3 \times 0.379 \, \text{mJy beam}_{\text{VLA}}^{1}$ as the limit.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 18. The same as Figure 16, but for the single-dish residual. The color coding is linear, with $3\sqrt{\sigma_{\text{int}}\sigma_{\text{sd}}}=3 \times 5.05$ mJy beam$^{-1}_{\text{FAST}}$ as the limit.</figcaption></figure>

ORCID iDs

Xuchen Lin (林旭辰) https://orcid.org/0000-0002-4250-2709
Jing Wang (王菁) https://orcid.org/0000-0002-6593-8820
Lister Staveley-Smith https://orcid.org/0000-0002-8057-0294
Suoqing Ji (季索清) https://orcid.org/0000-0001-9658-0588
Dong Yang (杨冬) https://orcid.org/0000-0002-5679-3447
Xinkai Chen (陈新凯) https://orcid.org/0000-0002-5016-6901
Fabian Walter https://orcid.org/0000-0003-4793-7880
Hsiao-Wen Chen https://orcid.org/0000-0001-8813-4182
Luis C. Ho https://orcid.org/0000-0001-6947-5846
Peng Jiang (姜鹏) https://orcid.org/0009-0003-3527-8520
Nir Mandelker https://orcid.org/0000-0001-8057-5880
Se-Heon Oh https://orcid.org/0000-0002-8379-0604
Bo Peng https://orcid.org/0000-0002-1605-0032
Céline Péroux https://orcid.org/0000-0002-4288-599X
Zhijie Qu (屈稚杰) https://orcid.org/0000-0002-2941-646X
Q. Daniel Wang https://orcid.org/0000-0002-9279-4041

References

Armillotta, L., Fraternali, F., & Marinacci, F. 2016, MNRAS, 462, 4157
Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167
Aung, H., Mandelker, N., Dekel, A., et al. 2024, MNRAS, 532, 2965
Bacchini, C., Fraternali, F., Iorio, G., & Pezzulli, G. 2019, A&A, 622, A64
Begelman, M. C., & Fabian, A. C. 1990, MNRAS, 244, 26P
Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717, 379
Berg, M. A., Lehner, N., Howk, J. C., et al. 2023, ApJ, 944, 101
Bordoloi, R., Lilly, S. J., Knobel, C., et al. 2011, ApJ, 743, 10
Borlaff, A. S., Lopez-Rodriguez, E., Beck, R., et al. 2021, ApJ, 921, 128
Bradley, L., Deil, C., Ginsburg, A., et al., 2023 astropy/regions, v0.8, Zenodo, doi:10.5281/zenodo.10144914
Braginskii, S. I. 1965, RvPP, 1, 205
CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 114501
Chen, H.-W., Qu, Z., Rauch, M., et al. 2023, ApJL, 955, L25
Chen, H.-W., Zahedy, F. S., Boettcher, E., et al. 2020, MNRAS, 497, 498
Chung, A., van Gorkom, J. H., Kenney, J. D. P., & Vollmer, B. 2007, ApJL, 659, L115
Clark, B. G. 1980, A&A, 89, 377

Colombo, D., Meidt, S. E., Schinnerer, E., et al. 2014, ApJ, 784, 4

Conselice, C. J. 2014, ARA&A, 52, 291

Cornwell, T. J. 1988, A&A, 202, 316

Cornwell, T. J. 2008, ISTSP, 2, 793

Cornwell, T. J., & Evans, K. F. 1985, A&A, 143, 77

Cowie, L. L., & McKee, C. F. 1977, ApJ, 211, 135

Cowie, L. L., & Songaila, A. 1977, Natur, 266, 501

Crain, R. A., & van de Voort, F. 2023, ARA&A, 61, 473

Di Teodoro, E. M., & Fraternali, F. 2015, MNRAS, 451, 3021

Di Teodoro, E. M., & Peek, J. E. G. 2021, ApJ, 923, 220

Dobbs, C. L., Theis, C., Pringle, J. E., & Bate, M. R. 2010, MNRAS, 403, 625

Duan, Q., Conselice, C. J., Li, Q., et al. 2024, MNRAS, arXiv:2407.09472

Ellison, S. L., Catinella, B., & Cortese, L. 2018, MNRAS, 478, 3447

Ellison, S. L., Fertig, D., Rosenberg, J. L., et al. 2015, MNRAS, 448, 221

Ellison, S. L., Wilkinson, S., Woo, J., et al. 2022, MNRAS, 517, L92

Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211

Ettori, S. 2015, MNRAS, 446, 2629

Faucher-Giguère, C.-A., & Oh, S. P. 2023, ARA&A, 61, 131

Ferland, G. J., Korista, K. T., Verner, D. A., et al. 1998, PASP, 110, 761

Field, G. B., Goldsmith, D. W., & Habing, H. J. 1969, ApJL, 155, L149

Fielding, D. B., Ostriker, E. C., Bryan, G. L., & Jermyn, A. S. 2020, ApJL, 894, L24

Fraternali, F. 2017, in ASSL, Gas Accretion Onto Galaxies, ed. A. Fox & R. Davé, Vol. 430 (Berlin: Springer), 323

French, D. M., & Wakker, B. P. 2020, ApJ, 897, 151

Gaspari, M., McDonald, M., Hamer, S. L., et al. 2018, ApJ, 854, 167

Ginsburg, A., Robitaille, T., & Beaumont, C., 2016 pvextractor: Position-Velocity Diagram Extractor, Astrophysics Source Code Library, ascl:1608.01010

Giovanelli, R., Haynes, M. P., Kent, B. R., et al. 2005, AJ, 130, 2598

Gooch, R. 1996, in ASPC, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes (San Francisco, CA: ASP), 80

Grand, R. J. J., van de Voort, F., Zjupa, J., et al. 2019, MNRAS, 490, 4786

Gronke, M., & Oh, S. P. 2020a, MNRAS, 492, 1970

Gronke, M., & Oh, S. P. 2020b, MNRAS, 494, L27

Gull, S. F., & Daniell, G. J. 1978, Natur, 272, 686

Gunn, J. E., & Gott, J. R., III. 1972, ApJ, 176, 1

Hafen, Z., Stern, J., Bullock, J., et al. 2022, MNRAS, 514, 5056

Han, J. L. 2017, ARA&A, 55, 111

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Heyer, M., Krawczyk, C., Duval, J., & Jackson, J. M. 2009, ApJ, 699, 1092

Hodges-Kluck, E. J., Miller, M. J., & Bregman, J. N. 2016, ApJ, 822, 21

Högbom, J. A. 1974, A&AS, 15, 417

Howard, S., & Byrd, G. G. 1990, AJ, 99, 1798

Huang, Y.-H., Chen, H.-W., Johnson, S. D., & Weiner, B. J. 2016, MNRAS, 455, 1713

Ji, S., Oh, S. P., & Masterson, P. 2019, MNRAS, 487, 737

Jörsäter, S., & van Moorsel, G. A. 1995, AJ, 110, 2037

Kacprzak, G. G., Muzahid, S., Churchill, C. W., Nielsen, N. M., & Charlton, J. C. 2015, ApJ, 815, 22

Karachentsev, I. D., Karachentseva, V. E., Huchtmeier, W. K., & Makarov, D. I. 2004, AJ, 127, 2031

Kauffmann, G., Li, C., Zhang, W., & Weinmann, S. 2013, MNRAS, 430, 1447

Klimenko, V. V., Kulkarni, V., Wake, D. A., et al. 2023, ApJ, 954, 115

Koda, J., Teuben, P., Sawada, T., Plunkett, A., & Fomalont, E. 2019, PASP, 131, 054505

Kooij, R., Grønnow, A., & Fraternali, F. 2021, MNRAS, 502, 1263

Kourkchi, E., & Tully, R. B. 2017, ApJ, 843, 16

Krumholz, M. R., Burkhart, B., Forbes, J. C., & Crocker, R. M. 2018, MNRAS, 477, 2716

Kuncic, Z., & Bicknell, G. V. 2004, ApJ, 616, 669

Lehner, N. 2017, in ASSL, Gas Accretion Onto Galaxies, ed. A. Fox & R. Davé, Vol. 430 (Berlin: Springer), 117

Lehner, N., Berek, S. C., Howk, J. C., et al. 2020, ApJ, 900, 9

Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2015, ApJ, 801, 25

Leroy, A. K., Sandstrom, K. M., Lang, D., et al. 2019, ApJS, 244, 24

Li, W., Nair, P., Rowlands, K., et al. 2023, MNRAS, 523, 720

Li, Z., Hopkins, P. F., Squire, J., & Hummels, C. 2020, MNRAS, 492, 1841

Mandelker, N., Nagai, D., Aung, H., et al. 2019, MNRAS, 484, 1100

Mandelker, N., Nagai, D., Aung, H., et al. 2020, MNRAS, 494, 2641

Manna, S., & Roy, S. 2023, ApJ, 944, 86

Martin, C. L., Ho, S. H., Kacprzak, G. G., & Churchill, C. W. 2019, ApJ, 878, 84

McClure-Griffiths, N. M., Stanimirović, S., & Rybarczyk, D. R. 2023, ARA&A, 61, 19

McKee, C. F., & Cowie, L. L. 1977, ApJ, 215, 213

Miller, M. J., & Bregman, J. N. 2015, ApJ, 800, 14

Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319

Moreno, J., Torrey, P., Ellison, S. L., et al. 2019, MNRAS, 485, 1320

Naab, T., & Ostriker, J. P. 2017, ARA&A, 55, 59

Nelson, D., Sharma, P., Pillepich, A., et al. 2020, MNRAS, 498, 2391

Oh, S.-H., Kim, S., For, B.-Q., & Staveley-Smith, L. 2022, ApJ, 928, 177

Oikawa, S., & Sofue, Y. 2014, PASJ, 66, 77

Ostriker, E. C., McKee, C. F., & Leroy, A. K. 2010, ApJ, 721, 975

Padnos, D., Mandelker, N., Birnboim, Y., et al. 2018, MNRAS, 477, 3293

Pan, H.-A., Lin, L., Hsieh, B.-C., et al. 2019, ApJ, 881, 119

Parsons, W., 3rd Earl of Rosse 1850, RSPT, 140, 499

Paturel, G., Petit, C., Prugniel, P., et al. 2003, A&A, 412, 45

Péroux, C., & Howk, J. C. 2020, ARA&A, 58, 363

Péroux, C., Nelson, D., van de Voort, F., et al. 2020, MNRAS, 499, 2462

Péroux, C., Weng, S., Karki, A., et al. 2022, MNRAS, 516, 5618

Pety, J., Schinnerer, E., Leroy, A. K., et al. 2013, ApJ, 779, 43

Plunkett, A., Hacar, A., Moser-Fischer, L., et al. 2023, PASP, 135, 034501

Prochaska, J. X., Werk, J. K., Worseck, G., et al. 2017, ApJ, 837, 169

Putman, M. E. 2017, in ASSL, Gas Accretion Onto Galaxies, ed. A. Fox & R. Davé, Vol. 430 (Berlin: Springer), 1

Quai, S., Hani, M. H., Ellison, S. L., Patton, D. R., & Woo, J. 2021, MNRAS, 504, 1888

Rahmati, A., Pawlik, A. H., Raičević, M., & Schaye, J. 2013, MNRAS, 430, 2427

Ramesh, R., Nelson, D., & Pillepich, A. 2023, MNRAS, 518, 5754

Richter, P. 2017, in ASSL, Gas Accretion Onto Galaxies, ed. A. Fox & R. Davé, Vol. 430 (Berlin: Springer), 15

Rots, A. H., Bosma, A., van der Hulst, J. M., Athanassoula, E., & Crane, P. C. 1990, AJ, 100, 387

Saintonge, A., Catinella, B., Tacconi, L. J., et al. 2017, ApJS, 233, 22

Sault, R. J. 1990, ApJL, 354, L61

Sault, R. J., Staveley-Smith, L., & Brouw, W. N. 1996, A&AS, 120, 375

Sault, R. J., Teuben, P. J., & Wright, M. C. H. 1995, in ASPC, Vol. 77, Astronomical Data Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes (San Francisco, CA: ASP), 433

Serra, P., Westmeier, T., Giese, N., et al. 2015, MNRAS, 448, 1922

Shangguan, J., Ho, L. C., Li, R., et al. 2019, ApJ, 870, 104

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Sorgho, A., Foster, T., Carignan, C., & Chemin, L. 2019, MNRAS, 486, 504

Sparre, M., & Springel, V. 2017, MNRAS, 470, 3946

Spitzer, L. 1962, Physics of Fully Ionized Gases (2nd edn.; New York, NY: Interscience)

Stanimirovic, S. 2002, in ASPC, Vol. 278, Single-Dish Radio Astronomy: Techniques and Applications, ed. S. Stanimirovic et al. (San Francisco, CA: ASP), 375

Steer, D. G., Dewdney, P. E., & Ito, M. R. 1984, A&A, 137, 159

Stein, Y., Dettmar, R.-J., Weżgowiec, M., et al. 2019, A&A, 632, A13

Stern, J., Fielding, D., Hafen, Z., et al. 2024, MNRAS, 530, 1711

Sun, J., Leroy, A. K., Ostriker, E. C., et al. 2020a, ApJ, 892, 148

Sun, J., Leroy, A. K., Schinnerer, E., et al. 2020b, ApJL, 901, L8

Sun, J., Leroy, A. K., Schruba, A., et al. 2018, ApJ, 860, 172

Sylos Labini, F., Benhaiem, F., Comerón, S., & López-Corredoira, M. 2019, A&A, 622, A58

Sylos Labini, F., Straccamore, M., De Marzo, G., & Comerón, S. 2023, MNRAS, 524, 1560

Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157

Tamburro, D., Rix, H.-W., Walter, F., et al. 2008, AJ, 136, 2872

Tan, B., Oh, S. P., & Gronke, M. 2021, MNRAS, 502, 3179

Theis, C., & Spinneker, C. 2003, Ap&SS, 284, 495

Thompson, A. R., Moran, J. M., & Swenson, G. W., Jr. 2017, Interferometry and Synthesis in Radio Astronomy (3rd ed.; Berlin: Springer), 89

Thorp, M. D., Ellison, S. L., Simard, L., Sánchez, S. F., & Antonio, B. 2019, MNRAS, 482, L55

Tumlinson, J., Werk, J. K., Thom, C., et al. 2011, ApJ, 733, 111

Turner, J. L., & Ho, P. T. P. 1994, ApJ, 421, 122

Vallat, R. 2018, JOSS, 3, 1026

van de Voort, F., Springel, V., Mandelker, N., van den Bosch, F. C., & Pakmor, R. 2019, MNRAS, 482, L85

Van Eck, C. L., Brown, J. C., Shukurov, A., & Fletcher, A. 2015, ApJ, 799, 35

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Vogelsberger, M., Marinacci, F., Torrey, P., & Puchwein, E. 2020, NatRP, 2, 42

Voit, G. M., Meece, G., Li, Y., et al. 2017, ApJ, 845, 80

Vugrin, K. W., Swiler, L. P., Roberts, R. M., Stucky-Mack, N. J., & Sullivan, S. P. 2007, WRR, 43, W03423
Walter, F., Brinks, E., de Blok, W. J. G., et al. 2008, AJ, 136, 2563

Wang, J., Catinella, B., Saintonge, A., et al. 2020, ApJ, 890, 63

Wang, J., Lin, X., Staveley-Smith, L., et al. 2024a, ApJ, 973, 15

Wang, J., Lin, X., Yang, D., et al. 2024b, ApJ, 968, 48

Wang, J., Yang, D., Lin, X., et al. 2025, ApJ, 980, 25

Wang, J., Yang, D., Oh, S.-H., et al. 2023, ApJ, 944, 102

Weng, S., Péroux, C., Karki, A., et al. 2023, MNRAS, 523, 676

Westmeier, T., Kitaeff, S., Pallot, D., et al. 2021, MNRAS, 506, 3962

White, M. C., Bicknell, G. V., Sutherland, R. S., Salmeron, R., & McGregor, P. J. 2016, MNRAS, 455, 2042

Wilde, M. C., Werk, J. K., Burchett, J. N., et al. 2021, ApJ, 912, 9

Yang, D., Wang, J., Qu, Z., et al. 2025, arXiv:2502.16926

Yang, Y., & Ji, S. 2023, MNRAS, 520, 2148

Yu, H., Zhu, M., Xu, J.-L., et al. 2023, MNRAS, 521, 2719

Yu, Q., Fang, T., Feng, S., et al. 2022, ApJ, 934, 114

Zhu, M., Yu, H., Wang, J., et al. 2021, ApJL, 922, L21