{'fraction_non_alphanumeric': 0.03675158202174266, 'fraction_numerical': 0.01509005354535129, 'mean_word_length': 4.251952847606874, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Pair-production multiplicities in rotation-powered pulsars

We discuss the creation of electron-positron cascades in the context of pulsar polar cap acceleration models and derive several useful analytic and semi-analytic results for the spatial extent and energy response of the cascade. Instead of Monte Carlo simulations, we use an integro- differential equation which describes the development of the cascade energy spectrum in one space dimension quite well, when it is compared to existing Monte Carlo models. We reduce this full equation to a single integral equation, from which we can derive useful results, such as the energy loss between successive generations of photons and the spectral index of the response. We find that a simple analytic formula represents the pair cascade multiplicity quite well, provided that the magnetic field is below 10^12 Gauss, and that an only slightly more complex formula matches the numerically-calculated cascade at all other field strengths. Using these results, we find that cascades triggered by gamma rays emitted through inverse Compton scattering of thermal photons from the neutron star's surface, both resonant and non-resonant, are important for the dynamics of the polar cap region in many pulsars.

The polar-cap pair-production model of pulsar emission has remained the chief theory of the region thought to give rise to the radiative emission for the past 30 years. Particles accelerate from the surface of the neutron star, drawn by the induced fields of the rotating magnetosphere. As first discussed by Sturrock (1971), these particles emit γ-rays which pair-produce in the high background magnetic field. These pairs then short out the magnetic field (Ruderman & Sutherland 1975), preventing further acceleration.
The physics of this pair production process has been partially explored, with progress ranging from the initial generational models of Tademaru (1973), the pair formation front structure calculations of Arons & Scharlemann (1979); Arons (1981), the detailed Monte Carlo simulations of Daugherty & Harding (1982), and the more recent "full-cascade" generational model of Zhang & Harding (2000). Our paper describes an improved analytical and theoretical understanding of this pair-production mechanism, not only reproducing the final γ-ray and pair spectra calculated by Daugherty & Harding (1982), but detailing the spatial extent of the pair-production process itself, giving a more thorough understanding of the pro-cess.
These straightforward analytic approximations to the pair production process were the basis of the model of Hibschman & Arons (2001), henceforth called Paper I. In this paper, we exhibit the full shower theory underlying those approximations and find a variety of formulae which represent the numerical solutions of the cascade equations, thus allowing users of the cascade physics to easily apply the results to other problems. We give a full description of the pair plasma emergent from pulsar polar caps. We reserve discussion of the emergent γ-ray luminosities and spectra to a separate paper, as these are of direct observational relevance.

gamma-ray opacity
The dominant opacity for pair production over pulsar polar caps is the γ-B process, wherein a high-energy photon interacts with the background magnetic field to form an electron-positron pair. The competing process of γ-γ pair production, wherein a high-energy photon interacts with a background X-ray to form a pair, is only competitive if the background X-ray flux is equivalent to a black-body with temperature nearly 10 7 K, well beyond the observed limits on stellar temperatures. According to Erber (1966), the opacity for γ-B pairproduction is where B q is the critical quantum magnetic field, B q = e/αλ 2 C = 4.41 × 10 13 Gauss, α F = e 2 / c ≈ 1/137 is the fine structure constant, λ C = /mc = 3.86 × 10 −11 cm is the reduced Compton wavelength, B is the local magnetic field strength, ψ is the pitch angle between the photon momentum and the local magnetic field, and where ǫ is the photon energy, in units of mc 2 . For the remainder of this paper, all energies will be quoted in terms of mc 2 , for convenience. This expression is accurate provided that ǫ sin ψ > 1, so that the created pair is in a high Landau level, and provided that B is small compared to B q ; B < B q /3 suffices. For pair production into low Landau levels, the full cross section must be used (Daugherty & Harding 1983), and for magnetic fields equal to or higher than B q , higher-order corrections become important, as discussed in Harding, et al. (1997). Most pulsars, however, pair produce into high Landau levels for the most significant photons and have B < B q /3. As this work intends to describe the "typical" pulsar, we neglect these high-energy and high-field effects.

Optical depth
Because the relativistic primary electron beam follows the magnetic field, the primary photons are beamed close to parallel to the local field. If the beam electrons have a Lorentz factor γ, the primary photons are beamed into a cone of opening angle ψ ∼ 1/γ. Typical beam Lorentz factors are of order 10 3 or higher, giving an initial pitch angle small compared to the pitch angles required for pairproduction. To a good approximation, then, we can treat the photons as if they are injected precisely parallel to the field lines.
As the photons propagate through the magnetosphere, the pitch angle between the photon momentum and the background magnetic field steadily increases. This pitch angle is given by sin ψ = |k ×B(r)|, wherek is the photon momentum direction, andB is the local magnetic field direction. If we assume ψ small enough that sin ψ ≈ ψ, the change in ψ along the photon path is then dψ/ds = |k × (k · ∇)B(r)|.
Since ψ is small,k remains close toB, so we can substitute, yielding dψ/ds ≈ |B × (B · ∇)B(r)| = ρ(r) −1 , where ρ(r) is the local magnetic field radius of curvature. This is accurate to second order in ψ.
For a dipole field, the photon pitch angle depends on radius (through second order in the colatitude, θ) as where r e is the radius of emission and r is the current radius, both measured from the center of the dipole, and ψ ∞ ≡ r e /ρ e , where ρ e = (4/3)(R * /θ * )f ρ (r e /R * ) 3/2 is the field line radius of the magnetic colatitude of the field line in question at the surface of the star and f ρ ≤ 1 is a factor used to crudely take into account the possibility of nondipolar components to the magnetic field at low altitude. Numerically this is ψ ∞ = 0.011(θ * /θ C )P −1/2 (r/R * ) 1/2 f ρ , with θ C the colatitude of the last closed field line, θ C = R * /R L and P is the period in seconds. If the dipole moment is displaced close to the surface of the star, GR bending of the photon paths would change this, but for almost all other cases GR perturbations are negligible.
The magnetic field and radius of curvature may then be rewritten as where B e and ρ e are respectively the magnetic field and radius of curvature at the emission point.
Using these relations to write the optical depth using t ≡ ψ/ψ ∞ as the integration variable, letting sin ψ ≈ ψ, yields and ǫ is the photon energy in units of mc 2 . The exponential above has a maximum at t = 1/4. Physically, this is the point where the increase in opacity due to the increasing pitch angle balances the decrease in opacity due to the steadily declining magnetic field.
Provided that Z ∞ is not too small, the integrand is a sharply peaked function of t, and can therefore be integrated using the steepest descents method, yielding When compared to numerical calculations, this expression is accurate to better than 10% in the neighborhood of τ = 1, for the range of magnetic fields and radii of curvature found in pulsars. Photons with τ ∞ < 1 will escape the magnetosphere, while those with τ ∞ > 1 will be absorbed. If τ ∞ (ǫ) is large, most of the absorption takes place on the leading edge of the exponential. Since the opacity is increasing exponentially in that regime, we can approximate the integral by expanding around the upper endpoint to find and where f (t) has been defined to be t −1 (1 − t) −3 for t < 1/4 and 256/27 for t > 1/4. (This result was previously derived in Arons & Scharlemann (1979).) Since the opacity decreases rapidly after t = 1/4, the optical depth saturates at that point.

Absorption peak
Given this opacity, we can find where any given photon will be absorbed. Due to the exponential dependence of the opacity, all photons are effectively absorbed at the maximum of α exp(−τ ), or approximately at τ = 1. Using (9) for τ , we find this peak occurs at where ln Λ ≡ ln(Λ 1 (ǫ)t 3 ). In the limit of small t a , this becomes Since the opacity saturates past ψ a = ψ ∞ /4, the minimum photon energy required to pair produce is This energy should be thought of as a critical scaling energy, not as the actual minimum energy which will pair produce, since the (1 − t) term was neglected. The actual minimum energy which will pair produce is As long as ǫ 5ǫ a , the small-t limit (12) is appropriate, and we find that a photon will be absorbed after propagating These expressions can be made tractable by treating ln Λ as a constant. Self-consistently evaluating ln Λ at t = 1/8 then yields (Arons & Scharlemann 1979). ln Λ = 16.2 + ln B 12 − 1 2 ln P.
3. synchrotron emission A charged particle with pitch angle ψ with respect to the magnetic field emits synchrotron radiation and rapidly spirals down to its lowest Landau level.
The basic rate of synchrotron emission is where ǫ B ≡ B/B q , K 5/3 is the modified Bessel function of order 5/3, and ǫ s is the characteristic synchrotron energy, This corresponds to a total power loss of (Jackson 1975) Due to synchrotron radiation, the Lorentz factor of the particle decreases according toγ = −P s . This energy loss occurs over a distance of 1.8 × 10 −6 γ cm, effectively instantaneous compared to stellar scales, where we have used equation (24) self-consistently to fix γψ and let sin ψ ≈ ψ. The final Lorentz factor may be found by noticing that the parallel component of the particle's velocity is conserved, as can be seen by transforming into a comoving frame. If the initial particle has a Lorentz factor of γ i = (1−β 2 i ) −1/2 and is moving at an angle ψ with respect to the magnetic field, then the final Lorentz factor is Integrating (17) over all time gives the number of photons generated during the decay. Using (19) to change the variable of integration from time to γ gives Changing variables to z = ǫ/ǫ s (γ) and re-arranging gives where Equation (22) gives the total number of photons of energy ǫ emitted by a particle with initial Lorentz factor γ i and initial pitch angle ψ. Nearly all of these photons are emitted in a cone with opening angle ψ, because the Lorentz factor of the particle remains high, so the angle ψ only changes once almost all of the particle energy has been radiated away.
In the case of a particle created by γ-B pair production from a photon of initial energy ǫ i , γ i = ǫ i /2 and ψ = ψ a (ǫ i ). Defining we find z i = ln Λǫ/ǫ i and z f = (1 + a 2 ) ln Λǫ/ǫ i . With these substitutions and counting the emission from both generated particles, the total number of synchrotron photons produced by one incident photon pair producing at ψ = ψ a is We can understand the synchrotron spectrum by examining its asymptotic limits. Figure 1 shows the function F (t) and its deviation from its zero-point value. The function is nearly constant for t < 1, with F (0) − F (t) ∝ t 5/6 , and exponentially decreases for t < 1. This divides the synchrotron response into three regimes. If ǫ < ǫ i /(1 + a 2 ) ln Λ, both z i and z f are less than one, and N ǫ (ǫ i ) ∝ ǫ −2/3 ; this is the usual synchrotron spectral index. If ǫ i /(1 + a 2 ) ln Λ < ǫ < ǫ i / ln Λ, z i < 1 but z f > 1, and N ǫ (ǫ i ) ∝ ǫ −3/2 ; this reflects the decline of the particle Lorentz factor. If ǫ > ǫ i / ln Λ, both z i and z f are greater than 1, and N ǫ (ǫ i ) decays exponentially. -Synchrotron response function. For t < 1 the function is roughly constant and F (0) − F (t) is a power law with index 5/6, while it exponentially decays at t > 1. This produces the synchrotron spectral indices of −2/3 at low energies and −3/2 at moderate energies, and the exponential tail at high energy.

Integral equation
Since particles are assumed to be beamed along the field lines, the specific number intensity at any point is where r e is a function of r and ψ via equation (3), Q γ is the volumetric photon emission rate, ρ is the radius of curvature of the field line, and τ (r e , r, ǫ) is the optical depth attained by a photon of energy ǫ propagating from r e to r. The radius of curvature enters since ρdψ is the path length where the field points in direction dψ.
The local volumetric rate of pair production events is then Using the synchrotron response, equation (25), the synchrotron emission rate is (28) Using equation (26) for I γ , dτ = α B ρdψ, and the sharppeaked nature of the absorption, we can exactly evaluate the ψ integral, yielding (29) wherer e = r e (r,ψ),ψ is the peak of the ψ integral, and τ max is the optical depth attained by a photon of energy ǫ propagating from the stellar surface to r.
Since, by equation (15), a photon emitted at r e travels at most 0.25r e before pair producing, we can treat the cascade as if it occurred on-the-spot (OTS). Photons with energies well above ǫ min pair-produce in a correspondingly shorter distance, so the photons near ǫ min set the spatial expanse of the cascade. This yields an integral equation for the synchrotron spectrum, where we have replaced τ max with τ ∞ and where This can be easily solved by matrix operator methods. If Q γ,syn and Q γ,src are represented by vectorsQ syn and Q src , and the integral operator above by a matrixK, the solution is simplyQ syn This gives an excellent approximation to the final γ-ray spectrum, and a close approximation to the final pair spectrum, as shown in Figure 2. The numerically calculated pair spectrum extends to lower energies than the on-thespot (OTS) spectrum, due to the change in magnetic field over the course of the cascade.

Moment equations
From the integral equation for the synchrotron photon distribution, we find several useful relations connecting the moments of the photon spectrum. Using equation (30) to calculate successive generations of synchrotron photons, we find where Q (i) γ is the i th -generation γ-ray spectrum and where we have replaced 1 − exp(−τ ∞ ) with a step function at ǫ min .
Multiplying by ǫ n and integrating the above expression over all energies gives where R n is the "reduced moment" of the spectrum, neglecting energies below ǫ min . Since K is known, the K n are known. Given a known source function, R n may be calculated. The moments of successive generations may then be calculated by assuming that R n ≈ Q n so that Q n . The limited number of generations truncates the sum of successive moments. Integrating the kernel, we find that the first two moments are In general, the n th moment is From K 0 and K 1 respectively, we can see that the total number of photons steadily increases while almost all of the energy remains in the photon spectrum, provided that φ = (1 + a 2 ) is large. The fraction of energy left in the pair spectrum is 1 − K 1 = φ −1/2 , which becomes significant for values of φ near 1, which occurs at magnetic fields greater than 3 × 10 12 Gauss.

Generations
A cascade of pair production then divides into several successive generations of photons converting into pairs which emit more photons. Following the work of Tademaru (1973) and later papers (Lu, et al. 1994;Wei, et al. 1997;Zhang & Harding 2000), we find, in our notation, that the characteristic synchrotron energies of successive generations are related via while the particles produced by generation j have final Lorentz factors where the approximate form assumes a is well greater than 1, which is true provided that B < 3 × 10 12 Gauss. These correspond exactly to the expressions cited by Zhang & Harding (2000), given the equivalence ln Λ = 4/3χ. From this, it is clear that for each absorbed photon, a fraction of the energy, is immediately re-emitted in lower-energy photons, with the remainder left in the electron-positron pair. As discussed by Zhang & Harding (2000), much of the the energy in the pairs is later radiated at energies near 2ǫ b γ f via resonant inverse Compton scattering, but this process generally does not create further pairs. These approximations effectively collapse the entire synchrotron emission of a created pair into a pulse of photons at the peak energy of the synchrotron spectrum. By using the moments of the photon spectra, we can derive similar approximations, but with wider applicability.
At each generation, the number of photons increases by a factor of K 0 , while the total energy in the photons decreases by a factor of K 1 , so the average energy of the j th generation isǭ whereǭ 0 is the average energy of the input spectrum. The total number of photons in the j th generation is where N 0 is the initial number of photons. Provided that ǫ j > ǫ a , most of the photons will convert into pairs. If we truncate the generations when the average energyǭ j reaches ǫ a , we find a power-law relationship between the initial energy and the final pair multiplicity, where the power-law exponent is given by The constant of variation depends on the behavior of the lowest-energy photons near ǫ min and will be discussed in Section 5.2. We can apply similar arguments to the spatial extent of the cascade. From equation (15), the photons of generation j, with mean energyǭ j , must propagate a distance ∆s = 0.25(ǫ a /ǫ j )R * before pair-producing. Relating the multiplicity of generation j, equation (44), to this propagation distance gives a power-law dependence between the multiplicity and the distance, M (s) ∝ s ν (47) where the power-law exponent, ν, is identical to that of the energy-multiplicity relation, equation (46).

numerical method
Since the OTS method cannot spatially resolve the pair cascade, we turn to a numerical calculation. Equation (29) provides the fundamental structure for the calculation. It gives the synchrotron spectrum injected into the magnetosphere at any point, as a function of the photons absorbed at that point.
To calculate the resultant spectrum, we simply evaluate the photon spectrum Q(ǫ) at logarithmically spaced points in ǫ. The relation between the synchrotron photons output and the absorbed spectrum at a given radius is then a simple matrix multiplication by N ǫ (ǫ ′ ). Due to the variation of magnetic field with altitude, this matrix will also depend on the altitude.
We then discretize the radial steps on a variable scale, with either purely logarithmic steps, in cases where we are primarily concerned with the γ-ray output and the location of the PFF, or a combination of logarithmic steps, up to a polar-cap radius, and linear steps beyond, to calculate a smooth pair spectrum. Since the final pair energies depend on the local magnetic field, the pair spectrum is more sensitive to the spatial variations in the field, requiring this scheme.
To begin the calculation, at each radial point we inject an initial spectrum. These spectra are either a deltafunction source at the lowest radius, for measuring the single-photon response, or the expected spectra emitted by a single beam particle as it passes through each radial bin, for a full physical cascade. In the second case, the spectrum is adjusted to prevent binning effects from affecting the total power injected.
Using the expression τ ∞ , equation (8), we find which photons are absorbed within the magnetosphere, and which escape. The escaping photons are simply accumulated. The absorbed photons are assumed to be absorbed at ψ(ǫ) = ψ a (ǫ), or at a radius of r ′ = r + ρ(r)ψ a . Each photon energy from each radial bin is absorbed at a different height; these final heights are then simply re-binned into the grid.
For each absorbed photon, equation (25) is used to determine the synchrotron spectrum injected at the absorption point. This process is then iterated, to find the secondaries produced by those synchrotron photons, and so on, until the energy remaining in synchrotron photons is negligible.
These calculations are streamlined by pre-calculating the synchrotron response matrices, N ǫ (ǫ i ; r), and the total absorption, τ ∞ (ǫ; r), at each altitude. In addition, the radial absorption matrix, r ′ (r; ǫ), is computed for each energy, reducing each step of the iteration to a series of matrix multiplications. Figures 4 and 5 show, respectively, the final γ-ray and pair spectra for both this model and the Monte Carlo model of Daugherty & Harding (1982). The predicted γray spectra match extremely well, while the pair spectra mismatch at low energies, due to the low resolution of the Monte Carlo method in that range, and at high energies, due to the approximation in this model of treating ln Λ as a constant, rather than allowing it to vary with photon energy.

single-photon results
To understand the cascade process itself, we first examine the response to a single photon of injected into the magnetosphere at the surface of the star. Physically, the resultant cascade depends on the energy of the input photon, the local magnetic field strength, and the local radius of curvature of the field lines.

Gamma-ray spectra
From the discussion of the optical depth, we know that only photons emitted at radius r with energies less than ǫ min (r) will escape the magnetosphere. Due to the decline of the magnetic field with r, ǫ min steadily increases with altitude.
However, the cascade produced by a single injected highenergy photon remains localized. If there are multiple generations of pair production, the photon energy must be greater than ln Λǫ a , from equation (40), which pairproduces in ∆r < r e /4 ln Λ. Since ln Λ is on the order of 20, this distance is short compared to r e , and the magnetic field remains effectively unchanged. Hence, we may treat the cascade process as if it occurred on-the-spot, and the energy cut-off remains at ǫ min (r e ).
The typical synchrotron response has a power-law exponent ν = −3/2, extending from energy ǫ i / ln Λ to ǫ i /(1 + a 2 ) ln Λ, and a power-law exponent of ν = −2/3, the typical value for synchrotron radiation, for lower energies. After processing through multiple generations, this translates to an output spectrum with ν = −2/3 up to ǫ min /(1 + a 2 ) ln Λ and ν = −3/2 from there to ǫ min . Figure 6 shows the gamma-ray spectra produced by various input energies. This demonstrates the expected cut-off at ǫ min , as well as the transition from a power-law exponent of −3/2 to −2/3 at ǫ min /(1 + a 2 ) ln Λ. The shape of the spectrum is different at low energies, due to the effects of the exponential tails of the opacity and synchrotron emission, but it quickly converges to the asymptotic form.

Pair multiplicity
Since the total photon energy at each pair-production event is roughly conserved, cascades in low-and mid-field pulsars effectively convert high-energy photons into photons with energy ǫ min . At the penultimate generation of an extended cascade, however, a significant fraction of the synchrotron photons have ǫ < ǫ min and escape. From the −3/2 synchrotron power law and using the expected multiplicity power-law dependence on energy, equation (45), we expect a final multiplicity of with ν as given in equation (46). For pulsars with moderate and low fields, B 3 × 10 12 , the quantity a is large, a negligible fraction of the energy of each absorbed photon is left in the generated pair, and ν ≈ 1, yielding a simple linear relation between input photon energy and final pair multiplicity.
In the opposite regime, very high fields, a is small, so the generated pair retains almost all of the initial photon energy, and there is little or no cascade; high-energy photons simply transmute into high-γ particles.
In Figure 7, we plot the numerically calculated energy dependence of the total multiplicity along with the theoretical results. The simplest multiplicity formula, equation (48) with the exponent ν set to 1, works well for magnetic fields below 10 12 Gauss, but over-estimates the multiplicity for higher fields, due to the increase in the amount of energy left in the particles in the higher-field cascades. For all magnetic fields, the spectrum has three well-defined regions. If ǫ < ǫ min , the initial photon is below threshold, and no pairs are produced. If ǫ min < ǫ < ln Λǫ min , the initial photon will pair produce, but all of the secondary synchrotron photons have energies below ǫ min and create no further pairs. Finally, for ǫ > ln Λǫ min , both the initial photon and the secondary synchrotron photons may create pairs, producing a smoothly increasing multiplicity. Figure 8 shows this change in the power-law exponent with magnetic field, along with the theoretical result, equation (46). The approximation closely matches the numerical results over the entire dynamic range.
From equation (47), we expect that the spatial development of the cascade should be a similar power-law in space. To find a closer approximation, we use an argument akin to that for the final multiplicity.
Since the first several generations of a high-energy cascade occur in the high-energy regime where equation (15) applies, the cascade develops as if ǫ a were the minimum energy, only deviating when the photon energy drops near ǫ min . Hence, a good approximation to the high-energy cascade is M (ǫ, s) = 1 + 1 √ ln Λ ǫ ǫ a 4(s − s 1 ) r e ν , s 1 < s < s 1 + 27 64 where s is the propagation distance from the emission point, s 1 is the point at which the initial photon pair produces, s 1 = 0.25(ǫ a /ǫ)r e , ν is given in equation (46)   The topmost solid line shows the theoretical prediction for B = 10 11 Gauss, while the other lines represent the numerical results for different magnetic field strengths: 1 × 10 11 , 1 × 10 12 , 2 × 10 12 , 3 × 10 12 , 5 × 10 12 , and 1 × 10 13 Gauss, respectively. The decrease in the power-law slope of the response with increasing magnetic field is clearly visible. the factor of 27/64 follows from the truncation of the cascade at ǫ min rather than ǫ a .
Examining the numerical results, shown in Figure 9, we found that, indeed, the pair production rate follows this rule very well. The total length of the cascade is approximately 0.11 R * , which is the pair production distance expected for a photon of the minimum energy capable of pair-producing, 64ǫ a /27, according to the ψ a model. The computed power law, shown in Figure 10 follows the predicted values almost exactly. Figure 11 shows the variation in the output pair spectrum, for different values of the input energy. The pair spectra are characterized by a sharp rise near ǫ a /a, followed by a ν = −3/2 power law up to ǫ 0 / ln Λ, where ǫ 0 is the energy of the injected photon, followed by an exponential decline. The spike corresponding to the input photon is also visible. Figure 12 shows the variation in the output pair spectrum, for a fixed ratio of input photon energy to ǫ a , ǫ 0 /ǫ a = 10 4 , for different values of the magnetic field. Here, we notice the similarity of the pair spectra at different field strengths. The lower limit of the pair spectra is at ǫ a /(1 + a 2 ) 1/2 , which is independent of B, provided that a ≫ 1. The upper limit is a fixed multiple of this value, and so remains constant as well.

Pair spectra
6. semi-numerical model The analytic models of Paper I did not take into account the variation of the power-law exponent, ν, with magnetic field. This does not affect the calculation of either the curvature or the resonant ICS cascades, since those mechanisms operate primarily on the primary photons, but the decline of ν with B limits the effectiveness of non-resonant ICS for magnetic fields above 10 12 Gauss.
Allowing for the change of the power law makes the analytic form impractical, but we can quantify the approximation by using a simple numerical model. If we still assume that all photons emitted by a given emission mechanism are emitted at the peak energy and use equation (49 for the spatial range of pair-production, we can calculate the location of the PFF.
In practice, we find that this method yields results which are nearly indistinguishable from the full numerical method, in the regime where the full numerical method predicts a finite PFF height. In Figure 13, we show the calculated PFF heights from the analytic model and in Figure 14 those from the semi-numerical model.

full cascade
The same numerical procedure used to calculate the expected response to one photon may be applied to a full model of a pulsar polar cap. We will ignore any variation across the polar cap, concentrating on a typical field line instead.
Given the pair production physics discussed above, a polar cap model only requires knowledge of the electric potential accelerating particles off the surface of the star, the relevant emission mechanisms which inject high-energy photons into the magnetosphere, and a model for the temperature of the star.
As a model of the accelerating potential, we use the cubic approximation to the Muslimov & Tsygan (1992) GR acceleration derived in Paper I. This model is easily computable and sufficiently accurate for our purposes.
The emission mechanisms included are curvature radiation and inverse Compton scattering, both resonant and non-resonant. The spectrum of curvature radiation is well-  -Final pair spectra produced by injecting a single photon at the stellar surface, evaluated at fixed voltage, Φcap = 10 14 V, and fixed input energy ratio, ǫ = 10 4 ǫa(B), for various values of the magnetic field. The magnetic fields used are 1 × 10 11 , 3 × 10 11 , 1 × 10 12 , 3 × 10 12 , and 1 × 10 13 Gauss, from highest peak amplitude to lowest. known, c.f. Jackson (1975), while we use equation (A17) for non-resonant scattering, with spatial dilution terms included, and equation (A2) derived from Dermer (1990) for resonant scattering. See Appendix A for descriptions of the precise emissivities used.
The stellar temperature is the most uncertain quantity in these detailed models. Not only is the precise cooling rate for neutron stars still a topic of debate, but the polar cap is likely heated by beam particles reversed and accelerated back down onto the cap. To minimize the uncertainties, we simply assume a fixed stellar temperature of 10 6 K.
From the emission mechanisms, we find the Lorentz factor of the beam as a function of altitude, and may then inject the appropriate emission spectrum at each radial bin in our grid. To ensure energy conservation, we slightly adjust the peak amplitude of each injected spectrum to compensate for the effects of the discrete energy grid.
We run the model multiple times. First, given an assumed polar cap temperature, we run it as described above to find the pair formation front altitude, namely the altitude at which κ g pairs have been produced per primary. At this point the pair plasma is sufficiently dense to short out the accelerating electric field.
We then recompute the Lorentz factor of the beam, with acceleration halting at the PFF, inject the new spectra, and follow the cascade to produce the final output γ-ray and pair spectra.
A typical output pair spectrum is shown in Figure 15, and a typical output γ-ray spectrum in Figure 16. These show many of the characteristics of the single-photon response. The pair spectrum shows the rapid rise, ν = −3/2 power-law, and exponential tail, but it also shows a lowamplitude tail of high-energy particles, due to the effects of inverse Compton scattering.
For most pulsars, the final pair and γ-ray spectra are similar to those in Figures (15) and (16).
In Figure 17, we plot the multiplicity as a function of pulsar period and cap potential for a fixed stellar temperature of 10 6 K. Given that fixed temperature, most pulsars have multiplicities ranging from 0.1-100, while most of the millisecond pulsars have lower multiplicities. In Figure 18, we set the radius of curvature of the field lines equal to the stellar radius, in order to model roughly the effects of an offset dipole. Here, most pulsars have multiplicities in the range of 10-1000, while the MSPs remain with significantly lower multiplicities. However, a parallel offset of the dipole moment would also increase the surface field, effectively increasing the cap potential. Given the sharp gradient in multiplicities at low-P , this could easily bring the multiplicities of the MSPs into parity with those of the other pulsars.

conclusions
This paper has given several approximations and descriptions of the pair creation process in pulsars, in the hopes that they will be useful for other researchers examining these objects. The development of the pair cascade in space is deceptively simple, a fact which has been used in Paper I to obtain several useful results about the regimes of dominance of the various emission mechanisms. A more detailed model, maintaining the variability of the power-law exponent, produces qualitatively identical results. Many pulsars have their PFF set by non-resonant ICS, with curvature dominating at high potential, and resonant ICS at high magnetic field.
In addition, due to the effects of NRICS, many pulsars seem to operate with comparatively low pair multiplicities. The total κ is in the range of 1-10 more often than it is in the 1000's predicted by curvature models.
The OTS approximation to the pair cascade describes the spectrum of γ-rays and particles produced by the pair production process very well, although the OTS pair spectrum underestimates the number of low energy pairs, due to the decline of the magnetic field with altitude.
The full cascade model may be used to predict the γray and pair spectra of individual pulsars. Although the pair spectrum produced will be modified by RICS, the high-energy γ-ray spectrum predicted should be observable. The γ-ray output from pulsars will be the subject of a subsequent paper.
We wish to acknowledge the assistance of Alice Harding, who provided the data for the comparisons of Figures 4 and 5.

emission mechanisms
For convenience, we describe the forms of the various emission mechanisms used in the numerical calculations.

Resonant ICS
For resonant ICS, we use the result from Dermer (1990), where H is the top hat function, defined before.

Non-resonant ICS
Although the curvature emission and resonant inverse-Compton scattering spectra are well-represented in the literature, the non-resonant ICS spectrum is not, to the accuracy we desire.
The chief difference between the physical situations considered by Blumenthal & Gould (1970) and conditions in the pulsar magnetosphere is that the background X-rays in pulsars are not isotropic. They originate from the hot surface of the star itself and only occupy a range in θ from θ = 0 to at most θ = π. In terms of µ = cos θ, thermal photons are emitted with angles from 1 to µ c , where where a is the radius of the polar cap, typically a = θ c R * , and z is the height above the stellar surface. For simplicity, we only consider emission along the z-axis, where the incoming photons are symmetric in φ. The range of µ populated with photons is then For the isotropic case, ∆µ = 2. If ∆µ < 2, the simplest consequence is that there are fewer photons, by a factor of ∆µ/2, potentially scattering off of the particle. Furthermore, these photons are concentrated near µ = 1, so the mean energy of the photons in the particle rest frame is reduced.
From (Blumenthal & Gould 1970), the full non-resonant inverse Compton cross section is The primed frame is the particle rest frame, the unprimed frame is fixed to the neutron star itself, and the subscript 1 indicates the scattered particles and angles. We have made the approximation that all of the incoming photons in the particle rest frame are beamed down the z-axis, so that the angle between the incoming and the scattered photon is simply equal to the latitude of the scattered photon, Θ ′ = π − θ ′ 1 . This approximation works very well for large Lorentz factors, but is less accurate in the Thomson regime. However, it gives negligible error in the Klein-Nishina regime.
The rest-frame scattered number emissivity is then Using the δ-function in the cross section to perform the ǫ ′ integral yields The integral over the source intensity is Changing variables to ǫ, using ǫ = γǫ ′ 0 (1 + βµ ′ ), this becomes The limits of the ǫ integration are set by which energies may be scattered such that they have energy ǫ ′ 0 in the rest frame. The minimum possible energy is ǫ min = ǫ ′ 0 /γ(1 − βµ c ) ≈ ǫ ′ 0 /γ∆µ, while the maximum is ǫ max = ǫ ′ 0 /γ(1 − β) ≈ 2γǫ ′ 0 , where the approximations assume that β ≈ 1.