{'fraction_non_alphanumeric': 0.04374253771844133, 'fraction_numerical': 0.01877781395853685, 'mean_word_length': 4.163077415739332, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Electrodynamics of Outer-Gap Accelerator: Formation of Soft Power-law Spectrum Between 100 MeV and 3 GeV

We investigate a stationary pair production cascade in the outer magnetosphere of a spinning neutron star. The charge depletion due to global flows of charged particles, causes a large electric field along the magnetic field lines. Migratory electrons and/or positrons are accelerated by this field to radiate gamma-rays via curvature and inverse-Compton processes. Some of such gamma-rays collide with the X-rays to materialize as pairs in the gap. The replenished charges partially screen the electric field, which is self-consistently solved together with the energy distribution of particles and gamma-rays at each point along the field lines. By solving the set of Maxwell and Boltzmann equations, we demonstrate that an external injection of charged particles at nearly Goldreich-Julian rate does not quench the gap but shifts its position and that the particle energy distribution cannot be described by a power-law. The injected particles are accelerated in the gap and escape from it with large Lorentz factors. We show that such escaping particles migrating outside of the gap contribute significantly to the gamma-ray luminosity for young pulsars and that the soft gamma-ray spectrum between 100 MeV and 3 GeV observed for the Vela pulsar can be explained by this component. We also discuss that the luminosity of the gamma-rays emitted by the escaping particles is naturally proportional to the square root of the spin-down luminosity.


Introduction
Recent years have seen a renewal of interest in the theory of particle acceleration in pulsar magnetospheres, after the launch of the Compton Gammaray Observatory (e.g., for the Vela pulsar, Kanbach et al. 1994, Fierro et al. 1998; for PSR B1706-44, Thompson et al. 1996;for Geminga, Mayer-Hasselwander et al. 1994, Fierro et al. 1998; for PSR B1055-52, Thompson et al. 1999). The modulation of the γ-ray light curves testifies to the particle acceleration either at the polar cap (Harding, Tademaru, & Esposito 1978;Daugherty & Harding 1982, 1996Sturner, Dermer, & Michel 1995), or at the vacuum gaps in the outer magnetosphere (Cheng, Ho, & Ruderman 1986a,b, hereafter CHR;Chiang & Romani 1992Romani & Yadigaroglu 1995;Higgins & Henriksen 1997. Both models predict that electrons and positrons are accelerated in a charge depletion region, a potential gap, by the electric field along the magnetic field lines to radiate high-energy γ-rays via the curvature process. However, there is an important difference between these two models: An polar-gap accelerator releases very little angular momenta, while an outer-gap one could radiate them efficiently. In addition, three-dimensional outer-gap models commonly explain double-peak light curves with strong bridges observed for γ-ray pulsars. On these grounds, the purpose of the present paper is to explore further into the analysis of the outer-gap accelerator. In the CHR picture, the gap is assumed to be geometrically thin in the transfield direction on the poloidal plane in the sense D ⊥ ≪ W , where D ⊥ represents the typical transfield thickness of the gap, while W does the width along the magnetic field lines. In this limit, the acceleration electric field is partially screened by the zero-potential walls separated with a small distance D ⊥ ; as a result, the gap, which is assumed to be vacuum, extends from the null surface to (the vicinity of) the light cylinder. Here, the null surface is defined as the place on which the local Goldreich-Julian charge density vanishes, where Ω refers to the angular frequency of the neutron star, B z the magnetic field component projected along the rotational axis, and c the velocity of light; s and z refer to the coordinates parallel and perpendicular, respectively, to the poloidal magnetic field. The star surface corresponds to s = 0; s increases outwardly along the field lines. The last-open field line corresponds to z = 0; z increases towards the magnetic axis (in the same hemisphere). If B z > 0 holds in the starward side of the null surface, a positive acceleration field arises in the gap. The light cylinder is defined as the surface where the azimuthal velocity of a plasma would coincide with c if it corotated with the magnetosphere. Its radius from the rotational axis becomes the so-called 'light cylinder radius', ̟ LC ≡ c/Ω. Particles are not allowed to migrate inwards beyond this surface because of the causality in special relativity.
It should be noted that the null surface is not a special place for the gap electrodynamics in the sense that the plasmas are not completely charge-separated in general and that the particles freely pass through this surface inwards and outwards. Therefore, the gap inner boundary is located near to the null surface, not because particle injection is impossible across this surface (as previously discussed), but because the gap is vacuum and transversely thin.
Then what occurs in the CHR picture if the gap becomes no longer vacuum? To consider this problem rigorously, we have to examine the Poisson equation for the electrostatic potential. In fact, as will be explicitly demonstrated in the next section, the original vacuum solution obtained in the pioneering work by CHR cannot be applied to a non-vacuum CHR picture. We are, therefore, motivated by the need to solve self-consistently the Poisson equation together with the Boltzmann equations for particles and γrays. Although the ultimate goal is to solve threedimensional issues, a good place to start is to examine one-dimensional problems. In this context, Hirotani and Shibata (1999a, b, c; hereafter Papers I, II, III) first solved the Boltzmann equations together with the Maxwell equations one-dimensionally along the field lines, extending the idea originally developed for black-hole magnetospheres by Beskin et al. (1992). In Paper I, II, and III, they assumed that the gap is geometrically thick in the transfield direction in the sense D ⊥ ≫ W .
There is one important finding in this second picture: The gap position shifts if there is a particle injection across either of the boundaries (Hirotani & Shibata 2001, 2002a hereafter Papers VII, VIII, IX). For example, when the injection rate across the outer (or inner) boundary becomes comparable to the typical Goldreich-Julian value, the gap is located close to the neutron star surface (or to the light cylinder). In other words, an outer gap is not quenched even when the injection rate of a completely chargeseparated plasma across the boundaries approaches the typical Goldreich-Julian value. Thus, an outer gap can coexist with a polar-cap accelerator; this forms a striking contrast to the first, CHR picture. It is also found in the second picture that an outer gap is quenched if the created particle density within the gap exceeds several percent of the Goldreich-Julian value. That is, the discharge of created pairs is essential to screen the acceleration field.
The purpose of this paper is to examine the second picture more closely. In all the previous works in the second picture, the particle energy distribution has been assumed to be mono-energetic in the sense that the particles attain the equilibrium Lorentz factor at each point, in a balance between the electrostatic acceleration and the radiation-reaction forces. In the present paper, we discard this assumption and explicitly consider the energy dependence of particles by solving the Boltzmann equations for positrons and electrons. We will demonstrate that the particle energy distribution cannot be represented either by a power law or by the mono-energetic approximation. We will further show that a soft power-law spectrum is generally formed in 100 MeV-GeV energies as a result of the superposition of the curvature spectra emitted by particles migrating at different positions.
In the next section, we describe the difficulties of electrodynamics found in the first picture. We then present the basic equations in § 3. and apply the theory to four γ-ray pulsars and compare the predictions with observations in § 4. In the final section, we discuss the possibilities of the unification of our picture with the CHR picture, as well as the unification of the outer-gap and polar-cap models.

Difficulties in Previous Outer-gap Models
To elucidate the electrodynamic difficulties in previous outer-gap models, we have to examine the Poisson equation for the electrostatic potential, Ψ. Since the azimuthal dimension is supposed to be large compared with D ⊥ in conventional outer-gap models, the Poisson equation is reduced to the following twodimensional form on the poloidal plane The true charge density, ρ e , is given by ρ e = ρ + + ρ − , where ρ + and ρ − represent the positronic and electronic charge densities, respectively. In the CHR picture, the magnetic field is supposed to have a singlesigned curvature. As a result, the particle number density (ρ + − ρ − )/e grows exponentially in z direction. Because of this exponential growth of the particle number density, it has been considered that most of the γ-rays are emitted from the higher altitudes (i.e., large z regions) in the gap.
To explain the observed γ-ray luminosities with a small D ⊥ , one should assume that the created current density becomes the typical Goldreich-Julian value in the higher altitudes. That is, the conserved current density per magnetic flux tube with strength B, should satisfy in the order of magnitude. However, such a copious pair production will screen the local acceleration field, −∂Ψ/∂s, as the Poisson equation (2) indicates. This screening effect is particularly important near to the inner boundary. Without loss of any generality, we can assume −∂Ψ/∂s > 0. In this case, because of the discharge, only electrons exist at the inner boundary. (We may notice here that external particle injections are not considered in the CHR picture.) Thus, we obtain ρ e /B = ρ − /B ∼ −Ω/(2πc) in the order of magnitude. In the vicinity of the inner boundary, we can Fourier-analyze equation (2) in z direction to find out that the −∂ 2 Ψ/∂z 2 term contributes only to reduce ∂(−∂Ψ/∂s)/∂s, the s gradient of the acceleration field. Thus, a positive −ρ GJ must cancel the negative ρ e to make the right-hand side be positive. That is, at the inner boundary, must be satisfied, so that the acceleration field may not change sign in the gap. It follows that the polar cap, where B z ∼ B holds, is the only place for the inner boundary of the 'outer' gap to be located, if the created particle number density in the gap is comparable to the typical Goldreich-Julian value (eq. [3]). Such a non-vacuum gap must extend from the polar cap (not from the null surface where ρ GJ vanishes) to the light cylinder. We can therefore conclude that the original vacuum solution obtained by CHR cannot be applied to a non-vacuum CHR picture when there is a sufficient pair production that is needed to explain the observed γ-ray luminosity.
To construct a self-consistent model, we have to solve equation (2) together with the Boltzmann equations for particles and γ-rays. We formulate them in the next section.

Poisson equation
Neglecting relativistic effects, and assuming that typical transfield thickness, D ⊥ , is greater than or comparable to the longitudinal width, W , we can reduce the Poisson equation for non-corotational potential Ψ into the one-dimensional form (Hirotani 2000b, Paper VI;see also § 2 in Michel 1974) As described at the end of § 3 in Paper VII, it is convenient to introduce the Debye scale length c/ω p , where and B in is the magnetic field strength at the inner boundary; e designates the magnitude of the charge on an electron, m e the rest mass of an electron. Thus, we can introduce the following dimensionless coordinate variable: where Ω 2 ≡ Ω/(10 2 rad s −1 ). By using such dimensionless quantities, we can rewrite the equation (5) into and the particle distribution functions are defined by N + (s, Γ) and N − (s, Γ) represent the distribution functions of positrons and electrons, respectively, at position s and Lorentz factor Γ. We evaluate the dimensionless Goldreich-Julian charge density B z /B in in equation (9) at each s, by using the Newtonian dipole field (see eqs. [9]- [12] in Paper VIII for details).

Particle Boltzmann Equations
On the poloidal plane, particles migrate along the magnetic field lines. Therefore, at time t, the distribution function N of particles obeys the following Boltzmann equation, where p refers to the particle momentum, F ext the external forces acting on particles, and S the collision terms. In the present paper, F ext consists of the Lorentz and the curvature radiation reaction forces.
Since the magnetic field is much less than the critical value (4.41 × 10 13 G), quantum effects can be neglected in the outer magnetosphere. As a result, curvature radiation takes place continuously and can be regarded as an external force acting on a particle. If we instead put the collision term associated with the curvature process in the right-hand side, the energy transfer in each collision would be too small to be resolved by the energy grids (see description below eq. [21] for details). We take the γ-ray production rate due to curvature process into account consistently in the γ-ray Boltzmann equations.

Stationary Boltzmann Equations
Imposing a stationary condition and neglecting the pitch-angle dependence, we can reduce equation (12) to (Appendix A) where Ω p designates the angular velocity of particles due to E ×B p drift. Since the drift motion due to the gradient and the curvature of B p is negligible for typical outer-gap parameters, Ω p coincides Ω provided that B φ = 0 and |E | ≪ |E ⊥ | = B p ̟/̟ LC hold, where ̟ refers to the distance from the rotation axis. Radiation-reaction force due to the curvature radiation, P CV /c, is given by where ρ c (s) refers to the curvature radius of the magnetic field line. For a justification of adopting the pure-curvature formula, see § 5.1.

Collision terms
We assume in this paper that γ-rays are either outwardly or inwardly propagating along the local magnetic field lines. For simplicity, we neglect the deviation of the γ-rays from the field lines due to magnetic curvature. In what follows, we assume that the soft photons are emitted from the neutron star and hence unidirectional at the gap. Then the cosine of the collision angle µ has a unique value µ + or µ − , for outwardly or inwardly propagating γ-rays, respectively; µ + and µ − are determined by the magnetic inclination α i at each position ξ. Under these assumptions, the source term can be expressed as where the dimensionless γ-ray distribution function g ± are defined as G ± refers to the number of γ-rays per unit volume per unit dimensionless energy ǫ γ = hν γ /m e c 2 . We will explicitly define the differential redistribution function for the pair production, ∂η p /∂Γ, by equation (30), after briefly describing the soft photon field in equations (23)-(29). If we multiply dΓ on both sides of equation (17) or (18), the first (or the second) term in the right-hand side represents the rate of particles disappearing from (or appearing into) the energy interval m e c 2 Γ and m e c 2 (Γ + dΓ) due to inverse-Compton (IC) scatterings; the third (or the fourth) terms does the rate of particles created via collisions between the outgoing (or the ingoing) γ-rays and ambient soft photons. The positrons (or electrons) are supposed to collide with the soft photons at the same angle as the outwardly (or inwardly) propagating γ-rays; therefore, the same collision angle cos −1 µ + (or cos −1 µ − ) is used for both IC scatterings and pair production in equations (17) and (18). The collisions tend to be head-on (or tailon) for inwardly (or outwardly) propagating γ-rays, as the gap approaches the star.
The IC redistribution function η γ IC (ǫ γ , Γ, µ) represents the probability that a particle with Lorentz factor Γ upscatters photons into energies between ǫ γ and ǫ γ + dǫ γ per unit time when the collision angle is cos −1 µ. On the other hand, η e IC (Γ i , Γ, µ) describes the probability that a particle changes Lorentz factor from Γ i to Γ in a scattering. Thus, energy conservation gives We introduce the following dimensionless IC redistribution functions, where the particle Lorentz factor is discretized as Γ = b l ; the γ-ray energy is divided into −l 0 + 1 + l m = 15 + l m bins as We normally use l m = 256 and ∆b = 1.5 × 10 5 . In general, η γ IC ±,i (l) can be defined by the soft photon flux dF s /dǫ s and the Klein-Nishina cross section σ KN as follows: where β ≡ 1 − 1/Γ 2 is virtually unity, Ω γ the solid angle of upscattered photon, the asterisk denotes the quantities in the electron (or positron) rest frame. In the rest frame of a particle, a scattering always takes place well above the resonance energy. Thus, the classical formula of the Klein-Nishina cross section can be applied for the present problem. The soft photon flux per unit dimensionless photon energy ǫ s [s −1 cm −2 ] is written as dF s /dǫ s . We assume in this paper that soft photons are supplied by the thermal radiation from the neutron star. If a soft emission component cannot be represented by a black body spectrum, we use the observed spectrum of dF s /dǫ s corrected for the interstellar absorption. If it can be represented by a blackbody component, we adopt where r and r NS refers to the distance from the star center and the neutron star radius, respectively; the Planck function [ergs s −1 cm −2 ster −1 ergs −1 ] is related with the blackbody temperature (δ ≡ kT /m e c 2 ) as Substituting equations (23) into (22), and executing the integration over ǫ γ , we obtain (Appendix B) where µ * and ǫ s * are related with µ and ǫ s by the Lorentz transformation as Since η γ IC ±,i (l) is defined at Γ = b l , Γ is replaced with b l in equation (28). Moreover, A j represents the observed emission area of the ith blackbody component normalized by 4πr 2 NS . For example, if a soft photon field consists of the whole neutron star surface emission with temperature T s and a heated polar cap emission with area πr 2 h and temperature T h , we obtain A 1 = 1.0, δ 1 = kT s /m e c 2 and A 2 = πr 2 h /(4πr 2 NS ), To define a dimensionless pair-production redistribution function, we consider the case when a γ-ray photon with energy m e c 2 b i materializes as a pair with Lorentz factor b l . Dividing the pair production rate [s −1 ] per Lorentz factor by ω p , we obtain where the pair-production threshold energy is defined by The collision angles cos −1 µ ± are determined by the geometry, in the same manner as η γ IC ±,i . The differential cross section for pair production is given by where σ T refers to the Thomson cross section and the center-of-mass quantities are defined as where µ = µ + (or µ − ) for η p+,l (or η p−,l ). Let us now discretize the differential equations. Denoting f k l as the quantity f evaluated at ξ = ξ k and Γ = b l (or ǫ γ = b l ), the source terms (eq. [17] and [18]) can be represented as The dimensionless γ-ray distribution functions are defined above the energy m e c 2 b l0 = 5.11 × 10 7 eV as where i = l 0 + 1, . . . , −1, 0, 1, 2, . . . , l m . On the other hand, the particle distribution functions are defined above Γ ≥ b 1 = 10 5 as where l = 1, 2, . . . , l m . Note that n ± are normalized by B at each point ξ k , while g ± are by B(ξ in ).
Since the IC scatterings do not alter the particle number, the sum of the first and the second terms in the right-hand side of equation (34) over l, vanishes identically. Therefore, holds, because the created number of positrons is identical to that of electrons.

Current Conservation
Let us, for the moment, consider the current conservation problem. Since we are considering a pure curvature process (i.e., neglecting the pitch-angle dependence of n ± ), and since Γ ≫ 1 holds, all the positrons (or electrons) migrate with velocity +c (or −c). As a result, particles at ξ = ξ k−1 migrate to ξ = ξ k after time (ξ k − ξ k−1 )/c. Figure 1 shows this situation, describing the particle propagation along the characteristics, whose gradient is given by (see eq. 14) in the phase space (ξ,Γ). Adding the produced pairs between ξ k−1 and ξ k , which are depicted by open circles in figure 1, we obtain for electrons. Adding both sides of these two equations, and utilizing equation (37), we obtain On these grounds, we can conclude that the current density is virtually kept constant for k (i.e., for position ξ = ξ k ), provided Γ ≫ 1. We solve the advection-type differential equations (14) and (15) on the two-dimensional phase space (ξ, Γ) by the Cubic Interpolated Propagation (CIP) method (e.g., Yabe & Aoki 1991, Yabe, Xiao, & Utsumi 2001.

Gamma-ray Boltzmann Equations
In this subsection, we consider the γ-ray Boltzmann equations. For simplicity, we assume that the density of outwardly (or inwardly) propagating γ-rays decreases (or increases) at the same rate with the magnetic field strength; this assumption gives a good estimate when W ≪ ̟ LC holds. Under this assumption, the γ-ray Boltzmann equations become (e.g., Paper VII)

±cB
∂ ∂s where η CV is the curvature radiation rate [s −1 ] into the energy interval m e c 2 (ǫ γ + dǫ γ ) by a particle migrating with Lorentz factor Γ. The first term in the right-hand side expresses the γ-ray absorption rate by pair production, while the second term does the emission rate via IC scatterings and curvature radiation.
Integrating equations (43) over ǫ γ between b i−1 and b i , we obtain where η CV+,i (l) refers to the curvature radiation rate η CV (b i , b l ) by a positron having Lorentz factor Γ = b l . In the same manner, η CV−,i (l) does that by an electron. Evaluating the right-hand side at ξ = ξ k−1 , we obtain the advanced distribution functions by the Euler method as follows:

Boundary Conditions
As we have described, the set of Maxwell and Botzmann equations consist of equations (8), (9), (14), (15), and (45). The Poisson equation and the γ-ray Boltzmann equations are ordinary differential equations, which can be straightforwardly solved by a simple discretization. On the other hand, the hyperbolictype partial differential equations (14) and (15) are solved by the CIP method.
In this section, we consider the boundary conditions to solve the set of Maxwell and Boltzmann equations. At the inner (starward) boundary (ξ = ξ in ), we impose (Paper VII) where the dimensionless positronic injection current across the inner boundary ξ = ξ in is denoted as j in .
As we described at the end of § 3.2, current density is conserved along the magnetic flux tube if particles are accelerated well above the initial Lorentz factor. Thus, we obtain at ξ = ξ in At the outer boundary (ξ = ξ out ), we impose The current density created in the gap per unit flux tube can be expressed as We adopt j gap , j in , and j out as the free parameters.
We have totally 2l m + 2(l m − l 0 + 1) + 4 boundary conditions (46)-(52) for 2l m + 2(l m − l 0 + 1) + 2 unknown functions n ±l , g ±i , Ψ, and E . Thus two extra boundary conditions must be compensated by making the positions of the boundaries ξ in and ξ out be free. The two free boundaries appear because E = 0 is imposed at both the boundaries and because j gap is externally imposed. In other words, the gap boundaries (ξ in and ξ out ) shift, if j in and/or j out varies.
In practice, we give n + (ξ in , Γ) and n − (ξ in , Γ) by a continuous function of Γ near the lowest Γ bin, b 1 , to avoid numerical oscillation due to the step function represented by equations (48) and (52). Nevertheless, detailed functional forms of the injected particle spectrum little affect the results.

Gap Position vs. Particle Injection
Let us briefly examine how the gap position changes as a function of j in and j out . For a detailed argument, see § 2.4 in Paper VII. Defining the particle spatial number density bȳ we can rewrite the Poisson equation (9) as where the created charge density in the gap is defined as and the two-dimensional screening effect (i.e., the −ψ/∆ 2 ⊥ term) is neglected in equation (55). The created charge density increases outward due to discharge and satisfies n gap (ξ in ) = −j gap and n gap (ξ out ) = j gap . We have n gap < 0 and f null > 0 in the inner part of the gap, while we have n gap > 0 and f null < 0 in the outer part. It follows from equation (55) that the created particles partially screen the original electric field. If j gap = j in − j out + B z (ξ in )/B in holds, we obtain dE /dξ = 0 at the inner boundary; that is, the gap is quenched due to the overfilling of space charges. Then, how about the injected particles? Let us neglect the created particles by imposing j gap ≪ 1, and consider the following three cases: • j in = j out In this case, E maximizes (i.e., f null vanishes) at the point where B z = 0 holds. That is, the gap is located around the null surface.
• j in ∼ 1, j out ∼ 0 In this case E maximizes at the point where B z = −B holds. That is, the gap is located near to the light cylinder.
• j in ∼ 0, j out ∼ 1 In this case E maximizes at the point where B z = B holds. That is, the gap is located near to the polar cap surface.
In short, created particles in the gap can quench the gap, while injected particles consisting of a single charge only shift it. The gap center is located at the 'generalized null surface' on which f null vanishes, provided that the created particle density is negligible (i.e., if n gap ≪ 1). Strictly speaking, injected particles emit γ-rays that can materialize as pairs; therefore, injected particles contribute for screening of E in a sense. Nevertheless, E is screened by the created particles discharged by E itself, not by the originally injected ones that do not return.

Application to Individual Pulsars
To solve the set of the Maxwell and Boltzmann equations, we must specify the X-ray field, dF s /dǫ s , which is necessary to compute the pair-production redistribution function (eq. [30]). In this paper, we use the X-ray fluxes and spectra observed for individual rotation-powered pulsars to specify dF s /dǫ s . In § 4.1, we summarize the observed properties of the X-ray field. Then we apply the theory to the Vela pulsar in § 4.2, to PSR B1706-44 in § 4.3, to the Geminga pulsar in § 4.4, and to PSR B1055-52 in § 4.5. We assume that the solid angle of the emitted γ-rays is 1 ster throughout this paper.

Input Soft Photon field
We consider the photons emitted from the neutron star surface as the seed photons for (γ-γ) pairproduction and IC scatterings. That is, we do not consider power-law X-ray components, because they are probably magnetospheric and beamed away from the accelerator. We evaluate the IR photon field, which is needed to compute the IC scattering rate, from the Rayleigh-Jeans tail of the surface (black-body) component. In table 1, we present the observed properties of the four γ-ray pulsars exhibiting surface X-ray components, in order of spin-down luminosity, L spin .
Vela (J0835-4513) From Chandra observations in 0.25-8.0 keV, the spectrum of this pulsar is turned out to consist of two distinct component: A soft component and a hard, power-law component. The hard component is probably a magnetospheric origin and hence will be beamed away from the gap. Thus, we consider only the soft component as the Xray field illuminating the outer gap. This component can be modeled either as a blackbody spectrum with kT = 1.49MK and A s = 0.045A * (d/0.25) 2 , or as a magnetic hydrogen atmosphere spectrum with effective temperature kT = 0.68MK (Pavlov et al. 2001). Based on high-resolution Ca II and Na I absorptionline spectra toward 68 OB stars in the direction of the Vela supernova remnant, Cha, Sembach, and Danks (1999) determined the distance to be 250 ± 30 pc. B1706-44 (J1710-4432) Gotthelf, Halpern, Dodson (2002) reported a broad, single-peaked pulsed profile with pulsed fraction of 23%, using the High Resolution Camera on-board the Chandra X-ray observatory. They fitted the spectroscopic data to find (at least) two components, e.g., a blackbody of kT = 143 eV with A = 0.129A * and a power-law component with photon index of −2.0. We consider that the former component illuminates the gap efficiently and neglect the second one. We adopt d = 2.5 kpc as a compromise between the smaller dispersion-measure distance of 1.8 kpc based on the free electron model by Taylor and Cordes (1993) and the larger H I kinematic distance of 2.4-3.2 kpc derived by Koribalski et al. (1995). Geminga (J0633+1746) The X-ray spectrum consists of two components: the soft surface blackbody with kT s = 50 eV and A s = 0.22A * (d/0.16) 2 and a hard power law with α = −1.6 (Halpern & Wang 1997). A parallax distance of 160pc was estimated from HST observations (Caraveo et al. 1996). B1055-52 (J1059-5237) Combining ROSAT and ASCA data, Greiveldinger et al. (1996) reported that the X-ray spectrum consists of two components: a soft blackbody with kT s = 68 eV and A s = 7.3A * (d/1.53) 2 and a hard blackbody with kT h = 320 eV and A h = 2.3 × 10 −4 (d/1.53) 2 . The distance is estimated to be 1.53 kpc from dispersion measure (Taylor & Cordes 1993).

The Vela pulsar
We first apply the theory to the Vela pulsar, using the X-ray field modeled by the blackbody spectrum in § 4.2.1-4.2.3, and by the hydrogen atmosphere spectrum in § 4.2.4.

Acceleration Field and Characteristics
Let us first consider the spatial distribution of the acceleration field, which is solved from the Poisson equation (9). We may notice here that s is related with ξ by equation (7). For the Vela pulsar, a small value j gap = 4.6 × 10 −5 gives the best-fit spectrum (see § 4.2.3 for details).
To compare the effects of particle injection, we present the E distribution for the three cases of j in = 0 (solid), 0.25 (dashed), and 0.50 (dash-dotted) in figure 2. The magnetic inclination is chosen to be α i = 75 • . We adopt j out = 0 throughout this paper, unless its value is explicitly specified. The abscissa denotes the distance along the last-open field line (i.e., the s coordinate used in eqs. [1] and [2]) normalized by ̟ LC .
As the solid line shows, E is located around the null surface when there is no particle injection across either of the boundaries. Moreover, E varies quadratically, because the Goldreich-Julian charge density deviates from zero nearly linearly near to the null surface.
As the dashed and dash-dotted lines indicate, the gap shifts outwards as j in increases. When j in = 0.5 for instance, the gap is located on the half way be-tween the null surface and the light cylinder. This result is consistent with Papers VII, VIII, and IX.
In figure 3, we present the characteristics of partial differential equation (14) for positrons by solid lines, together with E (ξ) when j in = 0.25 and α i = 75 • (i.e., the dashed line in fig. 2). We also superpose the equilibrium Lorentz factor that would be obtained if we assumed the balance between the curvature radiation reaction and the electrostatic acceleration, as the dotted line. It follows that the particles are not saturated at the equilibrium Lorentz factor in most portions of the gap.
In the outer part of the gap where E is decreasing, characteristics begin to concentrate; as a result, the energy distribution of outwardly propagating particles forms a 'shock' in the Lorentz factor direction. However, the particle Lorentz factors do not match the equilibrium value (dotted line). For example, near the outer boundary, the particles have larger Lorentz factors compared with the equilibrium value, because the curvature cooling scale is longer than the gap width. Thus, we must discard the mono-energetic approximation that all the particles migrate at the equilibrium Lorentz factor as adopted in Papers I through IX. We instead have to solve the energy dependence of the particle distribution functions explicitly.
The particles emit γ-rays not only inside of the gap but also outside of it, being decelerated by the curvature radiation-reaction force. The length scale of the deceleration is given by Since the typical Lorentz factor is a few times of 10 7 , l curv is typically much less than ̟ LC . Therefore, the escaping particles lose their energy well inside of the light cylinder.

Particle Energy Distribution
As we have seen in the foregoing subsection, the distribution function of the outwardly propagating particles forms a 'shock' in the outer part of the gap. In figure 4, we present the energy distribution of particles at several representative points along the field line. At the inner boundary (s = 0.184̟ LC ), particles are injected with Lorentz factors typically less than 4 × 10 6 as indicated by the solid line. Particles migrate along the characteristics in the phase space and gradually form a 'shock' as the dashed line (at s = 0.205̟ LC ) indicates, and attains maximum Lorentz factor at s = 0.228̟ LC as the dash-dotted line indicates. Then they begin to be decelerated gradually and escape from the gap with large Lorentz factors ∼ 2.5 × 10 7 (dotted line) at the outer boundary, s = s out = 0.241̟ LC .
Even though the 'shock' is captured by only a few grid points for the dash-dotted line in figure 4, the CIP scheme accurately satisfies the conservation of the total current density, For this case, j tot is accurately conserved at 0.25 within 0.2% errors even at the 'shock'.

Formation of Power-law Gamma-ray Spectrum
So far, we have seen that the outwardly propagating particles are not saturated at the equilibrium value and that such particles escape from the gap with sufficient Lorentz factors suffering subsequent cooling via curvature process. It seems, therefore, reasonable to suppose that a significant fraction of the γ-ray luminosity is emitted from such escaping particles.
We present the γ-ray spectrum emitted from outwardly propagating particles (i.e., positrons) for the Vela pulsar with α i = 75 • , j gap = 4.6 × 10 −5 , j in = 0.25 in figure 5. The dashed line represents the γ-ray flux emitted within the gap, while the solid one includes that emitted outside of the gap by the escaping particles. Therefore, the difference between the solid and the dashed lines indicates the γ-ray flux k-1 k Fig. 1.-Schematic picture of characteristics in the phase space (ξ,Γ), where the dimensionless coordinate ξ is related with s by equation (7). Positrons (or electrons) move along these characteristics when E > 0 (or E < 0). The filled circles indicate the particles moved from ξ = ξ k−1 to ξ k in a certain time interval. Created particles (open circles) are assumed to migrate immediately outwards after the birth.   emitted by the particles migrating outside of the gap. For comparison, we plot the phase-averaged EGRET spectrum, which is approximated by a power law with a photon index −1.7 (Kanbach et al. 1994). It follows from the figure that the γ-ray spectrum in 100 MeV-GeV energies can be explained by the curvature radiation emitted by the escaping particles. We adjusted the transfield thickness as D ⊥ = 0.16̟ LC = 2.8W so that the observed flux may be explained. The luminosity of the γ-rays emitted outside of the gap contribute 48% of the total luminosity 5.08 × 10 33 ergs s −1 between 100 MeV and 20 GeV. In another word, we do not have to assume a powerlaw energy distribution for particles (as assumed in some of the previous outer-gap models) to explain the power-law γ-ray spectrum for the Vela pulsar. This conclusion is natural, because a power-law energy distribution of particles will not be achieved by an electrostatic acceleration, and because magnetohydrodynamic shocks (i.e., real shocks) will not be formed in the accelerator.
Because the X-ray field is dense for this young pulsar, the pair-production mean free path, and hence the gap width becomes small (for details, see Hirotani & Okamoto 1998, Hirotani 2000a, Paper IV; Hirotani 2001, Paper V). As a result, the potential drop in the gap 2.24 × 10 13 V is only 0.81 % of the electromotive force (EMF) exerted on the spinning neutron star surface ∼ µ/̟ 2 LC = 2.79 × 10 15 V. Nevertheless, this potential drop is enough to accelerate particles into high Lorentz factors, 10 7.3 .

Hydrogen Atmosphere Model
We next briefly examine the γ-ray spectrum when the soft photon field is modeled as the hydrogen atmosphere spectrum (Pavlov et al. 2001), rather than the blackbody spectrum. In figure 6, we present the resultant spectrum for the same parameter set as in figure 5, except for the fitting model for the observed X-ray field. We fix s in instead of j gap for comparison purpose. It follows that IC flux increases 3 times while j gap does only 2 times. This is because the hydrogen atmosphere model gives much low-energy (IR) photons compared with the blackbody model, as can be understood from its small effective temperature and large emitting area. Nevertheless, the obtained GeV spectrum little differ from the blackbody model. This is due to the 'negative feedback effect' of the gap electrodynamics (Paper V), which ensures the existence of solutions in a wide parameter range and suggests

Solutions in a Wide Parameter Space
With both blackbody and the hydrogen atmosphere models, we can explain the observed γ-ray spectrum in a wide parameter space 45 • < α i < 75 • and 0.125 < j in < 0.25, by appropriately choosing j gap and D ⊥ . With increasing j gap (≪ 1), the νF ν [Jy Hz] peak energy increases because of the increased W , while the sub-GeV spectrum becomes hard because of the significant γ-ray emission within the extended gap (rather than outside of it) above GeV energies. On the other hand, D ⊥ affects only the normalization of the γ-ray flux; the γ-ray luminosity is proportional to D 2 ⊥ . Let us first fix j in at 0.25 and consider how the best-fit values of j gap and D ⊥ depend on α i . As we have seen, they are j gap = 4.6 × 10 −5 and D ⊥ = 0.155̟ LC = 0.84s in for α i = 75 • , where the gap inner boundary is located at s = s in = 0.1835̟ LC . However, the ratio D ⊥ /s in increases with decreasing α i and becomes 1.06 for 45 • (cf. 0.84 for α i = 75 • ). From geometrical consideration, we conjecture that D ⊥ should not greatly exceed s in ; we thus consider that α i > 45 • is appropriate for j in = 0.25. A large α i is preferable to obtain a small D ⊥ /̟ LC . However, since the radio pulsation shows a single peak, we consider that α i is not too close to 90 • . On these grounds, we adopted α i = 75 • as a moderate value in this section.
Let us next fix α i at 75 • and consider how the bestfit values of j gap and D ⊥ depend on j in . The ratio D ⊥ /̟ LC increases with decreasing j in and becomes 0.94 for j in = 0.125 (c.f. 0.84 for j in = 0.25). This is because the decreased flux of the outwardly migrating particles (due to decreased j in ) must be compensated by a large D ⊥ to produce the same γ-ray flux. Thus, we consider j in > 0.125 is appropriate for α i = 75 • . A large j in is preferable to obtain a small D ⊥ /̟ LC . However, for j in > 0.25, the gap is so extended that a significant γ-rays are emitted above GeV within the gap; as a result, the sub-GeV spectrum becomes too hard. On these grounds, we adopted j in = 0.25 as a moderate value.
For 45 • < α i < 75 • , 0.125 < j in < 0.25, and appropriately chosen j gap and D ⊥ , TeV flux is always less than 3 × 10 10 JyHz. Thus, one general point becomes clear: TeV flux is unobservable with current ground-based telescopes, provided that the emission solid angle is as large as 1 ster and that the surface thermal (not magnetospheric) X-rays are upscattered inside and outside of the gap. Since the magnetospheric X-rays will be beamed away from the gap and their specific intensity is highly uncertain, we leave the problem of the upscatterings of magnetospheric (power-law) X-rays untouched.

PSR B1706-44
We next apply the theory to a Vela-type pulsar, PSR B1706-44. To consider D ⊥ /s in as small as possible, we adopt a large magnetic inclination 75 • . We compare νF ν spectra for the three cases: j in = 0.4, 0.2, and 0.1.
To examine how the sub-GeV spectrum depends on j in , we fix the νF ν peak at the observed value, ∼ 2 GeV, by adjusting j gap appropriately. For the solid, dashed, and dash-dotted lines in figure 7, we adopt j gap = 2.2 × 10 −4 , 1.8 × 10 −4 , and 1.5 × 10 −4 , respectively, and d = 2.5 kpc. Moreover, the perpendicular thickness is adjusted so that the predicted flux may match the observed value (3.2 × 10 13 JyHz) at 1.4 GeV; for the solid, dashed, and dash-dotted lines, they are D ⊥ = 1.05̟ LC = 3.5s in , 0.61̟ LC = 4.1s in , 0.49̟ LC = 5.3s in , respectively. The sub-GeV spectrum becomes harder with increasing j in , because the ratio of the flux emitted outside of the gap and that emitted within it decreases as the gap extends with increasing j in . However, as the solid line indicates, the obtained sub-GeV spectrum is still too soft to match the observation. For a smaller j in (< 0.1), the sub-GeV spectrum becomes further soft. For a nonzero j out , the inwardly shifted gap is shrunk to emit smaller γ-ray flux, which results in a further greater D ⊥ /s in . On these grounds, the predicted sub-GeV spectrum becomes too soft or the γ-ray flux becomes too small (i.e., D ⊥ becomes too large) for any parameter set of j in and j out , if α i = 75 • and d = 2.5 kpc.
For a small inclination (α i < 75 • ), the gap is located relatively outside of the magnetosphere, because the null surface crosses the last-open field line at large distances from the star. Because the pairproduction mean-free path becomes large far from the star, the gap is extended for a small magnetic inclination. In such an extended gap, particles saturate at the equilibrium Lorentz factor at the outer part (Takata et al. 2002) and emit most of the γ-rays around the central energy of curvature radiation. As a result, a hard sub-GeV spectrum can be expected; however, we have to assume a large D ⊥ that exceeds Nevertheless, if d is much less than 2.5 kpc and α i < 75 • , we can explain the observed spectrum with moderate D ⊥ . For example, for d = 1 kpc, α i = 45 • , j in = 0.4, and j gap = 3.8 × 10 −3 , the gap exists in 0.51̟ LC < s < 0.74̟ LC and the resultant spectrum (dotted line in figure 7) matches the observation relatively well with a reasonable thickness, On the other hand, for a large inclination (α i > 75 • ), the sub-GeV spectrum becomes softer than α i = 75 • case for the same j in , j out , and j gap . Therefore, the spectrum will not match the observation whatever distance we may assume.
We can alternatively consider a CHR-like outer gap. Assuming a small D ⊥ (say, 0.05), we find that the gap extends along the field lines due to the screening effect of the zero-potential walls. In most portions of this extended gap, particles are nearly saturated at the equilibrium Lorentz factor. As a result, the sub-GeV spectrum becomes hard and match the observation with appropriate peak energy around 2 GeV. However, in this case, the γ-ray flux becomes too small to match the observed value, unless we adopt an unrealistic distance (say, 300 pc).
In short, the phase-averaged EGRET spectrum for this pulsar cannot be explained by our current model for any combinations of α i , j in , j out , j gap , and D ⊥ , if d = 2.5 kpc. Therefore, we suggest a small distance (e.g., 1 kpc) for this pulsar.

The Geminga Pulsar
Let us apply the theory to a cooling neutron star, the Geminga pulsar. For a small α i (e.g., 45 • ), the gap is so extended that the outer boundary exceeds the light cylinder. For a larger α i , on the other hand, not only the outer-gap emission, but also a polar-cap one could be in our line of sight. Since there has been no radio pulsation confirmed, we consider a moderate magnetic inclination α i = 60 • .
Since the soft photon field is less dense compared with young pulsars like Vela or B1706-44, the gap is extended along the field lines. For the set of parameters j in = 0.25, j out = 0, and j gap = 8.0 × 10 −5 , which gives the best-fit γ-ray spectrum, we obtain W ∼ 0.59̟ LC . We present the spatial distribution of E obtained for this set of parameters as the dashed line in figure 8, as well as the characteristics (solid lines) and the equilibrium Lorentz factor (dotted line). In the outer part of this extended gap, ρ GJ gradually increases with the distance, s. As a result, E (s) deviates from quadratic distribution and decrease gradually in the outer part of the gap. Because of this extended feature, particles are nearly saturated at the equilibrium Lorentz factor. In another word, the mono-energetic approximation adopted in Papers I-IX is justified for this middleaged pulsar. Particle distribution function forms a strong 'shock', which is captured only with one grid point, in 0.4 < s/̟ LC < 0.45. As a result, j tot fluctuates a little as figure 9 indicates. Nevertheless, it returns to the 0.2525 level, which is 1 % greater than the value it should be (0.250), as the characteristics begin to be less concentrated beyond s = 0.45̟ LC .
In figure 10, we present the resultant γ-ray spectrum. Because of the nearly saturated motion of the particles, they lose most of their energy within the gap. As a result, γ-ray luminosity associated with the escaping particles (5.3 × 10 31 ergs s −1 ), is negligibly small compare to that emitted within the gap (1.11 × 10 33 ergs s −1 ), which is represented by the dashed line.
It follows from the figure that the observed, soft spectrum between 200 MeV and 6 GeV can be explained by the present one-dimensional model. It is interesting to compare this result with what obtained for the Vela pulsar (figs. 5 and 6). Between 100 MeV and 1 GeV energies, both spectra are formed by the superposition of the curvature radiation emitted by the particles having different energies at different positions. The important difference is that the particles are saturated at the equilibrium Lorentz factor in the gap for the Geminga pulsar, while they are nearly mono-energetic but only decelerated via curvature process outside of the gap for the Vela pulsar. Because the particles are no longer accelerated outside of the gap, they emit γ-rays in lower energies compared with those still being accelerated in the gap. As a result, the γ-ray spectrum for the Vela pulsar becomes softer than that for the Geminga pulsar. Extending this consideration, we can predict that a γ-ray spectrum below GeV is soft for a young pulsar and tends to become hard as the pulsar ages.

PSR B1055-52
Let us finally apply the present theory to another middle-aged pulsar, B1055-52. To obtain a large γray flux for an appropriately chosen set of j in , j gap , and D ⊥ (< s in ), we adopt a large magnetic inclination,   -Total current density, j tot , for the same case as in figure 8. Even though the particle distribution function forms a strong 'shock' in the Lorentz factor direction in 0.4̟ LC < s < 0.45̟ LC (see fig. 8), j tot is conserved relatively accurately.
Since the acceleration field and the particle energy distributions are similar to the Geminga pulsar, we present only the computed γ-ray spectra for this pulsar in figure 11. The solid and dashed lines represent the spectra for the distances 0.50 and 1.53 kpc, respectively. For d = 0.50 kpc, j gap = 4.5 × 10 −3 is chosen so that the peak energy of curvature radiation may match the observed peak energy. In this case, the gap exists in 0.1450̟ LC < s < 0.3882̟ LC . For d = 1.53 kpc, j gap = 4.0 × 10 −2 is chosen and the gap exists in 0.1475̟ LC < s < 0.3818̟ LC .
It follows from the figure that the solid line matches the observed flux with a reasonable transfield thickness, D ⊥ = 0.61s in = 0.38W . On the other hand, for the dashed line, we have to choose D ⊥ = 0.25̟ LC = 1.7s in , which is unrealistic. The observed fluxes cannot be explained with acceptable gap width (e.g., D ⊥ < s in ) no matter what we may adjust j in , j out , and j gap if d = 1.53 kpc.
On these grounds, we conjecture that the distance 1.5 ± 0.4 kpc determined from the dispersion measure (Taylor & Cordes 1993) is too large and that a more closer distance, such as 500 pc derived from ROSAT data analysis (Ögelman & Finley 1993) or 700 pc estimated from a study of the extended nonthermal radio source around the pulsar (Combi, Romero, Azcárate 1997), is plausible. This conclusion is consistent with the results obtained in Paper IX under the monoenergetic approximation that all the particles are saturated at the equilibrium Lorentz factor, which is justified for middle-aged pulsars.

Discussion
In summary, we have quantitatively examined the stationary pair-production cascade in an outer magnetosphere, by solving the set of Maxwell and Boltzmann equations one-dimensionally along the magnetic field. We revealed that an accelerator (or a potential gap) is quenched by the created pairs in the gap but is not quenched by the injected particles from outside of the gap, and that the gap position shifts as a function of the injected particle fluxes: If the injection rate across the inner (or outer) boundary approaches the typical Goldreich-Julian value, the gap is located near to the light cylinder (or the star surface). It should be emphasized that the particle energy distribution is not represented by a power law, as assumed in some of previous outer-gap models. The particles escape from the gap with sufficient Lorentz factors and emit significant photons in 100 MeV-3 GeV energies via curvature radiation outside of the gap. The γ-ray spectrum including this component explains the phase-averaged EGRET spectra for the Vela pulsar, the Geminga pulsar and PSR B1055-52 between 100 MeV and 6 GeV. TeV fluxes are unobservable with current ground-based telescopes for these pulsars.
We show that synchro-curvature process can be approximated by the pure-curvature one in the next subsection. We then point out an implication to the γ-ray luminosity versus the spin-down luminosity in § 5.2, and discuss future extensions of the present method in § § 5.3-5.5.

Synchrotron vs. Curvature Processes
Let us first examine the pitch angles of particles injected across the inner boundary. Imposing that the synchrotron cooling length scale, is greater than the distance, ε̟ LC , for the particles to migrate before arriving the outer gap, we obtain sin χ < 7.2 × 10 −6 Ω 2 ε 0.1 Γ 2 where ε 0.1 ≡ ε/0.1, Γ 2 ≡ Γ/10 2 , Ω 2 ≡ Ω/(10 2 rad s −1 ), and B 10 ≡ B/(10 10 G). Particles having initial pitch angles greater than this value will lose transverse momentum to reduce the pitch angles below this value, while migrating the distance ε̟ LC . Therefore, if the particles are supplied from the polar cap with Γ ∼ 10 2 for instance, we can safely state that the the pitch angles are less than 10 −5 rad. If a particle having such a small pitch angle is accelerated to Γ = 10 7 in a weak magnetic field region B = 10 6 G, the ratio between the synchrotron and the curvature radiation rate becomes where B 6 ≡ B/10 6 G, Γ 7 ≡ Γ/10 7 . It follows that the synchro-curvature radiation can be approximated by a pure-curvature one for the particles injected across the inner boundary. It is noteworthy that the density of created particles is much less than that of injected ones (i.e., j gap ≪ j in ); therefore, the synchrocurvature effects for the created particles, which have much greater pitch angles compared with the injected ones, are negligibly small for the gap electrodynamics as well for the resultant γ-ray spectrum.

Gamma-ray vs. Spin-down Luminosities
It should be noted that the emission from the escaping particles attain typically 40% of the total γray luminosity for young pulsars. Thus, it is worth mentioning its relationship with the spin-down luminosity, where the braking index n is related to the spin-down rate asΩ = −kΩ n .
If the spin down is due to the magnetic dipole radiation, we obtain n = 3. The outwardly propagating particles escape from the gap with spatial number density (eq. [11]) where B out = B(ξ out ). Therefore, the energy carried by the escaping particles per unit time is given by where Γ esc (∼ 10 7.5 ) refers to the Lorentz factor of escaping particles. Note that Γ esc is essentially determined by the equilibrium Lorentz factor (dotted line in fig. 3) near the gap center. Since the equilibrium Lorentz factor depends on the one-fourth power of E , the variation of Γ esc on pulsar parameters is small. We can approximate B out as where r out refers to the distance of the outer boundary of the gap from the star center. Let us assume that the position of the gap with respect to the light cylinder radius, r out /̟ LC , does not change as the pulsar evolves; this situation can be realized if j in − j out is unchanged. Evaluating B at r = 0.5̟ LC , we obtain where n = 3 is assumed in the second line. To derive this conclusion, it is essential that the particles are not saturated at the equilibrium Lorentz factor. Thus, the same discussion can be applied irrespective of the gap position or the detailed physical processes involved. For example, an analogous conclusion was derived for a polar-cap model by Harding, Muslimov, and Zhang (2002). It is, therefore, concluded that the observed relationship L γ ∝ L spin 0.5 merely reflects the fact that the particles are unsaturated in the gap and does not discriminate the gap position.
Let us compare this result with what would be expected in the CHR picture. Since the gap is extended significantly along the field lines in the CHR picture, particles are saturated at the equilibrium Lorentz factor to lose most of their energies within the gap, rather than after escaping from it. We can therefore estimate the γ-ray luminosity as where D φ refers to the azimuthal thickness of the gap, P CV [ergs s −1 (particle) −1 ] represents the curvature radiation rate. Noting that the particle motion saturates at the equilibrium Lorentz factor satisfying P CV /c = e(−dΨ/ds), recalling that the acceleration field is given by −dΨ/ds ≈ ΩBD 2 ⊥ /4ρ c c in the CHR picture, and evaluating B at r = ̟ LC , we obtain where n = 3 is assumed again in the second line. Even though the escaping particles little contribute to the γ-ray luminosity in the CHR picture, it is worth mentioning the work done by Crusius-Wätzel and Lesch (2002), who accurately pointed out the importance of the escaping particles in the CHR picture, when we interpret L γ ∝ L 0.5 spin relation. As we have seen, the particles being no longer accelerated contribute for the γ-ray luminosity that is proportional to L 0.5 spin . Reminding that the particles migrate with larger Lorentz factors than the equilibrium value in the outer part of the gap (see fig. 3), we can expect roughly half of the γ-ray luminosity is proportional to L 0.5 spin (mainly between 100 MeV and 1 GeV), and the rest of the half to L spin (mainly above 1 GeV). As a pulsar ages, its declined surface emission results in a large pair-production mean free path, and hence W . Because |ρ GJ | ∝ r −3 becomes small in the outer part of such an extended gap, E (ξ) deviates from quadratic distribution to decline gradually in the outer part ( fig. 8). As a result, particles tend to be saturated at the equilibrium value. On these grounds, we can predict that the γ-ray luminosity tends to be proportional to L spin with age, deviating from L 0.5 spin dependence for young pulsars.
In the present paper, we have examined the set of Maxwell and Boltzmann equations one-dimensionally both in the configuration and the momentum spaces (i.e., only ξ and Γ dependences are considered.) In the next three sections, we discuss the extension of the present method into higher dimensions

Returning Particles
If we consider the pitch-angle dependence of particle distribution functions, we can compute the radiation spectrum with synchro-curvature formula (Cheng and Zhang 1996). Moreover, we can also consider the returning motion of particles inside and outside of the gap. The returning motion becomes particularly important when both signs of charge are injected across the boundary. For example, not only positrons but also electrons could be injected across the inner boundary from the polar-cap accelerator. If E > 0 for instance, the injected electrons return in the gap. This returning motion significantly affects the Poisson equation, if their injection rate is a good fraction of the Goldreich-Julian value.
It remains an unsettled issue whether an outer-gap accelerator resides on the field lines on which a polarcap accelerator exists. To begin with, let us consider the case when the plasma flowing between the polar cap and the outer-gap accelerator is completely charge separated. Such a situation can be realized, for instance, if only positively charged particles are ejected outwardly from the polar cap while there is virtually no electrons ejected inwardly from the outer gap. Neglecting the pair production, current conservation law gives the charge density, ρ e , per unit magnetic flux tube as where v refers to the particle velocity along the field line, and j tot the conserved current density per magnetic flux tube. At each point along the field line, ρ e should match ρ GJ . If the field line intersects the null surface, ρ e must vanish there; this obviously violates the causality in special relativity. Therefore, a stationary ejection of a completely charge-separated plasma from the polar cap can be realized only along the field lines between the magnetic axis and those intersecting the null surface at the light cylinder. On these grounds, it was argued that an outer-gap accelerator, which is formed close to the last-open field line, may not resides on the same field lines on which a polar-cap accelerator resides. This has been, in fact, the basic idea that an outer gap will not be quenched, because the particles ejected from the polar cap will flow along the different field lines. This idea was welcomed in outer-gap models, because a gap has been considered to be quenched if the external particle injection rate becomes comparable to the Goldreich-Julian value, which was proved to be incorrect in this paper.
In general, however, the plasmas are not completely charge separated and consist of both signs of charge (e.g., positrons and electrons). Such a situation can be realized, for instance, if both charges are ejected outwardly from a polar-cap accelerator, or if positively charged particles are ejected outwardly from the polar cap while electrons are ejected inwardly from the outer gap, or if there is a pair production between the two accelerators. In these cases, the velocities of both charges will be adjusted so that both the current conservation and ρ e = ρ GJ are satisfied at each point along the field lines. Therefore, it seems likely that a polar-cap accelerator and an outer-gap accelerator reside on the same field lines.
To examine if there is a stationary plasma flow between the polar cap and the outer gap, we must extend the present analysis into two dimensional momentum space in the sense that the pitch-angle dependence of the particle distribution functions is taken into account in addition to the Lorentz factor dependence. For example, if both charges are ejected from the polar-cap accelerator, electrons will return in the outer gap, screening the original acceleration field in the gap, and violating the original balance of ρ e = ρ GJ outside of the gap. Because the returning motion of particles can be treated correctly if we consider the pitch-angle evolution of the distribution functions, and because the pair production is already taken into account, our present method is ideally suited to investigate the plasma flows and E distribution self-consistently inside and outside of the gap.

Unification of Outer-gap Models
In addition to the extension into a higher dimensional momentum space, it is also important to extend the present method into a two-or three-dimensional configuration space. In particular, determination of the perpendicular thickness, D ⊥ , is important to constrain gap activities. There have been, in fact, some attempts to constrain D ⊥ in the CHR picture. Since −dΨ/ds is proportional to BD 2 ⊥ , particles energies, and hence the γ-ray energies increase with increasing D ⊥ (for a fixed B). Zhang and Cheng (1997) constrained D ⊥ , by considering the condition that the γrays cause photon-photon pair production in the gap. Subsequently, Cheng, Ruderman, and Zhang (2000) extended this idea into three-dimensional magnetosphere and discussed phase-resolved γ-ray spectra for the Crab pulsar. In addition, Romani (1996) discussed the evolution of the γ-ray emission efficiency and computed the phase-resolved spectra for the Vela pulsar, by assuming that BD 2 ⊥ declines as r −1 . However, in these works, screening effects due to pair production has not been considered; thus, the obtained D ⊥ , as well as the assumed gap position along the magnetic field, are still uncertain.
On the other hand, in our approach (picture), D ⊥ is not solved but only adjusted so that the γ-ray flux may match the observations. Therefore, the question we must consider next is to solve such geometrical and electrodynamical discrepancies between these two pictures. We can investigate this issue by extending the present method into higher spatial dimensions.

Unification of Outer-gap and Polar-cap Models
Electrodynamically speaking, the essential difference between outer-gap and polar-cap accelerators is the value of the optical depth for pair production. In an outer-gap accelerator, pair production takes place via γ-γ collisions and its mean-free path is much greater than the light cylinder radius. Therefore, a pair production cascade takes place gradually in the gap. In such a gap, E is screened out by the 'generalized Goldreich-Julian charge density', f null (eq. [57]), which increases outwards if Ω · B > 0. Since f null is negative (or positive) in the inner (or outer) part of the gap, there is a surface on which the righthand side of equation (55) vanishes, as long as the injected current density is less than the Goldreich-Julian value. The gap is located around this 'gener-alized null surface', which is explicitly defined in § 2 of Paper VII.
On the contrary, in a polar-cap accelerator, pair production takes place mainly via magnetic pair production, of which mean free path is much less than the star radius for a typical magnetic field strength (B ∼ 10 12 G, say). As a result, a pair production avalanche takes place in a limited region, which is called as the 'pair formation front', in the gap (Fawley, Arons, & Sharlemann 1977;Harding & Muslimov 1998, 2001Shibata, Miyazaki, Takahara 1998Harding, Muslimov, Zhang 2002). In the pair formation front, a small portion of the particles return to screen out E . Such a returning motion can be self-consistently solved together with E by our present method, if we implement the magnetic pair production and the resonant IC scattering redistribution functions in the source terms of the particles' and γ-rays' Boltzmann equations. We can execute the same advection-phase computation in CIP scheme; thus, all we have to do is to add these source terms in the non-advection-phase computation, which is not very difficult. Since analogous boundary conditions (e.g., E = 0 for a space-charge limited flow) will be applied, we expect the present method is also applicable to a polar-cap accelerator. This is an issue to be examined in our subsequent papers.
One of the authors (K. H.) wishes to express his gratitude to Drs. K. Shibata and A. Figueroa-Vinās for valuable advice on numerical analysis, and to Drs. K. S. Cheng and C. Thompson for fruitful discussion on theoretical aspects. He also thanks Canadian Institute for Theoretical Astrophysics for welcoming him as a visiting researcher.
On these grounds, we neglect z and φ dependence of the distribution functions and approximate equation (A7) as Neglecting the χ ⊥ derivative in equation (A5), and assuming |ṽ p | = c cos χ, we finally obtain c cos χB p ∂ ∂s When the particle is moving with Lorentz factor Γ = 1/ 1 − β 2 , we can relate the particle rest frame quantities (with asterisks) with the observer's frame (without asterisks) as dΩ * s =