††thanks: Corresponding author††thanks: Equal contribution.††thanks: Co-
corresponding author

#  _PyDaddy:_ A Python package for discovering
stochastic dynamical equations from time series data

Arshed Nabeel<EMAIL_ADDRESS>Ashwin Karichannavar  Shuaib Palathingal
Center for Ecological Sciences,
Indian Institute of Science, Bengaluru, India  Jitesh Jhawar  University of
Konstanz, Konstanz, Germany  Max Planck Institute of Animal Behaviour,
Konstanz, Germany  Danny Raj M  Department of Chemical Engineering,
Indian Institute of Science, Bengaluru, India  Vishwesha Guttal
<EMAIL_ADDRESS>Center for Ecological Sciences,
Indian Institute of Science, Bengaluru, India

###### Abstract

1\. Ecological dynamics, ranging from individual behaviours to ecosystem
dynamics, are inherently stochastic in nature. _Stochastic differential
equations_ (SDEs) are an important modelling framework to model dynamics with
intrinsic randomness. Here, we focus on the inverse question: If one has
empirically measured time series data from some system of interest, is it
possible to discover the SDE model that best describes the data?

2\. We present PyDaddy (Python library for Data Driven Dynamics), a toolbox to
construct and analyse interpretable SDE models based on time series data. We
combine traditional approaches for data-driven SDE reconstruction with an
_equation learning_ approach, to derive symbolic equations governing the
stochastic dynamics. The toolkit is presented as an open-source _Python_
library, and consists of tools to construct, visualise, and analyse SDEs.

3\. PyDaddy takes a time series as an input, and produces an interpretable,
analytical expression for the underlying SDE. The package includes
functionalities for visual examination of the stochastic structure of the
data, guided extraction of an analytical expression for the SDE, and diagnosis
and debugging of the underlying assumptions and the extracted model. Using
simulated time series datasets, exhibiting a wide range of dynamics, we show
that PyDaddy is able to correctly identify underlying SDE models. We
demonstrate the applicability of the toolkit to real-world data using a
previously published movement data of a fish school. The model discovered by
PyDaddy is consistent with the previous study.

4\. PyDaddy is an easy-to-use, open source package to discover SDEs that
govern the dynamics of time series data. The package, together with the
documentation and tutorial notebooks (which can be executed online even
without installation) are available at: https://pydaddy.readthedocs.io/

Data Driven Model Discovery, Langevin Dynamics, Self-organisation, Collective
motion, Mesoscale dynamics, Data Driven Dynamical Systems, Scientific Machine
Learning, Noise induced order

##  I Introduction

A central approach to modelling complex dynamical systems, across fields
including ecology, is via differential equations [1, 2]. Depending on the
importance of space and stochasticity in the system, these equations can be
ordinary, stochastic, or partial differential equations. Even when we begin
with a simple set of rules or interactions at a fine, local scale, one can
derive coarse-grained dynamical descriptions in the form of differential
equations. For example, in the context of ecological populations, the local
rules could be based on birth or death processes [3], movement and
interactions between organisms [4, 5, 6, 7], interactions of organisms with
the environment [8, 9], etc. These local interactions scale to shape the
dynamics at higher levels such as groups, populations or even ecosystems [10,
11, 12, 7, 9]. While such dynamical system based approaches are powerful and
have provided key biological insights, integrating empirical data with such
models continues to remain a challenge.

We are in the era of big data, where high-resolution data capturing the
dynamics of biological systems is increasingly available [13, 14]. Examples
span across scales of biological organisation, from individual (cells or
animals) movement trajectories [14] and group properties [5, 15, 16] to
population sizes [17, 18], fitness of populations [19] and ecosystem states
[20, 21, 22]. Importantly, data maybe available with high temporal
resolutions, sometimes even across space, thus opening avenues for better
integration between models and data.

To capture the dynamics of real-world systems accurately, one has to treat
state variables as stochastic rather than only accounting for the average
properties [8]. Not only is randomness intrinsic to biological systems, it may
also produce counter-intuitive effects that are absent in purely deterministic
systems [23, 11]. A suitable framework for analysing stochastic dynamics is
_stochastic differential equations (SDE)_ , which allows us to capture both
the deterministic and stochastic factors driving the dynamics. In the context
of integrating theory and data, we focus on the following central question:
can we construct stochastic dynamical models starting from observed time
series datasets?

Indeed, the answer is in the affirmative. Traditional approaches to stochastic
calculus, based on estimating so-called _jump moments_ , allow us to infer
stochastic differential equations from time series data [24, 25, 26, 27, 28].
More recently, techniques based on neural networks have also been developed
for SDE discovery [29, 30]. In the context of deterministic models, recently
developed _equation learning_ techniques allow us to discover simple,
interpretable differential equation models from time series data [31, 32].
Combining equation learning with the traditional methods of jump moments is
powerful, since it allows us to extract parsimonious and interpretable SDE
models directly from data [33, 34, 35]. However, these techniques are
scattered across physics and engineering literature and lack diagnostic tools.
Thus, there is a need for an end-to-end, user-friendly tool that takes in time
series data, applies these methods, outputs an SDE model and provides
diagnostic tools together with visualisations.

To address this gap, we introduce PyDaddy (Python library for Data Driven
Dynamics), a comprehensive and easy-to-use Python package for discovering
stochastic differential equations from empirical data. The package takes a
time series of a state variable (scalar or vector) as the input, and outputs
an interpretable SDE model. The package also implements a variety of
visualisation and diagnostic utilities. We illustrate the workflow of the
PyDaddy with the help of a real dataset from the collective animal motion
literature [15]. We also provide a detailed documentation along with ready-to-
use Python tutorial notebooks with a number of example datasets, to make the
package easy to use (see Data and Code availability section). We highlight
that the tutorial notebooks can be used without installation via Google Colab:
https://pydaddy.readthedocs.io/en/latest/tutorials.html.

##  II Discovering SDE models from time series data

###  II.1  Method: Combining jump moments and equation learning

Figure 1:  Sparse regression on the extracted instantaneous jump moments
reliably reconstructs governing equations from simulated data, for a wide
range of dynamics. The mathematical expressions for the actual and estimated
drift and diffusion functions are shown in the tables on the top (grey shaded)
and bottom (red shaded), respectively. Row _i_ shows a sample of the time
series $x(t)$ of the models Eqs 2-5, respectively, across four columns _(A-D)_
, while row _ii_ shows the corresponding histograms of $x(t)$. The arrows in
histograms in row _ii_ mark the deterministic stable equilibria of the
respective systems. Notice that the histograms in (A) and (B) are identical,
despite the underlying SDEs being very different. Similarly, the histograms in
(C) and (D) are similar. _(A)_ A simple Ornstein-Uhlenbeck process (Eq. 2),
with a linear drift and constant diffusion, giving rise to a unimodal
probability distribution of $x(t)$ that coincides with the deterministic
stable equilibrium. _(B)_ A system with a cubic drift and a quadratic
diffusion function (Eq. 3), also giving rise to a unimodal probability
distribution of $x(t)$. The deterministic part of this system has two stable
equilibria (two arrows in B-ii). However, due to the effect of stochasticity,
we observe only one mode of the time series which coincides with the unstable
deterministic equilibrium. _(C)_ A system with cubic drift and constant
diffusion (Eq. 4) giving rise to a bimodal probability distribution. The
bimodality here is explained entirely by deterministic equilibria, with noise
facilitating transitions between them. _(D)_ A system with a linear drift and
a quadratic diffusion (Eq.5), giving rise to a bimodal probability
distribution. The deterministic part of the system has only one stable
equilibrium at 0. Hence, the bimodality here is driven by stochasticity. Rows
_iii_ and _iv_ compare the ground-truth (black) drift and diffusion functions
with the estimated ones. The functions estimated using sparse regression are
shown as dark red dashed lines. In all cases, the proposed method accurately
recovers the drift and diffusion functions.

Consider a time series of some observed variable, such as population size or a
measure of group’s movement, denoted by $x$, captured at a high temporal
resolution. We wish to discover an equation governing the dynamics of $x$,
capturing both the deterministic and stochastic aspects. This can be done
using a stochastic differential equation of the form,

$\displaystyle\dot{x}=f(x)+g(x)\eta(t)$ (1)

(interpreted in the _Itô sense_) that describes the temporal dynamics of $x$.
Here, $f$ and $g$ are functions of $x$, and $\eta(t)$ is Gaussian white-noise
[36].

Eq. 1 describes how the rate of change of $x$ depends on the instantaneous
value of $x$. The average rate of change of $x$ over time is governed by $f$,
called the deterministic or _drift_ function. The fluctuation around the
average value is governed by $g^{2}$, called the stochastic or _diffusion_
function. When $g$ is a constant (i.e. independent of $x$), the strength of
the noise is the same for all values of $x$, and is called _additive_ or
_state-independent_ noise. When $g$ depends on $x$, the strength of the noise
depends on the instantaneous value of $x$, and is called _multiplicative_ or
_state-dependent_ noise [23].

Our goal is to discover the drift ($f$) and diffusion functions ($g^{2}$) from
the observed time series data. Specifically, we aim to find simple,
interpretable analytical expressions that describe $f$ and $g^{2}$, not just
the qualitative shape. Our approach to this problem involves two steps (see SI
Section S1 for the full mathematical details of the procedure):

  1. 1.

First, we extract the drift and diffusion components from the given time
series data, using the so-called _jump moments_ or the _Kramer-Moyal
coefficients_ [26]).

  2. 2.

Next, we use a technique based on sparse regression, sometimes referred to as
_equation learning_ , to find interpretable analytical expressions of the
extracted drift and diffusion functions [31, 33].

Conventional SDE estimation approaches usually stop with the first step. The
drift and diffusion functions are estimated as binned averages, based on user
input values of bin sizes and sampling time; although studies provide some
rules of thumb for selecting these parameter values, the results can
dramatically vary if incorrect values are given [27]. By combining the
conventional approach with sparse-regression-based equation discovery
techniques, our approach in PyDaddy precludes the need of arbitrary parameter
choices. In addition, the accuracy of the estimates is also considerably
improved, even when the available time series is short (see SI Section S3).

Modelling a given time series using an SDE is contingent on several
assumptions: for instance, the noise term $\eta$ is assumed to be uncorrelated
Gaussian noise. Therefore, the data-driven SDE discovery pipeline is
incomplete without sufficient diagnostic tests to make sure that the
assumptions are reasonably satisfied.

The following assumptions need to hold for the estimation to be accurate [37,
26]:

  * •

The residuals, defined as $r(t)=\dot{x}(t)-F(x(t))$, should be normally
distributed around 0. For vector time series, the residuals should also be
isotropic, i.e. the $\mathbf{x}_{i}$\- and $\mathbf{x}_{j}$-components of the
residuals should be uncorrelated for $i\neq j$.

  * •

The residuals should be uncorrelated across time.

  * •

_Pawula’s Theorem:_ Similar to the drift and diffusion components, one can
also estimate higher moments from the data. Pawula’s theorem states that, for
a drift-diffusion SDE, the third and higher moments should be zero [36]. In
practical scenarios, where we have a finite-length time series sampled at a
finite sampling time, the estimated higher moments will not go to zero in
general. In this case, alternative to Pawula theorem would be that
$K_{4}(x)\approx 3\cdot K_{2}(x)^{2}$, where $K_{2}$ and $K_{4}$ are the
second and forth Kramers-Moyal moments respectively. [26, 38]

In addition to this, there is also the issue of _model identifiability_ : how
can we be confident that the model we identify is the only (or at the very
least, the _‘best’_) model that describes the observed data? One way to
address this is to look at the self-consistency of the estimated model [27].
Suppose we estimate an SDE model using our approach. If we generate a
simulated time series using this estimated model, with identical length and
sampling interval as the original time series, then:

  * •

The simulated time series should have the same statistics—namely, histogram
and autocorrelation—as the original time series.

  * •

If we subsequently estimate a model from the simulated time series, this model
should match the originally estimated model.

###  II.2 Demonstration of the method from synthetic datasets

To demonstrate the estimation performance of the proposed method, we evaluated
the method on several simulated systems where the ground-truth governing
dynamics are known. We generated synthetic data from known SDEs, using several
different choices of $f$ and $g$. We then used our approach to estimate $f(x)$
and $g(x)$ from the simulated time series.

We consider the following two SDEs:

$\displaystyle\dot{x}$ $\displaystyle=-x+\sqrt{2}\cdot\eta(t)$ (2)
$\displaystyle\dot{x}$ $\displaystyle=x-x^{3}+2(\sqrt{1+x^{2}})\cdot\eta(t)$
(3)

Although Eq. 2 and 3 are quite different, the time series generated with these
equations as well as their histograms are very similar (Fig. 1 (A-i, A-ii,
B-i, B-ii): in fact, the steady-state distribution of $x(t)$ will be identical
for Eq. 2 and 3. Specifically, Eq. 2 has a single deterministic stable state
at $x^{*}=0$, which is reflected in the histogram as well (Fig. 1 (A-ii)). The
additive noise simply spreads the dynamics around the deterministic
equilibrium. On the other hand, the deterministic stable equilibria for Eq 3
are at $\pm 1$ but the mode of the histogram is at 0 (Fig. 1 (B-ii)). In this
case, the multiplicative noise-term completely alters the stability landscape,
annihilating the deterministic stable states and creating a new noise-induced
mode at $x=0$.

Similarly, consider another pair of SDEs:

$\displaystyle\dot{x}$ $\displaystyle=2x-3x^{3}+\frac{1}{2}\cdot\eta(t)$ (4)
$\displaystyle\dot{x}$ $\displaystyle=-x+\sqrt{2}(1-x^{2})\cdot\eta(t)$ (5)

Eq. 4 has two deterministic stable equilibria at $\pm\sqrt{2/3}$, and the
dynamics is spread around these equilibria by the additive noise (Fig. 1
(C-ii)). On the other hand, Eq. 5 has only one stable equilibrium at $x=0$;
yet, due to the effect of the noise term, the system exhibits the two modes in
the histogram away from the deterministic equilibrium (Fig. 1 (D-ii)).

Distinguishing between 2 and 3 (similarly 4 and 5) from either the time series
data or steady state distributions is challenging. Nevertheless, our approach
is able to distinguish these systems and extract the correct governing SDE in
each case (Fig. 1 rows iii and iv). We are also able to recover the analytical
expressions for the drift and diffusion functions almost exactly.

##  III The PyDaddy Package

Figure 2: Overall schematics of the PyDaddy package. PyDaddy takes as input
uniformly sampled 1D or 2D time series, and computes the drift and diffusion
components from the time series. Several functions are provided to visualise
data as time series or histograms, fit drift and diffusion functions, diagnose
whether underlying assumptions for drift-diffusion estimation are met, and to
export data.

We implement our technique for data-driven SDE discovery as a Python package,
PyDaddy (Python library for Data Driven Dynamics). In addition to the SDE
discovery functionality, based on combining sparse regression with the
computation of Kramer-Moyal coefficients (see SI Section 1), PyDaddy offers
many diagnostic and visualisation features that lets the user explore and
evaluate the results. PyDaddy takes either a scalar or 2-dimensional vector
data as the input, where the visualisations and diagnostics are most
effective. From a given uniformly sampled 1-D or 2-D time series, PyDaddy can
compute, visualise, and fit drift and diffusion functions. PyDaddy supports
time series with missing data. Diagnostic tools are provided to verify whether
the assumptions necessary for modelling the time series as an SDE are met.
Finally, PyDaddy can export data into a DataFrame or a CSV file. Fig. 2 shows
an overview of the package and its various functionalities.

For ease of use, we make the package and the associated tutorial notebooks
accessible online, even without installation, via Google Colab:
https://pydaddy.readthedocs.io/en/latest/usage.html#pydaddy-on-google-colab.

The package can also be installed on a local machine using PIP or Conda:

pip install pydaddy

or

conda install pydaddy -c tee-lab

PyDaddy features a single-command operation mode, which allows the user to run
a dataset through PyDaddy and generate a comprehensive HTML report. This can
be done using the command:

pydaddy <dataset-file-name>

The report contains interactive figures of the estimated drift and diffusion
functions (with default parameters and automatic threshold tuning—see below)
and diagnostic results about noise statistics and model self consistency.

PyDaddy also has a comprehensive Python API for advanced usage; which allows
the user to individually generate various plots, perform diagnostics, or fine-
tune the function fitting procedure.

###  III.1  SDE estimation workflow using a real dataset

Figure 3: The summary plot generated by PyDaddy when applied to schooling fish
data from the reference [15]. (A) Input time series (fish group polarisation),
a vector time series denoted by $\mathbf{m}=(M_{x},M_{y})$. One time step
corresponds to 0.12 s. (B) Autocorrelation functions of the time series,
showing that the autocorrelation decays with time. (C, D) Frequency histograms
of $\mathbf{m}$ and $|\mathbf{m}|$, depicting that the most likely state of
the fish group is schooling (i.e. highly polarised). (E, H) Preliminary
estimates of drift functions. (F, G, I, J) Preliminary estimates of diffusion
functions. The drift and diffusion functions shown here are bin-wise averages
of the jump moments (bin-wise averaging is done only for the visualisation).

We now demonstrate the PyDaddy estimation workflow, using a previously
published dataset of fish schooling [15]. Briefly, the dataset contains the
time series of the group polarisation vector, $\mathbf{m}=(M_{x},M_{y})$
(2-dimensional), for a group of 15 fish (_Etroplus suratensis_). The time
series is available at a uniform interval of 0.12 second, with many missing
data points. If the school of fish is highly ordered, with nearly all fish
moving in the same direction, the magnitude of the polarisation vector,
$|\mathbf{m}|$, will be closer to 1. On the other hand, if the group is
disordered, with fish moving in random directions, $|\mathbf{m}|$, will be
close to 0.

Authors [15] extracted an SDE that governs the group polarisation
$\mathbf{m}$, using the conventional approaches. Here, we will use this
dataset with PyDaddy to derive an SDE for $\mathbf{m}$.

####  III.1.1 Summary figure and estimation of drift and diffusion

The estimation workflow starts by creating a PyDaddy object dd, initialised
with a dataset.

dd = pydaddy.Characterize([x, y], t=0.12)

where x and y are Numpy arrays of the $m_{x}$ and $m_{y}$ time series
respectively. A summary figure is generated with basic statistics of the time
series such as the autocorrelation function and histograms of $\mathbf{m}$
(Fig. 3). The summary figure also shows a preliminary estimate of the drift
and diffusion functions, using binwise averaged conditional moments (see [15,
27] for more details on the binwise averaging method).

####  III.1.2 Fitting analytical expressions for drift and diffusion functions

Figure 4: Interpretable drift and diffusion functions estimated by PyDaddy
using the sparse regression method. The estimated intepretable functions
returned by PyDaddy are shown below each of the plots. $F_{1}$ and $F_{2}$ are
the components of the drift function, and $G_{11}$ and $G_{22}$ are the
components of the diffusion function. (The cross diffusion function
$G_{12}=G_{21}$ was estimated to be identically 0, and is not shown.)

The drift and diffusion plots above can be used to identify qualitative
features of the dynamics, such as the nature and positions of the
deterministic equilibrium points, nature of the state-dependent noise (i.e.
diffusion), presence or absence of cross diffusion, etc. PyDaddy can go a step
beyond and fit actual analytical expressions for the drift and diffusion
functions. This can be done using the fit() function of the dd object.

As mentioned in Section II.1, fitting is done using sparse regression. By
default, PyDaddy fits polynomial functions, but other function families can
also be fit using custom libraries. A _sparsification threshold_ , which
governs the number of terms in the fit (see SI Section S1B) can be either
chosen automatically using cross-validation (see SI Section S1C), or be set
manually.

For example,

dd.fit(’F1’, order=3, tune=True, plot=True)

fits a degree-3 polynomial (order=3) with sparsification threshold chosen
automatically (tune=True). The first argument names the function to be fitted:
F1 and F2 correspond to the two components of the drift; G11 and G22
correspond to diffusion, and G12 and G21 correspond to cross diffusion. The
plot=True argument tells the function to produce a cross-validation error plot
for a range of thresholds (see SI Fig. S1). The ideal threshold should strike
a balance between the error and number of terms in this plot. PyDaddy chooses
as high as a threshold as possible without increasing the cross-validation
error too much (see SI Section S1C for details).

The threshold can also be manually set by the user. For example,

dd.fit(’F1’, order=3, threshold=0.1)

produces a polynomial with degree up to 3, with only terms with coefficients
0.1 or higher.

Fig. 4 shows the drift and diffusion functions estimated by PyDaddy for the
example time series. In line with previous results [15], we find a linear
drift and a quadratic diffusion, suggesting that the group polarisation is
noise-induced. The cross-diffusion terms (not shown) were negligible and were
fit by the zero polynomial.

In the drift and diffusion functions in Fig. 4, coefficients with overlapping
error-bars can be combined (by replacing them with their averages), to get the
following vector SDE:

$\displaystyle\dot{\mathbf{m}}=-0.147\mathbf{m}+\sqrt{0.021+0.236(1-|\mathbf{m}|^{2})}\cdot\bm{\eta}(t)$
(6)

where we used ${\mathbf{m}}$ to represent the group polarisation in the
estimated equation. This equation is consistent with the result of [15].

We emphasize an important point here: although PyDaddy provides functionality
for automatic model selection, model selection works best when coupled with
some human input based on theory and understanding of the system. In
particular, we suggest that the maximum order of the polynomial, or more
generally the library of function to fit, needs to be specified based on
visual inspection of the drift/diffusion plots in conjunction with a
theoretical understanding of the system. The PyDaddy documentation [39]
contains an extended discussion about this aspect.

####  III.1.3  Diagnostics

Figure 5: Noise diagnostics figure generated by PyDaddy. (A) Distribution of
the residuals, which resembles a bivariate Gaussian with zero mean. The inset
matrix shows correlation matrix of the $x$ and $y$ components of the noise.
This should be close to the identity matrix. (B) Autocorrelation functions of
the residuals, with the black line marking the autocorrelation time and is
close to zero, as desired. (C, D) QQ-plots of the marginal distributions,
which nearly falls on a line of slope 1, again as desired by the assumptions
of the SDE model.

It is also essential to ensure that the assumptions involved in modelling the
time series as an SDE (SI Section S1E) are not violated. PyDaddy provides
diagnostic functions to verify that the assumptions are met. Diagnostic plots
in PyDaddy can be generated using

dd.noise_diagnostics()

Fig. 5 shows the noise diagnostic plots produced by PyDaddy. The diagnostic
plots show the distribution and autocorrelation of the residuals, and QQ-plots
of the marginals of the residuals against a theoretical Gaussian distribution.
These plots test for deviations of the residuals from an uncorrelated Gaussian
distribution (see Section II.1 for details).

In addition, PyDaddy can also perform a model self-consistency check (Section
S1E) using the following command:

dd.model_diagnostics()

This function uses the estimated SDE to generate a simulated time series, with
the same length and sampling interval as the original dataset. It then
estimates the drift and diffusion functions from the simulated time series,
using the same order and threshold used in the original estimation. Finally, a
summary figure is generated, comparing the histogram, drift and diffusion
functions of the original dataset to the simulated dataset (SI Fig. S4). If
the model is self-consistent, we expect a good match between the original and
the simulated time series.

##  IV Discussion

Here, we have presented a Python package PyDaddy—an easy-to-use tool for data-
driven discovery of dynamics of complex systems. Specifically, PyDaddy seeks
to find a data-driven stochastic differential equation that captures the
dynamics of the observed state variable. We combined the traditional
approaches of drift and diffusion estimation using jump moments, with
techniques based on sparse regression to recover SDE models from time series
data. The use of sparse regression helped us eliminate certain arbitrary
choices of parameters (such as bin-size and subsampling time) used in previous
methods [27, 28]. Equally importantly, we provide visual and diagnostic tools
to test the assumptions associated with data-derived functions, sparse fits,
as well as the noise.

We demonstrated the method with a time series data of schooling fish. However,
it is readily applicable to any ecological system, or more generally any
dynamical system, where sufficiently highly resolved and long data are
available. This raises the obvious question: how finely resolved and how long
a time series is necessary? As demonstrated in previous studies [27], the time
series must be at around an order of magnitude (or more) finely resolved than
the autocorrelation time of the state variable. When the time resolution is
comparable or more than the autocorrelation time, the dynamics of the system
can’t be captured; for such coarse time series, the estimated drift and
diffusion functions typically converge to a linear and upward-parabolic
functions, independent of the actual dynamics [28]. While we are not aware of
any studies that have looked at the minimum length required to estimate the
drift and diffusion reliably, based on our own numerical simulations of models
shown in this study, we speculate that the total length must be at least a few
orders of magnitude times the length of the autocorrelation time.

We expect many challenges, and hence opportunities for future work, when
applying these methods to complex ecological systems. First, despite the fact
that we have automated various steps, there are still a few choices the user
will make. These include the choice of the functions the user will chose to
fit and the sparsity threshold, since a less sparse function (i.e. one with
more terms) will always explain the variation of data better. We suggest that
one must choose the best fit functions using the _‘principle of parsimony’_ ,
and be guided by theory when possible. We have already alluded to the
importance of length of time series, and how it may impact on what we find as
best-fit functions. Secondly, in the current framework of SDEs, the structure
of noise is Gaussian and uncorrelated. This could be too simplistic since many
ecosystems are affected by different types of perturbations and stochasticity;
for example, a single large natural or human induced calamity is better
modelled as a shot noise rather than a Gaussian white noise [40]. Noise in
real systems may also be temporally correlated [41]. Finally, the key drivers
of the system may show seasonality (e.g. rainfall, nutrient flow), thus
violating the assumption of stationarity which is necessary for discovering
SDEs from time series data. Modifying and adopting the package described here
to account for realistic features presents interesting directions for future
work.

In summary, as the era of big data looms over much of biological sciences
including ecology, we hope that our presentation of a (relatively) easy to use
package that helps us characterise the governing dynamical equations from the
data will inspire usage of these relatively unknown methods in the field of
ecology and more broadly, to other complex dynamical systems.

## Data and Code Availability

The full code and documentation for the PyDaddy package is available at
https://github.com/tee-lab/PyDaddy and https://pydaddy.readthedocs.io
respectively.

We provide several tutorial notebooks to familiarise the user with the
package:

  * •

Getting started with PyDaddy

  * •

Getting started with vector data

  * •

Recovering SDEs from SDE-generated time series

  * •

Advanced function fitting

  * •

Exporting data

  * •

Fitting non-polynomial functions

The tutorial notebooks can be accessed here:
https://pydaddy.readthedocs.io/en/latest/tutorials.html

PyDaddy also comes with some sample datasets:

  * •

Real world data of fish schooling (taken from Jhawar et. al. [15])

  * •

Data generated from stochastic agent-based simulations, in 1-D and 2-D, based
on [27].

## Acknowledgements

We acknowledge support from Science and Engineering Research Board (VG),
CEFIPRA (VG), DST-FIST (VG), DST INSPIRE Faculty Award (DRM), Humboldt
Postdoctoral Fellowship (JJ) and MoE PhD Fellowship (AN abd SP). We thank
Ashrit Mangalwedhekar, Vivek Jadhav, Shikhara Bhat, Cassandre Aimon and
Harishankar Muppirala for comments on the package.

## Author Contribution

VG and DRM conceived the project. AN and AK developed the package, associated
methods and performed analysis with inputs from all coauthors. AN, DRM and VG
wrote the manuscript. All coauthors contributed to discussions and comments on
the manuscript.

## References

  * Strogatz [2018] S. H. Strogatz,  _Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering_ (CRC press, 2018).
  * Gotelli  _et al._ [2008] N. J. Gotelli _et al._ ,  _A primer of ecology_ , Vol. 494 (Sinauer Associates Sunderland, MA, 2008).
  * McKane and Newman [2004] A. J. McKane and T. J. Newman, Stochastic models in population biology and their deterministic analogs, Physical Review E 70, 041902 (2004).
  * Cheng  _et al._ [2014] H. Cheng, N. Yao, Z.-G. Huang, J. Park, Y. Do, and Y.-C. Lai, Mesoscopic interactions and species coexistence in evolutionary game dynamics of cyclic competitions, Scientific Reports 4, 1 (2014).
  * Yates  _et al._ [2009] C. A. Yates, R. Erban, C. Escudero, I. D. Couzin, J. Buhl, I. G. Kevrekidis, P. K. Maini, and D. J. T. Sumpter, Inherent noise can facilitate coherence in collective swarm motion, Proceedings of the National Academy of Sciences 106, 5464 (2009).
  * Lewis  _et al._ [2013] M. A. Lewis, P. K. Maini, and S. V. Petrovskii, Dispersal, individual movement and spatial ecology, Lecture Notes in Mathematics (Mathematics Bioscience Series) 2071 (2013).
  * Jhawar  _et al._ [2019] J. Jhawar, R. G. Morris, and V. Guttal, Deriving Mesoscopic Models of Collective Behavior for Finite Populations, in  _Handbook of Statistics_, Vol. 40 (Elsevier, 2019) pp. 551–594.
  * Black and McKane [2012] A. J. Black and A. J. McKane, Stochastic formulation of ecological models and their applications, Trends in ecology & evolution 27, 337 (2012).
  * Majumder  _et al._ [2021] S. Majumder, A. Das, A. Kushal, S. Sankaran, and V. Guttal, Finite-size effects, demographic noise, and ecosystem dynamics, The European Physical Journal Special Topics 230, 3389 (2021).
  * Loreau [2010] M. Loreau, From populations to ecosystems, in  _From Populations to Ecosystems_ (Princeton University Press, 2010).
  * Durrett and Levin [1994] R. Durrett and S. Levin, The importance of being discrete (and spatial), Theoretical population biology 46, 363 (1994).
  * Biancalani  _et al._ [2014] T. Biancalani, L. Dyson, and A. J. McKane, Noise-induced bistable states and their mean switching time in foraging colonies, Physical review letters 112, 038101 (2014).
  * Leyk  _et al._ [2019] S. Leyk, A. E. Gaughan, S. B. Adamo, A. de Sherbinin, D. Balk, S. Freire, A. Rose, F. R. Stevens, B. Blankespoor, C. Frye, _et al._ , The spatial allocation of population: a review of large-scale gridded population data products and their fitness for use, Earth System Science Data 11, 1385 (2019).
  * Nathan  _et al._ [2022] R. Nathan, C. T. Monk, R. Arlinghaus, T. Adam, J. Alós, M. Assaf, H. Baktoft, C. E. Beardsworth, M. G. Bertram, A. I. Bijleveld, _et al._ , Big-data approaches lead to an increased understanding of the ecology of animal movement, Science 375, eabg1780 (2022).
  * Jhawar  _et al._ [2020] J. Jhawar, R. G. Morris, U. Amith-Kumar, M. Danny Raj, T. Rogers, H. Rajendran, and V. Guttal, Noise-induced schooling of fish, Nature Physics 16, 488 (2020).
  * Tunstrøm  _et al._ [2013] K. Tunstrøm, Y. Katz, C. C. Ioannou, C. Huepe, M. J. Lutz, and I. D. Couzin, Collective states, multistability and transitional behavior in schooling fish, PLoS computational biology 9, e1002915 (2013).
  * Stenseth  _et al._ [1997] N. C. Stenseth, W. Falck, O. N. Bjørnstad, and C. J. Krebs, Population regulation in snowshoe hare and canadian lynx: asymmetric food web configurations between hare and lynx, Proceedings of the National Academy of Sciences 94, 5147 (1997).
  * Bjørnstad and Grenfell [2001] O. N. Bjørnstad and B. T. Grenfell, Noisy clockwork: time series analysis of population fluctuations in animals, Science 293, 638 (2001).
  * Lenski [2017] R. E. Lenski, Experimental evolution and the dynamics of adaptation and genome evolution in microbial populations, The ISME journal 11, 2181 (2017).
  * Carpenter  _et al._ [2020] S. R. Carpenter, B. M. Arani, P. C. Hanson, M. Scheffer, E. H. Stanley, and E. Van Nes, Stochastic dynamics of cyanobacteria in long-term high-frequency observations of a eutrophic lake, Limnology and Oceanography Letters 5, 331 (2020).
  * Xie  _et al._ [2008] Y. Xie, Z. Sha, and M. Yu, Remote sensing imagery in vegetation mapping: a review, Journal of plant ecology 1, 9 (2008).
  * Majumder  _et al._ [2019] S. Majumder, K. Tamma, S. Ramaswamy, and V. Guttal, Inferring critical thresholds of ecosystem transitions from spatial data, Ecology 100, e02722 (2019).
  * Horsthemke and Lefever [1989] W. Horsthemke and R. Lefever, Noise-induced transitions, Noise in nonlinear dynamical systems 2, 179 (1989).
  * Gradišek  _et al._ [2000] J. Gradišek, S. Siegert, R. Friedrich, and I. Grabec, Analysis of time series from stochastic processes, Physical Review E 62, 3146 (2000).
  * Friedrich  _et al._ [2011] R. Friedrich, J. Peinke, M. Sahimi, and M. R. R. Tabar, Approaching complexity by stochastic methods: From biological systems to turbulence, Physics Reports 506, 87 (2011).
  * Tabar [2019] R. Tabar,  _Analysis and data-based reconstruction of complex nonlinear dynamical systems_ , Vol. 730 (Springer, 2019).
  * Jhawar and Guttal [2020] J. Jhawar and V. Guttal, Noise-induced effects in collective dynamics and inferring local interactions from data, Philosophical Transactions of the Royal Society B 375, 20190381 (2020).
  * Rinn  _et al._ [2016] P. Rinn, P. G. Lind, M. Wächter, and J. Peinke, The Langevin Approach: An R Package for Modeling Markov Processes, Journal of Open Research Software 4, 10.5334/jors.123 (2016).
  * Dietrich  _et al._ [2021] F. Dietrich, A. Makeev, G. Kevrekidis, N. Evangelou, T. Bertalan, S. Reich, and I. G. Kevrekidis, Learning effective stochastic differential equations from microscopic simulations: combining stochastic numerics and deep learning, arXiv:2106.09004 [physics] (2021), arXiv: 2106.09004.
  * Evangelou  _et al._ [2022] N. Evangelou, F. Dietrich, J. M. Bello-Rivas, A. Yeh, R. Stein, M. A. Bevan, and I. G. Kevekidis, Learning Effective SDEs from Brownian Dynamics Simulations of Colloidal Particles, arXiv:2205.00286 [cs, math] (2022), arXiv: 2205.00286.
  * Brunton  _et al._ [2016] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the national academy of sciences 113, 3932 (2016).
  * Rudy  _et al._ [2017] S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz, Data-driven discovery of partial differential equations, Science advances 3, e1602614 (2017).
  * Boninsegna  _et al._ [2018] L. Boninsegna, F. Nüske, and C. Clementi, Sparse learning of stochastic dynamic equations, The Journal of Chemical Physics 148, 241723 (2018), arXiv: 1712.02432.
  * Callaham  _et al._ [2021] J. L. Callaham, J.-C. Loiseau, G. Rigas, and S. L. Brunton, Nonlinear stochastic modelling with langevin regression, Proceedings of the Royal Society A 477, 20210092 (2021).
  * Huang  _et al._ [2022] Y. Huang, Y. Mabrouk, G. Gompper, and B. Sabass, Sparse inference and active learning of stochastic differential equations from data, arXiv preprint arXiv:2203.11010 (2022).
  * Gardiner [2009] C. Gardiner,  _Stochastic methods_ , Vol. 4 (springer Berlin, 2009).
  * Kolpas  _et al._ [2007] A. Kolpas, J. Moehlis, and I. G. Kevrekidis, Coarse-grained analysis of stochasticity-induced switching between collective motion states, Proceedings of the National Academy of Sciences 104, 5931 (2007).
  * Lehnertz  _et al._ [2018] K. Lehnertz, L. Zabawa, and M. R. R. Tabar, Characterizing abrupt transitions in stochastic dynamics, New Journal of Physics 20, 113043 (2018).
  * pyd [2022] Pydaddy tutorials – advanced function fitting, https://pydaddy.readthedocs.io/en/latest/tutorials/3/advanced_fittiing.html (2022), [Online; accessed 1-May-2022].
  * Drury [2007] K. L. Drury, Shot noise perturbations and mean first passage times between stable states, Theoretical Population Biology 72, 153 (2007).
  * Johst and Wissel [1997] K. Johst and C. Wissel, Extinction risk in a temporally correlated fluctuating environment, Theoretical Population Biology 52, 91 (1997).
  * Shalev-Shwartz and Ben-David [2014] S. Shalev-Shwartz and S. Ben-David,  _Understanding machine learning: From theory to algorithms_ (Cambridge university press, 2014).


11institutetext: International Center of Relativistic Astrophysics Network
(ICRANet), Piazza della Repubblica 10, 65122 Pescara, Italy.
Dipartimento di Fisica, Università degli Studi di Roma “Sapienza”, P.le Aldo
Moro 5, 00185 Roma, Italy.

# Dark Matter Massive Fermions and Einasto Profiles in Galactic Haloes

I. Siutsou  C. R. Argüelles  R. Ruffini<EMAIL_ADDRESS>
(Received / Accepted )

On the basis of a fermionic dark matter model we fit rotation curves of The HI
Nearby Galaxy Survey THINGS sample and compare our 3-parametric model to other
models widely used in the literature: 2-parametric Navarro–Frenk–White,
pseudoisothermal sphere, Burkhert models, and 3-parametric Einasto model,
suggested as the new ”standard dark matter profile” model in the paper by
Chemin et. al., AJ 142 (2011) 109. The results from the fitting procedure
provides evidence for an underlying fermionic nature of the dark matter
candidate, with rest mass above the keV regime.

###### Key Words.:

Dark matter – Astroparticle physics – Galaxies: halos – Galaxies: kinematics
and dynamics

††offprints: I. Siutsou,

##  1 Introduction

The problem of the distribution of Dark Matter (DM) in galaxies, as usually
addressed in the literature, is mainly focused in the halo regions and
associated with the galaxy rotation curves obtained from the observations, see
e.g. Einasto (2013). A well-known approach used to deal with this problem is
the Navarro–Frenk–White (NFW) model (Navarro et al. 1997), expected to provide
a universal description of dark matter halos obtained under the following main
considerations: $1)$ N-body simulations in Cold dark matter (CDM) and
($\Lambda$CDM) cosmologies; $2)$ particles each of masses of $\sim
10^{9}M_{\odot}$ 111Modern numerical simulations can reach better resolution
down to particle masses of $\sim 10^{5}M_{\odot}$ (Gao et al. 2012).; $3)$
classical Newtonian physics.

Despite an indicated agreement of this model with the large scale structure of
the Universe, some problems remains at galactic scales, see e.g. Munshi et al.
(2013). A central characteristic of the NFW dark matter density profiles, is
that they show a cuspy and divergent behaviour through the center of the
configuration, while empirical profiles tend to show a core of constant
density, giving rise to the well-known core-cusp controversy, see e.g. de Blok
(2010).

Yet another important approach developed to understand the distribution of
matter in galaxies has been advanced by Einasto (1965) and Einasto & Haud
(1989). This is a phenomenological approach consisting in the proposal of an
empirical fitting function composed by three free parameters as detailed in
equation (1)

$\rho_{E}(r)=\rho_{-2}\exp{\left(-\frac{2}{n}\left[\left(\frac{r}{r_{-2}}\right)^{n}-1\right]\right)},$
(1)

where $\rho_{-2}$ and $r_{-2}$ are the density and radius at which
$\rho(r)\propto r^{-2}$, and $n$ is the Einasto index which determines the
shape of the profile.

Recent N-body simulations in $\Lambda$CDM cosmology by Navarro et al. (2004)
purported a novel dark matter halo model different from NFW. This model was
soon realized (Merritt et al. 2006) to be the same as the Einasto one as given
by equation (1).

After that, using the highest quality rotation curves available to date
obtained from The HI Nearby Galaxy Survey (THINGS) (Walter et al. 2008; de
Blok et al. 2008), the Einasto dark matter halo model has been proposed as the
standard model for dark matter halos by Chemin et al. (2011), as it provides
both cored and cusped distributions for different values of model parameters
(see Fig. 1). In that work, the fundamental core-cusp discrepancy is analyzed
in detail for the whole sample of galaxies under study. It is clearly shown
that for the majority of the galaxies considered in the sample, the cored
halos (compatible with near unity Einasto indexes) are preferred over the
cuspy ones (these instead compatible with higher Einasto indexes).

Figure 1: Comparison of the density profiles for different phenomenological
models of dark matter distribution.

We present here a novel approach focusing on galactic structures and an
underlying microphysical component of Dark Matter. The model is built upon the
following general considerations: $1)$ the Dark Matter component is assumed to
be chargeless spin-1/2 fermions; $2)$ the configurations are described by
General Relativity; $3)$ the particles are assumed to be isothermal in
thermodynamic equilibrium (i.e. without the need of pre-fixing any
cosmological history). The theoretical fundament of this new approach is
detailed in the model of semi-degenerate self-gravitating fermions first
introduced by Gao et al. (1990), and more recently with applications to
galactic dark matter by Ruffini et al. (2013); Argüelles et al. (2013). Our
model is based in the following main assumptions:

  1. 1.

the problem of galactic cores and halos have to be addressed unitarily;

  2. 2.

for definiteness we study the simplest problem of ”bare” massive particles,
neglecting at this stage all other interactions than the gravitational one and
fulfilling only the Fermi–Dirac statistical distribution

$f=\frac{1}{\exp\left({\frac{\epsilon-\mu}{kT}}\right)+1}=\frac{1}{\exp{\left(\frac{\epsilon}{\beta
mc^{2}}-\theta\right)}+1},$ (2)

where $\epsilon$ is kinetic energy of the particles, $\mu$ is chemical
potential, $T$ is the temperature, $k$ is Boltzmann constant and $c$ is the
speed of light. The mass of the particle ($m$), the temperature parameter
($\beta=kT/mc^{2}$) and the degeneracy parameter ($\theta=\mu/kT$) at the
center are the three free parameters of the model;

  3. 3.

we consider zero total angular momentum and also neglect any effect of
baryonic matter on the DM in the mathematical formulation.

It is shown that in any such system the density at large radii scales as
$r^{-2}$ independently of the values of the central density, providing the
flat rotation curve (Gao et al. 1990; Ruffini et al. 2013; Argüelles et al.
2013).

The dark matter halos obtained in the new dark matter approach proposed here
share a common or universal feature which shed more light on the core-cusp
discrepancy, while providing a new mass scale to the dark matter candidate.
Our density profiles always favor a cored behaviour (without any cusp) in the
observed inner halo regions, given the quantum nature of the fermionic
particles. Another fundamental outcome of our model is the range of the DM
particle mass, which must be $m\gtrsim 5$ keV in order to be in agreement with
typical halo sizes of the observed dwarf galaxies (Ruffini et al. 2013).

The theoretical formulation of Ruffini et al. (2013); Argüelles et al. (2013)
is based on the first principles physics and provides a physical complement to
Einasto phenomenological models. It also offers the necessity to approach the
Dark Matter distribution in galactic haloes with fermions with masses larger
than the above mentioned bound.

Figure 2: Semidegenerate density profile in dimensionless units for central
degeneracy parameter $\theta_{0}=15$ and central temperature parameter
$\beta_{0}=10^{-10}$.

The paper is structured as follows. We model the distribution of Dark Matter
as semidegenerate fully relaxed thermal self-gravitating general relativistic
fermionic solutions of Gao et al. (1990), see Sec. 2. The resulting density
profiles provide flat rotation curve at large distances, cored distribution of
dark matter in the halo, and a massive degenerate core at the very center, see
Sec. 3. We describe in Sec. 4 the actual procedure of fitting of rotation
curves, and then discuss the results in Sec. 5. Conclusions follow.

##  2 Model equations

We consider self-gravitating system of fermions in thermal equilibrium
following Gao et al. (1990) with occupation numbers given by

$f(\epsilon)=\frac{1}{e^{\frac{\epsilon-\mu}{kT}}+1}.$ (3)

Then equation of state reads

$\displaystyle\rho=m\frac{g}{h^{3}}\int\frac{1+\epsilon/mc^{2}}{e^{\frac{\epsilon-\mu}{kT}}+1}\,d^{3}p,$
(4) $\displaystyle
P=\frac{2}{3}\frac{g}{h^{3}}\int\frac{(1+\epsilon/mc^{2})^{-1}(1+\epsilon/2mc^{2})\epsilon}{e^{\frac{\epsilon-\mu}{kT}}+1}\,d^{3}p,$
(5)

where $g=2s+1$, $s$ is the spin of the particle, and integration is extended
over all 3-momentum space.

The Einstein equations for the spherically symmetric metric

$g_{\mu\nu}={\rm diag}(e^{\nu},-e^{\lambda},-r^{2},-r^{2}\sin^{2}\theta),$ (6)

where $\nu$ and $\lambda$ depend only on the radial coordinate $r$, together
with the thermodynamic equilibrium conditions of Tolman (1930) and Klein
(1949)

$e^{\nu/2}T=const\,,\qquad e^{\nu/2}(\mu+mc^{2})=const,$ (7)

can be written in the dimensionless form of Gao et al. (1990)

$\displaystyle\frac{d\hat{M}}{d\hat{r}}$
$\displaystyle=4\pi\hat{r}^{2}\hat{\rho},$ (8)
$\displaystyle\frac{d\theta}{d\hat{r}}$
$\displaystyle=-\frac{1-\beta_{0}(\theta-\theta_{0})}{\beta_{0}}\frac{\hat{M}+4\pi\hat{P}\hat{r}^{3}}{\hat{r}^{2}(1-2\hat{M}/\hat{r})},$
(9) $\displaystyle\frac{d\nu}{d\hat{r}}$
$\displaystyle=\frac{\hat{M}+4\pi\hat{P}\hat{r}^{3}}{\hat{r}^{2}(1-2\hat{M}/\hat{r})},$
(10) $\displaystyle\beta_{0}$
$\displaystyle=\beta(r)e^{\frac{\nu(r)-\nu_{0}}{2}},$ (11) $\displaystyle
e^{-\lambda}$ $\displaystyle=1-\frac{2\hat{M}(\hat{r})}{\hat{r}}\,.$ (12)

The following dimensionless quantities were introduced:

$\displaystyle\hat{r}=r/\chi,$ (13) $\displaystyle\hat{M}=GM/(c^{2}\chi),$
(14) $\displaystyle\hat{\rho}=G\chi^{2}\rho/c^{2},$ (15)
$\displaystyle\hat{P}=G\chi^{2}P/c^{4},$ (16)

where $\chi=2\pi^{3/2}(\hbar/mc)(m_{p}/m)$ is the characteristic length that
scales as $m^{-2}$, with $m_{p}=\sqrt{\hbar c/G}$ being the Planck mass, and
the temperature and degeneracy parameters, $\beta=kT/(mc^{2})$ and
$\theta=\mu/(kT)$, respectively. The constants of the equilibrium conditions
of Tolman and Klein have been evaluated at the center $r=0$, which we indicate
with a subscript ‘0’.

The system of coupled differential equations (8–12) is solved for initial
conditions $M(0)=\nu(0)=0$ and given set of free parameters $\beta_{0}$ and
$\theta_{0}$, $m$ for each galaxy under study as detailed below.

##  3 Properties of semidegenerate configurations

Galactic halos have to be necessarily composed from cold particles, so that
astrophysically relevant solutions will have temperature parameters $\beta\ll
1$. In this case, the general solution for semidegenerate configurations
($\theta_{0}\gtrsim 10$) present three different regions: an inner degenerate
compact core, an extended low-degenerate inner halo of almost constant density
and a non-degenerate outer halo with characteristic slope $\rho\propto r^{-2}$
(see Gao et al. (1990) and Fig. 2 for additional details). The infinite mass
of the configuration extended up to spatial infinity is not a problem, because
in reality it is limited by tidal interactions with other galaxies, which
introduce an energy cutoff into the distribution function, see e.g. Ingrosso
et al. (1992). However this is not important for the inner parts of the
configuration we are interested in.

In order to understand the crucial properties of this equilibrium
configurations we plot the circular velocity of a test body in the metric
fulfilling Eqs. (8–12) on Fig. 3. There are indeed four regions of the
solution for circular velocity, each with characteristic slope. The inner
region I correspond to the degenerate core of almost constant density, so that
$v_{\text{circ}}\propto r$. For increasing values of radial coordinate the
inner halo follows, itself composed of two different regions. In the region II
the density of dark matter sharply decreases and this Keplerian region is
dominated by the mass of the degenerate core, and as a result
$v_{\text{circ}}\propto r^{-1/2}$. For yet increasing values of the radial
coordinate the density of dark matter reach an almost constant value giving
rise to a plateau, see Fig. 2. As soon as the mass of the plateau prevails
over the mass of the core we have the region III where $v_{\text{circ}}\propto
r$. Finally in the region IV, after some oscillations the circular velocity
tends to a constant independent on $r$, corresponding to a pure Boltzmannian
regime and characteristic for the flat rotation curve of outer halo.

Figure 3: Dependence of $\hat{v}_{\text{circ}}=v_{\text{circ}}/c$ on
dimensionless radius $\hat{r}$, $\beta=10^{-8},\ \theta_{0}=50$.

We define the physical characteristics of each configuration as follows:

  * •

The characteristic radius of the core $\hat{r}_{c}$ is given by
$\hat{v}_{\text{circ}}(\hat{r}_{c})=\mathrm{max}$ in region I.

  * •

$\hat{M}_{c}$ is the mass of the core given by
$\hat{M}_{c}=\hat{v}^{2}_{\text{circ}}\hat{r}$ in the region II.

  * •

The characteristic radius of the inner halo $\hat{r}_{h}$ correspond to
$\hat{v}_{\text{circ}}(\hat{r}_{h})=\mathrm{max}$ in region III.

  * •

The characteristic mass of the inner halo $\hat{M}_{h}$ is given by
$\hat{M}_{h}=\hat{M}(\hat{r}_{h})$ just between the regions III and IV.

For the parameters in the region of
$\theta_{0}\in[0,200],\log\beta_{0}\in[-10,-5]$ we calculate a grid of models
and extracted numerically the physical characteristics mentioned above. Then
we fit the obtained values by different double parametric functions and find
out the best fitting formulae with the correspondent
($\beta_{0}$,$\theta_{0}$) dependence for the range of
$\theta_{0}\in[20,200]$. An interesting fact is that in our region of
parameters circular velocity $v_{\text{circ}}$ in the flat part of region IV
(i.e. $v_{\infty}$) is defined by temperature $\beta_{0}$ only. In the range
of astrophysically relevant parameters
$\theta_{0}\in[10,200],\log\beta_{0}\in[-10,-5]$ the scaling relation between
circular velocity and $\beta_{0}$ corresponds to the Boltzmannian relation
between $v_{\text{circ}}$ and one-dimensional dispersion velocity
$\sigma=\sqrt{kT/m}$ (see e.g. Binney & Tremaine (1987))

$\frac{v_{\infty}}{\text{km/s}}=\sqrt{2}c\sqrt{\beta_{0}}.$ (17)

We have here neglected general relativistic corrections which are very small
in these ranges of parameters, i.e. $e^{\nu(\infty)-\nu(0)}\approx 1$, and
then $\beta(\infty)\approx\beta_{0}$ by equation (11).

For the temperature and degeneracy free parameters in the range
$\log\beta_{0}\in[-10,-5]$, $\theta_{0}\in[20,200]$, respectively, we obtain
the following dimensionless scaling laws for core radius and mass222For radii
$r<r_{c}$ the configurations corresponds to region I, where $\rho(r)\approx
const$ and then from (18) and (19), taking into account (13) and (14), we can
write $\rho_{c}\propto M_{c}/r_{c}^{3}$ in terms of the chemical potential and
particle mass ($\beta_{0}\theta_{0}=\mu_{0}/mc^{2}$) as
$\rho_{c}\propto\mu_{0}^{3/2}m^{5/2}$. This dependence is precisely the one of
a fully degenerate non-relativistic Fermi gas in presence of an external
gravitational field, which is further coinciding with a polytrope of index
$n=3/2$ (see e.g. Shapiro & Teukolsky (1983)).

$\displaystyle\hat{r}_{c}=0.226(\beta_{0}\theta_{0})^{-1/4},$ (18)
$\displaystyle\hat{M}_{c}=0.234(\beta_{0}\theta_{0})^{3/4}.$ (19)

However, the core region is typically very small and is not constrained by
empirical data of the THINGS sample considered here. Moreover, the mass
contribution of regions I and II to the total mass $M_{h}$ at the end of
region III is $\lesssim 10^{-2}$ in our parameter range as shown in Ruffini et
al. (2013), indicating that only regions III and IV are the relevant ones to
be used in the fitting procedure against the data.

At this point it is important to emphasize that the theoretical treatment used
here to fit dark matter halos applies for any core size, even for the ones
which are close to its critical mass, $M_{cr}\sim 10^{9}M_{\odot}$, as studied
in Argüelles et al. (2014), where a relativistic treatment is mandatory. Even
though the regions I-II-III-IV corresponding to the sample considered here can
be well explained in terms of non-relativistic physics, the general
relativistic approach has been used for formal correctness, only giving
negligible corrections.

Dimensionless halo radius and mass have different scalings, they are
proportional not to $\theta_{0}^{\alpha}$, but to $\alpha^{\theta_{0}}$

$\displaystyle\hat{r}_{h}=0.953\beta_{0}^{1/4}(1.445)^{\theta_{0}},$ (20)
$\displaystyle\hat{M}_{h}=2.454\beta_{0}^{3/4}(1.445)^{\theta_{0}}.$ (21)

Formulas (20–21) represent perfect scalings in the region of parameters
considered above, which involves also the scaling of the whole rotational
curves in regions III and IV; formula (17) shows a perfect scaling in the flat
part of the rotation curve for region IV. Moreover, the Newtonian expression
for the dimensionless circular velocity
$\hat{v}_{circ}^{2}(\hat{r})=\hat{M}/\hat{r}$ is perfectly suitable in the
physical region under consideration. The expression for maximal rotation
velocity in the halo is thus obtained from (20) and (21) and reads

$\hat{v}^{2}_{h}(\hat{r}_{h}(\beta_{0},\theta_{0}))=2.575\beta_{0}^{1/2}.$
(22)

In next section we explain the fitting procedure with the use of the halo
scaling laws for regions III and IV obtained here.

##  4 Observed rotation curves and fitting procedure

In 2008, a sample of 34 nearby (closer than 15 Mpc) spiral and irregular
galaxies (Sb to Im) were observed with The HI Nearby Galaxy Survey (THINGS)
(Walter et al. 2008). These observations allowed to obtain the highest quality
rotation curves available to date due to the high spatial and velocity
resolution of THINGS. Then a sub-sample of these rotation curves,
corresponding to 19 rotationally dominated and undisturbed galaxies, were
combined with information on the distribution of gas and stars by de Blok et
al. (2008) to construct mass models for the dark matter component of the
sample. These models finally were used to quantify the dark matter
contribution for each galaxy by using the following formula

$V^{2}_{obs}=V_{gas}^{2}+\Upsilon_{*}V_{*}^{2}+V_{DM}^{2},$ (23)

which relates the observed input curves of $V_{obs}$, $V_{gas}$ and $V_{*}$,
defined below, with the dark matter rotation curve $V_{DM}$ to be determined
from the known input data once the the mass-to-light ratio $\Upsilon_{*}$ is
provided.

The total and gas observed rotation curves $V_{obs}$ and $V_{gas}$,
respectively, were both obtained from the THINGS data: the first was obtained
from velocity fields analysis and the second from the neutral hydrogen (HI)
distribution maps, as described in de Blok et al. (2008). Instead, each
stellar (light) rotation curve $V_{*}$ is obtained from the corresponding
stellar distribution observed in the K band (i.e. at $3.6\ \mu$m) by the
Spitzer Infrared Nearby Galaxy Survey (SINGS), independent of THINGS, and
described in de Blok et al. (2008) and references therein. Finally, the mass-
to-light ratio $\Upsilon_{*}^{K}$ was used to determine the rotation curve
associated with the stellar mass distribution from that of the measured light.

At this point it is relevant to further emphasize the underlying hypothesis to
which equation (23) is subject to. This is, each baryonic rotation velocity
$V_{gas}$ and $V_{*}$ was calculated from the correspondent observed baryonic
mass density distribution, and was defined as the velocity that each component
would induce on a test particle in the galactic plane as if they were isolated
of any external influence.

In de Blok et al. (2008) equation (23) was applied to test the cuspy Navarro-
Frenk-White and cored pseudo-ISO dark matter models against data as follows:
the (squared) rotation curves of the baryonic components (after appropriate
scaling with $\Upsilon_{*}^{K}$) were subtracted from the (squared) observed
rotation curve $V^{2}_{obs}$ to apply a reduced $\chi^{2}$ fitting procedure
in order to find the best fitting free parameters for each dark matter model.
Soon after, the same analysis was extended further to Einasto dark matter
profiles by Chemin et al. (2011), concluding that the Einasto model provides
the best match to the observed rotation curves when compared with NFW and
pseudo-ISO models with empirical fixed values for $\Upsilon_{*}^{K}$ for two
different stellar initial mass functions (IMFs).

Here we propose a different dark matter halo model, which is neither based on
numerical N-body simulations nor on phenomenological model proposals, but
relies on the underlying microphysical composition of the dark matter
candidate, as explained in former sections.

Thus, analogously to de Blok et al. (2008) and Chemin et al. (2011) we use the
HI high resolution observations of galaxies from THINGS survey (Walter et al.
2008). We analyze here the sample of 16 rotationally dominated and undisturbed
galaxies presented both in de Blok et al. (2008) and Chemin et al. (2011),
listed below in Table 1.

Regarding the rotation curves data of the baryonic components, we consider the
contributions of the gas, the stellar disk and a spherical stellar bulge
$V_{b}$ as given in Chemin et al. (2011). The halo rotation velocity
corresponding to the spherical dark matter component is taken from the two
parametric scaling, see (20–21), which we name from now on as the fermionic
dark matter velocity profile $V_{f}(r)$.

Once each component is provided, we make use of the equation analogous to (23)

$V^{2}_{obs}=V_{gas}^{2}+\Upsilon_{*}V_{*}^{2}+V_{b}^{2}+V_{f}^{2},$ (24)

With all the baryonic velocity terms ($V_{gas}^{2}$, $\Upsilon_{*}V_{*}^{2}$
and $V_{b}^{2}$) as observational inputs, we fit the HI observed rotation
curve $V^{2}_{obs}$ by Levenberg–Marquardt nonlinear least-squares algorithm,
in complete analogy as done in Chemin et al. (2011).

We did not take into account the contribution of molecular gas because the
total gas surface densities are dominated by atomic gas for the majority of
the sample, as explained in Chemin et al. (2011) and references therein. Total
rotation curve was taken from de Blok et al. (2008). We have not considered
models with free mass-to-light ratios, deferring it to a future paper (in
preparation). Instead following Chemin et al. (2011) we have adopted the fixed
mass-to-light ratios of stellar populations with a bursty star formation
history with a Kroupa IMF. We choose this IMF instead of the diet-Salpeter
IMF, also considered in Chemin et al. (2011) and de Blok et al. (2008), as it
generally provides better agreement with observations for rotation curves (see
Fig. 5 of Chemin et al. (2011)), and in some cases the Salpeter IMF leads to
rotational velocities due to stellar component only already in pronounced
excess over observed total rotational velocity (see, for example, cases of
NGC3521 and NGC5055 at Fig. 3 of Chemin et al. (2011)).

Together with the dark matter profile of semidegenerate configurations with
particle mass $m=10$ keV/$c^{2}$ and varying $\theta_{0}$ and $\beta_{0}$, and
for the sake of comparison, the following profiles were also used for fitting:

  * •

Cored profiles with central density $\rho_{0}$ and characteristic radius
$r_{0}$:

    * –

pseudo-isothermal sphere profile

$\rho_{DM}(r)=\rho_{0}\frac{r_{0}^{2}}{r^{2}+r_{0}^{2}},$ (25)

    * –

Burkert profile

$\rho_{DM}(r)=\rho_{0}\frac{r_{0}^{3}}{(r_{0}+r)(r_{0}^{2}+r^{2})}.$ (26)

  * •

Cusped profiles with characteristic radius $r_{-2}$ where the density profile
has a (logarithmic) slope of $-2$ (the ”isothermal” value) and $\rho_{-2}$ as
the local density at that radius. In the case of Einasto profiles a third
parameter is needed, the Einasto index $n$ which determines the shape of the
profile.

    * –

Navarro–Frenk–White profile

$\rho_{DM}(r)=4\rho_{-2}\frac{r_{-2}}{r}\left(\frac{r_{-2}}{r+r_{-2}}\right)^{2},$
(27)

    * –

Einasto profile333The family of Einasto profiles with relatively large indices
$n>4$ are identified with cuspy halos, while low index values $n<4$ presents a
cored-like behaviour (Chemin et al. 2011). The lower the $n$ the more cored-
like the halo profile.

$\rho_{DM}(r)=\rho_{-2}\exp\left\\{-2n\left[\left(\frac{r}{r_{-2}}\right)^{1/n}-1\right]\right\\}.$
(28)

##  5 Results and discussion

In this section, we compare the fits of rotation curves by the different
models considered in the last section. As we have to compare models with
different number of parameters, which are not nested into each other, we use
Bayesian Information Criterion (BIC) introduced by Schwarz (1978). It provides
a penalty to models with larger number of parameters to check what of them is
more likely to be correct. Model with minimum BIC value is preferred. For the
models with the same number of parameters, BIC is equivalent to $\chi^{2}$
criterion.

Table 1: Results of fitting

| Semidegenerate | Burkert
---|---|---
Galaxy | $\beta,10^{-8}$ | $\theta_{0}^{*}$ | $\chi^{2}_{r}$ | BIC |  $r_{0}$, kpc | $\rho_{0},\ 10^{-3}\ M_{\odot}/\text{pc}^{3}$ | $\chi^{2}_{r}$ | BIC
NGC2366 |  0.99${}\pm{}$ | 0.02 |  24.22${}\pm{}$ | 0.08 | 0.10 | _27_ |  2.2${}\pm{}$ | 0.2 |  43${}\pm{}$ | 10 | 0.12 | 35
NGC2403 |  7.10${}\pm{}$ | 0.06 |  27.43${}\pm{}$ | 0.07 | 2.7 | 902 |  4.08${}\pm{}$ | 0.06 |  83${}\pm{}$ | 3 | 2.3 | 866
NGC2841 |  10.8${}\pm{}$ | 0.4 |  32.04${}\pm{}$ | 0.17 | 2.6 | _366_ |  20.6${}\pm{}$ | 0.9 |  5.2${}\pm{}$ | 0.5 | 2.7 | 370
NGC2903 |  16.33${}\pm{}$ | 0.09 |  26.88${}\pm{}$ | 0.06 | 0.66 | _238_ |  2.89${}\pm{}$ | 0.06 |  388${}\pm{}$ | 18 | 1.1 | 283
NGC2976 |  9${}\pm{}$ | 3 |  28.8${}\pm{}$ | 0.4 | 0.44 | _93_ |  20${}\pm{}$ | 20 |  40${}\pm{}$ | 150 | 0.49 | 97
NGC3031 |  9.7${}\pm{}$ | 0.3 |  26.6${}\pm{}$ | 0.2 | 3.9 | 470 |  2.63${}\pm{}$ | 0.10 |  270${}\pm{}$ | 20 | 3.9 | 468
NGC3198 |  7.44${}\pm{}$ | 0.07 |  28.66${}\pm{}$ | 0.08 | 1.06 | 254 |  6.32${}\pm{}$ | 0.18 |  36${}\pm{}$ | 2 | 0.99 | _248_
IC2574 |  2.00${}\pm{}$ | 0.06 |  27.96${}\pm{}$ | 0.09 | 0.28 | 167 |  8.0${}\pm{}$ | 0.6 |  7.2${}\pm{}$ | 1.2 | 0.10 | 67
NGC3521 |  8.9${}\pm{}$ | 0.7 |  28.3${}\pm{}$ | 0.4 | 4.1 | _437_ |  5.4${}\pm{}$ | 0.4 |  60${}\pm{}$ | 8 | 4.2 | 437
NGC3621 |  7.35${}\pm{}$ | 0.11 |  28.70${}\pm{}$ | 0.08 | 3.0 | 496 |  6.48${}\pm{}$ | 0.12 |  34.3${}\pm{}$ | 1.3 | 2.5 | 475
NGC4736 |  3.01${}\pm{}$ | 0.13 |  22.4${}\pm{}$ | 0.5 | 2.1 | 302 |  0.84${}\pm{}$ | 0.07 |  870${}\pm{}$ | 160 | 1.9 | 296
DDO154 |  0.791${}\pm{}$ | 0.017 |  24.37${}\pm{}$ | 0.10 | 0.84 | 172 |  2.32${}\pm{}$ | 0.10 |  29${}\pm{}$ | 3 | 0.62 | 153
NGC5055 |  8.35${}\pm{}$ | 0.13 |  30.99${}\pm{}$ | 0.12 | 1.9 | 708 |  14.3${}\pm{}$ | 0.5 |  8.0${}\pm{}$ | 0.5 | 2.0 | 719
NGC6946 |  8.9${}\pm{}$ | 0.3 |  26.9${}\pm{}$ | 0.2 | 16.4 | 993 |  3.48${}\pm{}$ | 0.08 |  146${}\pm{}$ | 7 | 15.7 | 985
NGC7331 |  11.7${}\pm{}$ | 0.2 |  30.10${}\pm{}$ | 0.12 | 0.45 | 226 |  9.7${}\pm{}$ | 0.7 |  25${}\pm{}$ | 4 | 0.41 | 214
NGC7793 |  4.44${}\pm{}$ | 0.16 |  25.85${}\pm{}$ | 0.13 | 3.7 | 287 |  2.60${}\pm{}$ | 0.07 |  130${}\pm{}$ | 8 | 3.5 | 284

| Navarro–Frenk–White | Pseudo-ISO
---|---|---
Galaxy |  $r_{-2}$, kpc | $\rho_{-2},\ 10^{-3}\ M_{\odot}/\text{pc}^{3}$ | $\chi^{2}_{r}$ | BIC |  $r_{0}$, kpc | $\rho_{0},\ 10^{-3}\ M_{\odot}/\text{pc}^{3}$ | $\chi^{2}_{r}$ | BIC
NGC2366 |  200${}\pm{}$ | 1100 |  0.02${}\pm{}$ | 0.24 | 1.1 | 117 |  1.29${}\pm{}$ | 0.17 |  40${}\pm{}$ | 11 | 0.15 | 42
NGC2403 |  11.2${}\pm{}$ | 0.3 |  2.86${}\pm{}$ | 0.17 | 0.7 | 575 |  1.56${}\pm{}$ | 0.04 |  144${}\pm{}$ | 8 | 1.2 | 703
NGC2841 |  150${}\pm{}$ | 30 |  0.05${}\pm{}$ | 0.02 | 3.8 | 402 |  12.5${}\pm{}$ | 0.7 |  4.6${}\pm{}$ | 0.6 | 2.8 | 374
NGC2903 |  4.75${}\pm{}$ | 0.16 |  34${}\pm{}$ | 2 | 1.8 | 334 |  0.53${}\pm{}$ | 0.05 |  2300${}\pm{}$ | 500 | 3.9 | 406
NGC2976 |  900${}\pm{}$ | 40000 |  0.009${}\pm{}$ | 1 | 2.1 | 158 |  9${}\pm{}$ | 18 |  30${}\pm{}$ | 180 | 0.49 | 98
NGC3031 |  4.9${}\pm{}$ | 0.4 |  20${}\pm{}$ | 3 | 4.0 | 472 |  0.82${}\pm{}$ | 0.1 |  690${}\pm{}$ | 170 | 4.3 | 480
NGC3198 |  16.5${}\pm{}$ | 0.9 |  1.37${}\pm{}$ | 0.15 | 2.0 | 306 |  2.7${}\pm{}$ | 0.14 |  51${}\pm{}$ | 5 | 1.2 | 266
IC2574 |  500${}\pm{}$ | 1300 |  0.007${}\pm{}$ | 0.045 | 1.5 | 331 |  5.1${}\pm{}$ | 0.4 |  6.3${}\pm{}$ | 1.1 | 0.11 | 70
NGC3521 |  18${}\pm{}$ | 3 |  1.5${}\pm{}$ | 0.4 | 5.1 | 457 |  2.4${}\pm{}$ | 0.3 |  78${}\pm{}$ | 19 | 4.2 | 439
NGC3621 |  123${}\pm{}$ | 17 |  0.08${}\pm{}$ | 0.02 | 5.9 | 579 |  2.81${}\pm{}$ | 0.09 |  49${}\pm{}$ | 3 | 1.1 | 377
NGC4736 |  1.25${}\pm{}$ | 0.16 |  90${}\pm{}$ | 20 | 1.9 | 296 |  0${}\pm{}$ | 0.07 |  0${}\pm{}$ | ND | 2.3 | 311
DDO154 |  14${}\pm{}$ | 2 |  0.35${}\pm{}$ | 0.11 | 1.03 | 184 |  1.22${}\pm{}$ | 0.07 |  32${}\pm{}$ | 4 | 0.48 | 138
NGC5055 |  48${}\pm{}$ | 4 |  0.20${}\pm{}$ | 0.03 | 3.0 | 795 |  7.8${}\pm{}$ | 0.4 |  7.7${}\pm{}$ | 0.8 | 2.5 | 759
NGC6946 |  9.3${}\pm{}$ | 0.4 |  5.3${}\pm{}$ | 0.5 | 10.4 | 915 |  0.66${}\pm{}$ | 0.03 |  870${}\pm{}$ | 70 | 11.0 | 923
NGC7331 |  3000${}\pm{}$ | 40000 |  0.004${}\pm{}$ | 0.087 | 0.48 | 233 |  5.${}\pm{}$ | 0.5 |  28${}\pm{}$ | 6 | 0.32 | 190
NGC7793 |  17.0${}\pm{}$ | 1.9 |  1.4${}\pm{}$ | 0.3 | 4.1 | 294 |  1.47${}\pm{}$ | 0.05 |  126${}\pm{}$ | 10 | 4.0 | 293

| Einasto
---|---
Galaxy |  $r_{-2}$, kpc | $\rho_{-2},\ 10^{-3}\ M_{\odot}/\text{pc}^{3}$ | $n$ | $\chi^{2}_{r}$ | BIC
NGC2366 |  2.9${}\pm{}$ | 0.4 |  6.86${}\pm{}$ | 0.04 |  0.9${}\pm{}$ | 0.3 | 0.13 | 39
NGC2403 |  13.6${}\pm{}$ | 1.3 |  1.98${}\pm{}$ | 0.05 |  4.9${}\pm{}$ | 0.4 | 0.6 | _570_
NGC2841 |  24.5${}\pm{}$ | 0.6 |  1.091${}\pm{}$ | 0.006 |  0.54${}\pm{}$ | 0.08 | 2.5 | 367
NGC2903 |  5.33${}\pm{}$ | 0.15 |  28.983${}\pm{}$ | 0.005 |  2.9${}\pm{}$ | 0.2 | 1.6 | 326
NGC2976 |  70000${}\pm{}$ | ND |  0.019${}\pm{}$ | ND |  4.0${}\pm{}$ | 70 | 0.49 | 100
NGC3031 |  4.81${}\pm{}$ | 0.09 |  30.159${}\pm{}$ | 0.002 |  0.56${}\pm{}$ | 0.07 | 3.2 | _452_
NGC3198 |  11.5${}\pm{}$ | 0.4 |  3.029${}\pm{}$ | 0.006 |  1.80${}\pm{}$ | 0.17 | 1.1 | 257
IC2574 |  8.6${}\pm{}$ | 1.2 |  1.57${}\pm{}$ | 0.05 |  0.7${}\pm{}$ | 0.2 | 0.09 | _62_
NGC3521 |  9.3${}\pm{}$ | 0.5 |  6.27${}\pm{}$ | 0.02 |  1.2${}\pm{}$ | 0.3 | 4.3 | 443
NGC3621 |  37${}\pm{}$ | 10 |  0.3${}\pm{}$ | 0.4 |  6.1${}\pm{}$ | 0.9 | 0.57 | _296_
NGC4736 |  1.73${}\pm{}$ | 0.16 |  56.78${}\pm{}$ | 0.03 |  2.0${}\pm{}$ | 0.5 | 1.8 | _295_
DDO154 |  4.9${}\pm{}$ | 0.7 |  1.95${}\pm{}$ | 0.03 |  2.0${}\pm{}$ | 0.3 | 0.31 | _114_
NGC5055 |  21.4${}\pm{}$ | 0.4 |  1.464${}\pm{}$ | 0.002 |  0.39${}\pm{}$ | 0.04 | 1.2 | _626_
NGC6946 |  50${}\pm{}$ | 40 |  0.2${}\pm{}$ | 4 |  15${}\pm{}$ | 4 | 8.5 | _883_
NGC7331 |  1000${}\pm{}$ | 8000 |  0.002${}\pm{}$ | 900 |  11${}\pm{}$ | 21 | 0.23 | _158_
NGC7793 |  3.65${}\pm{}$ | 0.14 |  19.217${}\pm{}$ | 0.007 |  0.97${}\pm{}$ | 0.09 | 3.1 | _279_

ND means not constrained model parameters

The results of fitting are presented in the Table 1 and on Fig. 4. Due to the
scaling laws recalled in Sec. 3 the fits obtained here for particle mass
$m=10$ keV/$c^{2}$ can be also transferred to another particle mass by
changing the fitted central degeneracy parameter from $\theta_{0}^{*}$ (see
Table 1) according to the relation

$\theta_{0}(m)=\theta_{0}^{*}+12.52\log\frac{m}{10\text{ keV}/c^{2}},$ (29)

provided that $\theta_{0}(m)$ is larger than 20 and the influence of the
degenerate core on rotational velocity is negligible in the observed radial
range. From the values of $\theta_{0}^{*}$ obtained for the fitting of the
sample listed in Table 1, and the lower value of $\theta_{0}(m)$ from which
the scaling laws (18–21) are valid, it is possible to obtain from (29) a
preliminary lower limit for the particle mass $m\gtrsim$ few keV/$c^{2}$.
Nonetheless this limit should not be considered as an absolute lower limit for
the fermion mass of the model because the bound $\theta_{0}(m)\geqslant 20$ in
the formula above is a numerical limit, and no underlying physics has been
specified here for it. The formal way of providing an absolute lower limit for
the particle mass of our model when applied to typical spiral galaxies has
been found in Ruffini et al. (2013), and yields roughly an order of magnitude
less than the one inferred here.

From the 16 galaxies analyzed, our model has minimum BIC value in 5 cases
(NGC2366, NGC2841, NGC2903, NGC2976, NGC3521), Einasto model in 10 cases
(NGC2403, NGC3031, IC2574, NGC3621, NGC4736, DDO154, NGC5055, NGC6946,
NGC7331, NGC7793), and in the case of NGC3198 Burkert model is the best one,
marginally better than ours, which in turn is marginally better than Einasto.
Besides this general comparison in which apparently Einasto model is preferred
against our model, there is a more relevant comparison which must be made
considering that the semi-degenerate model provides cored halos only. For
this, we compare the Einasto model against the semi-degenerate one for the
sub-set of galaxies which are cored-like (i.e. with Einasto index $n\lesssim
4$), and then the same comparison is made for the sub-set of galaxies which
are cuspy-like (i.e. with Einasto index $n>4$). The important outcome of this
new BIC comparison is that our model is equivalently as good as Einasto for
the cored-like sub-sample, that is connected to the fact that Einasto profile
with $n\sim 1$ provides rotational curves that in a wide range of radii is
quite close to the one of ours at transition from region III to region IV. For
cored-like galaxies in 6 cases our model has lower BIC number than Einasto
model (NGC2366, NGC2841, NGC2903, NGC2976, NGC3198, NGC3521), and inversely in
other 6 cases Einasto is better (NGC3031, IC2574, NGC4736, DDO154, NGC5055,
NGC7793). Instead, for the cuspy-like sub-sample (NGC2403, NGC3621, NGC6946,
NGC7331), in all the cases Einasto model has lower BIC numbers as logically
one may expect due to the cored nature of the semi-degenerate halos.

If we take only the models with significant fits, i.e. with reduced $\chi^{2}$
less than one at least for one fit, then our model is preferred in 3 cases
(NGC2366, NGC2903, NGC2976), Einasto one in 5 cases (NGC2403, IC2574, NGC3621,
DDO154, NGC7331), and Burkert fit for NGC3198 is also significant. It is
remarkable that neither Navarro–Frenk–White model, nor pseudo-ISO model is
preferred against others in the THINGS sample.

If we take into account only two-parametric models, then we have the same 5
aforementioned cases for our model to be the best plus the case of NGC5055,
NFW model is preferred in 2 cases (NGC2403, NGC6946), pseudo-ISO one in 3
cases (NGC3621, DDO154, NGC7331), and Burkert model in 5 cases (NGC3031,
NGC3198, IC2574, NGC4736, NGC7331). Taking only significant fits, we get the
best performance of our model in 3 cases (NGC2366, NGC2903, NGC2976), NFW in 1
case of NGC2403, pseudo-ISO in 2 cases (DDO154, NGC7331), and Burkert in 2
cases (NGC3198, IC2574). From this we can conclude that our model is the best
two-parametric model of the set considered.

It is interesting to make pair comparison of our model with Burkert profile:
it is preferred statistically in 6 cases (NGC2366, NGC2841, NGC2903, NGC2976,
NGC3521, NGC5055) and is disfavored in 10 cases. However, in all these cases
besides IC2574 the preference is only marginal.

_NGC2366_ | NGC2403 | _NGC2841_
---|---|---
|  |
_NGC2903_ | _NGC2976_ | NGC3031
|  |
_NGC3198_ | IC2574 | _NGC3521_
|  |
NGC3621 | NGC4736 | DDO154
|  |
NGC5055 | NGC6946 | NGC7331
|  |
| NGC7793 |
|  |
Figure 4: Rotational curves $v_{r}(r)$ for some galaxies from THINGS survey
together with fits. Blue thick curves show best fits by dark matter
distributions of semidegenerate configurations, while magenta thin curves show
Einasto profile best fits. Galaxy names where semidegenerate profile fits
rotation curves better than Einasto profile are emphasized, and names where
profiles are comparable in fit quality are bolded. See digital version for
colored plots.

It should be mentioned that the velocity profile of the semidegenerate
configuration used for fitting is exact only in the case of Dark Matter
domination and thermal equilibrium at all radii, that can be not the case of
real galaxies. However, it is especially interesting that even such a
simplified model provides good correspondence to empirical rotational curves.

##  6 Conclusion

It follows from the results of fitting that the semidegenerate fermionic
distributions can fit dark matter in the THINGS sample of galaxies at least as
well as other profiles considered in the literature, with the important
”revenue” that this profile is theoretically motivated, and is not
phenomenological as most of the others. The cases when Einasto profile fits
rotational curve much better than semidegenerate profile show the general
cuspy behaviour of dark matter distribution, possibly representing a special
class of galaxies that are still not completely relaxed.

While Einasto profile is a pure phenomenological one based on best fit of the
observational and numerical simulation data, our profile is derived from the
first principles and based on the General Relativistic treatment of self-
gravitating neutral fermions. There is the distinct possibility that our
treatment gives the conceptual physical motivation for the existence of the
cored Einasto profile directly from the structure of the microphysical
constituents of dark matter.

###### Acknowledgements.

The authors wish to thank Prof. Herman Mosquera Cuesta for critical reading of
the manuscript and valuable suggestions.

## References

  * Argüelles et al. (2014) Argüelles, C. R., Fraga, B., & Ruffini, R. 2014, Journal of the Korean Physical Society (submitted)
  * Argüelles et al. (2013) Argüelles, C. R., Ruffini, R., Siutsou, I., & Fraga, B. 2013, Journal of the Korean Physical Society (submitted)
  * Binney & Tremaine (1987) Binney, J. & Tremaine, S. 1987, Galactic dynamics (Princeton University Press)
  * Chemin et al. (2011) Chemin, L., de Blok, W. J. G., & Mamon, G. A. 2011, AJ, 142, 109
  * de Blok (2010) de Blok, W. J. G. 2010, Advances in Astronomy, 2010
  * de Blok et al. (2008) de Blok, W. J. G., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2648
  * Einasto (1965) Einasto, J. 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87
  * Einasto (2013) Einasto, J. 2013, Brazilian Journal of Physics, 43, 369
  * Einasto & Haud (1989) Einasto, J. & Haud, U. 1989, A&A, 223, 89
  * Gao et al. (1990) Gao, J. G., Merafina, M., & Ruffini, R. 1990, A&A, 235, 1
  * Gao et al. (2012) Gao, L., Navarro, J. F., Frenk, C. S., et al. 2012, MNRAS, 425, 2169
  * Ingrosso et al. (1992) Ingrosso, G., Merafina, M., Ruffini, R., & Strafella, F. 1992, A&A, 258, 223
  * Klein (1949) Klein, O. 1949, Reviews of Modern Physics, 21, 531
  * Merritt et al. (2006) Merritt, D., Graham, A. W., Moore, B., Diemand, J., & Terzić, B. 2006, AJ, 132, 2685
  * Munshi et al. (2013) Munshi, F., Governato, F., Brooks, A. M., et al. 2013, ApJ, 766, 56
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Navarro et al. (2004) Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039
  * Ruffini et al. (2013) Ruffini, R., Argüelles, C. R., Rueda, J. A., & Siutsou, I. 2013, Phys. Rev. Lett.(to be submitted)
  * Schwarz (1978) Schwarz, G. E. 1978, Annals of Statistics, 6, 461
  * Shapiro & Teukolsky (1983) Shapiro, S. L. & Teukolsky, S. A. 1983, Black holes, white dwarfs, and neutron stars: The physics of compact objects (Wiley-Interscience)
  * Tolman (1930) Tolman, R. C. 1930, Physical Review, 35, 904
  * Walter et al. (2008) Walter, F., Brinks, E., de Blok, W. J. G., et al. 2008, AJ, 136, 2563