\section{Introduction}
\label{sec:int}

Dwarf spheroidal galaxies (dSphs) and globular clusters (GCs) have
similar luminosities but are strikingly different in terms of their
dark-matter (DM) content, as inferred from the kinematics of their
stars: while dSphs (and their ``ultra-faint'' extension to low luminosity) 
galaxies are the most DM dominated systems, there is
virtually no evidence that GCs have DM halos
\citep[][]{Heg96,Mas05a,Bra11}.  This finding might represent a
problem for standard cold dark matter cosmology, in which GCs are
expected to form within their own DM halos \citep{Pee84}. However it
is quite possible that GCs were originally embedded in massive DM
halos that evolved during the cluster lifetime as a consequence of
interaction with the cluster stars \citep{Bau08} and stripping by the
tidal field of the host galaxy \citep{Mas05b,Sai06}. Therefore, the
central parts of GCs might be left relatively poor in DM, because the
DM at the present time either populates the outer region of the
cluster or has been stripped \citep{Bek12}.  However, this is unlikely to be the
case in GCs observed to have stellar tidal tails, which argue against
very extended massive DM components \citep{Moo96}.  While it is not
well established to what extent two-body relaxation and tidal
stripping are effective in shaping hypothetical DM halos of GCs, it is
clear that these effects are minimized in clusters with long two-body
relaxation time and that orbit the outer parts of the host galaxy,
where tidal effects are small. The prototype of these systems is the
remote massive GC NGC~2419, located in the halo of the Milky Way at a
distance of $d=87.5\, {\rm kpc} $ \citep{DiC11} (Galactocentric distance of $\sim 95 \, {\rm kpc} $), 
which has half-mass
relaxation time $\simeq 43\, {\rm Gyr} $ (\citealt{Har10} --- 2010 edition). If NGC~2419 formed
within a massive DM halo, it is reasonable to expect that such a halo
has maintained its properties almost unaltered up to the present time.

In this paper we use the recently published kinematic data set of 166
stars \citep[][hereafter Paper I]{Iba11a} and the analysis of the
stellar mass-to-light ratio $M_*/L$ \citep{Bel12} of the globular cluster
NGC~2419 to study the DM content of this system.  
Paper~I showed that the present kinematic and structural data
on NGC~2419 can be fit acceptably well with models containing no DM; 
but this does not necessarily imply that models with DM cannot be as good.
Previous studies
have tried to put upper limits on the DM mass of NGC~2419
\citep{Bau09,Con11}. \citet{Bau09} have argued that the observed
velocity dispersion profile of the cluster (based on 40 stars) is not
consistent with a DM halo more massive than $M_{\rm dm}=10^7M_\odot$ (assuming
an isotropic stellar velocity distribution and NFW DM halo).
\citet{Con11} have put more stringent constraints ($M_{\rm dm}<10^6M_\odot$),
based on the observed outermost light profile of the cluster. However,
the argument of \citet{Con11} relies on the hypothesis that the
cluster has relaxed via two-body relaxation, which is not thought to
be the case for NGC~2419.

Given the collisionless nature of NGC~2419 it is natural to explore,
besides models with isotropic velocity distributions, 
also models with radially or tangentially biased stellar orbits. 
Of course it is important to exclude very radially anisotropic systems, 
which are radial-orbit unstable. However, in
this respect it must be noted that stellar systems with massive DM
halos can sustain relatively more radial kinetic energy than
corresponding purely baryonic systems \citep{Nip02,Nip11}.

In Paper~I we have shown that it is very problematic to explain
the observed structure and kinematics of NGC~2419 if modified
Newtonian dynamics (MOND; \citealt{Mil83}) is the correct theory of
gravity (but see \citealt{San12} and \citealt{Iba11b}).  It must be stressed that this
result does not necessarily imply that NGC~2419 does not contain a
substantial amount of DM. In fact, what we can infer on the basis of
the \citet{Iba11a} results is that it is unlikely that NGC~2419 has a
MOND-equivalent DM halo, i.e. a DM halo with density distribution such
that the overall gravitational potential of the cluster is the same as
predicted by MOND \citep[see][and references therein]{Mil86,Nip11}. 
Of course, there is no reason to expect that the putative DM halo of
this globular cluster must mimic a MOND gravitational potential, so
it is quite possible that there is a significant amount of DM
in NGC~2419, but that this DM is not distributed as in a
MOND-equivalent halo\footnote{It must be noted that stellar systems with
MOND-equivalent halos can sustain more radial anisotropy in the
central regions than MOND stellar systems \citep{Nip11}, so, at least in
principle, it is also possible that the cluster is well represented by
very radially anisotropic models within a MOND-equivalent halo.}.

Throughout the paper we adopt the above distance to NGC~2419,
which implies an angular scale of
$25.452\, {\rm pc} $ per arcmin.  Our conclusions are virtually independent of
this specific choice, because the estimated uncertainty on the
distance is small ($\sim 4\%$; \citealt{DiC11}). We note that our
results on the dark-to-luminous mass ratio of NGC~2419, obtained for
given $d$ and stellar mass to light ratio ($M_*/L_V$), hold for all the combinations of $d$ and
$M_*/L_V$ such that $(M_*/L_V) \times d$ has the same value. In the present
work we account for the uncertainty on $M_*/L_V$, which is of the
order of $30\%$ \citep{Bel12}, so it is clear that our results are not
affected by the smaller uncertainties on $d$. Other physical
properties of the cluster are listed in Table~1 of Paper~I.


\section{Models without dark matter}
\label{sec:models_without_DM}

\subsection{Distribution function based models}
\label{sec:dynamical_models_noDM}

Before analysing the case for the presence of dark matter in NGC~2419,
it is instructive to consider models that possess only stars. 
In Paper~I (see also \citealt{Zoc12}), we 
presented several Michie model fits to the system in Newtonian gravity, and we 
showed that one of these models gives a statistically excellent representation of the 
cluster. However, due to high computational cost, only a small grid of 11 models 
were investigated in that contribution. To probe the parameter space much more finely,
we re-wrote our Michie model code to allow it to be run on a parallel computer,
thereby enabling us to calculate a large ensemble of Michie model fits, using a Markov-Chain Monte Carlo 
algorithm. This has the advantage of allowing us to 
properly incorporate the measurement uncertainties and also yields posterior probability
distribution functions for the fitted parameters. Furthermore, we are now able to use the 
stellar mass to light ratio measurement of \citet{Bel12} to constrain the solutions.

The integration of the Michie models (optionally residing within within a dark matter
potential) is outlined in Appendix A. 
We initially set the dark matter potential $W_{\rm dm}=0$ so as to suppress that component, leaving
5 free fitting parameters: the Michie model stellar central density, the core radius,
the anisotropy radius, the central potential and the stellar mass to
light ratio\footnote{In Paper~I, we analysed the effect on the line of sight 
velocity distribution of the presence of binaries in the cluster. The most likely binary fraction
was found to be modest ($\lower.5ex\hbox{\ltsima} 10\%$), so we decided to neglect this effect in the
present contribution.}. For a given set
of these 5 parameters, we construct a simulated surface brightness
profile and the line of sight velocity distribution as a function of
radius, which can be compared directly to the observations.

As in Paper~I, we convolve the model velocity distributions
projected onto the line of sight by the estimated Gaussian errors
associated to each kinematic datum $i$, to give a probability
distribution $f_p(v,R_i)$ (we drop the model superscript $j$ used in
our earlier contribution, since we now consider a continuum of models
defined by the above 5 parameters). The likelihood of a model then
becomes:
\begin{equation}
\ln L = \sum_{i} \ln[f_p(v_i,R_i)] - \sum_{k} {{(\Sigma_p - \Sigma_k)^2}\over{2 \delta\Sigma_k^2}} + \ln \Big[P\Big({{M_*}\over{L_V}}\Big)\Big]  \, ,
\label{eqn:likelihood}
\end{equation}
where the second summation over the ($k=15$) surface brightness
profile data points is included to take account of the likelihood of
the density profile, and the third term is the contribution from the
observed stellar mass to light $(M_*/L_V)$ ratio. 

In their analysis of recent WFC3 data of this cluster, \citet{Bel12}
discuss in detail why the stellar mass to light ratio
$(M_*/L_V)$ must be constant, independent of radius. This is deduced 
from the uniformity of the luminosity function throughout the cluster,
confirming an earlier analysis \citep{Dal08} that reached the same conclusions
based on the lack of radial segregation of the Blue Straggler Star (BSS)
population\footnote{The argument that justifies a constant stellar 
$M/L$ with radius is as follows: the BSS stars, which are 
significantly more massive than any evolving cluster star, are not observed
to be radially segregated in NGC~2419, i.e. they have the same radial
distribution as single stars. This should not be observed if mass
segregation is at work, hence 2-body relaxation is ineffective in this cluster.
Hence the stellar mix is the same at any radius, and therefore the stellar
mass-to-light is the same at any radius. This argument is valid
because mass segregation is expected (and observed)
in most classical globular clusters, where 2-body encounters 
are a relevant mechanism of energy-exchange (relaxation). In
less dense systems such as dwarf galaxies, the BSS population \citep{Map09} 
would not allow us to draw this conclusion, as they are fully non-collisional.}.
This justifies the simplicity of the second term on the right hand side of Equation~\ref{eqn:likelihood}:
the surface brightness profile should trace the projected mass density.


\begin{figure*}
\begin{center}
\vbox{
\hbox{
\includegraphics[bb= 20 20 360 350, angle=0, clip, width=8.5cm]{N2419_DM_fig01a.eps}
\includegraphics[bb= 20 20 360 350, angle=0, clip, width=8.5cm]{N2419_DM_fig01b.eps}}
\hbox{
\includegraphics[bb= 20 20 360 350, angle=0, clip, width=8.5cm]{N2419_DM_fig01c.eps}
\includegraphics[bb= 20 20 360 350, angle=0, clip, width=8.5cm]{N2419_DM_fig01d.eps}}
}
\end{center}
\caption{The surface brightness profile (a) and velocity dispersion profile (b) of the 
best-fit models. As described in the
text, the observed data were first fit with a pure Michie model (5 free parameters) using
an MCMC algorithm that explored this parameter space. 
The blue line in (a) shows the resulting
excellent fit to the surface brightness profile. The model
simultaneously also fits the kinematics, which are also shown in blue on the
right-hand panel. The big hexagonal markers show a binned
representation of the kinematic data, presented here {\em only} for
visual purposes. The likelihood calculation does not use these
binned values, but rather the individual kinematic measurements, the
radial locations of which correspond to the radial locations of the
small dots. The corresponding velocity dispersion of the model
provides the ordinate value of these points.
The best two-component model consisting of a Michie 
model to represent the stellar component and a double
power law model to describe the dark matter is displayed in red. These values
are almost coincident with the single (pure stellar) Michie model solution.
For this composite model, a total of nine free parameters were fitted by
the MCMC algorithm. The green and turquoise dots show the best non-parametric 
solutions derived with the MCMC Jeans equation solver. The model in green
has no DM, while the turquoise model allows for a DM component.
Panel (c): the surface brightness residuals between the data and the models
 (a small shift in radius has been applied to avoid the superposition of the
points). Panel (d) shows the $\beta$ anisotropy profile of the models.
The best non-parametric model without DM (green) was derived from $1.2\times10^7$ MCMC trials, 
probing a vast parameter space of solutions; the fact its $\beta$ profiles is very similar
to that of the pure Michie model (blue) lends strong support to the use of the Michie 
model as a good representation of this globular cluster.}
\label{fig:Michie}
\end{figure*}


Having fixed the shape of the density
profile, we allow for variation of its normalisation (and that of the
total stellar mass $M_*$) by leaving 
$M_*/L_V$ as a free parameter. \citet{Bel12} found a best-fit value of
$M_*/L_V=1.5\pm 0.1$, but they argue that uncertainty in stellar
evolution models allows for a range of $1.2 \la M_*/L_V \la 1.7$ in the
most likely value. For the probability density function $P$, we therefore
adopt a flat distribution between $1.2 \le M_*/L_V \le 1.7$, tapered as
a Gaussian of sigma 0.1 beyond these edges.

We investigated this model parameter space with a Markov-Chain Monte
Carlo (MCMC) approach, using the Parallel Tempering technique
\citep{Gre05}. The Parallel Tempering procedure explores the parameter
space at different ``temperature'' (i.e. smoothing) levels, which
ensures that the solutions do not get stuck in local maxima. In our
implementation we used 4 parallel MCMC chains, allowing swaps between
the chains at random intervals of approximately 100 steps (see
\citealt{Gre05} for a detailed description of this method). To ensure rapid 
convergence, the proposal step sizes of all the parameters together were 
refined (in lock-step) every 1000 iterations, while every 10000 steps the proposal step size 
of all the parameters were refined separately. We set a target acceptance ratio of 25\%,
so that on average one in four trial steps was accepted (this is considered to be optimal
when fitting a large number of parameters, see \citealt{Gre05}).

It is well known that systems with strong radial anisotropy
are prone to the radial-orbit instability.  Typically, stability
criteria are expressed in terms of the Fridman-Polyachenko-Shukhman
parameter $\xi \equiv 2 T_{\rm rad}/T_{\rm tan}$, where $T_{\rm rad}$ and $T_{\rm tan} \equiv
T_{\vartheta}+T_{\varphi}$ are the radial and tangential components of the kinetic
energy tensor, respectively \citep[][]{Pol81,Fri84}.  Here we find it
convenient to use the parameter $\xi_{\rm half}$ which measures the same
ratio, but within the half-mass radius \citep{Nip11}.  In
paper~I we found that, in the absence of DM, $\xi_{\rm half}=1.5$ can
be taken to separate stable and unstable models of NGC~2419. This
condition is used as a prior to reject unstable solutions.

The resulting best fit from 100000 iterations of the lowest temperature chain is
shown with the blue line and blue dots in Fig.~\ref{fig:Michie}, and is compared to the observational 
data. The corresponding Michie model parameter values and their uncertainties are listed in Table~\ref{tab:Michie}.
The solution is fully compatible with the best solutions found in Paper~I and \citet{Iba11b}.


\subsection{Jeans equation based models}
\label{sec:kinematic_models_noDM}

While the Michie model gives a statistically very good description of the cluster (as we have shown in Paper~I), 
one is nevertheless left wondering how unique this solution is. In that earlier contribution we showed that 
isotropic Michie models (commonly referred to as King models) cannot account for the 
observations, but one could attempt at this stage to fit alternative models such as a Plummer, a polytrope or some other model.
However, Nature need not follow these simple analytic forms.
To attempt to survey the vast space of more general solutions, we have developed an
algorithm that takes a different approach: we will assess instead the statistical properties
of non-parametric solutions to the Jeans equation.

For a spherically-symmetric system the Jeans equation \citep{Bin08}:
\begin{equation}
{{G M(r)}\over{r }} = - \sigma_r^2 \Bigg[ {{\mathrm{d} \ln \rho}\over{\mathrm{d} \ln
        r}} + {{\mathrm{d} \ln \sigma_r^2}\over{\mathrm{d} \ln r}} + 2\Big(1 -
    {{\sgth^2}\over{\sigma_r^2}} \Big) \Bigg] \, ,
\label{eqn:Jeans}
\end{equation}
relates $M(r)$ the cumulative mass inside a radius $r$ to the square of the velocity dispersion profiles
in the radial ($\sigma_r^2$) and tangential ($\sgth^2$) directions and to the 
logarithmic derivatives of $\sigma_r^2$ and of the density ($\rho$).

The algorithm we have developed explores these functions from the cluster centre out to the tidal radius.
To sample the functions in an economical manner we found it convenient to assign one third of the radial
bins to probe the inner region within two core radii with logarithmic bins. Beyond two core radii the 
remainder of the radial points used linear intervals. We set the value of the central velocity dispersion
in the radial direction to be a model parameter. The radial 
velocity dispersion profile is then defined by the sequence $\Delta \ln \sigma_r^2 (r_i)$, and the density profile is defined
by the sequence $\Delta \ln \rho (r_i)$, where $r_i$ are the locations of the radial bins. To interpolate these profiles into
logarithmic derivatives, we assume that the functions can be approximated locally (between adjacent radial 
bins) by parabolas in $\ln \sigma_r^2$ (or $\ln \rho$) versus $\ln r$. It is worth noting that it is
highly undesirable to define the profiles directly in $\sigma_r^2$ and $\rho$, since this requires
extremely high accuracy to avoid huge interpolation errors in the logarithmic derivatives.
The remaining parameters for the 
dark matter free case, are the total cluster mass $M_*$ and the stellar mass to light ratio $M_*/L_V$.

Our model therefore has 3 parameters unrelated to the chosen binning ($\sigma_r^2(0)$, $M_*$, $M_*/L_V$), plus
$(n-1)$ parameters that define the dispersion profile and $n$ that define the density profile. 
We found it convenient to set the number of radial bins to $n=1+2^m$, and certainly with $n=129$ the profiles
are smooth and we could reconstruct a Michie model very accurately from artificial data. 

The algorithm projects the 3-dimensional functions onto the line of sight, to produce a velocity dispersion profile
and a surface brightness profile, from which the model likelihoods are measured using the observed data 
(Eqn.~\ref{eqn:likelihood}). 
One further assumption is made in the calculation of the likelihood: we assume that the velocity distributions
are Gaussian. This assumption is not necessary for the algorithm to work, and we could have used
the velocity moments of the kinematic data directly, or indeed we could have probed other parametrisations
of the velocity distribution. However, it is computationally faster and avoids the errors deriving 
from the calculation of the moments from the data. Furthermore, 
the assumption is supported by the data: as we showed in 
Paper~I, the observed distribution is consistent with a Gaussian function at all radii probed.
The line of sight velocity distributions are then convolved with the uncertainty distribution appropriate for each data
point.

We explore this model ($260$ free parameters) using a MCMC 
scheme\footnote{The Bayesian method we are using yields probability distribution functions
for the fitting parameters, as well as information on the correlations between these parameters. In this context it
is not a problem that we have more fitting parameters than data points.}. As before, we use the Parallel
Tempering technique with 4 chains to probe progressively more smoothed versions of the likelihood surface.
Following \citet{Gre05}, this is achieved by multiplying $\ln(L)$ by a factor $1/m^2$, where $m$ is the chain number.
The $m=1$ chain corresponds to the lowest ``temperature'' simulation, which is the one of primary interest. 

We implemented an improvement on this method, reducing the dimensionality of the high temperature
chains by redefining their radial profiles at only $n=1+2^m$ points. Cubic interpolation (in $\ln \sigma_r^2$ and $\ln \rho$ 
versus $\ln r$) between these anchor points was used to fill in intermediate values to maintain 129 radial 
bins in all chains. Although we attempted a simulation including a $m=5$ chain (i.e. with 9 radial anchor points), 
it was clear that the interpolation was too crude to be useful.

Some regularisation is required to tell the algorithm that it should avoid exploring very irregular profiles. This was implemented
by using a mild smoothing prior on the profiles of the logarithmic derivatives. 


An initial simulation with $2\times10^6$ iterations for each chain was run using a Metropolis-Hastings
acceptance criterion. As before, all parameters were tuned in lock-step every 1000 iterations, and then 
refined independently every $5\times10^5$ iterations, so as to ensure an acceptance criterion of 25\%.
Finally, the 1000 best solutions from this Metropolis-Hastings algorithm were used as inputs for another
$2\times10^6$ iteration run using an affine-invariant ensemble sampler algorithm. Our algorithm
is based on the ``stretch move'' algorithm proposed by \citet{Goo10}, which we expand to use the 
above-described multiple ``temperature'' parallel chains. The simulation uses a population of 1000
``walkers'' (which are non-coincident parameter space positions) which evolve from one iteration
to the next, and so sample efficiently the likelihood function.
The convergence of the method was checked by running simulations from 3 very different starting positions.

While the models built from the distribution function (described in
Section~\ref{sec:dynamical_models_noDM}) are consistent by construction, 
this Jeans-equation based approach can yield solutions
that do not have physically possible distribution functions. This is a
fundamental drawback of the technique, and cannot be avoided. 
However, some unphysical models are eliminated by ensuring that the solutions
obey the Global Density Slope-Anisotropy Inequality \citep[GDSAI;][]{Cio10a,Cio10b},
which is the requirement that at each radius the negative of the
logarithmic slope of the stellar density distribution
$\gamma_*\equiv -{\rm d}\ln \rho_*/{\rm d}\ln r$ cannot be less than twice the
local value of the anisotropy parameter, $\gamma_*(r)\geq 2\beta(r)$.
We remind the reader that the anisotropy parameter is defined as:
\begin{equation}
\beta(r) \equiv 1 -{\sgthsq+\sgphsq\over 2\sgrsq},
\end{equation}
where $\sgr$, $\sgth$ and $\sgph$ are, respectively, the $r$,
$\vartheta$ and $\varphi$ components of the velocity-dispersion
tensor. As in Section~\ref{sec:dynamical_models_noDM}, we also
impose the requirement that $\xi_{half} < 1.5$ to ensure stability of the resulting models.

In Fig.~\ref{fig:Michie} we display, in green, the surface brightness profile and 
velocity dispersion profile of the best solution to the Jeans equation found with this 
algorithm. The anisotropy profile of this solution is also shown (panel d); 
it is striking how closely this non-parametric solution matches
the behaviour of the best-fit Michie model from Section~\ref{sec:dynamical_models_noDM}
(shown in blue). 
That our Jeans equation solving algorithm should have selected this model, which is
so similar to the previously-fitted Michie model, out
of the vast space of possible models, lends strong support to the use of the Michie distribution function
to represent the stellar component in this system.
The small differences with respect to the Michie model fit in Fig.~\ref{fig:Michie} can be
attributed to the differences in the anisotropy profile but the
driving factor is the different stellar mass to light ratio (${M_*/L_V}=1.89$).


\section{Models with dark matter}
\label{sec:models_with_DM}

Following the layout of the previous section, we consider first adding dark matter
to the Michie model, and then freeing $M(r)$ in the Jeans
equation so as to also allow for dark matter in those solutions.


\subsection{Distribution function based models}
\label{sec:dynamical_models_DM}

We first consider the possibility that the putative dark-matter halo can be described by the double power law density profile
\citep{Bin08}:
\begin{equation}
\rho_{\rm dm} (r)=\frac{\rho_{\rm dm,0}}{\left(\frac{r}{r_{\rm s}}\right)^\gamma\left(1+\frac{r}{r_{\rm s}}\right)^{\delta-\gamma}},
\label{eq:twopower}
\end{equation}
where $\gamma$ and $\delta$ are the two power-law indices
characterising the shape of the profile, $r_{\rm s}$ is a scale radius, and
$\rho_{\rm dm,0}$ is a reference density.  In particular, when $\gamma=1$
and $\delta=3$ we get the \citet[][NFW]{Nav96} model; when $\gamma=0$
and $\delta=3$ we get a profile with a flat inner density
distribution, which is very similar to the profile proposed by
\citet{Bur95}. For combinations of $\gamma$ and $\delta$ such as those
mentioned above the integrated mass of the model diverges at large
radii, so it is useful to introduce as an additional parameter a
truncation radius $r_{\rm vir}$ such that $\rho_{\rm dm}=0$ for $r> r_{\rm vir}$.  The DM
mass contained within a radius $r$ is
\begin{equation}
M_{\rm dm}(r)=4\pi \int_0^{r}\rho_{\rm dm}(r')r'^2dr',
\end{equation}
so the total mass of the DM halo is $M_{\rm dm,vir}=M_{\rm dm}(r_{\rm vir})$.  If the
stellar density distribution is truncated at $r_{\rm t}$, the kinematics of
the cluster is influenced only by the mass distribution within 
the stellar tidal radius $r_{\rm t}$.
It follows that our results are not sensitive to the value of $r_{\rm vir}$,
provided that $r_{\rm vir}\geqr_{\rm t}$, which we assume throughout the paper (in
other words, our results are not sensitive to the halo concentration
$C\equivr_{\rm vir}/r_{\rm s}$). As a consequence, we define the quantity
$M_{\rm dm,t}\equivM_{\rm dm}(r_{\rm t})$ (the DM mass contained within $r_{\rm t}$), which is
more useful than $M_{\rm dm,vir}$ in order to characterise our models.  In
summary, our model halo is determined by specifying the following four
parameters: $\gamma$, $\delta$, $r_{\rm s}$ and $M_{\rm dm,t}$.


For the DM component, for any choice of the DM profile parameters ($\gamma$ and $\delta$)
the adimensional potential can be calculated by direct integration of equation~(\ref{eq_pois}) as
$$W_{\rm dm}(r)=9~
\frac{\rho_{\rm dm,0}}{\rho_{*,0}}\left(\frac{r_{\rm s}}{r_{\rm c}}\right)^{2}\int_{r/r_{\rm s}}^{r_{\rm t}/r_{\rm s}}\frac{1}{r'^{2}}\int_{0}^{r'}r''^{2}\rho_{\rm dm}(r'')dr''~dr'$$
In the particular cases of a NFW ($\gamma=1;~\delta=3$) and a cored
($\gamma=0;~\delta=3$) halo the above integral admits the analytical solutions:
\begin{equation}
\label{notev_eq}
W_{\rm dm}(r)=9~\frac{\rho_{\rm dm,0}}{\rho_{*,0}}\left(\frac{r_{\rm s}}{r_{\rm c}}\right)^{2}~\frac{\ln\left(1+\frac{r}{r_{\rm s}}\right)}{\frac{r}{r_{\rm s}}} 
\end{equation}
when $(\gamma=1,\delta=3)$, and
\begin{equation}
W_{\rm dm}(r)=-\frac{9}{2}~\frac{\rho_{\rm dm,0}}{\rho_{*,0}}\left(\frac{r_{\rm s}}{r_{\rm c}}\right)^{2}~\frac{\frac{r}{r_{\rm s}}-2\left(1+\frac{r}{r_{\rm s}}\right)~\ln\left(1+\frac{r}{r_{\rm s}}\right)}{\frac{r}{r_{\rm s}}\left(1+\frac{r}{r_{\rm s}}\right)}
\end{equation}
when $(\gamma=0,\delta=3)$.  

By incorporating this dark matter potential into the procedure
described in \S\ref{sec:dynamical_models_noDM}, one
can construct spherically symmetric
self-consistent stellar systems, truncated at a tidal radius
$\tilde{r}_{\rm t}$, in equilibrium in the presence of a DM halo \citep[see
  also][]{Amo11}. In these models the parameters $W_0$, $r_{\rm s}/r_{\rm c}$,
$\tilde{r}_{\rm a}$, $\rho_{\rm dm,0}/\rho_{*,0}$, $\gamma$ and $\delta$ drive the
shape of the density and velocity dispersion profiles, while the
values of $\rho_{*,0}$ and $r_{\rm c}$ determine their normalisation to
physical units (see Appendix).  

The stellar components of the models considered here are self-consistent by
construction, but in principle it is not guaranteed that there is a
positive distribution function generating the DM components. For
a wide class of two-component models, a necessary condition for
consistency is that also the DM component satisfies the GDSAI \citep{Cio10b}.
Clearly the specific form of the halo anisotropy
profile is not important for the purpose of the present investigation,
because the halo influences the observables only through its
gravitational potential. However, we impose that in our model the halo
density never increases with radius (${\rm d}\ln \rho_{\rm dm}/{\rm d}\ln r\leq0$ at
all $r$, i.e. $\gamma,\delta > 0$). This is a necessary condition for consistency only when the
halo is isotropic or radially anisotropic, but we impose it as a
general constraint because we think it is an astrophysically motivated
assumption. However, we noted that it is not excluded that
there are self-consistent tangentially biased models with
${\rm d}\ln\rho_{\rm dm}/{\rm d}\ln r<0$ at some $r$.

Even if consistent, some models can be excluded because they are
unstable. It has been shown that DM halo has a stabilising effect against the
radial-orbit instability \citep{Nip02,Nip11}, so at least some models
with DM halo are expected to be stable even if they have
$\xi_{\rm half}>1.5$.  However, a full stability analysis of two-component
models of NGC~2419 is beyond the purpose of the present paper, so here
we conservatively assume $\xi_{\rm half}\leq 1.5$ as a stability criterion
for all the models considered in the present work.  With this choice
we should effectively exclude unstable models.


We proceed to examine the relative likelihoods of the
two-component models built from a power-law dark matter halo and a
Michie model stellar component. As in \S\ref{sec:dynamical_models_noDM}, we
use a ``Parallel Tempering'' MCMC algorithm with 4 chains to explore
the parameter space. The composite model has 9 free
parameters: the stellar central density, the Michie model core radius,
the anisotropy radius, the central potential, the stellar mass to
light ratio, the dark matter central density, the DM characteristic
radius and the DM inner and outer power-law indices. For a given set
of these 9 parameters, we construct a simulated surface brightness
profile and the line of sight velocity distribution as a function of
radius, which can be compared directly to the observations.

Several (very weak) priors were imposed on the solutions. The stellar core radius
must lie within a factor of 10 of $10.88\, {\rm pc} $; the anisotropy radius
has to lie between 0.5 and 100 core radii;
the central potential parameter has to lie between 2 and 20; and the
stellar mass to light ratio is required to stay within the bounds of
$0.8<{M_*/L_V}<2.1$ (i.e. within $4\sigma$ of the flat region 
$1.2<{M_*/L_V}<1.7$ stated above in Section~\ref{sec:dynamical_models_noDM}). 
Furthermore, we impose that the characteristic
radius of the dark matter has to lie between $25\, {\rm pc} $ and $2.5\, {\rm kpc} $;
and that the inner and outer power law exponents lie in the ranges $0<
\gamma< 1.5$ and $2<\delta<4$. While certain black hole growth models
suggest that the inner power-law slope may be higher than $\gamma = 1.5$
(e.g., \citealt{Qui95}), this is likely to affect only the very central region of the
cluster where our data do not constrain the density profile.

Figure~\ref{fig:Michie} shows (in red) the most likely model found by the
MCMC algorithm in $10^6$ iterations. The solution has $r_{\rm c}=11.4\, {\rm pc} $,
$r_{\rm a}=1.6r_{\rm c}$, a dark matter characteristic radius of $28\, {\rm pc} $, and
$\gamma=0.06$ and $\delta=3.7$. However, we were surprised to find that
this model requires very little dark matter:
$M_{\rm dm}=5.0\times10^4M_\odot$ inside of 10\arcmin\ (the approximate limit
of the stellar structure), which is effectively insignificant compared
to the stellar mass of $8.1\times10^5M_\odot$. Marginalising over all
other parameters, we show the posterior distribution of $M_{\rm dm}$ in
Fig.~\ref{fig:composite_posterior}; clearly the data appear to favour
models without substantial dark matter.

To confirm this curious result, we re-ran the MCMC experiments forcing
certain astrophysically-motivated parameter sets: cored DM models
($\gamma=0$) and the NFW model ($\gamma=1$, $\delta=3$). In an almost
identical manner to the more general solution above, the model
likelihoods are found to decrease strongly with increasing DM
fraction. In an attempt to explore whether the MCMC algorithm was able
to escape from solutions with low dark matter, we also performed an
experiment forcing a massive DM component, by setting $M_{\rm dm}>
5\times10^6M_\odot$ as a strong prior. But again, the distribution of
solutions simply clustered near this low mass limit.

\begin{figure}
\begin{center}
\includegraphics[bb= 30 20 570 570, angle=0, clip, width=8.5cm]{N2419_DM_fig02.eps}
\end{center}
\caption{Unshaded histogram: the posterior distribution of the mass of the dark component
  inside of $10\arcmin$ (the approximate size of the stellar
  structure) derived by the MCMC algorithm that fits a Michie model
  inside a dark matter halo modelled as a double power-law. The hatched
  histogram shows the posterior distribution of the luminous component 
  of the cluster. Thus if we impose a
  Michie model on the luminous component (which as we have shown in
  Paper I, gives an excellent fit to the observations), the most
  likely solutions for the dark matter have very little mass relative
  to the stars. It would be natural to conclude that the cluster needs
  no dark component.}
\label{fig:composite_posterior}
\end{figure}

We conclude from this analysis that a simple Michie model where mass
traces light is actually preferred by the data over more complicated
models that include the possibility of having dark matter in the form
of a double power-law profile. The preference for a cored distribution
turns out to be very slight - the best models with $0.9 < \gamma < 1.1$
are only a factor of 1.4 less likely than the best model.
However, solutions with some DM seem better than those without DM. 
This is a manifestation of the tension between stellar and dynamical 
mass to light ratio that was noted in \citet{Bel12} (and also in \citealt{Sol12} 
in other clusters) and which is discussed further below.



\begin{figure}
\begin{center}
\includegraphics[bb= 30 20 570 570, angle=0, clip, width=8.5cm]{N2419_DM_fig03.eps}
\end{center}
\caption{Unshaded histogram: the posterior distribution of the mass of the dark component
  inside of $10\arcmin$, but this time calculated from an MCMC
  algorithm that solves the Jeans equation in a non-parametric way.
  A total of $1.2 \times 10^7$ MCMC iterations are presented. 
  The difference to Figure~\ref{fig:composite_posterior} is striking:
  with the freedom from analytic form, the preferred solutions
  possess a non-negligible dark matter component. The
  stellar mass of the models (hatched histogram) is approximately half 
  that of the dark matter.}
\label{fig:Jeans_posterior}
\end{figure}


\subsection{Jeans equation based models}
\label{sec:kinematic_models_DM}

The paucity of dark matter implied by our MCMC survey of the likelihood
surface of Michie models is compelling, since the Michie models form
a fairly general family of structures, but one is again left wondering to what 
extent this conclusion depends on the choices of the models
selected to examine both the cluster and the dark matter. 

For this reason we extended our analysis to a
more general set of spherical models where both the dark matter density
profile and the stellar density profile are free non-parametric profiles. 
We reuse the non-parametric MCMC tool
described in \S\ref{sec:kinematic_models_noDM} which solves the Jeans 
equation, finding solutions that are consistent with the data and our priors.


\begin{figure*}
\begin{center}
\hbox{
\includegraphics[bb= 30 20 570 570, angle=0, clip, width=8.5cm]{N2419_DM_fig04a.eps}
\includegraphics[bb= 30 20 570 570, angle=0, clip, width=8.5cm]{N2419_DM_fig04b.eps}
}
\end{center}
\caption{Panel (a): The cumulative mass profile of the dark matter component of
  the most likely model derived by the MCMC algorithm using the non-parametric Jeans
  equation method is shown with the continuous line. The dotted line
  shows the cumulative mass profile of the stellar component of the same solution.
  To first approximation the dark matter profile appears fairly similar to the stellar
  profile, and stands in stark contrast to cosmological expectations: the dashed line shows
  an NFW profile with $M_{vir}=10^7M_\odot$, concentration $c=23.5$ and $r_{vir}=4.5\, {\rm kpc} $. The dot-dashed 
  line shows an attempt to mimic the continuous line with an NFW profile; this however
  requires an extremely high concentration ($c=250$), even if $r_{vir}=13\arcmin = 0.33\, {\rm kpc} $.
  Panel (b): the corresponding density profiles.}
\label{fig:mass_profile}
\end{figure*}


However, given that we now allow for the possibility of dark matter, 
we require some additional information to close equation~\ref{eqn:Jeans}.
To this end we have chosen to free the tangential velocity dispersion profile (in addition
to $\sigma_r^2(0)$, $M_*$, ${M_*/L_V}$, and the profiles of 
$\Delta \ln \rho$ and $\Delta \ln \sigma_r^2$).
We found it convenient to tabulate this profile as $\Delta \ln \sgth^2$, adding
a further $n-1$ parameters to the MCMC fit, plus $\sgth^2(0)$. The model
thus has a total of 389 free parameters.

The MCMC simulations were run in an identical manner to those
described in \S\ref{sec:kinematic_models_DM}, using again 
a hierarchy of 4 different ``temperature'' chains, each with progressively lower
spatial resolution. As before, after initial runs with $2\times10^6$ 
iterations with a Metropolis-Hastings acceptance algorithm, we
switched to using a population of affine-invariant ensemble samplers and
reran the simulation for a further $2\times10^6$ iterations.

While the models built from the distribution function (described in
Section~\ref{sec:dynamical_models_DM}) are consistent by construction (at least
for as far as the stellar component is concerned), we cannot exclude that some of the  two-component models
considered are not generated by positive distribution functions. However, as in Section~\ref{sec:kinematic_models_noDM}
we impose at least some necessary conditions for consistency. We additionally
imposed the prior that the dark matter density profile decreases with radius
and that its logarithmic slope becomes increasingly negative with increasing radius, i.e. if
the radial locations $r_1$ and $r_2$ are such that $r_1 < r_2$, the logarithmic slopes
at these locations obey the condition:
$$ 0 > {{{\rm d}\ln \rho_{\rm dm}}\over{{\rm d}\ln r}} \Big|_{r_1} > {{{\rm d}\ln \rho_{\rm dm}}\over{{\rm d}\ln r}} \Big|_{r_2} \,.$$
This property is shared by many commonly-used models such as the Plummer sphere,
the Michie (and hence King) models, the double power law models if $\delta > \gamma$
(of which the NFW is a member), among others.


\begin{figure}
\begin{center}
\includegraphics[bb= 30 20 570 570, angle=0, clip, width=8.5cm]{N2419_DM_fig05.eps}
\end{center}
\caption{For the most likely non-parametric solution to the Jeans equation 
(shown previously in Figure~\ref{fig:mass_profile}), we calculate the stellar 
mass to light ratio required to eliminate the need for a dark matter component. The region between
the two dashed lines marks the maximum ($4\sigma$) range consistent with the analysis of \citet{Bel12}.}
\label{fig:ML}
\end{figure}

\begin{figure*}
\begin{center}
\vbox{
\hbox{
\includegraphics[bb= 20 20 360 350, angle=0, clip, width=8.5cm]{N2419_DM_fig06a.eps}
\includegraphics[bb= 20 20 360 350, angle=0, clip, width=8.5cm]{N2419_DM_fig06b.eps}}
\hbox{
\includegraphics[bb= 20 20 360 350, angle=0, clip, width=8.5cm]{N2419_DM_fig06c.eps}
\includegraphics[bb= 20 20 360 350, angle=0, clip, width=8.5cm]{N2419_DM_fig06d.eps}}
}
\end{center}
\caption{
The most likely Jeans equation solution of NGC~2419 using the new non-parametric MCMC technique
allowing for both a stellar and a dark matter component.
Panels (a) and (b) show the fit to the observed surface brightness and kinematics, respectively.
The data are reproduced from Figure~\ref{fig:Michie}. The small green dots in (a) display
the excellent model fit. The non-parametric velocity dispersion profiles in the radial and tangential
directions ($\sigma_r$ and $\sgth$) along with their line of sight projection ($\sigma_{LOS}$)
are displayed alongside the velocity dispersion estimates. We stress again that the {\em individual}
star velocities were used in the likelihood calculation. The 
remaining profiles required for the Jeans equation solver are displayed in the bottom panels: the {\it stellar}
density and its derivative are shown in (c), and the velocity dispersion derivatives in (d). 
For reference, we also show in panel (d) the anisotropy profile derived from the velocity 
dispersion profiles shown in (b). Note that the logarithmic derivative 
profiles are significantly noisier than the physical profiles. This model has a stellar mass to 
light ratio ${M_*/L_V}=1.20$, 
which is within the acceptable range reported by \citet{Bel12}. }
\label{fig:Jeans_composite}
\end{figure*}


Figure~\ref{fig:Jeans_posterior} shows the
marginalised posterior distribution of the dark matter in the solutions (unshaded 
histogram). The distribution peaks at $\approx 9.7\times 10^5 M_\odot$, suggesting the
presence of a significant dark matter fraction within the central 10\arcmin. 
The posterior distribution of the mass of the stellar component, shown with the hatched histogram,
reveals that roughly half the amount in stars as dark matter appears to be present in the 
solutions. At first sight it would appear that we are facing a simple degeneracy between the total
mass of the cluster and the mass to light ratio of the stellar population. This suspicion 
is reinforced by Figure~\ref{fig:mass_profile}a, which shows the cumulative mass profile
in the most likely solution: the dark matter (solid line) and stars (dotted line) follow each other
fairly closely (but see also the density profile in panel b). 
This best solution has a stellar mass to light ratio ${M_*/L_V}=1.20$, 
within the acceptable
range found by \citet{Bel12}. This would therefore suggest that the ``dark matter'' component is simply 
the result of the algorithm attempting to force a constant stellar ${M_*/L_V}$ over the whole radial range of
the cluster.

However, this option does not appear to be possible, since it would require a ${M_*/L_V}$ profile that is
inconsistent with observations and other astrophysical priors, as we show in Figure~\ref{fig:ML}. 
Here we have calculated the stellar ${M_*/L_V}$ that would be required to account for the total density at 
a given radius.
The constraints on the uniformity of the radial profile of the ${M_*/L_V}$ \citep{Bel12}, as well as
the implausible ${M_*/L_V}$ values beyond $\sim 5\arcmin$ make it unlikely that the dark matter component in our
Jeans equation solutions are simply due to erroneous assumptions in the ${M_*/L_V}$ profile.

The best Jeans equation solution is displayed in Fig.~\ref{fig:Jeans_composite}; the panels 
also illustrate the working of the method by showing the non-parametric profiles 
that are refined at each Markov-Chain Monte Carlo
iteration. 
As can be seen, the solution gives an almost perfect reproduction of the surface brightness
data, and interestingly, is able to fit the full range of the kinematic data, giving a
higher (and better) central velocity dispersion than any of the Michie models studied here or
in Paper~I. Beyond the central $\sim 1\arcmin$, the model is highly radially 
anisotropic\footnote{Recently, \citet{Zoc12} have studied NGC~2419 (among other clusters)
using so-called $f^{(\nu)}$ models, which can represent violently relaxed systems. In their
analysis (which assumes that no dark matter is present), they also manage
to produce a steep central velocity dispersion profile (see their Fig.~3).  Note
however that their best-fit $f^{(\nu)}$ model has ${M_*/L_V}=2.4$ and is quite
anisotropic $(\xi=1.77)$.}.


\section{Discussion}
\label{sec:Discussion}

With current data the answer to the question of whether the globular cluster NGC~2419
contains dark matter or not depends on the model (and hence one's priors) used to analyse the structure.
The choice of using a Michie model for the stars plus a double power law dark matter halo,
is appropriate for those who value highly dynamically self-consistent methods and who adhere 
strongly to cosmological expectations. 
The most likely solution of this composite model had $M_{DM}(<10\arcmin)=5\times 10^4M_\odot$, 
which would imply in standard
$\Lambda$CDM cosmology (following \citealt{Nav00})
a virial mass of $\sim 10^5M_\odot$, roughly an order of magnitude lower than the
mass of the stellar component in the cluster. 
The upper 99\% confidence limit for the dark mass inside of $10\arcmin$
according to the probability distribution function in Fig.~\ref{fig:composite_posterior} is $7.2\times 10^5M_\odot$, 
and this corresponds to a virial mass of $\sim 4\times 10^6M_\odot$.
Such a minuscule halo is some two orders of magnitude lower than expectations for the halo mass of dwarf 
spheroidal galaxies \citep{Wal12}, which have similar stellar content. 
Thus with this composite model the data prefers virtually no dark matter and one can rule out 
the possibility that the cluster is embedded at the centre of a dwarf-spheroidal-like mini-halo.

The issue of whether the cluster is located at the centre of the putative DM halo is important
for the validity of the analysis presented here. However, the case should in principle be rather similar
to that of the globular cluster M54 in the Sagittarius dSph, where the dynamical friction 
timescale for this similarly-massive globular cluster and its host halo to coincide is estimated to 
be $\sim 1.5\, {\rm Gyr} $ \citep{Bel08}. The expectation therefore is that NGC~2419 should 
reside at the DM halo centre, if it possesses one.

The alternative analysis, using our new Jeans equation solver yields a slightly
different answer. The strength of this method is that it does not require one
to force a particular analytic functional form on the data.
The disadvantage however, is that the resulting models are simply solutions to the Jeans 
equation, and so many of them are likely
to be unphysical (corresponding to a non-physical distribution function) or possibly unstable, 
but examining the ($\sim 10^7$)
models individually in the detail necessary to ascertain their
validity would be an immense task, far beyond the scope of this contribution. 

Nevertheless, taking the distribution in Figure~\ref{fig:Jeans_posterior} at face value gives a
99\% upper bound to the dark matter inside of $10\arcmin$ of $1.1\times10^6M_\odot$, which is somewhat
higher than the previous method. The strong difference is that the Jeans
equation solution has a 99\% lower limit at $7.5\times 10^5M_\odot$, effectively demanding more
dark matter than stars within $10\arcmin$. The difference with respect to
the DM-free models is that that the solutions are better able to follow the central enhancement
of the velocity dispersion. In order to confirm this suspicion, we re-ran the Jeans 
equation simulations from \S\ref{sec:kinematic_models_DM} (that include dark matter),
this time ignoring the most central 25 radial velocity measurements. The result was a solution very similar
to that found in \S\ref{sec:kinematic_models_noDM}, with only a small dark matter
fraction (mean total dark matter mass 
of $1.6\times10^5M_\odot$ inside of $10\arcmin$). This experiment shows that the difference in total dark matter
content between the distribution function model solution of \S\ref{sec:dynamical_models_DM} 
and the Jeans equation solution of \S\ref{sec:kinematic_models_DM} is due to the latter 
attempting to fit better the central velocity dispersion, and that it is those central
stars that are driving the requirement for a (small) dark matter fraction.
The MCMC algorithm accommodates this mass by changing primarily the anisotropy profile,
which can be seen to be negative at $r < 0.5\arcmin$ in Figure~\ref{fig:Jeans_composite}d (green dots).

It is interesting to consider which model is preferable: the best
composite Michie plus double power-law model, or the best model
resulting from the Jeans equation analysis? 
The traditional frequentist $\chi^2$ approach is unfortunately not of much
use due to the large number of parameters in the Jeans equation
solver (we would simply find a negative reduced-$\chi^2$, indicating that we have too
many parameters). However, answering this question properly with Bayesian
evidence is very challenging because the parameters used in the two
analysis methods are very different. Indeed, it would require
integrating the models over all parameter values, which would be
extremely computationally costly, and is beyond the scope of the present work.

In this paper we have used kinematics measurements of NGC~2419 to test
the hypothesis that this globular cluster is embedded in a massive DM
halo. Another dark component that can in principle be detected via
kinematics measurements is a central intermediate mass black hole
\citep{Bau05}.  However, the current kinematic sample does not allow one
to set significant constraints on the mass of this putative central
black hole, because it affects only the very central ($\lower.5ex\hbox{\ltsima} 0\scnd1$)
velocity-dispersion profile, which is barely probed by current
observations.


\subsection{Comparison with previous work}
\label{sec:comp}

\citet{Bau09} assume that the putative halo around NGC~2419 has a NFW profile (with
$r_{\rm s}=500\, {\rm pc} $) and use their data to limit the dark mass inside of 500\, {\rm pc} \ to
$M_{\rm dm}(r<500\, {\rm pc} )<10^7M_\odot$. Given their
model parameters, this corresponds to $M_{\rm dm}(r<10\arcmin)<3.9\times 10^6 M_\odot$.  
In the work of \citet{Con11}, they also assume
that the halo has a NFW profile (with $r_{\rm s}=250\, {\rm pc} $ and $r_{vir}=1\, {\rm kpc} $)
and conclude that $M_{\rm dm}(r<1\, {\rm kpc} )<10^6M_\odot$, which, for their model
parameters, gives a very tight $M_{\rm dm}(r<10')<2.4\times 10^4M_\odot$.  Our results for 
the two component models are in reasonable agreement with 
these previous results. The limits provided above, 
are more stringent than those of \citet{Bau09}, and unlike \citet{Con11}, our
method does not require the cluster to have suffered strong two-body relaxation, which as we have mentioned, 
should not have occurred in NGC~2419.

\begin{table*}
\begin{center}
\caption{Model parameters and derived quantities for the four families of
models investigated in this contribution. Where possible, we report the mean 
value and root mean square dispersion of the MCMC solutions. Due to
the construction of our Jeans equation algorithm, the stellar central density is expensive
to calculate at every MCMC iteration, so we report only the values of
the most likely solution. Furthermore, the posterior probability distribution for the 
dark mass in Fig.~2 is highly skewed, so we only
report the most likely value for that model.}
\label{tab:Michie}
\begin{tabular}{ccccc}
\hline\hline
  & Michie & Jeans without DM & Michie + DM & Jeans with DM \\
\hline
stellar central density ($M_\odot \, \, {\rm pc} ^{-3}$) & $57.3\pm4.8$ & $74.4$ (best value) & $59.2\pm10$ & $66.9$ (best value)\\
stellar core radius ($\, {\rm pc} $) & $12.6\pm0.6$ & ... & $10.8\pm0.5$ & ... \\
anisotropy radius ($\, {\rm pc} $) & $14.2\pm2.3$ & ... & $18.3\pm2.7$ & ... \\
central potential   & $3.6\pm0.4$ & ... & $4.7\pm1.1$ & ... \\
stellar $M/L$ & $1.72\pm0.11$ & $1.89\pm0.07$ & $1.62\pm0.17$ & $1.15\pm0.12$\\
Stellar mass ($10^5 M_\odot$) & $8.1\pm0.2$ & $9.0\pm0.1$ & $7.6\pm0.8$  & $5.5\pm0.6$ \\
Dark mass ($10^5 M_\odot$) & $0$ & $0$ & $0.5$ (best value)  & $9.7\pm0.8$ \\
\hline\hline
\end{tabular}
\end{center}
\end{table*}



\section{Conclusions}
\label{sec:con}


We have undertaken an analysis of current data on the kinematics,
structure and stellar mass to light ratio of the
Galactic globular cluster NGC~2419. This distant halo cluster is a
very interesting specimen to study, since it is massive and as such
may have been able to maintain a residual dark matter halo over the
course of its assimilation into the Milky Way. Thus NGC~2419 allows us
to examine the possibility of whether (some) globular clusters formed
along with (their now disrupted?) dwarf galaxy hosts within mini dark
matter halos.

Our earlier analysis (Paper~I) showed that a simple
Michie model gives an excellent fit to the state-of-the-art
measurements of the kinematics and surface brightness profile of the
cluster presented in that contribution. Here we extended that study to investigate whether composite
models with an additional dark matter component in the form of a
generic double power-law could also be accommodated. To this end we
implemented an MCMC parameter search for solutions of this composite
system. The Michie distribution function is found to be
highly restrictive; while it is as yet not possible to completely rule out 
the presence of dark matter in NGC~2419, the preferred dark matter
fraction is 6\% of the stellar mass within 10\arcmin\ (the approximate limit of the cluster),
and there is virtually no room for astrophysically
interesting amounts of dark matter that could hint at the formation
of clusters within dark matter halos.

In a major effort to assess the model-dependence of these results we developed
a new Markov-Chain Monte Carlo algorithm that solves the spherical Jeans equation using 
non-parametric profiles for the density and velocity dispersion profiles.
A very high number of free parameters is required to follow the gradients with sufficient resolution.
The code uses a hierarchy of resolutions so as to probe the parameter space
as efficiently as possible. These ideas (and the software) should be applicable in the future
to analyse other approximately spherical systems of interest, such as dwarf galaxies, elliptical galaxies, or galaxy clusters.

The more general Jeans equation solutions prefer a small dark matter
fraction, approximately twice the stellar mass content, and at face value
rule out the possibility that there is no dark matter. This last inference 
depends, however, on the validity of the conclusion by \citet{Bel12} that the stellar mass to
light ratio is invariant with radius.
For the most likely solution, the inferred dark matter profile has a very different
shape to an NFW model, so it is possible that this putative dark component, if real, 
is related to unexpected stellar remnants (the expected remnants were 
taken into account by \citealt{Bel12}) rather than cosmological dark matter
(note however, that the resolution of present-day simulations that include the relevant baryonic physics
do not allow reliable predictions for the profile of dark matter mini-halos of this mass).

In their analysis, \citet{Bel12} found that assuming a 100\% retention of all the dark remnants 
would have changed the best fit ${M_*/L_V}$ from 1.49 to 1.6. Hence in a
standard scenario, changing the fraction of remnant stars will not have a sufficiently
large effect to account for Fig.~\ref{fig:ML}.
On the other hand, as is pointed out in that contribution,
all the currently available (and highly uncertain) models invoked to explain the 
multiple populations in GCs envisage (a) a progenitor of the present-day cluster with a 
baryonic mass a factor $\lower.5ex\hbox{\gtsima} 10$ larger than today, (b) a preferential loss 
of first generation stars with respect to the second generation, (c) a strong 
difference in the radial distribution of stars of the two generations 
(the second one being more concentrated, as is generally observed; see \citealt{Lar11}).
This complex star formation history and strong variations in the cluster potential must have 
occurred in the first $<500\, {\rm Myr} $ of the cluster life. It can be conceived that at this stage 
both (i) an anomalous production and retention of dark remnants occurred, 
and (ii) a radial distribution of dark remnants different from that of the light 
was produced and then remained un-modified in that non-collisional environment.

Therefore if the observed dark matter is of baryonic origin it is likely that 
it has something to do with the origin and early enrichment of the multiple 
populations \citep{Coh11, Muc12} inhabiting NGC2419. On the other hand one may imagine that 
we see the small remnant of a larger non-baryonic cosmological halo that was 
nearly completely stripped away from the cluster at early times (see, e.g. \citealt{Pen08}). During the first 
approach to its peri-galacticon, it may have provided the potential well that 
allowed the cluster to retain ejecta from first generation stars.

Thus all the models we have explored turn out to have very little dark matter in their central
regions, and we rule out any connection with mini dark matter halos of the mass that
is associated with those in which dwarf galaxies are embedded. This conclusion,
though consistent with earlier work on this particularly promising cluster, and consistent
also with a large literature on more nearby systems (but which on the whole are clearly 
not simple isolated systems), remains puzzling. The connection between globular clusters 
and dwarf galaxy remnants is becoming more and more compelling \citep{Mac10},
so it is natural to expect that some clusters retained some of the dark matter of their former
host. It is possible that this problem will only be resolved conclusively by obtaining a larger
sample of distant halo clusters, but that will require observing further afield in M31,
and to be feasible will require the next generation of instruments.


\section*{Acknowledgements}
R.I. gratefully
acknowledges support from the Agence Nationale de la Recherche though
the grant POMMME (ANR 09-BLAN-0228).
C.N. is supported by the MIUR grant PRIN2008. 
AS acknowledges the support of INAF through the 2010 postdoctoral fellowship grant.
M.B. acknowledges the financial support of INAF through the PRIN-INAF 2009 grant assigned to the project {\em Formation and evolution of massive star clusters}, P.I.: R. Gratton. 

\section*{Apendix A}

The distribution function proposed by \citet{Mic63} has the following phase-space form:
\begin{eqnarray}
\label{eq_df}
f(E,L)&=&f_0~\exp\left(-\frac{L}{2\sigma_{\rm K}^{2}r_{\rm a}^{2}}\right)\left[\exp\left(-\frac{E}{\sigma_{\rm K}^{2}}\right)-1\right],\nonumber\\
f(r,v_{\rm r},v_{\rm t})&=&f_0~\exp\left[-\frac{v_{\rm t}^{2}}{2\sigma_{\rm K}^{2}}\left(\frac{r}{r_{\rm a}}\right)^{2}\right]~\times\nonumber\\
& &\left[\exp\left(-\frac{v_{\rm r}^{2}+v_{\rm t}^{2}}{2\sigma_{\rm K}^{2}}-\frac{\psi}{\sigma_{\rm K}^{2}}\right)-1\right],
\end{eqnarray}
where $E$ and $L$ are the energy and angular momentum per unit mass,
respectively, $\psi$ is the effective potential defined as the
difference between the cluster potential (including the contribution
of both stars and the dark matter --- the latter component will be added 
in Section~\ref{sec:dynamical_models_DM}) 
at a given radius $r$ and the potential at the
cluster tidal radius ($\psi=\phi-\phi_{\rm t}$), $f_0$ is a
normalization coefficient factor, $\sigma_{\rm K}$ is a reference velocity
dispersion, $v_{\rm t}$ and $v_{\rm r}$ are the radial and tangential component of
the velocity and $r_{\rm a}$ is the anisotropy radius: the velocity
distribution tends to be isotropic at $r<r_{\rm a}$ and radially anisotropic
at $r>r_{\rm a}$.

The density and the radial and tangential components of the velocity
dispersion can be derived by integrating the above distribution
function:
\begin{eqnarray}
\rho_*(r)&=& 4 \pi \int_{0}^{\sqrt{-2\psi}}\int_{0}^{\sqrt{-2\psi-v_{\rm r}^{2}}}
v_{\rm t}~f(r,v_{\rm r},v_{\rm t}) ~\mathrm{d} v_{\rm t} \mathrm{d} v_{\rm r},\nonumber\\
\sigma_{\rm r}^{2}(r)&=&\frac{4
\pi}{\rho_*(r)}\int_{0}^{\sqrt{-2\psi}} v_{\rm r}^{2} \int_{0}^{\sqrt{-2\psi-v_{\rm r}^{2}}} v_{\rm t}~ f(r,v_{\rm r},v_{\rm t})
~\mathrm{d} v_{\rm t} \mathrm{d} v_{\rm r},\nonumber\\
\sigma_{\rm t}^{2}(r)&=&\frac{4 \pi}{\rho_*(r)}\int_{0}^{\sqrt{-2\psi}}
\int_{0}^{\sqrt{-2\psi-v_{\rm r}^{2}}} v_{\rm t}^{3}~f(r,v_{\rm r},v_{\rm t})
~\mathrm{d} v_{\rm t} \mathrm{d} v_{\rm r}.\nonumber
\end{eqnarray}
It is convenient to work with the dimensionless quantities
\begin{eqnarray}
\zeta=\frac{v_{\rm r}^2}{2\sigma_{\rm K}^{2}}, & \eta=\frac{v_{\rm t}^2}{2\sigma_{\rm K}^{2}},\nonumber\\ 
\tilde{\rho}_*=\frac{\rho_*}{\rho_{*,0}}, & \tilde{r}=\frac{r}{r_{\rm c}},\nonumber\\
W=-\frac{\psi}{\sigma_{\rm K}^{2}}, & \tilde{r}_{\rm a}=\frac{r_{\rm a}}{r_{\rm c}},\nonumber
\end{eqnarray}
where $\rho_{*,0}$ is the central cluster density and
$$r_{\rm c}\equiv\left(\frac{9 \sigma_{\rm K}^{2}}{4 \pi G \rho_{*,0}}\right)^{1/2}$$
is the core radius (King 1966). The above equations can be then written 
as
\begin{eqnarray}
\label{eq_models}
\tilde{\rho}_*(\tilde{r})&=&\frac{\int_{0}^{W}\zeta^{-\frac{1}{2}}\int_{0}^{W-\zeta}e^{-\frac{\eta\tilde{r}^{2}}{\tilde{r}_{\rm a}^{2}}}(e^{W-\eta-\zeta}-1) 
~\mathrm{d}\eta ~\mathrm{d}\zeta}{\int_{0}^{W_0}\zeta^{-\frac{1}{2}}\int_{0}^{W_0-\zeta}e^{-\frac{\eta\tilde{r}^{2}}{\tilde{r}_{\rm a}^{2}}}(e^{W_0-\eta-\zeta}-1) 
~\mathrm{d}\eta ~\mathrm{d}\zeta},\nonumber\\
\sigma_{\rm r}^{2}(\tilde{r})&=&\frac{2\sigma_{\rm K}^{2}\int_{0}^{W}\zeta^{\frac{1}{2}}\int_{0}^{W-\zeta}
e^{-\frac{\eta\tilde{r}^{2}}{\tilde{r}_{\rm a}^{2}}}(e^{W-\eta-\zeta}-1) 
~\mathrm{d}\eta ~\mathrm{d}\zeta}{\int_{0}^{W}\zeta^{-\frac{1}{2}}\int_{0}^{W-\zeta}e^{-\frac{\eta\tilde{r}^{2}}{\tilde{r}_{\rm a}^{2}}}(e^{W-\eta-\zeta}-1) 
~\mathrm{d}\eta ~\mathrm{d}\zeta},\nonumber\\
\sigma_{\rm t}^{2}(\tilde{r})&=&\frac{2\sigma_{\rm K}^{2}\int_{0}^{W}\zeta^{-\frac{1}{2}}\int_{0}^{W-\zeta}
\eta e^{-\frac{\eta\tilde{r}^{2}}{\tilde{r}_{\rm a}^{2}}}(e^{W-\eta-\zeta}-1) 
~\mathrm{d}\eta ~\mathrm{d}\zeta}{\int_{0}^{W}\zeta^{-\frac{1}{2}}\int_{0}^{W-\zeta}e^{-\frac{\eta\tilde{r}^{2}}{\tilde{r}_{\rm a}^{2}}}(e^{W-\eta-\zeta}-1) 
~\mathrm{d}\eta ~\mathrm{d}\zeta},\nonumber\\
\end{eqnarray}
which can be solved once the potential $W(\tilde{r})$ is known.  This
last quantity can be written as the sum of the stellar and DM
contributions $W=W_{\rm *}+W_{\rm dm}$. Both contributions must obey
independently the Poisson equation
\begin{eqnarray}
\label{eq_pois}
{\bf \nabla}^{2}\phi_{i}&=&4 \pi G \rho_{i},\nonumber\\
\frac{{\rm d}^{2}W_{i}}{{\rm d}\tilde{r}^{2}}+\frac{2}{\tilde{r}}\frac{{\rm d} W_{i}}{d\tilde{r}}&=&-9\tilde{\rho_{i}} \, . \nonumber\\
\end{eqnarray}
Here, the subscript $i$ refers to either stars or DM. 
For the stellar component we leave
$W_0$ as a free parameter and solve equation~(\ref{eq_pois})
adopting the boundary conditions at the center
\begin{eqnarray}
\label{bound_eq}
W_{\rm *,0}&=& W_0-W_{\rm dm,0} \, , \nonumber\\
\frac{{\rm d} W_{\rm *,0}}{{\rm d}\tilde{r}}&=&0 \, .\nonumber
\end{eqnarray}

Equations~(\ref{eq_models}) and (\ref{eq_pois}) can be integrated
numerically to derive the three-dimensional density profile, and the
radial and tangential velocity dispersion profiles.  Finally, the
above profiles have been projected on the plane of the sky to obtain
the surface mass density
\begin{equation}
\Sigma_*(R)=2\int_{R}^{\tilde{r}_{\rm t}}
\frac{\tilde{\rho}_*\tilde{r} \mathrm{d} \tilde{r}}{\sqrt{\tilde{r}^{2}-R^{2}}}\nonumber
\end{equation}
and the line-of-sight velocity dispersion
\begin{equation}
\sigma_{\rm V}^{2}(R)=\frac{1}{\Sigma_*(R)}\int_{R}^{\tilde{r}_{\rm t}}
\frac{\tilde{\rho}_*\left[2\sigma_{\rm r}^{2}\left(\tilde{r}^{2}-R^{2}\right)+\sigma_{\rm t}^{2}R^{2}\right] \mathrm{d} \tilde{r}}{\tilde{r}\sqrt{\tilde{r}^{2}-R^{2}}} \, .\nonumber
\end{equation}

As described in Section~\ref{sec:dynamical_models_noDM}, the full
likelihood analysis performed here requires for each model the
knowledge of the distribution of projected velocities (not only its
second-order moment).  This function can be calculated as
\begin{eqnarray}
&&f_p(v,R)=\int_{0}^{r_{\rm t}}\frac{\rho_*(r) r}{\sqrt{r^{2}-R^{2}}}\int_{0}^{\sqrt{-2\psi-v^2}}\nonumber\\
&&\times\int_{0}^{\sqrt{-2\psi-v^2-v_{x}^2}} f(r,v_{\rm r}',v_{\rm t}')~{\rm d} v_{y} {\rm d} v_{x} {\rm d} r \, ,
\end{eqnarray}
adopting the following change of variables
\begin{eqnarray}
v_{\rm r}'&=&v~\sin\alpha +v_{x}~\cos\alpha,\nonumber\\
v_{\rm t}'&=&[(v~\cos\alpha -v_{x}~\sin\alpha)^{2}+v_{y}^{2}]^{1/2} \, ,\nonumber
\end{eqnarray}
where $\alpha=\arccos(R/r)$.



