Constraining gammaray pulsar gap models with a simulated pulsar population
Key Words.:
stars: neutron, pulsars: general, gamma rays: stars, radiation mechanisms: non thermal, methods: numerical, surveysWith the large sample of young ray pulsars discovered by the Fermi Large Area Telescope (LAT), population synthesis has become a powerful tool for comparing their collective properties with model predictions. We synthesised a pulsar population based on a radio emission model and four ray gap models (Polar Cap, Slot Gap, Outer Gap, and One Pole Caustic). Applying ray and radio visibility criteria, we normalise the simulation to the number of detected radio pulsars by a select group of ten radio surveys.
The luminosity and the wide beams from the outer gaps can easily account for the number of Fermi detections in 2 years of observations. The wide slotgap beam requires an increase by a factor of of the predicted luminosity to produce a reasonable number of ray pulsars. Such large increases in the luminosity may be accommodated by implementing offset polar caps. The narrow polarcap beams contribute at most only a handful of LAT pulsars. Using standard distributions in birth location and pulsar spindown power (), we skew the initial magnetic field and period distributions in a an attempt to account for the high Fermi pulsars. While we compromise the agreement between simulated and detected distributions of radio pulsars, the simulations fail to reproduce the LAT findings: all models underpredict the number of LAT pulsars with high , and they cannot explain the high probability of detecting both the radio and ray beams at high . The beaming factor remains close to 1.0 over 4 decades in evolution for the slot gap whereas it significantly decreases with increasing age for the outer gaps. The evolution of the enhanced slotgap luminosity with is compatible with the large dispersion of ray luminosity seen in the LAT data. The stronger evolution predicted for the outer gap, which is linked to the polar cap heating by the return current, is apparently not supported by the LAT data. The LAT sample of ray pulsars therefore provides a fresh perspective on the early evolution of the luminosity and beam width of the ray emission from young pulsars, calling for thin and more luminous gaps.
1 Introduction
After the radio detection of the first pulsar signal in 1967 (Hewish et al., 1968), a pulsar magnetosphere model was formulated by Goldreich & Julian (1969). A direct consequence of the Goldreich & Julian model is the establishment of a magnetospheric charge density that creates a forcefree pulsar magnetosphere. However, such a magnetosphere has no electric field along the magnetic field to accelerate charges and produce rays.
The detection, a few years later, of pulsed emission at ray energies from the Crab (McBreen et al., 1973) and Vela (Thompson et al., 1975) pulsars, and the detection of four more ray pulsars by Thompson et al. (1994) established that pulsars accelerate particles to energies of at least a few TeV suggesting that there are magnetospheric regions where the charge density departs from that of Goldreich & Julian, locally violating the forcefree condition and allowing particle acceleration. These regions were identified in two magnetospheric zones. In the inner magnetosphere, acceleration can take place both above the polar cap and in the slot gap, which extends to highaltitude along the last open magnetic field lines. In the outer magnetosphere, the outer gap extends from the null charge surface to the light cylinder. These gap regions correspond to three models: the lowaltitude slotgap model, hereafter Polar Cap (PC, Muslimov & Harding (2003)), the Slot Gap model (SG, Muslimov & Harding (2004)), and the Outer Gap model (OG, Cheng et al. (2000)).
In the polarcap model the emission comes from a region close to the neutron star (NS) surface and well confined above the magnetic polar cap. Charged particles from the neutron star are initially accelerated in the strong electrostatic field generated by a departure from the GoldreichJulian charge density (Arons & Scharlemann, 1979). Aided by inertial frame dragging (Muslimov & Tsygan, 1992), pulsars emit high energy photons by curvature radiation (CR) and inverse Compton scattering (ICS). The most energetic of these photons reach threshold for electronpositron pair production in the strong magnetic field at a Pair Formation Front (PFF), above which the secondary pairs can screen the electric field in a short distance. The pairs, produced in excited Landau states, emit synchrotron photons which trigger a pair cascade with high multiplicity. A small fraction of the pairs is actually accelerated. The pair plasma likely establishes forcefree conditions along the magnetic field lines above the PFF, as well as radiate rays. Over most of the polar cap, the PFF and ray emission occurs well within a few stellar radii of the NS surface. The main contribution to the ray emission comes from CR from the pairs moving upward. Since the CR intensity scales with the magnetic field lines curvature, it decreases from the polar cap edge toward the magnetic axis, conferring to the emission beam the structure of an hollow cone.
The slotgap emission is generated from the same polar cap electromagnetic pair cascade near the boundary of the closed magnetic field lines region where the parallel electric field and the PFF rises to higher altitude. Here electrons are accelerated over longer distances to produce the pair cascade. A narrow gap, the slot gap, is formed along the closed magnetic field surface where the PFF is never established, and electrons continue to be accelerated and radiating rays by selflimited curvature radiation into the outer magnetosphere. The resulting hollow beam is much broader and less collimated near the magnetic axis than the loweraltitude PC emission (see Section 7).
The outer gaps are vacuum regions characterised by a strong electric field along the magnetic field lines (Holloway, 1973; Cheng et al., 1976) above the null charge surface. Two outer gap regions (Cheng et al., 1976; Romani & Yadigaroglu, 1995; Cheng et al., 2000; Hirotani, 2006) can exist in the angular velocitymagnetic momentum plane, one for each pole. In the physical OG model, in the case of a nonaligned rotator, the gap region closer to the pulsar surface is more active than the other gap further away from the surface due to the pair production screening operating more efficiently at lower altitude. In the OG model a chargedeficient region forms in the outer magnetosphere above the null charge surface where a chargeseparated flow is formed. The induced electric field accelerates pairs radiating rays in a direction tangent to the lines. The ray photons interact with thermal Xrays from the NS surface to produce pairs on field lines interior to the last open field line. The pair formation surface screening the electric field defines the interior surface of the gap.
More than 2000 pulsars are listed in the ATNF database (Manchester et al., 2005), most of which were first observed at radio wavelength. We employ the following ten selected pulsar radio surveys in this study: Molonglo2 (Manchester et al., 1978), Green Bank 2 & 3 (Dewey et al., 1985; Stokes et al., 1985), Parkes 2 (70 cm) (Lyne et al., 1998), Arecibo 2 & 3 (Stokes et al., 1986; Nice et al., 1995), Parkes 1 (Johnston et al., 1992), Jodrell Bank 2 (Clifton & Lyne, 1986), Parkes Multibeam (Manchester et al., 2001) and the extended Swinburne surveys (Edwards et al., 2001; Jacoby et al., 2009). For these, the survey parameters are known with a high accuracy and they cover the largest possible sky surface while minimising the overlapping regions.
The advent of the LAT telescope on the Fermi satellite (Atwood et al., 2009) led to a drastic increase in the number of ray pulsars. After three years of observations the LAT detected about 106 pulsars, more than doubling the number of detections listed in the first pulsar catalog (Abdo et al., 2010) leading to the discovery of two well defined ray pulsar populations consisting of 31 millisecond pulsars, and 75 young or middle aged isolated, normal pulsars. To study and compare the collective properties of the LAT normal isolated pulsars and investigate the emission mechanisms that best explain the observed emission, we synthesised a pulsar population incorporating four important highenergy radiation gap models. The simulation takes into account the axisymmetric structure of our Galaxy and is designed to match the known characteristics of the group of older radio pulsar population than the younger group of pulsars sampled in rays. Four ray emission gap models have been assumed: the previously described Polar Cap (PC), Slot Gap (SG), and Outer Gap (OG), and a variation of the OG, hereafter the One Pole Caustic (OPC) (Romani & Watters, 2010; Watters et al., 2009) that differs from the OG in the energetics. We model the radio emission at two different frequencies, 1400 MHz and 400 MHz (Gonthier et al., 2004; Harding et al., 2007), comparing simulated radio fluxes with the flux thresholds of existing surveys.
The outline of this paper is as follows. In Sections 2 and 3, we describe the neutron star characteristics and evolution. In Sections 4, 5, and 6, we give a brief overview of the radio luminosity computation, ray gap widths, and ray luminosities computations. Sections 7 and 8 describe the pulsar lightcurve and flux computation. Section 9 reviews the radio and ray pulsar visibility calculations. We present the results in the final Section 10.
2 Neutron star characteristics
The neutron star mass, radius, and moment of inertia used in this paper have been chosen according to the experimental mass measurements in binary NSNS systems, Xray binaries, and NSwhite dwarf binaries shown in Figure 3 of Lattimer & Prakash (2007).
The assumed NS mass and radius are and km. The mass value lies between the weighted average and average values of Xray and white dwarfNS binaries estimates and, with the km, represent, a possible solution for the EOS that describe the NS interior (Figure 2 of Lattimer & Prakash (2007)).
The moment of inertia of a NS is evaluated by Equation 35 of Lattimer & Prakash (2007). For the 13 km radius and the 1.5 mass of our standard NS, we obtain kg m. Because of the uncertainty on the mass and radius estimates, this value has an uncertainty of about 70%.
For each simulated NS we have generated a value of the magnetic obliquity (angle between the pulsar rotation and magnetic axes) and of the observer line of sight (angle between the pulsar rotation axis and the observer line of sight). After the supernova explosion that generates the neutron star, the magnetic axis has equal probability to point in any direction of a 3 dimensional space. This is also true for the observer line of sight direction with respect to the pulsar rotational axis. The and distributions are isotropic.
The spindown power is defined as the rate with which the pulsar loses rotational kinetic energy, as
(1) 
The latter equation is based on the NS structure assumptions through the moment of inertia . Since mass and radius are chosen inside intervals of allowed values, the estimate is affected by an uncertainty of at least a factor 3.
The choice of mass, radius, and moment of inertia formulation yields a moment of inertia value that is 1.5 times higher than in the ATNF catalog. This helps to reduce the discrepancy found between the simulated and observed distributions (Section 10.3), while remaining well within the range of parameters allowed by the binary data and equations of state in Lattimer & Prakash (2007). The choice of different values for mass and radius would also impact the range of the and distributions of the evolved pulsar population (Section 10.2).
3 Neutron stars at birth and their evolution
We synthesised NSs with mass, radius, and moment of inertia as described in Section 2, and assuming a constant birth rate over the last 1 Gyr. It yields isolated ordinary pulsars to the left of the radio death line (see below). In order to match the observed radio pulsars and distributions, an exponential magnetic field decay with a time scale of yr has been assumed (Gonthier et al., 2004). The choice of such a short timescale decay is justified by the need to slow down the birth population enough to reproduce the characteristics of the observed radio sample. It provides a simple mathematical solution to a more physical model of the rotational evolution of the NS, yet to be developed. For our study, since we are dealing with young ordinary pulsars, this choice has been checked not to affect the obtained results.
The radio death line we used is defined as
(2) 
It is composed by three different segments (Story et al., 2007; Zhang et al., 2000), each one refers to a specific period interval characterised by the following and values
(3) 
3.1 Birth spinning and magnetic characteristics
The distribution of period at birth, , plotted in the right panel of Figure 1, follows a single gaussian of width 50 ms, centred at 50 ms, and truncated at 0 to avoid negative periods. The same distribution was adopted by Watters & Romani (2011) on the basis of radio luminosity arguments but it differs from the choice of Takata et al. (2011) who selected the birth period randomly in the range ms.
The magnetic field birth distribution shown in the left panel of Figure 1 has been built as the sum of two gaussians in [Tesla], both 0.4 in width, respectively centred at 8.5 and 9.1, and with an amplitude ratio of 1:7/12. Our choice represents a compromise between that of Watters & Romani (2011), a single gaussian centred at 8.65 and width 0.3, and the Takata et al. (2011) one, a single Gaussian centred at 8.6 and width 0.1. The high Gaussians provide energetic pulsars when evolved.
Both the and distributions have been optimized a posteriori to obtain, after evolution, a simulated pulsar sample as close as possible to the observed one by minimizing the observed lack of high objects (Section 10.3). The birth distribution has been derived from and by using the equation
(4) 
This formulation includes no dependence on the magnetic obliquity alpha, as proposed by Ruderman & Sutherland (1975) for the spinning down of magnetospheres carrying current flows. More recently, Spitkovsky (2006) numerically showed for forcefree magnetospheres that the spin down of orthogonal rotators is twice that of aligned rotators. In the nonideal case of a magnetosphere accelerating charges to produce pulsed emission the impact of alpha on is still under discussion, so we chose for this paper the alpha independent prescription of Ruderman & Sutherland (1975). Hereinafter, all luminosities are given as a function of to judge how the uncertainty on the spindown rate propagates.
3.2 Birth location and velocity in the Galactic plane
To follow the dynamical evolution of the pulsars in the Galactic reference frame, we synthesised their birth position , , in the Galaxy as well as their kick velocity and direction.
We emulated the distribution of the NS progenitors by using the location of the HII regions in the Galaxy. The latter are good tracers of massive stars because OB stars are required to ionise the hydrogen bubbles. For the number density of pulsars at birth as a function of Galactocentric distance, we used the HII region profile recently obtained by Bania et al. (2010) from radio observations that can probe HII regions to large distance with little absorption. Figure 2 shows the comparison between the Paczyński (1990) birth distribution used in earlier publications (Gonthier et al. (2004); Takata et al. (2011)) and the HII region profile used here. Both distributions extend from the Galactic centre up to 40 kpc and have been normalised to 1 over the Galaxy.
We assume that all the NSs are born in the Galactic disk, with an exponential thin disk distribution with a scale height of 50 pc (consistent with Watters & Romani (2011) that adopted an exponential thin disk with a 75 pc scale height) and with a surface density distribution defined in Figure 2. Due to the large supernova kick velocity, the neutron stars evolve quickly out of the plane of the Galaxy. The assumed kick velocity distribution is the same as in Watters & Romani (2011) and Takata et al. (2011). It is described by a Maxwellian distribution, characterised by a mean of 400 km s and a width of 256 km s (Hobbs et al., 2005).
3.3 Evolution
We have evolved both the pulsar position and velocity in the Galactic gravitational potential (described in Paczyński (1990) and Gonthier et al. (2002) Equations 17, 18, and 19, and Takata et al. (2011)). The spin characteristics have been evolved to the present time assuming a magnetic dipole.
The simulated pulsar population at birth is shown, in red, in the  diagram of Figure 3. Following Gonthier et al. (2002), by knowing the analytical expression for , it is possible to follow the evolution of the spin parameters from the birth time to the present time . The magnetic decay is described by
(5) 
where = 2.8 Myr is the decay timescale, and is the birth magnetic field in units of Tesla.
4 Radio emission model
After evolving the neutron stars in the Galactic frame, values of the radio dispersion measure (DM), and the radio scattering measure (SM), are assigned to each star using the NE2001 model (Cordes & Lazio, 2001). The sky temperature at 408 MHz () for each star is obtained using the allsky map from the study of Haslam et al. (1982).
The empirical radio emission model we have implemented in our simulations follows the work of Gonthier et al. (2004) and Harding et al. (2007). We assume that the radio beam is composed of a core component originating relatively near the neutron star surface and a conical component radiated at higher altitude, both centered on the magnetic axis in the corotating frame. The adopted form of this model is similar to that proposed by Arzoumanian et al. (2002), based on the work of Rankin (1983) and Kijak & Gil (2003) and modified to include frequency dependence by Gonthier et al. (2004). The total flux at a given frequency from the two components seen at angle to the magnetic field axis is
(6) 
where
(7) 
The index refers to the core or cone, is the spectral index of the total angle and frequency integrated flux for each component, is the component luminosity, and is the distance to the pulsar. The total solid angles of the Gaussian beams describing the core and cone components are
(8)  
(9) 
where the latter can be written as
(10) 
where is the frequency expressed in GigaHertz. The factor represents the portion of that is independent of the frequency and used later in Equation (14). The width of the Gaussian describing the core beam is
(11) 
where is the pulsar period in seconds. The annulus and width of the cone beam are
(12) 
(13) 
where (Gonthier et al., 2006), and
(14) 
is the radius of the open field volume at the emission altitude derived by Kijak & Gil (2003), and
(15) 
is in units of stellar radius. The ratio of the coretocone peak flux is expressed as
(16) 
and requires , where . Gonthier et al. (2006), who carried out a study of 20 pulsars having three peaks in their averagepulse profiles at frequencies 400, 600 and 1400 MHz, found a coretocone peakflux ratio
(17) 
It is consistent with the ratio of Arzoumanian et al. (2002) at periods above about 1 s, and predicts that pulsars with s are cone dominated. Such a picture is supported by the study of Crawford et al. (2001) who measured the polarisation of a number of pulsars younger than 100 kyr, finding that they possess a high degree of linear polarisation and very small circular polarisation, typical of cone beams. The luminosities of the core and cone components are
(18) 
where
(19) 
is evaluated from Equations 16 and 17, , , MHz and MHz, and
(20) 
as modified from Arzoumanian et al. (2002).
5 Pc & Sg: particle luminosity and gap width
5.1 Particle luminosity
The slot gap region is defined between the last open magnetic field line, defined by the colatitude , and the magnetic field line with a colatitude value where is the SG width expressed in units of the dimensionless colatitude of a PC magnetic field line, .
It is possible to define the emission component from the PC pair cascades along the PFF that forms on the inside surface of the SG by assuming that monoenergetic radiation is emitted tangent to field lines (Muslimov & Harding, 2003). is a function of pulsar period, , and surface magnetic field, (Muslimov & Harding, 2003). The photons from the polar cap pair cascade are emitted in the region defined by . The luminosity of the SG from each pole is
(21) 
where and are the primary charge density and potential as a function of the emission altitude and of , in units of NS radius, and is the magnetic azimuthal angle. Using the expressions for and for from Muslimov & Harding (2003), the PC particle luminosity (from the lowaltitude SG) is
(22) 
where is the spindown power, , is the NS moment of inertia in unit of 10 kg m, is the NS radius in unit of m, is a relativistic correction factor of order 1, is the correction factor for the dipole component of the magnetic field in a Schwarzschild metric, and is the pulsar obliquity (Muslimov & Tsygan, 1992; Harding & Muslimov, 1998).
Using the equations for and for from Muslimov & Harding (2004), the highaltitude SG particle luminosity from each pole can also be determined from Equation (21) as
(23) 
where and . The parameters and , are defined as:
where , , and . According to Muslimov & Harding (2004), the energies of the primary electrons in the SG quickly become radiationreaction limited, with the rate of acceleration balancing the curvature radiation loss rate, resulting in 100% efficiency with in this case.
5.2 Gap width
In the SG model, the width of the slot gap can be estimated as the magnetic colatitude where the variation in height of the curvature radiation PFF (in units of stellar radius) becomes comparable to a fraction of the stellar radius (Muslimov & Harding, 2003):
(24) 
In Equation 24, represents the dimensionless altitude, above the polar cap, of pair formation due to curvature radiation
(25) 
where s, and is the magnetic field in units of 10 Tesla. By solving numerically Equation 24 with defined in Equation 25, one obtains for a specific pulsar. The gap width value is then obtained as
(26) 
The parameter constrains both the energetics and emission pattern of the SG emission and impacts both the SG and PC luminosity (Section 10.6) and lightcurve sharpness and shape. For large values the lightcurve peaks appear too sharp compared with the observed LAT profiles, therefore the slot gap is too narrow and not energetic enough to explain the observed LAT fluxes. On the other hand, smaller values imply wider slot gaps, sufficiently luminous when compared with the observations, but lightcurve peaks too broad when compared with the observed ones.
As a result we compromise between the narrow lightcurve structures and the ray luminosity through a reasonable radiation efficiency . We tried two different approaches to constrain : one based on energetic arguments, and one based on the optimisation of the expected lightcurves for some of the LAT pulsars.
Since scales as and since we want the luminosity to be close to (Abdo et al. (2010), First pulsar catalog) we need to have to obtain a reasonable agreement with the LAT data. The luminosity remains close to for all the tested values, but favours to explain the bright LAT pulsars. A good compromise is found for .
One can calculate numerically the pair formation front shape for the P and B values of some of the best known pulsars, Crab, Vela, CTA1, and Geminga, to obtain an approximate value (Muslimov & Harding (2003)). The results yield =0.03, =0.1, =0.16, =0.3 for values between 0.020.6.
In order to investigate how the pulsar lightcurve changes as a function of , we performed a fit to some LAT lightcurves with the SG phaseplots (see Section 7), evaluated for a set of values obtained for different values. We studied the behaviour of the bestfit likelihood value as a function of for Vela, Crab, J10285820, J10485832, J2021+3651, and J2229+6114. For all the studied pulsars, in the range that allows bright enough pulsars, the maximumlikelihood value presents a local maximum between 0.2 and 0.4. This result is consistent with the estimate obtained from the luminosity study and the pair formation front evaluation from Crab, Vela, CTA1, and Geminga.
In this paper, we set =0.35. This value reproduces the bulk of the lightcurve structure of the observed objects and yields a reasonable estimate of the SG luminosity. In choosing , we put more emphasis on matching the sharply peaked lightcurves often recorded by the LAT than on achieving bright luminosities. This selection of was driven by the need to preserve realistic beam patterns (thus their brightness and visibility across the beam) and is a key assumption that contains the results of our population studies. We mitigated the low SG gammaray luminosities by using a radiative efficiency greater than 1 as discussed in section 8.2.
6 OG and OPC: particle luminosity and gap width
6.1 Gap width
To determine the gap width, we consider two different prescriptions. The first one (Watters et al., 2009) simply assumes that the gap width is equal to the ray radiation efficiency. Because of the relation observed in the first LAT pulsar catalog (Abdo et al. (2010)) , the gap width should follow as
(27) 
Our second prescription follows the calculations of the self sustaining OG model presented in Zhang et al. (2004). In this formulation, the Xrays that trigger the pair production come from the bombardment of the NS surface by the full return current from the OG. The bright Xray luminosity allows active OGs and ray emission for many old pulsars. The outer gap width across magnetic field lines is determined by computing the location of the pair formation surface. From Kapoor & Shukre (1998), the polar angle corresponding to the magnetic field line tangent to the light cylinder is:
(28) 
with the light cylinder radius given by
(29) 
Here is the distance between the pulsar and the point where the light cylinder is tangent to the magnetic field line corresponding to . The lower boundary of the outer gap is estimated from the nullcharge surface, , that in two dimensions is described by . By definition, the polar angle at the inner edge of the outer gap is
(30) 
The computation of is obtained from the relation
(31) 
which results in
(32) 
The relation that defines the fractional OG size in this case is:
(33) 
where is a factor that is numerically solved for each pulsar by taking into account the average distance at which primary rays are produced and along which magnetic field line they pair produce when they interact with an Xray coming radially from the NS surface. The average distance is defined in Zhang et al. (2004) as
(34) 
where and is the radius at which the fractional size of the outer gap stops to grow: )=1.
A full calculation of the width of the OG radiating layer is complicated (Hirotani, 2006, 2008) since both the screening and the radiation occur in the same location. For this paper, we assume that this is an infinitely thin layer on the gap inner edge and that it is uniform in azimuth around the magnetic axis whereas Hirotani (2006) finds a significant azimuthal dependence.
6.2 Particle luminosities
The assumed gap width defined in section 6.1 is not based on any physical prescription and is very different from the usual dependence luminosity (gap width) (both SG and OG) based on the electrodynamics.
The gap width luminosity is evaluated as
(35) 
In the OG case, from Zhang et al. (2004) and previous papers dealing with OG gap geometry, the total ray luminosity is
(36) 
where is the fractional width of the gap at the average gap radius .
7 Phaseplot and lightcurve generation
7.1 Assumptions and photon distributions
To provide the ray emission pattern for each emission mechanism, we used the geometric emission model from Dyks et al. (2004) based on the following assumptions: (i) the pulsar magnetic field is dipolar and sweptback by the pulsar rotation (retarded potentials) (Deutsch, 1955), (ii) the ray emission is tangent to the magnetic field line and oriented in the direction of the accelerated electron velocity in the star frame. Relativistic aberration and time of flight delays are taken into account.
In the computation of the emission pattern, the first step consists in localising the position of the magnetic field line from which the radiation is produced. Each field line is then divided into segments and for each segment the tangent direction and height with respect to the NS surface is evaluated. Since the emission gap is located, for each model, in a different magnetospheric region, the emission patterns are obtained by selecting the segments corresponding to the gap position in each model. The ray emission is assumed to be uniform along the field lines in the corotating frame. The phase of the pulsar emission is defined by the direction of the emitted photons with respect to the corotating frame. The result of this computation is the twodimensional emission pattern in the plane (,), shown for each implemented emission model in Figure 4, which we refer to as a phaseplot. Figure 4 also shows the evolution of the emission pattern as a pulsar ages.
To incorporate the radio emission geometry we modulate the field lines with the flux given by Equation 6. The differential flux radiated from a bundle of field lines centred at openvolume coordinates (, ) (see Dyks et al. (2004)) is
(37) 
where is the number of azimuthal divisions of each polar cap ring, and are the lower and upper boundary of the emission region, and (= 0.1 for the radio phaseplot) is the spacing of the rings on the PC in open volume coordinates. For the SG model is adjusted to have 20 rings within the gap. The flux is assumed to be emitted at an altitude of for the core component and at an altitude given by Equation 15 for the cone component.
In the PC model, the emission profile in colatitude is infinitely thin along the inner edge of the slot gap (with defined in section 5.2) and the intensity of the emission along the field line, , exponentially decreases from the polar cap edges to the magnetic pole
(38) 
Here , and is the curvilinear distance along the field line starting from the NS surface. Both and are in unit of .
To model the emission component from primary electrons in the SG model, we assume that radiation is emitted along the field lines in the slot gap, up to altitude (where ). We assume an emissivity distribution across the SG as:
(39) 
where at the center of the SG and
(40) 
Such a distribution, that peaks in the centre of the gap and decreases to zero at the gap edges, follows from the distribution of the SG potential (Muslimov & Harding, 2004).
For the OG/OPC model, we describe the emitting region as an infinitely thin layer along the inner surface of the gap. The radio and PC phaseplots show the hollow cone patterns centered on the magnetic pole, while the SG and OG phaseplots show the caustic emission patterns characteristic of outer magnetosphere emission (Dyks et al., 2004).
7.2 Lightcurve generation
The (,) phaseplot space has been sampled in bins.
Each bin of the phaseplot gives the number of photons per solid angle per primary particle that can be observed in the direction at the rotational phase :
(41) 
Each phaseplot is obtained for a specific set of pulsar parameters that define its magnetospheric structure: the spin period , the surface magnetic field , and the magnetic obliquity . For the studied models, the phaseplot has the following dependencies:
For each emission model, we have evaluated phaseplots for values, from 5 to 90, with a step of . For each value, the phaseplots have been evaluated for 2 magnetic field values and 9 spin period values for the PC and radio models, and for 16 gap width values in the SG and OG/OPC cases. The complete set of sampled parameters is listed in Table 1.
To obtain the lightcurve of a given NS, with a particular set of , , , and gap width parameters, we interpolated the phaseplots noted in Table 1. When comparing phase profiles for a different set of parameters, typically one profile will be narrower than another one making it nontrivial to interpolate between them. We adopted a nonlinear interpolation which expands the narrower lightcurve covering the smallest phase range up to the phase extent of the wider profile, then applies a linear interpolation, and contracting the expanded profile back down to the extent of the original parent profiles. This strategy preserves the thin peaks and high degree of modulation that characterises the pulsar emission profiles at radio and ray wavebands.
B  P  Gap Width values  

Tesla  milliseconds  Degrees  
PC  ,  30, 40, 50, 75, 100  590  0.04, 0.06, 0.08, 0.1, 0.13, 0.16, 0.2, 0.225 
300, 500, 750, 1000  5  0.25, 0.275, 0.3, 0.34, 0.38, 0.42, 0.46, 0.50  
SG  590  0.04, 0.06, 0.08, 0.1, 0.13, 0.16, 0.2, 0.225  
none  none  5  0.25, 0.275, 0.3, 0.34, 0.38, 0.42, 0.46, 0.50  
OG/OPC  590  0.01, 0.025, 0.04, 0.05, 0.067, 0.084, 0.1, 0.2  
none  none  5  0.3, 0.4, 0.5, 0.53, 0.56, 0.59, 0.62, 0.65  
Radio  ,  30, 40, 50, 75, 100  590  
300, 500, 750, 1000  5  none 
8 Flux calculations
8.1 Phaseplot normalisation and energy flux
Let us define as the radiative luminosity from each pole, either in the rays or in the radio. Assuming a value for the primary particle production rate and the energy of each photon, one obtains a radiation luminosity per phaseplot bin:
(42) 
where is a proportionality constant. One can normalise the phaseplot to the total radiation luminosity over the two poles according to:
(43) 
We define the specific intensity as
(44) 
It is now possible to obtain the average energy flux observed by an Earth observer for a line of sight :
(45) 
Here, is the pulsar distance.
From the Equations 44 and 45, we can write the average energy flux observed at the Earth as:
(46) 
The latter Equation establishes the relation between the luminosity derived in the framework of a given model, , and the integral of the pulsar lightcurve , obtained, from the phaseplot, for . This is related to the beaming factor discussed in Section 10.5
8.2 Computations: gap width and energy flux
We calculated the ray and radio lightcurve for each pulsar of the sample, storing the value of the integral for the flux computation. The width of the emission gaps is computed using Equation 26 for the SG, and equations 33 and 27 for the OG and OPC. Because the PC and SG models do not apply when the gap becomes too large (pairstarved gaps should then be used), the flux for gap widths larger than 0.5 has been set to 0. Because no emission remains visible from the thin inner edge of OG/OPC gaps when the gap width exceeds 0.7 the flux for gap widths larger than 0.7 has been set to 0. So all the pulsars with a gap width above these threshold levels are assumed to not produce any ray emission.
For the radio luminosity computation, , we have used Equation 20 to evaluate the total radio luminosity and Equation 18 to evaluate the luminosities of each core and cone component. The ray luminosity has been obtained by scaling the particle luminosities derived in Equations 22, 23, 36, and 35, for the PC, SG, OG and OPC models respectively by using a radiative eficiency . The latter has been chosen to provide a good agreement between the observed and simulated distributions as a function of characteristic age ( is the photon flux and the pulsar distance). This distribution involves only readily observable quantities. The solution adopted for each model is shown in Figure 5. The choice of radiative efficiencies are: , , , and . The high value of needed for the SG requires either a super GoldreichJulian current or a stronger value of the accelerating electric field in the gap compared to the original calculation by Muslimov & Harding (2004). This is quite possible if the polar cap is slightly offset, i.e., nonsymmetrical around the magnetic axis, as one expects from the shape of the magnetic field lines distorted by the stellar rotation. Harding & Muslimov (2011) show that this distortion leads to a larger pair multiplicity as well as an increased electrical field along the field lines, thus an enhanced luminosity. Offset polar caps can sustain the modest increase in particle energy that is required in the present population study to account for the flux and pulsar counts observed by the LAT without invoking a radiation efficiency larger than one. The offset polar cap prediction was not available at the time of the population synthesis work, so we keep here the original polar cap luminosity and .
8.3 Gammaray energy and photon flux
In order to evaluate the photon flux from the energy flux and to compare it with the LAT sensitivity in photon flux, we need to assume an emission spectrum for the pulsars. The typical photon spectrum of a LAT pulsar is well fitted by a powerlaw with an exponential cutoff, like
(47) 
In the first LAT pulsar catalog (Abdo et al., 2010), the distribution of the parameters and can be described by two Gaussians:
(48) 
The spectral index and the are defined as:
(49) 
The Gaussian widths, centroids, and correlation angle had been derived from the analysis of the spectral parameters measured for the 1st LAT pulsar catalogue. We took here the same values.
9 Gammaray and radio visibilities
9.1 ray pulsar visibility
In order to select the simulated pulsars that could be detected by the LAT during two years of observation, we made use of the 6 month pulsar visibility map published for the 1st LAT pulsar catalogue (Abdo et al., 2010) and of the 1 year pulsar visibility of blind pulsar searches (Dormody et al., 2011). The two maps have been used to estimate the ray detectability of the radioloud pulsars (corresponding to the LAT radio selected objects) and the radioquiet ones (corresponding to the LAT blind search objects) respectively. The maps give the minimum visible photon flux and have been obtained taking into account the real LAT observation time in the sky, the photon energy, and the effective collection area corrected for the different incidence directions. Since the sky survey mode for LAT observations has been continued after 6 months, the maps have been scaled to 2 years as the square root of time for the radioselected sensitivity map and linearly with time for the blind search sensitivity map. Very few pointed observations were programmed that would significantly alter the shape of the visibility map. Photons collected in survey mode largely dominate and the flux threshold for detectability is primarily limited by the intensity of the interstellar background.
9.2 Radio visibility
The synthesis of the population is not based on any assumed NS birth rate; we assume a flat star formation rate over the last 1 Gyr. Instead, the simulated sample was scaled to the real number of pulsars detected in the Galaxy. The scaling factor has been evaluated by selecting all the ATNF radio pulsars present in a select group of ten surveys and comparing this number with simulated radio pulsars visible in the same region. We have generated a large enough population to reduce the Poisson fluctuations and to improve the statistics in the analysis results.
(K Jy)  (K)  (MHz)  (s)  (ms)  (MHz)  (MHz)  

Molonglo 2  5.100  5.4  225  408  40.96  40.0  4.0  4.0 
Green Bank 2  0.886  7.5  30.0  390  136  33.5  16.0  2.0 
Green Bank 3  0.950  8.0  30.0  390  132  2.0  8.0  0.25 
Parkes 2  0.430  8.0  50.0  436  157.3  0.6  32.0  0.125 
Arecibo 2  10.91  8.0  100  430  39.3  0.4  0.96  0.06 
Arecibo 3  13.35  8.5  70.0  430  67.7  0.5  10.0  0.078 
Parkes 1  0.256  8.0  45.0  1520  157.3  2.4  320.0  5.0 
Jodrell Bank 2  0.400  6.0  40.0  1400  524.0  4.0  40.0  5.0 
Parkes MB  0.460  8.0  21.0  1374  2100.0  0.25  288.0  3.0 
Swinburne  0.427  10.0  21.0  1374  265.0  0.25  288.0  3.0 

Computed using .

Computed using

Computed using .

Computed using .
We selected the simulated pulsars within the visibility criteria of 10 radio surveys from the ATNF
database
(50) 
This is the factor we used to scale the simulated pulsar sample.
Radio pulsar selection
During a radio survey, the edges of the survey region are defined by the position of the radiotelescope beam centre. Nevertheless, because of the solid angle extension and complexity of the beam, it is possible to observe a pulsar slightly out of the declared survey region. Thus, to say that all the pulsars observed during a survey fall inside the declared survey coordinates edges is not totally correct. The first parameter we reevaluated for each survey is the number of pulsars seen inside a given region.
Molonglo 2          85.0  20.0  0.62  0.4  0.03 

Green Bank 2          18.0  90.0  0.32  0.7  0.03 
Green Bank 3  15.0  230  15  15      0.41  0.75  0.03 
Parkes2          90.0  0.00  0.90  0.75  0.03 
Arecibo 2  40  66  10  10  9.50  25.0  0.54  1.0  0.05 
Arecibo 3  38  66  8.1  8.2  5.00  26.5  0.66  0.7  0.05 
Parkes 1  92  20  4  4      0.41  0.6  0.03 
Jodrell Bank 2  5  105  1.3  1.3      0.50  0.8  0.03 
Parkes MB  105  52  6.03  6.35      0.98  0.9  0.05 
Swinburne  100  50  4.5  30      0.87  1.0  0.05 
The second important parameter is the survey efficiency . It is defined as a filling factor, e.g. the ratio between the actual solid angle covered by the radio telescope beam during the observations, and the area within the declared survey boundaries. The survey efficiency can be considered as the probability of observing a pulsar present in the survey region only if the parent spatial distribution is uniform. To evaluate the boundaries of the survey region and to define the survey efficiency we decided:

to slightly extend the sky survey boundaries in order to include the largest number of pulsars actually detected by a survey, without changing too much the original boundaries

to evaluate the detection flux threshold for each pulsar within a survey by scaling the Dewey formula (Dewey et al., 1985) with a free parameter, , to match the observations.
(51) where the threshold flux is expressed by the Dewey formula
(52)
The Dewey formula, or radiometer formula, takes into account the characteristics of a given radio telescope and detector as well as a pulsar period and direction to give the minimum flux the survey would be able to detect. In Equation 52 is the minimum signal to noise ratio taken into account, is the receiver temperature, is the sky temperature at 408 MHz, is the ratio between the radio telescope gain and the dimensionless factor that accounts for system losses, is the number of measured polarisations, is the total receiver bandwidth, is the integration time, and is the pulsar period. is the effective pulse broadening, defined as
(53) 
Here, is the intrinsic pulse width (Duty Cycle), is a lowpass filter time constant applied before sampling (when this parameter is unknown, a value equal to twice the sampling time has been used), is the pulse smearing due to the DM over one frequency interval , and is the pulse broadening due to interstellar scattering (Dewey et al., 1985). The dispersion broadening time, (ms), across one frequency channel, , is related to the dispersion measure (DM) as
(54) 
where is the mass of the electron, is the speed of the light, and , are the edges of the frequency channel. The dispersion measure, DM (pc cm), is obtained using the Cordes & Lazio (2001) NE2001 model. The same model provides the scattering measure, SM (kpc m), which allows to estimate the broadening time due to interstellar scattering as
(55) 
where is the pulsar distance in kpc (Johnston et al., 1992; Sturner & Dermer, 1996). The last term of equation 53, , is an additional time broadening added when the sampling is performed for a DM value different from the real one. It corresponds to the fourth term of equation 2 in Dewey et al. (1985) and becomes important just for low period pulsars.
The sky temperature at frequencies other than 408 MHz is obtained as:
(56) 
Tables 2 and 3 list all the radio telescope and detector characteristics of the surveys we took into account. Some of the survey parameters in the literature listed as average values have been reevaluated using the above mentioned prescription.
The scaling of the radiometer equation was motivated by the uncertainties related to the Dewey formula, because of flux oscillations due to scintillation. The scintillation is caused by the turbulent variation of the interstellar medium that the pulsar light has to cross before reaching the observer. The consequence is an oscillation (scintillation) of the pulsar flux that can introduce spurious detections of pulsars with a flux lower than the survey threshold or that can cause the nondetection of pulsars with a flux higher than the survey threshold. So we scaled the level in order to take into account possible spurious detections or missed detections due to scintillation. should not be lower than the flux of the weakest pulsar of the survey. A reasonable estimate is to employ the average of the lowflux tail of the pulsars of the survey.
In the ATNF database we can count how many pulsars fall within a survey boundary, how many would match the survey visibility criterion (flux ), and how many of these pulsars have really been observed by the survey. The comparison of the ratios of the radio flux recorded for each pulsar to the minimum visible flux in its direction provides an estimate of the Dewey scaling factor. The scaling values are given in Table 3 and the distribution of the flux ratios is shown in Figures 19, 20, and 21 (right plots) for each survey (only the ratios below 10 are displayed to focus near the visibility threshold).
Then, for each survey, we derived the ratio between the number of pulsars really detected by the survey and the total number of observable ATNF pulsars (the sum of the detected ones plus those that match the position and flux survey criteria but were not detected). We consider this last ratio as the new survey efficiency, , i.e. the percentage of pulsars detected by the survey with respect to all the detectable ATNF pulsars in the survey region. The new efficiency is listed, for each Survey, in Table 3.
By using the newly estimated survey parameters listed in Tables 2 and 3, and by using the radiometer equation 52, the number of real pulsars that meet the visibility criteria of our surveys is 1430 (ATNF database, January 2012). We use this number and the number of simulated pulsars that match the same criteria to scale the visible component of the simulated ray pulsar population in Equation 50.
PC  1431  3.6  0.5  0.87 
SG  1508  75  79  0.49 
OG  1496  123  66  0.65 
OPC  1524  107  94  0.53 
LAT  25  30  0.45 
10 Results
10.1 Detection statistics
Table 4 indicates, for each model, the numbers of NSs that passed the radio and/or visibility criteria and their comparison with the LAT detections after 2 years of observations. The number of radio visible pulsars in the simulation has been scaled to the 1430 ATNF radio pulsars that passed the same selection criteria.
The scale factor of 0.136 is required to match the simulated and observed radio samples and has been applied to all star counts quoted hereafter, in particular to the ray simulated samples. This scale factor implies a NS birth rate of NS/century over the last 1 Gyr. The choice of radiative efficiencies driven by a reasonable agreement in the evolutions with characteristic age shows that the wide beams produced in the intermediatehigh (SG) and outer models provide enough detections to account the LAT findings. The lowluminosity narrow PC beam fails in predicting the LAT detection number and the fraction of radioquiet objects because of the large overlap between the ray and radio beams.
10.2 Comparison of the total simulated and observed samples
Figures 6 and 7 show the comparison between the simulated distributions in the diagram for the radio visible component and for the ray visible population for each model. The simulated distributions reasonably describe the observed samples and are in nice agreement with the same distributions obtained by Takata et al. (2011). The simulated radio population is able to describe the observed distribution for the fastest rotators that are likely to sustain substantial ray emission and represent the LAT pulsar population.
The PC model reproduces poorly the observed population. Both the SG and OG models over predict the number of middle aged ray pulsars and under predict the number of young ray objects. Of those, the OG shows the poorer description of the data; the core of the distribution is too close to the pulsar death line and it lacks energetic pulsars. The OPC visible population best describes the observed and of the LAT population with a core centred on the observed objects and tails that cover the overall dispersion.
Figure 8 compares the total simulated populations and its ray subsample to the observed total sample of radio and/or ray visible objects for key characteristics: period, period first time derivative, surface magnetic field, and distance.
The simulated spin period distribution is too broad to describe the observed proportion between the number of intermediate period objects ( ms) and the wings of the distribution. The range of spin periods is well covered and well centred, but we lack simulated objects in the 0.31.0 second range. The simulated distributions in , , and are all shifted to an excess of young, energetic, and nearby pulsars compared to the observed ones. This results from the choice of birth characteristics and NS intrinsic characteristics (, , and formulation) that emphasised nearby and high objects while preserving the bulk of the radio distributions. This choice has been made a posteriori to minimise the lack of high objects discussed in section 10.3. Nevertheless, the discrepancies observed in Figures 6 and 8 are not only due to the choice of birth distributions, but also to a radio model ill adapted to explain the observed radio population at the highest s. Whereas this would be problematic to study radio beam models, the reasonable representation at ms and the excess of objects with s/s, T, and kpc, where most of the LAT pulsars are found, supports the study of ray models. The necessity of an improved radio model is a result of this paper and its formulation, beyond the purpose of this study, will be the subject of future work. In the histograms shown in Figure 8, the total distributions are dominated by the radio sample since the ray pulsars’ contribution is much smaller.
10.3 The spindown power
Figures 9 and 10 compare the distributions in spindown power and characteristic age for the LAT and the visible simulated pulsars. All models are significantly lacking simulated pulsars with spindown power W and characteristic age kyr. Additionally, all the models overpredict the number of low pulsars and favour a pulsar population older than that observed by the LAT.
The difference in shape of the observed and simulated histograms suggests that the inconsistency is not due to a simple scale mismatch, but to a deficiency in modelling the pulsar evolution: even by scaling the spin down power upward none of the models would be able to describe the observed distribution.
Even though the PC model fails to produce enough visible gammaray pulsars because its narrow beam is under luminous and rarely visible, its evolution with or age is less skewed to old age than the highaltitude SG or the outergap models.
The OG model provides the poorest description of the ray evolution. A strong evolution with age is predicted by the classical formulation of the OG because the gap size is controlled by the amount of Xrays emitted by the stellar surface heated by the backflow of primary charges returning from the gap (selfsustaining OG model). The strong evolution driven by this feedback is apparently not supported by the LAT data. The OPC model gives slightly better results but still fails to predict enough high pulsars. The similarity of the profiles obtained for the SG and OPC models shows that the relative lack of energetic ray pulsars is not related to
the number of visible hemispheres (two pole caustic SG, one pole caustic OPC), or to the evolution of the emission region with age within the open magnetosphere (the emitting layer moves closer to the magnetic axis with increasing age in the OPC case while it remains near the last closed B line, but widens with age in the SG model).
The underprediction of high visible ray pulsars is rather puzzling since they are the intrinsically brightest objects (high particle power and large ) with the widest beams (large open magnetosphere and thin gaps emitting near the closed field lines) sweeping widely across the sky. The problem affects all the models, so its origin does not depend much on the emission pattern or the luminosity trend with . For instance, the luminosity evolution of the OPC model was constructed to agree with the LAT data, yet the deficit of energetic