Probing gamma-ray emissions of -LAT pulsars with a non-stationary outer gap model
We explore a non-stationary outer gap scenario for gamma-ray emission process in pulsar magnetosphere. Electrons/positrons that migrate along the magnetic field line and enter the outer gap from the outer/inner boundaries activate the pair-creation cascade and high-energy emission process. In our model, the rate of the particle injection at the gap boundaries is key physical quantity to control the gap structure and properties of the gamma-ray spectrum. Our model assumes that the injection rate is time variable and the observed gamma-ray spectrum are superposition of the emissions from different gap structures with different injection rates at the gap boundaries. The calculated spectrum superposed by assuming power law distribution of the particle injection rate can reproduce sub-exponential cut-off feature in the gamma-ray spectrum observed by Fermi-LAT. We fit the phase-averaged spectra for 43 young/middle-age pulsars and 14 millisecond pulsars with the model. Our results imply that (1) a larger particle injection at the gap boundaries is more frequent for the pulsar with a larger spin down power and (2) outer gap with an injection rate much smaller than the Goldreich-Julian value produces observed GeV emissions. Fermi-LAT gamma-ray pulsars show that (i) the observed gamma-ray spectrum below cut-off energy tends to be softer for the pulsar with a higher spin down rate and (ii) the second peak is more prominent in higher energy bands. Based on the results of the fitting, we describe possible theoretical interpretations for these observational properties. We also briefly discuss Crab-like millisecond pulsars that show phase-aligned radio and gamma-ray pulses.
keywords:pulsars:general– radiation mechanisms:non-thermal–gamma-rays
The Fermi gamma-ray telescope (hereafter, Fermi) launched in 2008 has facilitated the study of gamma-ray emission process in the pulsar magnetosphere. The Large Area Telescope on-board the Fermi (hereafter Fermi-LAT) has measured the gamma-ray emissions from more than 150 pulsars (Fermi collaboration 2015), and has measured the spectra and the pulse profiles above 1GeV with unprecedented sensitivity. For example, Fermi-LAT found that the gamma-ray flux above the cut-off energy at around GeV decays slower than pure exponential function (Abdo et al. 2010a, 2013). This cut-off behaviour favours the emissions from the outer magnetosphere (e.g. slot gap, outer gap and annular gap) and rules out the classical polar cap scenario, which predicted a super exponential cutoff feature in the GeV spectrum because of the magnetic pair-creation process. Among Fermi-LAT pulsars, 20 pulsars are found to show pulsed emissions in the energy range GeV, including 12 up to GeV (Ackermann et al. 2013) and their spectra clearly indicate sub-exponential cut-off features above the cut-off energy (Ackermann et al. 2013). The pulsed gamma-ray emissions from the Crab pulsar show single power law spectrum above cut-off energy (GeV) and extends to TeV energy bands (Aleksi et al. 2011, 2012, 2014; Aliu et al. 2008, 2011, Abdo et al. 2010b). The GeV/TeV emissions from the Crab pulsar disagree with the spectra of the standard curvature radiation scenario (e.g. Cheng et al. 2000; Takata & Chang 2007; Harding et al. 2008), and will originate from the inverse-Compton scattering process in the outer magnetosphere (Lyutikov et al. 2012; Harding and Kalapotharakos 2015) or pulsar wind region (Aharonian et al. 2012). We (Leung et al. 2014) reported the detection of the pulsed emissions above 50GeV from the Vela pulsar, and showed that the previous models (e.g. Hirotani 2007; Takata et al. 2008) predicted a smaller flux level at 50-100GeV energy bands than the observed flux. A study of sub-exponential spectrum above cut-off energy will discriminate among emission models.
In addition to sub-exponential cut-off behaviour, the Fermi-LAT observations have revealed several interesting relations between the gamma-ray emission properties and the spin down characteristics; (1) the gamma-ray emission efficiency, which is the luminosity divided by the spin down power, decreases with the spin down power, and (2) the spectrum between 100MeV and the cut-off energy at around GeV tends to be softer for a larger spin down pulsars (Abdo et al. 2013), (3) the second peak in the light curve is in general more prominent in higher energy bands (e.g. Crab, Vela and Geminga pulsars, Abdo et al. 2013), and (4) Fermi-LAT millisecond pulsar with a higher spin down power and a larger magnetic field strength at the light cylinder tends to have Crab-like pulse profiles, in which radio/X-ray/gamma-ray pulses are in phase (Ng et al. 2014). Explanations for these observed properties with a model will advance in understanding of the nature of the high-energy emission process in the pulsar magnetosphere.
The cause of the formation of the non-exponential cut-off decay is still in debate. Abdo et al. (2010b) and Vigan and Torres (2015) argued that a sub-exponential cut-off in the observed spectrum could be understood as the superposition of several power law plus exponential cut-off functions with varying the photon index and the cut-off energy, for which the different components are produced at the different region of the pulsar magnetosphere cutting across our line of sight. The contribution of the inverse-Compton scattering process (likewise the Crab pulsar) is one of the proposed models to explain the high-energy tail of the Vela pulsar (e.g. Lyutikov et al. 2012). However, the required soft-photon number density in the magnetosphere to explain the observed GeV flux level will be much larger than one inferred from the optical/UV/IR observations of the Vela pulsar (Takata et al. 2008).
We (Leung et al. 2014) discussed the formation of the spectrum of the Vela pulsar within framework of the outer gap model, and proposed a non-steady model. In this new outer gap model, the electrons and positrons that enter the gap from outer and inner boundaries, respectively, control the gap structure (size, particle distribution and electric field strength etc.) and a smaller rate of the particle injection produces thicker outer gap and harder spectrum. The model suggested that the injection rate much smaller than Goldreich-Julian value produces the observed gamma-ray emissions above 10GeV. We argued that the rate of the particle injection at the gap boundaries could fluctuate with time and the observed gamma-ray spectrum is superposition of the emissions from different stationary gap structures with different injection rates.
In this paper, we will discuss a detail of the three-dimensional calculation method for this new outer gap model, since we did not provide it in our previous observational paper (Leung et al. 2014). In section 4, we present the predicted spectrum and light curve of the Vela pulsar. We will discuss the observed energy dependent light curve. In section 5, we will apply our model to other gamma-ray emitting pulsars, and will discuss how our model interprets the observed relation between the spectral softness below cut-off energy and the spin down power. We will also discuss the Crab-like millisecond pulsar and the limit of our model.
2 Theoretical model
2.1 Pulsar magnetosphere with outer gap accelerator
The global simulations have been developed to investigate structure of the magnetosphere with the high-energy emission region. Earlier particle simulations showed that the magnetosphere with no-pair-creation process settles down into a quiet state with electron cloud above the polar caps, a positively charged equatorial disc and vacuum gaps in the middle latitudes (Krause-Polstorff & Michel 1985; Smith et al. 2001; Wada & Shibata 2007). Recent particle-in-cell simulations have shown the pulsar magnetosphere with the discharged particles created by the pair-creation process. Chen & Beloborodov (2014) discussed that if pair-creation multiplicity is very high at outer magnetosphere around the light cylinder, the outer gap around the light cylinder was quenched and the magnetosphere is similar to the force-free solution with a super Goldrecih-Julian current sheet and the Y-point near the light cylinder, where are main high-energy emission region (Spitkovsky 2006). On the other hand, it is also suggested that if the pair-creation process in the outer magnetosphere is low, the outer gap can survive from the fill of discharge particles and it can be high-energy emission regions (Wada & Shibata 2007; Yuki & Shibata 2012). It is still under debate for the structure of pulsar magnetosphere as well as the high-energy emission region, since the current global simulations are difficult to deal with the realistic pair-creation process by taking into account the position dependent mean free path and soft-photon density.
In this paper, we assume that the pulsar magnetosphere has an outer gap and the high-energy gamma-rays are produced by the curvature radiation process of the discharge pairs inside the outer gap. Our local model precisely calculate the pair-creation rate in the outer magnetosphere. As we will see in section 4.1, the optical depth of the photon-photon pair-creation process around light cylinder is of order of for most of pulsars, and most of the gamma-rays emitted from outer gap can escape from the light cylinder.
2.2 Particle injection at the gap boundaries
To activate the gamma-ray emissions and subsequent pair-creation cascade in the outer gap, the charged particles (electrons and/or positrons) should enter the gap along the magnetic field line from outside the gap; the outer gap will be inactive without the injection of the particles at the gap boundaries. In this paper, we use terminology “injected current” to refer the electric current component carried by the electrons/positrons that enter the outer gap from the gap boundaries. The outer gap thickness in the poloidal plane affects to the magnitude of the accelerating electric field and therefore hardness/luminosity of the curvature emissions; a thinner outer gap produces a smaller accelerating electric field and a softer/fainter gamma-ray emissions. From electrodynamical point of view, we expect that the outer gap has a thickness so that the pair-creation cascade in the gap produces an electric current of order of the Goldreich-Julian value and hence the gap structure will be affected by amount of the particles (i.e. injected current) that enter the gap from the inner and/or outer boundaries. Takata et al. (2006) calculated two-dimensional outer gap structure and investigated the dependency of gamma-ray spectra on injection rates of the particles at the inner and outer boundaries. They demonstrated that a larger injection produces in general a thinner outer gap and a softer gamma-ray spectrum. For the inclination angle less than degree, the positrons and electrons can enter the gap from inner and outer boundaries, respectively.
The physical origin of the injected particles at the gap boundaries are argued as follows. As suggested by Shibata (1991, 1995), the polar cap accelerator, outer gap region, and the pulsar wind region, where the electric current crosses the magnetic field lines, should be connected by the current circulating the magnetosphere. As shown in global simulations (e.g. Yuki & Shibata 2012), we expect that the pair-creation process in the polar cap accelerator will make the current that flows higher latitude around the magnetic pole, while the discharged particles in the outer gap accelerator are main current carriers at lower-latitude region around the last-open field lines. The polar cap accelerator model usually assumes a particle injection from the neutron star surface. For the inclination angle , the electrons from the stellar surface are injected into the polar cap accelerator and initiate the pair-creation cascade process through the magnetic pair-creation and/or photon-photon pair-creation processes (Daugherty & Harding 1996). The discharged pairs form the current that flows higher altitude. Most of particles from the polar cap region will flow out from the magnetosphere and will form the pulsar wind. But it has been suggested that some of negative particles (for degree) cross the magnetic field lines towards equator due to drift (Wada & Shibata 2011l; Yuki & Shibata 2012), where is the radiation drag force, and they eventually return to the star along the magnetic field lines at the lower latitude. It is probable that on the way from the light cylinder to the star, the returning electrons enter to the outer gap along the magnetic field line from the outer boundary. The high-energy emissions by the returning electrons and subsequent pair-creation cascade processes produce the discharged pairs that also contribute to the current flowing lower latitude around the last-open field lines. In section 3.2 we discuss the current conservation along the magnetic field line.
We can argue several possibilities for the physical origin of the positrons that enter the gap from the inner boundary. In the polar cap accelerator, the discharged positrons will return to the polar cap region. If the star continuously absorbs the positrons more than electrons, it would be charged up positively. To keep the charge of the star at constant, the positrons should be re-emitted from stellar surface along the magnetic field lines outside polar cap accelerator. Such positrons could enter the outer gap from the inner boundary and contribute to the electric current flowing along the magnetic field lines that penetrate the outer gap.
Moreover, the gamma-rays produced in the outer gap will create more pairs around the inner boundary (c.f. Figure 2), and residual electric field could separate the charge particles. These discharged pairs could effectively become the injection current at the inner boundary, because the main emission region of the outer gap is beyond the null charge surface. Takata et al. (2010) also argued that the gamma-rays emitted towards the stellar surface by the incoming particles may generate new pairs via the magnetic pair-creation process near the stellar surface, and some new pairs could be returned to the outer gap due to complicated surface magnetic field structure. These returning positrons also could enter the gap from the inner boundary.
2.3 Outer gap with time dependent particle injection
Although the pulsed radio emissions averaged over longer time-scale is stationary, there is a wide range of variability in a shorter times scale in the radio emissions from the pulsar (e.g. Kramer et al. 2002; Lyne et al. 2010; Keane 2013). The micro-second variations seen in single pulse could be produced by spatial fluctuation in the emission region. The pulse-to-pulse variations on the timescale of millisecond to second likely represent timescale of the temporal variation of the structure of the emission region (e.g. time dependent pair-creation process/particle emissions from the stellar surface). The longer timescale (second to year) variations associated with the mode switching and nulling, which sometimes accompany the variations of the spin down rate, could be related with the changes of entire magnetosphere. These observations suggest that the switching between different states of magnetosphere is probably a general feature of the pulsars
In this paper, we assume that the outer gap structure is temporal variable and that the observed gamma-ray emissions are superposition of the different outer gap structures. We argue that the non-stationary behaviour of the outer gap is caused by the time variation of the rate of particle injections at the gap boundaries. We expect that the time-scale of variations is of order of or longer than the crossing timescale of the light cylinder, , where is the light cylinder radius and is the pulsar spin period. For example, as we discussed above, the discharged pairs produced around the inner boundary could effectively become the origin of injected particles at inner boundary. In such a case, the temporal variation of the outer gap will be related to the variation of the pair-creation rate around the null charge surface. Since the pair-creation rate depends on the gamma-ray intensity and surface X-rays intensity, which is affected by the returning particles (c.f. section 3.4), around the light cylinder, the expected time-scale of the variation will be of order of the light-cylinder crossing time-scale . The variation of the electrons returning from the pulsar wind region, which will enter the gap from the outer boundary, will be of order of or longer than the light-cylinder crossing time-scale, since the time-scale shorter than the crossing time-scale may be smoothed out during the travel around global magnetosphere.
We assume that the observed gamma-ray spectrum is a superposition of the emissions from various stationary gap structures with various particle injection rates at the gap boundaries, and the stationary outer gap structure for an injection rate forms with the crossing time-scale . For a fixed particle injection rate, our stationary solution will be stable for a small perturbation. For example, if the accelerating electric field increase from the stationary solution, the curvature photon energy and hence pair-creation rate increase from the stationary solution. The increase of the number of pairs try to screen the perturbed electric field.
3 Basic Equations
In this section, we describe our basic equations for solving the gap structure with a fixed injection rate at the boundary. By using vacuum rotating dipole magnetic field, we solve the Poisson equation to obtain the accelerating electric field (section 3.1), for which the charge density in the gap is obtained by solving the continuity equations for the electrons and positrons (section 3.2) with the curvature radiation process and pair-creation process (section 3.3). In section 3.2, we will discuss the conservation of the electric current along the magnetic field line.
3.1 The accelerating electric field
We investigate the outer gap structure under the steady condition that with being spin angular frequency. The electric field along the magnetic field line arises in the charge depletion region from so called Goldreich-Julian charge density, and it accelerates the positrons and electrons to an ultra-relativistic speed. The Poisson equation for the accelerating electric field is written as
where is the space charge density and is the Laplacian. In addition, is the Goldreich-Julian charge density, where is the component of the magnetic field projected to the rotation axis. The accelerating electric field along the magnetic field line is computed from , where is the distance along the magnetic field line.
To solve the Poisson equation (1), we adopt coordinate system based on the distance along the field line, , from the star () and the magnetic coordinates, and , which are angles measured from and around the magnetic axis, respectively (Hirotani 2006). We define at the north magnetic pole and (magnetic meridian) at the plane that includes the rotation axis and north magnetic pole for inclined rotator. The coordinates relate with the canonical spherical coordinates , for which axis coincides with the rotation axis, as
where is the stellar radius, is the local magnetic field strength and are components of the magnetic field, respectively. We define that corresponds to the rotation axis and is the magnetic meridian. We can relate between () and () as and with being the inclination angle.
3.1.1 Boundary conditions
For 3-D outer gap, there are six boundaries, that is, inner (stellar side), outer (light cylinder side), lower, upper, leading side and trailing side boundaries. For inclined rotator, the charge deficit region at the azimuthal angle is in general less active, because the null charge surface is located close to the light cylinder, and because the electric field is too small to boost the charge particles up to ultra-relativistic speed that can produce the high-energy gamma-rays. In this paper, therefore, we put the numerical boundaries on the magnetic field lines labelled by for the leading side (positive sign) and the trailing side (negative sign) of the gap, and impose the mathematical boundary conditions that .
For fixed azimuthal angle , the lower and upper gap boundaries lay on the magnetic field lines. We fix the lower boundary at the last-open field lines and impose on it. We also impose on the upper boundary and solve it’s position, for which the gap can create an assumed electric current density (c.f. section 3.5). In the calculation, we set the outer boundary near the light cylinder and impose on the boundary. We initially apply the numerical boundary at and solve the gap dynamics. If the electric field changes its sign around the given outer boundary, then we set new outer boundary at the location where the solved electric field changes its sign, because we anticipate that the outer gap should be unstable if the field-aligned electric field changes its sign inside the gap.
Finally, let us consider the inner boundary (stellar side). Because we assume that there is no potential drop between the stellar surface and the inner boundary, we impose the conditions and . Since arbitrary given boundary does not satisfy both the Dirichlet-type and the Neumann-type conditions, we seek for the appropriate boundary by moving the boundary step by step. With two-dimensional analysis, Takata et al. (2004) discussed that the inner boundary of the outer gap starts from the position where the charge density of the current carriers is equal to the Goldreich-Julian charge density. For example, the outer gap starts from the null-charge density of the Goldreich-Julian charge density, if the gap is vacuum. On the other hand, the inner boundary will locate on the stellar surface, if the electric current created inside the gap is in units of the Goldreich-Julian value, .
3.2 Continuity equations
In this paper, we assume that the inclination angle of the magnetic pole is less than 90 degree. In such a case, the positive electric field along the magnetic field line accelerates the positrons towards light cylinder and electrons towards the stellar surface, respectively. In the outer gap, we can anticipate that new born pairs in the gap are immediately charge separated and are boosted to ultra-relativistic speed by the electric field along the magnetic field line (Hirotani & Shibata 1999); that is, we can assume that all positrons and electrons in the outer gap move towards the light cylinder and towards the stellar surface, respectively, with the speed of light. Under these conditions, the continuities of the number density of the positrons (plus sign) and of the electrons (minus sign) may be written as
where is the source term due to photon-photon pair-creation process. The electric current density per magnetic flux tube in the gap is given by . Fixing (), the continuity equation (5) satisfies the current conservation along the field line, that is,
where we normalized the current density by the local Goldreich-Julian value. We define the normalized current densities carried by the positrons and electrons as
and . With the equation (5), the solutions for can be written as
respectively, where , and represent the positions of the inner boundary and outer boundary, respectively, and the injection current (or ) represents number of positron (or electron) that enters the gap from the inner (or outer) boundary per unit time per unit area and per magnetic flux tube. The origin of the particles injected into the gap were discussed in section 2.2. In terms of (), the conservation of the electric current along the magnetic field line becomes
which represents the current component carried by the created pairs in the outer gap (hereafter we use terminology “gap current” to refer ). Equation (7) tells us that the electric current along the magnetic field line is sum of the injection currents at gap boundaries plus gap current. We note that as long as the current flows along the magnetic field line that penetrates the outer gap, the magnitude of current density per magnetic flux tube is equal to both outside and inside the gap. Hence there is no current discontinuity along the magnetic field line. To close the current circuit, the trans-field current flow should appear in somewhere beyond the light cylinder (Shibata 1991, 1995). In this paper, since the structure of the magnetosphere outside the light cylinder is beyond out of scope, we just assume that the cross-field region is far from the outer boundary, and that the injected electrons cross the outer boundary along the magnetic field line.
Actual values for the total current , injected currents and should be solved with the complicated physics (e.g. energy-angular loss relation among the polar cap accelerator, outer gap and pulsar wind region, Shibata 1991) of the global pulsar magnetosphere. For example, the injection current might be solved together with the outer gap activity and positron re-emission from the neutron star surface, which is related to the charge redistribution over the polar cap region. As we discussed in section 2.2, the injection current at the outer boundary will be related to the physics of the formation of the pulsar wind. The total current running through the outer gap should be solved with global pulsar magnetosphere including the polar cap, outer gap and pulsar wind region. Because of the large theoretical uncertainties of the global structure, however, our local model treats () or as a set of the free parameters. In section 3.5, we describe how our model assumes the values of .
3.3 Curvature radiation and pair-creation processes
To calculate the source term in equation (5), we compute the pair-creation process between the gamma-rays emitted by the curvature radiation and thermal radiation from the stellar surface. We calculate the Lorentz factor of the accelerating electric field by assuming force balance between the acceleration force and the back reaction force of the curvature radiation process as
where is the curvature radius of the magnetic field line. The number of curvature photons emitted per unit time from the particle with a Lorentz factor is
The spectrum of the curvature radiation from the particle is described by
where is the modified Bessel function of the order of .
The emitted curvature photons may convert into new electron and positron pairs through the pair-creation process. The mean free path of the pair-creation is
with being the X-ray number density between energies and , the collision angle between an X-ray photon and a gamma-ray photon, and the threshold X-ray energy for the pair creation. In addition, the pair creation cross-section is given by
and is the Thomson cross-section. In this paper, we consider the thermal X-ray photons from the stellar surface. At the radial distance from the centre of the star, the photon number density between energy and is given by
where is the effective radius and refers to the surface temperature.
3.4 X-ray emissions from NS surface
In this paper, we consider two types of the surface emissions, namely, the neutron star cooling emissions and heated polar cap emissions. Both surface emission processes could contribute to the observed thermal X-ray emissions from the young/ middle-age pulsars (e.g. Caraveo et al. 2004 for the Geminga pulsar), while only heated polar cap emission should be important for the millisecond pulsars (Zavlin 2007; Takata et al. 2012)
For the neutron star cooling model, the temperature as a function of the age really depends on the neutron star model (Yakovlev & Pethick 2004), as shown in Figure 1. Although the observations of surface temperature could exclude the high mass neutron star with -condensate core model, it is still under debate for the neutron star model. In this paper, therefore, we use the cooling curve predicted by the standard model (thick line in Figure 1), which will provide typical surface temperature for fixed age of the neutron star. We assume the spin down age , where is the time derivative of the spin period, for the true age of the young/middle-age pulsars
For the heated polar cap emissions, we apply the model by Takata et al. (2012), in which the X-ray emissions from the heated polar cap region are composed of two components, namely, a core component and a rim component. The core component shows a higher temperature but a smaller effective radius (cm), and the rim component has a lower temperature but larger effective radius (cm). In their model, the bombardment of the returning pairs on the polar cap region causes the core component, while the irradiation of MeV gamma-rays that are emitted near the stellar surface by the returning particles heats up the surface and produce the rim components. We expect that with a smaller effective area (cm), the core component does not illuminate the outer gap, while the rime component with the effective radius cm has a greater likelihood of illuminating the outer gap. This model predicts the temperature of the rim component as (c.f. equation 21 in Takata et al. 2012)
In this paper, we assume cm to match with typical observed temperature and effective radius of the millisecond pulsars.
3.5 Model parameters
In the present model, we treat the injection currents and the gap current as the model parameters and they are relate to the total current as (c.f. section 3.2). In addition, the inclination angle is model parameter and we will fix the inclination angle in this paper. Since we focus on the observed phase-averaged spectrum, we do not introduce the Earth viewing angle. By integrating the emissions from whole outer gap region, we will compare the model spectrum with the observed phase-averaged spectrum.
For the injection currents at the gap boundaries, we assume constant over the boundaries and we impose also for reducing the model parameters, that is, we assume same particle injection rates at the inner and outer boundaries. Choice of the equal injection rates at the gap boundaries is arbitrary, and it is not necessary for the real case. For the photon-photon pair-creation process with the X-rays from the neutron star surface, the gap structure is more sensitively to the choice of the injection current at the outer boundary. This is because the pair-creation process between the gamma-rays emitted by the inward migrating electrons and the surface X-rays are head-on collision process, while the pair-creation process of the outward propagating gamma-rays from the positrons is tail-on collision process. Hence, most of the pairs are created by the inward propagating gamma-rays. We expect that if we assume no injection current at the outer boundary (i,e, ), the outer gap size will become thicker than the solutions discussed in this paper. We will study the gamma-ray emissions from the outer gap with in the subsequent papers. Here we define total injection current , as
which is time variable quantity in our model. In the model, we will apply (see section 3.6).
The gap current, is limited as follows. Figure 2 represents an example of the calculated gap structure in the plane defined by ; the left panel shows the photon-photon pair-creation rate in the gap and the right-hand panel shows the trans-field distribution of the total current density (). We can see in the figure that the calculated gap current () increases as increase of the height from the lower boundary (last-open field line), and that there is a maximum value of the gap current (). This is because the gamma-rays propagate in the convex side of the magnetic field lines. Around the upper boundary, the gap current decreases because the electric field decreases there and because the gamma-rays emitted at lower region do not illuminate upper region of the gap. In the Figure 2, the gap current therefore becomes maximum at around 70-80% of the gap thickness, and the maximum gap current is .
In the paper, we teat as the model parameter. We can find in Figure 2 that the position of the inner boundary on the magnetic field line that has a large gap current shifts towards stellar surface from the GJ null charge surface. As discussed in Takata et al. (2006), the inner boundary of the outer gap will locate near the stellar surface, if the created current is of order of . Within the framework of the calculation method, however, it is difficult to obtain such a stable solution, in which the field aligned electric field does not change its direction, if the gap current approaches to . In the calculation, hence, we assume the maximum gap current with a value slightly smaller than . In the model calculation, we assume the inclination angle , and we solve the location of the upper boundary so as to create the gap current of , which does not depend on the azimuthal angle (but see section 3.6 for large ). As long as , the exact value of does not affect much on the calculated gamma-ray spectra.
3.6 Calculation Process
With the specified injection current, , and fixed maximum gap current , we self-consistently solve the outer gap structure, as follow. We start the calculation by solving the Poisson equation (1) for a vacuum outer gap with a very thin thickness. Using the calculated electric field along the magnetic field line, we calculate the terminal Lorentz factor (8) at the each calculation grid. Given the injection current, , we solve the continuity equation (5) with the curvature radiation and pair-creation processes, and then we obtain new distribution of the charge density inside the gap. With the new charge distribution, we solve the Poisson equation to update the electric field, which subsequently modifies the charge density distribution. We iterate this procedure until all physical quantities converge.
If the gap thickness is too thin, the magnitude of the electric field is not enough high to boost the electrons/positrons up to ultra-relativistic speed, and the pair-creation cascade does not initiate in the gap. As a next step, therefore, we increase the thickness of the gap. For a fixed magnetic azimuth, the gap current has a distribution in the direction of the latitude , as Figure 2 indicates. If the maximum current density at fixed magnetic azimuth is smaller than (for ), we slightly increase the gap height. The increase of the gap height produces the increase of the accelerating electric field and results in the increase of the gap current. We solve the outer gap dynamics with new upper boundary and obtain new distribution of the gap current. Updating the gap upper boundary step by step, we finally obtain the desired gap structure for the fixed injection current .
In the model, the trans-field thickness of the outer gap is a function of the magnetic azimuth , and the gap thickness is the minimum at around the magnetic meridian . This is because the gap thickness relates to the radial distance to the null charge point on the last-open field line from the stellar surface. At around , the radial distance to the null charge point becomes minimum (c.f. Figure 3), and hence the number density of the surface X-rays around the inner boundary of the gap becomes maximum. Because the pair-creation rate increases as increasing of the number density of the surface X-rays, the pair-creation rate inside the outer gap becomes maximum around the magnetic meridian. As a result, the gap thickness becomes minimum around the magnetic meridian.
If the pair-creation rate is very low, the outer gap can become very thick. For example, on the magnetic field lines labeled by the azimuthal angle , since the null charge point at the last-open field lines are located near the light cylinder (c.f. Figure 3), the pair-creation rate is very low and therefore the outer gap can become very thick. In the calculation, we set possible maximum thickness at , namely 80% of the open field line region for the fixed . For such an azimuthal angle, the maximum gap current is smaller than . We note that there is critical magnetic field line for fixed , above which the null charge points on the magnetic field lines locate outside the light cylinder, and therefore a part of the outer gap in the calculation would locate outside right cylinder. In the present calculation, we ignore the radiation process and the pair-creation process outside the light cylinder, (1) because the emissivity of the curvature radiation and the pair-creation rate will be very low compared to those inside the light cylinder and (2) because the special relativistic effect (e.g. retarded electric potential) and magnetic field bending due to the magnetospheric electric current should be taken into account to obtain the correct gap structure outside the light cylinder. For the very thick outer gap, the pair-creation rate around the upper boundary is negligibly low and most emissions are produced in lower gap region.
Our model assumes that the observed gamma-ray spectrum is a superposition of the emissions from various gap structures with the various particle injection rates at the gap boundaries. Since our local model cannot determine the distribution of the injection rate, which will relate to the physics in the source region (e.g. polar cap), we assume a power-law distribution with
where is the normalization factor and it is calculated from . We fix the minimum injection rate at , because the solved outer gap for the most pulsars has the maximum thickness, , for . We set the maximum injection rate at so that the injection rate is smaller than created current in the gap . With the function form of equation (17), a larger (or smaller) injection current dominates in the distribution for the power-law index (or ). The superposed spectrum becomes
where is the gamma-ray spectrum for a fixed injected rate.
4 Application to the Vela pulsar
In this section, we will apply the model to the Vela pulsar (PSR J0835-4510) and will discuss the general properties of the outer gap structure and the gamma-ray spectrum predicted by the model.
4.1 Pair-creation in the gap
Figures 2 and 3 show the created gap current (in units of the Goldreich-Julian current density) per unit length, , in the gap and (magnetic meridian). In Figure 2, the bottom () represents the lower boundary and top () is the upper boundary of the gap. The results are for the injection rate , that is, in the present assumption.
We can see in the figures that the photon-photon pair-creation process creates a more gap current at around the inner boundary. This is because (1) the inward propagating gamma-rays mainly produce the pairs and (2) the mean free path is shorter at the inner magnetosphere. Our result confirms the results of our previous calculations (Cheng et al. 2000; Takata et al. 2006) and recent 3-D calculation (Hirotani 2015).
The field aligned electric field separates the electrons and positrons created inside the gap, which migrate inward and outward directions, respectively, for the inclination angle . Since most of pairs are created near the inner boundary, the positrons will feel most of full electric potential drop before escaping from the gap outer boundary, while the electrons will feel smaller potential drop between the inner boundary and the pair-creation position. As a result, the radiation power from the positrons is about factor of ten larger than that from the electrons. Figure 5 shows the spectra of the gamma-rays emitted by outward (solid line) and inward (dashed line) migrating particles. This result is also consistent with the previous studies (Takata et al. 2006; Hirotani 2015).
Figure 4 shows the distribution of the electric current carried by the positrons (solid line) and electrons (dashed line) along a magnetic field line for stationary outer gap. The electric field in the gap discharges the electrons and positrons and increases the electric current. In the figure, the current is constant below low inner boundary, which is located at , since we assume there is no electric field along the magnetic field line between the stellar surface and the inner boundary of the outer gap. In the outer magnetosphere around the light cylinder, the optical depth of the photon-photon pair-creation process is so low that the current density is almost constant along the magnetic field line.
One can estimate the pair-creation mean-free path and multiplicity around the light cylinder. The mean free path of the photon-photon collision for a gamma-ray may be estimated from
where is the collision angle, with being the temperature of the neutron star surface, stellar radius, Stephan-Boltzmann constant and Boltzmann constant. In addition, we assume the cross-section as with being Thomson cross-section. Optical depth is
which is much smaller than unity for the middle age pulsars and the millisecond pulsars.
A charge particle emit the gamma-rays with a rate of
A pair multiplicity by a charge particle accelerated inside the gap may be estimated as
4.2 Gap structure
Different injection rates produce different outer gap structures and the gamma-ray spectra. Figures 6 and 7 show the distribution of the electric current and the gap thickness, respectively, as a function of the magnetic azimuth. In addition, Figure 8 shows the gamma-ray spectrum for fixed injection rate. The solid line, dashed line and dotted line in the figures show the results for the different injection rate , and , respectively. We can see in Figure 6 how the azimuthal width of the “active” gap region (large current region) depends on the injection rate; the active region is wider for a larger injection rate. Figure 7 and Figure 8 show that as the injection rate increases, the averaged gap thickness becomes thinner and therefore the gamma-ray spectrum becomes softer.
We can find in Figure 8, the different injection rates produce similar amount of the gamma-ray luminosity. The pulsar electrodynamics implies that the gamma-ray luminosity is of order of , where is the total current flowing into the gap. As the injection rate increases, the total current becomes larger, while the potential drop, which depends on the gap thickness as , becomes smaller. Since these two effects compensate each other, the gamma-ray luminosity is insensitive to the injection rate.
Figure 8 also shows that gamma-ray spectrum for a fixed injection rate does not fit the observed spectrum in 0.1-100GeV of the Vela pulsar. With a small injection rate , the calculated spectrum explains the observed flux level above 10GeV, but the predicted spectral slope below 10GeV is steeper than the observed one. For the large injection , on the contrary, the predicted flux above the cut-off energy decays faster than the observed one, and it is difficult to reconcile with the observed flux above 10GeV. With the present framework of the 3-D calculation, we would expect that the superposition of the emissions from the different outer gap regions is not the main reason for the sub-exponential cut-off behaviour of the Vela pulsar.
4.3 Gamma-ray spectrum
We assume that the observed gamma-ray spectrum is a superposition of the emissions from various stationary gap structures with various injection currents () at the boundaries, for which we assume power-law distribution of the injection current (17), . We integrated the emissions of entire outer gap regions and used minimum -square method to find the best fitting index, , and normalization for the observed phase averaged spectrum.
Figure 9 compares the best fitting model spectrum with the phase-averaged spectrum for the Vela pulsar; the solid line shows the calculated spectrum with using the best-fitting index , implying that a larger injection rate dominates in the distribution. The dashed, dotted and dashed-dotted curves in Figure 9 show the contributions for the injection , and , respectively. As we can see in the figure, our model suggests that the emissions from the outer gap with smaller injection rates , produces observed spectrum above 10GeV, although it’s integrated flux in 0.1-100GeV energy bands is much smaller than the observed one. We also see that the emissions from the outer gap with a larger injection rate mainly contributes to the observed integrated flux, but it’s spectrum (e.g. dashed line) above 10GeV decays faster than the observed spectrum . The cut-off feature of the model spectrum (solid line) is in good agreement with the sub-exponential decay of the observations.
Around 100MeV, the model spectrum is steeper than the observed spectrum. This may imply that the distribution of the particle injection rates deviates from the simple power-law function, which has been assumed in the present calculation. Moreover, we have ignored the contributions of inward emissions, because the luminosity of the outward emissions is about one order of magnitude larger than that of the inward emissions, as Figure 5 shows. Around 100MeV, however, the flux level of the inward emissions is only several factor lower than that of the outward emissions and could contribute to the observed emissions. The inward emissions probably contribute to the non-thermal X-ray emissions from the Vela pulsar, which shows the multiple (three or four) peaks in the X-ray light curve (Harding et al. 2002; Takata et al. 2008).
Wang et al. (2010) proposed a two layer outer gap model, which divides the outer gap into two regions, namely, main acceleration region and thin screening region around the upper boundary. In the main acceleration region, the electric current density is smaller than Goldreich-Julian value and the curvature emissions produce GeV gamma-rays. In the screening region, a super Goldreich-Julian current screens the field aligned electric field, and the curvature radiation produces 100-500MeV gamma-rays. Within the framework of the present calculation method, it is difficult to reproduce the stationary outer gap, in which the field aligned electric field is positive-definite, with a super Goldreich-Julian current density. A more detail investigation will be necessary to explain the observed emissions around 100MeV of the Vela pulsar.
4.4 Light curves
The top and bottom panels in Figure 10 present the calculated light curves for 1GeV and GeV, respectively, of the Vela pulsar. The inclination angle and the viewing angle are and , respectively. In the figure, the rotation phase 0 and 0.5 correspond to the times when the south magnetic pole () and north magnetic pole (), respectively, point towards the observer.
We find in Figure 10 that the second peak is more prominent in higher energy bands, which is consistent with the observations (Leung et al. 2014). In present model, the outer gap emissions with a smaller injection rate explain the observed emissions above 10GeV of the Vela pulsar, as Figure 9 shows. For a smaller injection, the pair-creation process in the outer gap produces the pairs only on the magnetic fields around , as Figure 6 shows. As a result, the gamma-rays from the outer gap with a smaller injection rate are observed at around orbital phase in the light curve, where is the position of the second peak.
We would like to note that this tendency of energy dependent light curve predicted by the present model is general behaviour for all gamma-ray pulsars, since higher energy photons originate from the outer gap with a smaller injection current. Our model could provide a reason why the second peaks of the Fermi-LAT gamma-ray pulsars (Abdo et al. 2012) are more prominent in higher energy bands.
The model fitting for various Fermi-LAT pulsars may enable us to discuss how the power-law index relates to the spin down parameters. Since it is time consuming task to investigate for all Fermi-LAT pulsars (), we applied the model to 43 young/middle-age pulsars (YPSRs) and 14 millisecond pulsars (MPSRs). We chose those pulsars from the first Fermi-LAT 10GeV source catalog (Ackermann et al. 2013) and from the first Fermi-LAT pulsar catalog (Abdo et al. 2010b); we exclude the Crab pulsar, since its emission process is more complicated. We also apply the model to original millisecond pulsar J1939+2138 and black widow J1959+2048 (Guillemot et al. 2012) though they are not listed in two catalogs. Tables 1-3 summarize the spin down parameters for pulsars fitted in this study, and Figures 11-13 compare the best fitting model spectra with the Fermi-LAT spectra, for which we re-analyzed about 6 years (from 2008 August to 2014 August) Fermi data for the pulsars listed in first Fermi-LAT 10GeV source catalog, while we referred the published Fermi spectra (Guillemot et al. 2012; Abdo et al. 2013) for other pulsars. To obtain the best fitting model spectra with minimum -square method, we use the data points at the center of the errors (namely, the values at the filled boxes in each panel of Figures 11-13). We ignored the data at the lowest energy bin for fitting process since Fermi data at MeV may contain a larger uncertainty. The last columns in Tables 1-3 summarize the best fitting power-law index for the distribution of the injection rate
5.1 Injection current and spin down parameters
We investigate how the fitting power-law index relates to the spin down parameters. In Figure 14, we plot the fitting power-law index as a function of the spin down parameters, namely, rotation period (top-left panel), surface dipole magnetic field (top-right), spin down age (bottom-left) and spin down power (bottom-right). In each panel, and are factors of the linear correlation for the young/middle-age pulsars (open circles) and millisecond pulsar (filled boxes), respectively. We find no correlation between the fitting power-law index and the surface magnetic field (top right). For the rotation period, the correlation is strong for the millisecond pulsars but it is relatively weak for the young pulsars.
With the current samples, the correlation between the fitting power-law index and spin down power (bottom right in Figure 14) is relatively stronger for both young pulsars and the millisecond pulsars. There is a tendency that the fitting index increases with increasing of the spin down parameter, implying a larger current injection is more frequent for the pulsar with a larger spin down power. This tendency of the fitting index actually reflects the fact that the observed spectra below cut-off energy tends to be softer (i.e. larger photon index) for the pulsar with a larger spin down power (Abdo et al. 2013). In other words, our model provides an explanation why the observed spectrum below cut-off energy is softer for higher spin down pulsars.
The reason why typical amount of the injection rate increases with increasing of the spin down power would relate to the increasing of the available potential drop with the spin down power. The spin down power relates to the available potential drop on the polar cap region , namely,
We speculate that the injection currents originate from the pair-creation process in the acceleration region outside outer gap. The pair-creation rate should depend on the available potential drop , since the potential drop in the acceleration region is proportional to the available potential drop. For a larger available potential, a larger accelerating electric field will arise in the acceleration region, and hence more pairs that eventually migrate towards the outer gap will be created inside and outside acceleration region. Hence, we expect that the pulsar with a larger spin down power tends to produce a larger particle injection into the outer gap.
5.2 Class II millisecond pulsars
Pulsars with gamma-ray peak lagging, aligned with, and preceding the radio peak are divided into classes I, II, and III, respectively (Venter et al 2012). For young/middle-age pulsars, only Crab and Crab-twin (PSR J0540-6919) in LMC (Fermi-LAT collaboration, 2015) show the class II radio/gamma-ray pulse profiles. For Fermi-LAT millisecond pulsars, the sources with higher spin down power and stronger magnetic field at the light cylinder in general belong to class II pulsars (Ng et al. 2014). Non-thermal X-ray pulse profiles of the class II pulsars show similar peak structure and generally align with the gamma-ray and radio peaks. In Table 3, the symbol “*” indicates the class II millisecond pulsars.
Observed pulsed radio wave from the class II pulsars probably originates from the plasma process relating to the outer gap accelerator. We speculate that radio emission could be generated above outer gap region when copious amount of the pairs are created in outer magnetosphere, since this could give a shorter time-scale of plasma instability (Ng et al. 2014). As we can see in Table 3, our fitting suggests that a larger amount of the particles () are injected into the outer gap of the class II millisecond pulsars. This tendency would be preferable for our speculation for the radio emission process of the class II millisecond pulsars, since a larger injection rate can produce a more pairs outside the outer gap. We also note that the class II millisecond pulsars accompany the giant pulses (Romani & Johnston 2001; Knight et al. 2006), whose phase positions are in phase or close to the pulse peaks of normal pulsed radio emissions. Hence the mechanism of the giant pulses could relate to the large injection of the particles into the outer gap, which results in creation of large amount of the pairs. Hence our model expects a correlation between the giant radio pulses and X-ray/gamma-ray emission properties.
5.3 Middle age pulsars; J0633+1746 and J1836+5925
The middle-age pulsars J0633+1746 (known as Geminga) and J1836+5925 show that the cut-off behaviours above 10GeV decay slower than pure exponential cut-off, as Figure 11 shows. Within the present framework of the calculation it is difficult to produce 10GeV emissions from the outer gap for those two middle-age pulsars, and therefore there is a large discrepancy between the calculated and observed spectra. The typical potential drop and the accelerating electric field inside the gap are of order of and , respectively, where is the ratio of the gap thickness to the light cylinder radius at the light cylinder and it is of order of unity for middle-age pulsars. The typical Lorentz factor of the particles inside outer gap for middle-age pulsars is , which yields the typical energy of the curvature photons,
where we assumed . The spin down parameters of J0633+1746 and J1836+5925 provides GeV and GeV, respectively. This curvature photon energy can explain the position of the spectral cut-off around 2GeV, but it is difficult to reproduce the observed emissions above 10GeV for these middle-age pulsars.
Vigan and Torres (2015) fit the observed spectrum of the J0633+1746 by parameterizing the magnitude of the accelerating electric field in the outer magnetosphere. They argued that the observed flux peak position around 2GeV requires an accelerating electric field of order of V/m, which corresponds to a potential drop of order of V, namely of order of the available potential drop of the Geminga pulsar. Their model phase averaged spectrum also decays faster than the observation above 10GeV.
Takata & Chang (2009) argued if the last-open field lines could be different from the conventional one that is tangent to the light cylinder. Since the magnetic field must be modified by the rotational and the plasma effects in the vicinity of the light cylinder, the size of the polar cap could be larger than that of the pure dipole magnetic field (Romani 1996; Contopoulos et al. 1999; Gruzinov 2005). Since the available potential drop is proportional to the square of the polar cap radius, the model flux above 10GeV could increase as increasing of the polar cap size.
5.4 Effects of viewing geometry; PSR J0659+1414
We find in Figure 12 that spectral peak energy (GeV) of PSR J0659+1414 is significantly smaller than the model peak at 1GeV, which corresponds to the minimum curvature photon energy for the pair-creation process, . For PSR J0659+1414, it is likely that the Earth viewing angle cuts through the edge of the gamma-ray beam. The outer gap model predicts that more gamma-ray power is released in the direction of measured from the rotation axis, and hence the Fermi-LAT has preferentially detected pulsars with a larger inclination angle and larger viewing angles (Watters & Romani 2011; Takata et al. 2011). The observation bias would explain the double peak structure in the light curves for most of the Fermi-LAT pulsars. For PSR J0659+1414, the gamma-ray light curve shows single peak and furthermore the gamma-ray luminosity divided by the spin down power is (Abdo et al. 2012), which is one or two order of magnitudes smaller than those of the pulsars with similar spin down power. Hence, we expect that Earth viewing angle of PSR J0659+1414 greatly deviates from .
It could be also possible that we observe the inward emissions for PSR J0659+1414. There are very soft gamma-ray pulsars, which are dim in the Fermi-LAT bands but bright sources in hard/soft gamma-ray bands (e.g. PSR B1509-54, Wang et al. 2014; Kuiper & Hermsen 2015). We (Wang et al. 2014) suggested that the GeV-quiet soft gamma-ray pulsars are peculiar cases of the viewing geometry with the Earth viewing angle of , for which the outward emissions from the outer gap is out of line of sight. In the model, the inward emissions from the outer gap produce the observed spectrum of the GeV-quiet soft gamma-ray pulsars. With single pulse light curve and low efficiency of the gamma-ray emissions, PSR J0659+1414 could be another candidate for which we observe the inward emissions.
In summary, most of Fermi-LAT pulsars show that the spectral behaviour above cut-off energy decays slower than pure-exponential function. We discussed this sub-exponential cut-off feature with non-stationary outer gap accelerator. For the outer gap accelerator, the electrons and/or positrons that enter the gap from the inner and/or outer boundary control the gap structure and characteristic of the gamma-ray spectrum. We found that if rate of the particle injection at the gap boundaries fluctuate with time, the gamma-ray spectrum from the outer gap forms a sub-exponential cut-off feature. This model predicts that the emissions above 10GeV originate from a thicker outer gap with a small injection current, which also provides a theoretical explanation why the second peak is more prominent in higher energy bands. The observed gamma-ray spectrum below cut-off energy tends to be softer for the pulsar with a higher spin down rate. This observed tendency is explained if a larger rate of the particle injection is more frequent for a higher spin down pulsar. The class II millisecond pulsars are very unique gamma-ray emitting pulsars. Observed pulsed emission in radio/X-ray/gamma-ray bands and giant radio pulses probably originate from single site or neighbouring regions in outer magnetosphere. A large injection into the outer gap and subsequent pair-creation cascade of the class II millisecond pulsar will create copious pairs outside outer gap, which would enable to develop a plasma process for the radio emission. We expect that future studies for the evolution of the gamma-ray emission properties with the spin down power and the correlation of the radio/X-ray/gamma-ray emissions of the class II millisecond pulsars advance in understanding for nature of the multi-wavelength emission processes in the pulsar magnetosphere.
We express our appreciation to an anonymous referee for useful comments and suggestions. We thank C.-Y. Ng, A. H. Kong, C. Y. Hui, P. H. T. Tam, M. Ruderman, and S. Shibata for the useful discussions. J.T. and K.S.C. are supported by a GRF grant of Hong Kong Government under HKU17300814P. J.T. is supported by the NSFC grants of China under 11573010.
- abdo (2010) Abdo et al., 2010a, ApJ, 708, 1254
- abdo (2010) Abdo et al., 2010b, ApJ, 713, 154
- abdo (2010) Abdo et al., 2010c, ApJS, 187, 460
- abdo (2013) Abdo, A.A. et al., 2012, ApJS, 208,17
- acker (2013) Ackermann, M., et al., 2013, ApJS, 209, 34
- aharonian (2012) Aharonian, F.A., Bogovalov, S.V., Khangulyan, D., 2012, Nature, 482, 507
- aleksic (2011) Aleksi et al. 2011, ApJ, 742, 43
- aleksic (2012) Aleksi et al. 2012, A&A, 541, 13
- aleksic (2015) Aleksi et al. 2014, A&AL, 565, 12
- aliu (2008) Aliu, E., et al., 2008, Sci, 322, 1221
- aliu (2011) Aliu, E., et al., 2011, Sci, 334, 69
- cavero (2004) Caraveo, P.A., De Luca, A., Mereghetti, S., Pellizzoni, A., Bignami, G.F., 2004, Sci, 305, 376
- chen (2014) Chen, A.Y., Beloborodov, A.M., 2014, ApJL, 795, 22
- cheng (2000) Cheng, K.S., Ruderman, M., Zhang, L., 2000, ApJ, 537, 964
- cont (1999) Contopoulos, I., Kazanas, D., Fendt, C., 1999, ApJ, 511, 351
- dau (1996) Daugherty, J.K., Harding, A.K., 1996, ApJ, 458, 278
- fermi (2015) Fermi-LAT collaborations, 2015, Sci, 13, November, 801-805
- crun (2005) Gruzinov, A., 2005, PhRv Letter, 94, 1101
- gullemot (2012) Guillemot, L. et al. 2012, ApJ, 744, 33
- harding (2002) Harding, A.K., Strickman, M.S., Gwinn, C., Dodson, R., Moffet, D., McCullock P., 2002, ApJ, 576, 376
- harding (2008) Harding, A.K., Stern, J.V., Dyks, J., Frackowiak, M., 2008, ApJ, 680, 1378
- harding (2015) Harding, A.K., Kalapotharakos, C, 2015, ApH, 811, 63
- hirotani (1999) Hirotani, K., Shibata, S., 1999, MNRAS, 308, 54
- hirotani (2006) Hirotani, K., 2006, ApJ, 652. 1475
- hirotani (2007) Hirotani, K., 2007, ApJ, 662, 1173
- hirotani (2015) Hirotani, K., 2015, ApJ Letter, 798, 40
- Keane (2013) Keane, E. F., 2013, Proceedings of the International Astronomical Union, 291, 295
- Knight (2006) Knight, H.S., Bailes, M., Manchester, R.N., Ord, S.M., 2006, ApJ, 653, 580
- kramer (2002) Kramer, M., Johnston, S., Van Straten, W. 2002, MNRAS, 334, 523
- krause (1985) Krause-Polstorff, J., Michel, F.C., 1985, A&A 144, 72
- kuip (2015) Kuiper, L., Hermsen, W., 2015, MNRAS, 449. 3827
- leung (2014) Leung, Gene C. K., Takata, J., Ng, C.W., Kong, A.K.H., Tam, P.H.T., Hui, C.Y., Cheng, K.S., 2014, ApJ Letter, 797, 13
- lyne (2010) Lyne, A., Hobbs, G., Kramer, M., Stairs, I., Stappers, B., 2010, Sci, 329, 408
- lyu (2012) Lyutikov, M., Otte, N., McCann, A., MNRAS Letter, 422, 311
- ng (2012) Ng, C.-Y., Takata, J., Leung, G.C.K.m Cheng, K.S., Philippopoulos, P., 2014, ApJ, 787, 167
- ozel (2013) zel, F, 2013, RPPh, 76, 6901
- romani (1996) Romani, R.W., 1996, ApJ, 470, 469
- romani (2001) Romani, R.W., Johnston, S., 2001, ApJ Letter, 557, 93
- shibata (1995) Shibata, S., 1995, MNRAS, 276, 537
- shibata (1991) Shibata, S., 1991, ApJ, 378, 239
- smith (2001) Smith, I.A., Michel, F.C., Thacker, P.D., 2001, MNRAS, 322, 209
- spitkovsky (2006) Spitkovsky, A., 2006, ApJL, 648, 51
- taka (2004) Takata, J., Shibata, S., Hirotani, K., 2004, MNRAS, 354, 1120
- taka (2006) Takata, J., Shibata, S., Hirotani, K., Chang, H.-K., 2006, MNRAS, 366, 1310
- taka (2007) Takata, J., Chang, H.-K., 2007, ApJ, 670, 677
- taka (2008) Takata, J., Chang, H.-K., Shibata, S., 2008, MNRAS, 386, 748
- taka (2009) Takata, J., Chang, H.-K., 2009, MNRAS, 392, 400
- taka (2010) Takata, J., Wang, Y., Cheng, K.S. 2010, ApJ, 715, 1318
- taka (2011) Takata, J., Wang, Y., Cheng, K.S., 2011, MNRAS, 414, 2173
- taka (2012) Takata, J., Cheng, K. S., Taam, Ronald E., 2012, ApJ, 745, 100
- vebter (2012) Venter, C., Jphnson, T.J., Harding, A.K., 2012, ApJ, 744, 33
- vigano (2015) Vigan, D., Torres, D.F., 2015, MNRAS, 449, 3755
- wada (2011) Wada, T, Shibata,S, 2011, MNRAS, 418, 612
- wada (2005) Wada, T, Shibata,S, 2007, MNRAS, 376, 1460
- wang (2014) Wang, Y., Ng, C.W., Takata, J., Leung, Gene C. K., Cheng, K. S., 2014, MNRAS, 445, 604
- wang (2010) Wang, Y., Takata, J., Cheng, K.S., 2010, ApJ, 720, 178
- watters (2011) Watters, K., Romani, R.W., 2011, ApJ, 727, 123
- yakovlev (2004) Yakovlev, D. G., Pethick, C. J., 2004, ARA&A, 42, 169
- yuki (2012) Yuki, S., Shibata, S., 2012, PASJ, 64, 43
- zavilin (2007) Zavlin, V.E., 2007, Ap&SS, 308, 297