A FAN BEAM MODEL FOR RADIO PULSARS. I. OBSERVATIONAL EVIDENCE

Rather than developing a purely geometric beam model, we are interested in exploring how the emission mechanisms and propagation effects, secondary plasma flow, and magnetic field geometry shape the emission beam. The central question is what the beam geometry will be when the secondary relativistic particles produce broadband and coherent emission as they flow along a dipolar magnetic flux tube. Since there is no consensus on pulsar emission mechanisms and propagation effects, we prefer to use a phenomenological parameter q to describe the total effect (this will be defined in the assumption (3) below), which is adjustable according to specific emission and propagation mechanisms. The following basic assumptions are made to simplify the problem.

• 1.  
  The global magnetic field is dipolar.
• 2.  
  Assumptions on the particle flow: (1) A large amount of secondary electrons and positrons are generated in polar gaps (IVG or SCLF gap) with the multiplicity factor M ≫ 1, i.e., the total number density of electrons and positrons is approximately M times the local Goldreich–Julian (Goldreich & Julian 1969) number density n_gj. (2) Being quasi-neutral, the plasma flows freely along the magnetic flux tubes.
• 3.  
  Assumptions on emission mechanisms and propagation effects: (1) The emission is coherent in an elementary volume λ³ (λ is the wavelength of emission), where the emission power is proportional to the square of the number of particles therein. The emissions from separated elementary volumes are incoherent. (2) The emission from λ³ is broadband. The shape of the averaged radiation spectrum of a single particle is the same everywhere in the emission region. (3) The averaged emission power of a single particle is a power-law function of the emission altitude, viz. p_e∝r^q, where the altitude r is measured from the stellar center in this paper.

Some of the assumptions need to be defined more clearly. Firstly, the term "broadband" means that either the emission from an elementary coherent volume is intrinsically broadband or it is a collective effect due to the assembling of many narrowband emissions, which span a wide range of frequencies. Secondly, the "averaged power" and "averaged spectrum" for a single particle means that the quantities are averaged over a timescale for obtaining the mean pulse profile; therefore, they are stable and useful to model the mean structure of radio beam.

Although the assumption regarding the emission power and spectrum is phenomenological, it can be adapted to a variety of emission mechanisms and propagation effects by adjusting the power-law index q. For example, an electron or a positron with constant kinetic energy will emit less power at higher altitudes through the ICS process, roughly in a power law with q ∼ −2.⁴ Considering the curvature radiation instead, the index will be q ∼ −0.5.⁵ When a particular absorption effect correlated with the plasma density is considered, the plasma may be more transparent to the emission at higher altitudes; therefore, the relationship becomes flatter and may even have a positive index. In the case that the real relationship is not a power law, e.g., an exponential or some other function, the index q itself needs to be adjusted as a function of the coordinates of emission location to mimic the real relationship. It is a practical treatment for the emission and propagation mechanisms to simplify our model. In this paper, q is assumed to be a constant for a single pulsar. The statistical value will be constrained with a sample of pulsars in Section 3.2.2.

We first describe the general picture of the emission beam qualitatively with the above assumptions. As the secondary particles flow out along a flux tube, the number density keeps decreasing because of the divergent nature of the dipolar field. This normally leads to decreasing emission intensity, except that if it is too large a positive index q is assumed (q > 3 according to Section 2.2). Since the emission at a higher altitude points further from the magnetic pole, the intensity will decrease with increasing beam radius when q < 3, or on the contrary when q > 3, where the beam radius ρ is defined as the angle between the emission direction and the magnetic pole (see Figure 1). In this paper, the first case is called the radial limb-darkening intensity distribution. Figure 1 shows the schematic diagram for the formation of such a sub-beam generated from a flux tube.

Zoom In Zoom Out Reset image size

The azimuthal (transverse) intensity distribution depends on further assumptions about the particle distribution across the polar cap. In the simplest case, if the secondary particles are uniformly distributed, the intensity will be constant in any circular ring around the magnetic pole. Such a model can only account for single pulse profiles. More practically, the secondary particles may be generated in a number of separated flux tubes, and each flux tube will form a wedge-shaped sub-beam, which widens with increasing radius because of the divergence of the flux tube. There is an intensity valley between a pair of sub-beams due to a lack of particles in the region between two neighboring flux tubes. Then, the whole beam looks like an axial fan with a set of wedge-shaped sub-beams. This kind of beam is called the fan beam. When only a few flux tubes are assumed to be active in emission, the fan beam will show a very patchy pattern. Therefore, this model may give a reasonable explanation for the origin of patchy beams.

The above qualitative analysis has revealed a new type of emission beam pattern strikingly different from the conal beam structure. Unlike the conal beam, which has a circular or an elliptical boundary, the fan beam has no boundary. This is because we choose the broadband assumption, while the conal beam model chooses the narrowband assumptions. In the following subsections, we will derive the radial intensity–radius relationship, the transverse intensity–azimuth relationship, and an important prediction—the relationship between pulse width W and impact angle β. The beam structure and resultant profiles at various viewing geometries are simulated.

In this paper, the intensity is defined as the emission power over a unit solid angle around an emission direction. Since the emission from a secondary relativistic particle is beamed into a very narrow solid angle with a half opening angle of ∼1/γ, where the Lorentz factor γ is typically larger than 100, we assume that the emission direction is aligned with the tangent of the field line for the convenience of calculation. In order to derive the intensity, we first select an arbitrary subregion in a flux tube and calculate its solid angle formed by the tangents of the boundary field lines. Then, we calculate the volume, the particle number density, the total emission power, and finally, the emission intensity from the region.

Regarding the axial symmetry of magnetic dipole fields, it is convenient to describe an emission point with the polar angle θ, the azimuth angle φ, and the altitude r (counted from the stellar center) in the spherical coordinate system where the polar axis is the magnetic pole (Figures 1 and 2). The open field lines can be distinguished by two parameters: the azimuth φ between a field line plane and the meridian plane containing the rotation and magnetic axes, which is counted anticlockwise from the equatorial side of meridian plane, and the parameter f ≡ R_e/R_c, where R_e is the maximum altitude for an open field line and R_c = Pc/(2π) is the light cylinder radius, with P as the pulsar period. The parameter f indicates how close the field line is with respect to the magnetic pole. For the last open field lines, R_e varies from R_c to ∼3R_c depending on inclination and azimuth angles. For simplicity, we neglect this difference and assume R_e = R_c for all the last open field lines. Therefore, we have f = 1 for last open field lines and f > 1 for inner open field lines.

Zoom In Zoom Out Reset image size

Figure 2(a) shows a typical subregion in a slice of a flux tube. The slice confined by the field lines "1," "2," "3," and "4" is divided into a number of subregions of this kind, starting from the polar cap surface and extending to high altitudes. The field lines marked with "1" and "2" ("3" and "4") stand for the outer and inner (with respect to the magnetic pole) boundary field lines of the slice. Each pair of boundary lines has the same azimuth angle but different f parameters, where f₂ > f₁ and f₄ > f₃, meanwhile, f₁ ≃ f₃ and f₂ ≃ f₄. The whole flux tube can be divided azimuthally into a number of slices. The slices may have different boundary parameters f₁ and f₂, depending on the shape of the cross section of flux tube. We will consider two kinds of flux tube geometry in Section 2.4, in which one has a sector-shaped cross section (model A) and the other has a circular shape (model B). In model A, all the slices have the same pair of f₁ and f₂, while in model B, f₁ and f₂ depends on the azimuth φ of each subregion. The cross sections of flux tube and the corresponding slices are illustrated by Figure 2(b) for these two models.

For clarity, the poloidal and toroidal cross sections of the subregion are presented by Figures 2(c) and (d), respectively. In panel (c), the two straight lines with the polar angles θ and θ + dθ stand for the lower and upper boundaries (with respect to the polar cap) of the subregion. At a pair of points where a boundary field line intersect the straight lines, e.g., A and B, or C and D, the tangents form a small angle dρ, where the beam radius ρ is the angle between the tangent of a field line at r and the magnetic pole. In panel (d), the two azimuthal boundary field lines, marked with "1" and "3," are separated from each other by a small azimuth angle dφ. The tangents of field lines at two boundary points A and E open an angle dτ. The dρ and dτ should be much larger than the 1/γ beaming angle of a single particle so that the angular power distribution in the single-particle emission beam can be ignored. Whenever γ is more than a few hundreds, dρ ∼ 1° and dτ ∼ 1° would be large enough to ensure that. In this circumstance, the particle number density varies little poloidally between A and B (or C and D) as well as toroidally between A and E, enabling the following derivation to be valid.

The tangents at eight vertices of the subregion form a solid angle dΩ, as projected in the celestial sphere centered on the star (Figure 3). It is easy to find dΩ = dρdτ, where dτ ≃ sin ρdφ, which can be derived by using the law of cosines in the spherical triangle ▵MST. In the dipole field geometry, one has ρ ≃ 3/2θ when θ is not so large (see Appendix A). Then, the solid angle reads

$$\begin{equation} {d}\Omega \simeq \frac{9}{4}\theta {d}\theta {d}\varphi. \end{equation} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

When deriving the emission power dP for this subregion, it should be noted that the arc length ds of the inner boundary line would be significantly larger than that of the outer boundary line when f₂ ≫ f₁, and the particle number density would vary remarkably with f as well. Therefore, one needs to divide the whole subregion within f₁ and f₂ into a number of smaller ones with an interval df. The volume of such an elementary region is (see Appendix B)

$$\begin{equation} {d}V_{\rm f}\simeq R_{\rm c}^3\theta ^7f^2{d}f{d}\theta {d}\varphi. \end{equation} \tag{ 2 }$$

The next step is to find the particle number density in dV_f and calculate the emission power dP_f, where the subscript "f" indicates that these quantities are for the elementary region within f and f + df. Then, the total emission power dP is the integration of dP_f over f.

The altitude dependence of particle number density can be derived using the dipole field geometry and the assumption of free flow. According to the assumption, the initial total number density of secondary electrons and positrons is

$$\begin{equation} n_0=Mn_{\rm gj0}, \end{equation} \tag{ 3 }$$

with the Goldreich–Julian number density

$$\begin{equation} n_{\rm gj0}\simeq |\mathbf {\Omega }\cdot \mathbf {B_{\rm s}}|/(2\pi ce)\simeq B_{\rm s}\cos \alpha /Pce, \end{equation} \tag{ 4 }$$

where Ω is the rotation angular velocity, B_s is the surface magnetic field, α is the inclination angle, and c and e are the light velocity and the electron charge, respectively. In the case of free flow, when neglecting the trivial energy loss due to radiation, the particle number density follows the flux conservation law nS = n₀S₀, where S and S₀ are the cross section areas of a flux tube at altitudes r and R (stellar radius), respectively; S and S₀ are related to each other by the law of magnetic flux conservation, i.e., BS = B₀S₀. Combining these relations with the dipole magnetic field strength B ≃ B_s(R/r)³, one has

$$\begin{equation} n\simeq n_{\rm 0}(R/r)^3. \end{equation} \tag{ 5 }$$

Under the basic assumptions of emission mechanisms and propagation effects, the emission power from a coherent volume λ³ is

$$\begin{equation} \delta P_{\rm c}= (n\lambda ^3)^2i_{\rm e0}(r/R)^{ q}, \end{equation} \tag{ 6 }$$

where i_e0, the averaged emission power of a single particle, is modified by a power-law term (r/R)^q to represent the possible dependence on altitude. Note that i_e0 may be a function of emission frequency, but in the following calculation we ignore such a dependence. It is possible influence on the frequency dependence of pulse profiles will be studied in a subsequent paper.

In a larger volume dV_f that consists of many coherent volumes, the emissions are incoherently superposed, thus the emission power is

$$\begin{equation} {d}P_{\rm f} = \frac{{d}V_{\rm f}}{\lambda ^3}\delta P_{\rm c}={d}V_{\rm f} n^2{\lambda ^3}i_{\rm e0}(r/R)^{q}. \end{equation} \tag{ 7 }$$

Integrating over all the elementary regions from f₁ to f₂, one has the total emission power from the subregion

$$\begin{equation} {d}P = \int _{f_1}^{f_2}{d}P_{\rm f}. \end{equation} \tag{ 8 }$$

Eventually, the emission intensity reads

$$\begin{equation} I = \frac{{d}P}{{d}\Omega } . \end{equation} \tag{ 9 }$$

Substituting Equations (1)–(8) into Equation (9), and using the dipole field relations $B_{\rm s}\simeq 3.2\times 10^{19}(P\dot{P})^{1/2}({\rm Gauss})$ and r ≃ 4/9R_eρ² (for r ≪ R_e), we have the emission intensity (see Appendix B for details)

$$\begin{equation} I=AP^{q-4}\dot{P}\cos ^2\alpha \left(f_{1}^{q-3}-f_{2}^{q-3}\right)\rho ^{2q-6}, \end{equation} \tag{ 10 }$$

where the coefficient is

$$\begin{equation} A={1\over {3-q}}\left({{6.4\times 10^{19}} \over {3ce}}\right)^2 \left({{2c}\over {9\pi }}\right)^{q-3} R^{6-q} \lambda ^3 M^2 i_{\rm e0}. \end{equation} \tag{ 11 }$$

The last term on the right-hand side, ρ^{2q − 6}, indicates a radius dependence of emission intensity resulting from the assumptions of coherent emission, plasma free flow, and index q, which have been qualitatively discussed before. When q < 3, it predicts a limb-darkening relationship.

Beside this, the term $(f_{1}^{q-3}-f_{2}^{q-3})$ also affects the I − ρ relationship. It may cause different trends in the central and outer parts of emission beam. To see this, without losing generality, we consider a flux tube located between the last open field lines (f_om = 1) and a layer of the innermost open field lines with a fixed f_im, where the initial particle density n₀ is constant on the polar cap and the coherent emission starts from the polar cap (relaxation of these assumptions will be discussed soon). In order to see the emission from all the field lines between f_om and f_im, there must be ρ ⩾ ρ_t, where the transition radius ρ_t equals the opening angle of the polar cap boundary ρ_pc ≃ (3/2)(R/R_c)^1/2. In this case, namely for the outer beam with ρ ⩾ ρ_t, the term $(f_{\rm om}^{q-3}-f_{\rm im}^{q-3})$ is a constant, and the intensity follows I∝ρ^{2q − 6}. However, in the inner beam with ρ ⩽ ρ_t, for any emission direction with ρ, the outermost visible field line is determined where the emission direction from its foot point on the polar cap is aligned with the LOS, which leads to f₁ = (9/4)(R/R_c)ρ⁻². Therefore, the visible emission region is confined between f₁ and f₂ = f_im, which obviously shrinks, because f₁ > 1. Then, Equation (10) reduces to

$$\begin{eqnarray} I_{\rm inner} &=& AP^{-1}\dot{P}\cos ^2\alpha \nonumber \\ &&\times \left[1-\left(\frac{9\pi R}{2c}\right)^{3-q} f_{\rm im}^{q-3}P^{q-3}\rho ^{2q-6}\right]. \end{eqnarray} \tag{ 12 }$$

When q < 3, it shows an opposite trend in which the intensity increases with increasing radius in the inner beam.

In the above analysis, the polar cap boundary plays the role of the transition location for intensity distributions, because we assume a uniform n₀, set the lowest coherent emission altitude to be R, and select the last open field lines as the outer boundary of flux tube. Any deviation from these three assumptions will lead to a different transition location. For example, if a layer of inner open field lines is selected as the outer boundary of a flux tube, namely f_om > 1, or n₀ is peaked at an inner open field line within the flux tube, the transition location will move inward to the magnetic pole, thereby ρ_t < ρ_pc. The latter case can be seen in Figures 4 and 9 (models B1 and B2). On the contrary, if the lowest coherent emission altitude is high above the polar cap, i.e., r_L > R, while the other two assumptions remain, the transition location will have $\rho _{\rm t}\simeq \sqrt{r_{\rm L}/R}\rho _{\rm pc}>\rho _{\rm pc}$.

Zoom In Zoom Out Reset image size

To summarize, the radial intensity distribution in the sub-beam formed by a flux tube is twofold: starting from the magnetic pole (or a place near the magnetic pole), the intensity increases first, reaches its maximum at ρ = ρ_t, and then fades with increasing ρ in a form of power law. This feature can be seen in the top panels of Figures 5–12.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The transverse (azimuthal) intensity distribution in the beam depends on the azimuthal distribution of number density of secondary particles. Being enlightened by the idea of sparks in RS75 and GS00, it is assumed that there are some discharging flux tubes on the polar cap, but their transverse sizes are not necessary to be the same as those of sparks. The number density of secondary particles, or equivalently the multiplicity factor M, is probably a function of f and φ within a flux tube. Then, Equations (10)–(12), which only apply to the homogeneous distribution of M, must be revised. For simplicity, the lowest coherent emission height is set as r_L = R in the following derivation. The effect of a higher r_L will be discussed by the end of this section.

We start from an elementary region within f − f + df at a given φ in a flux tube. Its intensity should be dI∝(3 − q)M(f, φ)²f^{q − 4}df (by replacing f₁ and f₂ with f and f + df in Equation (10), respectively). Then, for the outer part of the sub-beam, namely, $\rho >(3/2)\sqrt{R/(R_{\rm c}f_1)}$, the total intensity from the whole flux tube reads

$$\begin{eqnarray} I_{\rm outer}(\varphi)&=&A_1P^{q-4}\dot{P}\cos ^2\alpha \rho ^{2q-6}\nonumber\\ &&\times \int _{f_1}^{f_2}M(f,\varphi)^2f^{q-4}{d}f, \end{eqnarray} \tag{ 13 }$$

where the coefficient is

$$\begin{equation} A_1=\left({{6.4\times 10^{19}} \over {3ce}}\right)^2 \left({{2c}\over {9\pi }}\right)^{q-3} R^{6-q}\lambda ^3 i_{\rm e0}. \end{equation} \tag{ 14 }$$

Note that the boundaries f₁ and f₂ may also be a function of φ, depending on the shape of the cross section of a flux tube on the polar cap.

For the inner part where $\rho \lesssim (3/2)\sqrt{R/(R_{\rm c}f_1)}$, the outmost visible field line corresponding to a given ρ shrinks inward, which satisfies f_o = 9R/(4R_cρ²). Thus, the lower boundary of integration should be determined by $f_1^{\prime }=\min [f_{\rm o}, f_1]$, where f₁ is the real outer boundary of the flux tube at a given azimuth φ. Then the total intensity is

$$\begin{eqnarray} I_{\rm inner}(\varphi)&=& A_1P^{q-4}\dot{P}\cos ^2\alpha \rho ^{2q-6}\nonumber\\ &&\times \int _{f_1^{\prime }}^{f_2}M(f,\varphi)^2f^{q-4}{d}f. \end{eqnarray} \tag{ 15 }$$

To explore how the geometry of flux tubes and the secondary particle distribution affect the beam geometry, we consider the two kinds of models described below.

In the first case (model A), M is assumed to be constant along the colatitude dimension but follows a Gaussian distribution in the azimuthal dimension, peaked at the central azimuth φ₀ of the flux tube, i.e., M = M₀exp [ − (φ − φ₀)²/(2σ²)]. The flux tube is assumed to be within [f₁, f₂] for all the azimuth angles, where f₁ and f₂ are constant. For the intensity in the outer part, Equation (13) is modified as

$$\begin{eqnarray} I_{\rm outer}^{\rm A}(\varphi)&=&A^\prime P^{q-4}\dot{P}\cos ^2\alpha {\rm e}^{-(\varphi -\varphi _0)^2/\sigma ^2} \nonumber \\ &&\times \left(f_{1}^{q-3}-f_{2}^{q-3}\right)\rho ^{2q-6}, \end{eqnarray} \tag{ 16 }$$

where the coefficient is

$$\begin{equation} A^\prime ={1\over {3-q}}\left({{6.4\times 10^{19}} \over {3ce}}\right)^2 \left({{2c}\over {9\pi }}\right)^{q-3} R^{6-q} \lambda ^3 M_0^2 i_{\rm e0}. \end{equation} \tag{ 17 }$$

For the intensity in the inner part, Equation (15) is modified as

$$\begin{eqnarray} I_{\rm inner}^{\rm A}(\varphi)&=&A^\prime P^{-1}\dot{P}\cos ^2\alpha {\rm e}^{-(\varphi -\varphi _0)^2/\sigma ^2} \nonumber \\ &&\times \left[1-\left(\frac{9\pi R}{2c}\right)^{3-q}f_{\rm im}^{q-3}P^{q-3}\rho ^{2q-6}\right]. \end{eqnarray} \tag{ 18 }$$

In the second case (model B), it is assumed that the cross section of a flux tube is circular and M follows a Gaussian distribution around the center of the flux tube. For convenience, a dimensionless colatitude χ ≡ θ/θ_pc is used below to describe the latitudinal position of the foot point of a field line on the polar cap, where θ and θ_pc are the polar angles of the foot points of an inner open field line and the last open field lines, respectively. The χ and f can be converted to each other by

$$\begin{equation} \chi ={3\over 2}\left(\frac{R}{fR_{\rm c}}\right)^{1/2}\rho _{\rm pc}^{-1}. \end{equation} \tag{ 19 }$$

Given the position of the center (χ₀, φ₀) and the dimensionless radius ℜ₀ of the flux tube cross section (ℜ₀ ⩽ 1 − χ₀), the Gaussian distribution of M is M = M₀exp [ − ℜ²/(2σ²)], where ℜ is the dimensionless angular distance from an arbitrary point to the flux tube center. Then the intensity of the outer part reads

$$\begin{eqnarray} I_{\rm outer}^{\rm B}(\varphi)&=&A_2P^{q-4}\dot{P}\cos ^2\alpha \rho ^{2q-6} \nonumber \\ &&\times \int _{f_1}^{f_2}{\rm e}^{-\frac{\chi _0^2+\chi ^2-2\chi _0\chi \cos (\varphi -\varphi _0)}{\sigma ^2}}f^{q-4}{d}f, \end{eqnarray} \tag{ 20 }$$

where

$$\begin{equation} A_2=\left({{6.4\times 10^{19}} \over {3ce}}\right)^2 \left({{2c}\over {9\pi }}\right)^{q-3} R^{6-q}\lambda ^3 i_{\rm e0} M_0^2, \end{equation} \tag{ 21 }$$

with the boundary

$$\begin{eqnarray*} &&f_{1,2}=(9/4)(R/R_{\rm c})[\rho _{\rm pc}\Re _0\sin \vartheta _{1,2}/\sin (\varphi -\varphi _0)]^{-2}, \end{eqnarray*}$$

where $\vartheta _{1,2}=\arccos [\sin (\varphi -\varphi _0)^2 \chi _0/\Re _0\mp \cos (\varphi -\varphi _0)\sqrt{1-(\chi _0/\Re _0)^2\sin (\varphi -\varphi _0)^2}]$ (see Appendix C for derivation). In the inner part of the sub-beam, there is

$$\begin{eqnarray} I_{\rm inner}^{\rm B}(\varphi)&=& A_2P^{q-4}\dot{P}\cos ^2\alpha \rho ^{2q-6} \nonumber \\ &&\times \int _{f_1^\prime }^{f_2}{\rm e}^{-\frac{\chi _0^2+\chi ^2-2\chi _0\chi \cos (\varphi -\varphi _0)}{\sigma ^2}}f^{q-4}{d}f. \end{eqnarray} \tag{ 22 }$$

In fact, Equations (16)–(18) and (20)–(22) give the full description for the radial and transverse intensity distribution of a sub-beam in models A and B, respectively. The beam shape and resultant profiles will be modeled in the following subsection.

In the above derivation it is assumed that the coherent emission starts from the polar cap. If the lowest emission altitude is higher, the transition radius for the outer and inner parts of beam should be $\rho _{\rm t}=(3/2)\sqrt{r_{\rm L}/(R_{\rm c}f_1)}$.

The above formulae only apply to a single flux tube; to form a global beam pattern, one needs additional assumptions of the configuration of discharging flux tubes across the polar cap, including the locations of flux tubes and the multiplicity factors therein (may be inhomogeneous). The following simulations for models A and B are presented as simple examples. There is a vast variety of configurations of discharging areas and particle density distributions, whose observational consequences differ in details. However, these simple examples are still useful to infer the necessary configurations that can account for observations.

The simulations follow two steps. In the first step, we designate the geometry of flux tubes. In model A, 8 discharging areas with sector-shaped cross sections are assumed to exist on the polar cap, each occupying an azimuthal range of 45° (see Figure 4). Their central azimuths are located every 45° starting from φ = 0. In model B, the discharging areas are assumed to be circular, which are all centered at χ₀ = 0.7 with a radius of ℜ₀ = 0.3. The centers are placed every 515 in the azimuth dimension starting from φ = −5°, so there are a total of seven discharge areas (see Figure 4).

In the second step, we assign the multiplicity distribution across the polar cap. It is of particular interest to see how the inhomogeneity of the multiplicity among flux tubes affects the beam pattern and average pulse profiles. Therefore, we investigate two subclasses for each model: one is the homogeneous case where all the discharging areas have the same multiplicity pattern, which are called models A1 and B1, respectively, and the other is the inhomogeneous case where two areas are dominant in pair production over the others, which are called models A2 and B2. The dominant areas are centered at φ = 0 and 180° in model A2 and centered at φ = −5° and 149° in model B2. In both of these sub-models, the maximum multiplicity M₀ in the dominant areas (at the central azimuth in model A2 and at the center of the area in model B2) are 10 times of those in the other areas.

In the following simulation the index q is fixed as 2, close to the values constrained from a sample of pulsars in Section 3.2. The lowest emission altitude is set as r_L = R for simplicity.

Figures 5–12 show the beam patterns of models A1, A2, B1, and B2 for inclination angles of α = 10° and 80°, respectively, together with the average pulse profiles at a number of impact angles for each model. In the contour maps of intensity distribution (the top left panels in each figure), we plot the equi-radius circles from 1 to 6 times of ρ_pc. This is especially helpful for viewing the intensity structure in the inner beam (the top right panels). A set of LOSs are also plotted in the contour maps, marked with their impact angles. We take ρ_pc ≃ 33 in our calculation, corresponding to a pulsar period of P = 0.15 s.

The above models assume a few flux tubes in the polar cap. To compare with the case of a large number of flux tubes, we simulate a model B3, where the single flux tube follows the pattern of model B, but a total of 90 flux tubes are evenly spaced along the polar cap boundary and the size of each flux tube is very small. We use the following parameters and assumptions in the simulation: χ₀ = 0.9, ℜ₀ = 0.05, and the same peak multiplicity factors for all the flux tubes. The diagram of the flux tubes and the simulated beam and profiles for α = 30° are presented in Figures 13 and 14.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Below we summarize the main features and the inference that can be drawn from the simulated results.

• 1.  
  Both kinds of models present a limb-darkening feature in both the radial and the azimuthal dimensions, except that the intensity trend undergoes transition near or within the polar cap boundary (mainly because of the assumption r_L = R). The radial limb-darkening feature in the outer beam follows a power law of ρ^{2q − 6}, which is a consequence of the cooperation of coherent emission, particle free flow, and altitude dependence of single-particle emission power. The transition of the radial intensity trend in the inner beam is caused by the shrinkage of the emission region for very small ρ and the attenuation of secondary particle density toward the magnetic pole (if it exists, e.g., in model B). The transverse limb-darkening phenomenon is caused by the continuous attenuation of particle density toward the edge of the flux tube.
• 2.  
  Discarding the detailed intensity distribution in the sub-beams, the shapes of sub-beams in models A and B look very similar. This is because the shape is determined by the divergence nature of dipolar flux tubes.
• 3.  
  In the outer beam, the pulse width broadens with increasing |β|. This is also a natural consequence of the divergence nature of dipolar flux tubes. However, this trend does not hold when the LOS sweeps across the central part of the beam, because the intensity distribution therein is more complex and the LOS may sweep across the bright parts of a few sub-beams. This break of the pulse width–impact angle relationship is indeed seen in observational data, as shown in Section 3.2.
• 4.  
  Models A1 and B1 tend to predict more complex and wider profiles than models A2 and B2. This is because all the flux tubes have the same activity in emission in models A1 and B1, thus the emission from the flux tubes outside the meridian plane are still strong enough to be observed, even though the radial limb-darkening effect has caused more attenuation in the observed intensity for them. In particular, models A1 and B1 predict a growing complexity with increasing |β| in the cases of small inclination angles, which is not supported by the data (see Table 1 for α and β and the number of pulse components N_c). Since normal pulsars have simple profiles (e.g., Karastergiou & Johnston 2007), these features strongly suggest that only a limited number of discharging flux tubes should be dominant in pair production and emission activity.
• 5.  
  Ignoring the spikes, model B3 can only produce single-component profiles for relatively large impact angles. In order to predict double- and multiple-component profiles, one has to assume that the peak multiplicity factors in the flux tubes are modulated azimuthally and form a few large-scale structures. General speaking, this kind of model with small-size flux tubes and large-scale multiplicity modulation is equivalent to the models invoking large-size flux tubes in the capability of predicting various kinds of pulse profiles.

Table 1. Parameters of Pulsars

PSR J Name	α	β	Freq.	W10	$W^1_{10}$	$W^2_{10}$	Freq.	Nc	P	$\lg \dot{P}$	L400	L1400	$\lg \dot{E}$	d
	(°)	(°)	(GHz)	(°)	(°)	(°)	(GHz)		(s)	(s/s)	(mJy kpc2)	(mJy kpc2)	(erg s−1)	(kpc)
0034−0721a	3 ± 1	3 ± 1	1.17 1	43.2 ± 2.3	40.9 2	45.4 2	0.4/1.4	2	0.943	−15.39	55.17	11.67	31.28	1.03
0133−6957	80 ± 90	−3.0 ± 0.5	0.66 3	13.0 ± 1.3	13.0 3		0.66	2	0.464	−15.92	29.28		31.68	2.42
0134−2937	14 ± 3	14.0 ± 0.9	0.44 3	24.0 ± 6.0	30.0 3	17.9 4	0.44/1.4	2	0.137	−16.11	28.52	7.60	33.08	1.78
0248+6021	$46^{+64}_{-21}$	$6.0^{+2.3}_{-2.5}$	2.10 5	19.0 ± 2.0		19.0 5	2.10	1	0.217	−13.26		54.80	35.32	2.00
0304+1932	162.4 ± 11.8	0.96 ± 0.63	1.42 6	17.3 ± 1.6	18.9 2	15.7 2	0.4/1.6	2	1.388	−14.89	24.37	2.71	31.28	0.95
0332+5434	51 ± 28	$-4^{+3}_{-2}$	0.41 7	25.4 ± 0.7	26.0 2	24.7 2	0.4/1.6	3	0.715	−14.69	1500.00	203.00	32.35	1.00
0528+2200	116.8 ± 4.6	−1.50 ± 0.08	1.42 6	19.8 ± 1.0	20.8 2	18.8 2	0.4/1.6	2	3.746	−13.40	296.31	46.79	31.48	2.28
0540−7125	42 ± 16	−1.2 ± 0.2	0.43 3	14.0 ± 1.4	14.0 3		0.43	4	1.286	−15.09	36.18		31.18	2.69
0627+0706m	86 ± 0.2	8.7 ± 0.3	1.40 8	7.1 ± 0.7	7.1 8		1.40	2	0.476	−13.53	234.46	130.92	34.04	7.88
0627+0706i	94 ± 0.3	0.7 ± 0.1	1.40 8	14.4 ± 1.4	14.4 8		1.40	2	0.476	−13.53	138.11	48.53	34.04	7.88
0630−2834	70 ± 50	$-12^{+9}_{-2}$	1.41 9	35.8 ± 2.4	38.2 2	33.4 2	0.4/1.6	1	1.244	−14.15	21.09	2.36	32.18	0.32
0659+1414	29 ± 23	8.9 ± 6.1	1.42 6	36.5 ± 5.8	42.2 2	30.5 2	0.4/1.6	1	0.385	−13.26	0.51	0.29	34.58	0.28
0738−4042	62	13	0.69 10	37.6 ± 4.7	42.3 10	32.9 10	0.69/3.1	3	0.375	−14.79	486.40	204.80	33.08	1.60
0814+7429b	8.7 ± 0.2	−4.7 ± 0.2	0.33 11	26.1 ± 1.1	25.0 2	27.2 2	0.4/1.6	1	1.292	−15.77	14.74	1.87	30.49	0.43
0826+2637m	98.9 ± 0.7	−3.03 ± 0.01	1.42 6	7.6 ± 0.9	8.4 2	6.6 2	0.4/1.6	1	0.531	−14.77	7.48	1.02	32.66	0.32
0826+2637i	81.1 ± 0.7	14.77 ± 0.01	1.42 6	23.0 ± 2.3		23.0 6	1.4	1	0.531	−14.77	0.07	0.01	32.66	0.32
0835−4510	60	−6.2 ± 0.6	2.30 12	17.8 ± 0.9	18.7 13	16.9 12	0.4/2.3	2	0.089	−12.90	392.00	86.24	36.84	0.28
0837+0610	50	2.6	0.33 14	10.0 ± 0.5	9.6 2	10.5 2	0.4/1.4	2	1.274	−14.17	46.14	2.07	32.11	0.72
0908−4913m	96.1 ± 0.4	−5.9 ± 0.6	1.40 15	15.0 ± 1.5		15.0 15	1.40	2	0.107	−13.82	20.14	7.19	35.69	1.00
0908−4913i	83.9 ± 0.2	6.3 ± 0.4	1.40 15	17.0 ± 1.7		17.0 15	1.40	2	0.107	−13.82	7.86	2.81	35.69	1.00
0934−5249	13 ± 3	−4.3 ± 0.3	0.66 3	10.8 ± 0.3	11.0 3	10.5 16	0.66/1.4	1	1.445	−14.33	154.53	10.30	31.78	2.93
0942−5657	90	2.7 ± 0.9	0.66 3	8.1 ± 0.9	9.0 3	7.1 16	0.66/1.4	2	0.808	−13.40	327.61	18.14	33.47	5.02
0953+0755m	105.4 ± 0.5	22.1 ± 0.1	1.42 6	30.4 ± 1.0	29.4 2	31.4 2	0.4/1.6	1	0.253	−15.64	26.72	5.61	32.75	0.26
0953+0755i	74.6 ± 0.5	52.9 ± 0.1	1.42 6	47.0 ± 4.7		47.0 6	1.42	2	0.253	−15.64	0.53	0.11	32.75	0.26
1015−5719	101 ± 5	20 ± 5	1.37 17	155 ± 16		155.0 17	1.37	3	0.140	−13.24		21.35	35.92	4.87
1036−4926	70 ± 50	2.6 ± 0.5	0.66 3	15.0 ± 1.5	15.0 3		0.66	3	0.510	−14.78	682.78		32.69	8.71
1042−5521	32 ± 9	3.4 ± 0.5	0.66 3	12.5 ± 0.5	13.0 3	12.0 16	0.66/1.4	1	1.171	−14.17	682.09	30.21	32.22	6.98
1057−5226m	75.2 ± 0.4	36.1 ± 0.6	1.37 8	37.8 ± 3.8		37.8 8	1.37	4	0.197	−14.23	107.02		34.48	1.53
1057−5226i	104.8 ± 0.4	6.1 ± 0.6	1.37 8	47.7 ± 4.8		47.7 8	1.37	3	0.197	−14.23	80.25		34.48	1.53
1136+1551	88 ± 5	7 ± 2	1.40 18	11.8 ± 1.3	13.0 2	10.5 2	0.4/1.6	2	1.188	−14.43	31.48	3.92	31.94	0.35
1141−6545c	$20^{+16}_{-8}$	2.3 ± 1.4	1.40 19	18.0 ± 1.8		18.0 19	1.4	2	0.394	−14.37		29.70	33.45	3.00
1253−5820	35 ± 10	6.6 ± 0.9	0.44 3	22.7 ± 4.4	27.0 3	18.3 16	0.44/1.4	3	0.255	−14.68	172.87	30.25	33.70	2.94
1328−4357	55 ± 25	3.0 ± 0.6	0.44 3	14.4 ± 1.6	16.0 3	12.8 20	0.44/1.6	2	0.533	−14.52	94.39	10.49	32.90	2.29
1513−5908	80	35 ± 10	1.35 21	95.2 ± 9.5		95.2 21	1.35	2	0.151	−11.82	29.04	18.20	37.23	4.40
1527−3931	70 ± 55	−1.2 ± 0.3	0.44/0.66 3	7.0 ± 0.7	7.0 3		0.66	2	2.418	−13.72	44.44		31.73	2.01
1543+0929	131.0 ± 5.7	−20.2 ± 2.3	1.42 6	109 ± 20	88.5 2	129.2 2	0.4/1.4	3	0.748	−15.36	2715.18	205.38	31.61	5.90
1549−4848m	92.5 ± 0.2	−3.5 ± 0.2	1.40 8	17.2 ± 2.8	19.9 3	14.4 8	0.66/1.4	2	0.288	−13.85	28.99	0.8	34.37	1.54
1549−4848i	87.5 ± 0.5	1.5 ± 0.2	1.40 8	13.6 ± 1.4		13.6 8	1.40	3	0.288	−13.85	11.31	0.31	34.37	1.54
1550−5418	160	14	8.36 22	117 ± 12		117.0 22	8.36	6	2.070	−10.63		313.06	35.00	9.74
1559−4438	37 ± 12	−3.1 ± 0.5	0.66/1.5 3	21.5 ± 4.5	17.0 3	26.0 3	0.66/1.5	2	0.257	−14.99	581.90	211.60	33.37	2.30
1603−5657	26 ± 6	6.3 ± 0.6	0.66 3	13.7 ± 4.3	18.0 3	9.4 16	0.66/1.4	2	0.496	−14.55	580.72	38.47	32.96	8.52
1622−3751	40	4.0	1.37 23	15.5 ± 1.6		15.5 23	1.37	2	0.731	−14.59		3.04	32.42	11.55
1651−5222	50 ± 20	−3.7 ± 0.6	0.66 3	14.3 ± 0.7	15.0 3	13.6 16	0.66/1.4	2	0.635	−14.74	939.14	118.41	32.45	6.39
1705−1906m	94 ± 2	−12 ± 2	4.85 24	17.7 ± 0.5	17.3 2	18.2 2	0.4/1.6	2	0.299	−14.38	35.13	9.65	33.79	1.18
1705−1906i	86 ± 2	−4 ± 2	4.85 24	8.1 ± 0.2	7.9 2	8.3 2	0.6/1.6	1	0.299	−14.38	5.27	1.45	33.79	1.18
1709−1640	90 ± 65	−2 ± 1	4.85 24	11.9 ± 1.0	12.8 2	10.9 2	0.4/1.6	1	0.653	−14.20	75.81	6.45	32.95	1.27
1710−2616	30	20	1.37 23	148 ± 15		148.0 23	1.37	4	0.954	−16.70		9.46	29.95	4.26
1722−3712m	90.7 ± 0.1	5.4 ± 0.3	1.40 8	14.9 ± 1.2	16.0 3	13.7 16	0.66/1.4	1	0.236	−13.96	145.46	18.63	34.51	2.51
1722−3712i	89.3 ± 0.1	6.7 ± 0.5	1.40 8	14.2 ± 1		14.2 8	1.40	2	0.236	−13.96	12.04	1.53	34.51	2.51
1739−2903m	84.2 ± 0.3	3.3 ± 0.2	1.40 8	22.7 ± 1.2	23.8 2	21.5 8	0.6/1.6	2	0.323	−14.10		14.23	33.97	3.19
1739−2903i	95.8 ± 0.4	−8.2 ± 0.4	1.40 8	16.6 ± 4.3	20.8 3	12.3 8	0.66/1.4	1	0.323	−14.10		6.12	33.97	3.19
1740+1311	90 ± 86	−2 ± 2	4.85 24	14.3 ± 1.4		14.3 2	1.4	3	0.803	−14.84	546.07	88.74	32.05	4.77
1750−3503	11 ± 2	−5.5 ± 0.5	0.66 3	56 ± 5.6	56.0 3		0.66	1	0.684	−16.42	745.44	20.31	30.67	5.07
1751−4657	90	−1.7 ± 1.2	0.43 3	11.5 ± 0.5	12.0 3	11.0 25	0.43/0.95	2	0.742	−14.89	74.26	10.61	32.10	1.03
1825−0935m	95	−7	0.69 10	25.7 ± 1.1	26.7 10	24.6 10	0.69/3.1	2	0.769	−13.28	2.95	0.98	33.66	0.30
1825−0935i	85	3	0.69 10	15.0 ± 1.0	16.0 2	14.0 2	0.6/1.6	2	0.769	−13.28	0.29	0.1	33.66	0.30
1828−1101m	97.3 ± 0.6	−10.6 ± 1.5	3.10 8	11.0 ± 1.1		11.0 8	3.10	1	0.072	−13.83		121.32	36.19	7.26
1828−1101i	82.7 ± 0.6	4.0 ± 1.4	3.10 8	10.0 ± 1.0		14.0 8	3.10	1	0.072	−13.83		31.53	36.19	7.26
1841+0912	86.1 ± 11.4	2.30 ± 0.04	1.42 6	14.1 ± 1.5	12.6 2	15.5 2	0.4/1.4	3	0.381	−14.96	124.00	10.54	32.89	2.49
1852−2610	17 ± 3	−8.1 ± 0.6	0.43 3	39.0 ± 3.9	39.0 3		0.43	2	0.336	−16.06	61.29	7.15	31.96	2.26
1900−2600	25	1.8	0.33 26	40.5 ± 1.3	41.8 2	39.2 2	0.4/1.6	5	0.612	−15.69	64.19	6.37	31.54	0.70
1906+0746m	$81^{+1}_{-66}$	8.9 ± 2.6	1.4 27	6.0 ± 0.6		6.0 28	1.4	1	0.144	−13.69	16.35	9.99	35.43	4.53
1906+0746i	$99^{+66}_{-1}$	−9.1 ± 2.6	1.4 27	21.0 ± 2.1		21.0 28	1.4	1	0.144	−13.69	2.12	1.3	35.43	4.53
J1915+1606	157 ± 5	5.0 ± 1.2	1.4 29	53.0 ± 2.0	55.0 30	51.0 29	0.43/1.4	2	0.059	−17.06	203.35	45.75	33.23	7.13
1917+1353	73 ± 19.4	5.4 ± 0.5	1.42 6	18.6 ± 1.1	19.7 2	17.5 2	0.4/1.4	2	0.195	−14.14	1075.00	47.50	34.59	5.00
1918+1444	118.0 ± 33.4	−1.00 ± 0.33	1.42 6	8.4 ± 0.8	9.1 2	7.6 2	0.4/1.6	3	1.181	−12.67	3.18	1.99	33.71	1.41
1932+1059m	35.97 ± 0.95	25.55 ± 0.87	1.42 6	19.5 ± 1.4	20.9 2	18.1 2	0.4/1.6	3	0.227	−14.94	29.12	3.46	33.59	0.31
1941+0121	138 ± 32	8 ± 4	0.82 31	60.3 ± 6.0	60.3 31		0.82	2	0.217	−15.72	14.30	3.50	32.87	2.89
2043+2740	67 ± 15	10 ± 2	1.36 32	16.2 ± 1.6		16.2 32	1.36	2	0.096	−14.90	19.15		34.75	1.13
2048−1616	34	−1.6	0.69 10	17.1 ± 1.6	18.6 2	15.5 2	0.4/1.6	3	1.962	−13.96	104.69	11.73	31.76	0.95
2053−7200	27 ± 8	1.7 ± 0.3	0.66 3	40.6 ± 1.6	39.0 3	42.2 33	0.66/1.4	2	0.341	−15.71	31.97	6.62	32.29	1.05
2113+4644	38 ± 37	−4 ± 4	4.85 24	74.3 ± 1.2	73.1 2	75.5 2	0.6/0.9	3	1.015	−15.15	3680.00	304.00	31.43	4.00
2144−3933	90	0.4 ± 0.1	0.66 3	6.0 ± 0.6	6.0 3		0.66	1	8.510	−15.30	0.41	0.02	28.50	0.16
2240+5832	$108^{+42}_{-98}$	$15^{+0.8}_{-12.3}$	1.41 5	13.4 ± 1.3		13.4 5	1.41	3	0.140	−13.81		595.41	35.34	14.85
2337+6151	90 ± 88	−8 ± 7	4.85 24	25.8 ± 0.6	25.2 2	26.3 2	0.4/1.4	2	0.495	−12.71	61.01	8.54	34.80	2.47
2346−0609	28 ± 9	−1.3 ± 0.3	0.44 3	21.0 ± 2.1	21.0 3		0.44	3	1.181	−14.87	42.26	7.68	31.51	1.96

Notes. The symbol "m" stands for main pulses, while "i" means interpulses. ^aα and β were constrained by fitting the multi-frequency pulse widths of average profiles, the width of the average drift profile, the fractional drift intensity, the interval between successive sub-pulses in the same period and parts of PPA (without the orthogonal mode jumps). ^bα and β were constrained by fitting the PPA sweep rate and multi-frequency pulse widths. ^cα and β were constrained by fitting multi-epoch pulse profiles and PPA variations with a precessional model. References. (1) Smits et al. (2007); (2) Gould & Lyne (1998); (3) Manchester et al. (1998); (4) Johnston et al. (2008); (5) Theureau et al. (2011); (6) Everett & Weisberg (2001); (7) Gil & Lyne (1995); (8) Keith et al. (2010); (9) Becker et al. (2005); (10) Johnston et al. (2007); (11) Rankin et al. (2006); (12) Krishnamohan & Downs (1983); (13) Taylor et al. (1993); (14) Rankin & Wright (2007); (15) Kramer & Johnston (2008); (16) Hobbs et al. (2004); (17) Johnston & Weisberg (2006); (18) Gangadhara et al. (1999); (19) Manchester et al. (2010); (20) Wu et al. (1993); (21) Crawford et al. (2001b); (22) Camilo et al. (2008); (23) Tiburzi et al. (2013); (24) von Hoensbroech et al. (1998); (25) van Ommen et al. (1997); (26) Mitra & Rankin (2008); (27) Desvignes et al. (2013); (28) Kasian (2012); (29) Weisberg & Taylor (2002); (30) Taylor & Weisberg (1982); (31) Boyles et al. (2013); (32) Noutsos et al. (2011); (33) Qiao et al. (1995).

Download table as:  ASCIITypeset images: 1 2

The above simulation shows that the pulse width W increases with increasing impact angle |β| in the outer beam. Below we derive the relationship, which will be used in Section 3.2 to test our model.

The above points (2)–(5) suggest that probably only a few flux tubes are active and visible, and we simply assume that the effective emission region is confined within an azimuth range of [ − φ, φ] around the meridian plane (note that this range may contain more than one flux tube). According to the spherical geometry relationships, we have

$$\begin{equation} \cos \zeta =\cos \alpha \cos \rho -\sin \alpha \sin \rho \cos \varphi, \end{equation} \tag{ 23 }$$

$$\begin{equation} \cos \rho =\cos \alpha \cos \zeta +\sin \alpha \sin \zeta \cos \Phi, \end{equation} \tag{ 24 }$$

and

$$\begin{equation} \frac{\sin \rho }{\sin \Phi }=\frac{\sin \zeta }{\sin \varphi }, \end{equation} \tag{ 25 }$$

where ζ = α + β is the viewing angle between the LOS and the rotating axis, and Φ is the pulse longitude.

Substituting Equations (24) and (25) into Equation (23), and using Φ = W/2, the formula is simplified as

$$\begin{equation} \cos \left(\frac{W}{2}+C\right)=\frac{\sin \alpha }{\tan (\alpha +\beta)(\cos ^2\alpha +\tan ^{-2}\varphi)^{1/2}}, \end{equation} \tag{ 26 }$$

where $C=\arctan (\sec \alpha /\tan \varphi)$. It can be found that regardless of whether β > 0 or β < 0, W will increase with |β|, which is shown by Figure 15(a). Note that Equation (26) only applies to the outer beam.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Deriving the W − |β| relationship for the inner beam is difficult because of the reversal intensity distribution. In some cases, the LOS may see emissions from other less active flux tubes far away from the meridian plane (see β = 1° and 2° in Figure 12 for examples). Since one tends to see more parts of the polar cap, ρ_b ≃ 1 ∼ 2ρ_pc may be a rough approximation for the beam size when |β| ≲ ρ_t.⁶

We must point out that the above mentioned pulse width is completely determined by the boundary of the flux tube, irrespective of the intensity distribution in the emission beam, hence it can be called the geometric pulse width. The most often used pulse width, W₅₀ or W₁₀, which is measured at the 50% or 10% level of pulse peak, obviously depends on the intensity distribution. It is possible that the intensity drops so abruptly as the LOS sweeps toward the lateral beam edges that W₁₀ becomes smaller than the geometric width W. Will the increasing trend of the W − |β| relationship still hold? Below we demonstrate that even in the extreme case such a trend still exists.

We consider an extreme case that the intensity is constant for any circle around the beam center. In the outer part of the beam, the intensity follows a simple limb-darkening relation, I = I₀ρ^δ (δ < 0). Given an LOS with β, one will see the maximum intensity when ρ reaches its minimal value, i.e., ρ = |β|; therefore, I_pk = I₀β^δ. For a level ηI_pk (η < 1), at which the pulse width is defined, the boundary intensity drops to I_b = ηI_{, pk}. Using the limb-darkening relation, one finds the boundary radius ρ_b = η^1/δβ. Then the corresponding pulse phase can be figured out via cos Φ_b = (cos ρ_b − cos αcos ζ)/(sin αsin ζ) (from Equation (24)), and the pulse width will be

$$\begin{equation} W=2\Phi _{\rm b}=2\arccos \left[\frac{\cos (\eta ^{1/\delta }\beta)-\cos \alpha \cos \zeta }{\sin \alpha \sin \zeta }\right]. \end{equation} \tag{ 27 }$$

Again, regardless of whether β > 0 or β < 0, W will increase with |β|, which is shown by Figure 15(b). This tendency can also be found in the profiles in Figure 14.

In the following section, we will use Equation (26) to calculate the pulse width when |β| > ρ_pc, but we fix the beam size as 2ρ_pc when |β| ⩽ ρ_pc.

It is known that binary pulsars may undergo relativistic spin-precession due to coupling between the spin and orbital angular momenta (Damour & Ruffini 1974; Barker & O'Connell 1975). As a result, the spin axis of the pulsar rotates around the total angular momentum vector, changing the viewing and impact angles. This effect enables us to "scan" the emission beam. So far, efforts to construct the beam structure have been made for six pulsars, PSR B1913+16 (e.g., Weisberg & Taylor 2002), B1534+12 (e.g., Arzoumanian 1995), PSR J1141−6545 (Manchester et al. 2010), PSR J1906+0746 (Desvignes et al. 2013), PSR J0737−3039A (Ferdman et al. 2013; Perera et al. 2014), and PSR J0737−3039B (e.g., Perera et al. 2012; Lomiashvili & Lyutikov 2014). For half of them, PSR J1141−6545, PSR J1906+0746, and PSR J0737−3039B, the two-dimensional beam structures were constructed with multi-epoch absolute flux density data; therefore, they can be used to probe the models via both the evolution of the pulse width and the radial distribution of emission intensity. But because the magnetosphere of PSR J0737−3039B is distorted by the wind from its companion PSR J0737−3039A, it is not an ideal case for the model test. For the other three pulsars, the profile at each epoch is first normalized by a pulse peak and then used to study the profile evolution or construct the beam structure. This method focuses on the evolution of relative intensity of different pulse components rather than the radial distribution of absolute intensity in the pulse beam. Therefore, only the evolution of the pulse width is useful to test the fan beam model for these 3 pulsars. Below we first present clear evidence from PSR J1141−6545 and PSR J1906+0746, and then we discuss the remaining pulsars, which provide less prominent evidence.

PSR J1141−6545 is a young binary pulsar with a precession rate of 1.4 deg yr⁻¹. Manchester et al. (2010) found that the average pulse profiles of this pulsar show remarkable variation at 1.4 GHz between 1999 and 2008. In the first a couple of years, the profile was dominated by the trailing part. The leading part grew stronger with time and later became comparable with the trailing part, leading to a bump-shaped profile. Despite the profile evolution, the pulse width at a very low intensity level, e.g., 1% of the peak, is nearly constant. The authors used a precessional beam model to fit the pulse profile and the absolute central PPAs obtained by fitting PPA data with the RVM. The derived impact angles varied from about −37 in 1999 to −09 in 2007, meaning that the LOS was moving toward the magnetic pole. With the data of flux density, the two-dimensional beam intensity structure was inferred, which shows that the maximum intensity was reached when β ≃ −2°. Beyond this angle, the intensity decreases as our LOS goes further to the magnetic pole, while below this angle, the intensity decreases as the LOS moves toward the magnetic pole, as indicated by two opposite arrows in the bottom left panel of Figure 16. The most striking feature is that the beam is quite asymmetric and is partially filled without any core or conal structure. On the basis of this point, the authors concluded that the beam is patchy.

The twofold radial intensity distribution is generally consistent with our model. The observed transition radius of the intensity distribution, ρ_t ≃ 2°, is close to the opening angle of the polar cap ρ_pc ≃ 2° for this pulsar with a period of 0.394 s. This suggests that the lowest coherent emission altitude should be close to the polar cap. The nearly constant pulse width does not conflict with our model. We have discussed in Section 2.4 that the pulse width near the transition radius may not follow the same increasing trend with impact angle as it does for the outer beam.

PSR J1906+0746 is a young binary pulsar with a higher precession rate of 2.2 deg yr⁻¹. The radio profile at 1.4 GHz has a narrow main pulse and a weak interpulse, which are separated by about 180°. The PPA data are smooth, following a simple RVM (Desvignes 2010). The RVM fitting at each epoch revealed that α ≃ 81°, while β varied from 68 to about 11° between 2005 and 2009 for the main pulse and from 112 to 7° for the interpulse (Desvignes et al. 2013).

The derived beams for both poles show asymmetric structures that are quite patchy. In both of the beams, the intensity decreases as the LOS moves away from the magnetic pole (see the bottom right panel of Figure 16 for the main pulse). Since the LOSs in both poles are much further from the polar cap boundary (ρ_pc = 33 for P = 0.144 s), this limb-darkening feature is consistent with the radial intensity trend for the outer part of the fan beam. The pulse width is nearly constant for the main pulse, but it is not clear for the interpulse because of poor signal-to-noise ratios.⁷ The near constant pulse width may reflect the change of the transverse intensity distribution at different radii, suggesting that the real physical process could be more complicated than our assumptions.

It is worth noting that the absence of emission in the other parts along our LOS, e.g., the black point A in the bottom right panel of Figure 16 that has the same radius as the white point B in the bright beam, is unlikely due to the limb-darkening effect. We suggest that perhaps only one flux tube is active in each pole for PSR J1906+0746; if there are more, they must be far away from the median and less active. As for PSR J1141−6545, two active flux tubes seem to be responsible for the asymmetric beam.

In order to compare with the observation, we simulate the model beams with model B2 for PSR J1141−6545 and the main pulse of PSR J1906+0746, as shown by the top panels in Figure 16. The observational beams of PSR J1141−6545 and PSR J1906+0746, shown in the bottom panels, are taken from Manchester et al. (2010) and Desvignes et al. (2013), respectively.

For PSR J1141−6545, two flux tubes are used in the simulation. Both of them have χ₀ = 0.7, ℜ₀ = 0.3, and the same peak multiplicity factors, but with different central azimuths, viz. φ_{0, a} = 0 and φ_{0, b} = 35°, as indicated by short and a long white arrows in the top left panel. The multiplicity is assumed to follow the same two-dimensional Gaussian distribution across the cross section of both flux tubes. Another important difference is the radial limb-darkening index, which is assumed to be δ = −4 for the flux tube at the meridian plane and δ = −2 for the lateral one. Although the simulated beam is not a satisfactory reproduction of the observed beam, the general intensity structure is similar.

For the main pulse beam of PSR J1906+0746, only one flux tube with χ₀ = 0.7, ℜ₀ = 0.3, and φ₀ = 40° is used in the simulation. The arrows in the simulated and observed beams represent approximately the direction of the radial intensity gradient. The simulations presented here are for comparing the general intensity distribution pattern between the simple models and observations. Some complex local structures in the observed beams are not reproduced and require detailed models.

PSR B1913+16 has a precession rate of 1.2 deg yr⁻¹. The radio pulse profile at 1.4 GHz consists of double peaks and a shallow bridge. Although the peak separation gradually decreases with time, the pulse width measured at a level below 50% peak intensity is roughly constant between 1981 and 2003 (Weisberg & Taylor 2005; Clifton & Weisberg 2008). In fact, the pulse width at very low intensity levels, e.g., below 10%, after keeping constant for about 10 years, has undergone a subtle increase since 1992. By fitting the profiles from 1981 to 2001 with a geometrical model of precessing radio beams and comparing the slope rate of the PPA curve expected by the RVM and the observed data, Weisberg & Taylor (2002) determined that the impact angle varied from about −36 in 1981 to −66 in 2001. Combining with the best-fit inclination angle α ∼ 158°, the negative β angles mean that our LOS was sweeping across the beam between the magnetic pole and the equatorial plane and was moving further from the magnetic pole.

Obviously, a circular conal beam model cannot account for the constant and subtle increasing trend of pulse width in such viewing geometry. The authors proposed that the beam should be elongated in the latitudinal direction and pinched in longitude near the center, forming a hourglass-shaped beam (Weisberg & Taylor 2002, 2005; Clifton & Weisberg 2008), which works well in modeling the historical data. Although the subtle pulse broadening is not as significant as the above simple fan beam model predicts, it is possible to be explained by taking into account some modifications, e.g., the aberration and retardation effects and the rotating dipole field and intensity modulation that may depend on emission altitude and azimuth. The narrowing of the peak separation and the increasing bridge emission may also be explained by this modified version if an appropriate number of flux tubes, e.g., three (corresponding to the leading, bridge, and trailing components), with different altitude-dependent intensity distributions are assumed. We notice that the current hourglass model predicts decreasing pulse width between 2003 and 2020 (Clifton & Weisberg 2008) and thus data in the near future will hopefully provide a clear test for the hourglass beam model and our fan beam model. The other quick test is to examine the evolution of flux density, if they exist in previous observations. We did not find useful information in regards to this matter, because the multi-epoch profiles have been normalized before constructing two-dimensional beam maps in published papers.

PSR B1534+12, a 37.9 ms millisecond pulsar (MSP) with a broad main pulse and interpulse, locates in a binary system orbiting with another neutron star (Wolszczan 1991). The most recently measured precession rate is 059 yr⁻¹ (Fonseca et al. 2014). The main pulse width of PSR B1534+12 at a 3% level of peak intensity is roughly increasing at a rate of 05 ∼ 1° per year as estimated from the profiles in Arzoumanian (1995) and Stairs et al. (2000). This tendency of profile broadening is confirmed by later observations (Stairs et al. 2004; Fonseca et al. 2014), accompanied by a secular decreasing trend in intensity of the central component with respect to the leading and trailing wings. The impact angle is constrained by fitting the observed PPA with RVM, which changed from about −33 in 1993 to −68 in 2009, indicating that the LOS was moving away from the magnetic pole (Fonseca et al. 2014). The tendency that the main pulse width increases with increasing impact angle (absolute value) is generally consistent with our model. It would be interesting to check whether the flux density decreases with time, as the fan beam model predicts, but the flux density data were not presented in the above references because of the normalization of multi-epoch profiles.

PSR J0737−3039A, a 22.7 ms MSP in the double pulsar system, has a predicted precession rate of 48 yr⁻¹ (Burgay et al. 2003), which is comparable to that of PSR J0737−3039B. However, the multi-epoch pulse profiles of PSR J0737−3039A were found to be fairly stable (Manchester et al. 2005; Ferdman et al. 2013), which is explained as a consequence of a very small misalignment between the pulsar spin axis and the orbital momentum. The broad main pulse and interpulse were modeled by circular cones, where the inclination angle was constrained to be nearly 90° and the two cones should originate from two opposite poles (Ferdman et al. 2013; Perera et al. 2014). We would suggest that the profile can be alternatively modeled by the fan beam model. Since our LOS is stable for this pulsar, it is unable to distinguish the beam models directly. Nevertheless, we notice that PSR J0737−3039A and PSR B1534+12 have some features similar to many other MSPs, e.g., a broad pulse width and complex profile shape. In the coming paper II, we will discuss how the fan beam model, in the context of emission flux tubes, has advantages for explaining these features for MSPs.

By numerically modeling the distortion of PSR J0737–3039B's magnetosphere induced by PSR J0737–3039A, Lomiashvili & Lyutikov (2014) determined that the emission beam is horseshoe shaped, which is somewhat similar to the arclike structure derived by Perera et al. (2010, 2012), and the emission region is located at about 3750 stellar radius (∼30% light cylinder radius), which is greater than the maximum altitude of 2500 stellar radius given by Perera et al. (2012). Considering that the orbital-phase dependent distortion of PSR J0737–3039B 's magnetosphere may influence the intrinsic emission beam structure, we choose to discard this pulsar.

In Section 2 we have derived two major relationships for the fan beam model, i.e., the radial limb-darkening relationship and the W − |β| relationship. A third relationship between the upper boundary of |β| and pulsar distance d will be derived in this section. In order to test these relationships, we collected a sample of pulsars with known pseudo radio luminosity L_app = F_νd², pulse width W₁₀, d, and the impact angle, where F_ν is the flux density at a particular frequency.

Unlike the other parameters, the impact angle cannot be directly measured. In literature, at least four kinds of methods have been proposed to constrain the impact and inclination angles. (1) When a pulsar presents a smooth S-shaped PPA curve, α and β can be obtained by fitting the PPA data with the following relation given by the RVM:

$$\begin{eqnarray} \tan (\Psi -\Psi _0)&=&\sin \alpha \sin (\Phi -\Phi _0)[\sin (\alpha +\beta)\cos \alpha \nonumber \\ &&-\cos (\alpha +\beta)\sin \alpha \cos (\Phi -\Phi _0)]^{-1}, \qquad \end{eqnarray} \tag{ 28 }$$

where Ψ is the PPA, and the subscript "0" denotes the values at the phase where the slope rate of PPA curve reaches its extremum. (2) The parameters can be constrained by fitting the pulse width and some properties of subpulse drifting for a few pulsars, e.g., the interval between successive subpulses in the same period. (3) For a few binary pulsars with a considerable precession rate, the parameters can be derived by fitting the profiles and PPA variation in terms of a precession model. (4) For gamma-ray pulsars, the parameters can be constrained by fitting the radio and gamma-ray profiles with gamma-ray emission and radio conal beam models. In this paper, we prefer the first three geometrical methods to avoid the dependence of emission models, and the related references can be found in Table 1.

We collected a sample of 64 normal pulsars⁸ from literature, of which the inclination and impact angles are derived with method (1) for 62 pulsars (including the binary pulsar PSR J1906+0746) and with methods (2) or (3) for the other 2 pulsars. A handful of pulsars that have too large of uncertainties for the impact angle are not included. Twelve pulsars in the sample have interpulse emission and hence contribute beam information from double poles. Table 1 gives the data for the total of 76 beams. The second column includes the inclination angles α, impact angles β, frequencies at which α and β are derived (with marked numbers for references), 10% peak pulse widths W₁₀ with uncertainties, pulse widths $W_{10}^1$ at lower frequencies and $W_{10}^2$ at higher frequencies, pairs of frequencies at which $W_{10}^1$ and $W_{10}^2$ are measured, number of profile components N_c that are identified by eye, pulsar periods P, period derivatives $\dot{P}$, pseudo luminosities at 400 MHz L₄₀₀ and at 1400 MHz L₁₄₀₀, rotation energy loss rates $\lg \dot{E}$, and pulsar distances d.

For most pulsars, the uncertainty of the pulse width induced by the frequency dependence of profiles is larger than the observational error at a single frequency. To count this major error source, we collected the pulse widths at two well-separated frequencies ($W_{10}^1$ and $W_{10}^2$), mostly at 0.4 MHz and 1.6 MHz. The uncertainty is then figured out by $|W_{10}^2-W_{10}^1|/2$. For those pulsars with only one frequency observation in literature, we assume an uncertainty of 10% for W₁₀. The data in the last six columns are taken from the ATNF Pulsar Catalog (Manchester 1995).

3.2.1. Test of the W − |β| Relationship

According to Equation (26), the pulse width depends on α, β, and the azimuthal width of a flux tube Δφ = 2φ. These parameters are different for pulsars, causing dispersion of pulse widths. In order to compare the model with observations, we perform a Monte Carlo simulation by randomly assigning α, β, and Δφ to a sample of ∼50,000 pulsars. If the projection of the magnetic pole is uniformly distributed in the celestial sphere, then the probability density function of inclination angle is P(α) = sin α/2.⁹ The probability density function of viewing angle is also assumed to be P(ζ) = sin ζ/2, and Δφ is assumed to be uniformly distributed between Δφ₁ and Δφ₂, where the boundaries are free parameters that can be estimated by comparing the simulated results with observational data.

Since we assume that the effective beam radius can be treated as ρ_b = 1 ∼ 2ρ_pc when the LOS sweeps across the inner part of the fan beam, the pulsar period will also affect the pulse width through its influence on ρ_pc. We assume a lognormal distribution for the pulse period, i.e.,

$$\begin{equation} p=\frac{1}{P\sigma \sqrt{2\pi }}\exp \left [-\frac{(\ln P-\mu)^2}{2\sigma ^2}\right ], \end{equation} \tag{ 29 }$$

where p is the probability density function of pulsar period P. The best-fit parameters μ and σ are obtained by fitting the ATNF data of the pulsar period longer than 50 ms, which are μ = −0.48 (corresponding to a peak period 0.62 s) and σ = 0.90. The period is also randomly assigned to each pulsar together with the other three parameters.

Once a pulsar is assigned with a group of α, ζ, P, and Δφ, the pulse width is calculated separately for two circumstances: using Equation (26) when |β| = |ζ − α| > 2ρ_pc and fixing it as W = 2Φ, which is figured out by substituting the assigned α, beam radius ρ_b = 2ρ_pc and ζ = α into Equation (24) when |β| ⩽ 2ρ_pc. Selecting |β| = 2ρ_pc as the criterion to separate the two circumstances and assuming the beam radius as 2ρ_pc are phenomenological choices to make the simulation generally coincide with the lower boundary of the observed W₁₀ for small impact angles (see Figure 17). But it can be reasonably explained in the context of the fan beam model, as will be shown below.

Zoom In Zoom Out Reset image size

Given the other parameters, smaller inclination angles tend to produce wider profiles. To avoid possible contamination of this effect, we divide the sample into three groups with different α, i.e., group A with α ⩽ 30°, group B with 30° < α ⩽ 60°, and group C with 60° < α ⩽ 90°. The inclination angles in Table 1 that are larger than 90° are converted to 180° − α before grouping, and the corresponding β are converted to −β. Each group has more than a dozen of pulsars, which are shown by the black dots in Figure 17. The simulated pulsars are also divided into three groups and are plotted by the gray dots.

One can see clearly the twofold W − |β| relationship from the observed data in Figure 17: W is almost independent to |β| when |β| < 6 ∼ 8°, while it increases with increasing |β| when |β| > 6 ∼ 8°. When selecting Δφ₁ and Δφ₂ as 40° and 160°, respectively, we find that the observed distribution is well reproduced. For comparison, the W − |β| relationships of Equation (26) under various groups of α and φ are plotted as dashed and dotted curves, e.g., curve "1" for α = 5° and φ = 80°. Obviously, such a simple equation only accounts for the W − |β| relationship for the outer beam.

These curves are helpful for understanding why we select |β| = 2ρ_pc as the criterion to distinguish the outer beam and the inner beam and set the effective beam radius as ρ_b = 2ρ_pc. We have tried some other criteria, but the simulations do not fit the lower part of the observed data. For instance, when selecting |β| = ρ_pc and ρ_b = ρ_pc, many simulated data points extend along these curves down to the small pulse width region, which is well below the lower boundary of the observed data points. When selecting |β| = 3ρ_pc and ρ_b = 3ρ_pc, the lower boundary of simulated widths will be higher than that of the observed data. That is the reason why we made the final choices of 2ρ_pc. A possible interpretation is that in the statistical sense, the lowest coherent emission altitude may be a few times of R so that the transition radius of the twofold intensity distribution is ρ_t ∼ 2ρ_pc, which then leads to twofold W − |β| relationships separated by |β| ∼ 2ρ_pc.

Two possible selection effects can be excluded for the increasing W₁₀ − |β| trend. First, we have checked the data of W₁₀ and pulsar period P. No trend is found between these two parameters for the sample, indicating that the increasing W₁₀ trend is not induced by possible dependence of the pulsar period. Second, it is known that the inclination and impact angles are hardly constrained with the RVM fitting method for very narrow profiles, because the PPA curves for a number of α-β pairs that produce the same maximum slope rate k = sin α/sin β can be barely distinguished in a very narrow phase interval around the flection point of the PPA curve. Then one may suspect whether this limitation can cause bias and hence the observed trend cannot be regarded as a robust evidence for the fan beam model. Especially when one considers an opposite case that the pulse width decreases with increasing |β| as predicted by the conal beam model, this limitation might cause the lack of samples with large impact angles (due to narrow pulse widths) and would be unfavorable for testing the prediction of the conal beam model. However, the limitation can be overcome by the large phase separation between the main pulse and interpulse when both of them are observed, even though the pulses are narrow. For the 25 main and interpulse beams in our sample (only the main pulse width is measurable for PSR J1932+1059), 11 cases have |β| > 8° (6 with |β| ⩾ 12°). With so many large-impact-angle samples, no violation of the increasing trend of W₁₀ is observed; therefore, this selection effect can be ruled out.

Finally, one may notice that the W − |β| relationships predicted by Equations (26) and (27) are actually different for positive and negative impact angles for a given α and other parameters and that the difference grows with decreasing inclination angle, as shown by the solid and dashed curves in Figure 15. However, the predicted difference is not seen in the distribution of current data in Figure 17, where the positive and negative impact angles are represented by black dots and open circles, respectively. The reason is probably that the scatter due to other parameters, e.g., α, φ, or δ, overtakes the difference caused by the sign of the impact angle (see the dashed curves in Figure 17 for the scatter induced by α and φ). Perhaps this difference can only be tested for a large sample of pulsars.

3.2.2. Test of the Intensity–Radius Relationship

The observed mean flux density, F = F_peakW_e/P, is a quantity averaged over the whole pulse period, where F_peak is the flux density at the maximum peak of the averaged pulse profile, and W_e/P is the pulse duty cycle. Combining with the pulsar distance and other parameters, F_peak can be converted to the emission intensity I_peak at ρ_peak in the radio beam, where ρ_peak is the radius for the maximum pulse peak. The fan beam model predicts that such an intensity follows I_peak∝ρ^{2q − 6} in the outer beam. Therefore, the derived intensity can be used to test this intensity–radius relationship.

The peak flux at a frequency range of ν ± Δν/2 can be converted to intensity by I_{peak, ν} = F_{peak, ν}d²Δν, namely I_{peak, ν} = F_ν(W_e/P)⁻¹d²Δν. Using Equation (10), we have

$$\begin{equation*} I_{\rm peak,\nu }=AP^{q-4}\dot{P}\cos ^2\alpha \left(f_{1}^{q-3}-f_{2}^{q-3}\right)\rho _{\rm peak}^{2q-6}. \end{equation*}$$

In order to display the dependence of ρ_peak clearly, we defined an intensity indicator parameter,

$$\begin{equation} Y=I_{\rm peak,\nu } A^{-1}P^{4-q}\dot{P}^{-1}\cos ^{-2}\alpha \left(f_1^{q-3}-f_2^{q-3}\right)^{-1}, \end{equation} \tag{ 30 }$$

so that

$$\begin{equation} Y=\rho _{\rm peak}^{2q-6}. \end{equation} \tag{ 31 }$$

Inserting I_{peak, ν} into the above equation, there is $Y=F_\nu d^2(W_{\rm e}/P)^{-1}\Delta \nu A^{-1}P^{4-q}\dot{P}^{-1}\cos ^{-2}\alpha (f_1^{q-3}-f_2^{q-3})^{-1}$, where F_νd² is the pseudo luminosity at a particular frequency ν.

When trying to test the Y − ρ_peak relationship, it is better to use

$$\begin{equation} Y=K\rho _{\rm peak}^{2q-6} \end{equation} \tag{ 32 }$$

instead of Equation (31) to fit data, where K is a coefficient to be constrained, because there are several free parameters in Y, which may cause uncertainties. These parameters include the boundaries of flux tube f₁ and f₂, index q, multiplicity M, stellar radius R, and emission power of a single particle i_e0. In the following calculation we take M = 10³, R = 10⁶ cm, i_e0 = 10⁻¹⁵ erg s⁻¹, f₁ = 1, and f₂ = 100. In fact, choosing different values does not affect the power-law index and the I_outer − ρ_peak relationships recovered later.

The pseudo luminosities L₄₀₀ = F₄₀₀d² and L₁₄₀₀ = F₁₄₀₀d² at 400 MHz and 1400 MHz are taken from the ATNF pulsar catalog, and the frequency range Δν is set as 100 MHz and 500 MHz, respectively. The W_e/P is simply replaced by W₁₀/360°, which may introduce an maximum uncertainty by a factor around two or three (see data in Gould & Lyne 1998). However, this uncertainty is much smaller than the uncertainties caused by the free parameters. Notice that cos α ∼ 0 when α ∼ 90°, which means that the net charge number density n_gj0 ∼ 0 and hence n₀ = Mn_gj0 ∼ 0 on the polar cap surface. We believe that the real secondary particle density of orthogonal rotators should be much higher than this and may be comparable to other cases with moderate and small α; therefore, we simply replace the correction term cos ⁻²α by a factor of one.

Finally, the free parameter q needs special treatment, because it affects both sides of Equation (32). We let it vary from –3 to 3 by a very small step size.¹⁰ Given a q value, Y is calculated for individual pulsars, and the least square fitting is applied to the Y and ρ_peak data to find a best-fit index δ for the relationship $Y\propto \rho _{\rm peak}^\delta$. If the relationship of Equation (32) does apply to the sample, there should be a q value satisfying 2q − 6 = δ. Therefore, the fitting process is repeated for the whole range of q to search for the solution. Because Equation (32) is only valid for the outer beam, in the following test the pulsars with |β| ⩽ 2ρ_pc are excluded.

When calculating ρ_peak with the peak pulse phase Φ_peak, we have considered the effect of profile shape on Φ_peak. The peak phases are set as Φ_peak ≃ W₁₀/4, 3W₁₀/8, and 5W₁₀/12 for profiles with double, quadruple, and sextuple components, where the profile center is always assumed to have Φ = 0. While the component number is odd, we simply assume that the maximum peak occurs at the profile center. Then ρ_peak is figured out with Equation (24) for each pulsar. The number of components, as listed in Table 1 by the column N_c, are identified by eye for each pulsar according to its multifrequency profile shapes whenever they are available in literature.

The right panels of Figure 18 show the modeled index 2q − 6 and the fitted index δ as a function of q when the fitting process is applied to 400 MHz and 1400 MHz data, respectively. The intersections give the solution q = 1.5 for 400 MHz and q = 1.75 for 1400 MHz. The left panels of Figure 18 present the Y − ρ_peak diagrams when the above solutions of q are adopted. Despite the considerable scattering, the data show a trend in which the intensity decreases when LOS becomes further from the beam center. The best-fit limb-darkening relationship obtained with the Levenberg–Marquardt method is $Y=10^{-1.6\pm 1.4}\rho _{\rm peak}^{-3.0\pm 1.5}$ for 400 MHz and $Y=10^{-1.7\pm 1.5}\rho _{\rm peak}^{-2.5\pm 1.5}$ for 1400 MHz at the 95% confidence level, respectively. Combining with the results for 400 MHz and 1400 MHz, we have an approximate limb-darkening relationship I_out∝ρ^{−(2.8 ± 1.8)}, with the index q ≃ 1.6 ± 0.9.

Zoom In Zoom Out Reset image size

Using the best-fit Y − ρ_peak relationships, we can recover the I_outer − ρ_peak relationships with the following equation:

$$\begin{equation*} I_{\rm outer}=YA\left(f_1^{q-3}-f_2^{q-3}\right)\Delta \nu _{\rm MHz}^{-1}P^{q-4}\dot{P}_{-15}10^{-15}, \end{equation*}$$

where I_outer is in unit of erg/s/MHz/sr, P is in unit of seconds, ρ_peak is in unit of degrees, $\dot{P}_{-15}=\dot{P}/(10^{-15}\,{\rm s/s})$, and Δν_MHz = 100 and 500 for 400 MHz and 1400 MHz, respectively. Substituting into the above parameters, we have

$$\begin{eqnarray} I_{\rm outer}^{400\,{\rm MHz}} &=& 10^{27.2\pm 1.4}P^{-1.50\pm 0.75}\nonumber \\ &&\times \dot{P}_{-15}\rho _{\rm peak}^{-3.0\pm 1.5}\ {\rm erg\ s^{-1}\ MHz^{-1}\ sr^{-1}} \end{eqnarray} \tag{ 33 }$$

and

$$\begin{eqnarray} I_{\rm outer}^{1400{\rm MHz}} &=& 10^{25.7\pm 1.5}P^{-2.25\pm 0.75}\nonumber\\ &&\times\dot{P}_{-15}\rho _{\rm peak}^{-2.5\pm 1.5} \ {\rm erg\ s^{-1}\ MHz^{-1}\ sr^{-1}}. \end{eqnarray} \tag{ 34 }$$

Although the above I_outer stands for the intensity at ρ = ρ_peak in the emission beam, it is reasonable to believe that these empirical relationships, in the statistical sense, can be used to describe the radial limb-darkening relationship for pulsar radio beams.

3.2.3. Impact Angle versus Pulsar Distance

Given an observing sensitivity, luminous pulsars have more of a chance to be detected at large distances. According to the radial limb-darkening relation in our model, a higher luminosity requires a smaller impact angle. Then, for a sample of pulsars, the upper limit of impact angle should decrease with pulsar distance. In the nearby region to the earth, less luminous pulsars may still reach the sensitivity, thus the impact angles are expected to be more scattered than in further regions.

To find out the relation between the distance d and |β|, let us assume a fixed minimum detectable flux density $S_{\rm min}=C_{\rm s}\sqrt{W/(P-W)}$ (Lorimer & Kramer 2005), where W is the equivalent pulse width, and C_s is a constant determined by some observational parameters, such as bandwidth, total integration time, etc. The corresponding peak flux density is $S_{\rm peak}=S_{\rm min}P/W\simeq C_{\rm s}\sqrt{P/W}$ when P − W ≃ P. In the fan beam model, an intensity $I_{\rm min}\propto |\beta |_{\rm max}^\delta$ would be detectable if S_peak ∼ I_min/d². This leads to

$$\begin{equation} |\beta |\propto C_{\rm s}^{2/(1+2\delta)}d^{4/(1+2\delta)}, \end{equation} \tag{ 35 }$$

if we use W∝|β| (a viable approximation to Equations (26) and (27), see Figure 15).

Using the ATNF data of distance, we plot the |β| − d diagram for the sample (Figure 19). It does show that the impact angles are less scattered at larger distances and the upper boundary of |β| decreases with distance until ∼15 kpc, where no more pulsars exist in the sample. The sparse data points at d < 1 kpc are probably due to incompleteness of the data set, but this does not affect the trend of the upper limit. To compare the above relationship with the data, we plot two |β| − d curves with δ = −2.0 and −3.5, which correspond to |β|∝d^−1.3 (dashed) and |β|∝d^−0.7 (solid), respectively. The |β| is restricted within 90°, because a β larger than 90° means viewing into the other pole, and it can be replaced by its supplementary angle. Therefore, the theoretical upper boundary at small distances is replaced by |β| = 90°. As shown in the figure, the data points are all placed under these boundary lines, and the radial limb-darkening relationship with δ ∼ −3 generally matches the upper boundary of the data.

Zoom In Zoom Out Reset image size

In order to examine if this consistency can extend to a larger distance, we have reviewed the currently known most distant pulsars in the Magellanic Clouds. Unfortunately, there is little polarization observation for these weak sources. But it was found that their pulse profiles are much narrower than those of Galactic pulsars (e.g., Manchester et al. 2006). This is a sign of small impact angle, if we believe that these two populations follow the same physics. We estimated W₁₀ for all the 21 normal pulsars in the Large Magellanic Cloud (LMC)¹¹ and the 5 normal pulsars in the Small Magellanic Cloud (SMC) using the published profiles (McConnell et al. 1991; Crawford et al. 2001a; Manchester et al. 2006; Ridley et al. 2013). For a few pulsars with the signal-to-noise ratio lower than 10, the lowest level pulse width is measured. The W₁₀ of LMC pulsars varies from 6° to 45°, with only 4 larger than 22°. The averaged W₁₀ is 15° for the 17 small-width pulsars and 20° for all the 21 pulsars. The W₁₀ of SMC pulsars varies from 11° to 22°, with an average value about 15°. These are considerably smaller than the averaged W₁₀ of 27° of our sample. As shown by the data points in Figure 17, most pulsars with W₁₀ ⩽ 20° have |β| ⩽ 10°; therefore, if there is no difference in the statistical property of emission beams between MC and Galactic populations, it can be inferred that most LMC and SMC pulsars have impact angles of |β| ⩽ 10°. As to the 4 LMC pulsars with larger W₁₀, they may also have small impact angles if the inclination angle is not too large, say α < 60° (see the top two panels in Figure 17).

The expected range of impact angles are plotted in Figure 19 with two short lines at d = 50 kpc for LMC (Pietrzyński et al. 2013) and d = 60 kpc for SMC (Hilditch et al. 2005). Again, the combined upper boundary of Galactic and MC samples is consistent with the limb-darkening relation with the index around −3.

Some other effects may influence the |β| − d distribution, but they cannot account for the trend of the upper boundary lines. The following analysis will strengthen the conclusion that the impact-angle-distance upper boundary line is evidence for our model.

• 1.  
  Given the luminosity and pulsar distance, a narrow pulse profile is easier to be detected than a wide pulse profile. To estimate the impact of this selection effect quantitatively, one has to assume a relationship between the pulse width and the impact angle. A possible estimation can be made as follows. Given a sensitivity of $S_{\rm min}=C_{\rm s}\sqrt{W/(P-W)}$, a pulsar with luminosity L will be detectable at d when its flux density is comparable to S_min, namely, L/d² ∼ S_min. A simple derivation leads to W∝d⁻⁴. If still using the approximation in the fan beam model, W∝|β|, we would have |β|∝d⁻⁴, which is too steep to account for the flat upper boundary line in the |β| − d diagram. In order to account for the apparent power-law index of the upper boundary line (∼−1), one has to assume W∝|β|^∼4, which is not reliable.
• 2.  
  The disparity in the sensitivity of different pulsar searching projects is unlikely to cause severe bias to the upper boundary line. For example, the limiting sensitivity is 0.14 mJy for Parkes multibeam survey (PMS, Manchester et al. 2001) and 0.08 mJy for 5% duty cycle for the survey of MC pulsars (Manchester et al. 2006). This difference can be denoted as C_{s, MC}/C_{s, PMS} = 4/7. Using Equation (35), the sensitivity of 0.08 mJy causes an increment of |β| by about 45% (or 25%) for MC pulsars compared with the |β| estimated with the sensitivity of 0.14 mJy, if δ = −2 (or −3) is used. Such an increment does not deviate much from the upper boundary line in Figure 19.
• 3.  
  Although the pulsar period is not written in Equation (35), it does affect that relationship in our model. However, no statistical distinction is found for P in all the distance ranges, including the MC pulsars. Therefore, the effect of P should be trivial.