galaxies: spiral — galaxies: dwarf — Local Group — galaxies: clusters: general
— dark matter

#  A universal scaling relation incorporating the cusp-to-core transition of
dark matter haloes

Yuka Kaneda11affiliation: Graduate School of Science and Technology,
University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki, 305-8577, Japan
Example: Present Address is xxxxxxxxxx<EMAIL_ADDRESS>Masao
Mori22affiliation: Center for Computational Sciences , University of Tsukuba,
1-1-1 Tennodai, Tsukuba, Ibaraki, 305-8577, Japan Koki Otaki 33affiliation:
Amanogawa Galaxy Astronomy Research Center, Graduate School of Science and
Engineering, Kagoshima University, 1-21-35 Korimoto, Kagoshima, Kagoshima
890-0065, Japan

###### Abstract

The dark matter haloes associated with galaxies have hitherto established
strong correlations within a range of observed parameters, known as scaling
relations of dark matter haloes. The origin of these scaling relations still
contains significant ambiguities and requires comprehensive exploration for
complete understanding. Utilising the correlation between the concentration
and mass of dark matter haloes inferred from cosmological $N$–body simulations
based on the cold dark matter paradigm ($c$–$M$ relation), we derive
theoretical scaling relations among other physical quantities such as the
surface mass density, the maximum circular velocity, and the scale radius of
the dark matter halo. By comparing theoretical and observed scaling relations
at various mass scales, it is found that the scaling relations observed in
dwarf galaxies and galaxies originate in the $c$–$M$ relation of the dark
matter halo. We predict that this theoretical scaling relation is also
established in galaxy clusters. Moreover, we propose a novel theoretical
scaling relation that incorporates the effects of the cusp–to–core transition,
which is supposed to occur in cold dark matter haloes. Our discussion
concludes with the exploration of potential observational verification of the
cusp–to–core transition process in dark matter haloes.

##  1 Introduction

Ninety years have passed since Zwicky (1933) postulated the existence of a
mysterious ’missing mass’ in galaxy clusters, later referred to as dark
matter. Since then, the evidence for dark matter has emerged from a variety of
astrophysical observations, including the rotation curves of galaxies and the
large-scale structure of the universe. These observations suggest that there
is more mass in the universe than can be accounted for by visible matter. The
nature of dark matter remains one of the most profound mysteries in
astrophysics and particle physics.

From an observational point of view, the scaling relations among dark matter
haloes play a vital role in exploring the nature of dark matter. The scaling
relations of dark matter haloes refer to the empirical relationships or
correlations between various properties of the cosmic structures. These
properties often include mass, size, concentration, and other characteristics
of the dark matter haloes. So far, they have been intensively explored,
coupled with the fact that they provide strong constraints on galaxy
formation. In particular, Burkert (1995) found a correlation between
characteristic radii and circular velocities at the characteristic radii of
dark matter haloes associated with galaxies based on the observed data of
seven dwarf galaxies. Following this, some studies have pointed out that the
central surface density of a dark matter halo, represented by the production
of the central density and the core radius, is constant over a wide range with
respect to the absolute $B$-band magnitude of nearby galaxies (Kormendy &
Freeman, 2004; Spano et al., 2008; Donato et al., 2009; Salucci et al., 2012;
Kormendy & Freeman, 2016; Salucci, 2019). On the other hand, recent studies
(Hayashi et al., 2017; Li et al., 2019; Del Popolo, 2023; Gopika et al., 2023)
have claimed that the central surface density weakly depends on the scale of
galaxies. In any case, these correlations should reflect the dynamical
processes which drive the formation and evolution of galaxies. However, the
origin of these empirical relations remains an open question.

The $\Lambda$ cold dark matter (CDM) theory has demonstrated exceptional
efficacy in elucidating and accounting for many of observational phenomena on
scales larger than galaxies. However, this has been considered to encounter
significant challenges when studying the galactic and sub-galactic scales. The
cusp–core problem (Moore, 1994) is one of the notable issues. Cosmological
$N$–body simulations based on the CDM theory yield a universal mass density
distribution that diverges at the centre of dark matter haloes (e.g. Navarro
et al. (1995); Navarro et al. (1996b); Fukushige & Makino (1997); Navarro et
al. (1997); Moore et al. (1999); Navarro et al. (2010); Ishiyama et al.
(2013)). The Navarro–Frenk–White (NFW) profile (Navarro et al., 1997) fits
well the density profiles of dark matter haloes formed in CDM simulations. In
addition, many past studies claim that the density distributions of haloes
formed in cosmological CDM simulations have a tight correlation between their
concentrations, their virial masses, and the redshift $z$, called
concentration–mass relation ($c$–$M$ relation, Prada et al. (2012); Klypin et
al. (2016)). This relation is an essential tool for characterising the mass
density distribution of dark matter haloes formed in CDM simulations.

On the other hand, observations of low–mass galaxies indicate that a
substantial fraction of their dark matter haloes have constant density at
small radii, a distribution called core profile (e.g. Moore (1994); Burkert
(1995); de Blok et al. (2001b); de Blok et al. (2001a); Oh et al. (2008); Oh
et al. (2011); Oh et al. (2015)). The Burkert profile (Burkert, 1995) fits
well the density structure of the observed dark matter haloes with cores. This
long–standing discrepancy between simulated and observed central density
distributions of the dark matter haloes is commonly referred to as the
’cusp–core problem’, and has been the subject of intensive analysis.

Several explanations of this discrepancy have been proposed. One of them is to
consider dark matter particles which intrinsically produce a core, such as
self–interactive dark matter (Spergel & Steinhardt, 2000) and fuzzy dark
matter (Weinberg, 1978; Hui et al., 2017; Burkert, 2020). The other is a
baryonic solution within the CDM framework. The mechanisms that transform a
central cusp into a core have been actively investigated, taking into account
the effect of baryon feedback, which is not included in the cosmological
$N$–body dark matter only simulations. One mechanism is the dynamical
friction, caused by infalling massive gas clumps, heats up dark matter, and
the distribution of dark matter expands, leading to diffuse distribution at
the centre (e.g. El-Zant et al. (2001); Goerdt et al. (2010); Inoue & Saitoh
(2011)). Another and most promising one is the stellar feedback, such as
supernova–driven outflow that alters the gravitational potential of a system
(Navarro et al. (1996a); Gnedin & Zhao (2002); Read & Gilmore (2005); Ogiya &
Mori (2011); Ragone-Figueroa et al. (2012); Ogiya & Mori (2014); for review
see, e.g. Bullock & Boylan-Kolchin (2017); Salucci (2019); Del Popolo & Le
Delliou (2022)). In the following, we refer to these baryonic effects that
alter a cusp into a core as the ”cusp–to–core transition” mechanism.

In the stellar feedback model, Ogiya & Mori (2011) suggested that only one
cycle of stellar feedback is not enough to change the density distribution of
a whole system irreversibly. Ogiya & Mori (2014) showed that recurrent change
of potential arising from the oscillatory feedback process is the primary
driving mechanism of the cusp–to–core transition which generates a dynamically
stable core. Furthermore, they found a simple relationship between the
oscillation period of the stellar feedback and the core radius of the dark
matter haloes assuming that the energy transport process is supported by
Landau resonance.

While Ogiya & Mori (2014) employed an idealised set-up to investigate the
effect of potential oscillation, the more realistic effects of stellar
feedback are investigated utilising zoom–in cosmological hydrodynamic
simulations (e.g. Madau et al. (2014); Fitts et al. (2017); Benítez-Llambay et
al. (2019); Lazar et al. (2020)). One of the important findings is that the
core formation depends on the amount of formed stars in one system. When a
galaxy forms a sufficient number of stars, the energy produced by supernovae
feedback is enough to redistribute dark matter and create a large core. On the
other hand, in the system with the ratio of stellar mass to virial mass less
than $10^{-4}$, the effect of baryon feedback is negligible and a cusp remains
(Di Cintio et al. (2014); Tollet et al. (2016); Lazar et al. (2020)). However,
the approximations used to implement the star–formation and feedback in
simulations have a significant impact on core formation. For example,
individual supernova remnants remain unresolved because their size ($\sim
10\,\mathrm{pc}$) is comparable to the resolution of simulations. The
densities of molecular clouds ($\sim 10^{3-4}\,\mathrm{cm^{-3}}$) are
unresolved and the simulations adopt artificial density thresholds (e.g. $\sim
100\,\mathrm{cm^{-3}}$) for star formation. Therefore, more sophisticated
methods for including star formation and stellar feedback in cosmological
simulations as sub-grid models are expected.

This study aims to investigate the origins of scaling relations in dark matter
haloes and assess the impact of the cusp–to–core transition on these
relations. We seamlessly integrate the scaling relation from dwarf galaxies to
clusters of galaxies introducing a simple and reasonable cusp–to–core
transition model that remains independent of the type of baryonic feedback.
Subsequently, the range of system mass in which the cusp–to–core transition
occurs is discussed. Furthermore, we evaluate how this process and the
critical mass for the transition can be observed using the latest and upcoming
facilities and instruments.

This paper is organised as follows: In section 2, the profiles and the
parameters of dark matter haloes used in this paper are summarised. In section
3, we relate the origin of the scaling relation for galaxy-scale dark matter
haloes to the $c$–$M$ relation. Then, we further investigate the scalings for
MW satellites in detail in section 4. In section 5, targets are extended to
groups and clusters of galaxies. We investigate whether the scaling relations
originating from the $c$–$M$ relation are reproduced in a wide mass range from
dwarf galaxies to clusters of galaxies. Until section 5, we use a cuspy
profile, which is the density distribution of haloes typically arising in
cosmological $N$–body simulations. However, cored profiles more accurately
replicate the observed density distributions. Therefore, we model a process
that transforms a cusp profile into a core profile in order to derive the
parameters of a cored profile, such as core radius and central density, after
the cusp–to-core transition, assuming the parameters of a cusp profile
according to the $c$–$M$ relation. In section 6, we predict a scaling relation
for cored haloes using this model. In section 7, we discuss the capability of
the new scaling relation and application of the cusp–to–core transition model.
To conclude, our work is summarised in section 8.

##  2 Definition of the dark matter halo properties

The density distribution of a CDM halo is well reproduced by the NFW profile:

$\rho(r)=\frac{\rho_{\mathrm{N}}}{\frac{r}{r_{\mathrm{N}}}\left(\frac{r}{r_{\mathrm{N}}}+1\right)^{2}},$
(1)

where $r$ is the distance from the centre and $\rho_{\mathrm{N}}$ and
$r_{\mathrm{N}}$ represent the scale density and scale radius of a halo,
respectively. The mass enclosed within the radius $r$ can be derived by
integrating equation (1) as

$M(<r)=4\pi\rho_{\mathrm{N}}r_{\mathrm{N}}^{3}f_{\mathrm{N}}\left(\frac{r}{r_{\mathrm{N}}}\right),$
(2)

where

$f_{\mathrm{N}}(x)=\ln\left(1+x\right)-\frac{x}{1+x}.$ (3)

Many previous studies (see section 1) have claimed that these haloes with
central cusps formed in cosmological CDM simulations present a tight
correlation between their concentrations of mass $c_{200}$, their masses
$M_{200}$, and the redshift $z$, called the $c$–$M$ relation. $c_{200}$ is
defined by

$c_{200}\equiv\frac{r_{200}}{r_{\mathrm{N}}}.$ (4)

The relational expression between the virial mass $M_{200}$ and the virial
radius $r_{200}$ for $z=0$ can be written by

$M_{200}\equiv\frac{4}{3}\pi 200\rho_{\mathrm{crit,0}}r_{200}^{3}.$ (5)

$\rho_{\mathrm{crit},0}=3H_{0}^{2}/(8\pi G)$ is the critical density of the
universe at $z=0$. $H_{0}$ is the Hubble constant and can be expressed as
$H_{0}=100\,h\,\rm{km\,s^{-1}\,Mpc}$ using the dimensionless Hubble constant
$h=0.674$ (Planck Collaboration et al., 2020).

The circular velocity is defined as

$V(r)=\left(\frac{GM(<r)}{r}\right)^{1/2}$ (6)

Then the maximum circular velocity $V_{\mathrm{max}}$ is given by
$V_{\mathrm{max}}=V(r_{\mathrm{max}})$ with the radius at which the circular
velocity becomes maximum, $r_{\mathrm{max}}$. $r_{\mathrm{max}}$ can be
derived by solving $dV/dr=0$. Using the NFW profile as a mass distribution, we
obtain

$r_{\mathrm{max}}=x_{\mathrm{N}}\,r_{\mathrm{N}}$ (7)

where $x_{\mathrm{N}}=2.16$.

There are two methods to estimate halo concentration from a resulting
distribution of particles in $N$–body simulations. Firstly, one can calculate
the spherically averaged density distribution and then fit it with an
appropriate function. In this case, the results may depend upon the selected
binning algorithm of particle position and the fit quality. Secondary, instead
of making fits to individual halo profiles, one can find the concentration by
solving the following equations numerically, having the
$V_{\mathrm{max}}/V_{200}$ ratio for each halo (Klypin et al., 2011; Prada et
al., 2012):

$\frac{V_{\mathrm{max}}}{V_{200}}=\left(\frac{c_{200}f_{\mathrm{N}}(x_{\mathrm{N}})}{x_{\mathrm{N}}f_{\mathrm{N}}(c_{200})}\right)^{1/2}.$
(8)

Here, we regard the virial velocity $V_{200}$ to be
$V_{200}=(GM_{200}/r_{200})^{1/2}$. In this paper, we use the concentrations
found by the later method.

The Burkert profile (Burkert, 1995) fits well the density structures of dark
matter haloes with cores,

$\rho(r)=\frac{\rho_{\mathrm{B}}}{\left(\frac{r}{r_{\mathrm{B}}}+1\right)\left(\frac{r^{2}}{{r_{\mathrm{B}}}^{2}}+1\right)},$
(9)

where $\rho_{\mathrm{B}}$ and $r_{\mathrm{B}}$ are a halo’s central density
and core radius, respectively. The mass enclosed within the radius $r$ can be
derived by integrating equation (9) as

$M(<r)=4\pi\rho_{\mathrm{B}}r_{\mathrm{B}}^{3}f_{\mathrm{B}}\left(\frac{r}{r_{\mathrm{B}}}\right),$
(10)

where

$f_{\mathrm{B}}(x)=\frac{1}{4}\ln\left(x^{2}+1\right)+\frac{1}{2}\ln(x+1)-\frac{1}{2}\arctan
x.$ (11)

$r_{\mathrm{max}}$ for the Burkert profile can also be derived by solving
$dV/dr=0$ which gives

$r_{\mathrm{max}}=x_{\mathrm{B}}\,r_{\mathrm{B}},$ (12)

where $x_{\mathrm{B}}=3.24$.

The pseudo–isothermal sphere (p-ISO) can also represent a cored profile,

$\rho(r)=\frac{\rho_{\mathrm{p-ISO}}}{1+\frac{r^{2}}{r_{\mathrm{p-ISO}}^{2}}}.$
(13)

where $r_{\mathrm{p-ISO}}$ is the core radius of p-ISO and
$\rho_{\mathrm{p-IS0}}$ is the central density. The mass enclosed within the
radius $r$ can be written as

$M(<r)=4\pi\rho_{\mathrm{p-ISO}}r_{\mathrm{p-ISO}}^{3}f_{\mathrm{p-ISO}}\left(\frac{r}{r_{\mathrm{p-ISO}}}\right),$
(14)

where

$f_{\mathrm{p-ISO}}(x)=x-\arctan x.$ (15)

##  3 Origin of the observed scaling relations

Original scaling relations and the relations converted to
$\bar{\Sigma}(<r_{\mathrm{max}})$ as a function of $V_{\mathrm{max}}$.  Name
Burkert (1995) Spano et al. (2008) original $\rho_{\mathrm{B}}=4.5\times
10^{-2}(r_{\mathrm{B}}/\mathrm{kpc})^{-2/3}\,M_{\odot}\,\mathrm{pc}^{-3}$
$\rho_{\mathrm{p-ISO}}r_{\mathrm{p-ISO}}\sim 150\,M_{\odot}\mathrm{pc}^{-2}$
converted to $\bar{\Sigma}(<r_{\mathrm{max}},V_{\mathrm{max}})$ [$*$]
$\bar{\Sigma}(<r_{\mathrm{max}})=2.5\\{V_{\mathrm{max}}/(\mathrm{km\,s}^{-1})\\}^{1/2}M_{\odot}\,\mathrm{pc}^{-2}$
$\bar{\Sigma}(<r_{\mathrm{max}})\sim 80\,M_{\odot}\,\mathrm{pc}^{-2}$
**footnotemark: $*$ Donato et al. (2009) Donato et al. (2009) dwarf Kormendy &
Freeman (2016)
$\rho_{\mathrm{B}}r_{\mathrm{B}}=140^{+80}_{-30}\,M_{\odot}\mathrm{pc}^{-2}$
$\rho_{\mathrm{B}}r_{\mathrm{B}}\sim 100\,M_{\odot}\mathrm{pc}^{-2}$
$\rho_{\mathrm{NI}}r_{\mathrm{NI}}=70\pm 4\,M_{\odot}\,\mathrm{pc}^{-2}$
$\bar{\Sigma}(<r_{\mathrm{max}})=37^{+21}_{-8}\,M_{\odot}\,\mathrm{pc}^{-2}$
$\bar{\Sigma}(<r_{\mathrm{max}})\sim 27\,M_{\odot}\,\mathrm{pc}^{-2}$
$\bar{\Sigma}(<r_{\mathrm{max}})=38\pm 2\,M_{\odot}\,\mathrm{pc}^{-2}$
{\dagger}{\dagger}footnotemark: ${\dagger}$ **footnotemark: $*$ The parameters
of p-ISO are converted to the ones of the Burkert profile utilising the method
described in appendix Matching the p-ISO to the Burkert profile since p-ISO
profile do not have maximum circular velocity.
{\dagger}{\dagger}footnotemark: ${\dagger}$ The parameters of the non-singular
isothermal sphere are considered to be the same as the parameters of p-ISO.

Figure 1: Left panel: a comparison between the scaling relations proposed by
Burkert (1995) (magenta line), Spano et al. (2008) (blue line), Donato et al.
(2009) (sky–blue line for galaxies, green line for dwarf galaxies), Salucci et
al. (2012) (yellow symbol), Kormendy & Freeman (2016) (orange line) and the
theoretical $\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ relation by
Moliné et al. (2023) (grey line). The coloured line segments start and end
with the domains of definitions of each scaling relation. The propagated
errors are shown in the shaded region in sky-blue (Donato et al., 2009) and in
orange, which is narrower than the thickness of the solid line (Kormendy &
Freeman, 2016). As a theoretical
$\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ relation, we adopt the
fitting curve of the $c_{V}$–$V_{\mathrm{max}}$ relation of subhaloes from the
Uchuu suite (Ishiyama et al., 2021) for $z=0$ proposed by Moliné et al. (2023)
(equation 18 with $C=c_{\mathrm{V}}(V_{\max})$ and
$S=V_{\max}\,\mathrm{km}\,\mathrm{s}^{-1}$, and the coefficients shown in
equation 19) converted to the
$\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ plane using equation
(23). Right panel: the median of $\bar{\Sigma}(<r_{\mathrm{max}})$ against
$V_{\mathrm{max}}$ derived for subhaloes of MW–like hosts in the high-
resolution cosmological $N$–body simulation, Phi-4096 (Ishiyama et al. (2021),
the thick solid line). The fitting curve of the $c_{V}$–$V_{\mathrm{max}}$
relation of MW subhaloes from Phi-4096 for $z=0$ (equation 18 with
$C=c_{\mathrm{V}}(V_{\max})$ and $S=V_{\max}\,\mathrm{km}\,\mathrm{s}^{-1}$,
and the coefficients shown in equation 24) converted to the
$\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ plane using equation (23)
is shown in black dotted line. The shaded region shows a 25-75th percentile
around it. The red symbols show observed MW and M31 satellites.

In this section, we attribute the origin of the empirical scaling relation for
dark matter haloes introduced in section 1 to the theoretical universality,
the $c$–$M$ relation, derived from cosmological dark matter only simulations.
We transform the empirical scaling relation into equations on the surface
density and $V_{\rm max}$ plane. Similarly, the $c$–$M$ relations by Moliné et
al. (2023), derived from the recent cosmological $N$–body simulation (Ishiyama
et al., 2021), is projected into relationships between surface density and
$V_{\rm max}$. The mutual connection between the empirical and theoretical
scaling relation, which govern the parameters of dark matter haloes, is
subsequently examined.

Burkert (1995) suggested the existence of a simple correlation between the
core radius of the Burkert profile $r_{\mathrm{B}}$ and the central density
$\rho_{\mathrm{B}}$. This relation is based on seven nearby dwarf galaxies
which are clearly dark matter dominated at their core radius. The parameters
of the dark matter haloes are determined by fitting the theoretically
predicted rotation curve to the observational $\mathrm{H_{I}}$ rotation
curves. They derived the following correlation: $\rho_{\mathrm{B}}=4.5\times
10^{-2}(r_{\mathrm{B}}\;\mathrm{kpc}^{-1})^{-2/3}M_{\odot}\,\mathrm{pc}^{-3}$.

Kormendy & Freeman (2004) firstly pointed out that the central surface
density, the product of a central density and a core radius is constant and
independent of galaxy luminosity. This relation was updated in Kormendy &
Freeman (2016). Kormendy & Freeman (2016) claimed that central surface
densities of dark matter haloes of Sc-Im galaxies and dwarf galaxies are
constant independent of their absolute $B$-band magnitude,
$\rho_{\mathrm{NI}}r_{\mathrm{NI}}=70\pm 4\,M_{\odot}\,\mathrm{pc}^{-2}$,
where $\rho_{\mathrm{NI}}$ is the central density and $r_{\mathrm{NI}}$ is the
scale radius of the nonsingular isothermal sphere. However, the central
surface densities of dwarf galaxies were artificially increased to match their
scaling relation with that of disk galaxies. Therefore we exclude the mass
range of dwarf galaxies from the domain of the Kormendy-Freeman relation shown
in the orange line in figure 1a.

Spano et al. (2008) studied the mass distribution of 36 spiral galaxies based
on rotation curve decomposition. p-ISO is adopted as a dark matter density
distribution. They suggested nearly constant central surface density around
$\sim 150\,M_{\odot}\mathrm{pc}^{-2}$ independent of galaxy luminosity.

Donato et al. (2009) extended this finding. Their sample comprised the
rotation curves of $\sim 1000$ spiral galaxies, the mass models of individual
dwarf irregular and spiral galaxies of late and early types with high-quality
rotation curves, and the galaxy–galaxy weak gravitational lensing signals from
a sample of spiral and elliptical galaxies. They determined that
$\rho_{\mathrm{B}}r_{\mathrm{B}}=140^{+80}_{-30}\,M_{\odot}\mathrm{pc}^{-2}$
for Burkert profile. They also showed that the observed kinematics of 6 Local
Group dwarf spheroidal galaxies indicate $\rho_{\mathrm{B}}r_{\mathrm{B}}\sim
100\,M_{\odot}\mathrm{pc}^{-2}$.

Furthermore, Salucci et al. (2012) extended this relation to the MW dSphs.
They performed Jeans analysis for eight MW dSphs and derived the parameters of
the Burkert profiles. They showed that the MW dSphs lie on the extrapolation
of the scaling law seen in spiral galaxies. Note that they did not fit the
eight dwarf galaxies with a function of constant central surface density; just
confirmed that the dSphs were consistent with the extrapolation of the spiral
galaxy scaling proposed by Donato et al. (2009).

To compare these observational scaling relations, assuming a core profile,
with halo properties from CDM simulations, we analyse the scaling relations
through a parameter that can be defined for both core and cusp profiles. Here,
we define the characteristic mean surface density of the dark matter halo
within the radius of the maximum circular velocity as follows:

$\bar{\Sigma}(<r_{\mathrm{max}})=\frac{M(<r_{\mathrm{max}})}{\pi
r_{\mathrm{max}}^{2}},$ (16)

where

$M(<r_{\mathrm{max}})=\int_{0}^{r_{\mathrm{max}}}4\pi\rho_{\mathrm{dm}}\left(r^{\prime}\right)r^{\prime
2}dr^{\prime},$ (17)

$\rho_{\mathrm{dm}}$ corresponds to the density profile of a dark matter halo
(see Hayashi & Chiba (2015); Hayashi et al. (2017)). Note that this quantity
differs from the usual surface density, i.e., the mass contained in the sphere
divided by the area, not the mass contained in the cylinder. If we assume the
NFW profile, the Burkert profile or p-ISO as $\rho_{\mathrm{dm}}$,
$\bar{\Sigma}(<r_{\mathrm{max}})$ is proportional to the product of scale
density and the scale radius or the central density and the core radius. This
means that $\bar{\Sigma}(<r_{\mathrm{max}})$ is several times the central
surface density that is reported to be constant for galaxies and dwarf
galaxies in the literature. The original observational scaling relations and
the ones converted to the $V_{\mathrm{max}}-\bar{\Sigma}(<r_{\mathrm{max}})$
plane are summarised in table 3. To derive
$\bar{\Sigma}(<r_{\mathrm{max}})-V_{\mathrm{max}}$ relation, the parameters of
p-ISO are converted to the ones of the Burkert profile using the method
described in appendix Matching the p-ISO to the Burkert profile since p-ISO do
not have maximum circular velocity. The parameters of the non-singular
isothermal sphere are considered to be the same as the parameters of p-ISO.
Figure 1a shows the scaling relations given in table 3 on the
$V_{\mathrm{max}}-\bar{\Sigma}(<r_{\mathrm{max}})$ plane. Solid lines with
different colours represent each scaling relation. The propagated errors are
shown in the shaded region in sky-blue (Donato et al., 2009) and in orange,
which is narrower than the thickness of the solid line (Kormendy & Freeman,
2016). The properties of dwarf galaxies estimated by Salucci et al. (2012) are
shown in the yellow symbols.

We explore the connection between these scaling relations derived from
observations and the theoretical properties of dark matter haloes. Theoretical
analyses based on cosmological simulations have widely recognized the
existence of universal properties of the cold dark matter haloes. They
manifest as a strong correlation between the total mass and central
concentration of the dark matter halo, often referred to as the $c$–$M$
relation. In the following, we will investigate the trajectory of this $c$–$M$
relation on the $V_{\mathrm{max}}-\bar{\Sigma}(<r_{\mathrm{max}})$ plane.

We adopt the fitting curve of the $c$–$M$ relation of subhaloes from the Uchuu
suite (Ishiyama et al., 2021) for $z=0$ proposed by Moliné et al. (2023):

$C(S)=c_{0}\left[1+\sum_{i=1}^{3}a_{i}\left[\log_{10}\left(S\right)\right]^{i}\right],$
(18)

where $C=c_{\mathrm{V}}(V_{\max})$, $S=V_{\max}/\mathrm{km}\,\mathrm{s}^{-1}$
and

$c_{0}=1.75\times 10^{5}\quad\mathrm{and}\quad a_{i}=\\{-0.9,0.27,-0.02\\}.$
(19)

These coefficients are also listed in table 3. Here, the definition of
$c_{\mathrm{V}}$ is

$c_{\mathrm{V}}=\frac{\bar{\rho}(<r_{\mathrm{max}})}{\rho_{\mathrm{crit,0}}}=2\left(\frac{V_{\mathrm{max}}}{H_{0}r_{\mathrm{max}}}\right)^{2},$
(20)

where $\bar{\rho}(<r_{\mathrm{max}})$ is the mean physical density attained
within the radius corresponding to the maximum circular velocity,
$r_{\mathrm{max}}$, expressed in units of the critical density of the Universe
(Diemand et al., 2007; Moliné et al., 2017). The relation between
$c_{\mathrm{V}}$ and $c_{200}$ is given by Diemand et al. (2007),

$c_{\mathrm{V}}=200\left(\frac{c_{200}}{x_{\mathrm{N}}}\right)^{3}\frac{f_{\mathrm{N}}(x_{\mathrm{N}})}{f_{\mathrm{N}}(c_{200})}.$
(21)

The maximum circular velocity is connected with the virial mass by equation
(8). In this way, a $c$–$M$ relation can be written as a
$c_{\rm{V}}$–$V_{\mathrm{max}}$ relation. Equation (18) with (19) is valid for
$7\mathrm{~{}km}\mathrm{~{}s}^{-1}\lesssim V_{\max}\lesssim
1500\mathrm{~{}km}\mathrm{~{}s}^{-1}$ due to the limitation of halo masses in
the simulation. The $c_{\rm{V}}$–$V_{\mathrm{max}}$ relation (equation 18 and
19) expressed on the $c_{200}-M_{200}$ plane is also described by equation
(18) with $C=c_{200}(M_{200})$, $S=M_{200}/M_{\odot}$, and

$c_{0}=48.29\quad\mathrm{and}\quad a_{i}=\\{-0.01863,-0.008262,0.0003651\\}$
(22)

The error of this parameterisation is less than 3 per cent. The
$c_{\rm{V}}$–$V_{\mathrm{max}}$ relation relation can also be expressed in
terms of $\bar{\Sigma}(<r_{\mathrm{max}})$ written as a function of
$V_{\mathrm{max}}$. Solving equation (20) for $r_{\mathrm{max}}$ and
substitute it into equation (16) yields

$\displaystyle\frac{\bar{\Sigma}(<r_{\mathrm{max}})}{M_{\odot}\,\mathrm{pc}^{-2}}$
$\displaystyle=1.10\times
10^{-12}\;\frac{h}{\pi}\left(\frac{G}{\mathrm{m^{3}\,kg^{-1}\,s^{-2}}}\right)^{-1}$
(23)
$\displaystyle\times\frac{V_{\mathrm{max}}}{\mathrm{km\,s^{-1}}}\left\\{c_{\mathrm{V}}\left(\frac{V_{\mathrm{max}}}{\mathrm{km\,s^{-1}}}\right)\right\\}^{1/2}.$

The grey dashed curve in figure 1a is the theoretical
$\bar{\Sigma}(<r_{\mathrm{max}})$ curve drawn on the basis of the
$c_{\rm{V}}$–$V_{\mathrm{max}}$ relation of the CDM haloes in equation (23).
This function is slowly monotonically increasing on this plane. The scaling
relations and the theoretical curve are consistent for the domain of the
definition of each scaling relation. The relations proposed by
donato_constant_2009 (donato_constant_2009, sky–blue line) and
kormendy_scaling_2016 (kormendy_scaling_2016, orange line) correspond to the
intermediate mass range around $V_{\mathrm{max}}\sim 100\,\mathrm{km\,s^{-1}}$
and they coincide with the theoretical
$\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ relation. The relation
proposed by spano_ghasp_2008 (spano_ghasp_2008, blue line) shows the largest
velocities, from $100$ to $1000\,\mathrm{km\,s^{-1}}$, and the highest
characteristic mean surface density, around
$100\,M_{\odot}\,\mathrm{pc^{-2}}$. The relation proposed by Spano et al.
(2008) is also consistent with the theoretical
$\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ relation. The empirical
scaling relations seem to be broken fragments of the theoretical
$\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ relation which is derived
from $c_{\rm{V}}$–$V_{\mathrm{max}}$ relation. In other words, what we have
thought of as empirical scaling relations between parameters of the dark
matter haloes is, actually, the $c$–$M$ relation, which roots in the nature of
the CDM haloes, seen from a different perspective. Observations over a limited
and narrow mass range make the observed scaling relations appear constant,
just as the theoretical $\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$
relation predicts. However, the theoretical
$\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ relation decreases as the
mass of the dark matter halo decreases. Therefore, the arguments that the
surface density depends on the scale of a galaxy (Hayashi et al., 2017; Li et
al., 2019; Del Popolo, 2023; Gopika et al., 2023) can be understood as a
natural consequence of the $c$–$M$ relation.

The scaling relation of dwarf galaxies has been poorly discussed in previous
studies. donato_constant_2009 (donato_constant_2009, green line) fitted only
six dwarf galaxies with a constant. salucci_dwarf_2012 (salucci_dwarf_2012,
yellow symbol) increases the number of dwarf galaxies to eight, however, they
have too large uncertainties to limit the scaling relation of dwarf galaxies.
In the next section, we discuss the scaling relation of dwarf galaxies
($V_{\mathrm{max}}\lesssim 100\,\mathrm{km\,s^{-1}}$) based on recent
estimations with improved observation.

Coefficients of equation (18).**footnotemark: $*$ $C$ $S$
environment{\dagger}{\dagger}footnotemark: ${\dagger}$
range{\ddagger}{\ddagger}footnotemark: ${\ddagger}$ profile\S\Sfootnotemark:
$\S$ $c_{0}$ $a_{1}$ $a_{2}$ $a_{3}$ \|\|footnotemark: $\|$$c_{\mathrm{V}}$
$V_{\rm{max}}\,\rm{km\,s^{-1}}$ global $7\sim 1500\,\rm{km\,s^{-1}}$ cusp
$1.75\times 10^{5}$ $-0.9$ $0.27$ $-0.02$ $c_{200}$ $M_{200}\,M_{\odot}$
global $10^{7}\sim 10^{15}\,M_{\odot}$ cusp $48.29$ $-0.01863$ $-0.008262$
$0.0003651$ $c_{\mathrm{V}}$ $V_{\rm{max}}\,\rm{km\,s^{-1}}$ MW subhalo $1\sim
40\,\rm{km\,s^{-1}}$ cusp $4.978\times 10^{4}$ $13.62$ $-17.33$ $5.439$
$c_{200}$ $M_{200}\,M_{\odot}$ MW subhalo $10^{6}\sim 10^{10}\,M_{\odot}$ cusp
$-2.947\times 10^{2}$ $-0.4375$ $0.05414$ $-0.002101$ $c_{\mathrm{V}}$
$V_{\rm{max}}\,\rm{km\,s^{-1}}$ global $7\sim 1500\,\rm{km\,s^{-1}}$ core
$4.925\times 10^{4}$ $-0.8505$ $0.2434$ $-0.02314$ $c_{\rm{B}200}$
$M_{200}\,M_{\odot}$ global $10^{7}\sim 10^{15}\,M_{\odot}$ core $41.21$
$-0.01431$ $-0.008294$ $0.0003589$ $c_{\mathrm{V}}$
$V_{\rm{max}}\,\rm{km\,s^{-1}}$ MW subhalo $1\sim 40\,\rm{km\,s^{-1}}$ core
$2.251\times 10^{4}$ $7.490$ $-10.02$ $3.220$ $c_{\rm{B}200}$
$M_{200}\,M_{\odot}$ MW subhalo $10^{6}\sim 10^{10}\,M_{\odot}$ core
$-2.386\times 10^{2}$ $-0.4412$ $0.05455$ $-0.002114$ **footnotemark: $*$Best-
ﬁt values of the parametrisations for the concentration parameter, $c_{200}$
or $c_{\rm{V}}$, as a function of $V_{\rm{max}}$ or $M_{200}$ (equation 18).
{\dagger}{\dagger}footnotemark: ${\dagger}$The third column, ”environment”,
shows the environments the parametrisations are for. ”MW subhalo” means that
the parameters are for the subhaloes of the MW-like host. ”Global” means that
the parameters are for the subhaloes without limitation on the mass of the
host halo. The ”global” relation is derived from
$c_{\rm{V}}$–$V_{\mathrm{max}}$ relation proposed in Moliné et al. (2023).
{\ddagger}{\ddagger}footnotemark: ${\ddagger}$The fourth column, ”range”,
shows the available range of these fittings. \S\Sfootnotemark: $\S$The fifth
column, ”profile”, shows which profile the parametrisations are for. ”Cusp”
means the parametrisation is the original relation derived from the results of
the cosmological $N$-body simulations and ”core” means the function is derived
by applying our cusp–to–core transition model to the results of the
simulations (see section 6 for detail).
\|\|footnotemark: $\|$The first line shows the coefficients for the original
$c_{\rm{V}}$–$V_{\mathrm{max}}$ relation in Moliné et al. (2023). Significant
figures are taken directly from Moliné et al. (2023). Other approximation
functions proposed in this paper are unified into four significant digits.

##  4 Theoretical $\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$
relation in the Milky Way–like environment

Recent developments in instruments and techniques have led to remarkable
progress in observing nearby galaxies. Moreover, near-future observations of
satellite galaxies in the MW and M31 using state–of–the–art instruments (such
as the Subaru Prime Focus Spectrograph on the Subaru Telescope and the James
Webb Space Telescope) will provide new insights into the nature of dark
matter. This situation motivates us to explore more precisely the theoretical
scaling relation for dwarf–galaxy–scale dark matter haloes in a Milky Way–like
environment, analysing a data set of one of the highest–mass–resolution
cosmological $N$–body simulations available described in the following
subsection. Here, we propose a new theoretical scaling relation for dwarf
galaxies based on the $c$–$M$ relation of the cosmological $N$–body simulation
which is of such high resolution that it allows us to study down to the scale
of dwarf galaxies.

###  4.1 Description of simulation data

The cosmological $N$–body simulation, Phi–4096, which has the highest
resolution among the Uchuu suite by Ishiyama et al. (2021), is utilised to
derive a $c$–$M$ relation for the dwarf–galaxy–mass haloes. Phi–4096 was
conducted to extend the halo mass range down at high redshift and evaluate the
resolution effects of the other simulations of the Uchuu suite. The basic
properties of Phi–4096 are listed in table 4.1. Since its particle mass is
$\sim 10^{3}\,M_{\odot}$, a subhalo with $10^{9}\,M_{\odot}$ consists of $\sim
10^{6}$ particles. Thanks to its relatively large volume and very high mass
resolution, it has a large number of MW–like haloes and subhaloes. Thus
Phi–4096 enables the statistical study and characterization of subhaloes with
a uniform resolution. The cosmological parameters used are the matter density
parameter $\Omega_{0}=0.31$, the baryon density parameter
$\Omega_{\mathrm{b}}=0.048$, the cosmological constant $\lambda_{0}=0.69$,
Hubble constant $h=0.68$, the spectral index of the density perturbation
$n_{\mathrm{s}}=0.96$, and the matter fluctuation amplitude $\sigma_{8}=0.83$.
These are consistent with the observational result of the Plank Satellite
(Planck Collaboration et al., 2020). Bound structures at each snapshot are
identified using the rockstar phase space halo$/$subhalo finder (Behroozi et
al., 2013a). The consistent trees code (Behroozi et al., 2013b) is used to
construct merger trees.

Properties of Phi–4096.**footnotemark: $*$ Name $N$ $\begin{array}[]{c}L\\\
(h^{-1}\mathrm{Mpc})\end{array}$ $\begin{array}[]{c}\varepsilon\\\
(h^{-1}\mathrm{kpc})\end{array}$ $\begin{array}[]{c}m_{\mathrm{p}}\\\
(h^{-1}\mathrm{M}_{\odot})\end{array}$ Phi–4096 $4096^{3}$ $16.0$ $0.06$
$5.13\times 10^{3}$ **footnotemark: $*$ $N$, $L$, $m_{\mathrm{p}}$, and
$\varepsilon$ are the total number of particles, box length, softening
parameter, and particle mass resolution, respectively.

We analyse nearly three hundred thousand subhaloes belonging to 27 MW–sized
host haloes virial mass ranging from $3.40\times 10^{11}$ to $2.04\times
10^{12}\,h^{-1}M_{\odot}$ at redshift $z=0$. The following two criteria are
imposed to select the subhaloes. First, to avoid the effects caused by the
mass resolution, the virial mass needs to be larger than $10^{6}M_{\odot}$ so
that a halo has at least 1,000 dark matter particles. Second, the subhalo has
to reside within twice the virial radius of a host halo at redshift $z=0$.

###  4.2 $\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ relation for MW
satellite-like haloes

Figure 2:  Top panel: median concentration (solid lines), 25-75th percentile
(shaded area), and fit (dashed lines) for a given
$c_{\rm{V}}$–$V_{\mathrm{max}}$ data of subhaloes of MW–like hosts in high-
resolution cosmological $N$–body simulation. Bottom panel: the relative
difference between fit and data.

In this section, the $c$–$M$ relation derived from the subhaloes of MW–like
hosts is compared with recent observations of MW satellites. We fit the median
of the concentrations inferred from the simulation by equation (18). The
fitting results are shown as the sky–blue solid (simulation) and dashed (fit)
lines in the top panel of figure 2. To make binned data, the exponent of the
mass $M_{200}$ is divided into increments of $0.1$. We exclude ranges
containing fewer than $10$ haloes per bin. Consequently, the data are valid
for $1<V_{\mathrm{max}}<40$ km s-1. The resultant fitting function, namely the
$c_{\rm{V}}$–$V_{\mathrm{max}}$ relation for subhaloes of MW–sized hosts, is
expressed as in equation (18) with $C=c_{\mathrm{V}}(V_{\max})$,
$S=V_{\max}/\mathrm{km}\,\mathrm{s}^{-1}$ and

$c_{0}=4.978\times 10^{4}\;\mathrm{and}\;a_{i}=\\{13.62,-17.33,5.439\\}.$ (24)

The sky–blue line in the bottom panel of figure 2 illustrates the relative
error in the fitting function for the cusp profile. It demonstrates that the
error remains below 10 per cent with $V_{\mathrm{max}}$ between
$1.5\,\mathrm{km\,s}^{-1}$ and $15\,\mathrm{km\,s}^{-1}$, indicating that is a
reliable approximation within this range. Furthermore, we fit the median
$c_{200}$ as a function of $M_{200}$ through equation (18) with
$C=c_{200}(M_{200})$, $S=M_{200}/M_{\odot}$ and find the coefficients to be

$c_{0}=-2.947\times
10^{2}\;\mathrm{and}\;a_{i}=\\{-0.4375,0.05414,-0.002101\\}.$ (25)

These coefficients are organised in table 3.

We adopted the observational estimation of the density profile of satellites
of MW and M31 derived by Hayashi et al. (2017). Hayashi et al. (2017)
estimated the dark matter halo profile with a higher precision than previous
studies by fitting the axisymmetric Jeans equation to the line–of–sight
velocities obtained from spectroscopic observations. They consider the
following density profile

$\displaystyle\rho(R,z)$ $\displaystyle=$
$\displaystyle\rho_{0}\left(\frac{m}{b_{\mathrm{halo}}}\right)^{\alpha}\left[1+\left(\frac{m}{b_{\mathrm{halo}}}\right)^{2}\right]^{-(\alpha+3)/2},$
(26) $\displaystyle m^{2}$ $\displaystyle=$ $\displaystyle R^{2}+z^{2}/Q^{2},$
(27)

where $\rho_{0}$ is a scale density, $b_{\mathrm{halo}}$ is a scale length in
the spatial distribution, and $Q$ is an axial ratio for the dark matter halo.

The comparison between the
$\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ relation of MW-satellite-
like haloes and observation of MW and M31 satellites is shown in figure 1b.
The solid grey line describes the median of $\bar{\Sigma}(<r_{\mathrm{max}})$
of MW-satellite-like haloes in the cosmological $N$-body simulation Phi-4096
and the shaded region illustrates the range between the 25th percentile and
the 75th percentile. The black dotted curve is the theoretical
$\bar{\Sigma}(<r_{\mathrm{max}})$ curve drawn on the basis of the
$c_{\rm{V}}$–$V_{\mathrm{max}}$ relation (equation 18 and 24) substituted in
equation (23). Note that the applicable range for the theoretical
$\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ relation, derived from
our $c$–$M$ relation, is $1<V_{\mathrm{max}}<40$ km s-1, which corresponds to
the typical range for subhaloes.

Figure 1b shows that the red symbols indicating the eight MW classical dwarfs
and five M31 dwarfs measured by Hayashi et al. (2017) are in excellent
agreement with this theoretical
$\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ relation. Because of the
limited number and significant errors of the observed dwarf galaxies in the
previous studies, the characteristic mean surface density inside
$r_{\mathrm{max}}$ has been proposed to be constant. However, the theoretical
$\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ relation we derived here
is consistent with recent observations of MW dwarf galaxies, showing that the
average surface density inside $r_{\mathrm{max}}$ depends on the system size.
We expect future observation using state-of-the-art telescopes and instruments
to estimate the density distributions of haloes of satellite galaxies with
$V_{\mathrm{max}}<20\,\mathrm{km\,s^{-1}}$.

##  5 Theoretical scaling relation from dwarf galaxies to galaxy clusters

In Section 3, we compared observations of the
$\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ relations on galaxy scale
with that derived from the $c$–$M$ relation predicted by the CDM model. The
result clearly shows that the observed scaling relations are deeply rooted in
the CDM prediction. In this section, we extend the theoretical scaling
relation not only to $V_{\mathrm{max}}-\bar{\Sigma}(<r_{\mathrm{max}})$ plane
but also to other physical quantities. Then, these relations are compared with
the observed quantities for a wide range of astronomical objects, from dwarf
galaxies to galaxy clusters.

We adopt the fitting function (equations 18 and 19) proposed by Moliné et al.
(2023). While the fitting function of theoretical $c$–$M$ relation we propose
(equation 24) is applicable for subhaloes associated with the MW–sized hosts
($10^{6}\sim 10^{10}M_{\odot}$), the fitting function suggested by Moliné et
al. (2023) is available for broader mass range $10^{7}\sim 10^{14}M_{\odot}$.
To compare the theoretical $c$–$M$ relation with observations ranging from
dwarf galaxies to galaxy clusters, we adopt the fitting function proposed by
Moliné et al. (2023) in the following subsections.

###  5.1 Description of employed observational data

This section describes employed observational data and selection criteria
adopted to ensure enough data quality.

We adopted the parameters which are derived assuming that haloes have an NFW
profile in order to simply compare them with the parameters of haloes inferred
from cosmological $N$–body simulations.

For dwarf galaxies, we employ data from Oh et al. (2015). Oh et al. (2015)
derived rotation curves from “Local Irregulars That Trace Luminosity Extremes,
The $\mathrm{H_{I}}$ Nearby Galaxy Survey” (LITTLE THINGS) and constructed
mass models of 26 galaxies. 17 $c_{200}$ and $V_{200}$ on their table 2 are
used.

For massive galaxies, we adopt the data from rotation curves obtained with
spectroscopic observations by Spano et al. (2008), de Blok et al. (2008), and
Sofue (2016). Spano et al. (2008) derived the mass models of 36 spiral
galaxies using rotation curve analysis. They combine the $\mathrm{H_{\alpha}}$
2D velocity field obtained using Fabry–Perot interferometry as part of the
GHASP survey and $\mathrm{H_{I}}$ data from the literature. $r_{\mathrm{N}}$
(columun 10) and $\rho_{\mathrm{N}}$ (columun 11) of their table 2 are used.
de Blok et al. (2008) presented the mass models of 19 galaxies using rotation
curve analysis from The $\mathrm{H_{I}}$ Nearby Galaxy Survey (THINGS).
$c_{200}$ (columun 8) and $V_{200}$ (columun 9) of their table 6 are used.
Sofue (2016) compiled rotation curves of more than one hundred spiral galaxies
from the literature. $r_{\mathrm{N}}$ (columun 7), $r_{200}$ (columun 9), and
$M_{200}$ (columun 10) of their table 1 are used.

Groups and clusters of galaxies are shown using the data from Gastaldello et
al. (2007), Merten et al. (2015), and Umetsu et al. (2016). Gastaldello et al.
(2007) presented radial mass profiles for 16 relaxed galaxy groups. Using the
X-ray data from Chandra and XMM-Newton, they solved the equation of
hydrostatic equilibrium to calculate the gravitating mass distribution.
$r_{\mathrm{N}}$ (columun 3) of their table 3, $c_{200}$ (columun 8), and
$M_{200}$ (columun 10) of their table 7 are used. Merten et al. (2015) have
determined the mass distribution of $X$–ray selected galaxy clusters from the
Cluster Lensing and Supernova Survey with Hubble (CLASH), spanning a redshift
range between $0.19$ and $0.89$. Data from HST are used to analyse both weak
and strong gravitational lensing and data from ground–based telescopes,
primarily Subaru Telescope, are utilised for weak–lensing analysis.
$\rho_{\mathrm{N}}$ (columun 2) and $r_{\mathrm{N}}$ (columun 3) of their
table 6, and $M_{200}$ (columun 6) and $c_{200}$ (columun 7) of their table 7
are used. Umetsu et al. (2016) provided comprehensive analysis of
strong–lensing, weak–lensing shear and magnification data for a sample of
$X$–ray–regular and high–magnification galaxy clusters at $0.19\lesssim
z\lesssim 0.69$ selected from CLASH. They constructed mass models of
individual clusters from a joint analysis of the full lensing constraints.
$M_{200}$ (columun 2), $c_{200}$ (columun 3), and $r_{\mathrm{N}}$ (columun 4)
of their table 2 are used. The object analysed in Umetsu et al. (2016)
overlaps with the objects analysed in Merten et al. (2015). We adopt the
updated values from Umetsu et al. (2016) whenever possible.

To ensure data quality, objects with unreliable estimation are excluded from
our sample. The exclusion criteria are as follows:

  1. 1.

Data points exhibiting uncertainties larger than the expected values are
excluded. 111 This criterion excludes the data of Salucci et al. (2012) and
Hayashi et al. (2017) already shown in figure 1.

  2. 2.

Spano et al. (2008) introduced artificial thresholds for the stellar mass-to-
luminosity ratios of the disc and bulge, $r_{\mathrm{N}}$, and
$\rho_{\mathrm{N}}$, to eliminate non-physical values. We exclude data
affected by these thresholds.

###  5.2 Comparison between theoretical scaling relations and observations

Figure 3: Comparisons between observations (symbols) and (a) the $c$–$M$
relation, (b) $r_{\mathrm{N}}$–$\rho_{\mathrm{N}}$ relation, (c)
$V_{\mathrm{max}}$–$r_{\mathrm{max}}$ relation, (d)
$\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ relation. The orange
symbols represent nearby dwarf galaxies from Oh et al. (2015). The yellow,
green, sky–blue symbols show bright galaxies by de Blok et al. (2008), Spano
et al. (2008), and Sofue (2016). The blue symbols are groups of galaxies from
Gastaldello et al. (2007). The purple and magenta symbols correspond to
clusters of galaxies from Merten et al. (2015) and Umetsu et al. (2016),
respectively.

Here, the theoretical scaling relations based on the $c$–$M$ relation of the
CDM model are compared with the observational data for four combinations of
parameters.

In figure 3a, the $c$–$M$ relation is compared with the observation from
dwarfs to clusters. The vertical axis shows $c_{200}$, and the horizontal axis
shows $M_{200}$. The black solid curve shows the fitting function proposed by
Moliné et al. (2023) (equation 18 and 22), and the symbols represent the
observed values of $c_{200}$ and $M_{200}$. The theoretical $c_{200}$ slowly
decreases as the mass increases. The theoretical $c$–$M$ relation is
consistent with the observations in the whole mass range. However, a
significant deviation is seen for the observed dwarf galaxies. It should be
noted that there is a significant ambiguity in the mass determination of dwarf
galaxies because the velocity dispersion is usually comparable level to the
rotation velocity in such systems. Therefore, spectroscopic observation with
higher resolution and dispersion and the development of mass models are
expected to constrain the $c$–$M$ relation models more accurately.

Observation and theory are compared on the
$r_{\mathrm{N}}$–$\rho_{\mathrm{N}}$ plane in figure 3b. To convert $c_{200}$
and $M_{200}$ to $r_{\mathrm{N}}$ and $\rho_{\mathrm{N}}$, we use the
following equations:

$\rho_{\mathrm{N}}=\frac{200\rho_{\mathrm{crit,0}}c_{200}^{3}}{3f_{\mathrm{N}}(c_{200})},$
(28)

and

$r_{\mathrm{N}}=r_{200}/c_{200},$ (29)

where

$r_{200}=\left(\frac{3}{4}\frac{M_{200}}{200\pi\rho_{\mathrm{crit,0}}}\right)^{1/3}.$
(30)

The observations indicate that $\rho_{\mathrm{N}}$ decreases as
$r_{\mathrm{N}}$ increases in agreement with the theoretical line. Dwarf
galaxies, which have the $r_{\mathrm{N}}$ of several $\mathrm{kpc}$, spread
around the theoretical curve due to the errors of their mass distribution. The
spiral galaxies are in good accordance with the theoretical prediction, with
radii from $1$ to $100\,\mathrm{kpc}$ and the scale densities range from
$10^{6}$ to $10^{8}\,M_{\odot}\,\mathrm{kpc^{-3}}$. On the other hand, a small
discrepancy is seen in the low mass and massive ends. One possible reason is
that some assumptions needed to estimate the mass of the stellar component,
such as maximum/minimum disk models or a mass–to–light ratio, to decompose the
rotation curve have a systematic effect on the parameters of the haloes. The
scale radii of galaxy clusters range from hundreds of $\mathrm{kpc}$ to
several thousands of $\mathrm{kpc}$ and the scale densities range from
$10^{5}$ to $10^{7}M_{\odot}\,\mathrm{kpc^{-3}}$.

Figure 3c shows the comparison on the $V_{\mathrm{max}}$–$r_{\mathrm{max}}$
plane. The vertical axis indicates $r_{\mathrm{max}}$, and the horizontal axis
corresponds to $V_{\mathrm{max}}$.

We convert $r_{\mathrm{N}}$ and $\rho_{\mathrm{N}}$ values to
$V_{\mathrm{max}}$ and $r_{\mathrm{max}}$ with equation (7) and

$V_{\mathrm{max}}=\sqrt{4\pi
G\rho_{\mathrm{N}}{r_{\mathrm{N}}}^{2}\frac{f_{\mathrm{N}}(x_{\mathrm{N}})}{x_{\mathrm{N}}}}.$
(31)

The theoretical curve shows that $r_{\mathrm{max}}$ increases as
$V_{\mathrm{max}}$ increases. For the lower mass range, $r_{\mathrm{max}}$ is
approximately proportional to $V_{\mathrm{max}}$ on a log-log scale. The
theoretical relation nicely reproduces the observed values for two orders of
magnitude in the maximum circular velocity or seven orders of magnitude in the
virial mass.

Figure 3d shows $\bar{\Sigma}(<r_{\mathrm{max}})$ as a function of
$V_{\mathrm{max}}$ as in figure 1. The theoretical relation nicely reproduces
the observations. According to the theoretical
$\bar{\Sigma}(<r_{\mathrm{max}})$–$V_{\mathrm{max}}$ relation, more massive
objects have larger characteristic mean surface density.

On all planes, there is a certain amount of dispersion on the observational
data set for the range of the galaxies and dwarf galaxies. This may be caused
by using NFW profiles when deriving the dark matter halo profiles, rather than
using a profile with a core more suitable to reproduce the observed profiles
(see section 1). Until now, we adopted the NFW profile fit for comparing with
the results of cosmological $N$–body simulations. In the next section, we
extend the $c$–$M$ relation to cored profiles and compare with observations
fitted with cored profiles.

##  6 Theoretical scaling relation after the cusp–to–core transition

In this section, we construct a scaling relation for cored haloes based on the
$c$–$M$ relation. In previous sections, we compared $\Lambda$CDM model
predictions with observations assuming the dark matter density profile to have
a central cusp. However, observations of disc galaxies, dwarf galaxies, and
low surface brightness disc galaxies generally indicate the existence of a
core in the centre of the dark matter haloes (see section 1). Most scaling
relations hold for the cored density distribution. While several solutions to
solve the cusp–core problem are introduced in section 1, we especially focus
on the baryonic solution in this section. We stand on the hypothesis that dark
matter haloes formed primarily with a cuspy density distribution and then some
dynamical processes lead to the formation of a core at the centre of the dark
matter haloes. We simplify the complicated physical processes and provide a
model to convert the parameters of a cuspy profile, namely, the scale density
and the scale radius, to the parameters of a core profile, the central density
and the core radius.

###  6.1 Cusp–to–core transition model

In this section, our cusp–to–core transition model is described. We adopt the
Burkert profile as a core profile and the NFW profile as a cusp profile. In
order to enhance applicability, we deliberately omit specific physical
processes responsible for the cusp–to–core transition when constructing a
model. Ogiya et al. (2014) proposed a crude transition model that can be
analytically solved, assuming mass conservation and that only the inner part
is changed during the cusp–to–core transition, considering $r_{200}\gg
r_{\mathrm{N}},\,r_{\mathrm{B}}$. They had $r_{\mathrm{B}}=\eta
r_{\mathrm{N}}$ and $\rho_{\mathrm{B}}=\zeta\rho_{\mathrm{N}}$ where
$\eta\approx 1.24$ and $\zeta\approx 0.526$ and regard these as $\eta\approx
1$ and $\zeta\approx 1$. However, in practice,
$c_{200}=r_{200}/r_{\mathrm{N}}$ is on the order of ten and cannot be regarded
as infinity. Therefore, the actual coefficients $\eta$ and $\zeta$ depend on
$c_{200}$ and a more careful approach is required.

We impose the same assumptions as Ogiya et al. (2014), however, adopt more
sophisticated conditions to formulate a transition model. First, only the
density distribution around the centre is affected, and the one in the
outskirts does not change. Second, the virial mass $M_{200}$ of a dark matter
halo is conserved. From the first condition, equating the density distribution
of the NFW profile and the Burkert profile at $r=r_{200}$, we have

$\frac{\rho_{\mathrm{N}}}{\frac{r_{200}}{r_{\mathrm{N}}}\left(\frac{r_{200}}{r_{\mathrm{N}}}+1\right)^{2}}=\frac{\rho_{\mathrm{B}}}{\left(\frac{r_{200}}{r_{\mathrm{B}}}+1\right)\left(\frac{r_{200}^{2}}{{r_{\mathrm{B}}}^{2}}+1\right)}$
(32)

As a requirement that meets the second condition, we adopt the conservation of
virial masses, in other words, we equate equations(2) and (10) at $r=r_{200}$
and yield

$\rho_{\mathrm{N}}r_{\mathrm{N}}^{3}f_{\mathrm{N}}\left(\frac{r_{200}}{r_{\mathrm{N}}}\right)=\rho_{\mathrm{B}}r_{\mathrm{B}}^{3}f_{\mathrm{B}}\left(\frac{r_{200}}{r_{\mathrm{B}}}\right).$
(33)

Eliminating $\rho_{\mathrm{N}}$ or $\rho_{\mathrm{B}}$ from equations (32) and
(33) and substituting $r_{\mathrm{B}}=\eta r_{\mathrm{N}}$, we obtain

$g(\eta)-\frac{f_{\mathrm{N}}(c_{200})}{f_{\mathrm{B}}(c_{200}/\eta)}=0,$ (34)

where

$g(\eta)\equiv\frac{(c_{200}+\eta)(c_{200}^{2}+\eta^{2})}{c_{200}(c_{200}+1)^{2}}.$
(35)

We solve equation (34) for $\eta$ numerically at each $c_{200}$ value of dark
matter haloes. Here, note that $r_{200}$ is determined from a given
combination of $r_{\mathrm{N}}$ and $\rho_{\mathrm{N}}$ using equation (2) and
(5). Then substituting $r_{\mathrm{N}}$ into equation (32) or equation (33),
we obtain $\rho_{\mathrm{N}}$.

In this way, we can derive the value of $\rho_{\mathrm{N}}$ and
$r_{\mathrm{N}}$ from a given combination of $\rho_{\mathrm{B}}$ and
$r_{\mathrm{B}}$.

###  6.2 Theoretical scaling relations for cored haloes

Combining the prediction for the cuspy haloes inferred from the cosmological
$N$–body simulation with our cusp–to–core transition model, we can predict the
theoretical relation for dark matter haloes which have experienced
cusp–to–core transitions. Here, we convert the properties of cuspy haloes
described by the theoretical $c$–$M$ relation to the properties of cored
haloes using our cusp–to–core transition model. Specifically, we derive
$\rho_{\mathrm{N}}$ and $r_{\mathrm{N}}$ from the theoretical $c$–$M$ relation
and convert them into $\rho_{\mathrm{B}}$ and $r_{\mathrm{B}}$ using our
cusp–to–core transition model. Then we recover $c_{\rm{V}}$ and
$V_{\mathrm{max}}$ from $\rho_{\mathrm{B}}$ and $r_{\mathrm{B}}$ and fit them
with equation (18). The fitting function of the properties of cored haloes
after this transition, in other words, the theoretical
$c_{\rm{V}}$–$V_{\mathrm{max}}$ relation for cored haloes from dwarf galaxies
to clusters of galaxies is given by equation (18) with
$C=c_{\mathrm{V}}(V_{\max})$, $S=V_{\max}/\mathrm{km}\,\mathrm{s}^{-1}$ and

$c_{0}=4.925\times 10^{4}\quad\mathrm{and}\quad
a_{i}=\\{-0.8505,0.2434,-0.02314\\}.$ (36)

We also perform a fitting for the median of $c_{200}$ against $M_{200}$ for
cored haloes from dwarf galaxies to clusters of galaxies. The concentration
for the Burkert profile is defined by

$c_{\rm{B}200}=\frac{r_{200}}{r_{\rm{B}}}.$ (37)

Adopting equation (18) as fitting function with $C=c_{200}(M_{200})$ and
$S=M_{200}/M_{\odot}$, coefficients become

$c_{0}=41.21\;\mathrm{and}\;a_{i}=\\{-0.01431,-0.008294,0.0003589\\}.$ (38)

The error of this parameterization is less than 3 per cent. The
$c_{\rm{V}}$–$V_{\mathrm{max}}$ relation for cored haloes of dwarf galaxies of
MW–sized hosts is given by equation (18) with $C=c_{\mathrm{V}}(V_{\max})$,
$S=V_{\max}/\mathrm{km}\,\mathrm{s}^{-1}$ and

$c_{0}=2.251\times 10^{4}\;\mathrm{and}\;a_{i}=\\{7.490,-10.02,3.220\\}.$ (39)

The fitting result is shown in the magenta line in figure 2. We also perform a
fitting for the median of $c_{200}$ against $M_{200}$ for cored haloes of
dwarf galaxies. Assuming equation (18) for fitting function with
$C=c_{200}(M_{200})$ and $S=M_{200}/M_{\odot}$, coefficients become

$c_{0}=-2.386\times
10^{2}\;\mathrm{and}\;a_{i}=\\{-0.4412,0.05455,-0.002114\\}.$ (40)

These coefficients are also listed in table 3 for a quick look.

We compare the $c$–$M$ relation for cuspy haloes and cored haloes in figure 4.
In the left column, the characteristic mean surface density inside a fraction
of $r_{\mathrm{max}}$ against $M_{200}$ is shown. The characteristic mean
surface density can be defined for both a cusp and a core profile, while the
central surface density can only be defined for core profiles. Moreover, the
characteristic mean surface density can be easily inferred from observables
compared to the inner slope of the dark matter profile since it requires only
the determination of the mass inside a certain radius and does not need the
recovery of a full-density profile. These enable us of robust discussions
using the characteristic mean surface density. The right column shows
three–dimensional density at a fraction of $r_{\mathrm{max}}$ against
$M_{200}$. The theoretical scaling relation for cuspy haloes is shown as the
solid line. The theoretical scaling relation for core haloes after the
cusp–to–core transition is shown as the dashed lines. The top two panels of
figure 4 show $\bar{\Sigma}(<r_{\mathrm{max}})$ against $M_{200}$ and
$\rho(r_{\mathrm{max}})$ against $M_{200}$, respectively. The difference
between the theoretical scaling relation of the cuspy profile and the one for
the cored profile is pretty small. To tell the difference between the cusp and
the core profile, we must see the inner density distribution since the density
distributions in the outer part are almost the same. The two lower rows in
figure 4 are similar to the upper row but at radii $Xr_{\mathrm{max}}$ with
$X<1$. The smaller the $X$ value is, the larger the difference between the
density of a cuspy halo and a cored halo becomes. Then cusps and the cores can
be distinguished using $\bar{\Sigma}(<Xr_{\mathrm{max}})$ or
$\rho(Xr_{\mathrm{max}})$. Although we examined various $X$, the results
reported here are only for $X=1,0.1$ and $0.01$.

The top panels in figure 4 indicate that if the spatial resolution of mass
determination by observations is equal to or worse than $r_{\mathrm{max}}$, it
is impossible to distinguish between the core profile and the cusp profile of
the dark matter halo. In the middle panels, the density at $X=0.1$ has a
difference of one order of magnitude between the theoretical prediction for
the cusps and the cores, hence, should be easy to distinguish observationally.
The density at $X=0.01$ has a difference of two orders of magnitude between
the theoretical prediction for the cusps and the cores.

The characteristic mean surface density $\bar{\Sigma}(<Xr_{\rm{max}})$ and the
density $\rho(Xr_{\rm{max}})$ can be derived from the $c$–$M$ relation
(equations 18 and 22). However, for convenience, it is desirable to provide an
approximation that describes the dependence of the characteristic mean surface
density and the density on the $X$ and the mass of the dark matter halo. We
propose the fitting function for the characteristic mean density and the
density:

$\log_{10}[D(M_{200},X)]=\log_{10}[F(X)]+\sum_{i=0}^{4}a_{i}\left[\log_{10}{\left(\frac{M_{200}}{M_{\odot}}\right)}\right]^{i}.$
(41)

Here, $D(M_{200},X)$ indicates the characteristic mean density or the density
at $Xr_{\rm{max}}$ and $F(X)$ undertakes the $X$ dependency. $D(M_{200},X)$,
$F(X)$ and coefficients $a_{i}$ are listed in table $*$ ‣ 6.2. These ﬁts work
well in the halo mass $10^{7}M_{\odot}<M_{200}<10^{15}M_{\odot}$, its accuracy
is better than 11 per cent at all $M_{200}$ values within this range and
$X=\\{1,0.1,0.01\\}$.

Coefficients of equation (41).**footnotemark: $*$ $D$
profile{\dagger}{\dagger}footnotemark: ${\dagger}$ $F$ $a_{0}$ $a_{1}$ $a_{2}$
$a_{3}$ $a_{4}$ $\bar{\Sigma}(<Xr_{\rm{max}})\,M_{\odot}\,\rm{pc}^{-2}$ cusp
$\frac{1}{X^{2}}\frac{f_{\mathrm{N}}(Xx_{\mathrm{N}})}{f_{\mathrm{N}}(x_{\mathrm{N}})}$
$0.1491$ $-0.3566$ $0.1122$ $-0.008745$ $0.0002310$
$\bar{\Sigma}(<Xr_{\rm{max}})\,M_{\odot}\,\rm{pc}^{-2}$ core
${\frac{1}{X^{2}}\frac{f_{\mathrm{B}}(Xx_{\mathrm{B}})}{f_{\mathrm{B}}(x_{\mathrm{B}})}}$
$-0.5081$ $-0.1924$ $0.08468$ $-0.006754$ $0.0001825$
$\rho(<Xr_{\rm{max}})\,M_{\odot}\,\rm{kpc}^{-3}$ cusp
$\frac{(x_{\mathrm{N}}+1)^{2}}{X(Xx_{\mathrm{N}}+1)^{2}}$ $9.635$ $-1.211$
$0.1975$ $-0.01525$ $0.0003990$
$\rho(<Xr_{\rm{max}})\,M_{\odot}\,\rm{kpc}^{-3}$ core
$\frac{(x_{\mathrm{B}}+1)(x_{\mathrm{B}}^{2}+1)}{(Xx_{\mathrm{B}}+1)(X^{2}x_{\mathrm{B}}^{2}+1)}$
$8.4911$ $-0.9302$ $0.1505$ $-0.01187$ $0.0003174$ **footnotemark: $*$ Best-
ﬁtting parameters of the characteristic mean surface density-mass and the
density-mass relation for the model equation (41) are summarised.
{\dagger}{\dagger}footnotemark: ${\dagger}$The second column, ”profile”, shows
which profile the parametrisations are for. ”Cusp” means the parametrisation
is the original relation derived from the results of the cosmological $N$-body
simulations and ”core” means the function is derived by applying our
cusp–to–core transition model to the results of the simulations (see section 6
for detail).

Figure 4: The left–hand column shows characteristic mean surface density
inside $X$ times $r_{\mathrm{max}}$ against $M_{200}$. The right–hand column
shows three–dimensional density at $X$ times $r_{\mathrm{max}}$ against
$M_{200}$. The theoretical scaling relation for the NFW dark matter profile is
shown as the solid lines (equivalent to those used in figure 3). The
theoretical scaling relation after the cusp–to–core transition is shown as the
dashed line. The symbols represent the estimated properties of dark matter
haloes associated with dwarf galaxies, galaxies, groups and clusters of
galaxies by observation. Objects for which a core profile is adopted are
plotted as filled symbols, and for which a cusp profile is adopted are plotted
as open symbols.

###  6.3 Observational verification of the cusp–to–core transition

The observational data points are added to the theoretical prediction in
figure 4. The symbols represent the estimated properties of dark matter haloes
associated with dwarf galaxies, galaxies, groups and clusters of galaxies by
observations. The colour assignments are the same as in the previous sections.
For papers that published both cusp profile and core profile fitting results
and $\chi^{2}$-value for each fitting (de Blok et al., 2008; Spano et al.,
2008), the fitting result with the better $\chi^{2}$-value is adopted. Oh et
al. (2015) listed both cusp profile and core profile fitting results, however,
did not list $\chi^{2}$-value. They mentioned that the decomposed dark matter
rotation curves of most sample galaxies are well matched in shape to those of
core–like haloes. Therefore, we adopt core profile fits for all galaxies. For
objects that are fitted by the p–ISO, we convert the p–ISO parameters to the
Burkert profile parameters using the method described in Appendix Matching the
p-ISO to the Burkert profile. In this way, it becomes possible to compare the
observations with the parameters after adopting our cusp–to–core transition
model, which assumes the Burkert profile as a cored profile. Objects for which
the core profiles are adopted are plotted as filled symbols, and for which the
cusp profile is adopted are plotted as open symbols.

Here, we can obtain hints of the system size dependence of the cusp–to–core
transition from figure 4. In the bottom two panels, less massive objects, such
as dwarf galaxies and galaxies are located between the $c$–$M$ relation of the
cusp profile and the core profile, while massive objects, such as clusters of
galaxies are well placed over the $c$–$M$ relation of the cusp. This indicates
that massive galaxies and galaxy clusters did not undergo cusp–to–core
transitions and remain cuspy. Moreover, it is seen that the transition mass
between the dark matter halo with a cusp and halo with a core is at about
$10^{11}\,M_{\odot}$, with the transition range extending from
$10^{10}\,M_{\odot}$ to $10^{12}\,M_{\odot}$ on either side. This is
consistent with the statement of the classical model advocated by Dekel & Silk
(1986) that $10^{11}\,M_{\odot}$ is the mass at which the energy that
cumulative supernovae inject into the surrounding gas and the energy required
to carry the gas to infinity against the gravitational potential of the dark
matter halo are comparable. In the less massive region
($<10^{10}\,M_{\odot}$), there exist both dark matter haloes with a cusp and
that with a core. This indicates that among less massive galaxies, some
objects experience significant cusp–to–core transitions while others weakly
experience it or do not. Similar trends have also been seen in recent galaxy
formation simulations. For example, NIHAO (Tollet et al., 2016) and FIRE-2
(Lazar et al., 2020) have shown that for objects that have smaller
stellar–to–halo mass ratio, the mass distribution near the centre of the dark
matter halo turns from cusp to core. What is even more interesting is that
these simulations demonstrate that dark matter haloes in faint dwarf galaxies
are little affected by stellar feedback and that the mass distributions in the
central parts of dark matter haloes retain cuspy structures (Tollet et al.
(2016); Lazar et al. (2020)). On the other hand, analysis of faint dwarf
galaxies around the MW by hayashi_diversity_2020 (hayashi_diversity_2020,
hayashi_dark_2022) shows a mixture of cusp–like dark matter haloes and
core–like dark matter haloes in those galaxies. We will postpone discussing
the properties of dark matter haloes of ultra–faint dwarf galaxies to our
upcoming paper.

For MW subhaloes, the 2.5th to 97.5th (the 25th to 75th) percentile range of
theoretical $\bar{\Sigma}(<Xr_{\mathrm{max}})-M_{200}$ relation for cuspy
haloes and cored haloes do not overlap when $X<0.021$ ($X<0.15$). If the
characteristic mean surface density can be measured for $X<0.021$ ($X<0.15$),
it will be possible to determine the critical transition mass at which the
cusp–to–core transition occurs in dark matter haloes. For example, if there is
a $10^{8}\,M_{\odot}$ object ($r_{\mathrm{max}}\sim 0.6\,\mathrm{kpc}$ from
the $c$–$M$ relation for MW–satellite–like haloes; equation 39),
$Xr_{\mathrm{max}}$ corresponds to $0.01\rm{kpc}$ for $X=0.021$
($0.09\rm{kpc}$ for $X=0.15$). To resolve $0.01\rm{kpc}$ ($0.09\rm{kpc}$) at a
distance of $10\rm{Mpc}$ (local universe), the required angular resolution is
at least $0.21$ arcseconds ($1.9$ arcseconds). Better than angular resolution
$\sim 0.21$ arcseconds is expected to be achievable with the Atacama Large
Millimeter/submillimeter Array (ALMA), in the near future, the Next Generation
Very Large Array (ngVLA), the Thirty Meter Telescope (TMT), and other
telescopes.

In this paper, the least massive sample of a galaxy we have is $\sim
10^{8}\,M_{\odot}$. Nearby ultra–faint dwarf galaxies are suitable to test our
theoretical prediction on the less massive side ($<10^{8}\,M_{\odot}$).
Ultra–faint dwarf galaxies offer a more accurate determination of the central
density distribution in dark matter haloes because they are uninfluenced by
baryonic effects. By measuring the spatial distribution of stellar velocity
dispersion and combining it with dynamical equilibrium models, one can derive
the dark matter density profiles. The latest and upcoming telescopes and
observational devices, for example, Subaru Prime Focus Spectrograph (PFS) is
anticipated to achieve an angular resolution $\sim 30$ arcseconds at its peak
performance. If $X=0.021$ ($X=0.15$) is adopted as the criterion for
determining the critical transition mass, dwarf galaxies with dark matter halo
masses exceeding $10^{8}M_{\odot}$ within 90 kpc (600 kpc) from us are
considered favourable targets for this telescope.

##  7 Discussions

In this section, we discuss the functional form of the dark matter halo
profiles used in the observations, and the possible extension of our
cusp–to–core transition model.

###  7.1 Functional form of the dark matter halo profile

The archival observation data used in our analysis assumed the cuspy
functional form or the cored functional forms of mass density distributions.
Actually, there are other profiles suggested to describe the observational
mass distribution of haloes, however, studying them goes beyond the scope of
this paper. Further study to include a function with varying inner slope, such
as the Einasto profile, is needed.

###  7.2 Application of the cusp–to–core transition model

####  7.2.1 Inner slope of the density and transition mass

Our cusp–to–core transition model assumes either a cusp or a flat core. In
other words, expressing the mass density distribution near the centre as
$\rho(r)\propto r^{\gamma}$, only $\gamma=-1$ or $\gamma=0$ are considered. In
reality, however, there is a possibility that a cusp may change to a shallower
cusp or moderate core ($-1<\gamma<0$). The $c$–$M$ relation of moderately
transformed halo should be between the $c$–$M$ relation for the cusp profile
(solid line) and the one for the core profile (dashed line) in figure 4.
Actually, the observed galaxies lie between the prediction for the cuspy
haloes and for the cores around $10^{11}\,M_{\odot}$ that is the mass thought
to be at which the energy that cumulative supernovae inject into the
surrounding gas and the energy required to carry the gas to infinity against
the gravitational potential of the dark matter halo are comparable (Dekel &
Silk, 1986). These galaxies might have intermediate profiles ($-1<\gamma<0$)
which we did not take into account in this paper. It is not clear how the
deviation caused by forcibly fitting a halo with slopes of $-1$ or $0$ affects
the conclusions in this paper.

####  7.2.2 Stability of a core

The stability of a core is closely connected with the balance between the time
scale for the relaxation of a dark matter halo and the duration of intensive
star formation. We assume that the halo retains a core profile once a halo
experiences the cusp–to–core transition. This assumption is verified if below
two conditions are satisfied.

Firstly, the cusp–to–core transition should occur after the halo reaches the
equilibrium state. Oñorbe et al. (2015) showed that a core formed in the early
epoch can be erased as a halo continues to accrete matter and experience
central mergers using zoom–in cosmological hydrodynamic simulations.

Secondary, stellar feedback should occur recurrently. Ogiya & Mori (2011)
numerically experimented the dynamical response of dark matter haloes to the
variance of the gravitational potential induced by instantaneous gas removal
from galaxy centres. They found the removal of gas can decrease the inner
density, however, after a few dozen dynamical times, a shallow cusp recovers.
On the other hand, Ogiya & Mori (2014) found that if there is a recurrent
change of gravitational potential driven by the cycle of star formation and
stellar feedback, the core retains for a few dozen dynamical times. Many
previous hydrodynamic simulation studies reveal periodic star formation in the
early epoch of galaxy formation, supporting a recurrent change of
gravitational potential. Consequently, it is reasonable to assume the core
retains for a few dozen dynamical times.

####  7.2.3 Epoch of the cusp–to–core transition

While we transform the scaling relation for cuspy haloes at $z=0$ to construct
a scaling relation for cored haloes, the actual cusp–to–core transition should
have taken place in the past. The epoch of the cusp–to–core transition depends
on the formation epoch of the dark matter halo and the star formation history.
Since the star formation history is strongly affected by stellar feedbacks and
radiative feedbacks, the ambiguity of the star formation model becomes very
large depending on the accuracy of the stellar feedback model and the
radiative feedback model. Therefore, as the first step, our study simply
applies the cusp–to–core transition model to the $c$–$M$ relation at $z=0$. As
a next step, it is interesting to consider the epoch of the cusp-to-core
transition as a function of redshift, taking into account the formation of the
dark matter halo and the star formation history (Ogiya et al., 2014). The
validity of our model can be discussed if we compare these results with
observations of high-redshift galaxies.

##  8 Conclusion

In this paper, we investigate the theoretical and observational scaling
relations from dwarf galaxies to clusters of galaxies. We summarise our main
findings as:

  1. 1.

The observational scaling relations between properties of galaxy-sized dark
matter haloes (Burkert, 1995; Spano et al., 2008; Donato et al., 2009; Salucci
et al., 2012; Kormendy & Freeman, 2016) are consistent with the theoretical
scaling relation based on the $c$–$M$ relation of dark matter haloes formed in
cosmological $N$–body simulations. In other words, the observed scaling
relations seem to originate from the $c$–$M$ relation of the CDM haloes.

  2. 2.

Utilising the results of Phi–4096, a cosmological $N$–body simulation whose
resolution is high enough to enable us to discuss the properties of MW
satellite–sized haloes, we derived the $c$–$M$ relation for the MW–satellite
sized subhaloes at $z=0$. The fitting functions are shown in equation (18)
with equation (24) or (25).

  3. 3.

The $c$–$M$ relation derived from the Uchuu suite by Moliné et al. (2023)
reproduces observations going from dwarf galaxies to clusters of galaxies in
four parameters planes, the $c_{200}$–$M_{200}$ plane, the
$r_{\mathrm{N}}$–$\rho_{\mathrm{N}}$ plane, the
$V_{\mathrm{max}}$–$r_{\mathrm{max}}$ plane, and the
$V_{\mathrm{max}}$–$\bar{\Sigma}(r_{\mathrm{max}})$ plane. However, the
dispersion of parameters is broad from dwarf galaxies to massive galaxies.
This might be because of forcibly fitting the NFW profile to the cored haloes.

  4. 4.

We provide the scaling relation for cored haloes for a wide mass range from
dwarf galaxies to clusters as equation (18) with equation (36) and for the
MW’s satellite mass range as equation (18) with equation (39), converting the
$c$–$M$ relation for a cuspy profile to the one for a cored profile using the
cusp–to–core transition model proposed in this work.

  5. 5.

Comparing the theoretical scaling relation for cored haloes with observations,
we find: a) Less massive objects, such as dwarf galaxies and galaxies lie over
both the theoretical scaling relation for cored haloes and cuspy haloes, while
massive objects, i.e., groups of galaxies and clusters of galaxies are well
placed over the theoretical scaling relation for cuspy haloes. There is a
critical mass around $10^{11}\,M_{\odot}$ for the transition between the dark
matter halo with a cusp and that with a core. This may be consistent with the
statement of Dekel & Silk (1986) that $10^{11}\,M_{\odot}$ is the mass at
which the stellar feedback and the gravitational potential of the dark matter
halo are comparable. This indicates that on smaller scales, some objects
experience strong cusp–to–core transitions while others weakly experience it
or do not. On the other hand, massive objects do not undergo cusp–to–core
transitions and remain cuspy. b) If the characteristic mean surface density
for $X<0.021$ (to discriminate 95 per cent of dark matter halos formed in
cosmological $N$-body simulation) or $X<0.15$ (50 per cent) is determined, it
is possible to identify the critical transition mass at which the cusp–to–core
transition of the dark matter halo occurs. The latest and upcoming telescopes
and observational devices, such as Subaru PFS, ALMA, TMT, and ngVLA would
achieve this resolution for the nearby dwarf galaxies.

We thank the anonymous referee for the detailed review and comments, and
Tomoaki Ishiyama for providing the results of the cosmological $N$–body
simulation Phi–4096. We are grateful for the helpful comments from Kohei
Hayashi and Masashi Chiba regarding observational aspects, especially about
local dwarf galaxies, and for valuable perspectives on future observations
from Takuya Hashimoto. We are grateful to Yudai Kazuno for many essential
discussions and Moemi Tanuma for providing valuable material. This work was
supported in part by JSPS KAKENHI Grant Number JP24K00669, JP24K07085,
JP23KJ0280, JP22KJ0370, and JP20K04022, and by Multidisciplinary Cooperative
Research Program in CCS, University of Tsukuba.

## Matching the p-ISO to the Burkert profile

In this section, how the parameters of the p-ISO are altered to the parameters
of the Burkert profile is explained. In section 6.2, we provide a prediction
for cored halo with the Burkert profile in the figure 4. On the other hand,
some observational studies modelled dark matter haloes as the p-ISO. However,
our cusp–to–core transition model cannot be applied to the p-ISO since it is
proportional to $r^{-2}$ at the outer region while the NFW profile is
proportional to $r^{-3}$. In other words, the core radius and the central
densities of the p-ISO which satisfy two requirements of our cusp–to–core
transition model do not exist. Thereby, we scale the parameters of the p-ISO
to those of the Burkert profile. We describe how we translate
$r_{\mathrm{p-ISO}}$ and $\rho_{\mathrm{p-ISO}}$ into $r_{\mathrm{B}}$ and
$\rho_{\mathrm{B}}$ below.

Since the accuracy of observations decreases as we see a more distant radius
from the centre of a galaxy, reproducing the density distribution of the
central region is essential. Therefore, we consider the central density
expressed as a parameter of the p-ISO to be identical to the central density
of the Burkert profile:

$\rho_{\mathrm{p-ISO}}=\rho_{\mathrm{B}}.$ (42)

Then, we expect that masses enclosed within the core radius are the same.
Equating equation (10) and (14) at core radius of the p-ISO, and supposing
$r_{\mathrm{B}}=\xi r_{\mathrm{p-ISO}}$, we have

$r_{\mathrm{p-ISO}}^{3}f_{\mathrm{p-ISO}}(1)-r_{\mathrm{B}}^{3}f_{\mathrm{B}}\left(\frac{1}{\xi}\right)=0.$
(43)

Solving equation (43) for $\xi$ by numerical procedure, we obtain
$r_{\mathrm{B}}$ from given $r_{\mathrm{p-ISO}}$.

## References

  * Behroozi et al. (2013a) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013a, ApJ, 762, 109
  * Behroozi et al. (2013b) Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., Busha, M. T., Klypin, A. A., & Primack, J. R. 2013b, ApJ, 763, 18
  * Benítez-Llambay et al. (2019) Benítez-Llambay, A., Frenk, C. S., Ludlow, A. D., & Navarro, J. F. 2019, MNRAS, 488, 2387
  * Bullock & Boylan-Kolchin (2017) Bullock, J. S., & Boylan-Kolchin, M. 2017, Annual Review of Astronomy and Astrophysics, 55, 343
  * Burkert (1995) Burkert, A. 1995, ApJ, 447, L25
  * Burkert (2020) Burkert, A. 2020, ApJ, 904, 161
  * de Blok et al. (2001a) de Blok, W. J. G., McGaugh, S. S., Bosma, A., & Rubin, V. C. 2001a, ApJ, 552, L23
  * de Blok et al. (2001b) de Blok, W. J. G., McGaugh, S. S., & Rubin, V. C. 2001b, AJ, 122, 2396
  * de Blok et al. (2008) de Blok, W. J. G., Walter, F., Brinks, E., Trachternach, C., Oh, S.-H., & Kennicutt, R. C. 2008, AJ, 136, 2648
  * Dekel & Silk (1986) Dekel, A., & Silk, J. 1986, ApJ, 303, 39
  * Del Popolo (2023) Del Popolo, A. 2023, Physics of the Dark Universe, 40, 101203
  * Del Popolo & Le Delliou (2022) Del Popolo, A., & Le Delliou, M. 2022, Galaxies, 9, 123
  * Di Cintio et al. (2014) Di Cintio, A., Brook, C. B., Macciò, A. V., Stinson, G. S., Knebe, A., Dutton, A. A., & Wadsley, J. 2014, MNRAS, 437, 415
  * Diemand et al. (2007) Diemand, J., Kuhlen, M., & Madau, P. 2007, ApJ, 667, 859
  * Donato et al. (2009) Donato, F., Gentile, G., Salucci, P., Frigerio Martins, C., Wilkinson, M. I., Gilmore, G., Grebel, E. K., Koch, A., & Wyse, R. 2009, MNRAS, 397, 1169
  * El-Zant et al. (2001) El-Zant, A., Shlosman, I., & Hoffman, Y. 2001, ApJ, 560, 636
  * Fitts et al. (2017) Fitts, A., et al. 2017, MNRAS, 471, 3547
  * Fukushige & Makino (1997) Fukushige, T., & Makino, J. 1997, ApJ, 477, L9
  * Gastaldello et al. (2007) Gastaldello, F., Buote, D. A., Humphrey, P. J., Zappacosta, L., Bullock, J. S., Brighenti, F., & Mathews, W. G. 2007, ApJ, 669, 158
  * Gnedin & Zhao (2002) Gnedin, O. Y., & Zhao, H. 2002, MNRAS, 333, 299
  * Goerdt et al. (2010) Goerdt, T., Moore, B., Read, J. I., & Stadel, J. 2010, ApJ, 725, 1707
  * Gopika et al. (2023) Gopika, K., Desai, S., & Paranjape, A. 2023, MNRAS, 523, 1718
  * Hayashi & Chiba (2015) Hayashi, K., & Chiba, M. 2015, ApJ, 803, L11
  * Hayashi et al. (2020) Hayashi, K., Chiba, M., & Ishiyama, T. 2020, ApJ, 904, 45
  * Hayashi et al. (2022) Hayashi, K., Hirai, Y., Chiba, M., & Ishiyama, T. 2022
  * Hayashi et al. (2017) Hayashi, K., Ishiyama, T., Ogiya, G., Chiba, M., Inoue, S., & Mori, M. 2017, ApJ, 843, 97
  * Hui et al. (2017) Hui, L., Ostriker, J. P., Tremaine, S., & Witten, E. 2017, Phys. Rev. D, 95, 043541
  * Inoue & Saitoh (2011) Inoue, S., & Saitoh, T. R. 2011, MNRAS, 418, 2527
  * Ishiyama et al. (2013) Ishiyama, T., et al. 2013, ApJ, 767, 146
  * Ishiyama et al. (2021) Ishiyama, T., et al. 2021, MNRAS, 506, 4210
  * Klypin et al. (2016) Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340
  * Klypin et al. (2011) Klypin, A. A., Trujillo-Gomez, S., & Primack, J. 2011, ApJ, 740, 102
  * Kormendy & Freeman (2004) Kormendy, J., & Freeman, K. C. 2004, 220, 377
  * Kormendy & Freeman (2016) Kormendy, J., & Freeman, K. C. 2016, ApJ, 817, 84
  * Lazar et al. (2020) Lazar, A., et al. 2020, MNRAS, 497, 2393
  * Li et al. (2019) Li, P., Lelli, F., McGaugh, S. S., Starkman, N., & Schombert, J. M. 2019, MNRAS, 482, 5106
  * Madau et al. (2014) Madau, P., Shen, S., & Governato, F. 2014, ApJ, 789, L17
  * Merten et al. (2015) Merten, J., et al. 2015, ApJ, 806, 4
  * Moliné et al. (2017) Moliné, Á., Sánchez-Conde, M. A., Palomares-Ruiz, S., & Prada, F. 2017, ApJ, 4974
  * Moliné et al. (2023) Moliné, Á., et al. 2023, MNRAS, 518, 157
  * Moore (1994) Moore, B. 1994, Nature, 370, 629
  * Moore et al. (1999) Moore, B., Quinn, T., Governato, F., Stadel, J., & Lake, G. 1999, MNRAS, 310, 1147
  * Navarro et al. (1996a) Navarro, J. F., Eke, V. R., & Frenk, C. S. 1996a, MNRAS, 283, L72
  * Navarro et al. (1995) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1995, MNRAS, 275, 720
  * Navarro et al. (1996b) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996b, ApJ, 462, 563
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Navarro et al. (2010) Navarro, J. F., Ludlow, A., Springel, V., Wang, J., Vogelsberger, M., White, S. D. M., Jenkins, A., Frenk, C. S., & Helmi, A. 2010, MNRAS, 402, 21
  * Ogiya & Mori (2011) Ogiya, G., & Mori, M. 2011, ApJ, 736, L2
  * Ogiya & Mori (2014) Ogiya, G., & Mori, M. 2014, ApJ, 793, 46
  * Ogiya et al. (2014) Ogiya, G., Mori, M., Ishiyama, T., & Burkert, A. 2014, MNRAS, 440, L71
  * Oh et al. (2011) Oh, S.-H., Brook, C., Governato, F., Brinks, E., Mayer, L., de Blok, W. J. G., Brooks, A., & Walter, F. 2011, AJ, 142, 24
  * Oh et al. (2008) Oh, S.-H., de Blok, W. J. G., Walter, F., Brinks, E., & Kennicutt, R. C. 2008, AJ, 136, 2761
  * Oh et al. (2015) Oh, S.-H., et al. 2015, AJ, 149, 180
  * Oñorbe et al. (2015) Oñorbe, J., Boylan-Kolchin, M., Bullock, J. S., Hopkins, P. F., Kerěs, D., Faucher-Giguère, C.-A., Quataert, E., & Murray, N. 2015, ApJ, 454, 2092
  * Planck Collaboration et al. (2020) Planck Collaboration, et al. 2020, A&A, 641, A6
  * Prada et al. (2012) Prada, F., Klypin, A. A., Cuesta, A. J., Betancort-Rijo, J. E., & Primack, J. 2012, MNRAS, 423, 3018
  * Ragone-Figueroa et al. (2012) Ragone-Figueroa, C., Granato, G. L., & Abadi, M. G. 2012, MNRAS, 423, 3243
  * Read & Gilmore (2005) Read, J. I., & Gilmore, G. 2005, MNRAS, 356, 107
  * Salucci (2019) Salucci, P. 2019, Astronomy and Astrophysics Review, 27, 2
  * Salucci et al. (2012) Salucci, P., Wilkinson, M. I., Walker, M. G., Gilmore, G. F., Grebel, E. K., Koch, A., Frigerio Martins, C., & Wyse, R. F. G. 2012, MNRAS, 420, 2034
  * Sofue (2016) Sofue, Y. 2016, PASJ, 68, 2
  * Spano et al. (2008) Spano, M., Marcelin, M., Amram, P., Carignan, C., Epinat, B., & Hernandez, O. 2008, MNRAS, 383, 297
  * Spergel & Steinhardt (2000) Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760
  * Tollet et al. (2016) Tollet, E., et al. 2016, ApJ, 456, 3542
  * Umetsu et al. (2016) Umetsu, K., Zitrin, A., Gruen, D., Merten, J., Donahue, M., & Postman, M. 2016, ApJ, 821, 116
  * Weinberg (1978) Weinberg, S. 1978, Phys. Rev. Lett., 40, 223
  * Zwicky (1933) Zwicky, F. 1933, Helvetica Physica Acta, 6, 110