{'fraction_non_alphanumeric': 0.031241052932953838, 'fraction_numerical': 0.011485928895026695, 'mean_word_length': 4.44652357365621, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Pulsar Polar Cap Heating and Surface Thermal X-Ray Emission I. Curvature Radiation Pair Fronts

We investigate the effect of pulsar polar cap (PC) heating produced by positrons returning from the upper pair formation front. Our calculations are based on a self-consistent treatment of the pair dynamics and the effect of electric field screening by the returning positrons. We calculate the resultant X-ray luminosities, and discuss the dependence of the PC heating efficiencies on pulsar parameters, such as characteristic spin-down age, spin period, and surface magnetic field strength. In this study we concentrate on the regime where the pairs are produced in a magnetic field by curvature photons emitted by accelerating electrons. Our theoretical results are not in conflict with the available observational X-ray data and suggest that the effect of PC heating should significantly contribute to the thermal X-ray fluxes from middle-aged and old pulsars. The implications for current and future X-ray observations of pulsars are briefly outlined.


INTRODUCTION
X-ray emission has been detected from several dozen pulsars in observations by ROSAT (Becker & Trumper 1997), ASCA (Saito 1998) and Chandra (Zavlin et al. 2000). In cases where good spectral measurements are available, the emission in the X-ray band seems to have both thermal and non-thermal components. While the non-thermal components are most likely of magnetospheric origin, especially when they connect smoothly to higher energy spectra, the origin of the thermal components is less clear because there are several mechanisms for production of thermal emission in the soft X-ray band. Pulsars younger than about 10 5 yr are expected to have significant neutron star (NS) cooling components and observed X-rays from these pulsars may be consistent with standard cooling scenarios (Page 1998). However, many older pulsars, including millisecond (ms) pulsars, beyond the age where cooling emission drops sharply, have strong thermal X-ray emission. In this case, some heating mechanism(s) must be operating. Two which have been proposed are heating by particles, accelerated in the magnetosphere, flowing back to the NS PCs (PC heating: Ruderman & Sutherland 1975, hereafter RS75;Arons 1981, hereafter A81) and internal heating by frictional forces in the crust (Shibazaki & Lamb 1989).
Because PC heating emission is intimately tied to the pulsar acceleration mechanism, X-ray emission of older pulsars provides very strong constraints on pulsar models. Models which predict X-ray luminosities from PC heating which are in excess of observed luminosities are not viable. The earliest PC models, based on vacuum gaps at the poles (RS75), predicted PC heating luminosities that were way above the first Einstein detections and upper limits on pulsars (Helfand et al. 1980). However, predictions of the space-charge limited flow (SCLF) model of A81 were within all observed values. Similarly, the PC heating predictions of the first outer-gap models (Cheng, Ho & Ruderman 1986), where particles are accelerated in the outer magnetosphere, were well above the observed X-ray luminosities. Later models allowed for a reduction of this huge PC X-ray emission by scattering in a "pair blanket" which forms near the surface (Wang et al. 1998). It is clear that models with PC vacuum gaps and outer gaps, where half of all the particle acceleration energy returns to the surface to heat the PCs, will have severe problems preventing or hiding the excess X-ray emission. In this respect, SCLF models have a distinct advantage in that the number of particles that are accelerated downward to heat the PC are only a small fraction of the number of primary particles that are accelerated upward. This is because the accelerating electric field (parallel to the magnetic field) arises not from a pure corotation charge, but from a slight imbalance between the corotation charge and the actual charge. The field grows with height above the surface until electron-positron pairs from the accelerated, upward moving electrons screen the field above a pair formation front (PFF). A81, in his study of SCLF acceleration, found that the screening is primarily due to positrons which turn around and accelerate backward to the stellar surface. The number of returning ("trapped") positrons required to screen the E is very small compared to the number of primary particles (electrons). In his studies, he assumed that pairs were produced only by curvature radiation (CR) photons.
We have begun an updated study of PC SCLF acceleration (Harding & Muslimov 1998; hereafter HM98) based on the work of Muslimov & Tsygan (1992;hereafter MT92), who found that including effects of general relativity produced a substantial difference in the parallel electric field induced above the surface of a NS. In particular, the effect of inertial frame dragging induces a field that is typically a factor of 50 times higher than that derived by A81, which is that induced by flaring of polar field lines alone. This relativistic effect of inertial frame dragging and the consideration of radiation from inverse Compton scattering (ICS) of thermal X-rays from the NS surface by the primary electrons introduces profound changes in SCLF models. We (HM98) found that pairs are produced by ICS photons from accelerating electrons at lower altitudes than are pairs from CR photons. Pairs from ICS radiation could therefore screen the E before any CR photons are produced. We further found that positrons accelerating downward from ICS PFFs could potentially cause unstable acceleration if ICS pairs from positrons screened the E above the surface. In our previous study, we computed self-consistently the height at which screening begins by assuming that the CR and ICS PFFs are located at the point where the first pair is produced. We also have given rough estimates for the returning positron flux and polar cap heating luminosity expected from CR PFFs (Zhang & Harding 2000). However, we did not model the structure of the PFF and thus could not determine the flux of returning positrons as a fraction of the flux of primaries (the returning positron fraction) or even whether full screening occurs. Many unanswered questions therefore remain, such as whether ICS photons can always screen the E and whether there are enough returning positrons to produce pairs which screen E near the surface.
To address these questions, as well as provide PC heating rates for the new SCLF models, we have made a detailed and self-consistent study of the screening at both CR and ICS PFFs. We have found that while ICS photons are able to produce PFFs in nearly all pulsars, ICS radiation is able to completely screen the E only in pulsars with higher surface magnetic fields. The magnetic field value required for ICS screening is strongly dependent on the surface temperature of the NS. For PC temperatures T < 10 6 K, pulsars with surface B 0 > ∼ 0.1 B cr , where B cr ≡ 4.4 × 10 13 G is the critical field strength, are capable of ICS screening, but for T > ∼ 10 6 K, the required field decreases with temperature. If the ICS photons are unable to completely screen the E , then the primary electrons will keep accelerating until they radiate CR photons which can produce a PFF. In some cases, such as ms pulsars, incomplete ICS screening via non-resonant scattering can make a significant contribution to PC heating, especially for those pulsars where CR photons cannot produce any pairs. The results of these studies will be presented in two papers. This first paper will discuss CR screening and PC heating in lower field pulsars, i.e. those below the threshold for ICS screening. In a second paper, we will discuss ICS screening in higher field pulsars, as well as partial ICS screening in older and ms pulsars, and will address the issues of lower PFFs and stability of ICS PFFs. In the present paper, we first discuss the SCLF electric field solution used and the conditions for the screening of this field. We then present analytic estimates and numerical calculations of the returning positron fraction and PC heating luminosity for CR screening and compare our predictions with ROSAT observations. The most recent relevant study is that of Hibschman & Arons (2001), who focus on the determination of theoretical radio pulsar death lines. Even though they use the frame-dragging electric field in their treatment of primary electron acceleration and incorporate both CR and ICS photons in the pair production, our approaches are essentially different. For example, in their calculation of the PC heating they use a very rough estimate of the returning positron fraction (which is generally inaccurate and is the same as that given by MT92). Also, they refer to the PFF as the location where the electric field is fully screened, whereas in this paper we define the PFF as the location where pair production begins. In the case of CR, this distinction is minor since the screening scale is small compared to the PFF height, but in the case of ICS the distinction is very important since the screening scale is comparable to the PFF height. Thus, the quantitative results of our studies are not directly comparable, although some of their conclusions are qualitatively similar. A more detailed comparison of our results with theirs will be given in our next paper, where we will present our theoretical radio pulsar deathlines.
In Section 2, we describe the model we employ to calculate the accelerating electric field in the PC region of pulsars. We outline our main assumptions and approximations. The electrodynamic boundary conditions at the stellar surface and upper PFF are specified. In Section 3, we present our calculations of the returning positron fraction. In this paper our analysis is limited to the parameter space where the CR is the primary source of photons producing pairs in the magnetic field. We give our analytic estimates first. Then, we discuss our numerical treatment of electron-positron pair dynamics and effect of partial backflow of positrons. In Section 4, we calculate both analytically and numerically the thermal X-ray luminosities from pulsar PCs produced by precipitating positrons. Summary and conclusions are given in Section 5.

Electric fields in the regime of space-charge limitation
In this Section we shall summarize some of the general solutions pertaining to the electrodynamics of open field line regions of rotating NSs derived earlier (MT92, MH97, and HM98) and which will be exploited throughout this paper. We should note that these solutions are based on the assumption that a regime known as space-charge limitation of current occurs in the acceleration of electrons ejected from the PCs of a NS. The effect of space-charge limitation is well-known since the very first experimental and theoretical studies of the vacuum diodes, and detailed discussion of this effect can be found in any textbook on electronic devices. In Section 2.2 we will briefly outline some of the features of the space-charge limitation specific to the physical conditions of the pulsar PC regions.
The general relativistic Poisson's equation describing the electric potential distribution (in the corotating frame of reference) within the region of open field lines in the magnetosphere of a rotating NS reads (MT92) where ρ is the space charge density, ρ GJ is the general relativistic analog of the Goldreich-Julian (GJ) space charge density, α = 1 − ǫ/η is the redshift function, ǫ = r g /R, r g is the gravitational radius of the NS, R is the stellar radius, η = r/R is the dimensionless radial coordinate, κ = ǫI/M R 2 , I and M are the moment of inertia and mass of the NS, respectively. The differential operators such as gradient, divergence, etc. should be taken in corresponding curvilinear coordinates. Eqn. (1) takes into account both the effect of dragging of inertial frames of reference and gravitational redshift. The former dramatically affects the electrodynamics of a NS, while the latter rather moderately modifies the corresponding electrodynamic quantities.
The study of particle acceleration within the region of open magnetic field lines in pulsars is grossly facilitated by the fact that for all known pulsars (including ms pulsars) the angular size of the PC is less than ∼ 0.2/ P/3 ms radian (P is a pulsar spin period). Thus, we can use a so-called small-angle approximation, which proves to be very satisfactory for most relevant problems. In this approximation, for general relativistic dipolar magnetic field, where Ω is the angular velocity of the NS rotation, B 0 is the surface value of the magnetic field strength at the magnetic pole and χ is the pulsar inclination angle. Here we use the dimensionless radial coordinate η (≡ r/R) and characteristics describing field-streamlines ξ = θ/θ(η) [which is a magnetic colatitude scaled by the half-opening angle of the polar magnetic flux tube, In this paper we continue our study of a regime of space-charge limitation of current, which allows us to derive an explicit expression for a steady-state distribution of space charge density (see MH97 for details) Using expressions (2)-(4) and appropriate boundary conditions we can solve Eqn. (1). In this paper we will be mostly interested in the solutions for the accelerating electrostatic potential (field) valid for altitudes less than one stellar radius from the PC surface i.e. for z = η − 1 < 1, even though their generalization for z > 1 proves trivial (HM98). Thus, for z < 1 the general solution for the potential can be expressed as where and Here Φ 0 ≡ (ΩR/c)B 0 R is the dimensional factor which is nothing but an order-of-magnitude maximum value of electrostatic potential generated by a magnetized globe (with dipolar external magnetic field) rotating in vacuo (see Deutsch 1955); k i andk i are the positive zeros of the Bessel functions J 0 and J 1 , respectively, indexed in ascending order; δ(η) is another correction factor [see Eqn. (A13)] accounting for the gravitational redshift effect, with H(1)δ(1) ≈ 1 (see e.g. HM98); and θ 0 is defined right after Eqn.
(3). The explicit expressions for the function F i are determined by the solution of Eqn.
(1) and appropriate electrodynamic boundary conditions. The solutions for the function F i relevant to our study will be presented in § 2.3 and the Appendix. In this paper we do not need and therefore do not present the explicit general-relativistic expressions for the magnetic field, but we shall refer an interested reader to our previous publications (MH97, HM98). Finally, our results will be based on the electrodynamic solutions of Eqn.
(1) subject to the standard Dirichlet boundary conditions.

Modification of space-charge limitation by returning positrons
In this Section we shall discuss the space-charge limitation of current taking into account the returning positron flux from the upper PFF. The idea of space-charge limitation of flow was in fact introduced into pulsar physics by Sturrock (1971): "What happens to the positrons produced by the electron-positron avalanche? If they were all returned to the surface, as one might expect by consideration of the electric field, the resulting space charge would reverse the sign of the electric field at the surface, cutting off the flux of primary electrons". Later on this idea was investigated by Michel (1974Michel ( , 1975, Tademaru (1973), and Cheng and Ruderman (1977), and then substantially quantified in the calculation of the electric fields produced by the relativistic electron beam within the region of open field lines in a NS magnetosphere (Arons & Scharlemann, 1978;Arons 1983, hereafter A83). They found that only a small fraction of positrons (e.g. ∼ θ 2 0 of the negative charge density already present in the electron beam, as estimated by Arons & Scharlemann) are returned relative to the primary electron beam. Note that the calculation of the fraction of positrons returning from the upper PFF is far from being a trivial problem. In this paper we attempt to address this problem again, taking into account most recent updates of the physics involved and revising some of the underlying assumptions.
As a result of injection of backflowing positrons into the initial electron beam, an excessive negative space charge unavoidably builds up which lowers the electrostatic potential above the PC surface. The lowering of potential above the surface and forcing its gradient to reverse sign at the very surface would inhibit electron emission, thus resetting the electric field to zero at the surface. The fundamental difference between the laboratory vacuum diode and PC voltage generator is that in the former the thermionic cathode is the only emitter of electrons, whereas in the latter the PC surface, serving as a cathode, both emits electrons and collects positrons from the upper PFF. Thus, for a vacuum diode the space-charge limitation of current is solely determined by the cathode, while for a NS PC this regime of operation depends on the upper PFF as well. Because of the intrinsic limitation of the total current by the GJ value, the injection of counterstreaming positrons suppresses the primary electron emission: the primary electrons should be ejected from the stellar surface with a density somewhat less than the GJ charge density to keep E = 0 at the surface. This effect provides a negative feedback between the primary electron ejection and flux of returning positrons, since the primary electrons ultimately determine (through the electrostatic potential drop they generate above the PC surface) the fraction of returning positrons. Now, using very basic arguments, we would like to demonstrate how the accelerating electric field is related with the fraction of returning positrons. In our previous paper (see HM98) we quantitatively explored the effect of rescaling of the lower boundary (e.g. stellar surface) caused by the formation of a lower pair front in the regime of space-charge limitation. The returning positrons add to the space charge density of primary electrons and lift up the boundary at which the ρ = ρ GJ condition is satisfied, similar to the effect produced by the formation of lower pair front. Let us use this reasoning and express the space charge density of a primary beam in terms of effective (rescaled) stellar radius, R E . For simplicity, we illustrate this by employing solution (4) at χ ≈ 0. The expression for the space charge density at the rescaled stellar surface then reduces to where η * = R/R E (≤ 1). If there were no returning positrons (η * ≡ 1) we would simply have The returning positrons perturb the primary electron beam causing extraction of less electrons than in the case of unperturbed electron beam and effectively mimic electrons missing in a beam.
Assuming that the regime of space-charge limitation is maintained, the following condition should be satisfied at the surface where δρ + is a small perturbation of (and of the same sign as) the electron space charge density caused by returning positrons. Introducing the dimensionless fraction of returning positrons, x + , so that δρ + = −σκx + , from the above equation we get that η * = (1 − x + ) 1/3 . Thus the positrons returning from the PFF reduce the accelerating electrostatic potential (field), so that the latter gets factorized by (see also HM98) It is now clear that the net effect of the injection (e.g. from the upper PFF) of counterstreaming positrons into primary electron beam is a weakening of the source of effective charge, or, equivalently, the transforming of general-relativistic parameter κ into (1 − x + )κ. Note that x + ≪ 1 in most cases. Thus the final result does not depend on the particular value of rescaled stellar radius and illustrates the dynamics of the feedback between accelerating electric field and fraction of returning positrons. Namely, the returning positrons tend to reduce the accelerating electric field and therefore suppress the ejection of primary electrons. The starving electron ejection is compensated by returning positrons, which mimic the outflowing electrons, thus recovering the original potential drop. Obviously, the increase in the fraction of backflowing positrons reduces the efficiency of pair production by primary electrons which may eventually cut off the supply of positrons themselves. One can expect some demand-supply balance to be established in the beam between the densities of primary electrons and backflowing positrons. Thus, we are now in a position to calculate acceleration of primary electrons and returning positrons produced at the PFF in a self-consistent way, taking into account the weakening of the primary beam by returning positrons. It is important that, after the appropriate parameters that favour the steady-state regime are found, we can use essentially the same expressions for the accelerating electric field/potential as in the case of a pure electron beam, simply because in a steady-state electrodynamic solution the returning positrons are undistinguishable from outflowing electrons. Then, the problem reduces to finding (through numerical iterations) a stable solution (see § 3.2.2 for details).

Screening condition at the upper PFF and accelerating electric field
In this paper we derive an appropriate solution for the accelerating electric field, and use this solution for a self-consistent analysis of a steady-state flow of charged particles (electrons and positrons) from the PC surface up through the PFF. We assume that both the PC surface and surface formed by the last open magnetic field lines are equipotential (Φ = 0). We must note that, within the context of the lower boundary condition, the PC surface is actually either the upper boundary of the stellar atmosphere or the bare stellar crust. Since the free emission of charges from the PC seems to be more favourable than charge starvation above the PC, the electron current should be consistent with the global magnetospheric current. We assume that in a steady-state regime the rate of ejection of electrons from the PC surface is approximately equal to the rate of inflow of electrons from beneath the surface. It is this basic condition that makes the ejection of electrons occur without producing significant polarization of charges at the PC surface. Note that, in a static model ignoring the global current closure the free emission of electrons from the PC would rather result in a filling of the polar magnetic flux tube with electron-ion plasma and any electron acceleration will be halted.
The upper boundary condition also needs to be discussed in detail, simply because it proves to be crucial for the calculation of the fraction of returning positrons from the upper PFFs. This condition was extensively discussed in the literature (see e.g. Arons & Scharlemann, 1979;A83) and it was suggested to set E = 0 and ∇ · E = 0 at the upper PFF. The problem is that these two conditions are formally inconsistent, as we shall demonstrate in the next paragraph. This inconsistency should not be ignored in any more or less realistic model simply because these two conditions cannot be satisfied not only at the same point but even at two points separated by a distance of order of the PFF thickness.
The formal condition assumed previously that ρ = ρ GJ at the upper PFF (at η = η 0 ) means that in the very vicinity of the PFF the radial part of electrostatic potential should satisfy the following equation [see e.g. MH97, Eqn. (48)] where 2/a = θ(η 0 )η 0 1 − ǫ/η 0 ∼ r pc /R (r pc is the PC radius). The non-trivial solution of this equation reads where C = const. The boundary conditions E = 0 and ∇ · E = 0 at η = η 0 are simply reduced to the following conditions for the function P ∂P ∂η = 0 and respectively. One can easily verify that the first of the above conditions is not satisfied. This means that it is incorrect to assume that both ∇ · E = 0 and E = 0 are satisfied at the same radial coordinate. The reason is that condition ∇ · E = 0 actually implies by Eqn. (13) that the electrostatic potential also vanishes at the upper PFF. However, the latter by no means requires that E vanishes at the same boundary. Moreover, as it is illustrated in Appendix [see Eqns (A2), (A3), and (A7)-(A10)], the vanishing of both ∇ · E and E is very unlikely even within entire PFF. In general, specifying both Dirichlet and Newmann conditions (e.g. ∇ · E = 0 and E = 0) overdetermines the problem and leads to there being no solution.
In our solution (which is essential for the kind of problem we discuss in this paper) we will only require that beyond the upper PFF, where Φ ∞ is the potential at altitudes well above the upper PFF. This condition implies that there is no perfect adjustment of the effective space charge density to the GJ space charge density at the very onset of the upper PFF. When E has dropped to a low enough value that positrons are no longer able to turn around, then there is still non-zero Φ but no additional returning flux. If, at the PFF, the GJ space charge were exactly compensated by the effective space charge of primary electrons and electron-positron plasma, then the electrostatic potential would drop to zero at the PFF. The vanishing of potential at the PFF suppresses the accelerating electric field and even causes it to reverse its sign well below the PFF. The grounding of potential at the PFF would therefore dramatically affect the acceleration of primary electrons: the electron beam will be cut off from the PFF by a layer with a reverse polarity of the electric field, and the primary electrons will no longer be able to maintain the PFF they would produce if there was no strong negative feedback under discussion. As a result, it is very unlikely that a stable upper PFF would be established at all, if the condition ρ = ρ GJ sets in at the PFF. We suggest that requirement (16) is more favorable for a non-disruptive regime of acceleration of primary electrons and formation of upper pair fronts. This requirement allows the parallel component of the electric field to partially penetrate into the relativistically moving pair plasma, with the bulk of the electron-positron pairs streaming out as a quasi-neutral beam. We therefore relax the intuitive requirement that relativistically moving electron-positron pairs completely screen out the electric field at or right above the PFF. A self-consistent treatment of the dynamics of primary electrons and pairs in a nearly screened electric field, as will be discussed below, may provide us with a reasonable estimate of the fraction of returning positrons.
Using general expression (5) we can derive the formulae for the accelerating potential (field) for 0 < z < 1. This solution has the property that E saturates (E → const) for r pc /R << z < 1, before it declines approximately as η −4 (see MT92, MH97, HM98). In this case (no explicit upper boundary) the solution for the function F i (z, γ i ) is and the corresponding formulae for the potential and field read where Simple expressions can be derived from the above formulae in some limiting cases. For example, for γz ≪ 1 (z ≪ r pc /R) we get In derivation of the above expressions we have used the formulae fitting the results of numerical summation with accuracy better than 2%.
Before we proceed with our modeling of field screening and calculation of the fraction of returning positrons, let us introduce a dimensionless minimum height z 0 (≡ η 0 − 1) at which a first pair is produced . Then for our modeling we will need the expression for the electric field applicable to the region of pair formation. Since the penetration depth of the electric field into the pair plasma (or the characteristic length scale of the electric field screening) is much less than the characteristic length scale of the accelerating electric field [e.g. Eqn. (25)], the electric field in the pair region can be satisfactorily described by the formula where E acc (z 0 ) is the accelerating electric field given by (19) and calculated at z = z 0 , ∆ s is the characteristic length scale of field screening.

Analytic estimate
Let us estimate very roughly the fractional density of positrons flowing back from the PFF. We shall restrict ourselves to the case where the pair creation is mostly determined by the CR photons. The condition that the positrons with energy ε + (in mc 2 ) turn around within the PFF can be written as where E is the accelerating electric field evaluated at the PFF, ∆ s is the same as defined in Eqn.
, and e = |e| is the elementary electron charge. The characteristic energy of turning around positrons can be estimated as where ε γ is the characteristic energy of pair-producing curvature photons, n e is the number density of primary electrons, n + is the number density of turning around positrons, and dN γ /dN e is the number of photons (per primary electron) producing those pairs whose positrons turn around and flow back to the surface. Note that where p(ε) = ( √ 3e 2 /2πρ c )γF (ε/ε cr ) is the spectral power of CR (see e.g. Landau & Lifshitz, 1975), ∆ε γ is the characteristic interval of energies in the CR spectrum from which the pair-producing photons generate positrons with the characteristic energy ∼ ε + . Using Eqns (27)-(29) we can write where ε cr = 3λc γ 3 /2ρ c (λc ≡h/mc = 3.9 × 10 −11 cm is the Compton wavelength) is the critical energy of the curvature spectrum. The factor Λ (≡ ∆ε γ /ε cr ) in Eqn. (31) takes into account the fact that only photons with approximately Gaussian distribution around some characteristic energy effectively produce pairs, with ∆ε γ being the spectral interval of photons that produce returning positrons. The factor F ( , if x >> 1, where K 5/3 is the modified Bessel function of order 5/3] accounts of the fact that returning positrons are produced by CR photons with energies greater than ε cr . Finally, factor 2eγ 4 /3ρ 2 c |E | ≤ 1 describes the efficiency of emission of CR by accelerating electron. It reaches a maximum in the so-called "saturation" regime where the power of electrostatic acceleration is almost exactly compensated by the CR losses. If we adopt in Eqn. (31) Λ ≈ 3, and ε γ ∼ (3 − 6)ε cr (this range for ε γ of photons producing returning positrons agrees with that resulting from our numerical simulations, with the characteristic value of ε γ being a factor of 2-3 smaller than the mean energy of pair-producing photons given by Eqn. [29] in HM98), we get It is interesting, that the fractional density of returning positrons given by the above rough formula is practically independent of pulsar parameters, except for the dimensionless factor in brackets. The latter also tends to be close to its maximum value of unity. The reason is that in most physically interesting situations γ should be well tuned to allow magnetic pair production, which occurs at or near the saturation part of the acceleration curve, where the further acceleration of electrons would be suppressed by the CR losses, i.e. where the condition e|E | ∼ 2e 2 γ 4 /3ρ 2 c is roughly satisfied.
Even though formula (32) gives the right range for the fractional density of returning positrons, it is based on a simplified mapping between the spectrum of pair-producing CR photons and the distribution function of returning positrons. Its main deficiency is that it requires a priori knowledge (say, based on numerical simulation) of parameters Λ and ε γ . Also, it does not reflect the fact that at the upper PFF the maximum density of returning positrons is limited by the value of |ρ GJ − ρ| and is therefore dependent on the altitude of the PFF. In a steady state one can estimate the maximum value of the density of returning positrons as one half of the difference |ρ GJ (z 0 ) − ρ(z 0 )|. Thus, using expressions (2) and (4) for ρ GJ and ρ we can come up with an alternative formula A similar formula has been used by Zhang & Harding (2000) in their study of the X-ray luminosities from the spinning-down pulsars.
Let us perform numerical estimates of z 0 based on formulae (21) and (25) for the accelerating electric field. To make the resulting estimates as compact as possible, we shall adopt the following parameters: ξ = 0.5, cos χ = 1, and κ = 0.15. This is still a justified simplification, since the resultant expressions will have rather weak dependence on these parameters. Then, Eqns (21) and (25) reduce to and respectively. Here E 6 ≡ E /10 6 esu, B 12 = B 0 /10 12 G, and P 0.1 = P/0.1 s.
The altitude of the PFF above the stellar surface can be estimated as [see e.g. HM98, Eqn.

Pair source function in screening region
The pair source function Q ± (γ ± 0 , x nr ) is the joint initial energy and spatial distribution of electron-positron pairs produced by the primary electrons, where x nr is the distance above z 0 R. Only the first generation of pairs is important to the screening process, as the higher pair generations are produced beyond a distance of z 0 + N s ∆ s (N s is the number of screening scales specified in the end of Section 3.2.2). It can be shown that the attenuation length of synchrotron photons from the first pair generation is much larger than ∆ s . Each primary electron is accelerated from its starting point at the stellar surface, radiates CR photons and the pair production attenuation lengths L of these photons are computed, as described in detail by HM98. The location of the first pair defines the PFF at z 0 . We compute the pair source function by accumulating the number of pairs per primary particle at each energy and height above the PFF as the particles accelerate up to the PFF, averaged over simulations for 10 -50 primary electrons. The primary particle trajectory is divided into discrete steps and at each step the CR spectrum at energy γ(s) is divided into discrete energy bins. A representative photon from each CR energy bin is propagated through the local field to determine whether it produces a pair or escapes (for details of such a calculation, see Harding et al. 1997). We accumulate a survival probability, where is the optical depth along the path. The PC angles of ms pulsars are large enough that the curvature of the photon trajectories in the strong gravitational field of the NS is important and we take this into account in computing the pair attentuation lengths (see Harding et al. 1997, for details). A random number ℜ is chosen to determine the pair production point of each test photon, when P surv (s = L) = ℜ. The momentum of each pair member is assumed to have half the energy and the same momentum, parallel to the local magnetic field, as the parent photon. The "number" of CR photons, n CR , represented by the test photon in each energy bin is estimated by dividing the energy radiated in the spectral interval, ∆ε i = ε i max − ε i min , by the average energy in that interval, ε i , where is the CR energy spectrum, ε cr is the critical energy and ρ c the field line radius of curvature [see also the formulae right after Eqn. (31) ] andγ CR is the CR loss rateγ CR mc 2 = 2e 2 cγ 4 /3ρ 2 c . The height above the PFF and the energy of the pairs are accumulated in a two-dimensional distribution, normalized by the total number of test primary photons and by n CR , to form a distribution Q ± (γ ± 0 , x nr ) of the number of pairs per energy per primary at each height and energy interval. An example of one of the computed pair source functions is shown in Figure 1. The first pairs have the highest energies because they have the shortest attenuation lengths. The number of pairs increases rapidly beyond the PFF and the mean energy of the pairs decreases.

Pair dynamics
The charge density along open field lines above the PFF will increase due to injection of pairs. The dynamical response of the pairs to the E sc , i.e. acceleration of electrons and deceleration of positrons, will determine the increase in effective space-charge density with height in the screening region. We find that the E sc is large enough that the major contribution to the space-charge density results from positrons that are turned around and accelerated back toward the NS surface. Above the PFF, the continuity equations for elections and positrons are where n + and n − are the first moments and β + and β − are the second moments of the positron and electron distribution functions, f ± (γ ± , x), respectively: Since the source function is electron positron pairs,ṅ + =ṅ − , we can subtract Eqn. (51) from Eqn. (50) to give so that the total charge density is The distribution functions of electrons and positrons f ± (v ± , x) are computed by dynamically evolving the pair source function, subject only to the force of the electric field in the screening region. Pairs are injected at discrete points, x nr , and with a discrete distribution of energies, γ ± 0 as described below. Each sign of charge in each of the n e energy bins is evolved in energy through a separate grid, x = (z − z 0 )R, of E sc (x) computed at discrete points x n . Using Eqn. (26), we find that with the resolution of the calculation, the region in which the positrons are non-relativistic at their turn-around points is not resolved. We are therefore justified in treating the particles as relativistic everywhere and may thus ignore the momentum equation. The energy of a particle at point x in the grid is Particles from the source function, Q ± (γ ± 0 , x nr ) are evolved through the grid until they either travel upward across the top boundary or downward across the lower boundary. At each grid height, the steady-state distribution function of electrons and positrons are computed Using Eqns (52), (53), (55) and (58), we can determine the density of upward-moving electrons, n − u (x), and positrons, n + u (x) and the downward-moving positrons, n + d (x) at each grid point. The total charge density due to pairs is then The downward-moving positrons will constitute an additional upward-moving negative current which will add to the current of primary elections (Section 2.2). Because the net current j is required to be the GJ current, we therefore need to adjust the primary current for the returning positron current. Thus, we can write where n − p is the number density of primary electrons and n + d (1) is the downward-moving positron density at the grid lower boundary (the PFF). The charge deficit, ∆ρ = ρ − p − ρ GJ at the PFF is thus also readjusted.
The calculation thus proceeds as follows. We first compute the location of the PFF above the surface, z 0 , and the pair source functions as described in Section 3.2.1. An initial screening scale length, ∆ s , is chosen as a fraction of z 0 , the distance to the PFF. The electric field in the screening region, above the PFF, is modeled by an exponential (Eqn. [26]) with scale length ∆ s . The pair dynamics determines a pair distribution function and a returning positron density, as described above. The primary electron flux and ∆ρ are adjusted for this returning positron flux, and the pair distribution function and a returning positron density are then recomputed using this adjusted charge deficit. A new screening scale length is determined from the point where ρ(N s ∆ s ) = ∆ρ, where N s is the number of E sc (x) scale heights required to guarantee that no positrons turn around above z 0 + N s ∆ s . The solution will then be self-consistent. We find that N s = 4.0 satisfies this condition and set it as a constant value. The scale height of the electric field in the screening region is then set to the new ∆ s . The iteration continues until the returning positron density and screening scale attain stable values.

Screening scale height and returning positron density
Convergence of the screening scale height ∆ s , the total charge density above the PFF, ρ(x)/ρ GJ , and the relative returning positron density, ρ + d (1), is achieved independently at each colatitude, ξ, across the PC. An example of a self-consistent solution for ρ(x)/ρ GJ and ∆ s is shown in Figure 2. The charge density increases rapidly immediately above the PFF, where the pair source function is growing exponentially and where E sc is high enough to decelerate and turn around all the produced positrons. At a distance above the PFF comparable to the screening scale height, E sc has decreased to a fraction of its value at the PFF and the rate of increase of the charge density begins to moderate, as fewer positrons are able to turn around. We find that the creation of the charge density and subsequent screening of E sc is due almost entirely to returned positrons rather than to velocity differences between electrons and positrons as the pairs dynamically respond to the electric field above the PFF. Due to the boundary condition imposed on the potential (see Eqn. [16]), E sc never decreases to 0, but only approaches 0 at infinity. Consequently, the self-consistent charge density never achieves the value ∆ρ but approaches it asympotically from below. Figure 3 shows solutions for the screening scale, ∆ s , and the returning positron density, ρ + /ρ GJ , across the PC as a function of the colatitude ξ, scaled to the PC opening angle at the surface, for various values of the pulsar period and surface magnetic field strength. Both ∆ s and ρ + /ρ GJ increase toward ξ = 0 (magnetic pole) and ξ = 1 (the outer edge of the PC). Since the value of the GJ density, ρ GJ is constant with ξ, the variation represents a true variation in ρ + across the PC. The increase in ∆ s and ρ + toward ξ = 0 is due to the increasing field line radius of curvature near the pole, causing the pair attenuation length to grow large. The pair source function grows more slowly with distance above the PFF, and the screening scale height increases. The returning positron density increases because as the height of the PFF increases, the charge deficit also increases, requiring more returning positrons for screening. Near the pole, there are no solutions for ∆ s and ρ + and screening is either incomplete or there is no screening at all. The increase in ∆ s and ρ + toward ξ = 1 is due to the decrease in E caused by the boundary condition Φ = 0 at ξ = 1 imposed on the solution to Poisson's equation. Consequently, the primary particles must accelerate over a longer distance to produce pairs, increasing the height of the PFF and thus ρ + . Also the positron turn-around distance becomes longer with a lower E , increasing the screening scale height.
However, the variation in ∆ s and ρ + (ξ) with ξ is small compared to their variation with pulsar period and surface field strength. Comparing Figures 3b,3c, and 3d, one can see that the screening scale height ∆ s increases by a factor of about ten as the period increases by only a factor of four, for a constant surface field. Comparing Figures 3a and 3c at a period of 0.1 s, ∆ s decreases by a factor of about six as the surface field increases by only a factor of two. Figure 4 shows the dependence of the screening scale length, ∆ s , on surface field and period. Generally, ∆ s decreases with increasing surface field strength because the pair attenuation length is shorter in higher fields. As a result, the pair density grows faster and E can be screened in a shorter distance above the PFF. Consequently, ms pulsars have relatively large screening scale lengths, but ∆ s is still small compared to z 0 . Figure 5 shows the dependence of ρ + on surface field and period. While ρ + is less sensitive to pulsar parameters than ∆ s , it decreases significantly as the surface field increases and as the period decreases. The numerical values of ρ + /ρ GJ are smaller than, but within a factor of two of, the analytic estimate of Eqn. (45). The variation of ∆ s and ρ + with pulsar parameters is primarily due to the dependence of E and the pair production attenuation length on period and surface field strength. By contrast, ∆ s and ρ + (ξ) vary by only factors of two or three across the PC of a single pulsar.
Since ρ + is always a small fraction of both the primary particle density and of ρ GJ , only the very first pairs of the full, multi-generation cascade are needed to screen the E . We therefore did not simulate the full cascade in computing the pair source functions needed in the present work. The vast majority of the cascade pairs are produced in the region where E ∼ 0 and they freely escape the magnetosphere with only radiation losses, but no significant acceleration. Previous simulations of the full cascade produced by the primary particles (Daugherty & Harding 1996) give multiplicities ∼ 10 3 − 10 4 pairs/primary. The fraction of the total number of secondary pairs that are returned to heat the stellar surface can thus be estimated by dividing the values of ρ + /ρ GJ in Figure 5 by the cascade multiplicity.

Analytic estimate
Let us derive an expression for the total power that can be deposited onto a single PC by precipitating positrons. The general expression can be written as where the integration is over the area of a sphere cut by the polar flux tube at the radial distance η 0 , and the factor α accounts for the general-relativistic correction to the current (j ∼ αcρ). Here dS = [S pc η 3 0 f (1)/πf (η 0 )]dΩ ξ , S pc = πΩR 3 /cf (1) is the area of the PC, dΩ ξ = ξdξdφ is an element of a solid angle in the PC region. Thus expression (61) reduces to After inserting expression for Φ [see Eqns (20), (24)] into (62) and normalizing ρ + by ρ GJ , we get where and are the fractions of pulsar spin-down power consumed by returning positrons, corresponding to the cases where Φ is given by Eqn. (20) and (24), respectively. Herė is the general-relativistic expression for the pulsar spin-down losses in vacuum (see MH97). In this formula B 0 /f (1) is the surface value of the magnetic field strength "as seen" by an infinitely remote observer, whereas B 0 is the value measured locally, at the NS surface. In a flat-space limit B 0 /f (1) simply transforms into B 0 , and we get a classical formula for the magneto-dipole losses. Note that in our derivation of formulae (64) and (65) we used only the components of potential and charge density that are proportional to cos χ and we assumed that z 0 is independent of ξ (flat PFFs), which is still a satisfactory approximation for rough analytic estimates. Our numerical calculations of f + take into account the curvature of PFFs (see further discussion of the numerical results), though. Now, using expression (41) for z 0 and (45) for ρ + /ρ GJ , we can combine Eqns (64) and (65) This expression illustrates that there is a rather strong correlation between the efficiencies of PC heating by precipitating positrons and pulsar spin-down age. The important implication of this correlation is that older pulsars should favour enhanced PC heating (see, however, the warning on use of Eqn. [67] in Section 4.2). It is interesting that the above order-of-magnitude analytic estimates are in good agreement with our numerical calculation.
A81 derived an expression for PC heating luminosity, taking into account the curvature of field lines but not general-relativistic effects in the derivation of E . Using his results we can derive the following expressions for f + The value of f A + in Arons' model also increases with pulsar age but is significantly lower than our result. This is because the frame dragging electric field and consequent acceleration energy in our model is higher.
During the photon cooling era, the effective temperature of a NS decreases with age according to a power law. For example, using the cooling calculations by Page and Sarmiento (1996), we can estimate that the cooling luminosity, L cool , roughly scales as τ −4 6 (τ 6 ≥ 0.3) and τ −3 6 (τ 6 ≥ 1) for the case with and without core neutron 3 P 2 pairing, respectively. By normalizing the luminosity of a cooling NS byĖ rot , we get the following expression for the cooling efficiency f cool = L cool E rot ≈ 10 −4 τ 6 0.15 (P 0.1 /τ 2 6 ) 2 if core neutron 3 P 2 − pairing and τ 6 ≥ 0.3, 2.5 P 2 0.1 /τ 3 6 if no core neutron 3 P 2 − pairing and τ ≥ 1.
¿From expressions (67) and (69) one can see that the pulsars with ages ≥ 10 6 yr are most likely candidates for those with luminosities dominated by the PC heating.

Numerical calculation of polar cap heating luminosity
The returning positrons accelerate through the same potential drop as the primary electrons and their energy heats the PC surface. We can calculate the luminosity heating the PC by using expression (62) and substituting the corresponding quantities from the numerical calculation. Figure 6 shows the dependence of L + (ξ) on ξ (from the integrand of Eqn. [62]), which reflects the distribution of heating across the PC. Although the distribution of ρ + (z 0 , ξ) has maxima at ξ = 0 and ξ = 1 (see Figure 3), γ max (ξ) has a strong maximum at ξ = 0 and decreases monotonically with ξ (see Figure 4 and 5 of HM98), causing L + to decrease with ξ.
The total positron heating luminosity, scaled with the spin-down luminosity as a function of characteristic pulsar age, P/2Ṗ is shown in Figure 7. The numerically computed L + /Ė rot increases nearly linearly with τ in the same way as the analytic estimate in the unsaturated regime (first expression of Eqn. [67]). L + /Ė rot also increases approximately linearly with period for a constant age. The PC heating luminosity is thus a negligible fraction of the spin-down luminosity in young pulsars with ages less than τ ∼ 10 5 yr, but becomes a significant source of heating in older pulsars. The numerical values are a factor of ∼ 6 lower than the estimated values of L + /Ė rot from the first expression of Eqn. (67) for the normal pulsars. For the ms pulsars, the analytic values are somewhat higher than the numerical values, but by less than a factor of ten. The reason that the analytic estimate, which simply assumes that the returning positron flux is half of the charge deficit at the PFF, is closer to the numerical values for normal-period pulsars may be that the screening scale heights are small and the returning positron fraction is not very sensitive to the structure of the electric field in the screening region. On the other hand the screening scale heights for the ms pulsars are about an order of magnitude larger, reflecting the larger (see Figure  4) growth scale of the pair source function. Consequently, the returning positron fraction for ms pulsars is more sensitive to the structure of the electric field in the screening region.
Although the analytic estimate of Eqn. (67) seems to be a very good estimate of the returning positron luminosity, we caution that it cannot be extrapolated to ages much beyond the numerical results, i.e. above τ ∼ 10 7 yr for normal pulsars and above τ ∼ 10 8 yr for ms pulsars. This is because the assumption of complete screening breaks down at age-period combinations where the pulsar cannot produce enough pairs. When the pair density grows too slowly and cannot reach a high enough level for complete screening, then the returning positron fraction drops below the charge deficit prediction. Beyond the upper end of each constant period line in Figure 7, numerical solutions of returning positron fraction and screening scale length do not exist. In fact, complete CR screening terminates on the constant period line at about the same point where PFFs are no longer produced. Beyond this point, the returning positron fraction and the heating luminosity will drop sharply. However, as we will discuss in our next paper, ICS pair fronts are produced by older pulsars, even when the ICS screening from these pairs is incomplete. This will result in a PC heating component due to returning ICS positrons from partial screening which can be quite significant for ms pulsars.
The results in Figure 7 have assumed a constant value of inclination angle, χ = 0.5. The expected variation of f + ∝ cos p χ (where p = 5/7 or 1/2) will cause f + (= L + /Ė rot ) to gradually decrease with increasing χ until χ approaches π/2. At this point, the second term in Eqn. (5) for the electrostatic potential (the relativistic equivalent of A83 solution) becomes important. For χ = π/2, our solution for the PC heating rate will be comparable to that of Eqn. (68). X-ray emission at energies 0.1 -2 keV has been detected from several dozen pulsars by ROSAT (Becker & Trumper 1997). There is a rough empirical correlation of X-ray luminosity with spin-down luminosity, giving L x ∼ 10 −3Ė rot , so that the observed level of X-ray emission in pulsars would be a constant line in Figure 7 at 10 −3 . The maximum calculated values of L + /Ė rot do reach the observed level, indicating the emission from PC heating in normal pulsars with τ > 10 5 yr and in ms pulsars is detectable.
In Table 1, we compare our computed values of flux and surface temperature from PC heating with measured values of some older pulsars in which hot thermal components have been detected. Our values are those at the NS surface and have not been corrected for gravitational redshift effects -i.e. L + and Φ + would be about 40% lower (a factor of α 2 ) and T + would be about 20% lower (a factor of α) for observers at infinity. Greiveldinger et al.(1996) have fit three-component (two thermal and one power-law) spectra to combined ROSAT and ASCA data from PSR0656+14 and PSR1055-52. These pulsars are young enough to have expected cooling as well as heating thermal components, so it is important to separate the emission from these two components in order to make comparisons with PC heating models. PSR0656+14, however, has an inferred surface dipole field of 9 × 10 12 G and so is not in the regime of CR pair fronts near the NS surface that we have treated in this paper. It may still have a CR pair front at higher altitude if the ICS pair front is unstable, a situation we will address in our next paper. The fits for the heated area A of the hot components in both cases are much smaller than the NS surface area, and are thus consistent with emission from heated PCs. Our predicted values of PC temperature from returning positrons agree fairly well with the measured values for these pulsars, although our predicted fluxes are significantly lower than those observed. Because these middle-aged pulsars are still dominated by cooling components and have power law components at high energies, extracting the relatively small heating component is very difficult.  have fit single-component blackbody spectra to ASCA data of PSR1929+10. This pulsar is too old to have detectable cooling component and the one-component fits of ASCA data indeed indicate emission from heated PC whose heated area is much smaller than even the standard PC area. We have thus assumed the measured value of area A in computing the theoretical PC temperatures, T + = (L + /σA) 1/4 , of these pulsars, instead of the canonical PC area, A pc = πR(ΩR/c). Our computed temperatures are in agreement with the measured values, within the uncertainties. Long period (older) pulsars have smaller PCs and thus the returning positron luminosity will heat a smaller area and we therefore would predict that pulsars with longer periods should have higher PC temperatures. PSR1929+10 does in fact have higher measured temperature and smaller heated area. Although Geminga is expected to have both cooling and heating components, it is more difficult to compare our theoretical results for positron heating with the observed values from a single-component thermal fit. The fact that our predicted flux is somewhat lower than the measured value is consistent with an additional cooling component. We have used the canonical PC area to compute our predicted temperature for Geminga, so we would expect it to be higher than a measured temperature which includes cooling from the whole NS surface.

SUMMARY AND CONCLUSIONS
In this paper we derived some practical formulae for the positron fluxes returning from the upper PFF in rotation-powered pulsars in cases where E screening is produced by pairs from CR. We presented the expected theoretical values for the PC X-ray luminosities due to the heating by precipitating positrons. The calculated efficiencies of PC heating explicitly depend on the pulsar spin-down age and spin period, and can be used in the analyses of thermal X-ray fluxes from pulsars. Our numerical calculation of returning positron flux and PC heating efficiencies show that the dependence on τ and P is the same as in analytic formulae, the actual values of the numerical quantities are smaller than the analytic estimates by a factor of 2-3. In summary, we have reached the following conclusions.
• The heating of the PC by returning positrons is possible if stable pair fronts develop, and E is not assumed to be fully screened at the onset of the PFF. In most cases we discuss in this paper the fraction of returning positrons is not affected by the details of our modeling of the screening of E and is mostly determined by the distribution of pairs beyond the PFF.
• In contrast to the results obtained earlier by A81, we find that the returning positrons may significantly contribute to the PC heating and therefore to the observed pulsar thermal X-ray fluxes, especially for older pulsars.
• Our theoretical model can be thoroughly tested against observations as soon as more data on the pulsed thermal X-ray fluxes from middle-aged and old pulsars become available. We anticipate that the pulsar X-ray light curve would have larger pulse fraction in the case of NS thermal emission produced by the PC heating than in the case of thermal emission during NS photon-era cooling.
• We predict that long-period pulsars will have higher surface temperatures from PC heating than those of short-period pulsars.
• Our calculations indicate that the PC heating increases toward the magnetic axis, which could manifest itself in a fine structure of the X-ray light curves and heated areas smaller than the PC area.
We caution that our analytic expression in Eqn. (67) for f + is not applicable for normal pulsars with τ > ∼ 10 7 yr or ms pulsars with τ > ∼ 10 8 yr. The reason is that the analytic formulae are based on constant magnetic field and do not take into account the cessation of pair formation by CR photons beyond a certain age (dependent on period). Thus, our plots (see Figure 7) for pulsars with 5 and 10 ms periods cannot be compared with the observational data for actual ms pulsars, most of which have ages τ > ∼ 10 8 yr where CR photons cannot produce pair fronts. Although our analytic formula for f + would give very large heating luminosities for these ms pulsars (beyond the region of validity of the formula), the actual heating efficiencies will be much lower. However, as we will address in our next study, ms pulsars are capable of producing pairs via non-resonant ICS and the resulting heating luminosities will be detectable. Furthermore, the results presented in this paper are not applicable to pulsars with surface magnetic fields above about 4 × 10 12 G, which are capable of ICS screening. The possibility of ICS pairs from those pulsars which are beyond their CR deathlines, as well as ICS screening and PC heating in higher-field pulsar, will be investigated in our next paper.