Assessment of Dark Matter Models Using Dark Matter Correlations across Dwarf Spheroidal Galaxies

Ahmad Borzou [email protected] EUCOS-CASPER, Department of Physics, Baylor University, Waco, TX 76798, USA

Abstract

The predicted size of dark matter substructures in kilo-parsec scales is model-dependent. Therefore, if the correlations between dark matter mass densities as a function of the distances between them are measured via observations, we can scrutinize dark matter scenarios. In this paper, we present an assessment procedure of dark matter scenarios. First, we use Gaia’s data to infer the single-body phase-space density of the stars in the Fornax dwarf spheroidal galaxy. The latter, together with the Jeans equation, after eliminating the gravitational potential using the Poisson equation, reveals the mass density of dark matter as a function of its position in the galaxy. We derive the correlations between dark matter mass densities as a function of distances between them. No statistically significant correlation is observed. Second, for the sake of comparison with the standard cold dark matter, we also compute the correlations between dark matter mass densities in a small halo of the Eagle hydrodynamics simulation. We show that the correlations from the simulation and from Gaia are in agreement. Third, we show that Gaia observations can be used to limit the parameter space of the Ginzburg–Landau statistical field theory of dark matter mass densities and subsequently shrink the parameter space of any dark matter model. As two examples, we show how to leave limitations on (i) a classic gas dark matter and (ii) a superfluid dark matter.

I Introduction

At the moment, the widely used cosmological model consists of a constant $\Lambda$ accounting for the vacuum energy, a cold collisionless dark matter (CDM), visible matter, and general relativity. Despite the popularity of $\Lambda$ CDM, there exist tensions between its predictions and observations. The discrepancy in the value of the Hubble constant [1], the impossibly early galaxy problem [2, 3], and the high-z quasar problem [4] can be mentioned as the challenges of the $\Lambda$ CDM model.

In the dark matter (DM) sector of $\Lambda$ CDM, at scales larger than 1 mega-parsecs, the observations are consistent with CDM. The cosmological model based on CDM provides a fairly accurate description of galaxy evolution, galaxy counts, and even galaxy morphology [5, 6]. Nevertheless, there are observations at the galactic scales that are hard to understand in the context of CDM. The observed mass densities of DM at the center of galaxies are (i) shallower [7, 8] and (ii) less steep [9, 10] than predicted by CDM cosmology. Therefore, CDM predictions of (i) the mass density of DM and (ii) the first derivative of the mass density are in disagreement with observations. Another class of observations that seem to contradict the predictions of CDM is related to the number of observed subhalos in galaxies such as our own Milky Way. While we have observed only $\sim$ 50 satellite or dwarf galaxies within the Milky Way, CDM predicts the number to be around 1000 [11]. Although some of these faint objects may have not been discovered, the difference between the observed and predicted counts is significantly high. Moreover, many of the observed satellite galaxies have total halo masses much less than the heaviest subhalos predicted by CDM. It is hard to understand why the heavier subhalos have failed to form galaxies while the less massive subhalos with lower efficiency in star formation have been observed [12].

Two classes of alternative scenarios have been introduced to solve the above mentioned small-scale problems. In the first category, baryonic feedback within the CDM framework accounts for the discrepancies [13]. In the second category, modified models of DM are suggested. Warm DM, self-interacting DM, degenerate fermionic DM, Bose–Einstein condensed models of DM, and superfluid DM are among the scenarios that are proposed for solving the small-scale problems [14].

Which one of these scenarios is better? By design, all of the proposals have a chance of explaining the mentioned small-scale problems while behaving more or less like CDM on larger scales. We need additional experiments or new analyses of the data from the existing experiments that can assess the DM scenarios in the domains that are independent of the mean density of dark matter. The present article is a contribution to the latter direction.

In this paper, we introduce a procedure to estimate the correlations between DM mass densities across dwarf spheroidal (dSph) galaxies and use that to lay limitations on DM models. We first learn the single-body phase-space density of the stars in the dSph galaxies from the motion of their stars. Next, we take a divergence of the Jeans equation and combine it with the Poisson equation to write the mass density distribution of DM in terms of the estimated single-body phase-space density of the stars. Estimation of DM mass density using the single-body phase-space density of stars has been reported in several publications. For example, see [15, 16] and the references therein. We use the estimated mass densities to explore the correlations between them. The estimated correlations shall be explained by any proposed DM scenario and results in a reduction of its parameter space. Since the correlations are independent of the mean DM mass distributions, as we show in Appendix A, the imposed restrictions are in addition to the requirement of explaining the small-scale observations.

There are two ways to use the measured correlations to leave limitations on a DM model. (i) In the case of sophisticated models, high resolution N-body simulations can be used to find the predicted correlations between mass densities across a small halo. The predicted correlations can then be compared with the measured correlations from observations. In this paper, we use the DM hydrodynamics simulations in the Eagle project to show that their CDM simulations produce no correlation between DM mass densities across small halos. (ii) In the case of simple analytic DM models, assuming that the halo exchanges DM with the surroundings to prevent the gravothermal catastrophe, we use the Ginzburg–Landau approach to construct the statistical field theory of the mass densities by expanding the free energy of the halo, in its general form, around the observed mean mass densities. The steadiness of the halo guarantees the smallness of the higher order terms. By neglecting these terms, we are able to straightforwardly derive an expression for the mass density correlations in terms of the coefficients of the free energy expansion. By comparing the correlations that are inferred from observations with the predicted correlations, one can place limitations on the coefficients of the free energy expansion. Since these coefficients are related to the underlying physics of a given model of DM, one can place bounds on the parameter space of the model.

As a showcase, we apply our procedure to the observations of the Fornax dSph galaxy, collected by the Gaia experiment. We observe no statistically significant correlation between DM mass densities that are apart by at least 100 (pc). We use this result to shrink the parameter space of (i) a classic DM gas and (ii) a superfluid DM as two examples of proposed DM models. It should be emphasized that the validity of these results depend on a few assumptions that are made to compensate the lack of observations of the z-components of positions and velocities of stars in dSph galaxies in the Gaia dataset. Therefore, this paper is more a presentation of a DM model assessing procedure than a full analysis of data. When the Gaia limitations are elevated, by, for example, integration of the results of other experiments, a reliable analysis will be possible. This full analysis is left for future works.

This paper is structured as follows. In Section II, we review the theoretical framework for deriving the correlations between mass densities of DM from the motion of stars in dSph galaxies. In the same section, we derive the theoretical form of the mass density correlations starting from the free energy of the model. In Section III, we estimate the DM correlations in a simulated small halo of the Eagle project. In Section IV, we retrieve the observations of the stars in Fornax dSph from Gaia and feed them into the theoretical framework and present the results. In the same section, a few DM models are assessed. A conclusion is drawn in Section V.

II Theoretical Layout

We begin with the widely used assumption that a given dSph galaxy has reached a steady state; see, for example, [17, 18, 19, 20]. In the case of the Fornax dSph, this assumption will be validated by data later in this article. Therefore, we start from the Jeans equation for the stars in the galaxy

\displaystyle\frac{1}{\rho^{*}}\partial_{j}\left(\rho^{*}\sigma^{2*ij}\right)+\partial^{i}\phi=0,

(1)

where an asterisk refers to the visible matter in the galaxy, $\phi$ is the gravitational potential, and the mass density and the dispersion velocity of the visible matter respectively read

	$\displaystyle\rho^{}=m^{}\int d^{3}vf^{*},$
	$\displaystyle\sigma^{2ij}=\int d^{3}vf^{}v^{i}v^{j},$		(2)

where $v$ is the velocity of stars, $f^{*}$ is the one-body phase-space density of stars, and $m^{*}$ is the granular mass of the stars, which will be canceled out later in the calculations.

After taking a divergence of Equation (1), and using the Poisson equation, the mass density of DM in the galactic halo reads

\displaystyle\rho=-\rho^{*}-\frac{1}{4\pi G}\partial_{i}\left(\frac{1}{\rho^{*}}\partial_{j}\left(\rho^{*}\sigma^{*ij}\right)\right),

(3)

where $G$ is the gravitational constant, and the mass $m^{*}$ will be canceled out in the second term. As far as the dwarf spheroidal galaxies are concerned, we can neglect the first $\rho^{*}$ on the right-hand side of this equation, and the DM mass density is approximately equal to the second term. Therefore, if we estimate the one-body distribution function $f^{*}$ from observations, the mass density of DM is known in terms of the positions in the galaxy.

II.1 Density Correlations from Observations

Assuming that DM mass density $\rho$ has been estimated, we first define the mean of DM mass density as

\displaystyle\bar{\rho}\equiv\frac{1}{\int d^{3}x}\int d^{3}x\,\rho(x),

(4)

where the integration is across a relatively large volume around the center of the halo. We define the DM mass density fluctuations as

\displaystyle\varphi(x)\equiv\frac{\rho(x)-\bar{\rho}}{\bar{\rho}}.

(5)

The correlations between the density fluctuations of DM, separated by distance $\Delta$ , can be estimated from observations

\displaystyle C(\Delta)=\frac{1}{\int d^{3}x\oint d\Omega_{\Delta}}\int d^{3}x\oint d\Omega_{\Delta}\,\varphi(\vec{x})\varphi(\vec{x}+\vec{\Delta}),

(6)

where $\oint d\Omega_{\Delta}$ means integration over all the directions of $\vec{\Delta}$ .

II.2 Density Correlations from DM models

We assume that the DM halo has reached a steady state after a long period of evolution. Moreover, we assume that the halo can exchange dark particles with its host and with the cosmic DM background and consequently obeys the statistics of a grand canonical ensemble with a partition function equal to ${\cal{Z}}=\sum_{E,N}\exp\left(-\beta\left(E-\mu N\right)\right)$ , where $E$ and $N$ are the energy and the number of particles of the halo, $\beta$ is the inverse of the temperature, and $\mu$ is the chemical potential due to the exchange of dark particles with the surrounding. The presence of $\mu$ in this statistics guarantees a state of minimum free energy and prevents the gravothermal catastrophe, which occurs in gravity dominant systems with a conserved number of particles. See Appendix A for more details. The partition function can be rearranged into the following form; see one such rearrangement for a trivial case in Appendix B,

\displaystyle{\cal{Z}}=\int D\rho\,e^{-\beta F[\rho]},

(7)

where $F[\rho]$ is the free energy functional, and $\int D\rho$ is the path integral over all possible DM density configurations.

Since the halo is in the steady state, i.e., $\frac{\delta F}{\delta\rho}|_{\bar{\rho}}=0$ , the free-energy functional can be expanded around the mean field $\bar{\rho}$ to write the partition function as (see [21, 22] for practical examples)

\displaystyle{\cal{Z}}=\int D\varphi\,\exp\left(-\frac{\beta}{2}\int d^{3}x\left(\varphi\,\frac{\delta^{2}F}{\delta\rho^{2}}|_{\bar{\rho}}\,\varphi+{\cal{O}}(\varphi^{3})+\varphi J\right)\right),

where $\frac{\delta^{2}F}{\delta\rho^{2}}|_{\bar{\rho}}=\gamma\partial^{2}+\nu^{2}$ is the inverse of the Greens’ function ${\cal{G}}^{-1}$ and $\gamma$ and $\nu^{2}$ are free parameters to be determined by the underlying physics. We have added an extra term $\varphi J$ to be set to zero later on and have set $D\rho=D\varphi$ after ignoring a normalization factor. In this paper, due to low data statistics, we ignore the higher order terms. In Appendix A, starting from the microscopic model of simple gases, the corresponding $\gamma$ and $\nu^{2}$ of that model are derived. Derivation of these coefficients in terms of the physics of other DM scenarios is left for the future. Since the higher order terms are ignored, the partition function can be calculated analytically and reads

\displaystyle{\cal{Z}}=\exp\left(\frac{\beta}{2}\int d^{3}x\int d^{3}yJ(\vec{x}){\cal{G}}(\vec{x}-\vec{y})J(\vec{y})\right),

(9)

where the normalization factor is dropped. To arrive at this equation, one needs to perform a linear transformation that diagonalizes the inverse Greens’ function in Equation (II.2). Next, the path integration can be separated into multiplications of independent one-dimensional Gaussian integrals with known answers.

The correlation between the DM fluctuations reads; see, for example, [23]

	$\displaystyle C(\Delta)$	$\displaystyle=$	$\displaystyle\beta^{-2}{\cal{Z}}^{-1}\frac{\delta^{2}{\cal{Z}}}{\delta J(\vec{x})\delta J(\vec{x}+\vec{\Delta})}\Big{\|}_{J=0}=\beta^{-1}{\cal{G}}(\Delta)$		(10)
		$\displaystyle=$	$\displaystyle 4\pi\left(\beta\gamma\Delta\right)^{-1}\exp\left(-\Delta/\sqrt{\beta\gamma/\beta\nu^{2}}\right).$		(10)

This correlation function should be compared with the one estimated from observations in Equation (6). To arrive at this equation, we note that the correlation is the weighted mean of the multiplication of the fluctuations $\langle\varphi(\vec{x})\varphi(\vec{x}+\vec{\Delta})\rangle\equiv{\cal{Z}}^{-1}\int D\varphi\,\varphi(\vec{x})\varphi(\vec{x}+\vec{\Delta})\,e^{-\beta F}$ . This expression of the correlation function is achieved through the second term in Equation (10) with the partition function given by Equation (II.2). If instead, we use the partition function in Equation (9) and take the two functional derivatives, the third term of Equation (10) would be the result. Since, by definition, $(\gamma\partial^{2}+\nu^{2}){\cal{G}}(\Delta)=\delta^{3}(\vec{\Delta})$ , the last term of Equation (10) has the form of a Yukawa potential.

It should be mentioned that, as far as the halo is in the steady state, the form of the partition function in Equation (II.2) is independent of the details of the DM model and the coefficients $\gamma$ , $\nu^{2}$ , and those of the higher order terms are the fingerprints of the model. In principle, if enough high precision data are collected, we would be able to estimate the coefficients up to sufficiently high orders and construct the true model of DM from data with no further assumption. The recipe would be to use the data to estimate the correlation function of Equation (6) as well as higher order correlation functions and then solve a system of equations to derive the coefficients of the free energy expansion. Unfortunately, the precision and low statistics of current experiments do not allow such a data-driven model building approach.

III Showcase I: Eagle Simulation

In this section, we treat the CDM simulation of the Eagle project [6, 24] as if it is the real data. Our goals are to (i) extract the DM correlations predicted by CDM such that we can test them later when they are understood via actual data and (ii) present the potentials of the theoretical method of Section II.

There are multiple simplifications when working with simulations. First, in actual data, we inevitably extract the mass density of DM from visible matter component. In simulation, the DM information is directly given. Second, an experiment like Gaia does not provide the z-components of the stars in dSph galaxies while these are known for all of the particles in simulations. Third, the systematic errors of star observations are large in an experiment like Gaia. When propagated to the DM sector, they become even larger. The systematic errors are absent in the simulations and we only need to deal with the statistical errors, which, as we see below, are quite small.

We use the particle data of the Eagle simulation with reference name RecalL0025N0752 at zero redshift and select a relatively small halo with GroupNumber = 1 and SubGroupNumber = 123 and halo mass of $\sim 10^{9}$ M_⊙. We retrieve the DM particle’s positions in the halo using the Python code snippets provided in [25]. In the following, we use the positions of DM particles to estimate the DM mass density $\rho$ and subsequently estimate the correlations. The whole process as well as statistical error estimation is implemented in a Python code that is publicly available at [26]. The code is intended to work with Gaia dataset, but the following procedures can be achieved with a minimal change.

We convert $(x,y,z)$ to dimensionless variables by dividing each of the components by the standard deviation of the corresponding component of all the DM particles in the halo. Next, we use the kernel density estimator [27, 28] to estimate $\rho$ in the three dimensional position-space. In this method, every particle contributes a Gaussian weight to a given point in the position-space. The sum of all the particles’ contributions to that point will be the probability of DM mass density there. In other words, the single-body phase-space density reads

\displaystyle\rho(\vec{x})=M{\cal{N}}\sum_{i}\exp\left(-\frac{|\vec{x}-\vec{x}_{i}|^{2}}{h^{2}}\right),

(11)

where $M$ is the total mass of the halo known through the simulation, ${\cal{N}}$ is the normalization factor and is set such that the position integral of the density is equal to $M$ , $i$ enumerates the DM particles in the halo, and $h$ is a free parameter to be determined such that the error is minimal. We use the implementation of the method in scikit-learn library, the neighbors class, and the KernelDensity method of the Python programming language. Ideally, we would like to explore the smallest lengths in a given halo, which requires small $h$ . Nevertheless, for a fixed number of particles in the halo, there is an optimal $h$ that minimizes the error in $\rho$ estimation but is greater than the ideal value. In general, better resolutions require higher number of particles in the simulation. To estimate the optimal $h$ , we use the GridSearchCV method of model_selection class of the scikit-learn library to explore the parameter space. See [29] for a description of the method. We find that the optimal $h$ is equal to $0.9$ , which, when scaled back to the position-space, is equivalent to $\sim$ 3 (kpc).

To estimate the statistical error, we use the estimated $\rho$ , which serves as the probability of finding DM particles at a given position, to draw a random sample dataset of the same size as the original one. Next, we estimate $\rho$ from the generated dataset by repeating the whole process above. We repeat the sampling and $\rho$ estimation until the standard deviation of the estimations of $\rho$ reaches steady values. The stable standard deviations are assigned as the statistical errors.

Finally, we use Equations (4) and (5) to compute the fluctuations $\varphi(\vec{x})$ , and substitute them into Equation (6) to estimate the correlations between fluctuations as a function of the distance $\Delta$ . The two-point correlation predicted by CDM can be seen in Figure 1. As the figure indicates, the statistical errors are relatively small and no significant correlation is predicted. Later, in Section IV, we analytically derive this result for a cold classic gas of DM, which is the underlying assumption of CDM simulations.

Refer to caption — Figure 1: We use the hydrodynamic simulations of the Eagle project and select a halo of mass $10^{9}$ M_⊙ at zero redshift. We estimate the DM mass density across the halo and use that to find the two-point correlations $C(\Delta)$ between DM mass density fluctuations as a function of the distance between the position of mass densities $\Delta$ . We observe no statistically significant correlation. The statistical errors are estimated and shown in the figure.

IV Showcase II: Gaia Data

In this section, we use the observations made by Gaia to learn the positions and velocities of the stars belonging to the Fornax dSph galaxy. After filtering and processing the data, we insert the results into the theoretical framework of Section II to derive the mass density distribution of DM and the correlations between them. The whole process is implemented in a Python code that is publicly available [26].

It should be emphasized that, in Section IV.3, we will compensate for the limitations of the Gaia experiment with an assumption whose validity is not known with certainty. Therefore, the results are only valid to the zeroth-order approximation for $f^{*}$ . In the future, one can integrate other observations with the Gaia data to remove the limitations. In that case, the assumptions made in this section are unnecessary and a thorough analysis will be possible.

IV.1 Selecting Stars

Billions of stars have been observed by Gaia. Among them, we need to extract those stars that belong to the Fornax dSph. For the Fornax dSph member selection, we use GetGaia, a Python code published in [30, 31], which performs a series of screening tasks on the dataset made available by Gaia. GetGaia starts with the stars within a large enough spherical region around the galaxy. Those stars with poor color and astrometric solutions are dropped. Next, those stars whose proper motion or parallax is not consistent with the bulk of the galaxy are filtered out. The selected stars have errors of less than 0.5 mas in parallax, less than 0.5 mas per year in proper motion, less than 0.1 mag in both integrated BP mean magnitude and integrated RP mean magnitude, and less than 0.01 mag in the G-band mean magnitude. Moreover, the stars with integrated BP/RP mean magnitudes of higher than 20.5 are removed. Since star selection and contamination removal in our work is entirely done by GetGaia, we refer to [30, 31] for a more detailed description of the selections above and the reasons behind them.

IV.2 Coordinate System

The calculations of Section II can be best carried out in a comoving Cartesian coordinate system $(x,y,z)$ whose origin is at the center of the Fornax dSph. The position and velocity of a star in this frame are

	$\displaystyle\vec{x}=\vec{r}-\vec{r}_{m},$
	$\displaystyle\vec{\dot{x}}=\vec{v}-\vec{v}_{m},$		(12)

where $\vec{r}$ refers to the position of the star in a spherical coordinate system $(r,\theta,\phi)$ attached to the sun with the polar and azimuthal angles defined as $\theta=\frac{\pi}{2}-\text{dec}$ and $\phi=\text{ra}$ , where dec and ra stand for the declination and right ascension of the equatorial coordinate system provided by Gaia. The subscript $m$ refers to the center of the galaxy in the spherical system, and $\vec{r}_{m}$ is defined as the mean of the position of the stars. Initially, the unit vectors of the co-moving coordinate system are set equal to the unit vectors of the spherical system at angles $\theta_{m}$ and $\phi_{m}$ that define $\vec{r}_{m}$

	$\displaystyle\hat{x}\equiv\hat{\theta}(\theta_{m},\phi_{m}),$
	$\displaystyle\hat{y}\equiv\hat{\phi}(\theta_{m},\phi_{m}),$
	$\displaystyle\hat{z}\equiv\hat{r}(\theta_{m},\phi_{m}).$		(13)

The components of the position and velocity of a star at $(r,\theta,\phi)$ , in the comoving coordinate system reads

	$\displaystyle x=\vec{r}\cdot\hat{x}\simeq D\,\hat{r}(\theta,\phi)\cdot\hat{\theta}(\theta_{m},\phi_{m}),$
	$\displaystyle y=\vec{r}\cdot\hat{y}\simeq D\,\hat{r}(\theta,\phi)\cdot\hat{\phi}(\theta_{m},\phi_{m}),$
	$\displaystyle\dot{x}=\vec{v}\cdot\hat{x}-\vec{v}_{m}\cdot\hat{x}\simeq-D\,\mu_{\delta}-\vec{v}_{m}\cdot\hat{x},$
	$\displaystyle\dot{y}=\vec{v}\cdot\hat{y}-\vec{v}_{m}\cdot\hat{y}\simeq D\,\mu_{\alpha*}-\vec{v}_{m}\cdot\hat{y},$		(14)

where we have used the orthogonality of the unit vectors to set $\vec{r}_{m}\cdot\hat{x}=\vec{r}_{m}\cdot\hat{y}=0$ , $D=147\pm 12$ (kpc) is the distance of the Fornax dSph from us [32], $\mu_{\delta}$ is the proper motion in declination, and $\mu_{\alpha*}$ is the proper motion in right ascension. We define $\vec{v}_{m}\cdot\hat{x}$ and $\vec{v}_{m}\cdot\hat{y}$ such that the mean of star velocities is zero, i.e., $\langle\dot{x}\rangle=0$ and $\langle\dot{y}\rangle=0$ . Moreover, some of the terms have been neglected knowing that for any of the stars, $r\simeq D$ , $\hat{\theta}\cdot\hat{\theta}_{m}\gg\hat{r}\cdot\hat{\theta}_{m}$ , $\hat{\theta}\cdot\hat{\theta}_{m}\gg\hat{\phi}\cdot\hat{\theta}_{m}$ , $\hat{\phi}\cdot\hat{\phi}_{m}\gg\hat{r}\cdot\hat{\phi}_{m}$ , $\hat{\phi}\cdot\hat{\phi}_{m}\gg\hat{\theta}\cdot\hat{\phi}_{m}$ . Figure 6 shows the histograms of these dot products for the stars in the Fornax dSph, confirming the validity of the inequalities. Finally, we rotate the comoving coordinate system around its z-direction until the new x and y components of the velocities are not correlated. This rotation results in the removal of the nondiagonal components of the dispersion velocity tensor, which we can confirm explicitly when the tensor is estimated from the data. The positions and velocities of the stars in the comoving coordinate system of the Fornax dSph are shown in Figure 8.

Unfortunately, the dSph galaxies are so far away that the z-components of the positions and velocities of their stars cannot be inferred using Gaia. To compensate for this limitation, we assume that the phase-space density $f^{*}$ has the following form $f^{*}\rightarrow f^{*}(x,\dot{x},y,\dot{y})f^{*}_{2}(z,\dot{z})$ . The probability distribution of stars in galaxies, i.e., $f^{*}$ , is often approximated as the Maxwell–Boltzmann distribution; see, for example, [33], which satisfies the assumed separation. Although this assumption is widely used in the community for the mean field analysis, its validity is not known for a next-order analysis such as the one presented in this paper. This restriction, hence the assumption, can be avoided when other datasets that have the z-components of the stars are combined with the Gaia dataset, when $f^{*}(x,\dot{x},y,\dot{y},z,\dot{z})$ can be estimated in the full six dimensions. We leave this for future work.

IV.3 Estimation of the Phase-Space Density

At this point, we would like to use the Cartesian components of the positions and velocities of the stars from the previous subsection and estimate $f^{*}$ of Equation (II). As mentioned above, we assume that $f^{*}$ is only a function of the four dimensional variable $\vec{q}\equiv(x,y,\dot{x},\dot{y})$ and the $z$ components have been integrated out. This assumption makes the following analysis a showcase of the method rather than a full analysis. The reason is that we are averaging the over- and underdensities along the $z$ -axis and then measuring the correlations between the z-averages across the x–y plane. Therefore, our discovery potential will be limited to special forms of DM substructures that are not washed out by integration over the $z$ components. In the future, when the z-information of the stars is collected from other experiments, the same procedure leads to thorough shrinkage of the parameter space of any proposed DM scenario, and the limitation on our analysis’ discovery potential will be elevated.

We convert $\vec{q}=(x,y,\dot{x},\dot{y})$ to dimensionless variables by dividing each of the components by the corresponding maximum value among all the stars. Next, we use the kernel density estimator [27, 28] to estimate $f^{*}$ in the four dimensional phase-space. The advantage of the kernel density estimator is that missing a few random member stars does not fundamentally change the estimation of the probability function but only increases its error. The method is the same as the one we used in Section III to extract the DM mass density probability distribution with the difference that it is used to estimate the probability of the four-dimensional dimensionless phase-space

\displaystyle f^{*}(\vec{q})=\sum_{i^{*}}\exp\left(-\frac{|\vec{q}-\vec{q}_{i^{*}}|^{2}}{h^{2}}\right),

(15)

where $i^{*}$ enumerates the stars of the dSph. We find that the optimal $h$ scaled back to the position-space is equivalent to $\sim$ 600 (pc) and scaled back to the velocities is equivalent to $\sim$ 38 (km/s) for Fornax dSph galaxy.

Figure 9 shows $f^{*}$ at the center of the Fornax dSph. The figure indicates that the velocity distribution of stars at that position is slightly different from the Maxwell–Boltzmann distribution, which is due to the fluctuations that we aim to compute. Nevertheless, the overall Maxwellian form of the distribution confirms our assumption in Section II that the Fornax dSph is in a steady state. We estimate one such distribution for every position on the x–y plane of the Fornax dSph. Having an estimation of single-body phase-space density, we insert it into Equation (II) to compute the density and dispersion velocity tensor at an arbitrary location $(x,y)$ . To carry out the derivatives of the two objects, we repeat the estimations of the density and the dispersion velocity tensor at $(x+dx,y),(x+2\,dx,y),(x,y+dy),(x,y+2\,dy)$ and use the finite difference method. Inserting the variables and their derivatives into Equation (3) and neglecting the first term on the right-hand side gives an estimation of the mass density of DM at $(x,y)$ . We repeat the same series of computations to estimate the DM mass densities over all the positions within a square region on the x–y plane of side length 800 (pc) centered at the origin of the Fornax dSph.

We carry out Equation (3) in a cylindrical coordinate system whose symmetry axis lies along the x-axis of the comoving coordinate system. Due to the rotation of x–y plane around the z axis, which was explained in Section IV.2, the x-axis is the symmetry axis of the system. Hence, the lack of observation of the z-component of the stars is compensated. The purpose of using a cylindrical over a spherical coordinate system is to avoid unnecessary additional assumptions about the x–y plane. More specifically, the coordinates of the cylindrical system, on the x–y plane, are $(y,\alpha,x)$ , where $x$ and $y$ are the components of the comoving Cartesian system and $\alpha$ is the polar angle on the y–z plane. In this system, the dispersion velocity tensor takes the following diagonal form $\sigma^{2i}_{~{}j}=(\sigma^{2}_{yy},\sigma^{2}_{yy},\sigma^{2}_{xx})$ , where the indices are raised and lowered by the following diagonal metric $g_{ij}=(1,y^{2},1)$ . Since the metric is different from unity, the partial derivatives in Equation (3) shall be replaced with the covariant derivatives whose connections are the Christoffel symbols of the metric.

We estimate both the statistical and the systematic errors of our estimation. To estimate the statistical error, we use the estimated $f^{*}$ to draw a random sample dataset of the same size as the original one. Next, we estimate $f^{*}$ from the generated dataset by repeating the whole process above. We repeat the sampling and $f^{*}$ estimation until the standard deviation of the estimations of $f^{*}$ reaches steady values. The stable standard deviations are assigned as the statistical errors. The systematic errors are estimated by propagating the errors of the variables in the Gaia dataset. We use the “uncertainties” package in the Python programming language [34] to carry out the propagation. We report that the statistical errors are negligible with respect to the systematic ones.

IV.4 Results

So far, we have found the DM mass density $\rho$ at every position within a square on the x–y plane with a side length of $\sim 800$ (pc) whose center is at $x=y=0$ . The estimated DM mass density distribution for the Fornax dSph can be seen in Figure 2.

We use Equations (4) and (5) to compute the fluctuations $\varphi(x,y)$ and substitute them into Equation (6) to estimate the correlations between fluctuations as a function of the distance $\Delta$ . The two-point correlations can be seen in Figure 3. We observe no meaningful correlations between the DM mass density fluctuations at distances $\Delta>100$ (pc). For $\Delta<100$ (pc), we observe that the correlation starts to deviate from zero and increases as $\Delta$ goes toward smaller distances. This correlation is induced by the kernel density estimator rather than being genuine. Unfortunately, given the current statistics of stars, the smoothing parameters $h$ are rather large. Therefore, $f^{*}$ and, as a result, the DM mass density at a given point $(x,y)$ are estimated using all the stars whose distances from $(x,y)$ are smaller than $h$ —the closer the stars are to the $(x,y)$ point, the more contribution to the estimation. The mass densities are expected to be correlated in distances sufficiently smaller than $h$ . At distances comparable to or larger than $h$ , no smoothing takes place and the estimated correlations are genuine.

At this point, we can use the results of Figure 3 to restrict the $\beta\gamma-\beta\nu^{2}$ parameter space of the free energy of Section II.2. We use the left-tailed $\chi^{2}$ method at a 5% significance level to place the limitations. The solid shaded region in Figure 4 is excluded for those models with $\nu^{2}>\gamma$ . Both the solid and the hatched shaded regions are excluded for models with $\nu^{2}\ll\gamma$ . The white region to the right of the Figure awaits future investigations. In the following, we discuss two DM models that belong to the two categories.

A Classic Weakly Collisional Gas of DM

We assume a gas of DM with weak gravitational collisions between its particles. The interaction between two dark particles is assumed to be a function of their distance $r$ and its potential energy is $u(r)$ . Defining $u_{0}\equiv\int d^{3}x\,u(r)$ and $u_{2}\equiv\int d^{3}x\,r^{2}u(r)$ , the coefficients of the free energy expansion in Equation (II.2) are

	$\displaystyle\beta\gamma=\beta\frac{\bar{\rho}^{2}}{m^{2}}u_{2},$
	$\displaystyle\beta\nu^{2}=\frac{\bar{\rho}}{m}\left(1+2\beta\frac{\bar{\rho}}{m}u_{0}\right),$		(16)

where $m$ is the mass of DM particles and the collisions are assumed weak enough that $\gamma<\nu^{2}$ . See Appendix A for the derivation.

To place a bound on $u(r)$ using the limitations of Figure 4, we note that at the center of small galaxies, the DM temperature to mass ratio is approximately equal to $10^{-6}~{}$ (K/eV) [35]. On the other hand, from Figure 2, the average mass density of DM at the center of the Fornax dSph is $\sim$ $0.5~{}(M_{\odot}/\text{pc}^{3})$ . Therefore, if DM mass is around 1 MeV, $\bar{\rho}/m$ $\sim$ $10^{59}(\text{pc}^{-3})$ . Substituting this number into Equation (IV.4), we can conclude that for a positive interaction, $\beta\nu^{2}$ is well above the excluded region in Figure 4. Therefore, a classical gas model of DM with no or positive $u_{0}$ can explain the nonexisting correlations of Figure 3. It should be emphasized that this simplistic model of DM by no means represents the sophisticated CDM model and the collisions due to baryonic feedback. However, the CDM model is still based on the evolution of a cold and classic gas. Therefore, the agreement between the analytic result presented here and the correlations predicted by the Eagle simulation of CDM presented in Figure 1 is not a surprise.

DM Models at the Critical Temperatures

Let us now focus on a scenario in which DM in small halos is in the superfluid state where the DM temperature is close to its critical point $T$ $\sim$ $T_{c}$ . It is known that close to the critical temperature, the correlation length diverges in the following form $\gamma/\nu^{2}$ $\sim$ $|T-T_{c}|^{-1}$ . As a result, the exponential in Equation (10) disappears and the correlation function takes the form of

\displaystyle C(\Delta)=4\pi(\beta\gamma\Delta)^{-1}.

(17)

One can observe that $\beta\nu^{2}$ , which was somewhat related to the number density of particles in Equation (IV.4), is absent in this equation. Using the $\chi^{2}$ test, we can conclude that only $\beta\gamma>10~{}(\text{pc}^{-1})$ is allowed in Figure 4.

V Conclusions

In this paper, we have used the simulations by the Eagle project and observations made by Gaia to estimate the mass density distribution of DM within the central part of (i) a simulated small halo and (ii) the Fornax dSph galaxy. For the two mentioned halos, we have computed DM mass density fluctuations as well as the correlations between them in the same regions. We have shown that the correlations between DM mass density fluctuations are not significantly different from zero when they are $>$ $100~{}$ (pc) apart in any of the two halos.

Our estimation of DM mass density correlations imposes restrictions on any proposed model of DM. Moreover, since correlations between density fluctuations are independent of the mean density, these limitations are in addition to those applied by observations of mass profiles of DM, through rotation curves, for example. We foresee two approaches for imposing restrictions on DM models using the estimation of mass density correlations from observations. In the first approach, provided that high-resolution N-body simulations exist, one can use Equation (6) to compute the predicted correlations using simulations and compare them with Figure 3. This route has been taken in this paper in Section III.

In the second approach, one writes down the general form of the statistical field theory of DM mass density. Assuming that the DM halo can exchange dark particles with its host and/or cosmic background, halo’s stability can be assumed, and gravothermal catastrophe is prevented by the induced chemical potential. Therefore, the free energy of the halo can be expanded around its steady state and the perturbation theory can be used to calculate the correlations between DM mass densities without knowing the underlying physics of the DM scenario. We have used Figure 3 to lay bounds on the coefficients of the free energy expansion. Since the coefficients are functions of the physics of the given DM model, the limitations can then be propagated into the model’s parameter space. We have used this approach to explore the parameter space of (i) a gas model of DM with weak collisions between its particles, and (ii) a superfluid DM at its critical temperature. The excluded regions of their parameter spaces have been presented.

The data analysis of this paper can be improved in several aspects. The radial velocities and distances of the member stars of dSph galaxies have been measured in other experiments. A combination of their datasets with Gaia observations can help avoid unnecessary assumptions. Such combinations might also help with identifying more member stars, which would subsequently increase the statistics and improve, or decrease, the smoothing parameter $h$ . The density of observed stars is not uniform across the halos. Especially, there are many more stars at the center of galaxies than at large distances from the centers. Hence, the analysis can enjoy an adaptive $h$ estimation with better resolution in regions where more stars have been observed. Since no public well-established implementation of this adaptive smoothing is available, we have left it for future work. Finally, new experiments targeting higher statistics and better precision can help to explore smaller length scales and tighten the limitations on DM models.

Appendix A Classic Gas Model of DM

The coefficients of the expansion of the free energy functional in Equation (7) are related to the underlying physics and serve as the fingerprints of a proposed DM scenario. In this section, we would like to use a simple example to demonstrate the map between the DM model and these coefficients.

We start by writing the total energy of a classic DM gas in the halo

\displaystyle E=\sum_{i}E_{i}+\sum_{ij}u_{ij},

(18)

where $i$ and $j$ enumerate dark particles, and $u_{ij}$ is the effective potential energy of the gravitational collisions between the particles. Since $\rho(x)=\sum_{i}m\delta^{3}(\vec{x}-\vec{x}_{i})$ , where $m$ is the mass of dark particles, the sum over the particles can be replaced by $\sum_{i}=\int d^{3}x\,\frac{\rho(x)}{m}$ . Therefore, in a grand canonical ensemble, the free energy functional reads

	$\displaystyle\beta F[\rho]$	$\displaystyle=$	$\displaystyle\beta F_{\text{free}}+\int d^{3}x\left(-\beta\frac{\rho}{m}\mu_{f}+J\rho\right)$		(19)
			$\displaystyle+\beta\int d^{3}x\,d^{3}x^{\prime}\frac{\rho(x)}{m}\frac{\rho(x^{\prime})}{m}u(\|\vec{x}-\vec{x}^{\prime}\|),$		(19)

where $\mu_{f}\equiv\left(\mu-m\phi\right)$ , and the terms related to an ideal gas read [36]

\displaystyle F_{\text{free}}=\int d^{3}x\bigg{(}\frac{\rho}{m}\text{Ln}\left(\frac{\rho}{m}\lambda^{3}\right)-\frac{\rho}{m}\bigg{)},

(20)

with $\lambda=\frac{h}{\sqrt{2\pi mkT}}$ where $h$ is the Planck constant. See Appendix B for the detailed derivation of $F_{\text{free}}$ and the second integral on the right-hand side of Equation (19).

At this point, we replace $\rho\rightarrow\bar{\rho}(1+\varphi)$ and separate the terms according to the powers of $\varphi$ . The zeroth-order terms appear as a normalization factor and will be canceled out later. The first and second order terms read

	$\displaystyle\beta F_{1}=\int d^{3}x\,\varphi\,\frac{\bar{\rho}}{m}\bigg{(}\text{Ln}\bar{\rho}+\text{Ln}\frac{\lambda^{3}}{m}-\beta\mu_{f}+2\beta\frac{\bar{\rho}}{m}u_{0}\bigg{)},$
	$\displaystyle\beta F_{2}=\int d^{3}x\,\bigg{(}\frac{1}{2}\frac{\bar{\rho}}{m}\varphi^{2}+\beta\frac{\bar{\rho}^{2}}{m^{2}}\varphi\int d^{3}x^{\prime}\,\varphi(x^{\prime})u(\|\vec{x}-\vec{x}^{\prime}\|)\bigg{)},$

where $u_{0}\equiv\int d^{3}x^{\prime}\,u(|\vec{x}-\vec{x}^{\prime}|)$ .

Our assumption regarding the steadiness of the halo implies that $\frac{\delta F[\varphi]}{\delta\varphi}|_{\varphi=0}=0$ , which means that the terms inside the integral in $F_{1}$ is zero. Hence, the steady condition requires that $\mu_{f}$ takes the following form

\displaystyle\mu-m\phi=2\frac{\bar{\rho}}{m}u_{0}+\frac{1}{\beta}\text{Ln}\frac{\bar{\rho}\lambda^{3}}{m}.

(22)

It should be noted that $\mu_{f}$ has to be negative for the free energy to have a minimum value corresponding to a stable steady state. If the number of particles in the system was conserved, the chemical potential $\mu$ on the left-hand side was absent and $\mu_{f}$ would be positive due to the negativeness of the gravitational potential. The latter is the origin of the gravothermal catastrophe in systems with dominant gravitational effects whose number of particles is not maintained by an external bath. The presence of the chemical potential on the left-hand side of Equation (22) allows us to assume that DM halos of dSph galaxies have negative chemical potentials that establish their stable steady states.

Assuming that $\varphi$ does not vary aggressively over the effective range of $u(|\vec{x}-\vec{x}^{\prime}|)$ , we use the so-called square gradient approximation, do an expansion up to the second order, and insert it into the integral in the second term of $F_{2}$

\displaystyle\int d^{3}x^{\prime}\,\varphi(x^{\prime})u(|\vec{x}-\vec{x}^{\prime}|)=\varphi(x)u_{0}+\frac{1}{2}\delta^{ij}\partial^{2}\varphi u_{2},

where $u_{2}\equiv\int d^{3}x^{\prime}\,(|\vec{x}-\vec{x}^{\prime}|)^{2}u(|\vec{x}-\vec{x}^{\prime}|)$ . Therefore, the effective free energy reads

\displaystyle\beta F[\varphi]=\frac{\beta}{2}\int d^{3}x\,\bigg{(}\gamma\varphi\partial^{2}\varphi+\nu^{2}\varphi^{2}\bigg{)},

(24)

where the $\gamma$ and $\nu^{2}$ are

	$\displaystyle\gamma=\frac{\bar{\rho}^{2}}{m^{2}}u_{2},$
	$\displaystyle\nu^{2}=\frac{\bar{\rho}}{m\beta}\left(1+2\beta\frac{\bar{\rho}}{m}u_{0}\right).$		(25)

Appendix B Free Energy Functional of an Ideal Gas

We begin with a sufficiently small volume interval $\delta V$ at position $\vec{x}$ of the halo with total energy $E(x)$ , the chemical potential $\mu(x)$ , that contains $N(x)=\frac{\rho(x)}{m}\delta V$ number of particles. We assume that this system is in the steady state such that its state’s probability is

\displaystyle{\cal{P}}_{x}=e^{-\beta\left(E(x)-\mu(x)N(x)\right)}.

(B1)

If the number of particles in the energy state $\varepsilon$ is $n_{\varepsilon}$ , the total energy of the system and total number of particles read

	$\displaystyle E(x)=\sum_{\varepsilon}n_{\varepsilon}\varepsilon+N(x)m\phi(x),$
	$\displaystyle N(x)=\sum_{\varepsilon}n_{\varepsilon}.$		(B2)

For the case of an ideal gas, no correlation does exist between different locations, making it a trivial case. Hence, the probability of finding the halo in a particular state reads

\displaystyle{\cal{P}}={\cal{P}}_{1}{\cal{P}}_{2}{\cal{P}}_{3}\cdots,

(B3)

where ${\cal{P}}_{i}={\cal{P}}_{x_{i}}$ . The partition function of the halo is the sum over all the probabilities

$\displaystyle{\cal{Z}}$	$\displaystyle=$	$\displaystyle\left(\sum_{N_{1},E_{1}}\sum_{N_{2},E_{2}}\cdots\right){\cal{P}}_{1}{\cal{P}}_{2}\cdots=\prod_{x}\sum_{N_{x}}\sum_{E_{x}}{\cal{P}}_{x}$	(B4)
	$\displaystyle=$	$\displaystyle\prod_{x}\sum_{N_{x}}e^{\beta N(x)\left(\mu(x)-m\phi(x)\right)}\frac{1}{N(x)!}$
		$\displaystyle\left(\frac{2\pi\delta V}{h^{3}}(2m)^{\frac{3}{2}}\int_{0}^{\infty}d\varepsilon\,\varepsilon^{\frac{1}{2}}e^{-\beta\varepsilon}\right)^{N(x)}$
	$\displaystyle=$	$\displaystyle\prod_{x}\sum_{N_{x}}e^{\beta N(x)\left(\mu(x)-m\phi(x)\right)}\frac{1}{N(x)!}\left(\frac{\delta V}{\lambda^{3}}\right)^{N(x)}.$

We note that the last two terms can be re-expressed in the following form

\displaystyle\frac{1}{N!}\left(\frac{\delta V}{\lambda^{3}}\right)^{N}=\exp\left(-N\,\text{Ln}\left(\frac{\lambda^{3}\rho}{m}\right)+N\right),

(B5)

where we have used $a=e^{\text{Ln}(a)}$ , the Stirling’s approximation for Ln $(N!)$ , and $\frac{N}{\delta V}=\frac{\rho}{m}$ . Hence, the partition function of the halo reads

	$\displaystyle{\cal{Z}}$	$\displaystyle=$	$\displaystyle\prod_{x}\int d\rho(x)\,\exp\left(-\delta V\frac{\rho}{m}\left(\text{Ln}\frac{\lambda^{3}\rho}{m}-1-\beta\mu_{f}\right)\right)$		(B6)
		$\displaystyle=$	$\displaystyle\int D\rho\exp\left(-\int d^{3}x\,\frac{\rho}{m}\left(\text{Ln}\frac{\lambda^{3}\rho}{m}-1-\beta\mu_{f}\right)\right),$		(B6)

where we have set $\sum_{x}\delta V\rightarrow\int d^{3}x$ .

Appendix C Additional Figures

This section presents a few additional figures.

References

Di Valentino et al. [2021] Di Valentino, E.; Mena, O.; Pan, S.; Visinelli, L.; Yang, W.; Melchiorri, A.; Mota, D. F.; Riess, A. G.; Silk, J. In the realm of the Hubble tension-a review of solutions. Classical and Quantum Gravity 2021, 38, 153001, https://doi.org/10.1088/1361-6382/ac086d.
Steinhardt et al. [2016] Steinhardt, C. L.; Capak, P.; Masters, D.; Speagle, J. The Impossibly Early Galaxy Problem. Astrophys. J. 2016, 824, 21, https://doi.org/10.3847/0004-637X/824/1/21
Yennapureddy and Melia [2019] Yennapureddy, M. K.; Melia, F. A comparison of the R_h=ct and $\Lambda$ CDM cosmologies based on the observed halo mass function. European Physical Journal C 2019, 79, 571, https://doi.org/10.1140/epjc/s10052-019-7082-z
Melia [2018] Melia, F. Model selection based on the angular-diameter distance to the compact structure in radio quasars. Europhysics Letters 2018, 123, 39001, https://doi.org/10.1209/0295-5075/123/39001
Vogelsberger et al. [2014] Vogelsberger, M.; Genel, S.; Springel, V.; Torrey, P.; Sijacki, D.; Xu, D.; Snyder, G.; Nelson, D.; Hernquist, L. Introducing the Illustris Project: Simulating the coevolution of dark and visible matter in the Universe. MNRAS 2014, 444, 1518–1547, https://doi.org/10.1093/mnras/stu1536.
Schaye et al. [2015] Schaye, J.; Crain, R.A.; Bower, R.G.; Furlong, M.; Schaller, M.; Theuns, T.; Dalla Vecchia, C.; Frenk, C.S.; McCarthy, I.G.; Helly, J.C.; et al. The EAGLE project: Simulating the evolution and assembly of galaxies and their environments. MNRAS 2015, 446, 521–554, https://doi.org/10.1093/mnras/stu2058.
Alam et al. [2002] Alam, S.M.K.; Bullock, J.S.; Weinberg, D.H. Dark matter properties and halo central densities. Astrophys. J. 2002, 572, 34–40, https://doi.org/10.1086/340190.
Oman et al. [2015] Oman, K.A.; Navarro, J.F.; Fattahi, A.; Frenk, C.S.; Sawala, T.; White, S.D.M.; Bower, R.; Crain, R.A.; Furlong, M.; Schaller, M.; et al. The unexpected diversity of dwarf galaxy rotation curves. MNRAS 2015, 452, 3650–3665, https://doi.org/10.1093/mnras/stv1504.
Simon et al. [2005] Simon, J.D.; Bolatto, A.D.; Leroy, A.; Blitz, L.; Gates, E.L. High-Resolution Measurements of the Halos of Four Dark Matter-Dominated Galaxies: Deviations from a Universal Density Profile. Astrophys. J. 2005, 621, 757–776, https://doi.org/10.1086/427684.
McGaugh et al. [2001] McGaugh, S.S.; Rubin, V.C.; de Blok, W.J.G. High-Resolution Rotation Curves of Low Surface Brightness Galaxies. I. Data. AJ 2001, 122, 2381–2395, https://doi.org/10.1086/323448.
Springel et al. [2008] Springel, V.; Wang, J.; Vogelsberger, M.; Ludlow, A.; Jenkins, A.; Helmi, A.; Navarro, J.F.; Frenk, C.S.; White, S.D.M. The Aquarius Project: The subhaloes of galactic haloes. MNRAS 2008, 391, 1685–1711, https://doi.org/10.1111/j.1365-2966.2008.14066.x.
Tollerud et al. [2014] Tollerud, E.J.; Boylan-Kolchin, M.; Bullock, J.S. M31 satellite masses compared to $\Lambda$ CDM subhaloes. MNRAS 2014, 440, 3511–3519, https://doi.org/10.1093/mnras/stu474.
Mashchenko et al. [2008] Mashchenko, S.; Wadsley, J.; Couchman, H.M.P. Stellar Feedback in Dwarf Galaxy Formation. Science 2008, 319, 174, https://doi.org/10.1126/science.1148666.
Feng [2010] Feng, J.L. Dark Matter Candidates from Particle Physics and Methods of Detection. Annu. Rev. Astron. Astrophys. 2010, 48, 495–545, https://doi.org/10.1146/annurev-astro-082708-101659.
Holmberg and Flynn [2000] Holmberg, J.; Flynn, C. The local density of matter mapped by Hipparcos. MNRAS 2000, 313, 209–216, https://doi.org/10.1046/ j.1365-8711.2000.02905.x.
Silverwood et al. [2016] Silverwood, H.; Sivertsson, S.; Steger, P.; Read, J.I.; Bertone, G. A non-parametric method for measuring the local dark matter density. MNRAS 2016, 459, 4191–4208, https://doi.org/10.1093/mnras/stw917.
Domcke and Urbano [2015] Domcke, V.; Urbano, A. Dwarf spheroidal galaxies as degenerate gas of free fermions. JCAP 2015, 2015, 002, https://doi.org/ 10.1088/1475-7516/2015/01/002.
Hayashi et al. [2016] Hayashi, K.; Ichikawa, K.; Matsumoto, S.; Ibe, M.; Ishigaki, M.N.; Sugai, H. Dark matter annihilation and decay from non-spherical dark halos in galactic dwarf satellites. MNRAS 2016, 461, 2914–2928, https://doi.org/10.1093/mnras/stw1457.
Strigari [2018] Strigari, L.E. Dark matter in dwarf spheroidal galaxies and indirect detection: A review. Rep. Prog. Phys. 2018, 81, 056901, https://doi.org/10.1088/1361-6633/aaae16.
Di Paolo et al. [2018] Di Paolo, C.; Nesti, F.; Villante, F.L. Phase-space mass bound for fermionic dark matter from dwarf spheroidal galaxies. MNRAS 2018, 475, 5385–5397, https://doi.org/10.1093/mnras/sty091.
Ohnesorge et al. [1991] Ohnesorge, R.; Löwen, H.; Wagner, H. Density-functional theory of surface melting. Phys. Rev. A 1991, 43, 2870–2878. https://doi.org/10.1103/PhysRevA.43.2870.
Hohenberg and Krekhov [2015] Hohenberg, P.C.; Krekhov, A.P. An introduction to the Ginzburg-Landau theory of phase transitions and nonequilibrium patterns. Phys. Rep. 2015, 572, 1–42, https://doi.org/10.1016/j.physrep.2015.01.001.
Tong [2017] Tong, D. Statistical Field Theory. Lecture Notes. 2017. Available online: http://www.damtp.cam.ac.uk/user/tong/sft.html (accessed on 20 July 2022).
Crain et al. [2015] Crain, R.A.; Schaye, J.; Bower, R.G.; Furlong, M.; Schaller, M.; Theuns, T.; Dalla Vecchia, C.; Frenk, C.S.; McCarthy, I.G.; Helly, J.C.; et al. The EAGLE simulations of galaxy formation: Calibration of subgrid physics and model variations. MNRAS 2015, 450, 1937–1961, https://doi.org/10.1093/mnras/stv725.
The EAGLE team [2017] The EAGLE team. The EAGLE simulations of galaxy formation: Public release of particle data. arXiv 2017, arXiv:1706.09899.
Borzou, A. [2019] Borzou, A. Python Code to Estimate the DM Mass Density in dSph Using Gaia Dataset, 2019. The Python Code Can Be Found Online. Available online: https://github.com/ahmadborzou/Gaia_dSph_DarkMatter (accessed on 20 July 2022).
Rosenblatt [1956] Rosenblatt, M. Remarks on Some Nonparametric Estimates of a Density Function. Ann. Math. Stat. 1956, 27, 832–837. https://doi.org/10.1214/aoms/1177728190.
Zanin Zambom and Dias [2012] Zanin Zambom, A.; Dias, R. A Review of Kernel Density Estimation with Applications to Econometrics. arXiv 2012, arXiv:1212.2812.
Heidenreich et al. [2013] Heidenreich, N.B.; Schindler, A.; Sperlich, S. Bandwidth selection for kernel density estimation: A review of fully automatic selectors. AStA Adv. Stat. Anal. 2013, 97, 403–433. https://doi.org/10.1007/s10182-013-0216-y.
del Pino et al. [2021] del Pino, A.; Fardal, M.A.; van der Marel, R.P.; Łokas, E.L.; Mateu, C.; Sohn, S.T. Revealing the Structure and Internal Rotation of the Sagittarius Dwarf Spheroidal Galaxy with Gaia and Machine Learning. Astrophys. J. 2021, 908, 244, https://doi.org/10.3847/1538-4357/abd5bf.
Martínez-García et al. [2021] Martínez-García, A.M.; del Pino, A.; Aparicio, A.; van der Marel, R.P.; Watkins, L.L. Internal rotation of Milky Way dwarf spheroidal satellites with Gaia Early Data Release 3. MNRAS 2021, 505, 5884–5895, https://doi.org/10.1093/mnras/stab1568.
McConnachie [2012] McConnachie, A.W. The Observed Properties of Dwarf Galaxies in and around the Local Group. AJ 2012, 144, 4, https://doi.org/10.1088/0004-6256/144/1/4.
Binney and Tremaine [2008] Binney, J.; Tremaine, S. Galactic Dynamics, 2nd ed.; United Kingdom: Princeton University Press: 2008.
[34] Lebigot, E.O. Uncertainties: A Python Package for Calculations with Uncertainties. Available online: http://pythonhosted.org/uncertainties/ (accessed on 20 July 2022).
Borzou [2021] Borzou, A. Estimation of the mass of dark matter using the observed mass profiles of late-type galaxies. JCAP 2021, 2021, 023. https://doi.org/10.1088/1475-7516/2021/08/023.
Lutsko [2006] Lutsko, J.F. First principles derivation of Ginzburg–Landau free energy models for crystalline systems. Phys. Stat. Mech. Its Appl. 2006, 366, 229–242. https://doi.org/10.1016/j.physa.2005.10.044.