# The SDSS-Gaia View of the Color–Magnitude Relation for Blue Horizontal-
branch Stars

Fabrícia O. Barbosa Universidade de São Paulo, Instituto de Astronomia,
Geofísica e Ciências Atmosféricas, Departamento de Astronomia,
SP 05508-090, São Paulo, Brazil  Rafael M. Santucci Universidade Federal de
Goiás, Instituto de Estudos Socioambientais, Planetário, Goiânia, GO
74055-140, Brazil  Universidade Federal de Goiás, Campus Samambaia, Instituto
de Física, Goiânia, GO 74001-970, Brazil  Silvia Rossi Universidade de São
Paulo, Instituto de Astronomia, Geofísica e Ciências Atmosféricas,
Departamento de Astronomia,
SP 05508-090, São Paulo, Brazil  Guilherme Limberg Universidade de São Paulo,
Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Departamento de
Astronomia,
SP 05508-090, São Paulo, Brazil  Angeles Pérez-Villegas Instituto de
Astronomía, Universidad Nacional Autónoma de México, Apartado Postal 106, C.
P. 22800, Ensenada, B. C., Mexico  Hélio D. Perottoni Universidade de São
Paulo, Instituto de Astronomia, Geofísica e Ciências Atmosféricas,
Departamento de Astronomia,
SP 05508-090, São Paulo, Brazil  Institut de Ciències del Cosmos (ICCUB),
Universitat de Barcelona (IEEC-UB), Martí i Franquès 1, E08028 Barcelona,
Spain

###### Abstract

We present an updated sample of blue horizontal-branch (BHB) stars selected
from the photometric and spectroscopic data from Sloan Digital Sky Survey and
its associated project Sloan Extension for Galactic Understanding and
Exploration (SEGUE). With this data, we selected candidates for A-type stars
in the color-color space and then a mixture modeling technique was implemented
in order to distinguish between BHB and main-sequence/blue-straggler stars
based on their surface gravity values ($\log\rm{g}$) estimated by the SEGUE
Stellar Parameter Pipeline. Our robust approach allows us to attribute
individual probabilities of each star truly being in the BHB stage. Hence, our
method is advantageous in comparison to previous SEGUE BHB selections that
adopted simple $\log\rm{g}$ cuts. We also revisit the color–magnitude relation
for these stars and propose two calibrations, based on updated distances for
Galactic globular clusters, to estimate absolute magnitudes with $(g-r)_{0}$
and $(u-r)_{0}$ colors.

Galaxy: stellar halo – stars: horizontal branch – stars: distances

##  1 Introduction

The Gaia mission (Gaia Collaboration et al., 2016) has provided a better
understanding of the Galaxy, in particular regarding the field of Galactic
Archaeology (Helmi, 2020; Brown, 2021). The astrometric information provided
for an unprecedented number of objects has dramatically changed the way we
study the Galactic halo (e.g., Belokurov et al., 2018; Myeong et al., 2018;
Koppelman et al., 2018; Malhan et al., 2018).

Despite the huge amount of direct measurements supplied by Gaia, distances
inferred from brightness are still of great value. At magnitude $G<15$, the
early third data release (EDR3) presents parallax uncertainties of
${\sim}0.02$ mas (Gaia Collaboration et al., 2021), and they increase
significantly for fainter stars. To overcome this limitation, we can use
various well-known distance tracers such as RR Lyrae (Shapley, 1916), Cepheids
(Leavitt & Pickering, 1912), and blue horizontal-branch (BHB) stars (Cacciari,
1999).

BHBs are metal-poor ([Fe/H]111[A/B]
$=\log\rm(N_{A}/N_{B})_{\star}-\log(N_{A}/N_{B})_{\odot}$, where $\rm N_{A}$
and $\rm N_{B}$ are the number density of atoms of the elements A and B,
respectively. $\star$ refers to the considered star, and $\odot$ refers to the
Sun. $\lesssim-0.5$; Santucci et al. 2015b) A or B-type stars that burn helium
in their cores. These evolved stars present a high and nearly constant
luminosity, making them perfect for investigating the outer regions of the
halo and the assembly history of our Galaxy (Xue et al., 2011, 2008; Deason et
al., 2011, 2017; Belokurov et al., 2014; Santucci et al., 2015b). In recent
works, BHBs were used to study dynamical substructures and stellar streams
(Yuan et al., 2019, 2020, 2022; Peñarrubia & Petersen, 2021; Li et al., 2022;
Wu et al., 2022), the connection between the apocenter pile-up of orbits and
the so-called “break-radius” of the stellar halo density profile (Deason et
al., 2018), the anisotropy of the halo velocity distribution (Lancaster et
al., 2019), the age gradient of the halo out to $\sim$35 kpc (Whitten et al.,
2019), to estimate the total dynamical mass of the Milky Way (Deason et al.,
2021; Bird et al., 2022), and even to demonstrate the influence of the Large
Magellanic Cloud in our Galaxy’s halo (Erkal et al., 2021; Petersen &
Peñarrubia, 2021).

The well-defined structure of the horizontal branch in the color-magnitude
diagram (CMD), a roughly constant luminosity, permits the development of a
distance calibration for these BHBs. The first approximation developed was a
linear fit, using the ($B-V$) color and absolute magnitude in the $V$ band,
for stars in globular clusters (Hayes & Philip, 1979). Likewise, Preston et
al. (1991) defined a smoother relation, a fourth degree polynomial, for the
same color-magnitude space. Two decades later, a widely used calibration was
presented by Deason et al. (2011, hereafter D11) based on magnitudes in the
$ugriz$ system for the Sloan Digital Sky Survey (SDSS; York et al., 2000)
eighth data release (DR8; Aihara et al., 2011), which had its color range
extended by Belokurov & Koposov (2016) afterwards. In the meantime, Fermani &
Schönrich (2013) argued that it is extremely important to take into account
the effect of the metallicity on the absolute magnitude estimation, proposing
a new calibration based on a statistical method. However, Santucci et al.
(2015a) and Utkin & Dambis (2020) showed that the differences between
considering or not the metallicity in the relations are negligible, with D11’s
estimates being $2.5\%$ higher, within (1-$\sigma$) errors of both
calibrations.

D11’s relation still remains the most used calibration for BHB stars (Santucci
et al., 2015a, b; Thomas et al., 2018; Whitten et al., 2019; Donlon et al.,
2020; Martin et al., 2022) even though photometric data have been updated
several times since then. Moreover, we can now compare photometric distances
of BHB stars with purely geometric estimates from Gaia’s parallaxes (e.g.,
Bailer-Jones et al. 2021) as well as new measurements for Galactic globular
clusters (Vasiliev & Baumgardt, 2021). These facts bring to light the
relevance of reviewing the D11’s calibration with recent data.

This paper is organized as follows. In Section 2, we describe the photometric
selection and revise a previous method to identify BHB stars. Section 3
presents the selection of stars in globular clusters and the method used to
define the absolute magnitude calibration. Finally, in Section 4 we discuss
our results.

##  2 Data

###  2.1 A-type stars

The initial selection of A-type stars was made using the photometry from the
sixteenth data release (DR16) of SDSS (Ahumada et al., 2020). For the
selection of BHB stars, we were specially interested in the spectroscopic data
obtained by the Sloan Extension for Galactic Understanding and Exploration
(SEGUE; Yanny et al., 2009) processed by the SEGUE Stellar Parameter Pipeline
(SSPP; Lee et al., 2008a, b)222Last run on DR9 (Ahn et al., 2012; Rockosi et
al., 2022)..

We implemented color cuts applying the following criteria:
$-0.3<(g-r)_{0}<0.1$ and $0.8<(u-g)_{0}<1.4$, similar to those used in
previous works (Sirko et al. 2004, D11). All the magnitudes were corrected
using the extinction coefficients ($A_{g}$, $A_{r}$, $A_{u}$) provided by the
SDSS catalog itself and we removed stars with relative errors in the $g$-band
magnitude greater than 1%.

The photometric selection is able to exclude several undesired objects, such
as white dwarfs, quasars and cooler spectral types (Yanny et al., 2009;
Vickers et al., 2012), but the major source of contamination, blue straggler
stars (BSSs), remains. The distinction between evolved and main-sequence
stars/BSSs is commonly made by investigating spectral features, specially
Balmer lines, whose depths are affected by effective temperature ($T_{\rm
eff}$) and widths by surface gravity ($\log\rm{g}$). With the output of SSPP,
we can directly inspect these stellar atmospheric parameters. Therefore, we
cross-matched the filtered sample with the SSPP catalog using
$5^{\prime\prime}$ radius. In addition to color filters, we restricted our
sample to stars with moderate signal-to-noise ratio ($S/N>10$) and $7500\,{\rm
K}<T_{\rm eff}<10000\,{\rm K}$ (Deason et al., 2012; Santucci et al., 2015a),
where $T_{\rm eff}$ is the estimate adopted by the pipeline. Duplicated stars
with the smallest $S/N$ were removed, which resulted in 16463 stars. The
restrictions above remove poor-quality data and cooler stars that could remain
after the color cut, which assures that contamination from other non-BSS stars
is minimal.

###  2.2 BHB stars

One of the techniques used to disentangle BSSs and BHB stars is the $f_{m}$
versus $D_{0.2}$ method (Pier, 1983), where $f_{m}$ is the minimum flux
relative to the continuum level and $D_{0.2}$ is a measurement of the line
width of Balmer lines, so it provides indirect information regarding both
$T_{\rm eff}$ and $\log\rm{g}$. Later, a different approach was proposed by
Clewley et al. (2002) based on the parameters of the Sérsic profile (Sersic,
1968), which describes the shape of the lines.

BSSs present a stronger $\log\rm{g}$ than those located in the horizontal
branch. Santucci et al. (2015a) showed that these stellar types are clearly
distinguishable for magnitudes $g_{0}<18$ with SEGUE/SDSS DR8 data, being
possible to classify them by fitting a combination of two Gaussian functions
to their $\log\rm{g}$ distributions. This method was proved to be in good
concordance with spectral analysis, with more than 90% of agreement. When
replicating this procedure with current SEGUE data, we noticed a change in the
peaks of both groups and a greater overlap in the $\log\rm{g}$ distribution as
presented in the left column of Fig. 1. This is observed even for relatively
bright A-type stars ($g_{0}<18$, top panel), which makes it more difficult to
separate these objects with that simple approach (dashed lines indicate the
Gaussian fits from Santucci et al. 2015a). The differences are probably due to
changes between the releases of SDSS on the $\log\rm{g}$ estimates considered
to obtain the final adopted parameter333See
https://www.sdss.org/dr16/spectro/sspp_changes/ for detailed information..

Figure 1: Histograms of $\log\rm{g}$. Left column: log g adopted by the
pipeline for stars with $g_{0}<18$ (top) and $g_{0}>18$ (bottom). Dashed lines
are the Gaussian distributions defined by Santucci et al. (2015a). Right
column: log g estimates provided by SSPP spectroscopically determined (top)
and from ANNRR method (bottom) for all stars.

###  2.3 Classification

Given the two-Gaussian-like morphology of the $\log\rm{g}$ distributions
observed for our sample of A-type stars (Fig. 1), we used a Gaussian Mixture
Model (GMM) unsupervised approach in order to distinguish BHBs and BSSs. For
this task, we utilize the scikit-learn (Pedregosa et al., 2011)
GaussianMixture444https://scikit-
learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html#sklearn.mixture.GaussianMixture.
package. In this GMM implementation, the expectation-maximization algorithm
(Dempster et al., 1977) is employed in the search for the best-fit model.

The GMM technique fits the data as a finite combination of $K$ Gaussian
distributions. As made previously by Santucci et al. (2015a), $K$ was defined
based on visual inspection of $\log\rm{g}$ estimates presented in Fig. 1 and
the assumption that the contamination is predominantly of main-sequence
stars/BSSs. Therefore, $K=2$ is an adequate value for the sample. Moreover,
GMM can be readily applied to data of arbitrary dimensionality. Therefore, we
take advantage of such flexibility and explore a suitable combination of
$\log\rm{g}$ estimates provided by SSPP (we refer the reader to Lee et al.
2008a for details about different approaches to determine $\log\rm{g}$ from
SEGUE spectra). We noticed that the distributions of both $\log\rm{g}_{\rm
ANNRR}$ and $\log\rm{g}_{\rm SPEC}$ exhibit clearly two peaks, as expected for
the BHBs/BSSs dichotomy, while it is not possible to observe this feature in
others. These two distributions are shown in the right column of Fig. 1.
Hence, we proceeded with the GMM separation within the two-dimensional space
defined by these $\log\rm{g}$ estimates. The final $\log\rm{g}$ adopted by the
pipeline was not considered an extra dimension as it consists of a weighted
mean of the valid estimates.

In order to guarantee the robustness of our method against uncertainties
reported by the SSPP, we constructed a set of $10^{4}$ realizations of each
star’s $\log\rm{g}$ estimates in a Monte Carlo framework. Then, we performed
the GMM classification for all iterations. Finally, the fraction of instances
that a star is attributed to a certain class (either BHB or BSS) is taken as
its membership probability for that given group. For this procedure, stars
without valid estimates of both $\log\rm{g}_{\rm ANNRR}$ and $\log\rm{g}_{\rm
SPEC}$ are removed. With this strategy, we achieved a sample of 5699/4590
stars classified as BHBs above 50%/99% probability555The full sample is
available at https://github.com/guilhermelimberg/bhb_dist..

The final classification obtained is shown in Fig. 2. The difference in the
uncertainties of the estimates greatly influences the classification, as the
$\log\rm{g}_{\rm ANNRR}$ presents more precise values ($\sim 0.06$) than
$\log\rm{g}_{\rm SPEC}$ ($\sim 0.21$). We cross-matched our sample with the
one from Santucci et al. (2015a) to evaluate the fraction of BSS
contamination. 10% of our BHB set was classified previously as BSSs, and,
among those with a probability greater than 99% of being BHB following the
method implemented here, 2% is possibly incorrectly assigned.

Figure 2: Distribution of classified stars with $g_{0}<18$ (left) and
$g_{0}>18$ (right) in the surface gravity space. Median errors are indicated
in the bottom right corner. Histograms show the distribution of log g from the
respective axis for stars classified as BHB and BSS. Colors indicate the
probability of being a BHB star.

##  3 Absolute magnitude calibration

###  3.1 Globular clusters stars

The procedure to construct an absolute magnitude relation follows previous
works (see Section 1), starting with the selection of BHB stars in globular
clusters. We used the photometric catalog from An et al. (2008), which
provides magnitudes for crowded fields observed by SDSS. The clusters
presenting a well defined horizontal branch were selected and their magnitudes
were corrected using the standard extinction ($E(B-V)$) from Schlegel et al.
(1998) along with the relative extinctions from Wang & Chen (2019) for $g-,u-$
and $r-$band. Vasiliev & Baumgardt (2021) attributed a membership probability
for stars in globular clusters based on proper motions and parallaxes from
Gaia EDR3. We selected stars from several globular clusters that were more
likely than 0.99 to belong to those clusters and we obtained their absolute
magnitude in the SDSS $g$-band ($M_{g}$) with the estimated distance for each
cluster given by these authors. The list of clusters, their heliocentric
distances and distance moduli are presented in Table 1.

Table 1: Heliocentric distances provided by Vasiliev & Baumgardt (2021) for each globular cluster, uncertainties, and their respective distance moduli. Cluster | D (kpc) |  $\sigma_{D}$ (kpc) |  $(m-M)_{0}$ (mag)
---|---|---|---
NGC2419 | 83.0 | 1.5 | 19.59
NGC4147 | 18.65 | 0.16 | 16.35
NGC5024, M53 | 18.59 | 0.15 | 16.35
NGC5053 | 17.30 | 0.14 | 16.19
NGC5272, M3 | 10.20 | 0.06 | 15.04
NGC5466 | 16.32 | 0.13 | 16.06
NGC5904, M5 | 7.49 | 0.05 | 14.37
NGC6205, M13 | 7.53 | 0.06 | 14.38
NGC6341, M92 | 8.60 | 0.05 | 14.67
NGC7078, M15 | 10.73 | 0.14 | 15.15
NGC7089, M2 | 11.62 | 0.13 | 15.33

To create the sample used to implement the calibration, we applied the limits
for colors as defined for the initial selection (see Section 2.1). Then, the
stars were selected in a single combined CMD, limiting the $M_{g}$ between
$-0.15$ and $1.15$. After this exercise, the remaining globular cluster
members were checked individually at the SIMBAD database (Wenger et al.,
2000), and those classified as variables, blue stragglers and other
undesirable types were removed. We also excluded stars with flags in the
magnitudes used, leaving us with 744 stars to derive the calibrations from.

###  3.2 Fitting the horizontal branch

Finding the best mathematical relationships to fit observable data is not an
easy task. In previous works, the absolute magnitudes for BHB stars have been
described as a high-degree polynomial (Preston et al., 1991; Deason et al.,
2011; Belokurov & Koposov, 2016). Instead of arbitrarily assuming that this
function is the best representation of the data, we explore the possible
combinations between colors and absolute magnitudes. For this task, we
employed the TuringBot software (Ashok et al., 2020), a code that performs
symbolic regression using a simulated annealing algorithm (Delahaye et al.,
2019; Chira & Plionis, 2019) in order to search for the best set of parameters
and mathematical operations to describe the data.

TuringBot is particularly interesting in this case, because it allows the
visualization of the estimated mathematical laws, allowing the user to choose
the most appropriate equations for their needs. Furthermore, the user is free
to choose the mathematical operations involved in the fitted functions, the
error metric for convergence, as well as the input variables. The best fits
are presented in a summarized box, combining the error and the complexity of
the equations. The complexity is defined by the sum of the “size” of the
mathematical operations, constants, and variables present in the solutions.
The program assumes that an input variable, constant, sum, subtraction, and
multiplication have size 1 each, division has size 2, and more complex
operations have higher sizes666More TuringBot details can be found in the
program documentation, available at:
https://turingbotsoftware.com/documentation.html.. We verified that the use of
very complex mathematical operations is unnecessary and does not improve the
average error of the equations presented by the software. Hence, we adopted
only the basic mathematical operations (sum, subtraction, multiplication, and
division) as input for the search for absolute magnitude calibrations. The
mean absolute error was used as a criterion for convergence and we tested the
dependence of all the most common available observable variables found in the
literature for estimates of this type, such as magnitudes, color indices, and
metallicity.

After evaluating all the combinations of input variables ($u_{0}$, $g_{0}$,
$r_{0}$, $(u-g)_{0}$, $(u-r)_{0}$, $(g-r)_{0}$, and [Fe/H]) presented in
Appendix A, we found that there is no significant dependence on metallicity in
the calibrations provided by the code, regardless of the mathematical
operations adopted and the algorithm convergence time. The colors $(u-g)_{0}$
and $(u-r)_{0}$ provided calibrations with smaller errors than the color
$(g-r)_{0}$, traditionally used in the absolute magnitude calibration of BHB
stars (D11; Belokurov & Koposov 2016), and also smaller than $(g-r)_{0}$ with
[Fe/H], which means we can achieve more accurate results that do not require
metallicity information. The observed improvement with $(u-g)_{0}$ and
$(u-r)_{0}$ color might be associated with the $u$ filter, whose transmission
curve is mostly between 3000Å and 4000Å, i.e., it is positioned in a region of
the spectrum where the Balmer discontinuity ($\sim$3645Å) is located, as well
as several Hydrogen lines from the Balmer series, which makes it a useful
indirect indicator of $T_{\rm eff}$ and $\log\rm{g}$ of the BHBs, atmospheric
parameters that are directly linked to the mass of the stars in the horizontal
branch (Valcarce & Catelan, 2008).

Fig. 3 shows the associated errors for each fit in the final BHB sample.
Clearly, $(u-r)_{0}$ presents a better performance than $(g-r)_{0}$ and
$(u-g)_{0}$. Using the same tool, we find that the relation proposed in D11 is
a function of complexity 33 and ${\rm error}=0.12$, whilst equations of lower
order present a much lower complexity with errors of ${\sim}0.10$.

We chose the first functions from which there is no significant decrease in
error, i.e., functions of complexity 6 in Fig. 3, as those that best describes
the data. Exists a singularity in the calibrations, however it is outside of
our color range. Hence, it does not imply an obstacle to their usage in the
context of this work.

Figure 3: Error comparison between fits for colors $(g-r)_{0}$, $(u-g)_{0}$
and $(u-r)_{0}$ from TuringBot.

$M_{g}=\frac{0.178}{0.537+(g-r)_{0}}$ (1)

$M_{g}=\frac{0.721}{(u-r)_{0}}-0.212$ (2)

###  3.3 Distances analysis

Left panels in Fig. 4 show the distribution of BHB stars in the CMD with color
$(g-r)_{0}$ (top) and $(u-r)_{0}$ (bottom). In the top left panel, we can
observe how the calibration proposed here (Eq. 1) provides magnitudes lower
than D11’s, which results in larger distances. The difference is minimal at
$(g-r)_{0}\sim-0.20$, where both equation come closer, and the smaller values
are a consequence of the inclusion of the cluster NGC7078, whose stars are
brighter and were not included in D11. On the other hand, the distribution
using $(u-r)_{0}$ has a lower dispersion (bottom left panel).

Figure 4: Left panels: color–magnitude diagrams showing BHBs used to define
the calibrations. Dash-dotted line represents the polynomial fit defined in
D11. Solid lines represent the calibration for $(g-r)_{0}$ and $(u-r)_{0}$
color presented in this work (Eq. 1 and 2, respectively, from top to bottom).
Median errors of the data are indicated in the bottom right corner of each
panel. Right panels: difference between distances calculated with D11’s
calibration and those presented in this work, respectively $(g-r)_{0}$ (top)
and $(u-r)_{0}$ (bottom).

In the right panels (Fig. 4), we also show the comparison between distances
estimated with the relation from D11 and each calibration defined in the
present work. For consistence with D11’s relation, only stars bluer than
$(g-r)_{0}=0$ were considered, and we rejected stars with BHB probabilities of
less than 0.99 to reduce the number of misclassified stars. The new
calibration using color $(g-r)_{0}$ provides distances about 5% larger than
D11’s for the reddest stars, while the other end of the color window attains a
relative difference of up to 9% (top right panel). For the color $(u-r)_{0}$,
the scatter is more uniform and much larger for the bluest stars (bottom right
panel).

When comparing with purely astrometric heliocentric distances, there is a
considerable scatter, even for stars closer than 5 kpc. For this comparison,
we selected stars with relative parallax uncertainty from Gaia EDR3777Gaia
Collaboration (2020). in the interval $0<\sigma_{\varpi}/\varpi<0.2$, re-
normalized unit weight errors within the recommended range
($\texttt{RUWE}<1.4$; Lindegren 2018), and also a BHB probability greater than
0.99 (${>}\ 300$ stars). In Fig. 5, we show the comparison between geometric
(left) and photogeometric (right) distances provided by Bailer-Jones et al.
(2021) and our calibration using $(u-r)_{0}$ color. For fainter stars, both
Bailer-Jones et al.’s (2021) distances are frequently underestimated. The gray
region indicates the interval within 20% of distance in the respective
horizontal axes, and we find 65% of the stars when using photogeometric
estimates and 54% with the purely geometric outside it. Gaia’s parallax
measurements potentially are not accurate enough for these BHBs and so the
final results are not representative of the sample (since the distances
inferred from Bayesian methods are strongly dependent on the measured
parallax). We also point out that this effect is unlikely the consequence of
an inappropriate classification since 94% of the stars possess
$\log\rm{g}_{\rm ADOP}<3.6$ and would receive the same label by Santucci et
al.’s (2015a) method. Finally, we cannot endorse the compatibility between
D11’s distances for BHB stars found in the Pristine survey (Starkenburg et
al., 2019) and Gaia DR2’s parallaxes. As the Pristine data are not publicly
available, it is not possible to evaluate whether the difference is due to the
BHB sample used.

Figure 5: Comparison between the distances estimated by Eq. 2
($D_{(u-r)_{0}}$), geometric ($D_{geo}$, left) and photogeometric ($D_{phot}$,
right) distances from Bailer-Jones et al. (2021). Colors indicate Gaia G
magnitude and gray area indicates the region of $1-\sigma$. Median errors are
indicated in the bottom right corner of the top panels.

Similar inconsistencies between photometric and astrometic-inferred distances
were also observed by previous works. Using OB stars, Shull & Danforth (2019)
noted an increase in the discrepancies at ${d>1.5}$ kpc, with B-type stars
showing smaller values of distances when considering parallaxes alone. Our
A-type stars sample seems to follow this same trend.

##  4 Discussion & summary

Since the last absolute magnitude calibration published for BHB stars we had
the advent of Gaia data, which allowed us to review the previous relationship
thanks to the better characterization of globular clusters (Vasiliev &
Baumgardt, 2021).

Using data from the SSPP catalog, we obtained a sample of ${\sim}5700$ BHB
stars implementing the GMM algorithm. This new approach is an alternative to
the previous individual Gaussian fits, as it is clear that they do not
distinguish the latest SSPP $\log{\rm g}$ distribution properly (Fig. 1).

To find which kind of function better describes the distribution of BHBs in
the CMD, we used a software for implementing symbolic regression. We suggest
two new color–magnitude calibrations based on photometry, including a relation
with $(u-r)_{0}$ color that has not been used before. This calibration
provides more accurate estimates than $(g-r)_{0}$ color in most cases. For the
bluest stars, the differences can exceed 10% of nominal values (Eq. 2).
However, this difference decreases for redder and more distant stars. Here, we
show that the calibrations can be simpler and achieve an acceptable result,
very similar to those from by D11’s relation.

We noted substantial differences between photometric and
geometric/photogeometric distances. A possibility would be inaccurate
estimates of $\log{\rm g}$ provided by SEGUE, as it was mainly designed for
cool stars ($T_{\rm eff}<7500\,{\rm K}$). However, Santucci et al. (2015a)
showed that this is unlikely to be the case. It could also be due to an
incorrect value of extinction for the SDSS photometry, yet it also does not
seem to explain the disparity. Most stars were observed in regions of low
extinction and we could not find any relation between the extinction values
and the inconsistency observed. The observed differences also lead us to
believe that measured parallaxes for these stars could be unreliable, as there
is no agreement between the distance estimates even for the closest stars. It
would be interesting to investigate whether the same discrepancy can be
observed with other halo tracers.

The new sample made available here can help to improve results already known
about the structures of the Milky Way stellar halo. For example, these BHBs
can be used to revisit the duality of the stellar halo or re-evaluate the
Galaxy’s mass estimate.

We thank the anonymous referee for the careful review and all the suggestions,
which greatly improved our work. This research was financed with public funds,
without which it would not have been possible. F.O.B. acknowledges CAPES
(PROEX; Proc. 88887.604787/2021-00). R.M.S. acknowledges CNPq (Proc.
306667/2020-7). S.R. acknowledges partial financial support from FAPESP (Proc.
2015/50374-0 and 2014/18100-4), CAPES, and CNPq. G.L. acknowledges FAPESP
(Proc. 2021/10429-0). A.P.-V. acknowledges the DGAPA-PAPIIT grant IA103122.
H.D.P. thanks FAPESP (Procs. 2018/21250-9 and 2022/04079-0). This work has
made use of data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. Funding for the Sloan
Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the
U.S. Department of Energy Office of Science, and the Participating
Institutions. SDSS acknowledges support and resources from the Center for
High-Performance Computing at the University of Utah. The SDSS web site is
www.sdss.org. SDSS is managed by the Astrophysical Research Consortium for the
Participating Institutions of the SDSS Collaboration including the Brazilian
Participation Group, the Carnegie Institution for Science, Carnegie Mellon
University, Center for Astrophysics — Harvard & Smithsonian (CfA), the Chilean
Participation Group, the French Participation Group, Instituto de Astrofísica
de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and
Mathematics of the Universe (IPMU) / University of Tokyo, the Korean
Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut
für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA
Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-
Institut für Extraterrestrische Physik (MPE), National Astronomical
Observatories of China, New Mexico State University, New York University,
University of Notre Dame, Observatório Nacional / MCTI, The Ohio State
University, Pennsylvania State University, Shanghai Astronomical Observatory,
United Kingdom Participation Group, Universidad Nacional Autónoma de México,
University of Arizona, University of Colorado Boulder, University of Oxford,
University of Portsmouth, University of Utah, University of Virginia,
University of Washington, University of Wisconsin, Vanderbilt University, and
Yale University. This research has also made use of RStudio (RStudio Team,
2022) and TOPCAT (http://www.starlink.ac.uk/topcat/, Taylor, 2005).

## References

  * Ahn et al. (2012) Ahn, C. P., Alexandroff, R., Allende Prieto, C., et al. 2012, ApJS, 203, 21, doi: 10.1088/0067-0049/203/2/21
  * Ahumada et al. (2020) Ahumada, R., Prieto, C. A., Almeida, A., et al. 2020, ApJS, 249, 3, doi: 10.3847/1538-4365/ab929e
  * Aihara et al. (2011) Aihara, H., Allende Prieto, C., An, D., et al. 2011, ApJS, 193, 29, doi: 10.1088/0067-0049/193/2/29
  * An et al. (2008) An, D., Johnson, J. A., Clem, J. L., et al. 2008, ApJS, 179, 326, doi: 10.1086/592090
  * Ashok et al. (2020) Ashok, D., Scott, J., Wetzel, S., Panju, M., & Ganesh, V. 2020, arXiv e-prints, arXiv:2010.11328.  https://arxiv.org/abs/2010.11328
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147, doi: 10.3847/1538-3881/abd806
  * Belokurov et al. (2018) Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611, doi: 10.1093/mnras/sty982
  * Belokurov & Koposov (2016) Belokurov, V., & Koposov, S. E. 2016, MNRAS, 456, 602, doi: 10.1093/mnras/stv2688
  * Belokurov et al. (2014) Belokurov, V., Koposov, S. E., Evans, N. W., et al. 2014, MNRAS, 437, 116, doi: 10.1093/mnras/stt1862
  * Bird et al. (2022) Bird, S. A., Xue, X.-X., Liu, C., et al. 2022, arXiv e-prints, arXiv:2207.08839.  https://arxiv.org/abs/2207.08839
  * Brown (2021) Brown, A. G. A. 2021, ARA&A, 59, doi: 10.1146/annurev-astro-112320-035628
  * Cacciari (1999) Cacciari, C. 1999, in Astronomical Society of the Pacific Conference Series, Vol. 167, Harmonizing Cosmic Distance Scales in a Post-HIPPARCOS Era, ed. D. Egret & A. Heck, 140–160
  * Chira & Plionis (2019) Chira, M., & Plionis, M. 2019, MNRAS, 490, 5904, doi: 10.1093/mnras/stz2885
  * Clewley et al. (2002) Clewley, L., Warren, S. J., Hewett, P. C., et al. 2002, MNRAS, 337, 87, doi: 10.1046/j.1365-8711.2002.05864.x
  * Deason et al. (2011) Deason, A. J., Belokurov, V., & Evans, N. W. 2011, MNRAS, 416, 2903, doi: 10.1111/j.1365-2966.2011.19237.x
  * Deason et al. (2017) Deason, A. J., Belokurov, V., Koposov, S. E., et al. 2017, MNRAS, 470, 1259, doi: 10.1093/mnras/stx1301
  * Deason et al. (2018) Deason, A. J., Belokurov, V., Koposov, S. E., & Lancaster, L. 2018, ApJ, 862, L1, doi: 10.3847/2041-8213/aad0ee
  * Deason et al. (2012) Deason, A. J., Belokurov, V., Evans, N. W., et al. 2012, MNRAS, 425, 2840, doi: 10.1111/j.1365-2966.2012.21639.x
  * Deason et al. (2021) Deason, A. J., Erkal, D., Belokurov, V., et al. 2021, MNRAS, 501, 5964, doi: 10.1093/mnras/staa3984
  * Delahaye et al. (2019) Delahaye, D., Chaimatanan, S., & Mongeau, M. 2019, in International Series in Operations Research & Management Science (ISOR), Vol. 272, Handbook of Metaheuristics, ed. M. Gendreau & J.-Y. Potvin (Springer), 1–35.ISBN 978–3–319–91085–7, doi: 10.1007/978-3-319-91086-4_1
  * Dempster et al. (1977) Dempster, A. P., Laird, N. M., & Rubin, D. B. 1977, J. R. Stat. Soc. B, 39, 1, doi: https://doi.org/10.1111/j.2517-6161.1977.tb01600.x
  * Donlon et al. (2020) Donlon, Thomas, I., Newberg, H. J., Sanderson, R., & Widrow, L. M. 2020, ApJ, 902, 119, doi: 10.3847/1538-4357/abb5f6
  * Erkal et al. (2021) Erkal, D., Deason, A. J., Belokurov, V., et al. 2021, MNRAS, 506, 2677, doi: 10.1093/mnras/stab1828
  * Fermani & Schönrich (2013) Fermani, F., & Schönrich, R. 2013, MNRAS, 430, 1294, doi: 10.1093/mnras/sts703
  * Gaia Collaboration (2020) Gaia Collaboration. 2020, Gaia Source Catalogue EDR3, IPAC, doi: 10.26131/IRSA541
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Hayes & Philip (1979) Hayes, D. S., & Philip, A. G. D. 1979, PASP, 91, 71, doi: 10.1086/130444
  * Helmi (2020) Helmi, A. 2020, ARA&A, 58, 205, doi: 10.1146/annurev-astro-032620-021917
  * Koppelman et al. (2018) Koppelman, H., Helmi, A., & Veljanoski, J. 2018, ApJ, 860, L11, doi: 10.3847/2041-8213/aac882
  * Lancaster et al. (2019) Lancaster, L., Koposov, S. E., Belokurov, V., Evans, N. W., & Deason, A. J. 2019, MNRAS, 486, 378, doi: 10.1093/mnras/stz853
  * Leavitt & Pickering (1912) Leavitt, H. S., & Pickering, E. C. 1912, Harvard College Observatory Circular, 173, 1
  * Lee et al. (2008a) Lee, Y. S., Beers, T. C., Sivarani, T., et al. 2008a, AJ, 136, 2022, doi: 10.1088/0004-6256/136/5/2022
  * Lee et al. (2008b) —. 2008b, AJ, 136, 2050, doi: 10.1088/0004-6256/136/5/2050
  * Li et al. (2022) Li, T. S., Ji, A. P., Pace, A. B., et al. 2022, ApJ, 928, 30, doi: 10.3847/1538-4357/ac46d3
  * Lindegren (2018) Lindegren, L. 2018.  http://www.rssd.esa.int/doc_fetch.php?id=3757412
  * Malhan et al. (2018) Malhan, K., Ibata, R. A., & Martin, N. F. 2018, MNRAS, 481, 3442, doi: 10.1093/mnras/sty2474
  * Martin et al. (2022) Martin, N. F., Venn, K. A., Aguado, D. S., et al. 2022, Nature, 601, 45, doi: 10.1038/s41586-021-04162-2
  * Myeong et al. (2018) Myeong, G. C., Evans, N. W., Belokurov, V., Sanders, J. L., & Koposov, S. E. 2018, ApJ, 856, L26, doi: 10.3847/2041-8213/aab613
  * Peñarrubia & Petersen (2021) Peñarrubia, J., & Petersen, M. S. 2021, MNRAS, 508, L26, doi: 10.1093/mnrasl/slab090
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Petersen & Peñarrubia (2021) Petersen, M. S., & Peñarrubia, J. 2021, Nature Astronomy, 5, 251, doi: 10.1038/s41550-020-01254-3
  * Pier (1983) Pier, J. R. 1983, ApJS, 53, 791, doi: 10.1086/190910
  * Preston et al. (1991) Preston, G. W., Shectman, S. A., & Beers, T. C. 1991, ApJ, 375, 121, doi: 10.1086/170175
  * Rockosi et al. (2022) Rockosi, C. M., Lee, Y. S., Morrison, H. L., et al. 2022, The Astrophysical Journal Supplement Series, 259, 60, doi: 10.3847/1538-4365/ac5323
  * RStudio Team (2022) RStudio Team. 2022, RStudio: Integrated Development Environment for R, RStudio, PBC, Boston, MA.  http://www.rstudio.com/
  * Santucci et al. (2015a) Santucci, R. M., Placco, V. M., Rossi, S., et al. 2015a, ApJ, 801, 116, doi: 10.1088/0004-637X/801/2/116
  * Santucci et al. (2015b) Santucci, R. M., Beers, T. C., Placco, V. M., et al. 2015b, ApJ, 813, L16, doi: 10.1088/2041-8205/813/1/L16
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772
  * Sersic (1968) Sersic, J. L. 1968, Atlas de Galaxias Australes (Universidad Nacional de Cordoba: Observatorio Astronomico)
  * Shapley (1916) Shapley, H. 1916, ApJ, 43, 217, doi: 10.1086/142246
  * Shull & Danforth (2019) Shull, J. M., & Danforth, C. W. 2019, ApJ, 882, 180, doi: 10.3847/1538-4357/ab357d
  * Sirko et al. (2004) Sirko, E., Goodman, J., Knapp, G. R., et al. 2004, AJ, 127, 899, doi: 10.1086/381483
  * Starkenburg et al. (2019) Starkenburg, E., Youakim, K., Martin, N., et al. 2019, MNRAS, 490, 5757, doi: 10.1093/mnras/stz2935
  * Taylor (2005) Taylor, M. B. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert, 29
  * Thomas et al. (2018) Thomas, G. F., McConnachie, A. W., Ibata, R. A., et al. 2018, MNRAS, 481, 5223, doi: 10.1093/mnras/sty2604
  * Utkin & Dambis (2020) Utkin, N. D., & Dambis, A. K. 2020, MNRAS, 499, 1058, doi: 10.1093/mnras/staa2819
  * Valcarce & Catelan (2008) Valcarce, A. A. R., & Catelan, M. 2008, A&A, 487, 185, doi: 10.1051/0004-6361:20078231
  * Vasiliev & Baumgardt (2021) Vasiliev, E., & Baumgardt, H. 2021, MNRAS, 505, 5978, doi: 10.1093/mnras/stab1475
  * Vickers et al. (2012) Vickers, J. J., Grebel, E. K., & Huxor, A. P. 2012, AJ, 143, 86, doi: 10.1088/0004-6256/143/4/86
  * Wang & Chen (2019) Wang, S., & Chen, X. 2019, ApJ, 877, 116, doi: 10.3847/1538-4357/ab1c61
  * Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9, doi: 10.1051/aas:2000332
  * Whitten et al. (2019) Whitten, D. D., Beers, T. C., Placco, V. M., et al. 2019, ApJ, 884, 67, doi: 10.3847/1538-4357/ab4269
  * Wu et al. (2022) Wu, W., Zhao, G., Xue, X.-X., Bird, S. A., & Yang, C. 2022, ApJ, 924, 23, doi: 10.3847/1538-4357/ac31ac
  * Xue et al. (2008) Xue, X. X., Rix, H. W., Zhao, G., et al. 2008, ApJ, 684, 1143, doi: 10.1086/589500
  * Xue et al. (2011) Xue, X.-X., Rix, H.-W., Yanny, B., et al. 2011, ApJ, 738, 79, doi: 10.1088/0004-637X/738/1/79
  * Yanny et al. (2009) Yanny, B., Rockosi, C., Newberg, H. J., et al. 2009, AJ, 137, 4377, doi: 10.1088/0004-6256/137/5/4377
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513
  * Yuan et al. (2020) Yuan, Z., Chang, J., Beers, T. C., & Huang, Y. 2020, ApJ, 898, L37, doi: 10.3847/2041-8213/aba49f
  * Yuan et al. (2019) Yuan, Z., Smith, M. C., Xue, X.-X., et al. 2019, ApJ, 881, 164, doi: 10.3847/1538-4357/ab2e09
  * Yuan et al. (2022) Yuan, Z., Malhan, K., Sestito, F., et al. 2022, ApJ, 930, 103, doi: 10.3847/1538-4357/ac616f

##  Appendix A Fitting the horizontal branch with [Fe/H]

For the analysis of dependence of metallicity, we have a smaller sample of
stars than the one used for the calibrations, about 40 stars from the original
sample were available in the SSPP data. The same procedure was done with both
the pure photometric sample and this one cross-matched with SSPP.

Each solution evaluation (considering different input variables) was taken in
a period of approximately 10 min, which is enough for the convergence of
several solutions, as TuringBot needs less time to converge than other similar
softwares (Ashok et al., 2020). The hardware involved in this process is
highly important for the convergence time. In our case, the program was
executed in a computer with an AMD Ryzen 7 2700 processor, with 16 threads,
all used at once. Tests were also made by running the program longer and no
significant improvement was observed.

Figure 6: Error comparison between fits for some combinations of magnitudes,
colors, and [Fe/H] from TuringBot.

Figure 6 shows how the error decreases with the increase of the complexity of
the functions (see Section 3.2 for the definition of “complexity”). Some fits
coincide since the program can create colors from magnitudes, if it is better
than the magnitudes alone, and not use all the variables provided. This is the
reason why the line for all the parameters is not visible. We can see that,
for equations up to complexity 25, there is no advantage of adding the
metallicity information as we can achieve similar errors using only
magnitudes.

##  Appendix B Globular cluster distances

Figure 7 presents a boxplot for the distances obtained with both calibrations
and the one from D11 compared to those provided by Vasiliev & Baumgardt
(2021). The size of each box represent the spread (25th and 75th quartiles)
and the center line indicates the median value. The number of stars in each
cluster is displayed above their NGC identifier.

The overall results indicate that Eq. 2 is in general more accurate than D11’s
relation. This conclusion is supported by better distance predictions where 10
of the 11 globular clusters were better constrained with Eq. 2. The
calibration with color $(u-r)_{0}$ revealed to be more accurate than using the
relation provided in this work for color $(g-r)_{0}$, as nine clusters present
lower dispersion and eight clusters show medians close to zero when using the
former color. NGC2416 was the only cluster where the color $(g-r)_{0}$
presented a better performance, however, this cluster was represented by only
9 members and therefore this result may be due to sub-sampling. The sample
considered has a bias to closer cluster, which is apparent by huge distance
gap between NGC2416, the farthest cluster, and NGC5466, the penultimate.
Therefore, we cannot ascertain whether or not the accuracy of distances
estimated using color $(u-r)_{0}$ varies with the distance.

Figure 7: Boxplot for relative distance difference using each calibration
presented in this work (Eq. 1 and Eq. 2) and that from D11. $D_{V21}$ is the
distance provided by Vasiliev & Baumgardt (2021), inferred from Gaia EDR3
data. The number of data points in each cluster is displayed at the bottom of
the panel.


# Hierarchical hub-filament structures and gas inflows on galaxy-cloud scales

J. W. Zhou  Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121
Bonn, Germany<EMAIL_ADDRESS>Timothy A. Davis  Cardiff Hub for
Astrophysics Research and Technology, School of Physics and Astronomy, Cardiff
University, Queens Buildings, Cardiff CF24 3AA, UK

###### Abstract

We investigated the kinematics and dynamics of gas structures on galaxy-cloud
scales in two spiral galaxies NGC5236 (M83) and NGC4321 (M100) using CO
(2$-$1) line. We utilized the FILFINDER algorithm on integrated intensity maps
for the identification of filaments in two galaxies. Clear fluctuations in
velocity and density were observed along these filaments, enabling the fitting
of velocity gradients around intensity peaks. The variations in velocity
gradient across different scales suggest a gradual and consistent increase in
velocity gradient from large to small scales, indicative of gravitational
collapse, something also revealed by the correlation between velocity
dispersion and column density of gas structures. Gas structures at different
scales in the galaxy may be organized into hierarchical systems through
gravitational coupling. All the features of gas kinematics on galaxy-cloud
scale are very similar to that on cloud-clump and clump-core scales studied in
previous works. Thus, the interstellar medium from galaxy to dense core scales
presents multi-scale/hierarchical hub-filament structures. Like dense core as
the hub in clump, clump as the hub in molecular cloud, now we verify that
cloud or cloud complex can be the hub in spiral galaxies. Although the scaling
relations and the measured velocity gradients support the gravitational
collapse of gas structures on galaxy-cloud scales, the collapse is much slower
than a pure free-fall gravitational collapse.

###### keywords:

Submillimeter: ISM; ISM: clouds; ISM: kinematics and dynamics; galaxies: ISM;
galaxies: structure; galaxies: star formation; techniques: image processing

##  1 Introduction

Measuring the dynamical coupling between density enhancements in giant
molecular clouds and gas motion of their local environment gives us a
perspective to understand the formation of hierarchical structures in high-
mass star formation regions McKee2007-45; Motte2018-56; Henshaw2020-4. High-
resolution observations of high-mass star-forming regions reveal the
structured arrangement of density enhancements within filamentary gas
networks, notably in hub-filament systems. Within these systems, converging
flows channel material into the central hub through the interconnected
filaments Peretto2013; Henshaw2014; Zhang2015; Liu2016; Yuan2018; Lu2018;
Issac2019; Dewangan2020; Liu2022-511; Zhou2022-514; Zhou2023-676. In
observations, gas inflow will naturally produce velocity gradients along
filaments Kirk2013; Liu2016; Yuan2018; Williams2018-613; Chen2019-875;
Chen2020-891; Pillai2020-4; Zhou2022-514; Zhou2023-676.

Henshaw2020-4 identified widespread velocity fluctuations spanning various
spatial scales and physical contexts within galaxies. They observed
oscillatory gas movements with wavelengths ranging from 0.3 to 400 parsecs,
intricately linked to regularly spaced density enhancements. These
enhancements are likely a result of gravitational instabilities
Henshaw2016-463; Elmegreen2018-863. Furthermore, the spatial correlation
between density enhancements and velocity gradient extrema may indicate
convergent motion due to gravitational collapse Hacar2011-533; Hacar2016-587;
Clarke2016-458; Misugi2019-881; Zhou2022-514; Zhou2023-676.

Zhou2022-514 investigated the physical properties and evolution of hub-
filament systems in $\sim$ 140 protoclusters, utilizing spectral lines
observed in the ATOMS (ALMA Three-millimeter Observations of Massive Star-
forming regions) survey Liu2020. Their findings indicate the presence of hub-
filament structures across a range of scales, spanning from 0.1 parsec to
several parsecs, in diverse Galactic environments. Additionally, slender
structures resembling filaments, such as spiral arms, have been identified at
scales below 0.1 parsecs in the vicinity of high-mass protostars Liu2015-804;
Maud2017-467; Izquierdo2018-478; Chen2020-4; Sanhueza2021-915. As proposed by
Zhou2022-514, self-similar hub-filament structures and filamentary accretion
seem to exist across scales, ranging from several thousand astronomical units
to several parsecs, within high-mass star-forming regions. This paradigm of
hierarchical/multi-scale hub-filament structures was generalized from clump-
core scale to cloud-clump scale in Zhou2023-676. Hierarchical collapse and
hub-filaments structures feeding the central regions are also described in
previous works, see Motte2018-56; Vazquez2019-490; Kumar2020-642 and
references therein.

Kinematically, the results in Zhou2023-676 also reveal the multi-scale hub-
filament structures in the G333 molecular cloud complex. The G333 complex
exhibits prevalent kinematic characteristics consistent with hub-filament
systems. Notably, the intensity peaks, acting as hubs, correlate with
converging velocities, indicating that the surrounding gas flows are
converging to dense structures. Specifically, there is a discernible increase
in velocity gradients at smaller scales. The filaments in the Large APEX sub-
Millimeter Array (LAsMA) and the Atacama Large Millimeter/submillimeter Array
(ALMA) observations show clear velocity gradients. The velocity gradients
fitted using the LAsMA and ALMA data exhibit consistent agreement over the
range of scales covered by ALMA observations in the ATOMS survey ($\textless$
5 pc). In Zhou2023-676, larger scale gas motions were investigated (the
longest filament $\sim$50 pc), yet similar results were obtained compared to
small-scale ALMA observations. Interestingly, the variations in velocity
gradients measured at distinct scales—small scales ($\textless$ 1 pc), medium
scales ($\sim$ 1-7.5 pc), and large scales ($\textgreater$ 7.5 pc)—align with
expectations from gravitational free-fall with central masses of $\sim$ 500
M⊙, $\sim$ 5000 M⊙ and $\sim$ 50000 M⊙ for the respective scales. This implies
that sustaining velocity gradients on larger scales necessitates higher
masses. The association of higher masses with larger scales suggests that the
inflow on a larger scale is driven by the larger-scale structure which may be
the gravitational clustering of smaller-scale structures. This observation
aligns with the hierarchical nature commonly found in molecular clouds and the
gas inflow from large to small scales. The large-scale gas inflow is likely
driven by gravity, indicating a state of global gravitational collapse within
the molecular clouds in the G333 complex. This supports the argument that
these molecular clouds serve as cloud-scale hub-filament structures. The
change in velocity gradients with the scale in the G333 complex indicates that
the morphology of the velocity field in position-position-velocity (PPV) space
resembles a ”funnel” structure. The funnel structure can be elucidated as the
acceleration of material converging towards the central hub and the
gravitational contraction of star-forming clouds or clumps. Large-scale
velocity gradients always associate with many intensity peaks, and the larger-
scale inflow is driven by the larger-scale structure, indicating that the
clustering of smaller-scale gravitational structures locally can serve as the
gravitational center on a larger scale. In a way, the funnel structure
provides insight into the gravitational potential well shaped by this
clustering.

In previous work, we have investigated gas dynamics and kinematics in
molecular clouds from cloud to core scales. The main task of this work is to
generalize the above physical pictures to galaxy-cloud scales, now the
smallest unit is the molecular cloud itself. In Zhou2022-514, dense core is
the hub in clump, in Zhou2023-676, clump is the hub in molecular cloud, here
we treat molecular cloud as the hub in galaxy. Self-similar or
hierarchical/multi-scale hub-filament structures and filamentary accretion
feeding the central regions exist from molecular cloud to dense core scales
Motte2018-56; Vazquez2019-490; Kumar2020-642; Zhou2022-514; Zhou2023-676, this
picture will be extended to galaxy-cloud scales in this work.

##  2 Data and target

We select two face-on spiral galaxies NGC5236 (M83) and NGC4321 (M100) from
the PHANGS-ALMA survey. We use the combined 12m+7m+TP PHANGS-ALMA CO (2$-$1)
data cubes to investigate gas kinematics and dynamics in the two galaxies,
which have a spectral resolution of 2.5 km s-1 and angular resolutions
$\sim$2.1 ′′ and $\sim$1.7 ′′, corresponding to linear resolutions $\sim$51 pc
and $\sim$123 pc for NGC5236 and NGC4321 at the distances 4.89 Mpc and 15.21
Mpc Leroy2021-257; Anand2021-501, respectively. The field-of-view (FOV) of
NGC5236 and NGC4321 in the ALMA observations are $\sim$ 10.5 Kpc $\times$ 10.8
Kpc and $\sim$ 17.1 Kpc $\times$ 15.5 Kpc, respectively. An overview of the
PHANGS-ALMA survey’s science goals, sample selection, observation strategy,
and data products is described in Leroy2021-257. A detailed description of
data calibration, imaging, and product creation procedures is presented in
Leroy2021-255. The PHANGS-ALMA CO (2$-$1) data cubes and other data products
(such as Moment maps) are available from the PHANGS team website
111https://sites.google.com/view/phangs/home. NGC5236 and NGC4321 were
selected for the following reasons:

1\. They are face-on galaxies. The inclination angles of NGC5236 and NGC4321
are 24o and 38.5o Lang2020-897, respectively. Thus, we can ignore the
projection effect in velocity gradient fitting in Sec.3.3;

2\. They have strong CO (2$-$1) emission, presenting large-scale continuous
structures (strong spiral arms), thus we can trace galaxy-scale gas motions;

3\. NGC5236 has similar morphology to NGC4321, but it is about 3 times closer
than NGC4321, thus we can resolve more detailed structures.

##  3 Results

###  3.1 Velocity component

Figure 1:  An overview of gas kinematics in NGC4321 and NGC5236. (a) and (d)
The velocity distribution of the entire galaxies in PPV space decomposed by
GAUSSPY+; (b) and (e) The main filaments of NGC4321 and NGC5236 identified by
FILFINDER algorithm overlaid on the Moment 0 maps; (c) and (f) The integrated
intensity distributions extracted from the Moment 0 maps along the main
filaments in panels (b) and (e); (c2) and (f2) The velocity distributions
extracted from the Moment 1 maps along the main filaments in panels (b) and
(e); (c1) and (f1) The velocity residual distributions along the main
filaments in panels (b) and (e) after subtracting the large-scale velocity
gradients due to the galaxy rotation.

Before studying gas kinematics, we need to ascertain the distribution of gas
components in the galaxy. We utilized the fully automated Gaussian decomposer
GAUSSPY+ Lindner2015-149; Riener2019-628, designed to break down intricate
spectra of molecular lines into multiple Gaussian components. The parameter
settings for the decomposition align with those employed by Zhou2023-676. For
the CO (2$-$1) data cube of NGC4321, 50071 spectra were fitted by GAUSSPY+,
98% spectra are single component, 2% spectra have two velocity components.
323059 spectra were fitted for NGC5236, 87% spectra show single component,
11.5% spectra have two components. The reduced $\chi^{2}$ values output in the
fitting of NGC4321 and NGC5236 are around 1, which means that both of them get
a good fit. Generally, multiple velocity components ($>$1 velocity component)
mainly concentrate on the central region of the galaxy. In this work, we
mainly focus on the gas kinematics on the spiral arms. Therefore, the analysis
can be carried out based on the Moment maps.

In Fig.1, we can see ubiquitous velocity fluctuations, which may be attributed
to gravitational instability, as discussed below. In the central region of the
galaxy, there is almost no velocity fluctuations, clear velocity fluctuations
are mainly distributed on the spiral arms, revealing the interconnection
between velocity fluctuations and star formation activities.

###  3.2 Scaling relations

####  3.2.1 Dendrogram

Figure 2: Dendrogram structures. (a) Leaves; (b) Branches. Red and orange eclipses represent type1 and type3 structures, respectively. Figure 3: Same as Fig.2, but for NGC5236. Table 1: Dendrogram structures in NGC 4321 and NGC 5236. | leaves | branches | type1 | type3 | total
---|---|---|---|---|---
NGC 4321 | 583 | 270 | 828 | 25 | 853
NGC 5236 | 1722 | 859 | 2339 | 242 | 2581

The issues of identifying structures in a PPV cube using the dendrogram
algorithm are discussed in Zhou2024-682-128. Consequently, we employed the
same approach as outlined in Zhou2024-682-128 and Zhou2024-682-173 to identify
dense gas structures. We conducted a direct identification of hierarchical
(sub-)structures based on the 2D integrated intensity (Moment 0) map of CO
(2$-$1) emission. Subsequently, we extracted the average spectrum of each
identified structure to delve into its velocity components and gas kinematics.
All the retained structures on the strictly masked Moment 0 map are reliable,
we therefore only require the smallest area of the identified structure larger
than 1 beam and do not set other parameters in the algorithm to reduce the
dependence of the identification on the parameter settings. Dendrogram
algorithm decomposes the intensity data (Moment 0 map) into hierarchical
structures. As illustrated in Fig.1 of Zhou2024-682-128, all identified
structures can be divided into three categories, i.e. i-leaves (isolated
leaves), c-leaves (clustered leaves) and branches. In this work, we merged
c-leaves and i-leaves. Finally, there are two kinds of structures, i.e. leaf
and branch structures. In Fig.2 and Fig.3, the structures identified by the
dendrogram algorithm exhibit a strong correspondence with the background
integrated intensity maps.

The algorithm characterizes the morphology of each structure by approximating
it as an ellipse. Within the dendrogram, the rms sizes (second moments) of the
intensity distribution along the two spatial dimensions define the long and
short axes of the ellipse, denoted as $a$ and $b$. As described in
Zhou2024-682-128, a smaller ellipse is obtained with $a$ and $b$,
necessitating a multiplication by a factor of two to appropriately enlarge the
ellipse. Then the effective physical radius of an ellipse is $R\rm_{eff}$
=$\sqrt{2a\times 2b}*D$, where $D$ is the distance of the galaxy. For a
structure with the area $A$, the total integrated intensity of the structure
is $I_{\rm CO}$, then the mass of the structure can be calculated by

$M=\alpha^{2-1}_{\rm CO}\times I_{\rm CO}\times A,$ (1)

here $\alpha^{2-1}_{\rm CO}\approx 6.7\rm M_{\odot}\rm pc^{-2}(\rm
K*kms^{-1})^{-1}$ Leroy2022-927.

The large-scale velocity gradients due to the galaxy rotation will contribute
to the non-thermal velocity dispersion. Before extracting the average spectra
of the identified structures, we subtracted the large-scale velocity gradients
in PPV data cube based on the constructed gas dynamical model in Sec.3.3.3.
According to their averaged spectra, the structures with absorption features
were eliminated firstly. Following Zhou2024-682-173, we only consider type1
(single velocity component) and type3 (blended velocity components) structures
in this work.

As detailed in Zhou2024-682-128, two screening criteria were applied for
structure refinement: (1) Removal of repetitive branch structures; and (2)
Exclusion of branch structures exhibiting complex morphology. Finally, there
are 853 and 2581 retained structures for NGC4321 and NGC5236, respectively, as
listed in Table 1 and marked in Fig.2 and Fig.3 by ellipses. NGC5236 is
$\sim$3 times nearer than NGC4321, the total number of the identified
structures in NGC5236 is also $\sim$3 times more than that in NGC4321.

####  3.2.2 Scaling relations

Figure 4: Scaling relations of leaf and branch structures of NGC 4321 and NGC
5236. (a)-(c) and (a1)-(c1) $\sigma-N*R$; (d)-(f) and (d1)-(f1) $\sigma-N$;
(g)-(i) and (g1)-(i1) $\sigma-R$. $\sigma$, $R$ and $N$ are the velocity
dispersion, effective radius, and column density of each structure,
respectively. Here, “P” represents the Pearson coefficient. Dashed vertical
line marks the boundary of two fittings.

The scaling relations provide insights into the physical states of the
structures. Although the identified structures in NGC5236 are much more than
that in NGC4321, the scaling relations in Fig.4 are comparable. Due to the
further distance, some faint structures were filtered out in NGC4321, which
can be seen in NGC5236. In Fig.4, the faint or low-density structures produce
dispersive tails, similar to the Fig.14 of Zhou2024-682-128 for the G333
molecular cloud complex in the Milky Way.

In Fig.4(g)-(i), it appears that only branch structures demonstrate a
discernible correlation between velocity dispersion and scale. Clumps or cores
characterized by high column density exhibit greater velocity dispersion
relative to those with lower column density, a phenomenon attributed to gas
motions associated with gravitational collapse Ballesteros2011-411;
Traficante2018-473; Ballesteros2018-479; Li2023-949; Zhou2024-682-128. As
depicted in Fig.4(d)-(f), the positive correlations observed between velocity
dispersion and column density suggest a gravitational origin for the velocity
dispersion. In the context of pure free-fall, an anticipated $\sigma-N*R$
relation would exhibit a slope of 0.5. For a more convenient comparison with
$\sigma-R$ and $\sigma-N$ relations, we transform the Heyer relation,
$\sigma/R^{0.5}\propto N^{0.5}$, to the form $\sigma\propto(R*N)^{0.5}$
(Ballesteros2011-411, Eq. 3 in). Both relations should ideally have a slope of
0.5. However, in Fig.4(a)-(c), the slopes of the $\sigma-N*R$ relations are
notably shallower than 0.5, suggesting a deceleration from the expected
behavior of pure free-fall gravitational collapse.

###  3.3 Velocity gradient

The gravitational collapse of the structures has been revealed by the scaling
relations. This section gives a more detailed discussion of the gravitational
collapse from the perspective of the velocity gradient.

####  3.3.1 Identification of filaments

The FILFINDER algorithm Koch2015 was used to identify and characterize
filaments in the galaxy. Following the method described in Zhou2022-514 and
Zhou2023-676, we used the Moment 0 maps of CO (2$-$1) emission from the data
products of the PHANGS-ALMA survey to identify filaments. Fig.1(b) and (e)
display the skeletons of the identified filaments superimposed on the Moment 0
maps. Notably, these skeletons exhibit a high degree of agreement with the gas
structures traced by CO (2$-$1) emission. In this work, we only consider the
filaments along the spiral arms (the main filaments), where have strong CO
(2$-$1) emission, which allows us to study large-scale continuous gas motions
in the galaxy.

Almost all algorithms have the issue of parameter settings, thus we try not to
discuss the identified filament itself in this work, such as its length and
width. Actually, we mainly regard the FILFINDER algorithm as a tool to draw
lines, these lines can make connections between discrete dense structures, or
put dense structures in gas networks. Then by extracting the velocity
information along these lines, we can study how the surrounding gas converges
to dense gas structures (gravitational centers). In Zhou2022-514 and
Zhou2023-676, the filaments traverse multiple local hub-filament structures,
characterized by fluctuations in velocity and density along the filaments.
These local dense structures, acting as gravitational centers, facilitate the
accretion of surrounding diffuse gas, culminating in the formation of local
hub-filament structures. In the PPV space, a hub-filament structure can be
represented as a funnel structure, as demonstrated in Fig.9 of Zhou2023-676.
The gradient of the funnel profile serves as an indicator of the gravitational
field’s strength. Consequently, examining the local velocity gradients along
the filaments provides insights into the strength of local gravitational
fields.

In Fig.1, although the velocity fluctuations at small scales can be seen
everywhere, they are hidden in global large-scale velocity gradients due to
the galaxy rotation. To derive the local velocity fluctuations, we must first
subtract the large-scale velocity gradients. We take two methods to solve this
issue.

####  3.3.2 Linear fitting

Figure 5: A segment of NGC5236’s main filament, used to demonstrate the
velocity gradient fitting. (a) Linear fitting of the large-scale velocity
gradient due to the bulk motion; (b) Velocity field of the segment in panel
(a) after subtracting the fitted large-scale velocity gradient. Local velocity
gradients are fitted in the ranges defined by the red vertical dashed lines in
panel (c), and straight lines show the linear fitting results. The color-
coding of straight lines are meaningless, but it is convenient for the reader;
(c) Blue and red dotted lines show the normalized velocity and intensity,
respectively. Orange box marks a gravitational coupling structure.

Fig.1(c2)&(c) and (f2)&(f) show the intensity-weighted velocity (Moment 1) and
integrated intensity (Moment 0) of CO (2$-$1) line emission along the main
filaments with intense velocity and density fluctuations. A segment of
NGC5236’s main filament was used to demonstrate the velocity gradient fitting
in Fig.5. Here, we eliminated the large-scale velocity gradient by employing a
simple linear function for modeling. Subsequently, we subtracted this model
from the velocity field obtained from the Moment 1 map along the main
filament. We divided the main filaments into many segments to ensure that each
segment is as straight as possible, thus increasing the accuracy of the linear
fit, as shown in Fig.5(a). The fitted straight line can be the baseline, then
we subtracted the baseline for each segment to leave only the residual local
velocity fluctuations. The good correspondence between velocity gradients and
intensity peaks in Fig.5(b) indicates that the surrounding gases are
converging to the center of the gravitational potential well. For each segment
of the main filament, we firstly estimated the global velocity gradients
between velocity peaks and valleys at two sides of the intensity peaks
(representing gravitational centers), and ignored local velocity fluctuations.
Then we also derived additional velocity gradients over smaller distances
around the local intensity peaks, as illustrated in Fig.5(b). Generally,
large-scale velocity gradients are associated with multiple intensity peaks.
Given that NGC4321 and NGC5236 are good face-on galaxies, here we ignored the
projection effect in velocity gradient fitting.

####  3.3.3 Gas dynamical model

Figure 6: (a) The velocity field (Moment 1) of NGC 4321 in observation; (b)
Gas dynamical model created by the Kinematic Molecular Simulation (KinMS)
package; (c) Velocity fluctuations after removing the created model.

As a comparison, we also removed the bulk motion or the galaxy rotation by
creating gas dynamical models using the Kinematic Molecular Simulation (KinMS)
package of Davis2013-429 based on the Moment-1 map, as shown in Fig.6(b). Then
we used the Moment 1 map to subtract the created model, and obtained the
velocity residual map. Finally, we extracted the velocity fluctuations from
the velocity residual map rather than the Moment 1 map. The subsequent
velocity gradient fitting is the same with Sec.3.3.2.

####  3.3.4 Statistical analysis of the fitted velocity gradients

Figure 7: (a) The correlation between the fitted velocity gradient and the
scale; (b) The mass distribution of the identified structures. Figure 8:
Statistical analysis of all fitted velocity gradients. (a) Velocity gradient
vs. the length. The color lines show the freefall velocity gradients for
comparison. For the freefall model, red, magenta, green, cyan, blue and black
lines denote masses of 105 M⊙, 105.5 M⊙, 106 M⊙ , 106.5 M⊙, 107 M⊙ and 107.5
M⊙, respectively. Panels (b), (c), and (d): Zoomed maps with lengths
$\textless$ 500 pc (small scale), $\sim$ 500 – 1500 pc (medium scale), and
$\textgreater$ 1500 pc (large scale) in panel (a).

In Sec.3.3.2 and Sec.3.3.3, the galaxy centers were excluded in the fitting.
As shown in Fig.7, two methods give consistent fitting results. In Fig.8, we
can find the same gas kinematic modes presented in Zhou2023-676 and
Zhou2024-682-173. The variations in velocity gradients at both small and large
scales (with a boundary around 500 pc) align with expectations from
gravitational free-fall, with central masses of $\sim$ 105–106.5 M⊙ and $\sim$
106–107.5 M⊙. This implies that sustaining velocity gradients on larger scales
requires correspondingly larger masses, and larger masses imply larger scales,
suggesting that the larger-scale inflow is driven by the larger-scale
structure which may arise from the gravitational clustering of smaller-scale
structures, in harmony with the presence of hierarchical or multi-scale hub-
filament structures within the galaxy and the gas inflows from large to small
scales. In the orange box marked in Fig.5(b) and (c), multiple peaks are
coupled together to form a gravitational potential well on larger scale, and
each peak itself is also a local gravitational potential well.

In Fig.7(a), almost all measured velocity gradients can be fitted in the mass
range $\sim$ 105–107 M⊙, which is consistent with the mass distribution of the
identified structures shown in Fig.7(b), indicating that local dense
structures and their complex as gravitational centers will accrete the
surrounding diffuse gas and then produce the observed velocity gradients at
different scales.

####  3.3.5 Deviation of the free-fall model

We can do a deeper analysis for the fitted velocity gradients with a simple
model. Assuming free-fall,

$\nabla v=-\frac{d}{dR}\sqrt{\frac{2GM}{R}}=\sqrt{\frac{GM}{2R^{3}}},$ (2)

in this equation, velocity gradient is more sensitive to scale. Thus we only
fit the correlation between velocity gradient and scale in Fig.7(a). In
equation.2, $\nabla v\propto R^{-1.5}$, however, the linear fitting in
Fig.7(a) gives $\nabla v\propto R^{-0.8}$, thus $\sim$2 times smaller than the
slope of the free-fall model, also indicating the slowing down of a pure free-
fall gravitational collapse presented in Sec.3.2.2.

##  4 Discussion

Gas kinematics on galaxy-cloud scales clearly deviate from the free-fall
model. The deviation may come from measurement or calculation biases. We only
measured the projection of the realistic velocity vector and hence
acceleration, thus tends to underestimate the observed acceleration. Moreover,
we only considered the molecular line CO (2$-$1), which only traces part of
the baryonic matters in the galaxy, thus underestimated the mass of the
gravitational center.

The shallower slope in Fig.7(a) may have physical correlation with the flat
rotation curve of the galaxy. Gas motions on galaxy-cloud scales may still
couple with the galactic potential, a comprehensive model should account for
the interplay between motions within the galactic potential and the self-
gravitational potential of the cloud. Dobbs2013-432; Meidt2018-854;
Meidt2020-892; Utreras2020-892. In the model proposed by Meidt2020-892, the
transition to free-fall collapse can only happen once the gas has completely
decoupled from the galactic potential. However, it is not clear down to which
spatial scales gas motions remain dynamically relevant to galactic potential
or start decoupling with galactic potential.

Tidal forces have been advocated in previous works to restrict or trigger star
formation. Previous studies by Ballesteros2009-393; Ballesteros2009-395 have
shown the impact of tidal forces arising from an effective Galactic potential
on molecular clouds. These forces have the potential to either compress or
disrupt molecular clouds, influencing the overall star formation efficiency.
Thilliez2014-31 examined the stability of molecular clouds within the Large
Magellanic Cloud (LMC) against shear and the galactic tide. However, their
findings indicate that star formation in the LMC is not impeded by either
tidal or shear instability. Moreover, in Ramirez2022-515, tidal stresses from
neighbouring molecular cloud complexes may increase interstellar turbulence,
rather than the galactic potential.

Gravity is a long-range force. A local dense structure evolves under its self-
gravity, but as a gravitational center, its influence can also affect
neighboring structures. At the same time, it also experiences the
gravitational pull from other nearby sources. The tidal and gravitational
fields are mutually interdependent. As described above, the
hierarchical/multi-scale gravitational coupling of gas structures means the
extensive tidal interactions in the galaxy. Whether the galactic potential has
effect on molecular clouds and their complexes or not, molecular clouds should
be affected by the cumulative tidal interactions of many nearby materials,
which may prevent gravitational collapse and growth of instabilities or star
formation in the cloud. Due to the diffuse and complex morphology of matter
distribution in the galaxy, a complete tidal calculation would be very
complex. One should derive the gravitational potential distribution from the
observed density distribution and then calculate the tidal field according to
the gravitational potential distribution, as presented in Li2024-528. Diverse
manifestations of gravitational effects on gas within molecular clouds were
unveiled in Li2024-528: Dense regions experience gravitational collapse, while
the surrounding gas is subject to significant tidal forces, suppressing
fragmentation. This gas, influenced by extensive tides, is directed towards
the dense regions, serving as fuel for star formation. The spatial
distribution of regions experiencing varying tidal influences elucidates the
observed hierarchical and localized pattern of star formation. Similar
mechanisms may also exist on galaxy-cloud scales, we will discuss this topic
in detail in future work.

Figure 9: The average number density distribution of the identified
structures.

In addition to the factors discussed above, magnetic fields may have a
significant impact on the gas kinematics of molecular clouds as suggested by
simulations Kim2021-911; Seifried2020-497; Ibanez2022-925; Ganguly2023-525 and
observations Crutcher2010-725; Li2011-479; Crutcher2012-50; Stephens2022-926;
Ngoc2023-953; Rawat2024-528. Especially, in simulations, the diffuse gas with
a number density ($n$) less than 100 cm-3 in the envelopes of molecular clouds
may be upheld against gravitational collapse due to magnetic support
Ibanez2022-925; Ganguly2023-525. Assuming the spherical geometry of the
identified clouds, the average number density can be calculated by

$\overline{n}=M/(\frac{4}{3}\pi R_{\rm eff}^{3})/(\mu m_{\rm H}),$ (3)

where $\mu$ = 2.37 is the mean molecular weight per ‘free particle’ (H2 and
He, the number of metal particles is negligible), mH is atomic hydrogen mass.
As shown in Fig.9, almost all of structures have the average number density ¡
100 cm-3. However, the real cloud shape should be more sheetlike instead of
spherical Shetty2006-647; Inutsuka2015-580; Arzoumanian2018-70; Kohno2021-73;
Arzoumanian2022-660; Rezaei2022-930; Zhou2023-519; Zhou2023-676;
Clarke2023-519; Ganguly2023-525. Therefore, the average number density of the
clouds should be significantly underestimated. Even if the average number
density estimates differ by an order of magnitude, what revealed by Fig.9
remains promising. The supportive role of the magnetic field may be a
significant factor contributing to the deviation from free fall motions, and
further detailed investigation is warranted in future research.

##  5 Summary

We investigated the kinematics and dynamics of gas structures on galaxy-cloud
scales in two spiral galaxies NGC5236 (M83) and NGC4321 (M100) using the CO
(2$-$1) line. The main conclusions are as follows:

1\. We directly identified hierarchical (sub-)structures according to the 2D
integrated intensity (Moment 0) map of CO (2$-$1) emission. Subsequently, we
extracted the average spectrum for each structure, delving into its velocity
components and gas kinematics. Considering that the large-scale velocity
gradients due to the galaxy rotation will contribute to the non-thermal
velocity dispersion, before extracting the average spectra of the identified
structures, we subtracted the large-scale velocity gradients in PPV data cube
based on the constructed gas dynamical model.

2\. In examining the scaling relation among velocity dispersion ($\sigma$),
effective radius ($R$), and column density ($N$) across all structures, it
becomes evident that $\sigma-N*R$ consistently exhibits a stronger correlation
compared to $\sigma-N$ and $\sigma-R$. The observed correlations between
velocity dispersion and column density suggest a potential link to
gravitational collapse, corroborated by the measured velocity gradients.
However, it is noteworthy that the slopes of the $\sigma-N*R$ relations appear
to be significantly shallower than the anticipated value of 0.5, implying a
deceleration from the characteristic behavior of a pure free-fall
gravitational collapse.

3\. We employed the FILFINDER algorithm to identify and characterize filaments
within the galaxy using integrated intensity maps. Observable velocity and
density fluctuations along these filaments enabled us to fit local velocity
gradients around intensity peaks, a process performed after removing the
global large-scale velocity gradients attributed to the galaxy’s rotation.

4\. Statistical analysis of the fitted velocity gradients on galaxy-cloud
scales shows the same gas kinematic modes presented on cloud-clump and clump-
core scales. The variations in velocity gradients at both small and large
scales (with a boundary around 500 pc) align with expectations from
gravitational free-fall, with central masses of $\sim$ 105–106.5 M⊙ and $\sim$
106–107.5 M⊙. This implies that sustaining velocity gradients on larger scales
requires correspondingly larger masses, and larger masses imply larger scales,
suggesting that the larger-scale inflow is driven by the larger-scale
structure which may arise from the gravitational clustering of smaller-scale
structures, in harmony with the presence of hierarchical or multi-scale hub-
filament structures within the galaxy and the gas inflows from large to small
scales.

5\. In free-fall model, the velocity gradient and scale satisfy $\nabla
v\propto R^{-1.5}$. However, in the observation, $\nabla v\propto R^{-0.8}$,
also indicating the slowing down of a pure free-fall gravitational collapse.

J. W. Zhou thanks V. Kalinova for the comments. It is a pleasure to thank the
PHANGS team, the data cubes and other data products shared by the team make
this work can be carried out easily. This paper makes use of the following
ALMA data:

ADS/JAO.ALMA#2013.1.01161.S,

ADS/JAO.ALMA#2015.1.00121.S,

ADS/JAO.ALMA#2015.1.00956.S,

ADS/JAO.ALMA#2016.1.00386.S,

ADS/JAO.ALMA#2017.1.00886.L.

ALMA is a partnership of ESO (representing its member states), NSF (USA) and
NINS (Japan), together with NRC (Canada), NSTC and ASIAA (Taiwan), and KASI
(Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA
Observatory is operated by ESO, AUI/NRAO and NAOJ.

## Data Availability

The PHANGS-ALMA CO (2$-$1) data cubes and other data products (such as moment
maps) are available from the PHANGS team website
222https://sites.google.com/view/phangs/home.