On the origin of the gamma-ray emission from Omega Centauri: Milisecond pulsars and dark matter annihilation

TL;DR

This paper investigates whether gamma-ray emission from Omega Centauri is due to millisecond pulsars or dark matter annihilation, analyzing observational data and theoretical models to identify the most plausible source.

Contribution

It compares MSP and dark matter scenarios for gamma-ray emission in Omega Centauri, providing constraints on dark matter properties and highlighting the need for further observational data.

Findings

MSP population could explain gamma-ray flux consistent with X-ray observations.

Dark matter annihilation into muons is incompatible with CMB constraints.

A narrow parameter space allows dark matter annihilation into light quarks.

Abstract

We explore two possible scenarios to explain the observed gamma-ray emission associated with the atypical globular cluster Omega-Centauri: emission from millisecond pulsars (MSP) and dark matter (DM) annihilation. In the first case the total number of MSPs needed to produce the gamma-ray flux is compatible with the known (but not confirmed) MSP candidates observed in X-rays. A DM interpretation is motivated by the possibility of Omega-Centauri being the remnant core of an ancient dwarf galaxy hosting a surviving DM component. At least two annihilation channels, light quarks and muons, can plausibly produce the observed gama-ray spectrum. We outline constraints on the parameter space of DM mass versus the product of the pair-annihilation cross section and integrated squared DM density (the so-called J-factor). We translate upper limits on the dark matter content of Omega-Centauri into lower limits on the annihilation cross section. This shows s-wave annihilation into muons to be inconsistent with CMB observations, while a small window for annihilation into light quarks is allowed. Further analysis of Omega-Centauri's internal kinematics, and/or additional information on the resident MSP population will yield much stronger constraints and shed light about the origin of this otherwise mysterious gamma-ray source.

Tables

Table 1. Table 1: Median values obtained by the inference for the MSP model parameters (see figure  3 ) and the corresponding Akaike information criterion parameter (AIC). The error bars for all the parameters are at a 95 % percent 95 95\% confidence.

$𝚪$	${\mathrm{𝐥𝐨𝐠}}_{\mathrm{𝟏𝟎}}\left({𝑬}_{𝒄𝒖𝒕}\right)$	${\mathrm{𝐥𝐨𝐠}}_{\mathrm{𝟏𝟎}}\left({𝑵}_{\mathrm{𝟎}}\right)$	${𝑵}_{𝑴𝑺𝑷}$
${1.58}_{-0.26}^{+0.26}$	${3.55}_{-0.14}^{+0.18}$	$-{6.89}_{-0.85}^{+0.71}$	${31.78}_{-7.97}^{+18.77}$
	${\mathrm{𝐥𝐨𝐠}}_{\mathrm{𝟏𝟎}}\left({𝑵}_{𝑺}\right)$	${𝜶}_{𝒔}$	AIC
	$-{11.43}_{-0.10}^{+0.08}$	${2.70}_{-0.13}^{+0.13}$	20.3

Table 2. Table 2: Median values found by the parameter inference for two different dark matter annihilation channels (see figure  4 ) and the corresponding Akaike information criterion parameter (AIC). The error bars are at a 95 % percent 95 95\% confidence.

	${𝝁}^{+}{𝝁}^{-}$	$\overline{𝒒}𝒒$
${𝒎}_{𝝌}$	${4.30}_{-0.08}^{+0.09}$	${9.10}_{-0.62}^{+0.69}$
${\mathrm{𝐥𝐨𝐠}}_{\mathrm{𝟏𝟎}}\left(⟨𝝈𝒗⟩𝑱\right)$	$-{4.34}_{-0.03}^{+0.03}$	$-{5.50}_{-0.03}^{+0.03}$
AIC	$179.4$	$47.4$

Equations

$$\mathcal{L}=\prod _i\frac{e^{-{\lambda }_i}{\lambda }_i^{d_i}}{d_i!},$$

$$\lambda =M(E,\Theta )+S(E,\eta )+B(E),$$

$$\frac{d\varphi }{dE}=N_0(E/\mathrm{MeV}{)}^{-\Gamma }{\mathrm{e}}^{-E/E_{\mathrm{cut}}},$$

$$M(E_i)={\int }_{E_i}^{E_{i+1}}{\int }_{\Delta \Omega }d\Omega ϵ(E)\mathrm{PSF}(\mathrm{E},\mathrm{\Theta })\frac{\mathrm{d}\varphi }{\mathrm{dE}}.$$

$$N_{\mathrm{MSP}}=\frac{L_{\gamma }}{⟨\dot{E}⟩\eta },$$

$$L_{\gamma }=4\pi S{d}^2.$$

$$\frac{d\varphi }{dE}=\frac{1}{8\pi }{J}_{\Delta \Omega }\frac{⟨\sigma v⟩}{m_{\chi }^2}\frac{dN}{dE},$$

$$\frac{dJ}{d\Omega }(\theta )$$

$$J_{\Delta \Omega }$$

$${\rho }_{\ast }(z,R)=\sum _j^M\frac{M_j}{(2\pi {\sigma }_j^2{)}^{3/2}{q}_j}\exp \left[-\frac{1}{2{\sigma }_{\mathrm{j}}^2}\left({\mathrm{R}}^2+\frac{{\mathrm{z}}^2}{{\mathrm{q}}_{\mathrm{j}}^2}\right)\right],$$

$$\rho (r)=\frac{{\rho }_s}{{\left(\frac{r}{r_s}\right)}^{\gamma }{\left(1+\left(\frac{r}{r_s}\right)\right)}^{(3-\gamma )}},$$

$${\frac{dN}{dE}}_{\mathrm{observed}}=\frac{1}{8\pi }\frac{⟨\sigma v⟩}{m_{\chi }^2}\frac{dN}{dE}{\int }_{\Delta \Omega }d\Omega \left(\frac{dJ}{d\Omega }\ast \mathrm{PSF}\right)(E,\theta ).$$

$$J_{\Delta \Omega }=\frac{{\int }_{\Delta \Omega }d\Omega \left(\frac{dJ}{d\Omega }\ast \mathrm{PSF}\right)(E,\theta )}{{\int }_{\Delta \Omega }d\Omega \mathrm{PSF}(\theta )}$$

$$f_{\mathrm{eff}}\frac{⟨\sigma v⟩}{m_{\chi }}<3.5\times 10^{-28}{\mathrm{GeV}}^{-1}{\mathrm{cm}}^3{\mathrm{s}}^{-1},$$

$$\mathrm{PSF}(E)={\int }_{\Delta \Omega }\mathrm{PSF}(\theta |E)d\Omega ,$$

$$\mathrm{PSF}(E)={\int }_{-\theta }^{\theta }dx{\int }_{-\sqrt{{\theta }^2-x^2}}^{\sqrt{{\theta }^2-x^2}}dy\mathrm{PSF}\left(\sqrt{(D-x{)}^2+y^2}|E\right),$$

$$f(x,y)=\int d{x}^{\prime }d{y}^{\prime }J(x^{\prime },y^{\prime })\mathrm{PSF}(x-x^{\prime },y-y^{\prime }).$$

$$\hat{f}(k)=\left({\int }_0^{\infty }2\pi r\frac{dJ}{d\Omega }(r){ȷ}_0(kr)dr\right)$$

$$\times \left({\int }_0^{\infty }2\pi r\mathrm{PSF}(r,E){ȷ}_0(kr)dr\right),$$

Full text

I Introduction

Omega Centauri ($\omega$-Cen ) is a globular cluster (GC) that has received significant attention over the last few decades due to its unique properties. Arguably the most massive GC in the Galaxy with a mass $\sim 10^{6}$ $\rm{M}_{\odot}$ 2013MNRAS.436.2598W ; Watkins:2015zza ; McLaughlin:2006mp ; 2015ApJ…812..149W , the nature of its stellar population is a rather hotly debated topic. Previous work suggests that $\omega$-Cen ’s star formation history is quite different from that of other GCs, making it challenging to explain its stellar population as an isolated object Bekki:2005rx ; Sollima:2005qy ; Joo:2012gr ; DAntona:2011ary ; 2017MNRAS.471L..31P . $\omega$-Cen has a high concentration of calcium and other heavy metals Lee:2009gg as compared to other GCs, which can be explained by supernovae metal enrichment Johnson:2009ku ; 1995ApJS..101..181W . However, $\omega$-Cen ’s present gravitational potential is not deep enough to retain the ejected material from a supernova explosion Johnson:2009ku . Such observations led to the hypothesis that this unique GC may be the remnant core of an ancient, tidally disrupted dwarf galaxy 1996ApJ…462..241N ; Hilker:2000hr ; Gnedin:2002un ; Carraro:2000ji ; Marcolini:2007kz ; Romano:2007kh ; 2013ApJ…762…36J ; Lee:2009gg ; Bekki:2005rx ; Rey:2003wk ; Fellhauer:2004gc ; Majewski:1999fw ; 1988IAUS..126..603Z ; 2004MNRAS.350.1141T . Unlike a typical GC, dwarf galaxies are dark matter (DM) dominated systems, as inferred from their dynamically-determined mass-to-light ratios Mateo:1998wg . If $\omega$-Cen is the vestige of one of such object it must have lost most of its external DM through tidal interactions with the Milky Way since its stellar kinematics has not yet shown any strong evidence for the presence of a DM halo. The recent discovery of tidal debris spread along $\omega$-Cen ’s orbit adds evidence to this scenario 2019arXiv190209544I .

The analysis of recent measurements of $\omega$-Cen ’s internal kinematics 2015ApJ…812..149W , using the Hubble Space Telescope (HSP) proper-motion catalogue, assuming a constant mass-to-light ratio $\Upsilon$, found $\Upsilon=2.66$. On the other hand, a value of $\Upsilon=1.87\pm 0.15$ was determined from stellar population-synthesis modeling McLaughlin:2006mp . Differences in mass-to-light ratios determined from the kinematics and stellar populations, could be explained due to variations in the initial mass function, or due to mass that can be attributed to the dark matter component, see the discussion in Cappellari:2011 for an example. Even though we do not claim here that the differences found for $\omega$-Cen per se suggests the presence of a relic DM halo, it can certainly be used to derive an upper limit to the amount of DM present on it.

Dark matter is currently one of the most pressing puzzles in physics and astronomy, making up nearly 26% Aghanim:2018eyx of the total energy density in the Universe. With no non-gravitational signals detected so far, its fundamental nature remains a mystery. One distinct possibility is that particle dark matter annihilates with its anti-particle to produce stable Standard Model particles, including photons. The typical energies of such photons would then be around and below the dark matter particle mass. For weak-scale dark matter particles, this would entail the production of gamma rays, a possible tell-tale sign of dark matter annihilation.

In 2010 the Fermi Large Area Telescope (Fermi-LAT) collaboration Abdo:2010nc ; Abdo:2010dk ; Ackermann:2012qk published an analysis of $\gamma$-ray emission from a collection of Milky Way GCs. They reported a non-zero flux from some GCs including $\omega$-Cen collaboration:2010bb . The conventional explanation is that the emission is due to populations of millisecond pulsars (MSPs) within the GCs. These ancient, short-period pulsars are thought to be “recycled” neutron stars “spun up” through the loss of orbital angular momentum to their binary companion. For $\omega$-Cen , the study concluded that a population of $19\pm 9$ MSPs are needed to reproduce the data, assuming all MSPs have the same values for emissivity power and efficiency collaboration:2010bb . Recent X-ray experiments show there are about thirty candidate MSPs in $\omega$-Cen  Henleywillis:2018ufu , but none of these has yet been confirmed as a pulsar.

The search for high-energy particles, in particular $\gamma$-rays, produced by DM annihilation (or decay) in astrophysical systems has been extensively pursued over the last decade (for a review see profumo2017introduction ; Conrad:2017pms ). Such an astronomical detection would shed light on the fundamental nature of DM and particle physics beyond the Standard Model. Given that $\omega$-Cen could be the remnant core of a disrupted galaxy and actually have a (subdominant) DM component, it is timely to explore under which conditions DM annihilation can reproduce Fermi-LAT observations. This can impose important constraints on the parameter space of DM particle properties as well as on the distribution of DM within $\omega$-Cen itself, e.g. Bartels:2018qgr ; Brown:2018pwq .

In this work we analyze $\sim$ 9.5 years of $\gamma$-ray observations by Fermi-LAT from the direction of $\omega$-Cen . Data reduction and analysis methods are outlined in Sec. II; In Sec. III we present the details, and the assumptions we make, for the MSP and DM scenarios for the origin of the $\gamma$-ray emission; in Sec. IV we describe the results of the analysis in the two scenarios of MSP and DM and we also provide an extended discussion on the selection of the annihilation channels and possible constraints on the thermal averaged cross section. Finally, we discuss our results and summarize our findings in Sec. V.

We emphasize that we do not claim here that DM annihilation is the most viable mechanism to explain the observed $\gamma$-ray emission from $\omega$-Cen . Our goal is, rather, to place conservative constraints on DM properties, including its abundance, in order to be compatible with observations. Finally we intend to compare this alternative explanation against a model where MSPs give rise to the $\gamma$-ray emission.

II Data selection and analysis

In this section we explain in detail the procedure used to reduce the data as well as the methodology to analyze it. The main purpose of the Fermi-LAT Abdo:2010nc ; Abdo:2010dk ; Ackermann:2012qk is to observe energetic $\gamma$-rays coming from Galactic and extra-galactic sources whose energy lies in the range from $\sim 100$ MeV up to $\gtrsim 300$ GeV. As stated before, our object of study is the globular cluster Omega Centauri ($\omega$-Cen ), with coordinates: $201.6970$ RA and $-47.4795$  DEC according to the Fourth Fermi Gamma-ray Source Catalog (4FGL) catalog Fermi-LAT:2019yla ($2.89\times 10^{-5}$ deg from $\omega$-Cen ’s position according to 2000A&AS..143….9W ). Using v10r0p5 of the Fermi Science Tools 111https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/, we select $\sim 9.5$ years of LAT data (between MET $239557417$s and $574181612$s). We adopt a region of interest (ROI) of $0.2\deg$ radius around the position of $\omega$-Cen in the 4FGL and select all Pass 8R3 Source events (evclass=128 and evtype=3), in the $0.2-50$ GeV range within the ROI using gtselect with a maximum zenith angle of $90\deg$. We use the standard filter "DATA_QUAL$>0$ && (LAT_CONFIG$==1$)" in gtmktime, we use dcostheta=0.025 and binsz=1 in gtltcube, and theta=300 and thetamax=1 deg in gtpsf to obtain the exposure, $\epsilon$, (observation time times effective area) and point spread function, PSF, both functions of energy. We extract the measured energy spectrum in 15 equally log-spaced energy bins between 0.2 and 50 GeV using gtbin with the PHA1 algorithm.

We analyze the data by fitting proposed models to the observed photon flux, see for instance Abdo:2009jfa ; Hooper:2016rap ; Cholis:2014noa in the context of LAT data. Since the the photon count coming from $\omega$-Cen is low the observations are most accurately described by a Poisson likelihood, defined by

[TABLE]

where the product is over energy bins $i$, $d_{i}$ is the number of counts observed in each bin, and $\lambda_{i}$ is the corresponding model prediction.

There is one known gamma-ray source (4FGL J1328.5-4727) within $1^{\circ}$ of our target according to the 4FGL222In the previous 3FGL catalog Acero:2015hja , there were two nearby sources within 1 degree of $\omega$-Cen . We report the analysis with only the one source present in the latest catalog. However, our conclusions are robust to the inclusion of the other one. (besides the source associated with $\omega$-Cen itself). According to the 4FGL, its energy spectrum corresponds to a power law with no cutoff. To take into account its contribution to the total number of events we will include it in the likelihood at the catalog position and use several nuisance parameters to describe its energy spectrum.

The general form for the model is given by

[TABLE]

where $M$ corresponds to the emission from $\omega$-Cen (to be defined in section III.1 for MSP, and III.2 for DM). $S$ corresponds to the nearby point source for which the flux is modeled with a simple power law with no cut-off: $d\phi/dE=N_{s}(E/E_{p})^{-\alpha_{s}}$. Here $N_{s}$ and $\alpha_{s}$ are the normalization and the spectral index, respectively, and will be treated as nuisance parameters that we marginalize over ($E_{p}$ is an arbitrary reference energy which we set to 4726.70 MeV according to the 4FGL catalog). To obtain the number of events contributed by this point source, we integrate the energy flux within the energy bins limits, as well as integrating the PSF (centered on the point source) over the $\omega$-Cen ROI as specified in Eq. (16).

Finally $B(E)$ corresponds to the diffuse and isotropic background components. We adopt the Fermi collaboration model333http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html consisting of interstellar galactic emission444gll_iem_v07.fits and an isotropic component555iso_P8R3_SOURCE_V2_v1.txt. The prediction $B(E)$ within each energy bin is found by multiplying the sum of these two component spectra by the LAT exposure and the solid angle of the ROI and then integrating over the energy bin.

With the likelihood defined, we perform a Monte Carlo Markov Chain (MCMC) analysis to generate posterior distributions for the model parameters. We use flat priors for parameters describing the emission from $\omega$-Cen , whereas Gaussian priors are used for the those describing the nearby source. These are specified in section IV. For the MCMC we use the emcee package, which is an affine invariant ensemble sampler 2013PASP..125..306F . Our base setup for the MCMC is to use 200 walkers with 1500 steps and a burn-in phase of 500 steps. To check for convergence we performed two checks: (a) we doubled the number of steps with a fixed number of walkers; (b) we doubled the number of walkers with the number of steps fixed as recommended by the authors of the sampler to reduce the noise in the chains. In both tests we find similar results to the base setup. Finally, to compare the different models we use the Akaike Information Criterion, $AIC=-2k-2\ln(\mathcal{Z})$, where $k$ is the number of free parameters in the model and $\mathcal{Z}$ is the evidence. Smaller values of the AIC indicate the model is a better fit to the data, after penalizing for the number of free parameters in the model. In this case evidences were computed using the pymultinest code Feroz:2007kg ; 2014A&A…564A.125B , that uses a nested sampling algorithm, for which one of the main features is to compute the evidence. The use of multinest also served as a cross check of the emcee results.

III Possible explanations for the $\omega$-cen observed $\gamma$-rays events

III.1 Milisecond Pulsars

Thus far, the only hypothesis explored to explain the observed flux of $\gamma$-rays from $\omega$-Cen is that they originate from MSPs collaboration:2010bb ; Cholis:2014noa ; Hooper:2016rap . In this scenario the photon flux is commonly modelled with an exponentially cut-off power law given by

[TABLE]

where $N_{0}$ is the normalization, $\Gamma$ is the spectral index and $E_{\rm cut}$ is the energy cutoff. Those three are the free parameters to be fitted in our the analysis, together with two nuisance parameters associated to the nearby point source. To compute the predicted number of events detected in an energy bin with limits $E_{1}$ and $E_{2}$ we need to include information from the instrument response, namely the exposure $\epsilon$ and the point spread function (PSF) which gives the probability that a photon originating from $\omega$-Cen is reconstructed within our defined ROI. Since these two quantities are energy-dependent, the expected number of events for the $i$-th energy bin is given by

[TABLE]

Note that the above equation takes $\omega$-Cen to be a point source of $\gamma$-rays. This is a valid assumption as the scale radius of the emission ought to be of the order the half-light radius Marcolini:2007kz , which is less than the width of the PSF ($\gtrsim 0.2^{\circ}$) at the energies relevant to our analysis.

Once the model parameters are determined, by means of the MCMC using the likelihood defined in Sec. II, the next step is to compute the associated number of MSPs giving rise to the emission. Following collaboration:2010bb , the number of MSPs is approximated by

[TABLE]

where $L_{\gamma}$ is the isotropic luminosity, $\eta$ is the estimated average spin-down to $\gamma$-ray luminosity conversion efficiency and $\langle\dot{E}\rangle$ is the average spin-down power. We adopt the values of $\langle\dot{E}\rangle$ and $\eta$ from Ref. collaboration:2010bb . The simplest approach is to approximate $L_{\gamma}$ directly from the observed energy flux $S$ and the distance to the source $d$ collaboration:2010bb ,

[TABLE]

Note that Ref. Hooper:2016rap estimates the number of MSPs in a globular cluster by correlating the number of binary systems that can be formed given a stellar density and $\gamma$-ray luminosity, whereas in Ref. Cholis:2014noa a correlation between the X-ray and $\gamma$-ray flux is used in order to compute how many MSPs could contribute to the $\gamma$-ray flux. These studies argue that the naive way to compute the number of MSPs used by collaboration:2010bb can lead to a wrong answer since, for a globular cluster, the $\gamma$-ray emission could be due to one really luminous MSP or the sum of several fainter ones. The details on how the MSPs are contributing to the $\gamma$-ray photon flux could lead to important differences in the analysis, specially when multiple components of the $\gamma$-ray source are considered. See for instance the analysis by Brown:2018pwq and Bartels:2018qgr on the 47 Tuc GC, which we discuss in section V. As mentioned earlier there are $\sim 30$ sources with X-ray-MSP-like emission in $\omega$-Cen Henleywillis:2018ufu , but there is no confirmation with radio observations that these sources are pulsars. Therefore, we consider the simple approach given by Eq. (6) to be sufficient to obtain the order of magnitude of the number of MSPs required to explain $\omega$-Cen ’s observed $\gamma$-ray photon flux.

III.2 Gamma rays from annihilating dark matter

As stated in Sec. I, one of our goals is to explore the possibility that $\omega$-Cen ’s energy spectrum is produced by the annihilation of DM into Standard Model particles. In this case the photon flux from $\omega$-Cen is given by666Notice that Eq. (7) assumes that the dark matter is its own antiparticle; if this is not the case, and there is an equal number of dark matter particles and antiparticles, an additional factor 1/2 would be needed.

[TABLE]

where $\left\langle\sigma v\right\rangle$ is the thermal averaged annihilation cross section, $m_{\chi}$ is the dark matter particle mass, $dN/dE$ is the photon spectrum produced for each annihilation, and $J_{\Delta\Omega}$ is the so-called astrophysical $J$-factor . The $J$-factor encodes information about the DM density distribution within $\omega$-Cen . It is the integral over a solid angle $\Delta\Omega$ of the $J$-profile $dJ/d\Omega$, which determines the spatial morphology of $\gamma$-ray emission from DM annihilation. The $J$-profile is defined as the integral of the square of the dark matter mass density $\rho(r)$ along a line of sight that makes an angle $\theta$ with the direction toward the center of $\omega$-Cen (e.g. Geringer-Sameth:2014yza ):

[TABLE]

One common approach is to use the observed stellar or gas kinematics to constrain the DM density distribution and then use this density profile to determine the $J$-factor for a given target. For nearby dwarf spheroidal galaxies this is usually done by fitting the free parameters of a proposed DM density profile using the Jeans equation (e.g. Strigari:2006rd ; Bonnivard:2015xpq ). For Galactic Center studies, the proposed Milky Way DM density profile is normalized so that it matches the inferred DM density in the solar neighborhood, and is then extrapolated to the very inner regions Abazajian:2012pn ; TheFermi-LAT:2017vmf . In dwarf spheroidals, propagating the uncertainties in the observed kinematics and the influence of different DM and velocity anisotropy profiles have yielded relatively constrained estimates for the $J$-factor  (e.g. 2015MNRAS.451.2524M, ; Geringer-Sameth:2014yza, ; Bonnivard:2015xpq, ).

We are not aware of any previous work pursuing this kind of analysis for the case of $\omega$-Cen , i.e. there are no specific constraints to the DM density profile that could be associated with it. Performing a Jeans (or similar) analysis of $\omega$-Cen kinematics including a DM component is by itself an interesting possibility which, however, is out of the scope of this paper. For this work we use available information on the different estimates of $\omega$-Cen ’s mass-to-light ratio $\Upsilon$ to set an upper limit on the dark matter mass possibly contained within $\omega$-Cen ’s half-light radius, $r_{h}$, and use that to pinpoint a reasonable range for the $J$-factor . We spell this out next.

Using population-synthesis models to compute stellar properties, a mass-to-light ratio of $\Upsilon_{\rm{syn}}=1.87\pm 0.15$ was found in Ref. McLaughlin:2006mp . On the other hand, in Ref. Watkins:2015zza , spherical, non rotating, dynamical models with constant mass-to-light ratio were fitted to the most recent Hubble Space Telescope proper motions 2014ApJ…797..115B , finding $\Upsilon_{\rm{dyn}}=2.66\pm 0.04$. The tension between these two values might be due to differences in the initial mass functions (IMF) Watkins:2015zza . However the possibility that a dark matter component can explain the difference between dynamical and stellar mass-to-light estimates has not been discussed. We do not attempt to demonstrate here that the difference is minimized by including DM as another $\omega-$cen component, but merely use the difference between the two determinations of $\Upsilon$ to set a possible maximum value for the DM contribution to the mass.

Following Ref. Watkins:2013fna , we compute the stellar mass using a cylindrical stellar density profile given by a Gaussian decomposition:

[TABLE]

where $M_{j}$ are the different mass components, $\sigma_{j}$ controls the extent along the major axis, and $q_{j}$ is the projected flattening. We set all parameters as in Watkins:2013fna . Using this decomposition we compute the stellar mass present in $\omega$-Cen for the two $\Upsilon$ values described above. Then we compute the difference between these two values which is approximately $\delta M(r_{h})=1.01\times 10^{5}\,{\rm M}_{\odot}$, where $r_{h}$ is the $\omega$-Cen half light radius equal to 7 pc Watkins:2013fna . We take $\delta M(r_{h})$ as the maximum total mass within $r_{h}$ that could be attributed to a DM halo. Note that according to the simulations of Bekki:2003qw , the dark matter mass that could remain in a $\omega$-Cen like object after a tidal disruption, assuming $\omega$-Cen is the remnant of a dwarf galaxy, would be $M_{DM}\sim 4\times 10^{5}M_{\odot}$. Since to date there are no observational constraints on $\omega$-Cen ’s DM component we can only point out that the total mass limit we adopt is similar to what is expected from predictions based on these numerical simulations.

To model the DM density in Eq. (8), we adopt a generalized Navarro-Frenk-White density profile;

[TABLE]

where $\rho_{s}$ and $r_{s}$ are the characteristic density and scale radius respectively.

Further, we make some assumptions on the properties of the DM halo. A simple picture of $\omega$-Cen being the remnant of a disrupted galaxy would require the scale radius to be larger than the observed luminous half mass radius, otherwise the stellar component would not survive the disruption. On the other hand, one could extrapolate mass-concentration relations such as the one proposed in Ref. Klypin2014:1411.4001v2 , which seems to hold even for small halos Pilipenko2017:1703.06012v2 , and compute the scale radius for a halo mass of $10^{6}M_{\odot}$, which results in a few parsecs. Therefore we can safely assume the scale radius to be larger or of the order of $r_{\rm h}$. For a given $r_{s}$ and $\gamma$ we determine $\rho_{s}$ so that the enclosed mass within the half-light radius is equal to the DM mass limit found from the difference between the values of the mass-to-light-ratios computed above: $M_{\rm DM}(r_{h})=\delta M(r_{h})=1.01\times 10^{5}\,{\rm M}_{\odot}$.

It is important to take into account that, just as in the calculation for the number of events using the MSP model, we must include the information regarding the instrument response. The case of DM is trickier than for MSPs since the angular extent of the dark matter emission profile may be comparable to or broader than the PSF. The way to incorporate the PSF information is through its convolution with the $J$-profile so that the observed number of events is given by

[TABLE]

Over a small patch of sky the integral over solid angle becomes a 2D rectangular integral which can be performed using a zeroth order one-dimensional Hankel transform Geringer-Sameth:2014yza . For this we use the python package mcfit which uses the FFTlog algorithm Hamilton:1999uv . We give the details of this procedure in Appendix B.

Formally, the density profile parameters $r_{h}$ and $\gamma$ should be sampled by the MCMC, but this would involve performing the convolution on every MCMC step, which would be computationally expensive. As a simpler approach, we will use a single value of the $J$-factor for all energy bins which is then multiplied by the energy-dependent PSF. Then we will give an estimate of the error that we might be introducing by this approximation. In Figure 1 we show, for three energy bins, how the convolution of the $J$-profile with the PSF integrated over the ROI (dashed lines), and divided by the integrated PSF, approaches to the unconvolved $J$-factor for high energies (solid lines) and substantially differs for low energies (dashed lines). That is, for the high-energy bins the equality

[TABLE]

holds. The relatively large difference, between the convolved and unconvolved $J$-factor in the low energy bins will lead to different results specifically for the analysis of the $\bar{q}q$ channel. As we discuss in Sec. V.2, the approximation made does not have a great impact on our conclusions. We also note that the analysis of deMenezes:2018ilq prefers a point-source model of $\omega$-Cen to an extended one, which is also consistent with our approximation of using a constant $J$-factor for all energy bins.

IV Analysis and Results

Using the models for $\gamma$-ray emission from $\omega$-Cen described in Sections III.1 and III.2, we now proceed to analyze the data using the statistical methods specified in Section II.

IV.1 MSPs

We model MSP emission in $\omega$-Cen as a power law spectrum with exponential cut-off given by Eq. (3) with three free parameters: $N_{0}$, $\Gamma$, $E_{\rm cut}$. The two nuisance parameters $N_{s}$ and $\alpha_{s}$ describe the energy spectrum of the nearby source. For the $\omega$-Cen parameters we use flat priors over the following ranges:

• •

$0<\Gamma<3$,

• •

$2<\log_{10}(E_{\rm cut}/{\rm MeV})<8$,

• •

$-18<\log_{10}(N_{0}/{\rm MeV})<0$.

We adopt Gaussian priors for $N_{s}$ and $\alpha_{s}$ based on the best fitting values and error bars reported in the 4FGL catalog Fermi-LAT:2019yla :

• •

$\dfrac{N_{s}}{10^{-14}}=3.81\pm 0.39\;{\rm MeV}^{-1}$

• •

$\alpha_{s}=1.98\pm 0.07$

The results for the fit are shown in Figure 3, where we have added a last panel for the derived quantity of interest, the number of MSPs in $\omega$-Cen . The flux from the MSP contribution that corresponds to the median values of the posterior is shown in Figure 2 (black dashed line). Also shown are the fluxes from the nearby source and the background components (green and blue symbols respectively). In Figure 3 the shaded gray region in the $N_{\rm{MSP}}$ panel corresponds to the range of expected number MSPs reported by Ref. collaboration:2010bb . We can see those results more or less agree with our finding, although the distribution peaks at a slightly higher value, and has a long tail towards a high number of MSP.

IV.2 Dark Matter

Under the dark matter annihilation scenario, the shape of the observed energy spectrum is determined by the mass of the DM particle and the model for final state of the annihilation. There exists some theoretical motivation for considering heavy fermion pairs as the dominant annihilation final state, based upon arguments such as helicity suppression (as in the case for Majorana fermion dark matter, which includes the supersymmetric neutralino Jungman:1995df ) or considering scalar mediators, typically featuring Yukawa-like interactions with fermions (see e.g. McDonald:1993ex ; Boucenna:2011hy ; Profumo:2010kp ; Profumo:2013yn ). However, since here we intend to fit for a $\gamma$-ray spectrum whose features indicate quite clearly a low-mass dark matter particle, it would be inappropriate to consider heavy fermions, such as top and bottom quarks, as well as pair-annihilation to gauge or Higgs bosons. Therefore we focus on the following final states: muon-antimuon ($\mu^{+}\mu^{-}$) and light quark-antiquark ($\bar{q}q$) pairs.

To compute the number of events for these specific annihilation channels we use the following spectra.

• •

Muons: $\frac{dN}{dy}=\frac{\alpha}{\pi}\left(\frac{1-(1-y)^{2}}{y}\right)\left(\ln\left[\frac{s(1-y)}{m_{\mu}^{2}}\right]-1\right)$, where $y=E/m_{\chi}$, $\alpha$ is the fine structure constant, $s$ is the center of mass Mandelstam variable and $m_{\mu}$ is the muon mass. This corresponds to the use of the Altarelli-Parisi splitting function (see e.g. Bartels:2017dpb ), which is appropriate as long as the center of mass energy is much larger than the lepton mass, as is definitely the case for the range of masses of interest to us.

• •

Light quarks: $\frac{dN}{dE}=\alpha_{1}\frac{E}{m_{\chi}}\left(\frac{E}{m_{\chi}}\right)^{-3/2}e^{-\alpha_{2}E/m_{\chi}}$, which is a good approximation to simulated $\gamma$-ray spectra, with $\alpha_{1}=0.95$ and $\alpha_{2}=6.5$ Bergstrom:1997fj ; Strigari:2006rd .

Before detailing the MCMC analysis for these annihilation channels we summarize existing constraints. First, for s-wave annihilation, there are strong constraints on DM with particle mass from 1 MeV up to 1 GeV coming from the angular power spectrum of the Cosmic Microwave Background (CMB). An extra injection of energy from DM annihilation can ionize the neutral hydrogen during the Dark Ages, changing the effective free electron fraction left after recombination, and thus affecting the power spectrum. A thorough analysis of CMB data leads to the constraint Aghanim:2018eyx

[TABLE]

where $f_{\rm{eff}}$ is the effective efficiency function which quantifies the total amount of deposited energy in the hydrogen gas. While $f_{\rm{eff}}$ is redshift-dependent, it was shown that the effect on the power spectrum is the same if we consider an effective value at $z\sim 600$ Slatyer:2012yq . We used a value of $f_{\rm{eff}}=0.3$ for both channels, as found in Slatyer:2012yq . From Eq. 14 we obtain an exclusion limit in the mass vs cross-section plane that is independent of the $\omega$-Cen $\gamma$-ray analysis. This corresponds to the solid black line in Fig. 5. We also show for reference the expected thermal averaged DM cross-section as a function of the particle mass that would lead to a correct relic density for a weakly interacting massive particle (WIMP), as computed by Steigman:2012nb .

To proceed further with the dark matter MCMC analysis we need to set priors for the DM mass and cross section. For the DM mass we choose a wide range that encloses the energies where the $\gamma$-ray flux is observed using the aforementioned constraints. On the other hand, whether or not we use an energy-bin dependent $J$-factor or a constant one in our model (Eq. 12), it is clear from Eq. 7 that there is a degeneracy between the $J$-factor and the cross section. When using a constant value for $J$, we can combine the two into a single parameter $\left\langle\sigma v\right\rangle J$. The prior on this parameter combination is not so evident. By choosing a $\log_{10}(J/\mathrm{GeV^{2}cm^{-5}})=21$, which we find to be a plausible maximum value for $\omega$-Cen based on Figure 1, we consider a prior for the cross section in the range $10^{-31}<\left\langle\sigma v\right\rangle/\mathrm{cm^{3}s^{-1}}<10^{-21}$. This would result in a range of $-10<\log_{10}(\left\langle\sigma v\right\rangle J)<0$. For smaller values of the $J$-factor , the limits of the $\left\langle\sigma v\right\rangle J$ product would be highly constrained by the CMB. With this estimate we chose a large enough prior to include a wide parameter space in the $\left\langle\sigma v\right\rangle J$ Vs mass plane:

• •

$0.5<m_{\chi}/{\text{GeV}}<20$,

• •

​ $-20<\log_{10}(\left\langle\sigma v\right\rangle J/\text{GeV}^{2}\text{cm}^{-2}\text{s}^{-1})<0$,

For the nuisance parameters of the nearby source we used the same Gaussian priors as in the MSP case. In Figure 4 we show the results of the MCMC for the $\mu^{+}\mu_{-}$ and $\bar{q}q$ annihilation channels.

In Figure 5 we present the CMB constraints (black solid line), the relic density (black dashed line), and the results of the MCMC analysis (gray vertical shaded region), where the lower limit on $\langle\sigma v\rangle$corresponds the reference value of $\log_{10}J=21$. As we can see, for our choice of the $J$-factor , annihilation into muons is excluded by Planck constraints, in the s-wave approximation, whilst the case of $\bar{q}q$ shows a small allowed window. A better determination of the $J$-factor will help reduce further the allowed parameter space. On the other hand, this could also be interpreted as an upper limit for the $J$-factor for a given $\langle\sigma v\rangle$.

As we claimed in Section III.2, the formal way to proceed in the analysis would be to include the convolution of the PSF with the $J$-factor . This can be done by introducing such computation into the MCMC, leading to an exploration of a six-dimensional parameter space, the thermal averaged cross section $\langle\sigma v\rangle$, the DM density profile slope $\gamma$, the scale radius $r_{s}$, the particle mass $m_{\chi}$, and the two nuisance parameters associated to the nearby source. The computation of the aforementioned convolution is computationally demanding, but the actual limitation is that there is a clear degeneracy between the first three parameters. Therefore, in this work we argue that using a constant $J$-factor multiplied by the PSF, instead of the convolution will not alter our conclusions. As seen on Figure 1, the convolution for the higher energy bins approaches the value the non-convolved $J$-factor . In order to quantify the error introduced by taking a constant $J$-factor we took a value of $r_{s}=7$ pc and $\gamma=0.5$, computed the convolution for all energy bins, and performed the MCMC analysis over the rest of the parameters ($\langle\sigma v\rangle$, $m_{\chi}$, $N_{s}$ and $\alpha_{s}$). We then compared the result with an identical procedure but fixing the $J$-factor to the non-convolved value for this combination of $r_{s}$ and $\gamma$, $log_{10}(J)=19.3$. We found that for $\bar{q}q$ the value of the inferred DM particle mass changes approximately $\sim 17\%$, while in the case of DM annihilating into $\mu^{+}\mu^{-}$ leaves the posterior unchanged. This can be attributed to the fact that the muon spectrum does not contribute strongly to the first energy bins, which have the largest discrepancy once the convolution is taken into account. In both channels the posteriors on $\langle\sigma v\rangle$are unaffected.

V Discussion

The motivation for this work is connected to the following observations:

(i) Given that $\omega$-Cen ’s formation is not entirely understood and that its present-day properties (such as metallicity) differentiate it from other GCs, it is possible that it is the remnant of a dwarf galaxy whose DM halo has been tidally disrupted 1996ApJ…462..241N ; Hilker:2000hr ; Gnedin:2002un ; Carraro:2000ji ; Marcolini:2007kz ; Romano:2007kh ; 2013ApJ…762…36J ; Lee:2009gg ; Bekki:2005rx ; Rey:2003wk ; Fellhauer:2004gc ; Majewski:1999fw ; 1988IAUS..126..603Z ; 2004MNRAS.350.1141T .

(ii) Recent studies on $\omega$-Cen ’s internal kinematics suggest a larger mass-to-light ratio than expected from a stellar population analysis.

(iii) There is a well-identified $\gamma$-ray emission associated to it that so far cannot be definitively identified with known astrophysical objects, MSPs in particular.

Under these premises, we studied $\omega$-Cen by fitting the $\gamma$-ray observed flux using two different models for its origin: unresolved MSPs and DM annihilation.

V.1 MSP

In the MSP scenario the quantity of interest is the number of MSPs inferred from the observed flux. Taking a simple approach, we concluded that a median value of $\sim 30$ MSPs are needed to explain the observed $\gamma$-ray emission, slightly higher than the expected number reported in collaboration:2010bb . This can be understood as current data indicates a higher energy cut off $E_{\rm{cut}}$ and harder spectral index $\Gamma$ than reported in Ref. collaboration:2010bb . Nonetheless, the number of MSPs obtained here is consistent with the previous work of deMenezes:2018ilq who estimated the number using a similar approach to collaboration:2010bb . However, it has been discussed whether the estimate of the MSP population computed through Eq. (4) is reliable Cholis:2014noa ; it could be that the luminosity of $\omega$-Cen comes from a single very luminous pulsar or it could be the contribution of many fainter ones. The fact that there are $\sim 30$ MSP candidates Henleywillis:2018ufu in X-ray points to the plausible explanation of multiple MSP that also produce the $\gamma$-ray emission. With stringent constraints on the number of MSPs a deeper analysis on the luminosity function will be needed.

V.2 Dark Matter

Another possible scenario is that $\omega$-Cen is the remnant of a dwarf galaxy, leading to the possibility that dark matter makes a (subdominant) contribution to $\omega$-Cen ’s mass. In this work we make a first attempt to use the observed $\gamma$-ray flux associated with $\omega$-Cen to constrain the DM particle properties, and overall DM content itself. Given the energy range of the observed flux, we explored the possibility of it being produced by dark matter particles annihilating into $\mu^{+}\mu^{-}$ and $\bar{q}q$. For the analysis we treated the product of the astrophysical $J$-factor and the thermal averaged cross section as a free parameter $J\left\langle\sigma v\right\rangle$ along with the DM particle mass $m_{\chi}$. The constraints on such parameters, found on the MCMC analysis, are reported in Table 2. Figure 5 summarizes our results.

With the posteriors in hand, we compare the analysis results with current DM annihilation constraints. In Figure 5, the light shaded region corresponds to the thermal-averaged cross section values that are excluded by the latest Planck CMB constraints (for s-wave annihilation). The dark gray vertical band corresponds to the allowed region at the 95% confidence level, for the DM mass, found in this work. The lower limit on the cross section is set by the largest $J$-factor we consider possible for $\omega$-Cen . For the case of DM annihilating into $\bar{q}q$, there is an allowed range of $J$-factor values that will be consistent with both constraints on the maximum amount of DM in $\omega$-Cen , and the $\gamma$-ray flux. Conversely, the case of DM annihilating into $\mu^{+}\mu^{-}$, constraints are more stringent. All of the parameter space is already excluded by CMB, meaning that it will require a very small cross-section and/or a very large $J$-factor (which translates into very compact objects) to explain the $\gamma$-ray flux. Note that taking a p-wave annihilation instead on s-wave relaxes CMB constrains due to the fact that the p-wave DM annihilation does not contribute importantly to the thermal history of the Universe until low redshift.

Despite the fact that only for the $\bar{q}q$ case there is a small window of the allowed parameter space, an assumption that the spectra is produced by multiple contributors can in principle change the values of the mass and cross-section needed to explain the data, possibly lowering the cross section needed to reproduce the photon flux. Recent work assessed the possibility that the $\gamma$-ray emission from another globular cluster, 47 Tuc, could be explained by both MSP and DM annihilation. According to Brown:2018pwq the observed photon flux is better fit by a model that includes DM than by a model of MSPs alone, excluding the latter with more than five sigma significance. While Bartels:2018qgr argue there is no significant evidence for a DM component if the variance in the spectrum of a population of MSPs is included in the analysis. This stresses the importance of taking into account the uncertainties in the MSP contribution to the photon flux, specially for mixture models (MSPs and DM). In this work, however, we focused on separate limits for each DM annihilation channel in an effort to gain insight not only on possible explanations to $\omega$-Cen ’s $\gamma$-ray emission, but also to obtain conservative limits on the DM content of $\omega$-Cen under the tidally disrupted dwarf hypothesis.

The results presented were analyzed by taking a $0.2\deg$ ROI. Nonetheless, for both the MSP and DM hypotheses we repeated our analysis using a $0.15\deg$ and a $0.5\deg$ ROI to check for consistency. The posteriors for all three angular regions overlap within $1\sigma$ with the three different data sets.

The predicted $\gamma$-ray flux for each DM model and for the MSP model is shown in Figure 2, based on the values presented in Tables 1 and 2. We see that, qualitatively, the three models are a relatively good fit to the data. A quantitative comparison using the AIC factor, shown in the aforementioned tables, leads to the conclusion that the MSP explanation provides a better fit than the DM models. As in Brown:2018pwq ; Bartels:2018qgr a joint analysis with both MSPs and DM could be performed. After an eventual confirmation of a population of MSPs in $\omega$-Cen , they could be included as an additional background in the DM annihilation scenario. In this case, the DM constraints would become much more stringent. Because of $\omega$-Cen ’s potentially large $J$-factor , the DM constraints may eventually be competitive with those from dwarf spheroidal galaxies. To reach this level we also need a better knowledge of the dark matter density distribution within $\omega$-Cen . This can be explored in future dedicated analysis of $\omega$-Cen ’s stellar kinematics.

To summarize, here we have addressed for the first time the possibility that the $\gamma$-ray emission from $\omega$-Cen can be explained by DM annihilation. Such models are quite constrained. Constraints can also be given in terms of the maximum content of DM in $\omega$-Cen . Further studies of $\omega$-Cen will shed light on the origin of this peculiar system. Dedicated analysis of the $J$-factor will help narrow down the possibilities for DM to explain the observed $\gamma$-ray flux and the confirmation of the existence of MSPs (or other astrophysical sources) will allow stronger constraints on the DM content and properties of $\omega$-Cen .

Acknowledgements

This work was partially funded by a UCMEXUS-CONACYT collaborative project and CONACYT project 28689. AXGM acknowledges support from Cátedras CONACYT. SP is partly supported by the U.S. Department of Energy grant number de-sc0010107.OV acknowledges support from the PAPIIT-UNAM grant IN112518.

Appendix A Integrated PSF for nearby point sources

As stated in Section III, to account for the number of events we must consider the instrument response and in particular the information contained in the point spread function (PSF). The energy-dependent probability of detecting a photon within the ROI is computed by integrating over the ROI,

[TABLE]

where $\theta$ is the angle between the point within the ROI and the $\gamma$-ray source.

For the small angles under consideration we can take the sky as a flat plane and integrate the PSF in rectangular coordinates. The standard integration for $\omega$-Cen is straightforward and has no major complication. However, when we perform the integration for the nearby source, we must integrate its PSF also in a circle centered on $\omega$-Cen , this can be achieved by

[TABLE]

where $\theta$ is radius of the ROI and $D=0.63^{\circ}$ is the angular separation between the source and $\omega$-Cen .

Appendix B Convolution of the Point Spread Function and the $J$-factor

To compute the differential flux we must convolve the $J$-profile and the PSF as in Eq. (12). We make use of the very efficient procedure introduced in Geringer-Sameth:2014qqa . Just as in the last section, we take the source to be in a 2D rectangular plane where the angular separation is approximated by the Euclidean distance. The convolution then has to be done through a 2D integral in the $xy$ plane. We compute the convolution very efficiently by going to Fourier space and taking advantage of the circular symmetry of both the $J$-profile and PSF. Let $(J*{\rm PSF})(x,y)$ be denoted by $f(x,y)$. Then

[TABLE]

The Fourier transform of $f$ is a function only of the magnitude of the $k$-vector $\hat{f}(k_{x},k_{y})=\hat{f}\left(k=\sqrt{k_{x}^{2}+k_{y}^{2}}\right)$ and takes the form

[TABLE]

where $\jmath_{0}$ is the Bessel function of zero order and $dJ/d\Omega$ is defined in Eq. 8 ($r$ being the Euclidean approximation to angle $\theta$). The two terms in (18) (except for a factor of 2$\pi$) are the zeroth-order Hankel transforms, which can be performed efficiently using the FFTlog algorithm Hamilton:1999uv and the python package mcfit. We followed the description of the algorithm described in Ref. Geringer-Sameth:2014qqa to avoid numerical issues when using mcfit and we used the default parameters given by the package functions.