Abstract

We have performed a Monte-Calro simulation for Galactic population
of pulsars and for the γ-ray observations.
We apply two-layer outer gap model, which has been developed
by Wang, Takata & Cheng, for
the γ-ray emission process, and study the radiation characteristics
as a function of the magnetic inclination angle (α)
and the Earth viewing angle (ζ). In our model,
the γ-ray flux and the spectral
cut-off energy tend to decrease as the inclination and
viewing angles deviate from 90∘.
The emerging spectrum above 100 MeV becomes soft with
a photon index p∼1.8−2 for ζ→90∘ and
p∼1.2−1.3 for ζ≪90∘.
Our simulation predicts that
the pulsars with larger inclination angles (α=70−90∘) and
larger viewing angles (ζ=70−90∘) have been preferentially detected
by the Fermiγ-ray telescope,
and hence the observed pulse profiles of the
γ-ray pulsars have the double peak structure rather
than single peak. In the simulation,
most γ-ray millisecond
pulsars are categorized as the radio-quiet γ-ray pulsars, because
its radio fluxes are under the sensitivities of the
major radio surveys. Even we drastically increase the radio
sensitivity by a factor of ten, the number of radio-selected
millisecond pulsars detected by the Fermi ten years observations is still
much less than the expected γ-ray-selected millisecond pulsars,
indicating the
radio-quiet millisecond pulsars must contribute to the Fermi unidentified
sources and/or the γ-ray background radiations.
We argue that γ-ray pulsars observed with
a smaller viewing angle (ζ≪90∘)
will appear as low-efficient γ-ray
pulsars. For example, unique radiation properties of
the low-efficient γ-ray pulsar, PSR J0659+1414,
can be explained by the present our gap model with a viewing geometry of
α∼ζ=40∘−50∘.

keywords:

The Large Area Telescope (LAT) on board the Fermiγ-ray telescope
has increased number of γ-ray pulsars, and the recent Fermi catalog includes more that 60 γ-ray pulsars, which includes
9 millisecond pulsars (Abdo et al. 2009a,b, 2010a3; Saz Parkinson et al. 2010).
Furthermore, the detection of radio millisecond pulsars associated with
about 20 unidentified Fermi point sources (e.g. Ray 2010; Caraveo 2010; Ranson et al. 2011; Keith et al. 2011) has been reported,
suggesting that the millisecond pulsar, as well as the canonical pulsar,
is one of the major Galactic γ-ray source. It can be expected that
more γ-ray pulsars will be added to the list over the Fermi mission.
The spectral shape and the pulse morphology measured by the
Fermi have been used to
discriminate particle acceleration and γ-ray emission models; the
polar cap model (Ruderman & Sutherland 1975; Daugherty & Harding 1982, 1996),
the slot gap model (Arons 1983; Muslimov & Harding 2004; Harding et al. 2008;
Harding & Muslimov 2011) and the outer gap model
(Cheng, Ho & Ruderman 1986a,b; Hirotani 2008; Takata, Wang & Cheng 2010a).
The polar cap model assumes the acceleration region near the stellar surface,
and the slot gap and outer gap models assume the emission region extending
to outer magnetosphere. The cut-off features of the γ-ray spectra
of the Crab and the Vela pulsars measured
by Fermi imply that the γ-ray
emission site of the canonical pulsars
is located in the outer magnetosphere rather than
near polar cap region, which produces a cut-off feature
steeper than the observed one (Aliu et al. 2009; Abdo et al. 2009c, 2010d).
Romani & Watters (2010) and Watters & Romani (2011) have studied
morphology of the pulse profiles of the young pulsars predicted
by the our gap and slot gap models, and they argued statistically
that the outer gap geometry is more
consistent with the Fermi observations than the slot gap model.
Venter, Harding & Guillemot (2009)
found that the observed pulse profiles of several millisecond pulsars
detected by the Fermi cannot be explained by the outer gap and/or the slot gap models, and
proposed a pair-starved polar cap model, in which the particles are
continuously accelerated up to high altitude because of the insufficient
multiplicity of the pairs.

The high quality data measured by the Fermi
enables us to perform a detail study
for population of the γ-ray pulsars.
Takata et al. (2010a) have studied the relation between
the emission properties (luminosity and spectral cut-off energy) and pulsar
characteristics (e.g. rotation period and magnetic field). They proposed
that the outer gap accelerator model controlled by
the magnetic pair-creation process can explain the observed population
statistics better than that controlled by the photon-photon
pair-creation process (Zhang & Cheng 1997, 2003).
Wang, Takata & Cheng (2010) have fitted the observed phase-averaged
spectra by using a two-layer outer gap model, in which the accelerator
consists of a wide but low charge density main region and a narrow
but high charge density
screening region. They suggested that the relation between the γ-ray
luminosity (Lγ) and the spin down power (Lsd)
can be expressed as Lγ∝Lβsd with β∼0 for
Lsd≥1036 erg/s, while β∼0.5
for Lsd≤1036 erg/s. This relation is consistent with
the theoretical expectation, i.e. the gap fractional size f
is determined by f=min(fzc,fm), where fzc and fm
are the gap size determined by photon-photon pair-creation process
(Zhang & Cheng 1997, 2003) and magnetic pair-creation process
(Takata et al 2010a), respectively (c.f. section 3.2).
Equating fzc and fm corresponds to Lsd∼1036 erg/s.

The study for population synthesis
of the radio pulsars has been developed by using the detailed modeling
of the radio emissions and the radio surveys
(e.g. Bailes & Kniffen 1992; Sturner and Demer 1996;
Faucher-Gigu`ere and Kaspi, 2006).
Also the population study of the γ-ray pulsars have
been developed by several authors (e.g.
Gonthier et al. (2002) for the polar cap model, Cheng & Zhang
(1998) for the outer gap model, and Story, Gonthier & Harding (2007)
for the slot gap model). For example,
Story et al. (2007) studied the
population of γ-ray millisecond pulsars with
the slot gap accelerator model, and predicted the Fermi observations.
They predicted that the Fermi will detect 12 radio-loud and
33-40 radio-quiet γ-ray millisecond pulsars.
With the Monte-Calro simulation of the outer gap,
Takata, Wang & Cheng (2011a,b) have explained
the observed distributions of the characteristics
of the γ-ray pulsars detected by the Fermi with the
six-month long observations. They predicted that at least 80 γ-ray
canonical pulsars have a γ-ray flux exceeding
the sensitivity of the six-month Fermi observations, suggesting
the present observations have missed many γ-ray emitting pulsars.
Watters & Romani (2011) simulated
the Galactic distribution of the young γ-ray pulsars, and compared
the simulated pulse morphology with the Fermi results.
They argued statistically that the outer gap model explains
the distributions of
the pulse morphology (e.g. the phase-separation of the two peaks)
measured by the Fermi better than the slot gap model.
The population studies (e.g. Kaaret & Philip 1996;
Cottam, Jean Faucher-Gigu`ere
& Loeb 2010; Takata et al. 2011b) have also pointed
out that unidentified pulsars, in particular the millisecond pulsars located
at high-Galactic latitudes will associate with the Fermi
unidentified sources (Abdo et al. 2010b),
and will contribute to the γ-ray background radiations.

In the our previous studies (Takata et al. 2011a,b), we
ignored the dependence of the radiation characteristics
on the viewing geometry (i.e. the magnetic inclination angle and
the Earth viewing angle measured from the pulsar’s rotation axis),
and focused on the distributions of the
properties (e.g. the rotation period and magnetic field)
of the γ-ray pulsars. However,
the observed radiation characteristics, such as the flux, spectral
cut-off energy and pulse profile,
must be affected by the viewing geometry. To perform a more solid study on
the population, which is compared with the high-quality Fermi data,
it is required a three-dimensional model that takes into account
the dependence of the radiation characteristics on the viewing geometry.

In this paper, we develop a Monte-Carlo study for the population
of the γ-ray pulsars with the radiation model that takes into account
the dependence on the viewing geometry.
In section 2, we review the Monte-Carlo simulation for
the population of the γ-ray pulsars. We discuss our outer gap model
in section 3.
In section 4.1, we discuss the
dependence of the radiation characteristics on the viewing geometry.
In section 4.2, we compare the results of the Monte-Calro simulation
with the Fermi six-month long observations.
We also show
the expected population of the γ-ray pulsars if the Fermi
observations continue five-years or ten-years. In section 5,
after summarizing our simulation results,
we discuss the viewing geometry of PSR J0659+1414, which is
known as a low efficient γ-ray pulsar.

In this paper, we denote the canonical pulsar and the millisecond pulsar
as CP and MSP, respectively.
We assume that the birth rates of the CPs and the MSPs
are ∼10−2yr−1 and
10−6∼10−5yr−1 (Lorimer et al. 1995, Lorimer 2008),
respectively, and we ignore the millisecond pulsars in globular clusters.
The birth location
is determined by the spatial distributions given by
(Paczynski 1990),

ρR(R)=aRe−R/RexpRR2exp,

ρZ(Z)=1Zexpe−|Z|/Zexp,

(1)

where R is the axial distance from the axis through the Galactic
centre perpendicular to the Galactic disk and Z is the distance from
the Galactic disk,
Rexp=4.5 kpc and
aR=[1−e−Rmax/Rexp(1+Rmax/Rexp)]−1 with Rmax=20
kpc. In addition, we apply Zexp=75 pc for the CPs
and Zexp=200 pc for the MSPs, respectively.

To obtain current position of each simulated pulsar, we
solve the equation of motion from its birth to the current time.
The equation of motion is given by

dR2dt2=v2ϕR−∂Φtot∂R,

(2)

dZ2dt2=−∂Φtot∂Z,

(3)

and

Rvϕ=constant.

(4)

Here vϕ is the azimuthal component of the velocity, Φtot=Φsph+Φdis+Φh is the total gravitational potential, where
Φsph, Φdis and Φh are spheroidal, disk
and halo components of the Galactic gravitational potential, and are
given by

Φi(R,Z)=−GMi√R2+[ai+(Z2+b2i)1/2]2,

(5)

where i=sph and dis, asph=0, bsph=0.277 kpc,
Msph=1.12×1010M⊙, adis=3.7 kpc, bdis=0.20 kpc,
and Mdis=8.07×1010M⊙, while for the halo
component

Φh(r)=−GMcrc[12ln(1+r2r2c)+rcrtan−1(rrc)],

(6)

where rc=6.0 kpc and Mc=5.0×1010M⊙ (c.f. Burton & Gordon 1978; Binney & Tremaine 1987; Paczynski 1990).
The Lagrangian in units of energy
per unit mass is given by

L=v2(R,Z,ϕ)2−Φtot(R,Z),

(7)

where v is the velocity.

For the initial velocity of each modeled pulsar,
we assume random isotropic direction of the velocity,
and assume magnitude drawn from
a Maxwellian
distribution with a characteristic width of
σV=265 km/s for the CP and σV=70 km/s for
the MSP (c.f. Hobbs et al. 2005), namely,

ρV(V)=√π2V2σ3Ve−V2/2σ2V.

(8)

2.1 Pulsar characteristics

In this section, we describe how we calculate the current value of the various
characteristics of the simulated pulsars.

Canonical pulsars

For canonical γ-ray pulsar,
which is born soon after the supernova explosion,
the true age is more or less equal to the spin down age τ=P/2˙P,
where P and ˙P are the rotation period and its time derivative.
This allows us to calculate the present distributions of the pulsar
characteristics from the initial distributions.

For the Crab pulsar, the initial period is estimated to be
P0∼19ms. A short birth rotation period is also
expected for the young 16ms pulsar PSR J0537-0691 (Marchall et al. 1998),
suggesting most pulsars were born with P0∼20−30 ms.
There is a good evidence
that some pulsars were born with a longer initial period (e.g.
P0∼62 ms of PSR J1911-1925, Kaspi et al. 2001).
However, we found that the distribution of initial rotation
period does not affect
much to the results of the following Monte Carlo simulation.
For example, we compared the results of two simulations with different
distribution of the initial period; one is that all simulated pulsars are
randomly distributed in the range P0=20−30 ms, and other is that
70 % of pulsars are distributed at P0=20−30ms and
30 % are P0=20−100 ms. In such a case, we could not see any significant
difference in the populations of the simulated γ-ray pulsars.
The difference in the simulated distributions of the γ-ray pulsars
becomes significant if we assume that >50% of new born pulsars
is distributed with the initial period of P0≥30 ms.
Because we do not know well the exact
distribution of the initial period, we randomly choose
the initial period in the narrow range P0=20−30 ms,
which provides a more consistent result of our simulation
with the Fermi observations.

We assume a Gaussian distribution in log10Bs for the initial
distribution of the stellar magnetic field measured at the magnetic equator,

ρB(log10Bs)=1√2πσBexp[−12(log10Bs−log10B0σB)2].

(9)

In this study, we apply log10B0=12.6 and σB=0.1.
Because the canonical γ-ray pulsars are younger than 10 Myr, we ignore
evolution of the stellar magnetic field, which may be important
for the neutron star with an age older
than 10 Myr (Goldreich & Reisenegger 1992; Hoyos, Reisenegger & Valdivia, 2008).

With a constant stellar magnetic field in time, the period evolves as

P(t)=(P20+16π2R6sB2s3Ic3t)1/2,

(10)

where Rs=106 cm is the
stellar surface and I=1045gcm2 is the neutron star momentum
of inertial. In addition, we artificially assumed that the magnetic
inclination angle does not evolve with the spin down age,
while a decreases on the spin down time scale has been pointed out
(e.g. Davis & Goldstein 1970; Michel 1991; Tauris & Manchester 1998).
The time derivative of the rotation period is calculated
from

˙P(t)=8π2R6sB2s3Ic3P.

(11)

Millisecond pulsar

We assume that all MSPs are born through the so called recycled process,
in which the accretion of the matter from the low mass companion star
spins up the neutron star. It implies that the true age of the binary system
is different from the spin down age of the MSP.
However, we expect that the Galactic distribution does not depend
on the spin down age of the MSPs. With the typical velocity of the
observed MSPs , V∼70 km/s, it is expected that the
displacement of the MSPs (or binary system)
with the typical age, ≥100 Myr,
becomes larger than the size of the Galaxy. With the relatively
slow velocity, V≤100 km/s, however,
the MSPs remain bound in the Galactic potential and hence
their Galactic distribution does not depend on the spin down age.

The initial rotation period of MSP is related with the
history of the accretion process after the decaying stage of the magnetic
field (Campana et al. 1999; Takata, Cheng and Taam 2010b).
However, the description of the
transition from an accretion powered to the rotation powered phase is
not well understood due to the complexities in the description
of the interaction between the magnetosphere
of a neutron star and its accretion disk (Romanova et al. 2009).
Furthermore, the true age of the MSP after the supernova explosion
in the binary system is different from the spin down age of the MSP.
These theoretical uncertainties make it difficult
to obtain the present distributions of the MSP
properties from the initial distributions.

In this paper, we assign the “current”
pulsar properties for each
simulated MSP, instead of modeling from the initial distributions; namely,
we (1) randomly select the age of the simulated
MSP up to 10 Gyr, (2) shifts the simulated MSP from its birth location
to the current location, and (3) assign the parameters of the MSP
following the observed distributions.
We assign the period time derivative (˙P)
and the stellar magnetic field (Bs) following the observed
˙P−Bs distribution (Manchester et al. 2005).
From the assigned period time derivative
and the stellar magnetic field, the current rotation period and the
spin down age are calculated from

P−3=0.97B28˙P−1−20ms

(12)

and

τ=1.5×109P−3˙P−120yr

(13)

respectively (Lyne & Graham-Smith, 2006).
Here B8 is the stellar magnetic field in units of 108 G,
P−3 and ˙P−20 are the rotation period in units
of 1 millisecond and its time derivative in units of 10−20, respectively.
In fact, above process can provide
a consistent Galactic distribution of the radio MSPs with the
observations (c.f. Takata et al. 2011b).

2.2 Radio emissions

Using the empirical relation
among the radio luminosity, rotation period, and period time derivative,
the distribution of the
radio luminosity at 400 MHz is expressed by (Narayan & Ostriker 1990)

ρL400=0.5λ2eλ,

(14)

where λ=3.6[log10(L400/<L400>)+1.8] with
<L400>=η106.64˙P1/3/P3, and
L400 is the luminosity in units of mJykpc2.
Here η is a scaling factor to adjust the observed distribution,
and η=1 (or 0.05) for the CPs (or MSPs).
The radio flux on Earth is given by S400=L400/d2,
where d is the distance to
the pulsar. We scale the simulated 400 MHz luminosity
to the observational frequency using
a typical photon index ∼2 for the CPs and ∼1.8 for the MSPs,
respectively (Kramer et al. 1997, 1998).

We take into account the beaming of the radio emission.
For the CPs, we apply the half-angle, which is measured from the magnetic axis,
of the radio cone studied by Kijak and Gil (1998, 2003),

ωCP∼1∘.24r1/2KGP−1/2,

(15)

where

rKG=40ν−0.26GHz˙P0.07−15P0.3,

where ˙P−15 is the period time derivative in units of 10−15, and
νGH is the radio frequency in units of GHz.
For the MSPs, the half-angle does not depend on the frequency and
is approximately described as (Kramer & Xilouris 2000),

ωMSP∼ω0(P/1s)−1/2,

(16)

where ω0 is randomly chosen in the range 2.75∘−5.4∘.
The radio emission can be detected by observer with a viewing angle
between max(0∘,α−ωi)
and min(α+ωi, 90∘),
where i=CP or MSP, and α is the inclination angle
between the rotation axis and the magnetic axis.

We use the ten radio surveys (Molongo 2, Green Band 2 and 3, Arecibo 2 and 3,
Parkes 1, 2 and MB, Jordell Bank 2 and Swinburne IL), whose
system characteristics are listed in table 1 of Takata et al (2011a)
and the references therein. To calculate the dispersion measure, we apply
the Galactic distribution of electrons studied by Cordes & Lazio (2002).

The pulsar rotation energy, which is thought to be essential energy source
of the γ-ray emission, can be released by both the current
braking torque and
the magnetic dipole radiation. According to the analysis
of the force-free magnetosphere done by Spitkovsky (2006),
the spin down power depends
on the inclination angles as Lsd∝1+sin2α; in other words,
Lsd changes only by a factor of 2 with the inclination angle.
In the present Monte-Carlo study, therefore, we ignore the dependence of
the spin down power on the inclination angle, and apply the conventional
expression that Lsd=4(2π)4B2sR6s/6c3P4.

3.1 Two-layer outer gap model

The outer gap dynamics is controlled by the pair-creation process, which
produces a charge distribution in the trans-field
direction (Cheng, Ho & Ruderman
1986a,b; Takata, Shibata & Hirotani 2004).
Wang et al. (2010, 2011) argued that the
outer gap should be approximately
divided into two layers, i.e. the main acceleration region
starting from the last-open field lines and the screening region lying at
the upper part of the gap. In the main acceleration region,
the charge density is ∼10%
of the Goldreich-Julian value (Goldreich & Julian 1969), and
a strong electric field accelerates the electrons and positrons up to
the Lorentz factor of Γ∼107.5.
The accelerated particles emit several GeV photons
via the curvature radiation process. In the screening region, the large
number of pairs created by the pair-creation process starts to screen out
the gap electric field.
The curvature radiation from the screening pairs produces mainly ∼100 MeV
photons.

A simple description of the electric field structure inside
the two-layer outer gap is discussed in Wang et al. (2010, 2011).
We denote x, z and ϕ
as the coordinates along the magnetic field line,
perpendicular to the magnetic field line in the poloidal plane and in the
magnetic azimuth, respectively. We expect that the particle number
density increases exponentially near the boundary (z=h1)
between the main acceleration and screening regions (Cheng et al. 1986a,b),
and that the charge density is almost
constant in the screening region (Hirotani 2006). Hence, we approximately
describe the distribution of the charge density in the z-direction with
the step function as follows

ρ(→r)={ρ1(x,ϕ),$if$0≤z≤h1(x,ϕ),ρ2(x,ϕ),$if$h1(x,ϕ)<z≤h2(x,ϕ),

(17)

where,|ρ1|<|ρGJ|<|ρ2|,
z=0 and z=h2 correspond to the last-open field line and
the upper boundary of the gap, respectively. For simplicity,
we define the boundary between main acceleration region and the screening
region, i.e. h1, with a magnetic field line
with such approximation that h1/h2 is constant along the
magnetic field line. The present model
predicts that the charge density in the screening region should be proportional
to the Goldreich-Julian charge density (Wang et al. 2010). This situation will
be satisfied because of a lot of the pairs created by the
pair-creation process in the screening region. In the main acceleration region,
because the number density is much smaller
than the Goldreich-Julian value, its
distribution along the magnetic field line does not important
for the electric field distribution. Therefore, we
approximate that ρ−ρGJ∼g(z,ϕ)ρGJ(→r) for both
main acceleration and screening regions, where

g(z,ϕ)={−g1(ϕ),$if$0≤z≤h1(x,ϕ),g2(ϕ),$if$h1(x,ϕ)<z≤h2(x,ϕ).

(18)

We assume that g1>0 and g2>0 so that |ρ|<|ρGJ|
for the main acceleration region and |ρ|>|ρGJ| for
the screening region.

To obtain the typical strength of the electric field in the gap,
we find the solution of the Poisson equation
for each azimuthal angle (c.f. Wang et al. 2010, 2011),

∂2∂z2Φ′(x,z,ϕ)|ϕ=fixed=−4π[ρ(x,z,ϕ)−ρGJ(→r)]ϕ=fixed,

(19)

where Φ′ is the electric potential of the accelerating field.
Here we assumed that the derivative of the Potential field in the z-direction
is much larger than that of x-direction, and of ϕ-direction.

In this paper, we neglect
the z-dependence of the Goldreich-Julian charge
density, and approximate the Goldreich-Julian charge density as
ρGJ(x,ϕ)∼−ΩBx/2πcRc (Cheng et al. 1986a,b),
where Ω and Rc are the angular frequency of the pulsar
and the curvature radius of the field line, respectively.
The boundary conditions on the lower (z=0) and upper (z=h2) boundaries
are given by

Φ′(x,z=0,ϕ)=0andΦ′(x,z=h2,ϕ)=0

(20)

respectively. Imposing Φ′ and ∂Φ′/∂z are continuous
at the boundary z=h1, we obtain the solution as

Unknown environment '%'

(21)

where

C1(x,ϕ)=−g1h1(h1−2h2)+g2(h1−h2)2h22,

D2(x,ϕ)=−g1h21+g2h22h22,

and z′≡z/h2(x,ϕ). The accelerating electric field,
E||=−∂Φ′/∂x, is writted as

Unknown environment 'array%

(22)

where we used the relations of the dipole field
that ∂(Bh22)/∂x∼0,
∂z′/∂x=∂(z/h2)/∂x∼0,
∂(h1/h2)/∂x∼0, and approximated that
∂Rc/∂x∼0.

On the upper boundary, we anticipate that the total potential field
(co-rotational potential + non co-rotational potential) in the
gap is continuously connected
to the co-rotational potential field outside the gap. This screening
condition is described by

∂Φ′∂z|z=h2=−E⊥(x,z=h2,ϕ)=0.

(23)

This condition gives the relation between (h1,h2) and (g1,g2) as

(h2h1)2=1+g1g2.

(24)

In this paper, we do not consider the azimuthal distribution of the
dimensionless charge density g1 and g2, because we discuss the general
properties of the γ-ray emissions. The azimuthal structure will be
important for explaining the detailed observed
properties, such as the existing of the third peak in the pulse profile of the Vela
pulsar measured by the Fermi, and the energy dependence of the pulse phase of the third peak
(c.f. Wang et al. 2011).

The typical Lorentz factor of the accelerated particles can be estimated by
force balance between the electric field and the curvature radiation drag
force as

Γ=(3R2c2eE||)1/4.

(25)

The spectrum of the curvature radiation emitted by the individual
particle is written as

Pc(Eγ,→r)=√3e2ΓhRcF(χ),

(26)

where χ=Eγ/Ec with Ec=3hcΓ3/4πRc and

F(χ)=χ∫∞χK5/3(ξ)dξ,

where K5/3 is the modified Bessel function of order 5/3.
A γ-ray spectrum measured on Earth may be expressed by
(e.g. Hirotani 2008)

dFγdEγ∼1d2∑→riN(→ri)Pc(Eγ,→ri)Rc(→ri)△Ai,

(27)

where N∼|ρGJ|/e is the particle number density, →ri
represents the radius to the emission point, from which the emission
is measured by the observer, △Ai is
the area of the calculation grid
in perpendicular to the magnetic field lines. The integrated energy flux
between 100 MeV and 300 GeV can be calculated from

Fγ,100=∫300GeV100MeVdFγdEγdEγ.

(28)

In the gap, the curvature photons are emitted in the direction of the
particle motion, which may be described as (Takata, Chang & Cheng 2007)

→v=vp→B/B+→r×→Ω,

(29)

where the first term represents the motion along the magnetic field line, vp
is calculated from the condition that |→v|=c, and the second term is
the co-rotation motion.
The emission direction measured from the rotation axis (i.e. viewing angle)
ζ and the pulse phase ψ are calculated from (Yadigaroglu 1997)

{cosζ=→n⋅→eΩψ=−ψn−→r/Rlc⋅→n,

(30)

where →n=→v/v, →eΩ is the unit vector in the direction
of the rotation axis, and ψn is the azimuthal angle of the emission
direction.

3.2 Outer gap geometry

In this paper, we adopt a rotating vacuum dipole field as the magnetosphere,
and we assume that a strong emission region extends between
the null charge surface
(→Ω⋅→B=0) of the Goldreich-Julian charge density and
the radial distance r=Rlc, where Rlc=2π/cP is the light
cylinder radius. The electrodynamic studies have pointed out
that the gap current
can shift the inner boundary toward the stellar surface
(e.g. Takata et al. 2004). However, because it is expected
that the curvature radiation below the null charge surface
appears with a emissivity
much smaller than that above the null charge surface (e.g. Hirotani 2006),
we ignore its contribution to the calculation.

We define the fractional gap thickness measured on
the stellar surface as,

f≡h2(Rs,ϕ)rp(ϕ),

(31)

where rp is the polar cap radius. Note because the electric field E||
is proportional to Bh22, it can be found that E||∝f2.

Zhang & Cheng (1997, 2003) have argued a self-consistent outer gap model
controlled by the photon-photon pair-creation process between
the curvature photons and the X-rays from the stellar surface.
They estimated the gap fraction as

fzc,CP=h2(Rs,ϕ)rp(ϕ)∼D⊥(Rlc)Rlc=5.5(P/1s)26/21(Bs/1012G)−4/7

(32)

for the CPs and

Extra open brace or missing close brace

(33)

for the MSPs. Here δr5 is the distance (in units of 105 cm)
from the stellar surface to the position where the local magnetic
field is comparable to
the dipole field, and it will be δr5∼1−10 cm.

We note that Zhang & Cheng (1997, 2003) estimated the gap fraction
by a completely vacuum electric field
E||=ΩBf2R2lc/cRc.
With the same gap fraction,
the solution described by equation (22)
gives an electric field at least a factor of four smaller
than that used in Zhang & Cheng (1997, 2003). This difference can
be important for the typical energy of the curvature radiation (Ec),
because Ec∝E3/4||. In other words,
if we derive the gap fraction from the pair-creation condition that
EXEc=(mec2)2, where EX is the X-ray photon energy,
the present model predicts a fractional gap thickness
larger than that of Zhang & Cheng (1997, 2003). In this paper,
reducing by a factor of four for the electric field in the model of
Zhang & Cheng (1997, 2003), we apply the gap fraction increased by a factor of
43/7∼1.8 from those in equations (32) and (33).

Takata et al. (2010a) proposed that the outer gap model can be
controlled by the magnetic pair-creation process taken place
near the stellar surface.
They argued that half of particles created in the gap or injected
particles into the outer gap at the outer boundary
return to the stellar surface, and these returning particles
emit ∼100 MeV photons near the stellar surface.
A good fraction of 100 MeV photons can make pairs by the magnetic
pair creation process. They argued that if the magnetic
field lines near the surface, instead of
nearly perpendicular to the surface, are bending side-wards due to the
strong local field, the pairs created in these local magnetic field
lines can have a pitch angle, which is defined by the angle between
the directions of the particle’s motion and of the field line,
smaller than 90∘, which results in an
outgoing flow of pairs, and hence the pairs control the size of the outer gap.
We also note that the outgoing flow may be
produced by the Compton scattering. With this model,
the fractional gap thickness is estimated as

fm∼D⊥(Rs)Rp=0.8KP1/2,

(34)

where K∼B−2m,12s7 is the parameter characterizing the
local parameters, i.e., Bm,12 and s7, which
are the local magnetic field
in units of 1012G and the local curvature radius in units
of 107cm, respectively. By fitting the emission characteristic of the
γ-ray pulsars observed by the Fermi,
they estimated as K∼2 for the CPs and
K∼15 for the MSPs.
When the fractional gap thickness fm is smaller (or larger)
than fzc, the magnetic pair-creation (or photon-photon pair-creation)
process controls the gap thickness.

We note that the outer gap should only exist between the last-open
field lines and the critical magnetic field lines that have
the null charge points at the light cylinder.
Figure 1 shows the polar angle,
which is measured from the magnetic axis, of
the polar cap rim (θp, solid line) and of the foot points of the critical field lines
(θc, dashed line) as a function of the magnetic azimuth. The
maximum gap fraction is defined by fmax=(θp−θc)/θp,
and is represented by the dotted-line in Figure 1. The
azimuthal angle ϕ=0 in Figure 1 corresponds to the
plain spanned by the rotation axis and the magnetic axis. In this paper,
we anticipate that the azimuthal expansion of the active gap is limited as
f≤fmax(ϕ) and −90∘≤ϕ≤90∘.

Figure 1: Latitudes of the polar cap rim θp (solid line) and of the
critical magnetic field lines θc (dashed line). The dotted line
shows the maximum gap fraction fmax=(θp−θc)/θp. The
magnetic azimuth ϕ=0 corresponds to the magnetic meridian, which includes
the rotation and magnetic axises. The results are for the inclination angle
α=45∘ and the rotation period P=0.1 s.Figure 2: Phase plot of the photons having energies larger than 100 MeV.
The results are for α=60∘,
B=3×1012 Gauss and fzc=0.1.

In this paper, we present the results of the two-layer outer gap model
by using the ratio of h1/h2=0.95 and
the dimensionless charge density in the main region of 1−g1=0.3.
The value h1/h2=0.95 is chosen because
Wang et al. (2010) fitted the phase-averaged
spectra of mature pulsars by using h1/h2∼0.95. The value 1−g1=0.3
is slightly larger than 1−g1∼0.1 used in Wang et al. (2010), but
reproduces a more consistent simulated population with the Fermi
observations. The gap height h2 and the charge density in screening region
g2 are calculated from equations (31) and (24),
respectively.

4.1 Dependence on the inclination and viewing angles

The dependence of the characteristics of the calculated
γ-ray spectra on the inclination angle and the viewing
angle are summarized in Figures 2∼8.
Figure 2 shows the phase plot of the photons having the energies
larger than 100 MeV. The result is for α=60∘,
B12=3×1012 G and fzc=0.1.

γ-ray flux

The left panels in Figures 3 and 4 show the
γ-ray flux (≥100 MeV) as a function of
the viewing angle ζ and of
the inclination angle α. The vertical line (y-axis) represents
the fractional γ-ray flux,
which is defined by the flux measured by f3Lsd/d2
(c.f. Takata et al. 2011). Here we calculated
the rotation period from Bs=3×1012G and
the gap fraction f=fzc. We can see that
the calculated flux tends to decrease as the line of sight approaches to
the rotation axis, where ζ=0∘.
In the sky map of Figure 2, we see that the large viewing
angle (ζ→90∘) can encounter more
intense emission region, whereas the small viewing angle (ζ≪90∘) will encounter the less intense region
or even miss the emission region.
With the outer gap geometry running from
the null charge surface to the light cylinder,
the observer with a smaller viewing angle (ζ≪90∘)
may miss the emissions from higher latitudes (z≥h2/2)
and measure the only
emission from the lower part of the gap (z∼0).
Because the accelerating electric field
vanishes at the lower boundary of the gap (z=0), the Lorentz factor
of the accelerated particles and the resultant emissivity of the
curvature radiation around lower boundary of the gap are significantly
decreased. Consequently, a smaller viewing angle tends to measure
a smaller flux of the curvature radiation.
In the left panels of Figures 3 and 4,
we also see the tendency that the decrease
of the fractional flux with the decrease of the viewing angle is more gradual
for larger inclination angle. This is because the emission from the outer gap
with larger inclination angle covers more wide region of the sky.
These dependences of the γ-ray flux on the viewing geometry
predict that the Fermi has preferentially detected the pulsars with
larger inclination angles and larger viewing angles near 90∘
(c.f. section 4.2.2).

The right panels in Figure 3 and 4 summarize
the dependence of the fractional γ-ray flux on
the gap fractional thickness. We can see that the factional flux
decreases as the gap fraction increases.
This dependence is related with the expansion
of the outer gap in the azimuthal direction.
In the present calculation, we have assumed that the azimuthal
expansion of the active gap is limited as
f≤fmax(ϕ) and −90∘≤ϕ≤90∘.
As the dotted line in Figure 1 shows, fmax acquires a
maximum value near the magnetic meridian, where ϕ=0∘, and tends to
decrease as the azimuthal angle deviates from ϕ=0. This
implies that the width of the active gap in the azimuthal direction narrows
with the increase of the gap fractional thickness, f. Consequently,
the fractional flux tends to decrease as the gap fraction increases.

Cut-off energy

The cut-off energy is calculated as the peak energy in the spectral
energy distribution of the emission from the main acceleration region.
Figures 5 and 6 show the
dependence of the cut-off energies for the CPs and for the MSPs, respectively,
on the viewing geometry. The vertical axis represents the cut-off energy
measured by 3hcΓ30/4πRlc, where
Γ0=(3R2lcE||,0/2e)1/4 with E||,0=f2B(Rlc).
Figures 5 and 6 show that the cut-off energy
decreases with the decrease of the viewing angle.
As we argued in section 4.1.1,
a large viewing angle ζ∼90∘ can encounter the emission
from the strong accelerating electric field region at the
middle of the gap (z∼h2/2), whereas
the observer with a small viewing angle (ζ≪90∘) measures
the emission from the small electric field region near
the lower boundary (z∼0).
As a result, the cut-off energy in the spectrum decreases with
the decrease of the viewing angle. We can also find in Figures 5
and 6 that the fractional cut-off energy does not depend much
on the fractional gap thickness, f.

Photon index

We applied the minimized-χ2 method to fit the spectrum
between 100 MeV and the cut-off energy with a single power low form.
Figures 7
and 8 show the
dependence of the photon index for the CPs and MSPs, respectively, on
the viewing geometry. Our model predicts that the spectral shape
is relatively soft with a photon index
p∼1.8−2 for larger viewing angle
(ζ→90∘) and hard with
p∼1.2−1.3 for smaller viewing angle (ζ≪90∘).
In Figures 7 and 8,
the transition from soft to hard spectra
occurs in a narrow range of the viewing angle.
In the present two-layer model,
the screening region and the main acceleration region produce
the γ-ray photons with a typical energy of ∼100 MeV and
∼1 GeV, respectively. For the viewing angle closer to
ζ∼90∘, because the observed γ-ray radiation
consists of the emissions from both main acceleration
and screening regions, the emerging spectrum becomes soft with
a photon index of p∼1.8−2 above 100 MeV.
For a smaller viewing angle,
on the other hand, because the emission
of the screening region is missing,
the emission from the only main acceleration region contributes to the spectrum.
In such a case, the emerging spectrum has
a photon index p∼1.2−1.3 above 100 MeV, which closes to
a mono-energetic curvature spectrum.

Figure 3: The γ-ray flux (>100 MeV) measured on the
as a function of the viewing angle. The results are for the canonical
pulsar with Bs=3×1012 G and f=fzc.
Left:The γ-ray flux for
the inclination angle is α=20∘ (solid line),
40∘ (dashed line) and 80∘ (dotted line). The results are
for the gap fraction of f=0.25. Right:The γ-ray flux for
f=0.1 (solid line), 0.2 (dashed line), 0.4 (dotted line) and
0.55 (dashed-dotted line). The results are for the inclination angle of
α=50∘. Figure 4: The same with Figure 3, but for the millisecond pulsars
with Bs=3×108 G. Figure 5: Dependence of spectral cut-off energy on the viewing geometry.
The vertical line represent the cut-off energy in
units of (3/4π)hcΓ30/Rlc, where
Γ0=(3R2lcE||,0/2e)1/4 with E||,0=f2B(Rlc).
The results are for the canonical
pulsar with Bs=3×1012 G and f=fzc.
The lines correspond to same cases as Figure 3.Figure 6: The same with Figure 5, but for the millisecond pulsars
with Bs=3×108 G. Figure 7: Dependence of photon index on the viewing geometry.
The results are for the canonical
pulsar with Bs=3×1012 G. The lines correspond to same cases
as Figure 3.Figure 8: The same with Figure 7, but for the millisecond pulsars
with Bs=3×108 G and f=fzc.

4.2 Results of the Monte-Carlo simulation

six-months

five-years

ten-years

CPs

Nr

Ng

Nr

Ng

Nr

Ng

Ra. Sen. (x 1)

40

39

56

138

59

182

Ra. Sen. (x 2)

51

34

77

123

81

163

Ra. Sen. (x 10)

72

28

130

90

145

120

Beaming

76

28

155

79

177

102

Table 1: Population of simulated radio-selected (Nr)
and γ-ray-selected (Ng) CPs for six-month, five-year and ten-year Fermi observations.
The first line; the results for ten radio surveys
listed in table 1 of Takata et al. (2010a). The second and third lines are
the results
with the sensitivities increased by a factor of two and ten, respectively, and the bottom is the populations associated with only beaming effects
of the radio emission.

six-months

five-years

ten-years

MSPs

Nr

Ng

Nr

Ng

Nr

Ng

Ra. Sen. (x 1)

10

52

14

200

16

284

Ra. Sen. (x 2)

16

48

26

190

29

274

Ra. Sen. (x 10)

45

32

82

152

94

227

Beaming

106

11

321

41

438

62

Table 2: The same with Figure 1, but for the MSPs

In this section, we present the results of the Monte-Carlo simulation, in
which we assume the birth rates of the CPs and of the MSPs are
∼0.015 per year and ∼9×10−6 per year, respectively.
We also assume that the viewing angle and the inclination angle are
randomly distributed. For the inclination angle close to α∼90∘, in which
the null charge points are located on or very close to the stellar
surface, it may cause a numerical error when
we search the last-open field lines of the rotating vacuum field.
To avoid it, we limit the inclination
angle below α≤85∘. Because a comparison
between the simulated and the observed
distributions for the various properties (e.g. rotation
period and magnetic field) of the γ-ray pulsars
was done in Takata et al. (2011a,b),
we focus on the characteristics
of the γ-ray radiations with the viewing geometry.

For the sensitivity of the observations, we refer the Fermi
first pulsar catalog (Abdo et al. 2010a). For the
radio-selected pulsars, the Fermi archived the sensitivity
F∼10−11erg/cm2s
(or ∼3×10−11erg/cm2s) for the galactic latitudes
|b|≥5∘ (or <5∘) with the six-month long observations.
For the CPs, the Fermi six-month data allow us to detect
the pulsed period by blind search,
if the γ-ray flux is larger than
F≥2×10−11 for |b|≥5∘ and F≥6×10−11
for |b|≤5∘. For the MSPs, because there are no
such detections of the γ-ray-selected pulsar so far,
we cannot simulate the Fermi
sensitivity of the blind search. In this paper, therefore, we simulate
the population of the γ-ray-selected MSPs with
the Fermi sensitivity of the blind search of CPs.

Viewing geometry

Figures 9 and 10 plot the inclination angle (α)
and the viewing angle (ζ) for the simulated 100
canonical and millisecond
γ-ray pulsars, respectively. The filled-circles and the boxes represent
the radio-selected and γ-ray-selected γ-ray
pulsars, respectively. For the CPs,
we can see in Figure 9 that the
radio-selected and γ-ray-selected pulsars distribute at
different region in α−ζ plane.
Specifically, the radio-selected CPs group together
around the line α=ζ. This is because we have assumed
that the magnetic axis is the centre of the radio cone,
indicating that the radio emission can be
detected by observer with a viewing angle ζ∼α.
The scattering from the line α=ζ is related with
the width of the radio cone. For example,
if we assume a narrower cone than that
of equation (15), the amplitude of the scattering is reduced.
Figure 9 also shows that the
detected γ-ray pulsars are mainly distributed with
a viewing angle ζ∼70−90∘. These results are consistent with
the results of Watters & Romani (2011).

Unlike with the CPs, the distributions of the radio-selected
and γ-ray-selected MSPs on the α-ζ
plane are overlapped each other. This is because the width of the radio beam
described by equation (16) covers almost whole sky. We note
that most γ-ray-selected MSPs in the simulation
irradiate the Earth with the radio emissions, but the radio fluxes are
lower than the sensitivity of the simulated radio
surveys (c.f. Tables 1 and 2).

In Figures 9 and 10, we can see that
no γ-ray pulsars are detected with a smaller inclination angles (α≪90∘)
and a smaller viewing angles (ζ≪90∘). This is because the γ-ray flux decreases
as the viewing angle and/or the inclination angle decrease, as we discussed
in section 4.1
(c.f. Figures 3 and 4). Hence,
our simulation results predict that the pulsars with
larger inclination and larger viewing angles
(α,ξ→90∘) have
been preferentially detected by the Fermi six-month long observations.

Population

Tables 1 and 2 show the simulated population of the CPs and MSPs, respectively,
with the Fermi six-month,
five-year and ten-year long observations. Here we scale
the sensitivity of the Fermi observations as ∝√T, where
T is the length of the observation time. In addition,
the second lines (“Ra. Sen. (x2)”) and the third lines (“Ra.Sen. (x10)”)
in the tables show the results for the ten
radio-surveys but we increase the sensitivities by a factor of two and of ten,
respectively, and
the fourth lines (“Beaming”) show the population associated with
the only beaming effect of the radio emission.

With the previous radio surveys (first line in Table 1), the present
simulation shows that
40 radio-selected and 39 γ-ray-selected CPs can be
detected by the Fermi six-month observations, indicating that the present
model predicts more γ-ray pulsars than the Fermi observations
(∼20 for both radio-selected and
γ-ray-selected γ-ray pulsars).
To explain the difference between the simulated and observed numbers,
the several reasons will be expected; (1) the Fermi
sensitivity will become much worse at the Galactic plane, and
(2) the γ-ray emissions from
the pulsars will be missed by source confusion with
the complex regions and the unresolved sources that are not modeled in the
diffuse backgrounds. The predicted number 10 of the radio-selected
γ-ray MSPs is consistent with the observed number 9.
The present model predicts 56 (or 59) radio-selected
and 138 (or 182) γ-ray-selected CPs
can be detected by the Fermi with five-year (or ten-year)
observations. For the MSPs, 14 (or 16) radio-selected γ-ray pulsars
will be detected by the Fermi with five-year (or ten-year) observations.

We see in the first lines of Tables 1 and 2 that the simulated
numbers of radio-selected γ-ray CPs and MSPs increase
only ∼20 and ∼10 sources, respectively,
over even ten-year Fermi observations.
This implies that most presently
known radio pulsars (∼2000 for CPs and ∼80 for MSPs)
might not be discovered by the Fermi.
For the γ-ray-selected pulsars, on the other hand,
the simulation predicts that the Fermi can
detect about 140 CPs for about five-year observations.
We note that the predicted number of the radio-loud (or radio-quiet)
γ-ray pulsars really depends on the sensitivities
of radio surveys, as Tables 1 shows; for example,
the predicted number of the γ-ray-selected CPs
after ten-year of the Fermi observations decreases from 182 to 120 if
the sensitivities of the radio surveys increases by a factor of ten.
This indicates that a deep radio search may find more
radio-emissions from the γ-ray-selected pulsars,
such as LAT PSRs J1741-2054 and J2032+4127 (Camilo et al. 2009).
Table 1 also shows if the sensitivities of the radio surveys increase
by a factor of ten, the ratio of the radio-loud and radio-quiet
γ-ray pulsars is almost determined by the beaming effect of
the radio emissions.

For the MSPs (Table 2), most simulated pulsars are
categorized as the γ-ray-selected pulsars
with the previous sensitivities of the radio surveys, although
the Fermi has not confirmed the radio-quiet MSPs.
We argue that it may be difficult to identify radio-quiet
MSPs, because the detection of the rotation period
by the Fermi blind search is much harder than that of the CP.
Furthermore, if the MSP is in binary system, the effects of the orbital
motion on the observed rotation period make it even harder to confirm
the millisecond rotation period by the blind search.
Takata et al. (2011b) discussed that
the γ-ray-selected MSPs in the simulation correspond
to the Fermi unidentified sources located at higher Galactic latitudes.

We note that the radio cones from MSPs are quite huge
so that most MSPs irradiate the Earth with the radio emissions, as
the bottom line in Table 2 shows. This implies that
more radio MSPs associated with the Fermi
unidentified sources will be detected by the future radio surveys,
such as the discovery of new
20 radio MSPs associated with the Fermi unidentified sources (Ray 2010;
Caraveo 2010; Ransom et al. 2011; Keith et al. 2011).
As third line in Table 2 shows, however, even we drastically increase the radio
sensitivity by a factor of 10, the number of radio-selected MSPs
detected by 10 year Fermi observations
can increase from 16 to 94, it is still
much less than the expected 227 γ-ray-selected millisecond pulsars.
Unless the Fermi sensitivity of the blind search is improved,
the most γ-ray MSPs will not be identified and will
contribute to the Fermi unidentified sources and/or the
γ-ray background radiations.

γ-ray radiation characteristics

Figures 11-13 show the
radiation characteristics (γ-ray luminosity, cut-off energy and photon index, respectively) versus the characteristics (the spin down power or
the magnetic field strength at the light cylinder) for the γ-ray
pulsars with the six-month observations. In those figures,
we plot the Fermi data with errors taken from Abdo et al. (2010a)
and Saz Parkinson et al. (2010), and present the simulated pulsars
in the sub-figures. For the simulated pulsars, we randomly
choose 200 simulated
γ-ray pulsars (except for the γ-ray-selected MSPs).

In Figure 11, the γ-ray luminosity is calculated from
Lγ=4πd2Fγ,100, where d is the distance,
and the histograms of the values of the γ-ray
luminosity are projected along the right-hand axis. The solid and
dashed histograms represent the distributions
for the simulated and observed γ-ray pulsars, respectively.
The present model predicts
that most γ-ray CPs have a spin down power of
Lsd∼1035−38erg/s and a γ-ray luminosity
of Lγ∼1034−36, while MSPs have a
Lsd∼1033−35erg/s and
Lγ∼1032.5−34.5erg/s, which are consistent with
the Fermi observations.

In Figure 11, the spin down power Lsd and
the γ-ray luminosity Lγ of the simulated pulsars can be related as Lγ∝Lβsd with β∼0 for Lsd≥1035−36erg/s
and β∼0.5 for Lsd≤1035−36erg/s. In
the present emission model, the γ-ray luminosity is
proportional to Lγ∝f3Lsd as
Figure 3 and 4 indicate (c.f. Takata et al. 2011b).
The change of the slope is caused by switching gap closure process
between the photon-photon pair-creation process and the
magnetic pair-creation process.
As the equations (32)-(34) show, the gap fraction depends
on the rotation
period and the magnetic field as fzc∝P26/21B−4/7
for the photon-photon pair-creation process
and fm∝P1/2 for the magnetic
pair-creation process. These relations imply that
the γ-ray luminosity depends on the spin down power as
Lγ∝L1/14sd for the photon-photon pair-creation process
and Lγ∝L5/8sd for the magnetic pair-creation process.
Equating fzc (32) and fm (34) corresponds to
Lsd∼1035−36erg/s. The change
of the slope β has been found by Wang et al (2010),
who used the two-layer outer gap model to fit the phase-averaged
spectrum of the mature γ-ray pulsars observed by the Fermi.

Although the general trend of the relation between the γ-ray
luminosity and the spin down power
is explained with a simple form Lγ∝Lβsd, we can
see in Figure 11 that
some simulation samples deviate from the relation. For example,
some simulated pulsars with a spin down power of
Lsd∼1036erg/s have a γ-ray luminosity of
Lγ∼1032−33erg/s, which
is about two or three order smaller than the typical
value of Lγ∼1035erg/s.
We emphasize that this low efficiency of the γ-ray emission
is mainly caused by the effects of viewing angle. The
pulsars lying on the relation Lγ∝Lβsd
are observed with viewing angles of ζ∼90∘, whereas
those (“apparently”) low-efficient γ-ray pulsars are
observed with smaller viewing angles.
Although the flux depends on also the inclination angle,
the flux is more sensitive to the viewing angle than the inclination angle,
as Figures 3 and 4 show.
We note that lower efficient γ-ray pulsars tend to locate
closer to the Earth.

Figure 12 shows the cut-off energy versus
the magnetic field at the light cylinder, Blc.
The symbols and the histograms correspond to same cases
as Figure 11. As sub-figure in Figure 12 indicates,
the present simulation predicts that most
γ-ray CPs and MSPs have
a magnetic field at the light cylinder of Blc≥103 G
and Blc≥104 G, respectively, and have a cut-off energy
smaller than ∼2 GeV. These features will be consistent
with the Fermi observations. In the figure, some simulated
CPs have a cut-off
energy significantly smaller than the typical value Ec∼2 GeV.
As well as the case of the γ-ray luminosity,
this deviation from the typical value is cause by the effects of the
viewing angle, as Figures 5
and 6 imply. In Figure 12, our simulation predicts
that the typical cut-off energy (1-1.5 GeV) of the MSPs
is smaller than that (∼2 GeV)
of the CPs. It is difficult to discuss the difference in the
observed cut-off energies between the CPs and the MSPs, because of
the large observational errors.

Figure 13 represents the photon index versus the
spin down power. We can find that the model
distribution of the photon index has
two peaks at p∼1.2−1.3 and p∼1.8−2, and
the observed distribution (dashed
histogram) may also have two peaks at p∼1.3 and p∼1.7.
With the present model, the hard component of p∼1.2−1.3 corresponds
to the spectrum associated with the emission from the only main acceleration
region, while the soft component p∼1.8−2 corresponds
to the spectrum composed of the emissions from
both main and screening regions, as we discussed in section 4.1.3.
The present model does not predict very hard spectrum with a index
p≤1, which has been indicated for some Fermi pulsars.

Figures 14 and 16 represent the averaged
apparent fractional thickness, which is defined by
fa≡(4πd2Fγ/Lsd)1/3, and the photon index, respectively, as a function of the spin down power, and Figure 15 shows
the averaged cut-off energy as a function of the magnetic field at the
light cylinder. In the figures,
the solid and dashed lines represent the results for the
radio-selected and γ-ray-selected CPs, respectively.
In addition, the thick and thin lines correspond to the results
of the simulation and the Fermi observations, respectively.
In Figure 14, our model predicts a tendency that
the apparent fractional thickness, fa, tends to decrease
with the increase of the spin down power. This is because the true
fractional thickness, fzc or fm,
tends to decrease with the increase of the spin down power, that is,
fzc∝L−13/21sdB1/21 from equation (32) and
fm∝L−1/8sdB1/4 from equation (34).
Because the photon index of the spectrum tends to increase
with the decrease of the fractional gap thickness
as Figure 7 shows, the averaged photon index
in Figure 16 increases with the spin down power. We find that these behaviors are qualitatively consistent with the Fermi observations.
In Figure 15, our model predicts that the typical
cut-off energy does not depend much on the magnetic field
strength at the light cylinder, while the present Fermi data
may have a tendency that the cut-off energy increases with the magnetic field
strength at the light cylinder. However,
because the observational errors are
so large, a more deep observation to reduce the errors
may be required to discuss the tendency.

In Figures 17 and 18, we summarize the distributions
of the radiation characteristics for the simulated canonical and millisecond
γ-ray pulsars, respectively, including both radio-selected and
γ-ray-selected pulsars. The solid and dashed lines are results
for the simulated six-month and ten-year Fermi observations, respectively.
Comparing the solid and dashed lines,
we find that the distributions do not depend much on the time span of
the Fermi observations. However, we note that a longer observation
enables to detect γ-ray pulsars with smaller fluxes, which include
γ-ray pulsars with the viewing angle close to the rotation axis
(c.f. Figures 3 and 4).
For the observer with a viewing angle close to the rotation axis, the photon
index tends to be p∼1.2−1.3,
as Figures 7 and 8 show.
Therefore, our model prediction is that
a longer Fermi observation will detect
more γ-ray pulsars with
photon indexes p∼0.12−0.14,
as right panels in Figures 17 and 18
indicate.

Pulse profiles

Figures 19-22 present the calculated pulse
profiles for the different type of the γ-ray pulsars.
For each type of the γ-ray pulsars,
64 samples are randomly chosen to present the model prediction on
the statistical distribution of the morphology of the pulse profiles.
From the left to right panels and from the top to bottom panels in the figures,
the inclination angle increases.
In principle, we can quantitatively compare the simulated
morphology of the pulse profiles in Figures 19-22
(e.g. phase-separation, number of peak)
with the Fermi observations. However, we would like to point out
that it is very difficult to quantify the morphology to compare with the
observations.
In Figure 19, for example, the pulse profile represented
in upper-left conner is indeed a double peak structure with the narrow
phase-separation between two peaks. In the observations, however,
the identification of the double peak structure with such narrow peak
separation really depends on source counts and timing ephemerides
of the pulsars. In Figure 20,
the simulated pulse profile presented
at the seventh-line and the first-column (α=74.2∘ and
ζ=67.5∘) has the first peak much smaller
than the second peak. In the observations, the detection
of such a small peak will depend on the strength of the
background radiation. Moreover,
because the intensity of pulse peak will depend on
the modeling of the gap structure, it is difficult to discuss
the distribution of the
intensity ratio of the first and the second peaks with the present simple
three-dimensional model. In the present paper, therefore, we avoid
a quantitative discussion for the morphology of the pulsed profile.
Instead, we can only present a qualitative comparison.

We note that the peaks emerging in the calculated pulse profiles
are caused by the so-called caustic effect (Romani and Yadigaloglu 1995;
Cheng, Ruderman and Zhang 2000; Dyks, Harding and Rudak 2004),
in which more photons are observed at narrow width of the rotation phase due to
the special-relativistic effects, i.e. the aberration of the
emission direction and photon’s travel time. In such a case,
the pulse phase and the number of the main peak are not affected much by
modeling of the gap structure, whereas
the intensity of the peak depends on the gap structure.

We can see in Figures 19 and 20 that
the pulse profiles of the simulated CPs (in particular,
γ-ray-selected CPs) are described by
double peak structure rather than single peak profile.
This result may explain the tendency that most canonical
γ-ray pulsars observed by the Fermi show the double peak structure.
In the present simulations, we have predicted that the γ-ray
pulsars measured from ζ∼90∘ are
preferentially detected in the simulation (c.f. section 4.1).
With the viewing angle closer to ζ∼90∘, the
calculated pulse profile tends to have the double peak structure.

For the radio-selected canonical γ-ray pulsars,
although the double peak structure are more or less common feature,
single pulse structure or the double peak
structure with a narrow phase separation stands out for pulsars
with inclination angles α∼40∘−50∘, as we can see in
Figure 19. This is because the
simulated radio-selected canonical pulsars
have the viewing angle similar to the inclination
angle (c.f. Figure 9), ζ∼α.
With the present outer gap geometry,
the pulse profile of the
smaller viewing angle (ζ≪90∘) tends to have
single peak or double-peak with a narrow phase separation, as the
phase plot of Figure 2 shows.

Comparing the pulse profiles of the CPs (Figures 19 and 20)
and of the MSPs (Figures 21 and 22),
one may see that the pulse profile
with single peak is more common for the MSPs than the CPs.
In the present model, the fractional gap thickness fZC,MSP
of equation (33) for the MSPs tends
to be larger than fZC,CP of equation (32) for the CPs.
With a larger fractional gap thickness,
the emission from the higher altitude
(larger z in units of the light radius) contributes to
the pulse profile. Because the emission from the higher altitude
produces two caustic peaks with a narrower phase separation, the calculated
pulse profiles of the MSPs have single peak or double
peak structure with the narrow phase separation more than that of the CPs.

In this paper, we have applied the so called two-layer outer gap model for
the γ-ray radiations from the pulsars, and have focused
on the dependence of the γ-ray radiation
characteristics on the inclination angle and the viewing angle. We showed
that the γ-ray flux and the spectral cut-off energy decreases as
the viewing angle deviates from ζ=90∘.
The spectrum above 100 MeV
becomes soft with a photon index p∼1.8−2 for the
observer with a larger viewing angle (ζ→90∘),
whereas it becomes hard with a index p∼1.2−1.3
for the observer with a smaller viewing angle (ζ≪90∘).
The spectrum with a
photon index p∼1.8−2 consists of the emissions
from both main acceleration and screening regions,
while p∼1.2−1.3 corresponds to the emission from the only
main acceleration region.

We have developed the Mote-Carlo simulation
for the population of γ-ray emitting pulsars. The our simulation
predicts that 56 (or 59) radio-selected, 138 (or 182)
γ-ray-selected canonical γ-ray pulsars and 14 (or 16)
radio-selected γ-ray MSPs can be detected
by five-year (or ten-year) Fermi observations.
Even we drastically increase the radio
sensitivity by a factor of ten, the most simulated MSPs are expected as
the γ-ray selected MSPs, which will contribute as the Fermi
unidentified sources and/or γ-ray background radiations.
For the viewing geometry, our simulation predicts
that the radio-selected CPs have the inclination angle (α)
and the viewing angle (ζ) of α∼ζ, and that most
γ-ray pulsars detected by the Fermi have been measured with
a viewing angle of ζ=70∘∼90∘. We demonstrated
that the spin down power Lsd and
the γ-ray luminosity Lγ of the simulated pulsars is
related as Lγ∝Lβsd with β∼0 for Lsd≥1035−36erg/s
and β∼0.5 for
Lsd≤1035−36erg/s. The change of the
slope β is associated with the switching of the gap closure mechanism
from the photon-photon pair-creation process
of Lsd≥1035−36erg/s to
the magnetic pair-creation process of Lsd≤1035−36erg/s.
We showed that the distribution of the photon index
of the simulated γ-ray pulsars
has two peak at p∼1.8−2 and p∼1.2−1.3, which may explain
the observed distribution. We expect
that more γ-ray pulsars with the hard spectrum of p∼1.2−1.3
will be detected by the future Fermi observations.
For the pulse profiles, the present simulation will explain
the observational tendency that
most canonical γ-ray pulsars detected by the Fermi
have the pulse profiles with the double peaks.
The our model expects that single pulse profile or
the double pulse profile with a narrower phase separation is more
common for MSPs than CPs.

The present model can be used to diagnose the viewing geometry
of the γ-ray pulsars. For example, the present model predicts that
most present γ-ray pulsars detected by the Fermi will have
the viewing angle, ζ∼70∘−90∘,
and hence the pulse profile with the double peaks is
more common than that with the single peak. On the other hand, we expect that
PSR J0659+1414, which is known as the extremely low efficient
γ-ray pulsar (Abdo et al. 2010c), has a unique viewing geometry.
PSR J0659+1414 is the radio-loud
γ-ray pulsar and its spin down power is
Lsd=3×1034 erg/s. Interestingly,
the γ-ray emissions from PSR J0659+1414 has been observed
with a γ-ray luminosity,
Lγ∼3×1032erg/s, which
is about two order smaller than the typical value for the pulsars
with a similar spin down power.
To explain the observed low efficiency of the γ-ray radiation,
our simulation predicts that the inclination angle and
the viewing angle are relatively small, say
α∼ζ∼40∘−50∘, and hence the apparent
efficiency is significantly reduced from
the typical value (c.f. Figure 3).
PSR J0659+1414 is also known that the spectral cut-off energy,
Ec∼0.7 GeV, is suggestively smaller than the typical value
of Ec∼2 GeV. This behavior is also able to be explained by the
effects of the viewing geometry, as Figure 5 indicates.
Furthermore, our simulation expects that the pulsars measured as the low efficient γ-ray
pulsars are located closer to Earth.
In fact, PSR J0659+1414 with d∼0.28 kpc is
one of the nearest γ-ray pulsars.
Finally, the observed broad single pulse profile of PSR J0659+1414 is also
expected by the present calculation. As discussed in section 4.2.4 (Figure 19),
the single pulse structure or double peak
structure with a narrow phase separation is more common for
the radio-selected canonical γ-ray pulsar with an
inclination angle α∼40∘−50∘.
On these ground,
the unique radiation properties of PSR J0659+1414 can be explained,
if the Fermi has measured the γ-ray emission from the our gap
with the viewing geometry α∼ζ∼40∘−50∘.
It is expected that the population of the
low-efficient γ-ray pulsars will be increased by the future
Fermi observations. The diagnose its γ-ray efficiency,
pulse profile and cut-off energy
will provide a strong constraint on the emission models.

We thank A.H. Kong, C.Y. Hui, B. Rudak,
Lupin-C.C. Lin, M.Ruderman, R.E. Taam and
S.Shibata for the useful discussions.
We express our appreciation to an anonymous referee for useful
comments. We also thank the Theoretical Institute
for Advanced Research in Astrophysics (TIARA) operated under the Academia
Sinica Institute of Astronomy and Astrophysics, Taiwan,
which enable author (J.T.) to use the PC cluster at TIARA.
KSC is supported by a 2011 GRF grant of the Hong Kong SAR
Government entitled “Gamma-ray Pulsars”.