High Energy Radiation from Spider Pulsars

The population of millisecond pulsars (MSPs) has been expanded considerably in the last decade. Not only is their number increasing, but also various classes of them have been revealed. Among different classes of MSPs, the behaviour of black widows and redbacks are particularly interesting. These systems consist of an MSP and a low-mass companion star in compact binaries with an orbital period of less than a day. In this article, we give an overview of the high energy nature of these two classes of MSPs. Updated catalogues of black widows and redbacks are presented and their X ray/$\gamma$-ray properties are reviewed. Besides the overview, using the most updated eight-year Fermi Large Area Telescope point source catalog, we have compared the $\gamma$-ray properties of these two MSP classes. The results suggest that the X-rays and $\gamma$-rays observed from these MSPs originate from different mechanisms. Lastly, we will also mention the future prospects of studying these spider pulsars with the novel methodologies as well as upcoming observing facilities.


What Are Millisecond Pulsars?
Rotation-powered pulsars, which act as unipolar inductors by coupling the strong magnetic field and fast rotation, radiate at the expense of their rotational energy. As a result, the rotation of a pulsar gradually slows down as it ages. When the rotation becomes too slow to sustain the particles' acceleration in their magnetospheres, the radiation beam shuts down. In such a case, we say a pulsar is "dead". While their magnetospheric particle accelerators have been turned off, other emission mechanisms (e.g., accretion) can still possibly work in these dead pulsars.
By the time of writing, there are 2702 pulsars in total, including radio pulsars, radio-quiet γ-ray pulsars and magnetars [1]. In Figure 1, we show the distribution of their period P and period derivativeṖ. Applying a clustering analysis in this 2D parameter space with k-means partitioning, the whole population can be divided into two parts. The largest one is displayed as black dots in Figure 1. This partition spans a range of P ∼0.02 − 23.5 s andṖ ∼ 5 × 10 −18 − 5 × 10 −10 s/s. This group includes canonical pulsars as well as magnetars. Assuming the surface magnetic field is dipolar and the rotational energy of the pulsar goes entirely to the dipolar radiation, one is able to estimate their surface magnetic field strength B s and their characteristic age τ as: and where the radius R and the moment of inertia I of the neutron star are taken to be 10 6 cm and 10 45 g cm 2 respectively.
These imply that the surface field strength and the characteristic age of this group are at the order of B s ∼ 10 10 − 10 15 G and τ ∼ 10 3 − 10 8 years (see Figure 1). As a canonical pulsar ages and spins down, it moves toward the lower right in the P −Ṗ diagram. One should note that there is an absence of systems at the lower right region in the diagram. Such a region is a so-called "graveyard" zone. The dashed line on Figure 1 is known as the death line and separates the neutron stars that can sustain particle acceleration in their inner magnetospheres from those that cannot [2]. On the other hand, there is a dispersed group at the upper right region of the P −Ṗ diagram. These objects with surface dipolar fields > 10 14 G are magnetars.
The other partition is clustered at the lower left corner of the P −Ṗ diagram, which spans a range of P ∼ 0.001 − 0.5 s andṖ ∼ 5 × 10 −22 − 10 −17 s/s. The spin parameters imply these pulsars have a weak magnetic field at the surface B s ∼ 10 7 − 10 11 G and as old as τ ∼ 10 9 − 10 10 years. These are commonly referred to as millisecond pulsars (MSPs). They are displayed as the grey dots in Figure 1.
It is a consensus that MSPs originate from evolved compact binaries. The standard scenario for the formation of MSP is that a dead pulsar is rejuvenated through accreting material from its binary companion [3]. During this "recycling" phase, the system appears as a low-mass X-ray binary (LMXB). While the details of the magnetic field decay is still a subject under discussion (e.g., [4]), the surface field of the neutron star can be somehow buried by the accretion. Approaching the end of this phase (on a Gyr timescale), the neutron star will gain sufficient angular momentum and resuscitate particle acceleration in the magnetosphere and hence the pulsed radiation again. This marks the birth of a rotation-powered MSP. This scenario was proposed shortly after the discovery of the first MSP PSR B1937+21 [5]. The fast rotations, weak surface magnetic fields and old ages of MSPs are all found to be consistent with this recycling scenario [3]. Period-period derivative (P-Ṗ) diagram of all currently known pulsars. The population can be divided into two groups (black and grey dots) based on k-means partitioning. The locations of black widow and redback millisecond pulsars (MSPs) in this parameter space are highlighted by the green and blue circles, respectively. Lines of constant dipolar surface magnetic field (Equation (1)) and characteristic age (Equation (2)) are shown. The dashed line illustrates the death line for radio pulsars by assuming a multipolar magnetic field configuration [2].

Fermi Gamma-Ray Space Telescope-A Game Changer
Since the discovery of the first pulsar in 1967 [6], radio observations have been the major drive for the progress in pulsar astronomy for a long time. For hunting pulsars, a number of extensive radio surveys have been carries out. For example, Parkes multibeam pulsar survey discovered ∼800 pulsars in 1997-2004 ( [7][8][9]). Once the sky positions of pulsars have been pin-pointed by the radio telescopes, multiwavelength investigations of their nature can be carried out. However, a survey of pulsars simply based on radio observations has certain limitations. For example, the ground-based radio observatories cannot provide a full sky coverage and a high cadence monitoring of their behaviour. compiled the updated catalogues of all known spider MSPs in the galactic field and globular clusters, separately. Currently, there are 44 confirmed black widows systems that have been detected with 27 in the galactic field and 17 reside in globular clusters. For redbacks, there are 26 in total with 14 systems confirmed in the the galactic field and 12 in globular clusters. The distributions of black widows and redbacks in the P −Ṗ diagram are illustrated in Figure 1. In Figure 2, we show their spatial distribution in our galaxy. In the last decade, there have been many interesting high energy phenomena that have been observed from these spider pulsars. These systems contain low-mass companion stars in tight orbits. The interactions between the MSPs and their companions results in many interesting high energy phenomena. In the following sections, we will give an overview of their emission properties in X-ray and γ-ray. (1) O = Significant optical orbital modulation; o = possible optical orbital modulation; X = significant X-ray orbital modulation; x = possible X-ray orbital modulation; g = possible gamma-ray orbital modulation. (2) R = radio pulsations; G = gamma-ray pulsations; X = X-ray pulsations; O = optical pulsations. (1) O = Significant optical orbital modulation; o = possible optical orbital modulation; X = significant X-ray orbital modulation; x = possible X-ray orbital modulation; g = possible gamma-ray orbital modulation.

Black Widows
While the recycling scenario has been proposed for more than 30 yrs, many details of the formation process of MSPs remain uncertain. As mentioned before, an MSP is expected to be the end product of the evolution of a compact binary. One may be puzzled to notice that ∼30% of the known MSPs in the galactic field are found to be isolated [96]. One proposed possibility to explain their solitude is that the high energy emission from these rejuvenated pulsars have ablated their companions [97]; and the companion stars will eventually evaporate entirely.
This proposed scenario was motivated by the discovery of PSR B1957+20 ( = PSR J1959+2048) [61]. It is a binary that contains a 1.6 ms MSP and a very low mass companion (M c ∼ 0.02M ) in a 9.2 h orbit. Eclipses of radio pulsations have been observed from this system during the phases called inferior conjunction when its companion lies between the MSP and us. Eclipses result from the absorption or scattering of the radio signals from the pulsar by the dense ionized gas streaming off from its evaporating companion ( [98][99][100]). Such system is dubbed black widows because the situation is similar to a female spider that devours its mate after mating.
While the pulsed emission is the defining characteristic of a pulsar, it only consumes a tiny fraction of the rotational energy of the neutron star. Even with the pulsed radiation across the whole electromagnetic spectrum summed up, it can hardly exceed a few tens of percents of a pulsar's spin-down power. Most of the rotational energy of the neutron stars is carried away by the relativistic pulsar wind outflows. While the pulsed emission originates within the magnetosphere, the wind region extends beyond the light cylinder. We can detect the presence of the wind through its interactions with surroundings. For the case of PSR B1957+20, the effect of its pulsar wind can be revealed in several ways. PSR B1957+20 has a tranverse velocity of v ∼ 220 km/s [1], which is found to be supersonic and a bow-shock nebula can be formed. The termination shock radius is determined by the balance between the wind particles and the interstellar medium at the head of the shock. Through the interactions with the shocked medium, the relativistic wind particles radiate synchrotron emissions as they trail behind the pulsar's motion. This results in a cometary-tail bow-shock morphology. Figure 3 shows the X-ray bow-shock nebula of PSR B1957+20 as observed by Chandra, which was firstly discovered by [101].
Systems of this kind can enable us to study the shock physics in the interstellar medium (ISM) [102]. The observed length l of the X-ray tail can be interpreted as the distance traversed by the pulsar within the electron synchrotron cooling timescale t c , i.e., l ∼ v p t c , where v p is the proper-motion velocity of the pulsar. The cooling time in the X-ray band is ∼10 8 B −3/2 mG (hν X /keV) −1/2 s, where B mG is the inferred magnetic field strength in the emitting region [102]. Adopting the inferred tail length of ∼9.7 × 10 17 cm at the distance of 2.5 kpc and a proper motion velocity of 220 km s −1 , the synchrotron cooling timescale is estimated to be ∼4.9 × 10 4 yrs. This yields a magnetic field of B ∼ 17.7 µG in the shock region. Considering a magnetic field strength of ∼ 2 − 6 µG in the ISM, the compression factor of the magnetic field in the termination shock is estimated to be ∼3 [102].
Spatially-resolved spectroscopy can also shed light on the properties of the nebula [63]. For a synchrotron nebula powered by the MSP, a softening of the spectrum of the X-ray tail as a function of the distance from the pulsar is expected if synchrotron cooling of the particles injected at the termination shock is dominated. By dividing the tail into two segments, the photon index of the segment closer to the pulsar is found to be Γ ∼ 1.6 and that of the other is Γ ∼ 2.1 [63], which is consistent with the aforementioned scenario.
Currently, there are >70 pulsar wind nebulae detected in the X-ray [103]. However, only three of them are confirmed to be associated with MSPs (PSRs B1957+20 [101], J2124-3358 [104] and J1911-1114 [105]). We therefore encourage dedicated search into the archival X-ray imaging data, which can possibly reveal unreported faint X-ray nebulae associated with MSPs. . X-ray image of the bow-shock nebula associated with black widow MSP PSR B1957+20 obtained from the public data with an effective exposure of ∼165 ks as acquired by Chandra (Observation ID: 9088). This is the same data as used in the study by Huang et al. [63].
The interactions between the pulsar wind and the ablated material from its companion can also form intrabinary shock. Because of the geometry of the shock, this can lead to the modulation of high energy emissions across the orbit. In Figure 4, we show the X-ray light curve of PSR B1957+20 folded at the orbital period [63]. The observed flux is found to peak just before and just after the pulsar enters the radio eclipse region (the blue shaded regions). This can be interpreted as the Doppler boosting effect caused by the bulk flow in the downstream region.
The shock geometry is controlled by the ratio of the momentum fluxes of the pulsar wind to the stellar wind. For PSR B1957+20, such ratio is estimated to be ∼10 which suggests the companion star is confined by the pulsar wind and the shock [63]. The opening angle of the cone-like shock is ∼50-60 • , which corresponds to a separation of ∼0.15 in orbital phase. Because the emission is concentrated in the forward direction of the flow, we therefore expect that double peaks due to the Doppler boosting effect appear at the ∼0.15 phase before and after the phase of the radio eclipse (i.e., inferior conjunction), which is consistent with the observed result (see Figure 4) [63].
Apart from X-ray, evidence for orbital modulation of this black widow MSP is also found in a γ-ray regime. The γ-ray emission of PSR B1957+20 at energies larger than 2.7 GeV has been found to be modulated at the orbital period [64] (see Figure 5). On the other hand, the emissions below 2.7 GeV are steady and are dominated by the pulsar emission. Figure 5 clearly shows that the enhanced emission above 2.7 GeV appears around the inferior conjunction. It has been speculated that the modulated GeV emission originates from the inverse-Compton scattering of the thermal radiation of the companion star off the "cold" (i.e., low energy of the leptons in the co-moving frame of the plasma) ultra-relativistic pulsar wind [64].
Besides PSR B1957+20, only two other black widow MSPs (PSR J1311-3430 [45] and PSR J2241-5236 [72]) have their possible γ-ray modulations reported, which can possibly be a signature of intrabinary shock. However, their significance is not very high, including the prototypical case of PSR B1957+20. As more than 10 years of data have been accumulated by Fermi LAT, one can reexamine and check these γ-ray modulations together with the updated instrumental responses and the background model.
Currently, 27 black widows have been discovered in the galactic field (see Table 1) and 17 are found in globular clusters ( Table 2). The companion stars in these systems are likely to degenerate and have masses 0.1M [17,36]. For those in the galactic field, a number of them have their X-ray counterparts identified [20]. In Section 5, we will provide a general discussion on the X-ray and γ-ray natures of this population as a whole.  Orbital modulation of PSR B1957+20 in X-ray as observed by Chandra [63]. The error bars correspond to 1σ uncertainties assuming Poisson noises. The eclipses of the radio pulses occur at the orbital phases of 0.2-0.3 and 1.2-1.3 which are highlighted by the blue regions. The grey regions represent the phases for extracting the X-ray spectrum of PSR B1957+20 in the eclipsing region in Huang et al. [63]. Two orbital cycles are shown for clarity. Figure 5. γ-ray orbital modulation of PSR B1957+20 observed by Fermi large area telescope (LAT). This is discovered by Wu et al. [64]. The error bars correspond to 1σ uncertainties assuming Poisson noises. The shaded regions correspond to the phase of radio eclipse (i.e., 0.2-0.3 and 1.2-1.3). Two orbital cycles are shown for clarity.

Redbacks
Redbacks are close relatives of black widows, which are the other class of pulsar binaries that show intense interactions between the pulsars and the companion stars. They are named after the species of the poisonous spiders which originate from Australia. While the range of the orbital period spanned by these systems is similar to that of black widow MSPs (P b ≤ 20 h), their companions are late-type non-degenerate stars with masses of M c ∼ 0.2 − 0.4M which are significantly higher than that of black widows (M c 0.1M ) [17]. Let us revisit the end stages of LMXB evolution. When a neutron star has been spun up sufficiently in the accretion-powered phase, particle acceleration can be re-initiated in its magnetosphere and the MSP is hence turned on [3]. However, the radio pulses cannot be observed unless the accretion disk and the dirty environment around it have been somehow cleared. Mechanisms involve processes such as pulsar wind ablation and radiation pressure from the MSP can possibly help to clear up the environment [3,5]. However, the exact transition between the accretion-powered LMXB and the rotation-powered MSP has not been witnessed until a decade ago.
The first identified redback MSP, PSR J1023+0038, was initially identified as an LMXB FIRST J102347.6+003841 [106]. The source clearly showed an accretion disk before 2002 [107]. Since then, the disk was found to have disappeared. Subsequently, the radio pulsations have been detected in 2007 and PSR J1023+0038 revealed itself as a newly born MSP [27]. It is a long-sought missing piece in the evolutionary picture of compact X-ray binaries. Shortly after its discovery, γ-ray emission as well as the X-ray orbital modulation (see Figure 6) from PSR J1023+0038 have been subsequently detected [108]. Interestingly, PSR J1023+0038 shows that MSP might not be just the end point in the evolution of compact binaries. Since 2013 late June, the radio pulsation of this system has disappeared ( [109,110]) and a new accretion disk has been formed as shown by the re-appearance of strong double peaked Hα emission [111]. All this evidence indicate that this system was switching from a rotation-powered state back to an accretion-powered state.
Changes in all other wavelengths have been found to accompany this transition (see Figure 7) ( [28,29]). First, its γ-rays suddenly brightened within a few days in 2013 June/July and have remained at a high γ-ray state. Second, both UV and X-ray fluxes have increased by roughly an order of magnitude. Moreover, the system does not show any X-ray orbital modulation after the transition. On the other hand, the inset box shows the detailed evolution of the γ-ray emissions from 6 June to 24 July. Each data point of UV/X-ray represents an individual observation taken by Swift. Each γ-ray data points in the main panel and inset corresponds to two weeks and three days, respectively. In the cases where the detection significances is ≤ 3σ, upper limits at 95% confidence are given instead [28].
Assuming the MSP is still active and the absence of radio pulsation is simply a result of absorption/scattering due to increased local charge density, the aforementioned multiwavelength behaviour of PSR J1023+0038 can be explained with the model illustrated in Figure 8 [28]. In this scenario, a new accretion disk was formed as a result of sudden increase of the stellar wind. This can provide an explanation for all the other observed phenomena after the transition.
The enhancement of the UV emission can be naturally accounted for by the newly formed disk. The increase of the γ-ray flux by an order of magnitude suggests that a new emission component emerges after the transition. We propose that the inverse-Compton scattering process of the cold-relativistic pulsar wind off the optical/UV photons from the accretion disk produces the additional γ-rays. The accretion disk extends beyond the light cylinder radius (R lc ). R s is the distance to the intra-binary shock from the pulsar. R c is the critical distance from the pulsar at which the γ-rays from its magnetosphere evaporate the disk matter at R < R c ∼ 3 × 10 9 cm. UV/Optical photons mainly originate from the disk at R ∼ 10 9−10 cm. Shock is formed through the interaction between the pulsar wind and the stellar wind. This produces the non-thermal X-ray emissions. The inverse-Compton process of the cold-relativistic pulsar wind off UV/Optical photons from the disk produces the additional γ-rays [28].
Before the transition in June 2013, X-ray modulation was observed from PSR J1023+0038 where the emission originated from the intra-binary shock ( [108,112]). It has been suggested that during the MSP phase, the emission region is closer to the companion star and its orbital variations are caused by the eclipse of the emission region by the companion star [112]. After the transition, it is possible that the size of the emission region is getting bigger than that of before late June 2013. In such a case, the companion star can only block a negligible fraction of this emission and hence this explains the disappearance of the orbital X-ray variation. It has been speculated that the increase in the mass transfer from the companion star pushes the emission region back toward the pulsar, and more fraction of the pulsar wind is stopped by the shock, resulting in an increase in the X-ray emissions from the system [28].
In Figure 9, the multiwavelength spectral energy distributions (SEDs) of PSR J1023+0038 before and after the transition in June 2013 are compared [28]. Moreover, the theoretical models based on the aforementioned scenario are overlaid on the observed data. It appears that this scenario proposed by [28] can explain the results quite well.  Figure 9. Multi-wavelength spectral energy distributions of a PSR J1023+0038 system before (left) and after (right) late June 2013. Calculations with a model consist of emission components from the pulsar magnetosphere (outer gap); shock and pulsar wind (PW) are compared with the observed data before and after the transition. For further details, please refer to [28]. PSR J1023+0038 is an archetypal examples of a growing group of transitional MSPs (tMSPs). This includes PSR J1824-2452I, which resides in the globular cluster M28 [91]. Its X-ray counterpart is a transient IGR J1824-2452. In 2013, it was detected as an X-ray outburst with rotational and orbital ephemeris the same as PSR J1824-2452I. Interestingly, after a month long X-ray burst, the system was found to be reactivated as an active radio MSP within a few days. Another tMSP XSS J1227-4859 was previously identified as an LMXB with positional coincidence with Fermi LAT source 1FGL J1227.9-4852/2FGL J1227.7-4853 [113]. In December 2013, a transition was accompanied with a brightness drop in opticals, X-rays and γ-rays [114,115]. Radio searches at the beginning of 2014 eventually detected the pulses at 607 MHz from the system [116].
tMSPs are not the only redback that show dramatic changes. Very recently, we observed another redback system PSR J1048+2339 with the Lulin 1 m telescope in Taiwan and 2 m Liverpool telescope in La Palma [34]. Its orbital modulation in the optical regime has been found to change drastically in a timescale of less than two weeks (see Figure 10). From the observations taken in March 2018, the optical light curve folded at the orbital period of 6 h resembles an ellipsoidal modulation of the companion star which shows two peaks in a cycle at φ = 0 and φ = 0.55. Ellipsoidal modulation is a consequence of the orbital motion for a tidally distorted companion, which is commonly seen among redback systems (e.g., PSR J2129-0429 [117]). In Figure 10, the minimum at φ = 0.25 corresponds to the inferior conjunction.
When PSR J1048+2339 was observed again in early April 2018, its light curve was found to be completely different from that obtained in March (see Figure 10). It became a single-peaked sinusoidal profile with a maximum at φ = 0.65 and a minimum at φ = 0.25. Furthermore, the brightness of the companion is increased by one magnitude. All this suggests the modulation is now dominated by pulsar wind heating of the companion [34]. The change from an ellipsoidal modulation to a brightened sinusoidal profile occurred in less than 14 days [34]. The irradiation power of the system has been found to increase by a factor of six, and it has been speculated that this is related to the activity of the companion star [34]. The magnetic field of the companion star could play an important role by connecting to the shock region and guides the pair plasmas to the companion star surface.

X-Ray Properties
As the number of MSPs have significantly enlarged over the last decade, it is possible to carry out a statistical analysis of their population as a whole and examine similarities and differences among various classes. Recently, the authors of Lee et al. [20] performed an X-ray census of all MSPs in the galactic field. By utilizing all the available X-ray data, 47 MSPs have their X-ray counterparts detected (see Table 1 in [20]). On the other hand, upper limits have been placed on the X-ray emission of another 36 MSPs (see Table 2 in [20]). Using this censored data, the authors of Lee et al. [20] examined the empirical relation between their X-ray luminosities L x and spin-down powerĖ (see Figure 11). The best-fit relation is found to be L x 10 31.05Ė1. 31 35 erg/s in 2-10 keV (i.e., the Akritas-Thiel-Sen (ATS) line in Figure 11 [20]), whereĖ is the spin-down power with units of 10 35 erg/s.
Concerning the spider pulsars, there are eight redbacks and 12 black widows in the X-ray selected sample adopted in the study by Lee et al. [20]. By comparing their physical properties with the non-spider MSPs, while the distributions of the surface magnetic field strengths among different classes are found to be similar, the magnetic field strength at the light cylinder as well as the spin-down powers of redbacks are found to be significantly higher than those of non-spider MSPs in binaries.
Comparing the rotational and orbital parameters of redbacks and black widows, there is no significant difference found [20]. However, L x of redbacks are found to be significantly higher than those of black widows [20] (see Figure 12). Moreover, there is an indication that the X-ray emission of redbacks is harder than that of black widows [20] (see Figure 12). The difference in their X-ray luminosities can be explained by the different contributions from their intrabinary shocks [20]. The shock luminosity is proportional to L X ∝ δĖ, where δ is the fraction of the pulsar wind blocked by the outflow from the companion star and/or the companion star itself. δ can be estimated by the fraction of the sky intercepted by the companion star with δ ∼ (R R /2a) 2 , where R R is the Roche-lobe radius of the companion star. Since the Roche-lobe radius is estimated as R R /a = 0.462[q/(1 + q)] 1/3 with q being the ratio of the companion mass to the neutron star mass, the fraction can be expressed as δ = 0.053[q/(1 + q)] 2/3 . With the typical values of the mass ratio for both classes of spider pulsars, Lee et al. [20] estimate δ ∼ 1% for redbacks (q = 0.1) and δ ∼ 0.2% for black widows (q = 0.01). This may explain their difference in L x .

γ-Ray Properties
While the aforementioned work has investigated the X-ray properties of spider MSPs, there is no systematic analysis of their γ-ray properties that has been reported elsewhere. For a complete review on the high energy nature of these MSPs, we here present a statistical analysis of their γ-ray properties.
First, we investigate the L γ −Ė relation for the most updated sample of γ-ray pulsars, where L γ is the γ-ray luminosity in 0.1-100 GeV. We make use of the Fermi LAT 8-Year Point Source Catalog (4FGL) of the version updated on 29 May 2019 [11]. Among 5065 sources detected at a significance above 4σ, 239 of them are classified as pulsars. For computing their L γ andĖ, we need the estimates of their distance as well asṖ. Searching for such properties of the 4FGL pulsars from the ATNF catalog [1], we obtain a sample of 169 pulsars with their L γ andĖ plotted in Figure 13. By fitting a linear model to log L γ and logĖ, we obtain a best-fit relation of log L γ = (0.63 ± 0.05) logĖ + (11.87 ± 1.78) with the uncertainties estimated by bootstraping. We plot this relation as a straight line in Figure 13. The data points for the spider pulsars are highlighted by the coloured symbols. We found they follow the general L γ −Ė trend inferred from the entire γ-ray population.
Since the luminosities and the spectral properties of black widows and redbacks in X-ray regimes are found to be rather different [20], it is instructive to compare the corresponding properties of these two classes in γ-ray. The γ-ray spectrum of a pulsar is typically characterized by an exponentially cutoff power law with the following form: where E 0 is a chosen energy scale expressed in MeV, K is the prefactor which describes the flux density of the pulsar, Γ is the photon index which describes the low energy spectral slope, a and b are the exponential factor and and the exponential index which describe the cutoff at high energy. In the adopted 4FGL catalog, except for the six brightest γ-ray pulsars which were fitted with Equation (4) with b as free parameter, all the other pulsars are fitted with this function with b fixed at 2/3 [11].
In Figure 14, we construct the cumulative distribution functions for the γ-ray properties of black widows (blue symbols) and the redbacks (green symbols). We first compare L γ of these two classes of spider MSPs (top left panel of Figure 14). In contrast to the case for X-rays, the distributions of L γ of these two groups appear to be similar to each other. Applying a two-samples Anderson-Darling (A-D) test to their distribution, which is the same method adopted by [20] in comparing the X-ray properties of different classes of MSPs, we obtain a p-value of 0.6 and hence we conclude that there is no difference in the γ-ray luminosity distributions between black widows and redbacks. We have also compared theĖ of this γ-ray selected sample (top right panel of Figure 14) and we also do not find any difference between black widows and redbacks.
In the lower-left and lower-right panels of Figure 14, we compare the distributions of Γ and a between the groups of 23 black widows and 10 redbacks from the 4FGL catalog. While these plots show some differences in their γ-ray spectral properties, A-D tests yield the p-values of 0.16 and 0.30 in comparing their distributions of Γ and a, respectively. Hence, we conclude that the γ-ray spectral properties of black widows and redbacks are consistent with each other.
In this analysis, we do not find any significant difference in the γ-ray properties between black widows and redbacks. The different contributions from the intrabinary shocks in these two MSP classes, which lead to their differences in the X-ray, do not result in any notable difference in the γ-ray regime. This suggests that the X-rays and γ-rays from these spider MSPs originated from different mechanisms. we conclude that the γ-ray emission of the spider pulsars are dominated by their magnetospheric radiation.

Future Prospects
We have seen that a vast of dramatic high energy phenomena have been observed from the spider MSPs. Here we would like to give a wish list for what kind of investigations we can do and what we can expect to learn from these targets in the near future.
1. Catch the changes in the act: PSR J1023+0038 is a classic example of a redback and demonstrates the swings between accretion-powered state and rotation-powered state. The emission properties of the system have been found to change drastically along with the transition (Figures 7-9). Studying these variabilities can allow better understanding of the evolution of an LMXB in the recycling phase. However, it is impossible to catch such a transition without a long-term monitoring campaign. Almost all the spider pulsars are γ-ray emitters. Therefore, the continuous all-sky surveying mode of Fermi LAT provides us with the ideal instrument to monitor their behaviours. In Figure 15, we show the long-term γ-ray light curve of PSR J1023+0038 as observed by Fermi LAT. The transition in June 2013 (MJD 56450) can be clearly noted. Once a similar jump in the γ-ray flux has been spotted from any redbacks, a fast response X-ray/UV follow-up with the Neil Gehrels Swift Observatory can hence be triggered. 54500 55000 55500 56000 56500 57000 57500 58000 58500 59000 2. Globular clusters vs. Galactic field: In Section 5, we highlighted the results from the population analysis of the X-ray selected MSPs in the galactic field. We mentioned that the X-ray emission from the redbacks in the galactic field are apparently brighter and harder than the other classes of MSPs, including black widows ( Figure 12). In Tables 1 and 2, we have shown that the current sample sizes of the spider pulsars in the galactic field and in globular clusters are comparable. The formation processes of MSPs in globular clusters is dynamical [119], which can be very different from those in galactic field. It is instructive to compare the rotational, orbital and emission properties between these two populations. As the properties of redbacks in the galactic field are found to be different from the other classes, it is interesting to see if their counterparts in globular clusters also share similar behaviours. One can also compare the redbacks/black widows in these two populations to explore if the dynamical formation processes in globular clusters can result in any difference.
3. Spider hunting: Thanks to Fermi, the MSP population has been significantly enlarged through the synergy of picking MSP-like candidates from the unidentified γ-ray sources and the multiwavelength identification campaign. In our previous MSP-like candidates selection, we adopt a simple set of selection criteria which was based on our current knowledge of MSPs (see Section 2). However, the accuracy of picking the right candidates by this conventional method is unlikely to be optimal because the patterns in the data that define an MSP can be overlooked by human investigators. To improve the efficiency and the accuracy of the MSP identification campaign, machine learning techniques should be employed in γ-ray source classification.
We note that some studies have used such techniques to classify γ-ray sources (e.g., [120]). However, since the features for building the classification model were selected manually in these works, it is unlikely that the power of machine learning has been fully exploited. Recently, we have developed a scheme by coupling the classification algorithms with automatic feature extraction methods [121]. With this scheme, one can improve the prediction accuracy and pick the MSP-like γ-ray sources with higher confidence. Using the known γ-ray selected MSPs for testing our algorithm, we can attain an accuracy > 95% in classifying MSPs while all the previous attempts can only achieve at most ∼90%. Therefore, we expect that the multiwavelength follow-up observations of the targets selected by this method can result in more spider pulsars.
Very recently, a new X-ray mission eROSITA was launched, on 13 July 2019. 1 eROSITA is equipped with seven identical Wolter-I mirror modules. Each module consists of 54 nested mirror shells. It is going to perform an all-sky survey in 0.3-10 keV with unprecedented sensitivity. Therefore, it will provide a full X-ray coverage of all unidentified γ-ray sources. On the other hand, the upcoming large synoptic survey telescope (LSST) 2 will have an 8.4 m mirror with an exceptionally wide field of 3.5 • . It is capable of surveying the entire observable sky in the optical regime from its site in only three nights. It can provide deep multi-epoch data for the periodicity searches which enable us to identify the companions of the MSP-like objects.
In Table 3, we summarize the properties of a number of promising spider pulsar candidates. They were identified by the periodicities detected in X-ray/optical which resemble those orbital modulations found in redbacks/black widows. Three of these targets, 3FGL J0427.9−6704, FL8Y J1109.8−6500 and 3FGL J1544.6−1125, have shown stepwise jumps in their long-term X-ray light curves that are similar to the transitions between accretion-powered and rotation-powered states that have seen in PSR J1023+0038. Furthermore, Hα emissions have been detected from these three sources, which indicates the possibility of accretion. This makes them the candidates of tMSPs. Follow-up observations of them are strongly encouraged. 3 .  [47,[136][137][138] O = Significant optical orbital modulation; o = Possible optical orbital modulation; X = Significant X-ray orbital modulation; x = Possible X-ray orbital modulation; g = Possible gamma-ray orbital modulation.
4. Push to very high energy regime: In a very high energy regime (VHE; >100 GeV), neither the pulsed emission nor the pulsar wind nebula has been detected from any MSP so far. However, the optical and X-ray modulations observed from spider MSPs indicate the presence of ablation and heating of their companions. The interaction of the pulsar and stellar winds can create a shocked region in which particles can be accelerated [139]. Such accelerated electrons can possibly Comptonize stellar radiation from the companion to the TeV energies. This process can produce orbital modulated VHE γ-rays.
On the other hand, black widow PSR B1957+20 has a bow-shock X-ray nebula (See Figure 3). Assuming the observed non-thermal X-ray emission up to ∼ 10 keV is produced in a synchrotron process, this implies the magnetic field of the nebula of the order of µG and the presence of high energy electrons with energies up to the order of a hundred TeV [102]. These electrons can also possibly Comptonize the ambient soft photon fields and result in VHE γ-rays [140]. It is interesting to note that the limiting fluxes in the TeV regime placed by MAGIC are very close to the theoretically predicted values (see Figure 8 in [141]).
Cherenkov telescope array (CTA) is the next generation VHE observatory which consists of more than 100 telescopes located in both the northern and southern hemispheres. Upon completion of the construction of the full array, CTA will have a much larger collecting area, wider energy coverage and a larger field-of-view than any existing VHE-observing facility. With its unprecedented performance at such high energies, it is not unreasonable to speculate that a number of spider pulsars can be detected in the VHE regime. This will open a new widow for us to study these MSPs and help us to further constrain the acceleration processes of leptons within the binaries and their surrounding.