EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 608 (December 2017) A&A, 608 (2017) A59 Full HTML
Free Access
Issue
A&A
Volume 608, December 2017
Article Number A59
Number of page(s) 13
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/201731084
Published online 07 December 2017

© ESO, 2017

1. Introduction

Gamma-ray binaries are systems composed of a massive star in orbit with a compact object and are characterized by broad non-thermal emission peaking (in νFν) at energies above 1 MeV. The latter feature distinguishes gamma-ray binaries from high-mass X-ray binaries (HMXBs), whose spectral energy distribution peaks in X-rays, whereas the former feature distinguishes these binaries from recycled binary millisecond pulsars that have a low-mass companion. The compact object in gamma-ray binaries is likely to be a young, rotation-powered neutron star with non-thermal radiation due to the interaction of energetic pulsar wind particles with the stellar wind and radiation field of the O or Be companion. There is ample indirect evidence for this binary pulsar wind nebula scenario even though scattering in the stellar wind prevents detection of the expected radio pulsar in most gamma-ray binaries (see Dubus 2013, for a review). Accordingly, we explicitly assume in the following that the compact object in gamma-ray binaries is a pulsar. However, many of our results are equally applicable if the gamma-ray emission is powered by non-thermal jet emission from an accreting black hole (e.g. Massi et al. 2017). Clear evidence for gamma-ray jet emission exists for the accreting sources Cyg X-1 and Cyg X-3, but these objects are not gamma-ray binaries according to our definition because they are orders of magnitude more luminous in X-rays than in gamma rays.

There are six gamma-ray binaries detected in high energy (HE; 0.1−100 GeV) or very high energy (VHE; >100 GeV) gamma rays. Of these, two were initially detected as HE gamma-ray sources in all-sky surveys (LS I +61°303, Gregory & Taylor 1978; and 1FGL J1018.6-5856; Ackermann et al. 2012), two were independently detected in HE gamma rays and X-rays before the association was made (LS 5039, Paredes et al. 2000; and LMC P3, Corbet et al. 2016), one was detected serendipitously in VHE observations of the Monoceros Loop (HESS J0632+057; Hinton et al. 2009), and one was detected in a radio pulsar survey (PSR B1259-63; Johnston et al. 1992). Follow-up observations established that these sources are binaries harbouring a massive star and that their non-thermal emission is modulated on the orbital period. In addition to those six gamma-ray binaries there are also four systems discovered in radio surveys with a young pulsar in orbit with a massive star, but where variable gamma-ray emission associated with the binary has yet to be detected because of low pulsar power, a long orbital timescale, and/or large distance. These four systems are PSR J0045-7319, PSR J1638-4725, PSR J1740-3052, and PSR J2032+4127 (see Stairs et al. 2001; Bassa et al. 2011; Madsen et al. 2012; Lyne et al. 2015, respectively).

Gamma-ray binaries are probably a short-lived phase in the evolution of massive star binaries, following the birth of the neutron star and preceding the HMXB phase, when the neutron star accretes material captured from the stellar wind instead of holding it back (see Tauris & van den Heuvel 2006, for a review on the formation of compact objects in binaries). Accretion occurs if the ram pressure from accreting matter is able to overcome the pulsar wind, turning off the pulsar mechanism (Shvartsman 1971; Illarionov & Sunyaev 1975; Lipunov et al. 1994; Campana et al. 1995). A gamma-ray binary can thus transition to a HMXB on the typical spin-down timescale of young pulsars, i.e. a few 105 yr. The evolution of the companion eventually leads to a second supernova with the formation of another compact object. Therefore, besides the unique opportunities gamma-ray binaries provide to understand the physics of pulsar winds, these binaries also offer a window into the pulsar and orbital parameters of systems that remain bound after a supernova and constrain the formation paths to double neutron stars and coalescing compact objects.

Achieving these goals depends on our ability to explore the population of gamma-ray binaries. The number of gamma-ray binaries in our Galaxy has been estimated from a few dozen to a few thousand systems from population synthesis studies of HMXB evolution (Meurs & van den Heuvel 1989; Iben et al. 1995; Portegies Zwart & Verbunt 1996; Portegies Zwart & Yungelson 1998). Gamma-ray binaries are more likely to stand out in gamma rays rather than in radio, optical, or X-ray surveys where they are usually inconspicuous. The discovery of LMC P3 in the Large Magellanic Cloud suggests that we may have already accessed most of the observable gamma-ray binary population in our own Galaxy (Corbet et al. 2016).

Here, we aim to provide the first detailed estimate of the number of gamma-ray binaries based on HE and VHE observations and to evaluate the prospects for further discoveries. To do this, we simulated observations of gamma-ray binaries to assess the probability of detections in mock gamma-ray surveys that are designed to follow as closely as possible those performed or planned with the Fermi Large Area Telescope (Fermi-LAT), the High-Energy Stereoscopic System (H.E.S.S.), the High-Altitude Water Cherenkov observatory (HAWC), and the Cherenkov Telescope Array (CTA; Sect. 2). One difficulty in assessing the detectability is that the gamma-ray flux can vary strongly with orbital phase. We use input gamma-ray orbital light curves based on templates constructed from observations (Sect. 3) or based on a radiative model (Sect. 4). The estimated population size and expectations for future detections are discussed in Sect. 5.

2. Simulating surveys

We simulate a measurement as the flux average of the gamma-ray binary light curve integrated over a certain duration and energy range. The duration of the measurement, number of measurements (visits), and their distribution throughout time vary according to the type of instrumentation. The observability and detectability of the system depend on the assumptions made for each type of survey that was simulated. The observability only depends on the part of the sky surveyed and the location of the binary system. The detectability depends on the sensitivity of the survey, cadence of the visits, and emission properties of the system.

We simulate five types of surveys with properties as close as possible to existing or envisioned surveys, without carrying out a full end-to-end simulation of the observations and data analysis chain. In our opinion, current knowledge about the radiative mechanisms in gamma-ray binaries does not justify performing such complex and costly end-to-end simulations. The level of detail in our mock surveys is appropriate for the basic emission model that we develop in Sect. 4, which represents gamma-ray binary spectra at 1 GeV and 1 TeV with mono-energetic electrons. In the GeV domain, we simulate the Fermi-LAT third catalogue (“3FGL-like”) and the Fermi All-sky Variability Analysis (“FAVA-like”). In the TeV domain, we simulate the H.E.S.S. Galactic Plane survey (“H.E.S.S.-like”), a “HAWC-like” survey, and the CTA Galactic Plane survey (“CTA-like”). A source is considered detected if its average flux in ph cm-2 s-1 during the observation exceeds the threshold of the survey as defined below. We do not address the issue of how the detected gamma-ray source is identified as a gamma-ray binary, presumably through multi-wavelength follow-up observations. In particular, we make no attempt to investigate how binaries can be identified through a period analysis, such as that performed on the Fermi-LAT catalogue by Ackermann et al. (2012). Here, the orbital modulation only intervenes as the source of flux variability between observation windows.

2.1. The 3FGL-like survey

The 3FGL-like survey tests whether the binary would have been included in the third Fermi-LAT catalogue (Acero et al. 2015). The whole Galactic plane is covered so the gamma-ray binary observability is 100%. The measurement is assumed to last four years, ignoring any time variation in exposure. The threshold for detection is taken to be 10-9ph cm-2 s-1 (1−100 GeV) based on the flux distribution of sources detected within 10° of the Galactic plane in Acero et al. (2015, see their Fig. 24). We set the energy threshold at 1 GeV because the HE component of binaries peak around this energy and because the Galactic diffuse emission, which we do not take into account, is weaker than at 100 MeV. The 3FGL catalogue includes LS 5039, LS I +61°303, and 1FGL J1018.6-5856. LMC P3 is also part of the catalogue but is confused with other sources in the Large Magellanic Cloud. We also consider, where indicated, the impact of continued Fermi-LAT observations in the future. This extended 3FGL survey assumes a detection threshold that is lowered by a factor 2 and an exposure that is increased by a factor 4 (16 yr of observations).

2.2. The FAVA-like survey

The FAVA-like survey is based on the search for 5.5σ deviations from a long-term average model of the GeV emission observed with the Fermi-LAT (Abdollahi et al. 2017). Following the FAVA procedure, we search for deviations on a weekly timescale, which sets the duration of the simulated measurement, over a time span of eight years. Again, any time variation of the exposure is ignored and the observability is 100%. We consider the system to be detected by this survey if its weekly average flux subtracted from its long-term average flux (over eight years) is greater than 10-6.5ph cm-2 s-1(>100 MeV). Although the exact threshold changes as a function of location in the Galactic plane and spectrum, this choice should be conservative based in Fig. 4 of Abdollahi et al. (2017). The FAVA survey is potentially more sensitive than the 3FGL survey to systems such as PSR B1259-63 with long orbital periods and short duty cycles for GeV emission. The second FAVA catalogue includes LS I +61°303 and PSR B1259-63.

2.3. The H.E.S.S.-like survey

The H.E.S.S.-like survey is based on the survey of the Galactic plane carried out by the H.E.S.S. collaboration since 2004 and is meant to be representative of what current imaging arrays of Cherenkov telescopes (IACTs) can achieve. The survey covers only part of the Galactic plane, i.e. −110° ≤ l ≤ 65°. The observability of a system is decided by checking that it is observable for at least 2 h at some point in the year at a zenith angle smaller than 45°, assuming the geographical location of the H.E.S.S. array, and that its longitude is within the surveyed area. To produce a schedule of observations, we randomly distribute 25 visits of 2 h over a time span of eight years, i.e. we assume a uniform survey exposure of 50 h is achieved. We take into account that observations occur at night, ensuring each binary has a preferred observation season. However, we do not account for Moon-less operations, which influence the distribution of observable time on a monthly timescale. The latter effect averages out over a timescale of a few years, whereas the former (prefered season) does not. The measured fluxes from each visit are then averaged and compared to a detection threshold of 3.6 × 10-13ph cm-2 s-1 ( >1 TeV). This threshold corresponds to a flux of 20 mCrab1. The exposure times and sensitivity are comparable to those of the H.E.S.S. survey2. We also consider whether a detection could be claimed from a single visit, scaling the threshold by a factor (50/2)1/2.

2.4. The CTA-like survey

The CTA-like survey is intended to test the potential performance of the CTA array in detecting new gamma-ray binaries. The guiding principles are identical to the H.E.S.S.-like survey. We assume that the survey is divided up into two blocks carried out in parallel during the first two years of operations, based on the plans for an initial Galactic Plane survey by the CTA Consortium (Vercellone 2017). The first block, carried out by the south array in Chile, covers longitudes −60° ≤ l ≤ 60° down to a sensitivity of 2.7 mCrab using six visits of 2 h. The second block, carried out by the north array in the Canary Islands, covers 60° ≤ l ≤ 150° down to 4.2 mCrab in four visits of 2 h. We also consider the full survey covering the entire Galactic plane and carried out over a time span of 10 yr (see Fig. 6 in Vercellone 2017, for details). The observability of each system is decided as for the H.E.S.S.-like survey using the planned locations for the arrays.

2.5. The HAWC-like survey

Finally, we test for the detection of binaries using the extended air shower array HAWC. The high duty cycle and full-sky monitor capacity of HAWC can make it more sensitive to flaring gamma-ray binaries than IACTs such as H.E.S.S. and CTA. Here, the binary is observable if it transits with a zenith angle smaller than 45° at the location of the HAWC array in Mexico. We then simulate one measurement per day at the time of transit and with a duration equal to transit duration. The HAWC sensitivity after five years of operation is comparable to that achieved in the H.E.S.S. Galactic Plane survey, i.e. 20 mCrab above 1 TeV, assuming a source transit duration across the sky of 6 h (Carramiñana 2016). The threshold for daily detection is close to 1 Crab for a 6 h transit, corresponding to the transit time of a source that passes close to zenith, i.e. with a declination close to +19°. The dependence of the threshold with source declination is taken into account using the curve showing sensitivity as a function of declination for a E-2.5 spectrum in Fig. 10 of Abeysekara et al. (2017). We test for detection in each daily measurement and in the accumulated exposure over five years of HAWC operations.

thumbnail Fig. 1

Template light curves for known gamma-ray binaries. Two orbits are shown except for PSR B1259-63 where the plot focusses on periastron passage. The various measurements are shown in grey, with error bars omitted for clarity. Arrows indicate measurement upper limits. Left: photon flux integrated above 1 GeV in units of 10-9 ph cm-2 s-1 based on Fermi-LAT measurements. Right: photon flux integrated above 1 TeV in units of 10-12 ph cm-2 s-1 based on IACT measurements. The GeV and TeV data are taken from Abdo et al. (2009) and Aharonian et al. (2006) for LS 5039,  where ΓGeV = 2.54; Hadasch et al. (2012) and Acciari et al. (2011) for LS I +61°303, where ΓGeV = 2.42 and ΓTeV = 2.6; Bordas et al. (2016) for PSR B1259-63, where ΓGeV = ΓTeV = 2.7; Ackermann et al. (2012) and Abramowski et al. (2015) for 1FGL J1018.6-5856, where ΓGeV = 3.1 and ΓTeV = 2.7. The GeV data for LMC P3 are from Corbet et al. (2016), where ΓGeV = 2.8. The TeV data for HESS J0632+057 are from Aliu et al. (2014).

3. Extrapolating from observed gamma-ray binaries

We currently have five binaries with measured orbital modulations at both GeV and TeV energies and one with a GeV modulation (LMC P3). We aimed to find out the sensitivity of the surveys to the detection of these binaries. We estimated this sensitivity here by constructing a template light curve for each of the known gamma-ray binaries and, after proper scaling for distance, deriving the probability for detection once the binary was randomly located in the Galaxy.

3.1. Template light curves

Figure 1 shows template light curves for each of the known gamma-ray binaries based on the GeV and TeV observations available at the time of writing. In most cases, we simply took the phase-folded measurements and interpolated these measurements using splines. For PSR B1259-63, LS I +61°303, and HESS J0632+057, the error bars, phase coverage, or orbit-to-orbit variations make it difficult to assess the mean orbital light curves. In these cases, our templates are meant to be representative of the behaviour of the system in that they roughly capture the amplitude and phase variations that have been observed. The template GeV and TeV light curves are given in ph cm-2 s-1 integrated above 1 GeV and 1 TeV, respectively. We converted to these units assuming a simple power law when the data were not directly available in this format. The source of the data and the power-law photon index Γ (with dN ∝ E− ΓdE) that we used can be found in the caption to Fig. 1. Given the low statistics, the GeV light curve of HESS J0632+057 (not shown in Fig. 1) is described as a two bin light curve (orbital phases 0.0−0.5 and 0.5−1.0) using the spectral parameters in Table 1 of Li et al. (2017), who have reported the first detection of this system at GeV energies3.

Table 1 lists the detected systems for each mock survey presented in Sect. 2, given the template light curves and locations of the known binaries in the Galaxy. LS I +61°303 is not detected in the H.E.S.S.-like survey owing to its location. LS 5039 and LS I +61°303 are observable with HAWC but are undetected in the HAWC-like survey because of their unfavourable declinations; this non-detection is consistent with the 18 months of HAWC observations that have been reported to date (Abeysekara et al. 2017). PSR B1259-63 is always detected in the FAVA survey. The system is just below the threshold of the 3FGL survey when the four-year survey time span includes only one periastron passage of the 3.4 yr orbit, as observed. PSR B1259-63 has a very small probability ( < 1%) of being detected in the H.E.S.S. survey because of its low duty cycle and flux, but this source a 50% chance of being detected in the CTA-like survey. HESS J0632+057 is outside the H.E.S.S. and initial CTA-like surveys; this source is detected in the full CTA survey. 1FGL J1018.6-5856 is only detected in the 3FGL-like survey and in the full CTA-like survey. LMC P3 is only detected in the 3FGL-like survey. These results are fully consistent with the actual 3FGL, FAVA, and H.E.S.S. survey observations.

Table 1

Detected systems in the mock surveys based on the template light curves in Fig. 1.

3.2. Galactic distribution

We assumed that gamma-ray binaries are located in or close to the spiral arms of our Galaxy, similar to the O and B stars and HMXBs to which they are directly related. The Galaxy is modelled as four one-dimensional spiral arms. We used the arm formula of Ringermacher & Mead (2009) with parameters adjusted to reproduce the Galactic structure in Fig. 5 of Russeil (2003). Our Sun is 8 kpc away from the Galactic centre. Binaries are spread out uniformly across the Galaxy disk (15 kpc), keeping only those within 1 kpc of a spiral arm and more than 3 kpc away from the Galactic centre to account for the older stellar population in the bulge (see Fig. 2). The binaries are assumed to reside in the Galactic plane (b = 0°). The model Galactic longitude distribution (Fig. 2) compares well to the HMXB longitude distribution plotted in Grimm et al. (2002) or Walter et al. (2015).

The ground-based surveys (H.E.S.S., HAWC, and initial CTA) access only part of the Galactic plane, hence only a fraction of the binaries are observable for them. These fractions are given in Table 2 for both a distribution along spiral arms, as described above, and a strictly uniform disk distribution. This makes little difference. In the following, we consider only the spiral arm distribution.

thumbnail Fig. 2

Top: map of randomly generated locations for gamma-ray binaries in the Galaxy. Bottom: longitude distribution of gamma-ray binaries under the assumptions described in Sect. 3.2.

Table 2

Fraction (%) of observable binaries in each survey.

Table 3

Fraction of detected systems in each survey using the light curves in Fig. 1 as templates (see Sect. 3.3).

3.3. Detectable fraction based on observed light curves

Table 3 reports the detectable fraction of gamma-ray binaries in the various surveys based on a sample of 104 systems randomly distributed in Galactic location and in starting orbital phase for each template light curve in Fig. 1 after scaling for distance. The light curves are then run through the simulated observations of each survey described in Sect. 2 to test for detection. For example, the 3FGL-like survey detects 77.9 ± 0.8% of the binaries that are assumed to have a GeV light curve identical to LS 5039. Here and in the remainder of the paper, the errors represent the 95% confidence interval derived from a likelihood analysis4.

1FGL J1018.6-5856 and LMC P3 have high enough luminosities that such systems are detectable anywhere in the Galaxy with the 3FGL survey. Unsurprisingly, the FAVA-like survey is best suited for low duty cycle light curves such as PSR B1259-63. Our analysis recovers that PSR B1259-63 is detected in FAVA, but would not include LS I +61°303. The latter is inconsistent with the FAVA catalogue and is the result of neglecting the orbit-to-orbit variations that are detected in this source (Ackermann et al. 2013). The FAVA-like analysis also fails to detect LS 5039 and 1FGL J1018.6-5856 since the amplitude of the flux variations on weekly timescales are insufficient to flag them. In this case, this is consistent with the actual FAVA results (Abdollahi et al. 2017).

The numbers remain small in the H.E.S.S. and HAWC-like surveys. The high sensitivity of the CTA array leads to detection fractions that are comparable to those achieved by the Fermi-LAT. Very few systems are detected as transients in the ground-based surveys (i.e. detected only in one visit); the flux increase compared to the average in a highly eccentric system rarely compensates the higher sensitivity threshold for observations on a shorter duration. Hence, the fraction of detected systems in VHE surveys decreases with longer orbital period even if the systems have comparable maximum TeV luminosities owing to a high eccentricity.

The fractions in Table 3 give an estimated detection probability from which we can constrain the maximum size of the underlying population5. Any other existing system with identical properties to 1FGL J1018.6-5856 or LMC P3 would have been detected since the detection probability is 1. For LS 5039, knowing that the system is detected in the 3FGL survey, the 77.9% probability implies with >95% confidence that there are at most three systems with comparable light curves in our Galaxy and most likely only one. The same conclusion is reached for LS I +61°303. For PSR B1259-63, given the FAVA detection, the 12.6% detection probability implies 7+26-6\hbox{$^{+26}_{-6}$} such systems in our Galaxy. Therefore, on average, there may be one more LS 5039 or LS I +61°303 type system, and six other PSR B1259-63-like systems in the Galaxy that could have escaped detection in the Fermi-LAT data.

The VHE detection probabilities are not as constraining as those of HE except for HESS J0632+057. The detection probability is only 0.8% in both the 3FGL and H.E.S.S.-like surveys. The lack of detection in those surveys (HESS J0632+057 is outside the H.E.S.S. survey area; Sect. 3.1) places an upper limit of <231 on the number of HESS J0632+057-like systems in the Galaxy. The initial CTA-like survey should detect 11-6+8\hbox{$11^{+8}_{-6}$} of those 231 systems or reduce their estimated number to 8-7+30\hbox{$8^{+30}_{-7}$} should it only detect HESS J0632+057 after the full 10 year Galactic Plane survey. The Cherenkov Telescope Array Galactic Plane survey will thus be able to strongly constrain the number of such systems.

4. A synthetic population

In the preceding section we estimated the number of existing gamma-ray binaries from the properties of the known systems. However, these systems represent only the upper end of the luminosity function of gamma-ray binaries. In this section we estimate this number from a synthetic population model. Building this population requires a model for the gamma-ray emission of binaries, which is a bold enterprise given the current knowledge. While there is general agreement that anisotropic inverse Compton scattering of photons from the star and γγ pair production at TeV energies must play a role, since these processes naturally lead to orbital modulations the details vary significantly from model to model. Modulated Doppler boosting is also very likely to intervene if the emission occurs in a pulsar wind bow shock. Reproducing the orbital phases of gamma-ray detections in systems with Be companions, such as PSR B1259-63, has proven particularly difficult; this is possibly because of the complex interaction between the pulsar and the circumstellar material surrounding its companion. In the following, we adopted a simple model with the intention of minimizing the number of parameters while still being able to produce orbital light curves comparable to the observed light curves.

4.1. Orbital parameters

The binary eccentricities e were assumed to follow the thermal distribution (Ambartsumian 1937) p(e)de = 2ede with the additional conditions that e <emax = 1−(Porb/ 2 days)− 2/3 to ensure that the companion does not fill more than 70% of its Roche lobe at periastron and that the binaries are circularized (e = 0) for Porb ≤ 2 days (see Moe & Di Stefano 2017, and references therein). The inclination of the system is derived by randomly picking a vector on a sphere. The argument of periastron and the orbital phase at the time of the first simulated observation are picked from a uniform distribution between 0 and 2π. Finally, we uniformly sampled the logarithm of orbital periods between 1 and 104 days in order to assess the fraction of detected systems as a function of Porb, except in Sect. 4.4 where this is slightly modified for a more realistic representation of the Porb distribution of HMXBs.

4.2. Radiation model

We assumed that the radiation is due to Compton upscattering of stellar photons with an initial energy ≈ 10 eV. The GeV (resp. TeV) emission then requires electrons with a Lorentz factor γ = 104 (resp. γ = 106). For simplicity, we assumed monoenergetic distributions at these energies. This is supported by the observed GeV spectra of gamma-ray binaries, which generally consist of a hard power law with an exponential cut-off around 1 GeV. This is also admissible in the TeV range where soft power laws are observed such that most of the photons have an energy close the threshold energy of the VHE observations. The true particle distributions are likely to be more complex, but the GeV and TeV emissions are dominated by electrons of these energies and assuming more complex distributions (power laws, see Sect. 5) does not have a significant impact on the results.

We computed the inverse Compton bolometric power radiated by these particles, assuming that they are located at the position of the compact object. If the electron distribution is isotropic, the light curve in the Thomson approximation for Compton scattering is given by Lγ=NeσTcU⋆(1−βμ)[(1−βμ)γ2−1],\begin{equation} L_\gamma=N_{\rm e} \sigma_{\rm T} c U_\star (1-\beta\mu) \left[(1-\beta\mu)\gamma^2-1\right] \label{star} , \end{equation}(1)where U⋆=(1/c)σSBT⋆4(R⋆/d⋆)2\hbox{$U_\star=(1/c) \sigma_{\rm SB} T^4_\star (R_\star/d_\star)^2$} and T⋆ is the star temperature, R⋆ its radius, d⋆ its distance to the particles, and μ = cosθ represents the angle between the line of sight and the binary axis. The angle θ varies from π/ 2 + i (superior conjunction) to π/ 2−i (inferior conjunction) with i the system inclination. The massive star was assumed to have a radius of 10 R⊙ and temperature of 33 000 K. The analytic formula is valid for γ = 104, where the Thomson approximation is acceptable. However, stellar photons scatter in the Klein-Nishina regime when γ = 106. Hence, we numerically integrated the Compton kernel to derive the anisotropic emitted power instead of using Eq. (1) (see Dubus et al. 2010).

The total number of electrons Ne is related to the injected power in particles Pinj by Ne=Pinjγmec2×min{τesc,τic},\begin{equation} N_{\rm e}=\frac{P_{\rm inj}}{\gamma m_{\rm e} c^{2}} \times \mathrm{min}\left\{ \tau_{\rm esc},\tau_{\rm ic}\right\} , \end{equation}(2)where τesc is the escape timescale of the particles from the gamma-ray emitting region (see below) and τic is the inverse Compton loss timescale, which in the Thomson regime is τic=γmec243σTcU⋆γ2·\begin{equation} \tau_{\rm ic}=\frac{\gamma m_{\rm e} c^2}{\tfrac{4}{3}\sigma_{\rm T} c U_\star \gamma^2} \cdot \end{equation}(3)Hence, ⟨ Lγ ⟩ = Pinj (integrated over all angles) if the particles radiate efficiently before they leave the vicinity of the star (τic ≤ τesc), otherwise the radiated power is reduced to the fraction of particles that are in the emission zone ⟨ Lγ ⟩ = (τesc/τic)Pinj. The latter can be rewritten using Kepler’s third law as Lγ∝1/d⋆∝Porb−2/3\hbox{$L_\gamma \propto 1/d_{\star} \propto P_{\rm orb}^{-2/3}$}, hence there is a break in the distribution of ⟨ Lγ ⟩ /Pinj as a function of Porb for the orbital period where τic = τesc. This can be seen in Fig. 3, in which the mean of the average orbital luminosity is plotted for a sample of 104 binaries with orbital periods ranging from 1 to 104 days and randomly sampled eccentricities. The break is at Porb ≈ 10 days because we decided to set τesc = d⋆/c. Such a fast escape timescale is reasonable in the context of gamma-ray binaries, where the accelerated particles flow away relativistically in a bow shock (Dubus et al. 2015). This assumption is also conservative in that it may underestimate the number of detections by minimizing the radiative efficiency. The influence of this choice on the results is further discussed in Sect. 5.2.

thumbnail Fig. 3

Mean orbit-averaged gamma-ray flux, normalized to the injected power, as a function of Porb (see Sect. 4.2).

Since TeV photons are likely to create pairs by interacting with photons from the star, we calculated the expected γγ absorption at 1 TeV in the point source limit following Dubus (2006). Absorption reduces the average TeV luminosity for short orbital period systems, where the stellar radiation density is highest, resulting in a strong decrease of LTeV with Porb below 10 days (Fig. 3). The GeV emission is not affected by γγ absorption. However, we also took into account eclipses of the (point-like) gamma-ray emission zone by the star. This results in a slight decrease of the average GeV power at short Porb, instead of the expected flat distribution ⟨ Lγ ⟩ = Pinj. Hence, this model predicts the radiative efficiency is maximum for systems with Porb ≈ 10 days.

Figure 4 shows the average fractional amplitude of the model TeV light curves, measured as (fmax−fmin)/(fmax + fmin), where f is the flux. The mean amplitude increases slightly from short to long orbital periods owing to the larger eccentricities permitted (see Sect. 4.1) but eclipses and γγ absorption strongly increase the amplitude at short Porb. The average TeV variability amplitude at long Porb is about 80%, implying fmin ≈ 0.11fmax. Figure A.1 shows examples of GeV and TeV light curves that are computed via the radiative model described in this section.

thumbnail Fig. 4

Mean fractional amplitude of the simulated GeV (dashed line) and TeV (dash-dotted line) gamma-ray light curves (see Fig. 3) and fraction of the systems showing eclipses as a function of Porb (thick line).

Table 4

Fraction of detected systems in each survey using synthetic light curves as templates (see Sect. 4.3).

4.3. Detectable fraction based on synthetic light curves

To check for consistency with the results of Table 3, we produced 104 synthetic light curves using the orbital period and eccentricity (when known) for each observed system, i.e. leaving the system orientation free. We then normalized the synthetic light curves to the maximum observed luminosity. The systems are distributed throughout the Galaxy. The detection fractions in Table 4 are within a factor 2 or less of those in Table 3, showing comparable trends when looking at objects, orbital period, or surveys. The exception is PSR B1259-63 where the detected fraction in the 3FGL-like survey is a factor 5 higher because the model typically produces a lower amplitude light curve than observed, hence a higher average flux (see below). Despite this shortcoming, our simple radiative model should still be able to yield realistic estimates of the average detection rate for a population of systems.

We then produced synthetic light curves for a sample of binaries with random orbital parameters and a given injected power. Figure 5 shows the fraction of systems detected in the mock HE and VHE surveys discussed in Sect. 2, as a function of Porb and Pinj. The FAVA-like survey is much less efficient at detecting systems than the 3FGL-like survey. At short orbital periods, the sensitivity is insufficient to detect systems on a time span of a week. At long orbital periods, the amplitude of the variations in the model light curves (Fig. 4) is insufficient to provide a significant advantage to this burst search strategy compared to the integration strategy employed in the 3FGL-like survey. The latter is extremely efficient when the injected power in HE-emitting particles exceeds 1035 erg s-1, even for long Porb compared to the integration time (four years).

The VHE surveys access only part of the Galactic plane so their maximum efficiency does not reach 100% even for high injected powers in VHE-emitting particles. The results show comparable efficiencies for the H.E.S.S.- and HAWC-like surveys. The design of these two surveys, notably the visit frequencies, does not appear to play a major role in the detectable fraction; the peak at Porb ≈ 10 to 100 days simply reflects the higher radiated luminosity expected for those orbital periods in the model (see Fig. 3). The CTA-like survey is much more sensitive, detecting nearly all accessible systems for Pinj ≥ 1035 erg s-1 regardless of orbital period. Again, the sensitivity at long Porb results from our model, which on average gives a minimum flux around 11% of the maximum flux (Fig. 4 and Sect. 4.2). This enables the detection of long orbital period systems even when the phases of maximum flux are not sampled by the visits.

thumbnail Fig. 5

Detected fractions in the HE (left panels) and VHE (right panels) surveys. Each panel contains four curves corresponding to Pinj = 1033,1034,1035,1036 erg/s (dark to light blue sequence in each panel). The detection fraction increases when going from small Pinj (darker blue) to high Pinj (lighter blue) until the detected fraction saturates at the value given by the observable fraction (Table 2).

thumbnail Fig. 6

Orbital period (left panels) and spin-down power distribution (right panels) of a random sample of 105 systems. Top panels show the full population as a light blue histogram with the line-filled fraction showing the distribution of the binaries that are accreting according to Eq. (4). The bottom panels zoom in to highlight the systems detected in any of the present-day surveys (combining the H.E.S.S., 3FGL, and FAVA-like surveys; medium blue histogram) or any of the future surveys (combining the full CTA, extended 3FGL, and HAWC-like surveys, dark blue histogram). The latter essentially shows the systems detected in any survey since a binary detected in one of the present-day surveys has a nearly >99% chance of being re-detected in one of the future surveys.

4.4. Full population model

A full population model requires assumptions on the injected power Pinj and how it relates to the total available power Ė measured by pulsar spin-down. The parameter Pinj is likely to be different for the GeV and TeV emitting particles, whether they arise from different populations or from the same power-law distribution. We used PSR B1259-63, which is the only system with a measured Ė = 8 × 1035 erg s-1, to estimate the power going to the GeV- and TeV-emitting particles. Simulating GeV and TeV light curves with the same orbital period and eccentricity as PSR B1259-63 (i.e. following the procedure described in Sect. 4.3), we found that injection fractions PGeV = 0.07Ė and PTeV = 0.01Ė are needed to reproduce, on average, the peak gamma-ray fluxes listed in Table 3, and we adopted these values in the following.

We built the mock population by randomly sampling probability distributions of Ė and Porb. Following Lutovinov et al. (2013), we took a flat distribution in log Porb tapered by Gaussian edges at log Porb (days) = 1.3 and 3.7. This probability distribution (Fig. 6, top left panel) results from the evolution of pre-HMXB binaries (Bhadkamkar & Ghosh 2012). For Ė, we took as input the distribution of spin-down powers extracted from the ATNF pulsar catalogue (Manchester et al. 2005)6, selecting only those pulsars with a pulse period >10 ms and a spin-down timescale < 107 yr to exclude recycled millisecond pulsars. The resulting Ė distribution is shown in the top right panel of Fig. 6.

Table 5

Detection fractions for the population shown in Fig. 6.

We calculated the detection fraction in the various surveys from a random sample of 105 systems (Table 5). The distributions of detected systems in present-day or future surveys are shown in the zoomed-in bottom panels of Fig. 6. The detection fractions are biased towards short Porb and high Ė, as expected from the results of Sect. 4.3. The population model naturally accounts for the existence of radio pulsars in binaries that remain undetected in gamma rays because of their long orbital periods and low spin-down powers. PSR J0045-7319, PSR J1638-4725, PSR J1740-3052 with Porb = 51,1941,231 days (resp.) and Ė = 0.2, 0.4, and 5 × 1033 erg s-1 (resp.) are examples of such systems, which we do not expect to be readily detectable by the gamma-ray surveys (Stairs et al. 2001; Bassa et al. 2011; Madsen et al. 2012).

The pulsar wind pressure can be insufficient to hold off accretion from the stellar wind of the companion for low Ė and short Porb. In such a case, we considered that the system is an accreting HMXB. The accreted material can reach the neutron star surface or the neutron star magnetic field can propel material out, depending on the respective locations of the co-rotation radius and magnetospheric radius, both of which are within the light cylinder. In both cases we considered that the system does not emit gamma rays. We used a simple criterion to test whether a system is accreting or not, assuming the massive star wind is isotropic, uniform with the same constant mass loss rate Ṁw = 10-6M⊙ yr-1, and the velocity vw = 1000 km s-1 for all systems. The system is accreting if the pulsar spin-down power is less than Ė <4×1033(Ṁw10-6 M⊙ yr-1)(103 km s-1vw)3(0.1 AUap)2erg s-1,\begin{equation} \dot{E} < 4\times 10^{33} \left(\frac{\dot{M}_{\rm w}}{10^{-6}\,M_\odot\,\rm yr^{-1}}\right)\left(\frac{10^{3}\rm\,km\,s^{-1}}{v_{\rm w}}\right)^3\left(\frac{0.1\rm\,AU}{a_p}\right)^2\rm erg\,s^{-1} \label{eq:acc} , \end{equation}(4)where ap is the binary separation at periastron (Illarionov & Sunyaev 1975). This criterion is simplistic in regards to the complex physics of wind launching, capture, Be circumstellar disks, etc. (Dubus 2013), but we chose values of Ṁw and vw that are likely to overestimate the fraction of accreting systems. We find about 23% of the sampled systems are accreting, mostly at short Porb and low Ė as shown by the line-filled histogram in the top panels of Fig. 6. Despite this, we find negligible overlap with the population of systems detected in the HE and VHE surveys because these select high Ė systems; < 0.4% of the detected systems are also flagged as accreting.

About 40% of the binaries that are detected in HE can be found in both the 3FGL and FAVA-like surveys, whereas less than 2% are detected only in the FAVA-like survey. The detection of PSR B1259-63 in FAVA without a concurrent detection in the 3FGL survey is therefore unlikely in our model, as the statistics in Sect. 4.3 already showed. However, PSR B1259-63 is close to our crude 3FGL detection threshold using the template light curve so details in the 3FGL detectability may come into play (orbit-to-orbit fluctuations, Galactic diffuse emission). The probability to detect a system in one of the VHE surveys (H.E.S.S., HAWC, and CTA-like) is 4.23 ± 0.13%, with most of the detections arising from the CTA-like survey. Altogether, the probability to detect a gamma-ray binary in any of the surveys is 5.32 ± 0.14% (3FGL-, FAVA-, H.E.S.S.-, HAWC-, or CTA-like). Only a very small number are detected in VHE surveys without a detection in the HE surveys with this model. For instance, all of the systems detected by the H.E.S.S.-like survey are also detected by the 3FGL-like survey.

5. Discussion

5.1. Estimated population of gamma-ray binaries

We explored two ways to estimate the number of gamma-ray binaries. The first (Sect. 3.3) employed the light curves of the known systems as templates to evaluate the detection probability of identical systems distributed throughout the Galaxy. The results show that systems such as LS 5039, LS I +61°303, LMC P3, and 1FGL J1018.6-5856 are already detectable throughout most of the Galaxy, therefore it is highly unlikely that more than one or two have escaped notice (Table 3). Future HE detections are more likely to be of low duty cycle systems, such as PSR B1259-63, the total number of which is estimated at 7+26-6\hbox{$^{+26}_{-6}$}. One such anticipated detection is that of PSR J2032+4127, which is a 2 × 1035 erg s-1 pulsar in an eccentric, >20 yr orbit around a Be star that will pass periastron in late 2017 (Lyne et al. 2015; Ho et al. 2017)7. The largest source of uncertainty is the number of HESS J0632+057-like systems with a ratio of TeV to GeV luminosity about two orders of magnitude higher than the other binaries (see Table 3 and Li et al. 2017). There may be as many as ≈230 such systems in our Galaxy; this is an upper limit that CTA will decrease to 8-7+30\hbox{$8^{+30}_{-7}$} if none are discovered in the full Galactic Plane survey (Sect. 3.3).

The known systems represent only the upper end of the luminosity function of gamma-ray binaries. Thus, our second estimate for the number of gamma-ray binaries employed a full population model based on a series of assumptions on the radiative process, distributions of orbital parameters, and injected power (Sect. 4). In the HE domain, with four systems in the 3FGL and FAVA surveys (Table 1, excluding LMC P3 since it is not in the Galactic plane), the total parent population is estimated at 82+108-56\hbox{$^{+108}_{-56}$} systems based on the detection fraction in Table 5. In the VHE domain, with only LS 5039 detected in the H.E.S.S. survey, the parent population is constrained to 85+290-81\hbox{$^{+290}_{-81}$} systems. Combining all the information in Table 1 into the likelihood function, i.e. assuming four systems in the HE surveys, one system in H.E.S.S., none in HAWC, and at least five in the full CTA survey, the population is estimated at 101+89-52\hbox{$_{-52}^{+89}$} gamma-ray binaries in our Galaxy. These numbers are consistent with the predictions from population synthesis of HMXBs (Sect. 1).

Gamma-ray surveys are ≥ 50% complete for Ė ≥ 1036 erg s-1 (Fig. 6), but they access only a handful of systems in a population of about a hundred binaries. A few additional systems, such as PSR J2032+4127, may be detected through their pulsed gamma-ray emission without showing binary-related gamma-ray emission. We have not attempted to take this into account. The spin-down distribution of detected Fermi-LAT pulsars peaks at log Ė = 35.5 (see Sect. 5.2 below), suggesting this is unlikely to make a difference to the number of systems detected in gamma rays. A couple dozen binaries may be visible as accreting X-ray sources, indistinguishable from other HMXBs except perhaps through their neutron star spin periods or through propeller-induced behaviour. SAX J0635+0533 (Cusumano et al. 2000) and A0538-66 (Skinner et al. 1982) are possible examples. These two systems clearly have much faster spin periods ( < 70 ms) than all the other known X-ray pulsars in HMXBs ( >1 to 1000 s), suggesting that the neutron star may not yet have spun down significantly from its birth period.

PSR J0045-7319, PSR J1638-4725, and PSR J1740-3052 are representative of the low Ė systems that represent the majority of the pulsar plus massive star population: ≈ 55% of the sampled systems have Ė ≤ 1034 erg s-1 and are not accreting. Adding in the 23% that are accreting, this implies that 78% of the population is inaccessible to gamma-ray surveys. Estimating their detection rate in radio (SKA) or X-ray surveys (eROSITA) is beyond the scope of this work, but we note that the long Porb, high eccentricity systems are clearly more susceptible to be detected as radio pulsars (Lipunov et al. 1994), providing a complementary way to access the pulsar plus massive star population.

5.2. Systematic uncertainties in the population synthesis

How dependent are our results on the assumptions of the model? The Galactic distribution and binary parameters should not be a major source of concern since these have already been scrutinized in population studies of high-mass X-ray binaries (Walter et al. 2015). The distribution of Ė for gamma-ray binaries is entirely unknown and taking as input the Ė distribution of young pulsars in the ATNF catalogue probably suffers from a variety of selection biases; such biases are notably because it is not obvious that the birth spin period and evolution should be identical in isolated pulsars and binaries – if only because mass loss and kick during the supernova are necessarily weaker if the newly born neutron star is to remain bound to its companion. Yet, our assumption on Ė is not likely to have a major impact on gamma-ray observations since these are mostly sensitive to the high end of this distribution. The Ė distribution of detected binaries (Fig. 5) actually resembles the Ė distribution of young pulsars detected in gamma rays with the Fermi-LAT8, which we find to be well approximated by a Gaussian centred at log Ė = 35.5 with a standard deviation σ = 1. Taking this distribution as input increases the detection fraction, decreasing the population size inferred from current observations without changing much the number of expected detections in future surveys. However, this distribution cannot account for the known radio pulsars in orbit around massive stars with low Ė. The strongest impact of our assumption on the Ė distribution is therefore on the relative numbers of pulsar plus massive star binaries that are found in radio and gamma-ray surveys.

The light curves of gamma-ray binaries have proven difficult to model, even in the cases in which we have the most information, questioning the validity of our radiative model. For instance, relativistic beaming of the emission is thought to be an important factor in shaping the light curves (e.g. An & Romani 2017). A refined light curve model is desirable but may not change our results much. First, despite its simplicity, the detection fractions inferred from the model are broadly consistent with those inferred from the observed light curves (Sect. 3.1). Its main shortcoming is that it predicts lower amplitudes than observed, overestimating detection rates for PSR B1259-63-like systems. However, these do not dominate the detected systems (Fig. 6). Second, we experimented with a more complex radiative model, using a power-law distribution of particles and including Doppler boosting (assuming a particle bulk velocity of c/3 directed away from the star as in Dubus et al. 2010). There was surprisingly little difference between the detection fractions shown in Fig. 5 and those obtained with the more complex model despite substantial changes to the light curves from relativistic Doppler boosting, indicating that the flux level is more important than the detailed shape of the light curve in setting the detection fractions.

Improvements to our radiative model should thus concentrate on the injected power Pinj and the radiative efficiency τesc/τic, both of which set the flux level. A longer escape timescale increases the radiative efficiency but this needs to be compensated by a lower fraction Pinj/Ė in order to match the maximum flux from PSR B1259-63 (Sect. 4.4). For example, taking τesc = 10d/c implies a decreased injection fraction (PGeV ≈ 0.02Ė and PTeV ≈ 0.002Ė). The combination yields an estimated population of 105+92-54\hbox{$_{-54}^{+92}$} systems, i.e. very close to our previous estimate of 101+89-52\hbox{$_{-52}^{+89}$}. The peak of the gamma-ray flux distribution is pushed to longer periods than in Fig. 3, leading to a flatter distribution in the fraction of detected systems as a function of Porb. In principle, the Porb distribution of detected gamma-ray binaries could thus be used to constrain τesc, assuming excellent knowledge of their parent Porb distribution. A distribution of Lγ/Ė as a function of Porb would narrow down possibilities for the radiative and injection efficiency, which are presently only constrained by observations of PSR B1259-63. Gamma-ray observations of PSR J2032+4127 at periastron passage provide a second constraint on these efficiencies. Inversely, future observations of gamma-ray binaries as a population also have the potential to constrain the relative efficiencies in the GeV and TeV range, as described below (Sect. 5.3).

Many of our results remain applicable even if gamma-ray binaries are not powered by pulsar spin-down. The results of Sect. 3, based on the template light curves, do not depend on this assumption. The results of Sects. 4.1−4.3 are also applicable as long as the emission arises from electrons located close to the compact object (e.g. at the base of a jet) upscattering stellar radiation. Even if the compact object is a black hole, the compact object mass remains much lower than the companion mass, so any difference in orbits is minor for the radiation model. However, differences can be expected in the full population model (Sect. 4.4) since we made use of the distribution of spin-down powers of pulsars. We would need some assumptions on the distribution of jet power to perform an equivalent calculation and deduce the parent population. However, as stated in the introduction, we consider it very unlikely that gamma-ray binary emission arises from accretion-powered jets (Dubus 2013).

5.3. Future gamma-ray observations

We have aimed to discover what future observations hold in store. With an estimated population size of 101 gamma-ray binaries, up to 8 new binaries might be detected in an extended 3FGL survey with a most likely value of 2 new detections beyond the known sample. New discoveries are less likely in the VHE surveys. Once the expected detections (Table 1) are taken into account – up to 3 new detections are predicted in the HAWC survey, 5 in the initial CTA survey, and 6 in the full CTA survey, with 95% confidence limits – the most likely outcome statistically is no new detection. The reason is that the detection probabilities remain small for these surveys.

Serendipitous discoveries in deep VHE observations of Galactic sources (e.g. HESS J0632+057) can complement the surveys. We find that the probability for a chance detection of a gamma-ray binary is 0.17 ± 0.03% in a 100 h CTA exposure towards the Galactic centre, covering 6° in Galactic longitude, and reaching 1 mCrab at 1 TeV. This is ≈1.7 times the detection rate from the Galactic Plane survey over a comparable area, i.e. there are roughly 7 previously undetected systems for every 10 systems detected in the Galactic Plane survey of this deep field. Having 20 such deep pointings, spread around the Galactic plane towards areas of special interest such as the Galactic centre, the Cygnus and Westerlund regions, or the Sagittarius-Carina spiral arm (see Fig. 2), adds 1.4% to the detected fraction with CTA. Combining surveys and deep pointings can thus yield a detection rate that is comparable to or greater than that in the Fermi-LAT survey.

Any discovery in a VHE survey would have a major impact on the estimated population number, raising it to higher values. A discrepancy could appear between the actual number of sources detected in the VHE and HE surveys since the model predicts that essentially all TeV sources should be detected at GeV energies. Some tension is already present in the model. The maximum likelihood Lm obtained by treating the HE and VHE surveys independently is ≥ 20% for both, with corresponding population numbers of 82+108-56\hbox{$^{+108}_{-56}$} (HE) and 132+268-86\hbox{$^{+268}_{-86}$} (VHE). Combining the HE and VHE numbers into a single likelihood gives the estimate of 101+89-52\hbox{$_{-52}^{+89}$} systems presented above, but Lm drops to 4%. This low probability indicates that the model has difficulty accounting for both the number of GeV and TeV detections when they are taken from the same underlying binary population. This can be resolved by increasing the injection fraction Pinj at 1 TeV, with the effect of raising the detection probability in VHE surveys and lowering the parent population size to a value that slackens the tension with the HE constrains, or by lowering it at 1 GeV with opposite effects on detection probability and population size. Hence, the relative numbers of HE and VHE detections can constrain the relative injection efficiencies. In any case, regardless of the value of Pinj, the population of VHE-emitting systems is unlikely to be greater than 230 systems, otherwise HESS J0632+057-like systems would be detected in the 3FGL survey or in the H.E.S.S. Galactic Plane survey (Sect. 3.3). This number is close to the upper limit on the population size estimated from synthetic light curves (190 systems). Both estimates thus converge to a maximum gamma-ray binary population of ≈ 200 systems. With 200 systems, up to 10 (resp. 14) new binaries could be detected in the initial (resp. full) CTA survey, where the most likely number is 4 (resp. 6) discoveries.

6. Conclusions

We have modelled the population of gamma-ray binaries and evaluated the fraction of systems that can be detected in various HE and VHE surveys, taking into account the variability of their gamma-ray emission. The number of gamma-ray binaries is constrained to 101+89-52\hbox{$_{-52}^{+89}$} systems in our Galaxy. This number matches expectations from HMXB population synthesis.

Gamma-ray binaries are rare systems and we do not expect a watershed of discoveries in the near future. Pursuing the Fermi-LAT survey to ≈ 2024 should lead to a handful of discoveries of mostly PSR B1259-63-like systems. At very high energies, combining Galactic Plane surveys and deep observations of Galactic sources with CTA should provide a comparable number of discoveries. However, the number of HESS J0632+057-like systems with very weak GeV emission is a major source of uncertainty. Observations already indicate that the GeV and TeV emission originate from different particle populations. A VHE survey could therefore reveal a population of binaries that cannot be seen with the Fermi-LAT. Such a population is limited to ≲230 systems based on the lack of HESS J0632+057-like systems in the Fermi-LAT 3FGL survey and the H.E.S.S. Galactic Plane survey. With 200 systems, four new gamma-ray binaries can be expected in the first two years of the CTA Galactic Plane survey. Of course, these numbers refer only to gamma-ray binaries and do not limit gamma-ray detections from other types of binaries such as novae, colliding wind binaries, binary millisecond pulsars, and microquasars.

Detecting a system depends more on its orbit-averaged flux than on the shape of the gamma-ray light curve. Thus, the scheduling of visits from ground-based instruments plays a minor role in setting the detected fraction. The average flux is set by the efficiency with which spin-down power is radiated in the HE and VHE bands. This is the most important source of uncertainty in our model. Ideally, this should be constrained by measuring the pulsar spin-down power and radiated luminosity for as many systems as possible. At present, this is limited to PSR B1259-63with the possible addition of PSR J2032+4127 in the near future. Alternatively, this relative efficiency in the HE and VHE bands can be constrained statistically by the relative number of sources detected in HE and VHE surveys.

About 55% of pulsars in orbit around massive stars are hardly accessible to gamma-ray observations, which are most sensitive to the high Ė, short Porb systems. Low Ė and long Porb binaries are likely to be more efficiently accessed by radio pulsar surveys, which are thus fully complementary to the gamma-ray observations. Another significant fraction, ≈ 23%, may actually be visible as accreting X-ray pulsars or propellers instead of binary pulsar wind nebulae. Future work should strive to combine detection probabilities in gamma rays with detection probabilities in radio (SKA) and X-ray (eROSITA) surveys.

1

For the VHE surveys, we converted Crab units to integrated flux above 1 TeV using 1 Crab ≡ 1.82 × 10-11ph cm-2 s-1, based on the Crab spectrum measured by Albert et al. (2008).

3

The Fermi-LAT detection reported by Malyshev & Chernyakova (2016) is compatible with the detection of the low energy end of the VHE spectrum rather than the detection of a distinct GeV spectral component as in the other gamma-ray binaries.

4

Here, we estimate the probability p to detect a binary in a survey. We find m detections after running a random sample of n systems through our mock survey procedure. The likelihood function is L(p)=Cnmpm(1−p)n−m\hbox{${L}(p)=C_n^m p^m (1-p)^{n-m}$}, where Cnm\hbox{$C_n^m$} is the binomial coefficient. The function has a maximum Lm for p = m/n. Defining the test statistic as T = 2log (Lm/L) and applying Wilks’ theorem, the 95% confidence interval on p is calculated from the interval, where T ≤ 3.84 (the cut-off value in a χ2 distribution with one degree of freedom beyond which the probability ≤ 5%). The chosen number of systems n to simulate is a compromise between computational time and statistical error.

5

We take the detection probability p derived by the simulation and find the population n that maximizes the likelihood (see footnote 4) to detect m observed systems. Here, m = 1 for each type of gamma-ray binary.

7

PSR J2032+4127 is a pulsed Fermi-LAT source and coincident with an extended, persistent VHE source. This source does not yet show evidence for variable gamma-ray emission related to binary motion, as seen in PSR B1259-63 and the other gamma-ray binaries.

Acknowledgments

We are grateful to Masha Chernyakova, Markus Böttcher, and Jamie Holder for their careful CTA internal review that helped improve this work. G.D. and P.O.P. acknowledge support from Centre National d’Etudes Spatiales (CNES).

References

  1. Abdo, A. A., Ackermann, M., Ajello, M., et al. (Fermi-LAT Collaboration) 2009, ApJ, 706, L56 [NASA ADS] [CrossRef] [Google Scholar]
  2. Abdollahi, S., Ackermann, M., Ajello, M., et al. (Fermi-LAT Collaboration) 2017, ApJ, 846, 34 [NASA ADS] [CrossRef] [Google Scholar]
  3. Abeysekara, A. U., Albert, A., Alfaro, R., et al. (HAWC Collaboration) 2017, ApJ, 843, 40 [NASA ADS] [CrossRef] [Google Scholar]
  4. Abramowski, A., Aharonian, F., Ait Benkhali, F., et al. (H.E.S.S. Collaboration) 2015, A&A, 577, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Acciari, V. A., Aliu, E., Arlen, T., et al. (VERITAS Collaboration) 2011, ApJ, 738, 3 [NASA ADS] [CrossRef] [Google Scholar]
  6. Acero, F., Ackermann, M., Ajello, M., et al. (Fermi-LAT Collaboration) 2015, ApJS, 218, 23 [Google Scholar]
  7. Ackermann, M., Ajello, M., Ballet, J., et al. (Fermi-LAT Collaboration) 2012, Science, 335, 189 [NASA ADS] [CrossRef] [Google Scholar]
  8. Ackermann, M., Ajello, M., Ballet, J., et al. (Fermi-LAT Collaboration) 2013, ApJ, 773, L35 [NASA ADS] [CrossRef] [Google Scholar]
  9. Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. (H.E.S.S. Collaboration) 2006, A&A, 460, 743 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Albert, J., Aliu, E., Anderhub, H., et al. (MAGIC Collaboration) 2008, ApJ, 674, 1037 [NASA ADS] [CrossRef] [Google Scholar]
  11. Aliu, E., Archambault, S., Aune, T., et al. (VERITAS and H.E.S.S. Collaborations) 2014, ApJ, 780, 168 [NASA ADS] [CrossRef] [Google Scholar]
  12. Ambartsumian, V. A. 1937, Astron. Zhurn., 14, 207 [Google Scholar]
  13. An, H., & Romani, R. W. 2017, ApJ, 838, 145 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bassa, C. G., Brisken, W. F., Nelemans, G., et al. 2011, MNRAS, 412, L63 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bhadkamkar, H., & Ghosh, P. 2012, ApJ, 746, 22 [NASA ADS] [CrossRef] [Google Scholar]
  16. Bordas, P., Dubus, G., Eger, P., et al., 2017, in 6th International Meeting on High-Energy Gamma-Ray Astronomy, AIP Conf. Proc., 1792, 040017 [Google Scholar]
  17. Campana, S., Stella, L., Mereghetti, S., & Colpi, M. 1995, A&A, 297, 385 [NASA ADS] [Google Scholar]
  18. Carramiñana, A. 2016, J. Phys. Conf. Ser., 761, 012034 [CrossRef] [Google Scholar]
  19. Corbet, R. H. D., Chomiuk, L., Coe, M. J., et al. 2016, ApJ, 829, 105 [Google Scholar]
  20. Cusumano, G., Maccarone, M. C., Nicastro, L., Sacco, B., & Kaaret, P. 2000, ApJ, 528, L25 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  21. Dubus, G. 2006, A&A, 451, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Dubus, G. 2013, A&ARv, 21, 64 [NASA ADS] [CrossRef] [Google Scholar]
  23. Dubus, G., Cerutti, B., & Henri, G. 2010, A&A, 516, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Dubus, G., Lamberts, A., & Fromang, S. 2015, A&A, 581, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Gregory, P. C., & Taylor, A. R. 1978, Nature, 272, 704 [NASA ADS] [CrossRef] [Google Scholar]
  26. Grimm, H.-J., Gilfanov, M., & Sunyaev, R. 2002, A&A, 391, 923 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Hadasch, D., Torres, D. F., Tanaka, T., et al. 2012, ApJ, 749, 54 [Google Scholar]
  28. Hinton, J. A., Skilton, J. L., Funk, S., et al. 2009, ApJ, 690, L101 [Google Scholar]
  29. Ho, W. C. G., Ng, C.-Y., Lyne, A. G., et al. 2017, MNRAS, 464, 1211 [NASA ADS] [CrossRef] [Google Scholar]
  30. Iben, I. J., Tutukov, A. V., & Yungelson, L. R. 1995, ApJS, 100, 217 [NASA ADS] [CrossRef] [Google Scholar]
  31. Illarionov, A. F., & Sunyaev, R. A. 1975, A&A, 39, 185 [NASA ADS] [Google Scholar]
  32. Johnston, S., Manchester, R. N., Lyne, A. G., et al. 1992, ApJ, 387, L37 [NASA ADS] [CrossRef] [Google Scholar]
  33. Li, J., Torres, D. F., Cheng, K.-S., et al. 2017, ApJ, submitted arXiv:1707.04280] [Google Scholar]
  34. Lipunov, V. M., Nazin, S. N., Osminkin, E. Y., & Prokhorov, M. E. 1994, A&A, 282, 61 [NASA ADS] [Google Scholar]
  35. Lutovinov, A. A., Revnivtsev, M. G., Tsygankov, S. S., & Krivonos, R. A. 2013, MNRAS, 431, 327 [NASA ADS] [CrossRef] [Google Scholar]
  36. Lyne, A. G., Stappers, B. W., Keith, M. J., et al. 2015, MNRAS, 451, 581 [NASA ADS] [CrossRef] [Google Scholar]
  37. Madsen, E. C., Stairs, I. H., Kramer, M., et al. 2012, MNRAS, 425, 2378 [NASA ADS] [CrossRef] [Google Scholar]
  38. Malyshev, D., & Chernyakova, M. 2016, MNRAS, 463, 3074 [NASA ADS] [CrossRef] [Google Scholar]
  39. Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993 [NASA ADS] [CrossRef] [Google Scholar]
  40. Massi, M., Migliari, S., & Chernyakova, M. 2017, MNRAS, 468, 3689 [Google Scholar]
  41. Meurs, E. J. A., & van den Heuvel, E. P. J. 1989, A&A, 226, 88 [NASA ADS] [Google Scholar]
  42. Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15 [Google Scholar]
  43. Paredes, J. M., Martí, J., Ribó, M., & Massi, M. 2000, Science, 288, 2340 [Google Scholar]
  44. Portegies Zwart, S. F., & Verbunt, F. 1996, A&A, 309, 179 [NASA ADS] [Google Scholar]
  45. Portegies Zwart, S. F., & Yungelson, L. R. 1998, A&A, 332, 173 [NASA ADS] [Google Scholar]
  46. Ringermacher, H. I., & Mead, L. R. 2009, MNRAS, 397, 164 [NASA ADS] [CrossRef] [Google Scholar]
  47. Russeil, D. 2003, A&A, 397, 133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Shvartsman, V. F. 1971, Soviet Astron., 15, 342 [NASA ADS] [Google Scholar]
  49. Skinner, G. K., Bedford, D. K., Elsner, R. F., et al. 1982, Nature, 297, 568 [NASA ADS] [CrossRef] [Google Scholar]
  50. Stairs, I. H., Manchester, R. N., Lyne, A. G., et al. 2001, MNRAS, 325, 979 [NASA ADS] [CrossRef] [Google Scholar]
  51. Tauris, T. M., & van den Heuvel, E. P. J. 2006, in Compact Stellar X-ray Sources, eds. W. H. G. Lewin, & M. van der Klis (Cambridge, UK: Cambridge University Press), Cambridge Astrophys. Ser., 39, 623 [Google Scholar]
  52. Vercellone, S. 2017, in 6th International Meeting on High-Energy Gamma-Ray Astronomy, AIP Conf. Proc., 1792, 030001 [Google Scholar]
  53. Walter, R., Lutovinov, A. A., Bozzo, E., & Tsygankov, S. S. 2015, A&ARv, 23, 2 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Example synthetic light curves

thumbnail Fig. A.1

Example lightcurves computed from the model described in Sect. 4.2 (full line: GeV emission; dashed line: TeV emission, taking into account γγ absorption). The orbital period and eccentricity of the binary system is indicated in the title of each plot. The lightcurves are normalised to the maximum value. The systems shown here are a random selection of the systems flagged as detected in Fig. 6. Periastron passage is at phase 0.

All Tables

Table 1

Detected systems in the mock surveys based on the template light curves in Fig. 1.

In the text
Table 2

Fraction (%) of observable binaries in each survey.

In the text
Table 3

Fraction of detected systems in each survey using the light curves in Fig. 1 as templates (see Sect. 3.3).

In the text
Table 4

Fraction of detected systems in each survey using synthetic light curves as templates (see Sect. 4.3).

In the text
Table 5

Detection fractions for the population shown in Fig. 6.

In the text

All Figures

thumbnail Fig. 1

Template light curves for known gamma-ray binaries. Two orbits are shown except for PSR B1259-63 where the plot focusses on periastron passage. The various measurements are shown in grey, with error bars omitted for clarity. Arrows indicate measurement upper limits. Left: photon flux integrated above 1 GeV in units of 10-9 ph cm-2 s-1 based on Fermi-LAT measurements. Right: photon flux integrated above 1 TeV in units of 10-12 ph cm-2 s-1 based on IACT measurements. The GeV and TeV data are taken from Abdo et al. (2009) and Aharonian et al. (2006) for LS 5039,  where ΓGeV = 2.54; Hadasch et al. (2012) and Acciari et al. (2011) for LS I +61°303, where ΓGeV = 2.42 and ΓTeV = 2.6; Bordas et al. (2016) for PSR B1259-63, where ΓGeV = ΓTeV = 2.7; Ackermann et al. (2012) and Abramowski et al. (2015) for 1FGL J1018.6-5856, where ΓGeV = 3.1 and ΓTeV = 2.7. The GeV data for LMC P3 are from Corbet et al. (2016), where ΓGeV = 2.8. The TeV data for HESS J0632+057 are from Aliu et al. (2014).
In the text
thumbnail Fig. 2

Top: map of randomly generated locations for gamma-ray binaries in the Galaxy. Bottom: longitude distribution of gamma-ray binaries under the assumptions described in Sect. 3.2.
In the text
thumbnail Fig. 3

Mean orbit-averaged gamma-ray flux, normalized to the injected power, as a function of Porb (see Sect. 4.2).

In the text
thumbnail Fig. 4

Mean fractional amplitude of the simulated GeV (dashed line) and TeV (dash-dotted line) gamma-ray light curves (see Fig. 3) and fraction of the systems showing eclipses as a function of Porb (thick line).
In the text
thumbnail Fig. 5

Detected fractions in the HE (left panels) and VHE (right panels) surveys. Each panel contains four curves corresponding to Pinj = 1033,1034,1035,1036 erg/s (dark to light blue sequence in each panel). The detection fraction increases when going from small Pinj (darker blue) to high Pinj (lighter blue) until the detected fraction saturates at the value given by the observable fraction (Table 2).
In the text
thumbnail Fig. 6

Orbital period (left panels) and spin-down power distribution (right panels) of a random sample of 105 systems. Top panels show the full population as a light blue histogram with the line-filled fraction showing the distribution of the binaries that are accreting according to Eq. (4). The bottom panels zoom in to highlight the systems detected in any of the present-day surveys (combining the H.E.S.S., 3FGL, and FAVA-like surveys; medium blue histogram) or any of the future surveys (combining the full CTA, extended 3FGL, and HAWC-like surveys, dark blue histogram). The latter essentially shows the systems detected in any survey since a binary detected in one of the present-day surveys has a nearly >99% chance of being re-detected in one of the future surveys.
In the text
thumbnail Fig. A.1

Example lightcurves computed from the model described in Sect. 4.2 (full line: GeV emission; dashed line: TeV emission, taking into account γγ absorption). The orbital period and eccentricity of the binary system is indicated in the title of each plot. The lightcurves are normalised to the maximum value. The systems shown here are a random selection of the systems flagged as detected in Fig. 6. Periastron passage is at phase 0.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.