\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)$.


\section{Introduction}
The series of cubic RX$_2$Zn$_{20}$ has been recently closely examined, where R is a rare earth atom and X is Fe, Co, Ru, etc.\cite{jia1,jia2,jia3,jia4,nasc,tian,wang,isi1}. R is located in the cubic symmetric site surrounded by 16 Zn atoms. The exchange interaction between two R atoms is weak because R atoms are diluted in this compound. In fact, in the case of X = Co or Ru, the Curie temperature is less than 10 K for any R. In the case of the Fe series, however, $T_c$ is significantly enhanced.
It has been suspected that the magnetism of Fe atoms is involved in this high $T_c$.
However, YFe$_2$Zn$_{20}$ is a nonmagnetic material, i.e., the Fe atoms in YFe$_2$Zn$_{20}$ are non-magnetic.
At present, what causes the high $T_c$ of RFe$_2$Zn$_{20}$ is not clear.
We closely examined the magnetization, the magnetic susceptibility and the specific heat of DyFe$_2$Zn$_{20}$ in various magnetic fields, and found experimentally unusual magnetic anisotropies, i.e., the large magneto-crystalline anisotropy at 2 K and a metamagnetic transition in the fields along the [100] direction,
which disappears in the temperatures above 30 K below $T_c$.
These anomalous features are analyzed based on the following hamiltonians; the crystalline electric field, the exchange interaction between two R atoms and the exchange interaction between Fe and R atoms.
We have taken into account the contribution of the 3d electrons of Fe atoms based on the simplified SCR theory.
In this paper, we show that the enhancement of $T_c$, the strong magnetic anisotropy, the metamagnetic behavior at low temperatures, the disappearance of the magnetic anisotropy below $T_c$ down to 30 K, and no anomaly of magnetic specific heat at $T_c$, etc. can be explained by our calculation.

\section{Experimental results}
Figure \ref{MvsH_3} shows the magnetization curves $M(H)$ in the magnetic fields $H$ along the three principal axes [100], [110] and [111] at 2 and 30 K. $M(H)$ at 2 K is significantly anisotropic, as has been shown by Jia {\it et al.}\cite{jia1}.
The easy axis of magnetization is the [111] direction.
At zero field, the values of residual magnetic moments $M$ along the [110] and [100] directions
approximately equal to the projection of $M$ to these directions;
the calculated ratio of $M$ along the  [111], [110] and [100] directions  is 1:$\sqrt{2/3}$:$1/\sqrt{3}$.
In the figure, a metamagnetic behavior of $M$ is seen in $H$ along the [100] direction in approximately 1.2 T.
The anisotropic behavior of $M(H)$ disappears at 30 K even though this temperature $T$ is much lower than $T_c$.
In the inset of Fig. \ref{MvsH_3}, the details of the field dependence of $M(H)$ are shown in $H$ along the [100] direction.
The field of the metamagnetic transition gradually decreases with increasing temperature,
and the metamagnetic behavior disappears at 25 K.
At 15 and 20 K, the metamagnetic behavior is not clearly observed,
but an unusual $M(H)$ remains in $H$ lower than 0.2 T.
\begin{figure}[t]
\centering
\includegraphics[width=75mm]{MvsH_3}
\caption{Magnetization curves $M(H)$ of DyFe$_2$Zn$_{20}$ at 2 K and 30 K in $H \parallel $  [100],  [110], and [111].
Inset shows the temperature dependence of $M(H)$ in $H \parallel $ [100]. }
\label{MvsH_3}
\end{figure}
\begin{figure}[t]
\centering
\includegraphics[width=75mm]{MvsT}
\caption{Temperature dependence of $M/H$ of DyFe$_2$Zn$_{20}$ in $H \parallel $  [100], [110], and [111].
Inset shows the temperature dependence of $H/M$. }
\label{MvsT}
\end{figure}

Figure \ref{MvsT} shows the temperature dependence of $M(T)/H$ in $H$ // [100], [110], and [111].
$H$ is applied at 0.1 T.
From this figure, it is apparent that the strong anisotropy of $M(T)$ at low $T$ rapidly shrinks at approximately 20 K and disappears above 30 K.
The Curie temperature $T_c$ is 53 K, which was determined by $M(T)$ curve in $H = 0.01$ T,
because the temperature change of $M(T)$ at $T_c$ is more clearly observed in low $H$ than that in higher $H$.
The inset of Fig. 2 shows the temperature dependence of $H/M$, which obeys the Curie-Weiss law.
The paramagnetic Curie temperature $\theta_p$ is 21.1 K, and the effective moment $\mu_{\rm eff}$ is 10.8 $\mu_{\rm B}$, which is close to the effective moment of the free ion of Dy atom, 10.6 $\mu_{\rm B}$.
The magnetism of DyFe$_2$Zn$_{20}$ seems to be originated from the magnetic ions of rare earth Dy atoms.
However the temperature dependence of $M(T)$ is unusual;
$M$ decreases linearly with increasing $T$, and the temperature change of $M$ at $T_c$ is not clear.

The specific heat $C$ of DyFe$_2$Zn$_{20}$ and YFe$_2$Zn$_{20}$ was measured.
The magnetic part $C_{\rm m}$ of DyFe$_2$Zn$_{20}$ was defined as
$C_{\rm m}=  C_{\rm DyFe_2Zn_{20}} - C_{\rm YFe_2Zn_{20}}$, where
any corrections due to the difference of molar weight
have not been made.
Our $C_{\rm m}(T)$ of DyFe$_2$Zn$_{20}$ was almost the same
as the one reported by Jia {\it et al.}\cite{jia1}.
The two characteristic features of $C_{\rm m}(T)$, which is not shown in the present paper,
are (1) a broad peak around 15 K, (2) an unusual and unclear change of $C_{\rm m}$ at $T_c$.

The specific heat $C_{\rm m}$ was measured in $H$ along the three principal axes,
and we have found anomalous behavior of $C_{\rm m}$ especially below 20 K.
Figure 3 shows the temperature dependence of $C_{\rm m}(T)$ in $H$ along the [100] direction.
In $H=0$, $C_{\rm m}(T)$ shows only a broad peak with the amplitude of 13 J/mol K at approximately 12 K.
In 0.2 T, however, a tiny thorn appears at 16 K,
and then the thorn grows into an apparent peak with increasing $H$, as shown in the figure.
Then the peak abruptly disappears in $H= $ 1.4 T.
The entropy $S_{\rm m}(T)$ reaches to 14.2 J/mol K ($R$ ln 5.6) at 20 K in $H=0$,
and  decreases to 9.0 J/mol K ($R$ ln 3) at 20 K in $H=8$ T, where $R$ is the gas constant.
\begin{figure}[t]
\centering
\includegraphics[width=75mm]{CmvsT100detail}
\caption{Magnetic part of the specific heat $C_{\rm m}(T)$ of DyFe$_2$Zn$_{20}$ in various $H \parallel $ [100].}
\label{fig:3}
\end{figure}
In the case of $H$ along the [111] direction,
contrary to $C_{\rm m}(T)$ in $H$ // [100],
$C_{\rm m}(T)$ steadily decreases without any anomalies by applying $H$,
and the broad peak is gradually weakened  with increasing $H$ and disappears in 8 T.

\section{Analyses and discussion}
We analyze the anomalous behavior of $M(T)$, $M(H)$ and $C_{\rm m}(T)$ of DyFe$_2$Zn$_{20}$
taking into account the magnetism of both the Dy and Fe atoms.
First, we estimate the contribution from the Fe atoms based on the experimental results of YFe$_2$Zn$_{20}$,
which is a so-called nearly-ferromagnet.
Figure 4  and the insets show the inverse magnetic susceptibility $1/\chi(T)$, the magnetization curves $M(H)$ and
the specific heat divided by $T$, $C(T)/T$.
Note that $M$ is the magnetization per one mole of YFe$_2$Zn$_{20}$, that is, per two Fe atoms.
\begin{figure}[t]
\centering
\includegraphics[width=85mm]{YFZ_2chivsT111_2}
\caption{Temperature dependence of the inverse magnetic susceptibility of YFe$_2$Zn$_{20}$ in $H \parallel $  [111]. Open circles and the solid line denote the experimental data and theoretical line, respectively.
Lower inset shows the temperature dependence of $M(H)$ curves at various $T$ in $H \parallel $ [111].
Solid lines denote the theoretical lines.
Upper inset shows the specific heat divided by $T$.
The line in this inset is a guided line for the eyes.}
\label{fig:4}
\end{figure}
$M(H)$ in the inset in Fig. 4 is experimental raw data,
but the $\chi$ is corrected by subtracting small ferromagnetic impurities less than 0.0005 $\mu_{\rm B}$/2Fe.
The electronic specific-heat coefficient $\gamma_0$ of YFe$_2$Zn$_{20}$ is 0.05 J/mol K, as shown in the inset.
We analyzed the magnetic properties of YFe$_2$Zn$_{20}$ based on the equation,
\begin{eqnarray}
H = \left(\frac{1}{\chi_0}-I+dT^2+b \langle m^2\rangle \right)M_{\rm Fe}+bM_{\rm Fe}^3,
\end{eqnarray}
where $M_{\rm Fe}$ denotes the magnetization of Fe atoms.
This equation is equivalent to the Stoner model when $\langle m^2\rangle=0$.
The parameter $I$ is a Stoner's enhancement factor,
and $\langle m^2\rangle$ is a correlation function of magnetization of Fe.
Thus, this equation is equivalent to the SCR theory
if this term is correctly taken into account\cite{mor1,mor2,taka}.
In this paper, however, to simplify the analysis, we assume that  $\langle m^2\rangle$ is approximately proportional to the temperature $T$,
that is, $b\langle m^2\rangle=cT$.
This assumption can hold in the paramagnetic region,
because the SCR theory provides the same form as the Curie-Weiss law in this region.
Therefore, experimentally, we can use the equation,
\begin{eqnarray}
H =  (a+cT+dT^2)M_{\rm Fe}+bM_{\rm Fe}^3,
\end{eqnarray}
where $a=1/\chi_0-I$.
In the limit of $H \rightarrow 0, M_{\rm Fe}  \rightarrow 0$, and thus eq. 2
corresponds to the Curie-Weiss law (CW) of the  localized model if $d=0$.
The term of $d$ is a correction term from the straight line of the CW law.
These parameters are experimentally estimated as; $a= 96$ mol/emu, $c= 2.18$ mol/emu K,
$d= 3.7 \times 10 ^{-4}$ mol/emu K$^2$.
The solid line in Fig. 4 is a calculated line using these parameters.
As shown in the inset in Fig. 4,
at low temperatures in finite magnetic fields,
$M(H)$ curve is not a straight line.
The second term in eq. 2 is a correction term to correct the straight line.
Thus, $b$ is experimentally obtained as 9.19$\times 10^{-5}$ Oe/(emu/mol)$^3$.
The solid lines of $M(H)$ curves in the inset of Fig. 4 are the calculated ones using these parameters.
The parameter $d$ is not effective at low temperatures.
The parameter $\chi_0$ is estimated as 6.85$\times 10^{-4}$ emu/mol
using the experimental value of $\gamma_0$.
The Stoner's enhancement factor $\alpha=\chi_0I$ is 0.93.
\begin{table}[b]
\caption{Characteristic parameters.
$a$ [mol/emu], $c$ [mol/emu K], $d$ [10$^{-4}$mol/emu K$^2$],
$b$ [10$^{-5}$Oe/(emu/mol)$^3$],
$\theta_p$ [K], $p_{\rm eff}$ [$\mu_{\rm B}$/Fe], $\alpha$: Stoner's enhancement factor.}
\label{table:1}
\begin{tabular}{ccccccccc}  \hline
sample &  a &    c  &  d    &  b   & $\alpha$ & $\theta_p$ & $p_{\rm eff}$  &\\ \hline
\#2  &    96 &  2.18 &$-$3.7 & 9.19 & 0.93   & $-$44.1    & 1.36  & \\ \hline
\end{tabular}
\end{table}

Next, we define the following hamiltonian to analyze the anomalous properties of DyFe$_2$Zn$_{20}$,
\begin{eqnarray}
H= H_{\rm CEF}+ H_{\rm exch} + H_{\rm Z} +H_{\rm FeR},
\end{eqnarray}
where $H_{\rm CEF}$ is the crystalline-electric-field (CEF) hamiltonain, $H_{\rm Z}$ the Zeeman hamiltonian,
and $H_{\rm exch}$ the exchange hamiltonian between the two Dy atoms.
The former three terms are hamiltonians for Dy atoms in DyFe$_2$Zn$_{20}$.
$H_{\rm CEF}$ is written as $W(\frac{x}{F_4}O_4+\frac{1-|x|}{F_6}O_6)$,
$H_{\rm Z}$ is $-\Mbf_{\rm R}\Hbf_{\rm ext}$,
and $H_{\rm exch} = -\Mbf_{\rm R}\Hbf_{\rm mol}$,
where $\Mbf_{\rm R}$ is the magnetization of Dy atoms, $\Hbf_{\rm ext}$ the external field,
$\Hbf_{\rm mol}$ the molecular field at Dy atom caused by surrounding Dy atoms,
$\Mbf_{\rm R} = -g_J\mu _{\rm B}\Jbf$, $\Hbf_{\rm mol} = n_{\rm RR}\Mbf_{\rm R}$,
and  $n_{\rm RR}$ the exchange constant between two Dy atoms.
The last term in eq. 3 is an exchange interaction between Dy and Fe atoms, which is written as
\begin{eqnarray}
H_{\rm FeR}= -n_{\rm FeR}  \Mbf_{\rm Fe} \Mbf_{\rm R},
\end{eqnarray}
where $\Mbf_{\rm Fe}$ is the magnetization of Fe atoms, and
$n_{\rm FeR}$ the exchange constant between Dy and Fe atoms.
Note that the physical quantities of $\Mbf_{\rm R}$, $\Jbf$, $\Hbf_{\rm mol}$, and $\Hbf_{\rm ext}$ are three-dimensional vectors.
The molecular field at the Dy atom caused by the surrounding Fe atoms is written as
\begin{eqnarray}
n_{\rm FeR}\Mbf_{\rm Fe}=n_{\rm FeR}M_{\rm Fe}
\left(\frac{\Hbf_{\rm ext}+n_{\rm FeR}\Mbf_{\rm R} }{|\Hbf_{\rm ext}+n_{\rm FeR}\Mbf_{\rm R}|}\right),
\end{eqnarray}
where $M_{\rm Fe}=|\Mbf_{\rm Fe}|$, and the term in the parenthesis is a unit vector parallel to $\Mbf_{\rm Fe}$.
In eq. 5, we assumed that
the magnetic moment of Fe is isotropic, i.e., aligning parallel
to the direction of $\Hbf_{\rm ext} + n_{\rm FeR}\Mbf_{\rm R}$,
which is the sum of $\Hbf_{\rm ext}$ and the molecular field at the Fe site.
The direction of $\Hbf_{\rm ext}$ is not always parallel to the direction of $\Mbf_{\rm R}$,
especially in the ferromagnetic temperature region.
Consequently,
$M_{\rm Fe}$ obeys the following equation,
\begin{eqnarray}
|\Hbf_{\rm ext} + n_{\rm FeR}\Mbf_{\rm R}|=  (a+cT+dT^2)M_{\rm Fe}+bM_{\rm Fe}^3.
\end{eqnarray}

We summarize the hamiltonian for Dy atom in the frame of the molecular-field approximation as follows:
\begin{eqnarray}
H &=& H_{\rm CEF}  + g_J\mu_{\rm B}(1+f_{\rm Z}) \Hbf_{\rm ext} \Jbf  \nonumber \\
&& -(g_J\mu_{\rm B})^2 (n_{\rm RR}+f_{\rm exch})\langle \Jbf\rangle\Jbf,
\end{eqnarray}
where
\begin{eqnarray}
f_{\rm Z} &=&  \frac{n_{\rm FeR}M_{\rm Fe}}{|\Hbf_{\rm ext}-g_J\mu_{\rm B}n_{\rm FeR}\langle \Jbf\rangle|},  \\
f_{\rm exch} &=&  \frac{n_{\rm FeR}^2M_{\rm Fe}}{|\Hbf_{\rm ext}-g_J\mu_{\rm B}n_{\rm FeR}\langle \Jbf\rangle|}.
\end{eqnarray}
The magnetization of Dy atom is calculated by
\begin{eqnarray}
\langle \Jbf\rangle= \frac{{\rm Tr}~ \Jbf\exp(-\beta H)}{{\rm Tr}~ \exp(-\beta H)},
\end{eqnarray}
where $\beta = 1/k_{\rm B}T$.
The specific heat $C$ is calculated by
\begin{eqnarray}
C = N_{\rm A}\frac{\partial}{\partial T}\left (\langle H\rangle -\frac{1}{2}\langle H_{\rm exch}\rangle \right),
\end{eqnarray}
where the second term in the parenthesis is a correction term to correct
the double-counting of the exchange energy between the two Dy atoms.

There are eight parameters
to fit the experimental $\Mbf(T)$, $\Mbf(\Hbf)$, $C(T)$, and $C(\Hbf)$ in this calculation:
 $x,~W,~n_{\rm RR},~n_{\rm FeR}, a,~ c,~d$ and $b$.
However, the parameters $a,~c,~d$ and $b$ are fixed as in YFe$_2$Zn$_{20}$.
Moreover, the parameters $x,~W$ and $n_{\rm RR}$ are estimated by adjusting
the magnetic crystalline anisotropy and $T_c$ for  DyRu$_2$Zn$_{20}$,
which is the compound without magnetic Fe atoms, where $T_c$ of DyRu$_2$Zn$_{20}$ is 4 K.
The final parameter is $n_{\rm FeR}$, which is determined by
adjusting $T_c$ of DyFe$_2$Zn$_{20}$.
The eight fitting parameters seem to be easily found.

\par
However, difficulties of the calculation exist
in the iteration procedures
when we solve simultaneously  both $\Mbf_{\rm R}$ of Dy atoms and $\Mbf_{\rm Fe}$ of Fe atoms
using eqs. 6 and 10.
In this paper, we have adopted a two-step procedure to solve these equations.
First, we ignore the anisotropy of $\Mbf_{\rm R}$.
When $H_{\rm CEF}$ is excluded from the total hamiltonian,
the directions of $\Mbf_{\rm Fe}$ and $\Mbf_{\rm R}$ are parallel to $\Hbf_{\rm ext}$ in the case of positive $n_{\rm FeR}$.
$M_{\rm R}$ is easily expressed using the Brillouin function $B_J(x)$ as
\begin{eqnarray}
\langle J_z\rangle = -JB_J(x),
\end{eqnarray}
where
\begin{eqnarray}
x =
\left(\frac{g_J\mu_{\rm B}J(H+n_{\rm RR}\langle J_z\rangle +n_{\rm FeR}M_{\rm Fe})}{k_{\rm B}T}\right).
\end{eqnarray}
In the case of negative $n_{\rm FeR}$,
the directions of $\Mbf_{\rm Fe}$ and $\Mbf_{\rm R}$ are not always parallel to $\Hbf_{\rm ext}$.
However, we used eq. 12 for both cases for simplicity.
Using eqs. 6 and 12, $M_{\rm Fe}$ and $\langle J_z\rangle$ are easily calculated.
As the results of this step,
the contribution of $M_{\rm Fe}$ is found to be  approximately 4 \% of $M_{\rm R}$.
Then, as the second step of the calculation,
the anisotropic $\Jbf$ of Dy atom is re-calculated based on eq. 10
by the iteration procedure.
In this step, we fixed $M_{\rm Fe}$ at the value calculated in the first step.

Figures 5, 6 and 7 show the calculated results of $M_{\rm R}(H)$, $M_{\rm R}(T)$ and $C(T)$, respectively.
The calculated $M$ in these figures is $M_{\rm R}$, because $M_{\rm Fe}$ is presumably smaller than $M_{\rm R}$.
The CEF parameters and the exchange interaction constants are listed in Table 2.
\begin{figure}[t]
\centering
\includegraphics[width=75mm]{cal_MvsH_3}
\caption{Calculated $M(H)$ curves at 2 and 30 K in $H \parallel $ [100], [110], and [111].
Inset shows the calculated $M(H)$ curves at various $T$ in $H \parallel $  [100].}
\label{fig:5}
\end{figure}
\begin{figure}[t]
\centering
\includegraphics[width=75mm]{cal_MvsT}
\caption{Calculated temperature dependence of $M/H$ in $H \parallel $ [100], [110] and [111].
Inset shows the calculated $H/M$ in $H \parallel $ [100]. The thin line in the inset is a guided line for the eyes.}
\label{fig6}
\end{figure}
\begin{figure}[t]
\centering
\includegraphics[width=75mm]{cal_CvsT100}
\caption{Calculated  specific heat $C_{\rm m}(T)$ in various $H \parallel $ [100].}
\label{fig:7}
\end{figure}
\begin{table}[b]
\caption{CEF parameters, $W$ [K] and $x$, and the exchange constants between R-R and R-Fe,
$n_{\rm RR}$ [T/$\mu_{\rm B}$] and $n_{\rm FeR}$ [T/$\mu_{\rm B}$].}
\label{table:2}
\begin{tabular}{ccccc}  \hline
 $W$     & $x$    & $n_{\rm RR}$ & $n_{\rm FeR}$ \\ \hline
$-0.036$ & $-0.2$ & 0.188        & $-14.8$       \\ \hline
\end{tabular}
\end{table}
The energy splitting by CEF is considerably small; the ground state is
$\Gamma_8=-7.5$ K, the first excited state is
$\Gamma_8=-3.5$ K, and the overall energy splitting is 16.3 K.

The calculated results of $M(H)$ and $M(T)$ are in good agreement with the experimental results in Figs. 1 and 2
both qualitatively and quantitatively.
The calculated results represent the following characteristic features of DyFe$_2$Zn$_{20}$:
the calculated anisotropy $M(H)$ at 2 K coincides with the experimental one,
the easy axis of magnetization is parallel to [111],
the metamagnetic transition occurs in $H$ along the [100] direction in approximately 2 T at 2 K,
the metamagnetic transition field decreases with increasing $T$ and goes down to zero near 20 K,
the magnetic anisotropy disappears for all the principal axes above 30 K,
Concerning $C(T)$ in $H$ along the [100] direction, the calculated results reproduce well the experimental features:
$C(T)$ in $H=0$ makes a broad peak at approximately 10 K without any anomalous thorn,
a small anomalous thorn appears by applying $H$ along the [100] direction
and it grows into an apparent peak with increasing $H$ and then disappears above 3 T.
In the cases of $H$ along the [110] and [111] directions,
we could also analyze the experimental results of $M$ and $C$ by our calculation.

We summarize the experimental and calculated results of DyFe$_2$Zn$_{20}$.
DyFe$_2$Zn$_{20}$ is a ferrimagnet with two magnetic moments of Dy and Fe atoms.
The former ion has 10 $\mu_{\rm B}$/Dy, the latter 0.4 $\mu_{\rm _B}$/2Fe.
This small magnetic moment of Fe atom  has been currently inferred
from the M$\ddot{\rm o}$ssbauer experiments\cite{tamu}.
The overall splitting of CEF is 16 K.
The exchange interaction $n_{\rm RR}$ between the two Dy ions is extremely weak compared with the exchange interaction $n_{\rm FeR}$ between the Dy and Fe atoms.
Thus the $T_c$ of DyFe$_2$Zn$_{20}$ is enhanced by this $n_{\rm FeR}$, because
it is 80 times larger in magnitude than $n_{\rm RR}$ in DyFe$_2$Zn$_{20}$.
This small $n_{\rm RR}$ is caused by the distance between Dy ions in DyFe$_2$Zn$_{20}$,
and the large value of $n_{\rm FeR}$ is presumably
due to the large overlapping of the wave function between the itinerant 3d and the localized 4f electrons.

\vskip 5mm
This study is partly supported by a Grant in Aid for Scientific Research (No. 21540356) from the Japan Society for the Promotion of Science.