Globular clusters are the oldest gravitationally bounded stellar systems in the Galaxy. Nearly 160 globular clusters are known today (Harris 1996, 2010 edition). These clusters form a spherical halo around the Galaxy, and many of them are located within the Galactic bulge. Owing to the high concentration of stars within globular clusters, these clusters host a large number of compact objects, including neutron stars and White Dwarfs (WDs) many such objects are found in binary systems, forming for example Low-Mass X-ray Binaries (LMXBs) and cataclysmic variables.

Since the 1970s, it has been known that the formation rate per unit mass of LMXBs (Alpar et al. 1982) is several orders of magnitude higher in globular clusters than in the rest of the Galaxy (Clark 1975; Katz 1975). Theoretical arguments have long asserted that the formation of LMXBs is made efficient through frequent stellar encounters. Using the X-ray populations in various globular clusters unveiled by the Chandra X-Ray Observatory, Pooley et al. (2003) and Gendre et al. (2003) found a positive correlation between the number of LMXBs in globular clusters and the stellar encounter rate, Γ_c, putting the dynamical formation scenario of LMXBs in globular clusters on an observational ground. Because MilliSecond Pulsars (MSPs) are generally believed to be descendants of LMXBs, it is quite natural that ~80% of the known MSPs are detected in globular clusters (cf. Manchester et al. 2005). Recent observations of a class of transitional MSPs further support a strong evolutionary link between LMXBs and radio MSPs (see, e.g., Stappers et al. 2014; Takata et al. 2014; Tam et al. 2014).

Using the cumulative luminosity distribution functions of radio MSPs in globular clusters as a probe of the MSP population in the clusters, Hui et al. (2010) found that the number of MSPs in a globular cluster is correlated with its stellar encounter rate, as well as its metallicity. This finding provides an observational evidence of the dynamical origin of MSPs in globular clusters.

Similar to young pulsars, MSPs can be efficient γ-ray emitters (Abdo et al. 2009c). Given the large numbers of MSPs that reside in globular clusters, it has been expected that globular clusters are γ-ray emitters as well, as was first suggested by Chen (1991). In this paper, we review the main results of high-energy γ-ray observations (100 MeV to ~100 GeV) of Galactic globular clusters, which were mainly contributed by data acquired by using the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope. Observations at even higher energies (i.e., >100 GeV) by using Cherenkov telescopes such as H.E.S.S. also complement the above mentioned GeV observations and will be presented at the end of this article. The interested reader is referred to an earlier review that presents a more theoretical point of view on γ-ray emission from globular clusters (Bednarek 2011).

Observations by using the Energetic Gamma-Ray Experiment Telescope (EGRET) aboard the Compton Gamma-ray Observatory in the 1990s failed to detect any globular clusters, including 47 Tucanae that hosts the largest number of MSPs known at that time (Michelson et al. 1994; Manandhar et al. 1996). Recent LAT observations revealed that the upper limit on the EGRET flux for 47 Tucanae is only twofold higher than the measured LAT flux.

The launch of the Fermi Gamma-ray Space Telescope in 2008 has enabled the discovery of γ-ray emission from globular clusters. The first globular cluster for which γ-ray emission was discovered was 47 Tucanae (Abdo et al. 2009b, see Fig. 1). When this discovery was announced, MSPs as a class were being established as γ-ray emitters (Abdo et al. 2009c), and MSPs were believed to be the only stable known γ-ray sources in globular clusters. The observed γ-rays from the direction of 47 Tucanae were thus naturally attributed to the MSPs in the cluster (Abdo et al. 2009b). This assertion was strengthened by the measured cutoff at ~2.5 GeV in the spectrum (Fig. 1) that is close to the average cutoff energy obtained for nearby, γ-ray emitting MSPs in the Galactic field (Abdo et al. 2009c). From the mean γ-ray luminosity of individual MSPs in the Galactic field, i.e., ~5×10³³ erg s^-1 (Abdo et al. 2009c), an average MSP would not be detected at distances of several kiloparsecs. It is therefore generally believed that ã-rays from globular clusters do not come from a single MSP, but represent collective emission that originates from the cluster population of MSPs as a whole. That being said, the very luminous PSR J1823−3021A and PSR B1821-24 are exceptional (Freire et al. 2011; Wu et al. 2013). Please also refer to Section 2.4 for more details. Assuming that the average γ-ray luminosity of MSPs in 47 Tucanae is similar to that of nearby MSPs, the number of MSPs that yield the observed flux from 47 Tucanae (~2.5×10^-11) erg cm^-2 s^-1; Abdo et al. 2010b) is ~50, nearly twofold higher than 23 MSPs in 47 Tucanae detected in radio surveys . Such discrepancy may stem from incomplete radio surveys of MSPs, e.g., owing to the viewing angle effect.

The second globular cluster that was found to emit γ-rays was Terzan 5 (Kong et al. 2010, see Fig. 1). Terzan 5 hosts the largest number of detected MSPs up to this date (i.e., 33). Its high-energy spectrum also exhibits a cutoff at ~3 GeV (Abdo et al. 2010b; Kong et al. 2010; see Fig. 1). The observed energy flux above 100 MeV is ~7.1×10^-11 erg cm^-2 s^-1 (Abdo et al. 2010b).

Taking the distance estimates of 47 Tucanae and Terzan 5 to be 4.0 kpc (McLaughlin et al. 2006) and 5.5 kpc (Ortolani et al. 2007), respectively, the γ-ray luminosity of Terzan 5 is roughly fivefold that of 47 Tucanae, indicating that either the actual number of MSPs is higher in Terzan 5 compared with 47 Tucanae, or the average γ-rays efficiency of the underlying MSPs in Terzan 5 is higher than that in 47 Tucanae, or a combination of both. There are some indications from cumulative radio luminosity distributions that Terzan 5 may indeed contain more MSPs than 47 Tucanae (Hui et al. 2010). If the number of MSPs that reside in Terzan 5 is smaller than ~5 times that of 47 Tucanae, environmental factors such as the background soft photons (whose number is much higher in Terzan 5 than in 47 Tucanae), as suggested by Bednarek & Sitarek (2007) and Cheng et al. (2010), may explain the higher γ-ray luminosity of Terzan 5.

After these two initial discoveries, the number of γ-ray emitting globular clusters has continued to grow (e.g., Abdo et al. 2010b; Tam et al. 2011; Zhou et al. 2015). These new findings have established globular clusters as a population of γ-ray sources, and in turn allow one to perform population studies (see Section 3). The Fermi LAT 4-year Point Source Catalog (3FGL) includes 15 sources that are categorized as globular clusters, which do not include two globular clusters for which pulsed γ-rays were observed (see Section 2.4). These 15 sources are listed in Table 1.

Attempts have been made to detect pulsed γ-rays from individual MSPs in globular clusters. In the first two years after the launch of the Fermi Gamma-ray Space Telescope, no γ-ray pulsation was found by using available ephemerides for MSPs in various globular clusters (Abdo et al. 2009b, 2010b).

The first successful effort to detect γ-ray pulsation from an individual MSP in a globular cluster was made by Freire et al. (2011). They found a significant, 7-σ detection of γ-ray pulsations above 100 MeV from PSR J1823-3021A in the globular cluster NGC 6624 (Fig. 2). The high measured flux in the pulsed component significantly constrained the total number of MSPs in NGC 6624 to be smaller than 32.

Another γ-ray pulsation was detected from the MSP PSR B1821-24 that resides in the globular clusters M28 (Wu et al. 2013 Fig. 3). It has a spin-down power of   erg s^-1, which makes it the most energetic MSP detected so far. The significance of pulsed γ-ray emission was reported to be ~4.3σ at energies > 200 MeV (Wu et al. 2013). The pulse profile is shown in Fig. 3. The γ-ray pulsation was finally established by Johnson et al. (2013). A source slightly offset from the M28 core (detected when using off-pulse data of PSR B1821-24) was reported, but its association with M28 remains uncertain (Johnson et al. 2013).

There have been attempts to relate the γ-ray luminosity of a globular cluster to other properties of the cluster. Such correlations shed light on the origin of γ-rays emitted from globular clusters.

The two-body stellar encounter rate has long been believed to relate to the binary formation rate and hence the number of MSPs in a globular cluster (see, e.g., Pooley et al. 2003; Hui et al. 2010). The first attempt to correlate γ-ray luminosity with the stellar encounter rate was performed by Abdo et al. (2010b). Using the first 8 known γ-ray globular clusters and the expression where ρ₀ is the central luminosity density and r_c is the core radius (Verbunt & Hut 1987), they found that the γ-ray luminosity of a globular cluster is correlated with the stellar encounter rate, with a linear correlation coefficient of 0.7. Assuming that each MSP emits a similar amount of γ-rays, the number of MSPs in these γ-ray emitting globular clusters is said to be consistent with the number obtained by other estimators (Abdo et al. 2010b).

Additional parameters besides the stellar encounter rate are thought to relate to the number of MSPs in a globular cluster. Hui et al. (2010) identified the metallicity of a globular cluster as another important indicator of the number of MSPs. A binary system in a globular cluster with higher metallicity is more efficient in orbital shrinkage by magnetic braking. This yields a higher likelihood of a successful Rochelobe overflow (Ivanova 2006) and in turn leads to a higher formation rate of MSPs. Hui et al. (2010) found a positive correlation between metallicity and the population of MSPs. Using the available γ-ray data from 15 globular clusters, Hui et al. (2011) found that the γ-ray luminosity of a globular cluster is correlated with the cluster metallicity, [Fe/H], at 99.9% level, as logL_γ= (0.59±0.15)[Fe/H]+(35.56±0.15).

On the other hand, primordially formed binaries are not related to the stellar encounter rate. If they form the majority of binaries, one would expect the binary population to be correlated with the cluster mass parameter, M_GC. Assuming a constant mass-to-light ratio, M_GC can be estimated from the absolute visual magnitude M_V according to M_GC=10^-4M_v) (e.g., Hui et al. 2010). This estimate of the cluster mass is advantageous because it is not correlated with the encounter rate, thus serving as an independent parameter in correlational studies. However, no correlation of M_V with the γ-ray luminosity was found (Hui et al. 2011).

Apart from the above three cluster parameters, the energy densities of optical (u_opt) and infrared (u_IR) photons at the cluster location are also important for predicting the γ-ray properties, because they serve as the seed photons in IC models (Cheng et al. 2010). Using the GALPROP code (Strong & Moskalenko 1998) to estimate the soft photon energy densities, the γ-ray luminosity was found to correlate with the optical and infrared photon densities, respectively, at >96% confidence level, as as well as (Hui et al. 2011).

With the updated sample of γ-ray globular clusters compiled in the 3FGL catalog (Acero et al. 2015; cf. Table 1), we have re-examined the above correlations by using the revised catalog of globular clusters (Harris 1996, 2010 edition). The corresponding soft photon densities were also estimated by using GALPROP. NGC 6624 and M28 were excluded from all these revised statistical analyses because their γ-ray emissions are dominated by the individual energetic pulsars (see Section 2.4). The results are summarized in Table 2. We found that the correlations remained robust. As indicated by non-parametric analysis (i.e., Spearman rank correlation coefficient), which is more robust than the analysis of linear correlations, the correlations of logL_γ - [Fe/H] and log L_γ - log u_IR are the strongest among the tested parameters. The best fits for all log L_γ data are displayed in Fig. 4.

As mentioned above, the total energy output of γ-rays indeed scales with two factors that are crucial to the formation of MSPs: the stellar encounter rate and metallicity. This suggests an intimate relationship between the observed γ-ray emission from globular clusters and the LMXB population, which are precedents of the population of MSPs. The new findings that optical and infrared photon energy densities are also correlated with L_γ conform with the IC model (Cheng et al. 2010), in which both the number of MSPs in a globular cluster and the soft photon energy densities play important roles in the resulting γ-ray output. This prompted Hui et al. (2011) to combine two of the above cluster parameters and to perform a 2-dimensional regression analysis.

In view of the updated sample of γ-ray globular clusters and their physical properties (see Section 3.1), we have revisited these fundamental planes relationships, and the results are given as follows:

The edge-on views of the best-fit fundamental-plane relationships are shown in Fig. 5. These relationships should be taken into account in any realistic model of γ-ray emission from globular clusters.

Two classes of models of γ-ray emission from globular clusters, both attributed to MSPs, have been discussed in the literature. In pulsar magnetosphere models, γ-ray emission is produced via curvature radiation within the magnetospheres of MSPs in the clusters (e.g., Venter & de Jager 2008; Venter et al. 2009). The pulsed emission of a number of MSPs is superposed at different frequencies such that unless a small number of luminous MSPs dominate the combined γ-ray emission, as in the cases for NGC 6624 and M28, pulsations from individual pulsars are difficult to detect.

In another class of models that involve IC emission, electrons are accelerated close to the MSPs or (re-) accelerated in colliding wind shocks and scatter off the optical, infrared, or cosmic microwave background photons (Bednarek & Sitarek 2007; Cheng et al. 2010), giving rise to the observed GeV emission. In this scenario, γ-rays are intrinsically unpulsed.

The γ-ray spectrum and especially the cutoff energy as observed for a number of globular clusters are generally similar to nearby MSPs in the Galactic field. Such coincidence may suggest that the observed γ-rays represent curvature radiation from pulsar magnetospheres. On the other hand, averaged spectra may not be good discriminators for testing various models, because IC models can explain the GeV spectra equally well (Cheng et al. 2010). Moreover, according to Cheng & Taam (2003), MSPs in globular clusters may be very different from those in the Galactic field based on their radio and X-ray properties. In fact, the former may possess complicated magnetic fields owing to frequent stellar encounters (Cheng & Taam 2003), thereby strongly affecting the polar/outer gap structures and in particular quenching the outer gap. This scenario is supported by the fact that the majority of the MSPs in 47 Tucanae are thermal X-ray emitters (Bogdanov et al. 2006).

The differences between the emission signatures of the two models should provide a diagnostic. In the pulsar magnetosphere models, γ-rays with energies of up to a few gigaelectronvolts are emitted, while IC processes may give rise to γ-rays with higher energies, reaching teraelectronvolts, depending on the energy of seed photons. A candidate globular cluster reported to emit >10 GeV emission is Liller 1 (Tam et al. 2011). However, this globular cluster does not exist in, e.g., the 3FGL catalog. In addition, no globular clusters have been detected in the 1FHL (Ackermann et al. 2013) and 2FHL catalogs, which include γ-ray sources detected by the LAT above 10 GeV and 50 GeV, respectively, thus constraining any strong emission above ~10 GeV as predicted by some IC models. Very High Energy (VHE) observations by using Cherenkov telescopes provide further constraints on models at >100 GeV energies (see Section 5).

The size of the emission region also differs between the two classes of models. The γ-ray emission region predicted by the pulsar magnetosphere models is more compact and that predicted by the IC models is more extended. The observed γ-ray emission region is consistent with a point source for almost all globular clusters and shows no sign of extended emission (e.g., Abdo et al. 2010b). While one would expect to see point sources in γ-rays in the pulsar magnetosphere models, the IC models predict a certain degree of extension. For example, Cheng et al. (2010) suggest that the IC emission size should be >10 pc, which in the case of 47 Tucanae corresponds to ≥8' (see also their Fig. 4), larger than the upper limit (i.e., 4.8' ) reported by Abdo et al. (2010b). Therefore, unless the actual diffusion coefficients are lower than those assumed in their calculations, the model of Cheng et al. (2010) may not be able to explain the compactness of the emission size. With the increased angular resolution made available by the Pass 8 data and better photon statistics, it is expected that better constraints will be provided on the emission size of the γ-ray emission.

Clapson et al. (2011) present radio data acquired by using the Effelsberg 100-m telescope in the vicinity of Terzan 5, which may impose some constraints on the IC models. The measured radio flux at 11 cm in the circle of radius 0.15° around the core of Terzan 5 (region 1 in their Fig. 1, reproduced here as Fig. 6) is (0.14±0.21) Jy (local background emission contributes about 30–40% of this flux). Assuming a magnetic field strength of 10⁶ Gauss in the cluster, and according to Eqs. (25)–(27) of Cheng et al. (2010) that assume that the same populations of electrons lose their energy by synchrotron and IC radiation, the modeled radio fluxes at 11 cm are ~10 Jy and ~6.8 Jy, within a few arc-minutes, when the seed photons in their IC model are infrared and optical photons, respectively. It may still be possible that low-energy electrons that are responsible for the synchrotron radiation below the peak at ~44 GHz diffuse further out into region 11, thereby explaining the enhanced radio emission in region 11, for which the flux is (3.86±0.34) Jy, while high-energy electrons with a much shorter diffusion length stay close to the core, up-scatter ambient photons, and give rise to the GeV emission from Terzan 5. More detailed modeling is needed to better understand the astrophysical conditions in that cluster.

Besides originating from MSPs, γ-ray emission from particles accelerated by non-accreting WDs was also suggested (Bednarek 2012). Owing to a large number of WDs, these particles may also contribute to the observed γ-rays from globular clusters. In this scenario, particles are accelerated in the inner magnetosphere of WDs, diffuse out, and produce γ-rays via IC scattering off various soft photon fields, e.g., star light, infrared radiation, and cosmic microwave background radiation.

In summary, current γ-ray and low-frequency data are discordant with some IC models, while cumulative pulsed emission from individual MSPs with fewer modeling uncertainties remains a more favorable model for MeV–GeV emission from globular clusters. That being said, the latter model does not naturally predict any correlation of γ-ray luminosity with optical and/or infrared light at the observed cluster positions (see Section 3). Whether or not both mechanisms are at work in some γ-ray-emitting globular clusters remains an open question. Detecting more MSPs in globular clusters with future radio arrays such as SKA may help to better understand the full population of MSPs, while continuous surveys by various γ-ray detectors are equally crucial, by further increasing the sample of globular clusters that emit γ-rays (for correlational studies), and by deeper observations of individual globular clusters (to constrain the γ-ray spectrum and the emission size).

Gamma-ray emission above 100 GeV can be produced via IC scattering off various soft photon fields by relativistic leptons accelerated in the pulsar magnetosphere (Venter et al. 2009), or leptons re-accelerated in shocks generated by pulsar winds (Bednarek & Sitarek 2007). The expected spectrum of >100 GeV emission is uncertain, as the spectrum of relativistic leptons and the nature of target photons remain uncertain in this scenario. Nevertheless, observations of nearly 20 globular clusters by using several γ-ray instruments (H.E.S.S., MAGIC, and VERITAS) above 100 GeV have been performed, but no detection has been reported yet (Kabuki et al. 2007; Aharonian et al. 2009; Anderhub et al. 2009; McCutcheon et al. 2009; Abramowski et al. 2013).

Terzan 5 is truly exceptional in this regard. The Abramowski et al. (2011) announced the 7-σ detection of >0.4 TeV emission in the vicinity of Terzan 5 (Fig. 7). They estimated the probability of chance coincidence of this new VHE source, HESS J1747−248, with an unknown active Galactic nucleus or pulsar wind nebula to be 10^-4. The reported photon flux of 1.2×10^-12 cm^-2 s^-1 in the 0.44–24 TeV range corresponds to 1.5% of the flux of the Crab nebula in this energy range. The spectrum can be fit by a simple power law of photon index Γ_γ = 2.5±0.3_stat±0.2_sys. Unlike the emission in the 0.1–10 GeV range, the TeV emission region is offset from the cluster core (at the 2-σ level) and shows sign of extended feature.

Bednarek & Sobczak (2014) discuss a scenario in which a bow shock forms ahead of where Terzan 5 is moving, and accelerated particles within the globular cluster escape preferentially in the direction opposite to that of the Terzan 5 motion, generating TeV γ-rays through IC processes. Another possible leptonic scenario is that the HESS J1747−248 originates from IC emission by electrons accelerated in colliding shocks between the collective pulsar wind from the cluster and the Galactic wind, which might also explain the fact that the TeV emission site is on the way between the cluster core and the Galactic center.

In leptonic scenarios, as noted by the Abramowski et al. (2011) and by Domainko (2011), there should be a Klein-Nishina cutoff at a few teraelectronvolts, which has not been observed. However, the lack of a cut-off may simply result from low photon statistics at these high energies, which can be tested by performing deeper TeV observations, e.g., by current IACTs or CTA.

The Abramowski et al. (2011) further discusses two hadronic scenarios: cosmic-rays accelerated by a past supernova or in a short GRB remnant. They have the advantage of being able to explain the observed simple power law spectrum up to 20 TeV. In the supernova interpretation, they argue that given the lack of molecular clouds and thus the low interstellar medium density at this location (n≈0.1 cm^-3 is assumed), the energy of cosmic rays needed to produce the observed TeV flux reaches 10⁵¹ erg, which is rather high for a supernova. The more exotic short GRB scenario was further developed by Domainko (2011), but a similar energy was required and it remains unclear why this particular short GRB can be so energetic (i.e., at least 10⁵¹ erg, which is at the high end of isotropic-equivalent energy of short GRBs; Nakar 2007) and efficient in transferring the energy to cosmic-ray particles. In this scenario, thermal X-rays at the flux level of 10^-12 erg cm^-2 s^-1 are produced in shocks caused by sub-relativistic ejecta expelled during the merger event that heat the interstellar medium (Domainko & Ruffert 2008). Future X–ray observations can help to discriminate various scenarios.

In summary, the origin of HESS J1747−248 remains unclear. The simplest leptonic model for Terzan 5, namely electrons IC up-scattering off longer-wavelength emission, is strongly inconsistent with upper limits derived for 15 other globular clusters if one assumes a similar leptonic model for all these globular clusters (Abramowski et al. 2013). Future X-ray observations at the center of the TeV emission, as well as deeper TeV observations, will be crucial for tackling the origin of the TeV emission.