Radiation by the superluminally moving current sheet in the magnetosphere of a neutron star

The mechanism by which the radiation received from obliquely rotating neutron stars is generated remains an open question half a century after the discovery of pulsars. In contrast, considerable progress has recently been made in determining the structure of the magnetosphere that surrounds these objects: numerical computations based on the force-free, magnetohydrodynamic and particle-in-cell formalisms have now established that the magnetosphere of an oblique rotator entails a current sheet outside its light cylinder whose rotating distribution pattern moves with linear speeds exceeding the speed of light in vacuum. Here we insert the description of the current sheet provided by the numerical simulations in the classical expression for the retarded potential and thereby calculate the radiation field generated by this source in the time domain. We find a radiation consisting of highly focused pulses whose (i) spectrum can extend from radio waves to gamma rays, (ii) brightness temperature can exceed 10^(40) K, (iii) linear polarization can be 100%, (iv) two concurrent polarization position angles are approximately orthogonal often and swing through 180 deg across the pulse profile in most cases, (v) circular polarization reverses sense across some components of the pulse profile, (vi) microstructure is determined by the thickness of the current sheet, and (vii) whose flux density diminishes with the distance D from the star as D^(-3/2) (rather than D^(-2)) in certain directions. The intrinsically transient radiation process analysed here is thus capable of generating an emission whose features are strikingly similar to those of the emissions received from pulsars and magnetars and from the sources of fast radio bursts and gamma-ray bursts.


Main
From the results obtained by numerical simulations of the magnetospheric structure of an obliquely rotating neutron star, [3][4][5] Tchekhovskoy et al. 7 have derived a semi-analytic de-scription of the distributions of the electric and magnetic fields that permeate the plasma surrounding these objects (see also ref. 14). The described fields in conjunction with Maxwell's equations provide an explicit expression for the space-time distribution of the density of magnetospheric charges and currents including that of the current sheet (see §1 of the Supplementary Information). The surface on which the current sheet is distributed spirals away from the light cylinder in the azimuthal direction at the same time as undulating in the latitudinal direction (see Fig. 1a). Its motion consists of a rotation with the angular frequency of rotation of the central neutron star, ω, and a radial expansion with the speed of light in vacuum, c. This is not incompatible with the requirements of special relativity because the superluminally moving distribution pattern of the current sheet is created by the coordinated motion of aggregates of subluminally moving charged particles. [15][16][17] In this paper we treat the distribution of charges and currents that make up the current sheet at any given time as a prescribed volume source whose density can be inserted in the retarded solution of the inhomogeneous Maxwell's equations to find the radiation field it generates in unbounded free space. The only role we assign to the rest of the magnetosphere, whose radiation field is negligibly weaker than that of the current sheet, is to maintain the propagation of this sheet. The multi-wavelength focused pulses emitted by the current sheet escape the plasma surrounding the neutron star in the same way that the radiation generated by the accelerating charged particles invoked in most current attempts at modelling the emission mechanism of these objects does. 6,18 The current sheet is described by charge and current densities whose space-time distributions depend on the azimuthal coordinate ϕ and time t in the combination ϕ − ωt only. The radiation field we are after can be built up, therefore, by the superposition of the fields of the uniformly rotating volume elements that constitute this source. Superluminal counterpart of the field of synchrotron radiation, which plays the role of such a Green's function for the present problem, entails intersecting wave fronts that possess a two-sheeted cusped envelope (see Figs. 1b, 1c and 1d). Outside the envelope only one wave front passes through the observation point at any given observation time; but inside the envelope three distinct wave fronts, emitted at three different values of the retarded time, simultaneously pass through each observation point. Coalescence of two of the contributing retarded times on the envelope of wave fronts results in the divergence of the Green's function on this surface. At an observation point on the cusp locus of the envelope all three of the contributing retarded times coalesce and the Green's function has a higher-order singularity (see §2 of the Supplementary Information).
Constructive interference of the emitted waves and formation of caustics thus play a crucial role in determining the radiation field of the current sheet. Not only the integral that defines the Green's function for the problem but every one of the repeated integrals in the classical expression for the retarded potential in the present case entails either singularities or nearby saddle points that coalesce and thereby result in further focusing of the radiation at certain observation points (see §3 of the Supplementary Information).

Results
The results reported below are consequences mainly of the shape and motion of the current sheet: two features of the magnetosphere which are the same not only both for a dipolar and a monopolar magnetic field at the surface of the star, but also both close to and far from the light cylinder. 7, 14

Pulse profiles and polarization position angles
Two examples of the longitudinal distributions of the Stokes parameters for the radiation generated by the current sheet are shown in Fig. 2. The total intensity I and the intensities of the linearly and circularly polarized components of the radiation, L and V , are plotted for given values of the inclination angle α and colatitude θ P of a far-field observer in units of where B 0 is the magnitude of the star's dipolar field at its magnetic pole, r s0 is the star's radius and R P is the distance of the observer from the star. In these examples, and in Fig. 3, the origin of the azimuthal coordinate of the observation point is shifted to place the pulse window next to longitude zero. The electric field of the radiation turns out to be the sum of two distinct parts with differing polarization position angles (see §4 of the Supplementary Information). Figure 2 also shows the longitudinal distributions of the position angles of these two parts which we refer to as P and Q modes. Here, those branches of the multi-valued function defining the position angle are adopted that yield continuous position-angle distributions across various components of a given pulse. The variable k −1 u , whose value determines the wavelength of the small-amplitude modulations (microstructure) of the distributions shown in Fig. 2, represents a lower limit to the thickness of the current sheet in units of the light-cylinder radius. The thickness assigned to the current sheet by the description in ref. 7 is zero. However, a superluminally moving source is necessarily volume-distributed: 17 it would give rise to a divergent field if it has no thickness. We have circumvented this shortcoming of the description given in ref. 7 by replacing the infinitely long integration range in the Fourier representation of the Dirac delta function describing the current sheet by the truncated wave-number interval |k| ≤ k u (see Method). Figure 3 illustrates an example of a radically different type of pulse: one detectable near those observation points for which two nearby saddle points of the integral over the latitudinal distribution of the current sheet coalesce, thus giving rise to a much tighter focusing of the emitted waves. Though their profiles over the entire pulse window look similar to those of other pulses (as in Fig. 3a), such pulses display extraordinarily large amplitudes and short widths once their peaks are resolved (as in Fig. 3b). We shall see below that the extraordinary values of the amplitudes and widths of such pulses, illustrated by the example in Fig. 3, are what underpin the high brightness temperatures and broad frequency spectra of the radiation generated by the current sheet.
Further examples of pulse profiles and position-angle distributions can be found in the Supplementary Information. It should be added that at any given value of the inclination angle α, the pulse observed at π − θ P differs from that observed at θ P only in that the intensity V of its circularly polarized part is replaced by −V and its longitude ϕ P is replaced by ϕ P + π. Moreover, the results for α > π/2 follow from those for α < π/2 by replacing θ P , ϕ P and V by π − θ P , ϕ P + π and −V , respectively.

Brightness temperature
By equating the magnitude of the Poynting flux of the radiation to the Rayleigh-Jeans law for the energy that a black body of the same temperature would emit per unit time per unit area into the frequency band ∆ν centred on the frequency ν, it can be shown that the brightness temperature of the radiation T b is related to the Stokes parameter I by (2) in which I is in units of I 0 ,B 0 is the value of the magnetic field B 0 in units of 10 12 Gauss, d is the value of the star's radius r s0 in units of 10 6 cm, D is the value of the observer's distance R P in kpc,ν and ∆ν are the values of the radiation frequency and its bandwidth in units of 10 8 Hz and 10 6 Hz, respectively, andP = 10 2 /ω is the value of the star's rotation period in seconds. Equation

Frequency spectrum
Given that the radiation field of the current sheet depends on the observation time t P and the azimuthal coordinate ϕ P of the observation point only in the combination ϕ P −ωt P , the frequency spectrum of this radiation is equally well described by the Fourier decomposition of its longitudinal distribution. In the present case, the content of this spectrum stems from two factors. One factor is the thickness of the current sheet manifested in the sharp small-amplitude modulations of the pulse profile whose wavelengths are proportional to k −1 u (see Fig. 1a). The other factor is the full width at half maximum (δϕ P ) of the component of the pulse profile with the highest peak (see Fig. 3b). While the spectral distribution of the former lies in the radio band when k u 10 5 , that of the latter ranges from radio waves to gamma rays: the width δϕ P = 1.21 × 10 −26 radian of the pulse depicted in Fig. 3b, for example, implies a frequency spectrum that extends as far as ν ω/(2πδϕ P ) 1.31 × 10 27P −1 Hz.
In the lower frequency range, the Stokes parameter I has the power-law dependence ν −β on frequency with a spectral index β that assumes the following values in various regimes: 2/3, 1, 4/3, 5/3, 2, 7/3 (see §4 of the Supplementary Information).

Flux density and its rate of decay with distance
In the present case where the length subtended by the longitudinal extent of the radiation beam embodying a high-frequency pulse can be shorter than 1 cm at D = 1 kpc, the flux density S of the radiation is related to the Stokes parameter I by in which I is in units of I 0 and δϕ P in radians. In the case depicted in Fig. 3 where α = 60 • , θ P = 90 • and k u = 10 7 , for example, the flux density S has the value 32.1 S 0 erg/(sec×cm 2 ) at D = 1 kpc. The separation between two nearby saddle points of the integral over the latitudinal distribution of the current sheet decreases as R −1/2 P with increasing distance R P of the observation point from the star. The enhanced focusing of the radiation that is caused by this shortening of the separation between the saddle points gives rise to both a narrowing of the width and an augmenting of the peak intensity of the emitted pulse. Because the peak intensity and the width of the pulse are modified with different rates, this effect results in a flux density that diminishes with increasing distance as R −3/2 P rather than R −2 P (or equivalently as D −3/2 instead of D −2 ). The latitudinal width δθ P of such nonspherically decaying pulses is of the order of (R P ω/c) −1 radians in most cases. But there are, in general, several latitudes near which such pulses can be observed; the number and locations of these latitudes are determined by the values of α and R P (see §4 of the Supplementary Information).
The violation of the inverse-square law encountered here is not incompatible with the requirements of the conservation of energy because the radiation process discussed in this paper is intrinsically transitive. The difference in the fluxes of power across any two spheres centred on the star is in this case balanced by the change with time of the energy contained inside the shell bounded by those spheres (see Appendix C of ref. 9 for a demonstration of this feature).
Given their limited latitudinal extent, the non-spherically decaying pulses generated by the current sheet of a neutron star are more likely to be observed when, as a result of the precession of the star's rotation axis, the radiation beams embodying such pulses sweep past the Earth. Using the decay rate R −2 P , we would over-estimate the power emitted by the sources of the bursts of radiation we would receive in this way by as large a factor as 10 14 if the neutron stars that generate the bursts lie at cosmological distances. The enormously high energetic requirements normally attributed to the sources of fast radio bursts and gamma-ray bursts 12, 13 could therefore be artefacts of the invariably made assumption that the radiation fields of all sources necessarily decay as predicted by the inverse-square law.

Method
Both the outcomes of the numerical simulations of the magnetosphere of an oblique rotator 3-6 and their semi-analytic description 7 have appeared in the literature as space-time distributions of the electric and magnetic fields. Once inserted in Maxwell's equations, the field distributions described in ref. 7 yield the following expressions for the charge and current densities, ρ (cs) and j (cs) , of the current sheet in terms of the spherical polar coordinates (r s , θ, ϕ) and time t: and where C = sin α sin θ cos (φ +r s −r s0 ) + cos α cos θ, h = sin α cos θ cos (φ +r s −r s0 ) − cos α sin θ, δ is the Dirac delta function,r s = rω/c,r s0 = r s0 ω/c and (ê rs ,ê θ ,ê ϕ ) are the base vectors of the coordinate system (see §1 of the Supplementary Information). For the purposes of the present analysis, it is essential that the finiteness of the duration of the source is taken into account (see Appendix B of ref. 9). If ρ (cs) and j (cs) are turned on at t = 0, then the coordinates t and ϕ in equations (4) and (5) both range over (0, ∞) but the values of the combinationφ = ϕ − ωt in which they occur has a limited range of length 2π, e.g., 0 ≤φ < 2π.
As can be seen from the alternative form ϕ =φ + ωt of equation (7),φ is a Lagrangian coordinate that labels the rotating volume elements of the current distribution on each circle r s =const, θ =const, by their azimuthal positions at the time t = 0. This coordinate cannot range over a wider interval because the aggregate of volume elements that constitute a rotating source in its entirety can at most occupy an azimuthal interval of length 2π at any given time (e.g., at t = 0).  Figure 1: a, Snapshot of a single turn of the current sheet about the light cylinder (the cylinder at which the speed of any entity co-rotating with the neutron star would equal the speed of light c). This surface undulates within the latitudinal interval (π/2 − α, π/2 + α) each time it turns about the rotation axis, thus wrapping itself around the light cylinder as it extends to the outer edge of the magnetosphere (α which denotes the angle between the magnetic and rotation axes of the star has the value π/3 in this figure). b, Cross sections of the wave fronts (the circles in light blue) emanating from a volume element of the current sheet with fixed radial and latitudinal coordinates. This figure is plotted for a source element the radius of whose orbit (the dotted circle) is 2.5 times the radius c/ω of the light cylinder (the green circle). Cross sections of the two sheets of the envelope of these wave fronts with the plane of the orbit (shown in dark blue) meet at a cusp and wind around the rotation axis, while moving away from it all the way to the far zone. c, Three-dimensional plot of the envelope of wave fronts emanating from a superluminally rotating source element and d, the cusp of this envelope along which the two sheets of the envelope meet tangentially. The cusp touches the light cylinder where it crosses the plane of the orbit and moves away from the axis of rotation as it and the envelope itself spiral into the far zone maintaining a symmetry with respect to this plane.  Figure 2: a, Longitudinal distributions of the intensities of the total (I) and circularly polarized (V ) radiation generated by the current sheet in units of I 0 for the inclination angle α = 5 • at a far-field observation point with the colatitude θ P = 5 • . Distribution of the intensity of the linearly polarized radiation (L) is essentially coincident with that of I in this case. The sharp modulations of the pulse profile (its microstructure) are discernible here because this figure is plotted for k u = 10 2 which corresponds to a relatively thick current sheet. The thinner the current sheet, the shorter are the wavelengths of these modulations. At its peak, the right-hand component of this pulse has the intensity I = 1.05 × 10 12 I 0 and the longitudinal width 6.76 × 10 −9 second when k u is 10 7 , i.e., when the thickness of the current sheet is of the order of 10 −7 c/ω. b, Position angles of the two polarization modes versus longitude for the pulse depicted in a. Note that not only does V reverse sign across the left-hand component of this pulse but also the position angle of the Q mode swings through 180 • across the profile of this component. Note also the approximate orthogonality of the two modes over most of the pulse window. c, The Stokes parameters I, L and V versus longitude for α = 80 • , θ P = 110 • and k u = 10 4 . This is an example of a pulse with several widely separated peaks for which the Stokes parameters are comparable in magnitude over some longitudinal intervals. Though hardly visible, small amplitude sharp modulations pervade also the distribution here. d, Position angles of the two polarization modes versus longitude for the pulse depicted in c exemplifying position-angle distributions that undergo differing variations across different components of the pulse profile. Part d is plotted with k u = 10 7 for clarity.  Figure 3: a, Distributions of the Stokes parameters I, L and V in units of I 0 for the inclination angle α = 60 • , an observation point in the equatorial plane θ P = 90 • , a distance of 10 13 light-cylinder radii and k u = 10 7 . In this example, the colatitude of the observation point lies close to one of the critical colatitudes for the given values of the inclination angle and distance (the colatitudes at which two nearby saddle points of the integral representing the retarded potential coalesce). b, The right-hand component of the pulse depicted in a is here plotted over a sufficiently short longitudinal interval to resolve its peak and width. The values of α, θ P , k u and distance in b are the same as in a but the origin of longitude is shifted in b for clarity. The shape of the pulse depicted in b is the same in all cases of this type.
As pointed out in the text, we have replaced the range of integration in the following Fourier representation of the Dirac delta function that appears in equations (4) and (5), by the finite interval (−k u , k u ) to circumvent the divergence that arises from the vanishing of the thickness of the current sheet. In this paper, the retarded potential that arises from the charge and current densities described by the modified versions of equations (4) and (5) is inserted in to obtain the corresponding expression for the generated radiation field (E, B), where (x, t) and (x P , t P ) are the space-time coordinates of the source points and the observation point P , respectively, and R is the magnitude of the separation R = x P − x (see, e.g., ref. 8).
To satisfy the required boundary conditions at infinity the free-space radiation field of an accelerated superluminal source has to be calculated (in the Lorenz gauge) by means of the retarded solution of the wave equation for the electromagnetic potential. There is a fundamental difference between the classical expression for the retarded potential and the corresponding retarded solution of the wave equation that governs the electromagnetic field. While the boundary contribution to the retarded solution of the wave equation for the potential that appears in Kirchhoff's surface-integral representation can always be rendered equal to zero by means of a gauge transformation that preserves the Lorenz condition, the corresponding boundary contribution to the retarded solution of the wave equation (or any other equation) for the field cannot be assumed to be zero a priori. Not to exclude emissions whose intensity could decay more slowly than predicted by the inverse-square law, it is essential that the radiation field is derived from the retarded potential (see § 3 of ref. 9 where this point is expounded).
Since the problem discussed in this paper entails the formation of caustics we cannot proceed to the far-field limit |x P | → ∞ before evaluating the radiation field, as is customarily done in radiation theory. The far-field approximation of the argument of the delta function in (13) would replace spherical wave fronts by planar wave fronts thereby relinquishing the possibility of their constructive interference. Moreover, given the exceptionally short scales of the longitudinal (or equivalently temporal) variations of the present radiation, it would be intractably more difficult to obtain the time-domain results reported in this paper by means of a frequency-domain analysis (see, e.g., ref. 19). Further technical reasons why a conventional approach to the present problem does not work are discussed in Appendix B of ref. 9.
We calculate the retarded potential (13) in §3 of the Supplementary Information by first performing the integration with respect to t at a fixed (r s , θ,φ), i.e., by superposing the contributions of the uniformly rotating volume elements of the current sheet illustrated in Figs. 1b-d. A uniform asymptotic approximation to the value of the integral over t, which entails two nearby saddle points, 20 is obtained by the time-domain version of the method of Chester, Friedman and Ursell. 21,22 The divergences of the resulting expression (which acts as the Green's function for the problem) on the envelope of the emitted waves and its cusp stem from the relativistic restrictions inherent in Maxwell's equations and reflect the fact that no superluminal source can be point-like. 17 Once the product of the charge-current density with this Green's function is integrated over theφ extent of the current sheet with the aid of Hadamard's regularisation technique, 23,24 an expression is obtained that is singular only on the hyperboloid swept by the cusp loci of the envelopes of various source elements (see Fig. 1d). This remaining singularity is also removed when the integration with respect to the radial extent of the current sheet is performed: an integration that receives its main contribution from the intersection of the hyperboloid in question with the current sheet. The last integral (that with respect to θ) is singularity-free and once again entails two nearby saddle points that coalesce for certain values of the colatitude θ P of the observation point.
The five-dimensional integration (with respect to t, r s , θ,φ and k) required for evaluating the radiation field (14) has here been carried out analytically. The only assumption made in the analysis presented in the Supplementary Information is that the radiation frequency appreciably exceeds the rotation frequency ω/2π.
A final remark is in order: it is often presumed that the plasma equations used in the numerical simulations of the magnetospheric structure of an oblique rotator should, at the same time, predict any radiation that the resulting structure would be capable of emitting. 3,4 Irrespective of the formalism on which they are based (whether force-free, MHD or particle-in-cell), the plasma equations used in these simulations are formulated in terms of the electric and magnetic fields (as opposed to potentials). It has already been demonstrated in §3 of ref. 9, however, that the gauge freedom offered by the solution of Maxwell's equations in terms of potentials plays an indispensable role in the prediction of the characteristics of the present radiation. The absence of high-frequency radiation (and, specifically, the type of radiation described in this paper) is hardwired into the numerical simulations that have been performed to determine the magnetospheric structure of an oblique rotator by the imposition of the standard boundary conditions on the fields in the far zone (see §3 of ref. 9).

Supplementary Information
Supplementary Information for this article is included below in the form of a paper, a paper that presents the mathematical derivations of the results reported in the present article.
Supplementary Information for "Radiation by the superluminally moving current sheet in the magnetosphere of a neutron star" 1 Semi-analytic description of the magnetosphere In a reference frame marked by the spherical polar coordinates (r s , θ, ϕ) for which θ = 0 coincides with the axis of rotation, the distributions of the electric and magnetic fields in the magnetosphere of an oblique rotator are, according to [1], given by where t (≥ 0) is time, ω is the angular frequency of rotation of the star, c is the speed of light in vacuum, α is the angle between the rotation and magnetic axes of the star (henceforth referred to as the inclination angle), r s0 is the radius of the star and B 0 is the magnitude of the star's dipolar field at its magnetic pole. The caret on r s and r s0 (and any other variable with the dimension of a length) is used in this paper to designate a variable that is rendered dimensionless by being measured in units of the light-cylinder radius c/ω.
Inserting these fields in Maxwell's equations, we obtain the corresponding distributions of the charge density ρ and the current density j in the magnetosphere: the base vectors (ê rs ,ê θ ,ê ϕ ) are those of the spherical coordinate system (r s , θ, ϕ) and δ is the Dirac delta function. In deriving (1.12), we have made use of the fact that here E depends on ϕ and t only in the combination ϕ−ωt and r 2 s B rs depends on ϕ andr s only in the combination ϕ +r s , i.e., that ∂E/∂t = −ω∂E/∂ϕ and ∂(r 2 s B rs )/∂r s = (ω/c)∂(r 2 s B rs )/∂ϕ. The magnetospheric current sheet shown in figure 1.1 occurs where the argument of the delta function in (1.12), i.e., the function C defined by (1.6), vanishes. Since the radiation whose frequency appreciably exceeds the rotation frequency ω/2π is due entirely to this current sheet, we base our analysis of the radiation field in the following sections only on those terms in the above expressions for ρ and j that involve the Dirac delta function. We will see in § 4 that the radiation generated by this sheet is in addition considerably more intense than that generated by the rest of the magnetosphere. Disregarding the other terms in (1.11) and (1.12), we obtain and for this sheet's charge and current densities, wherê and When the inclination angle α lies in the interval [0, π/2], the above source densities are non-zero only in as can be seen from the argument of the delta function that appears in (1.18) and (1.19). It turns out that the symmetry of the current sheet with respect to θ and α enables us to infer the results for π/2 < α < π from those for 0 < α < π/2 (see § 3.4).
For the purposes of quoting the results derived in [2], it is convenient to express the base vectors of the spherical polar coordinates used above in terms of the base vectors (ê r ,ê ϕ ,ê z ) of the cylindrical polar coordinates (r, ϕ, z) defined by i.e., to replace (ê rs ,ê θ ,ê ϕ ) viâ e rs = sin θê r + cos θê z ,ê θ = cos θê r − sin θê z .
Equations (1.18) and (1.19) can then be written as describe the dimensionless charge densityρ and the cylindrical components of the dimensionless current densityj of the magnetospheric current sheet as functions of (r s , θ,φ) (see (1.10) and (1.20)).

Formulation of the problem 2.1 The free-space solution of Maxwell's equations that satisfies the required boundary conditions at infinity
To satisfy the required boundary conditions at infinity the free-space radiation field of an accelerated superluminal source has to be calculated (in the Lorenz gauge) by means of the retarded solution of the wave equation for the electromagnetic potential. There is a fundamental difference between the classical expression for the retarded potential and the corresponding retarded solution of the wave equation that governs the electromagnetic field. While the boundary contribution to the retarded solution of the wave equation for the potential that appears in Kirchhoff's surface-integral representation can always be rendered equal to zero by means of a gauge transformation that preserves the Lorenz condition, the corresponding boundary contribution to the retarded solution of the wave equation (or any other equation) for the field cannot be assumed to be zero a priori. Not to exclude emissions whose intensity could decay more slowly than predicted by the inverse-square law, it is essential that the radiation field is derived from the retarded potential (see § 3 of [2] where this point is expounded). Accordingly, we base the analysis in this paper on the classical expression for the retarded potentials that arise from the charge and current densities described by (1.26) and insert this in to obtain the corresponding expression for the generated fields: where (x, t) and (x P , t P ) are the space-time coordinates of the source points and the observation point P , respectively, R is the magnitude of the separation R = x P − x,n = R/R is a unit vector along the radiation direction and δ denotes the derivative of the Dirac delta function with respect to its argument (see, e.g., [3]). Since the problem we will be analysing entails the formation of caustics we cannot proceed to the far-field limit |x P | → ∞ before evaluating the above integral, as is customarily done in radiation theory. The far-field approximation of the argument of the delta function in (2.3) would replace spherical wave fronts by planar wave fronts thereby relinquishing the possibility of their constructive interference. As we will be applying our results to astronomical objects we can, however, approximate the unit vectorn by its far-field valuê at this stage, so that the magnetic field can be obtained from B =n ∞ ×E once the electric field is known. In (2.4), (ê r P ,ê ϕ P ,ê z P ) are the cylindrical base vectors at the observation point P .
For the purposes of the present analysis, it is essential that the finiteness of the duration of the source is taken into account (see appendix B of [2]). If ρ (cs) and j (cs) are turned on at t = 0, then the coordinates t and ϕ in (1.26) both range over (0, ∞) but the values of the combinationφ = ϕ − ωt in which they occur has a limited range of length 2π, e.g., 0 ≤φ < 2π. (2.5) As can be seen from the alternative form ϕ =φ+ωt of (1.20),φ is a Lagrangian coordinate that labels the rotating volume elements of the current distribution on each circle r =const, z =const, by their azimuthal positions at the time t = 0. This coordinate cannot range over a wider interval because the aggregate of volume elements that constitute a rotating source in its entirety can at most occupy an azimuthal interval of length 2π at any given time (e.g., at t = 0). In this section, we mark the spatial coordinates of the source points by cylindrical polar coordinates and eliminate t in favour ofφ. Thus changing the variables of integration in (2.3) to (x, t) = (r, ϕ, z,φ) and introducing the dimensionless coordinatesr = rω/c andẑ = zω/c, we obtain the function g(r, ϕ,ẑ;r P , ϕ P ,ẑ P ) is defined by g ≡ ϕ − ϕ P +R, (2.8) and the variable φ in the argument of the delta function stands for We have expressed the range of ϕ integration as a sum of the intervals of length 2π that the element initially located atφ traverses during each of its individual rotations: m is a positive integer enumerating successive rotation periods (the first rotation period being designated by m = 1) and the summation extends over the set of rotations executed by the source over its lifetime.

The Green's function for the problem and its loci of singularities
To put the current density j (cs) = j (cs) zêz into a form suitable for performing the integration with respect to ϕ, we need to express the ϕ-dependent base vectors (ê r ,ê ϕ ,ê z ) associated with the source point (r, ϕ, z) in terms of the constant base vectors (ê r P ,ê ϕ P ,ê z P ) at the observation point ( Once the resulting expression, is inserted in (2.6) and δ (g − φ) is written as −∂δ(g − φ)/∂φ (see (2.9)), we arrive at denotes the outcome of the remaining integration with respect to ϕ and ). Note that sinceφ + 2(m − 1)π designates the same source element aŝ ϕ + 2mπ the dependence onφ of the limits of integration in (2.13) does not contribute toward the values of the derivatives of G j with respect toφ. The function G j (r,φ,ẑ;r P ,φ P ,ẑ P ) here acts as the Green's function for the present problem. It describes the Liénard-Wiechert field that arises from an individual volume element of the rotating distribution pattern of the source. If we specialize the current distribution to a rotating point charge q, i.e., let j r = j z = 0 and j ϕ = r ωqδ(r −r )δ(φ)δ(z) with a constant r , then (2.13) at an observation point in the far zone would describe the familiar field of synchrotron radiation when r < c/ω and a synergic field combining attributes of both synchrotron andČerenkov emissions when r > c/ω.
Depending on the value of for a given source point (r,φ, z) with rω > c, the ϕ-dependence of the function g that appears in the definition of the Green's function G j in (2.13) has one of the generic forms shown in figure 2.1. As can be seen from the curve labelled ∆ > 0 in this figure, there are values, of the retarded position of the source point at which vanishes and so G j diverges. These turning points of g occur at source points for which ∂(R| ϕ=φ+ωt )/∂t = −c, i.e., the source points that approach the observer, along the radiation directionn, with the speed of light at the retarded time. The inflection point of g (see the curve labelled ∆ = 0 in figure 2.1), at which in addition vanishes, occurs at source points that approach the observer not only with the wave speed but also with zero acceleration at the retarded time, i.e., for which both ∂(R| ϕ=φ+ωt )/∂t = −c and ∂ 2 (R| ϕ=φ+ωt )/∂t 2 = 0 at the time when g| ϕ=φ+ωt = φ and ∂g/∂ϕ = ∂ 2 g/∂ϕ 2 = 0. In (2.21), is the value ofR at the extrema ϕ ± of g.
The envelope of the wave fronts emanating from a given rotating source element (r,φ,ẑ), on which ∂g/∂ϕ vanishes, consists of the rigidly rotating two-sheeted surfaceφ −φ P = g(ϕ ± ) in the space (r P ,φ P ,ẑ P ) of observation points. This surface, which is shown in figures. 2.2 and 2.3, is described by (2.19) and (2.22)). The two sheets of this surface tangentially meet along a cusp on which ∂ 2 g/∂ϕ 2 as well as ∂g/∂ϕ vanishes (see figures 2.3 and 2.4). Three distinct wave fronts, emitted at three differing values of the retarded time, pass through any given observation point inside the envelope. At an observation point located on the envelope or its cusp, respectively two or all three of these waves coalesce and interfere constructively (see figure 2.2).

Bifurcation surface of an observation point
Reciprocally, the locus in the space of source points (r,φ,ẑ) on which ∂g/∂ϕ vanishes is a two-sheeted cusped surface issuing from the fixed observation point P (see figure 2.5). We refer to this locus, which is described by (2.23) for fixed values of (r P ,φ P ,ẑ P ) rather than fixed values of (r,φ,ẑ), as the bifurcation surface of the observation point P . The two sheets φ = φ + and φ = φ − of this surface meet along the following cusp: (2.24) where m is the same integer as that appearing in (2.13). In this paper we refer to both C and its projection onto the (r s , θ) plane as the cusp locus of the bifurcation surface; whether it is C itself or its projection that is referred to will be clear from the context. The source points inside the bifurcation surface, close to its cusp, make their contributions toward the observed value of the field at three distinct retarded positions in their trajectory (where a horizontal line g = φ in figure 2.1 intersects the curve ∆ > 0 between its extrema), while those outside the bifurcation surface make their contributions at a single retarded position (where the curve ∆ < 0 is intersected by g = φ in figure 2.1). For the source points on the bifurcation surface (i.e., those for which g = φ ± in figure 2.1), two of the contributing retarded positions coalesce at the extrema of the curve ∆ > 0 in figure 2.1 1: Generic forms of the function g(ϕ) for source points whose (r,ẑ) coordinates lie across the boundary ∆ = 0 delineating the projection of the cusp curve of the bifurcation surface onto the (r,ẑ) plane. Depending on whether φ lies outside or inside the interval (φ − , φ + ), contributions are made toward the observed field (i.e., the argument g(ϕ) − φ of the Dirac delta function in (2.13) vanishes) at either one or three retarded positions of the source. For a horizontal line g = φ that either approaches an extremum of g(ϕ) from inside the interval (φ − , φ + ) or passes through an inflection point of g(ϕ), two or all three of the retarded positions in question coalesce and so their contributions interfere constructively to form caustics. This figure is forr = 3 and only shows two rotation periods. At higher speeds, the difference between the values of φ + and φ − can be large enough for a horizontal line g = φ to intersect g(ϕ) over more than one rotation period (see figure 36 in [2]). Contributions toward the observed field can thus arise, not only from one or three, but from any odd number of retarded positions of the source. There are contributions from more than three retarded times whenever the rotation period of the source is shorter than the time taken by the collapsing sphere |x − x P | = c(t − t P ), centred on the observation point P , to cross the orbit of the source. Three-dimensional view (in the space (r P ,φ P ,ẑ P ) of observation points) of the envelope of wave fronts emanating from the rotating source point (r,φ,ẑ). This envelope consists of two sheets that tangentially meet along a cusp (see figure 2.4). The singular sheet, i.e., the sheet that issues from the source point with an initial conical shape, is that described byφ P =φ − φ − (r P ,ẑ P ;r,ẑ). The cusp along which the two sheets of the envelope of wave fronts meet and are tangent to one another. This cusp touches and is tangent to the light cylinderr P = 1 on the planeẑ P =ẑ and spirals outward into the far field on the hyperboloid generated by the revolution of the curve ∆(r P ,ẑ P ;r,ẑ) = 0.
giving rise to a divergent value of the Green's function at P (see figures 9 and 10 of [2]). For the source points located on the cusp locus C of the bifurcation surface (i.e., those for which ∆ = 0 in figure 2.1), all three of the contributing retarded positions coalesce at the inflection point of the curve ∆ = 0 in figure 2.1 giving rise to a higher-order singularity in G j .

A uniform asymptotic approximation to the value of the Green's function near the cusp locus of the bifurcation surface
The time-domain version [4] of the method of Chester et al. [5] can be employed to derive a uniform asymptotic approximation to the value of G j for the source points close to the cusp C of the bifurcation surface. The result, which corresponds to that derived in §4.5 of [2] for the case n = 1, is The two sheets φ = φ ± of the bifurcation surface issuing from the observation point P , the cusp C of this surface and the light cylinderr = 1. In contrast to the envelope of wave fronts which resides in the space of observation points, the surface shown here resides in the space (r,φ, z) of source points: it is the locus of source points that approach P , along the radiation direction, with the speed of light at the retarded time. The two sheets of this surface meet along a cusp that tangentially touches the light cylinder atẑ =ẑ P and moves outward spiralling around the rotation axis on the hyperboloid generated by the revolution of the curve ∆(r,ẑ;r P ,ẑ P ) = 0 (see figure 2.4). The source points on this cusp approach the observer along the radiation direction not only with the speed of light but also with zero acceleration at the retarded time. If the position of the observation point is such that the cusp shown here intersects the current sheet, there will be wave fronts with differing emission times that are received simultaneously: while the source points outside the bifurcation surface make their contributions toward the value of the observed field at a single instant of retarded time, the source points inside this surface make their contributions at 3 (or 5, 7, · · ·) distinct instants of retarded time.
and H(x) denotes the Heaviside step function. The two-dimensional loci χ = ±1 across which the above expression for the Green's function G j changes form correspond to the two sheets φ ± of the bifurcation surface, respectively. As a source point (r,φ, z) in the vicinity of the cusp C approaches the bifurcation surface from inside, i.e., as χ → 1− or χ → −1+, G in j diverges. However, as a source point approaches either one of the sheets of the bifurcation surface from outside, the numerator and the denominator in (2.27) vanish simultaneously and G out j tends to a finite limit: Note that c 1 , and hence p j and q j , are independent of m (see (2.19), (2.23) and (2.30)). Thus the Green's function G j is singular only on the inner side of the bifurcation surface (see figures 9 and 10 of [2]).

Hadamard's finite part of the divergent integral representing the field
It follows from (2.25) and (2.26) that the factor ∂G j /∂φ in the integrand of the integral (2.12) diverges as (1 − χ 2 ) −3/2 and so has a non-integrable singularity on the bifurcation surface where χ 2 equals 1. This singularity has arisen because we differentiated the retarded potential (2.1) under the integral sign when calculating the field. Had we evaluated the integral in (2.1) prior to differentiating it we would have found a singularity-free expression.
Interchanging the orders of integration and differentiation is mathematically permissible when the integrand is discontinuous only if one treats the resulting integral as a generalized function and so one handles any non-integrable singularities that consequently arise by means of Hadamard's regularization technique (see [6], [7] and the illustrative example in appendix A of [2]). Hadamard's procedure consists of performing an integration by parts and discarding the divergent (integrated) term in the resulting expression. The remaining finite part is the value that Hadamard's regularization assigns to the integral; in the present case, it is the value we would have obtained if we had first evaluated the finite integral representing the retarded potential and had differentiated the result [Φ(x P , t P ), A(x P , t P )] of that evaluation subsequently. (The more direct approach, in which the potential is first evaluated and then differentiated, cannot of course be carried out for any realistic source distribution analytically.) Theφ coordinatesφ ± of the two sheets of the bifurcation surface depend on the observation time t P [see (2.31) and (2.9)], so that these two sheets move across theφ extent of the source distribution as t P elapses. If the position of the observation point is such that the cusp locus of the bifurcation surface intersects the source distribution, the two sheets of this surface (which tangentially meet at the cusp) will divide the volume of the source into a part that lies inside and a part that lies outside the bifurcation surface. The Lagrangian coordinatesφ designating the initial azimuthal positions of the constituent volume elements of a source that fully occupies an annular region range over the interval 0 ≤φ < 2π. The (r,ẑ) coordinates of these source elements either fall in ∆ ≥ 0 or in ∆ < 0. The elements in ∆ ≥ 0 are always divided into two sets: a set inside the bifurcation surface for whichφ − ≤φ ≤φ + and so the Green's function G j has the form G in j and a set outside for whichφ lies either in (0,φ − ) or in (φ + , 2π) and so G j has the form G out for all values of (r,ẑ) within the magnetosphere, then the source lies entirely outside the bifurcation surface and G j has the form G sub j . Note that, for certain space-time coordinates of the observation point P , the values ofφ − andφ + that lie in the interval (0, 2π) could correspond to different rotation periods, i.e., to different values of m [see (2.19), (2.22) and (2.23)]. To simplify the notation, here we adopt an observation time t P at which the values ofφ − andφ + that lie in the interval (0, 2π) correspond to the same rotation period m.
Breaking up the volume of integration in the expression for the radiation field E into the domains of validity of G in j , G out j and G sub j , we can therefore write theφ-integral over u j in (2.12) as If we now integrate every term of the above expression by parts, recall thatφ = 0 labels the same source point as doesφ = 2π, and use the fact that the exact version of G j given in (2.13) is periodic inφ as well as in ϕ (with the same period 2π), we arrive at an expression that reduces to once the integrals over G in j , G out j and G sub j are combined in the light of (2.25). We have seen in the last paragraph of § 2.4 that G in nj diverges atφ =φ ± . The Hadamard finite part of Iφ is therefore given by the right-hand side of (2.38) without the divergent terms involving G in j |φ =φ − and G in j |φ =φ + : where Fp{Iφ} denotes the Hadamard finite part of the divergent integral Iφ (see [6,7]).
Once the integral with respect toφ in (2.12) is equated to the expression on the righthand side of (2.39), we find that and The term E v constitutes the contribution from the entire volume of the source while the terms E b ± denote the contributions from the discontinuities of the Green's function on the two sheets φ = φ ± of the bifurcation surface, respectively. We will see that the terms E b ± describe an unconventional radiation field with characteristics that turn out to differ from any previously known radiation fields.
3 Radiation field of the current sheet

The contribution from discontinuities of the Green's function
Since the description of the current sheet given in § 1 is in terms of spherical polar coordinates, it is more convenient, for the purposes of evaluating the contributions E b ± to the radiation field E, to change the integration variables in (2.42) from (r,ẑ) to (r s , θ) while continuing to designate the orientations of the current density and the field by means of the cylindrical base vectors (ê r P ,ê ϕ P ,ê z P ). Equations (1.27)-(1.30), (2.14), (2.35) and (2.42) then jointly yield in which ∆, c 1 , p j , q j andφ ± are expressed as functions of (r s , θ) by the insertion of ( i.e., whenφ assumes one of the following valueŝ ϕ l = (−1) l arccos(cot α cot θ) −r s +r s0 + (2n + 1)π, (3.3) where l is either 1 or 2 and n is the integer for which the requirement 0 ≤φ l ≤ 2π set by (2.5) is met. Hence, an alternative form of this delta function, suitable for first performing the integration with respect tor s in (3.1), is and v l3 = [r s cos α cos θ − (−1) l sin α sin θ(1 − cot 2 α cot 2 θ) 1/2 ]ê z P −r s cos αn ∞  4)). The step function H ∞ in (3.5) ensures that the contribution from the mth rotation cycle reaches a far-field observer at (R P , θ P , ϕ P ) during the intervalR P + 2(m − 1)π ≤ ωt P ≤R P + 2mπ of observation time (see (2.34)). If the observation time is set at midpoint of this interval, i.e., ωt P =R P + (2m − 1)π, (3.9) the argument of the delta function in (3.5) assumes the form  3)). Note that there is no loss of generality in fixing the observation time: because t P only appears in the combination ϕ P − ωt P , the temporal variation of the radiation is equally well described by its dependence on ϕ P .
Replacing the delta function in (3.5) by its Fourier representation we therefore obtain in which we have interchanged the orders of integration with respect tor s and k. Note that the range of integration overr s is set by ∆ ≥ 0. A pair of integration variables more suitable than (r s , θ) for evaluating (3.12) are η and τ defined by (see (2.18)) and τ = arccos(csc α cos θ) (3.14) since the Jacobian of this transformation removes the singularities of the integrand that occur on the boundaries (θ = π/2 ± α and ∆ = 0 where c 1 = 0) of the integration domain. If we now expresŝ r s and θ everywhere in terms of η and τ by inverting (3.13) and (3.14) we obtain θ = arccos(sin α cos τ ), (3.16) and so can rewrite (3.12) as wherer s and θ henceforth stand for the functions of (η, τ ) given by (3.16) and (3.17). We shall see below that the limiting value of the ratio η/c 1 at η = 0, where c 1 vanishes, is finite. The task of the rest of this section is to evaluate the right-hand side of (3.18) by treating it as a repeated integral.

Evaluation of the asymptotic values of the integrals over η by the method of stationary phase
It follows from (3.10) in conjunction with (3.16) and (3.17) that ∂f ± l /∂η vanishes at η = 0, so that the main contribution towards the asymptotic value for large |k| of the integral over η in (3.18) comes from the vicinity of the cusp locus C of the bifurcation surface at which η is zero (see, e.g., [8]). We can therefore approximate the functions f ± l in the phase of the integrand of (3.18) by the following leading terms in their Taylor expansions in powers of η: f denotes the value ofr s at the cusp locus C (see (3.17)). Note that the third-order term has to be included in the above expansion to take account of the difference between f + l and f − l . Note, moreover, that according to (2.30a) and (3.10), c 1 can be written as [ 3 4 (f + l − f − l )] 1/3 which implies, in conjunction with (3.19), that c 1 2 −1/3 ξ near η = 0.
If we now apply the method of stationary phase to the η-integral, i.e., insert (3.19) in (3.18) and replace the amplitude of the exponential in the integrand of (3.18) by its value at the stationary point η = 0 (see, e.g., [8]), we obtain the following expression for the leading term in the asymptotic expansion of this integral for large |k|:

26)
and use has been made of the fact that η/c 1 = 2 1/3 (r 2 Pr 2 sC sin 2 θ − 1) 1/2 in the limit η → 0. (Note that b tends to 1 asR P tends to infinity.) The above asymptotic approximation for large |k| is justified since the frequency of the radiation we are interested in has a much higher value than the rotation frequency ω/2π of the central neutron star.
Hence, the unconventional contribution E b + − E b − to the radiation field E that was encountered in (2.40) can be written as once the integration with respect to ξ has been carried out analytically to obtain in terms of the Airy function Ai and in terms of the generalized hypergeometric function 3 F 4 (see, e.g., [9]). Mathematica has been used to perform the above integrations by first employing a change of integration variable to cast I Q and I P in the forms of sine and cosine Fourier transforms and using the relations in §5.5 of [9] to simplify the gamma functions in the resulting expressions.

Dominance of the contribution from large values of |k|
In this section we assess the expectation that the main contribution toward the value of the radiation field should come from large values of |k|. We will do this by evaluating the k-integral in (3.28) once with the exact value of its integrand and another time with the asymptotic value of its integrand for large |k| and comparing the outcomes of these two evaluations.
The functions P l , Q l and b in the integrand of (3.28) are independent of k. Once the Airy function that appears in the expression for I Q is expressed in terms of Bessel functions, the k-integral multiplying Q l assumes the form of a tabulated Fourier transform (cf., [10]) and so can be evaluated exactly to obtain The k-integral multiplying P l in (3.28) also consists of the sum of five Fourier transforms each of which can be evaluated explicitly by means of Mathematica. The result is  Note that repeated use has to be made of the relations in §5.5 of [9] to cast the result in the above simplified form.
The factor a that multiplies k in the arguments of the Airy function in (3.29) and the generalized hypergeometric functions in (3.30) is non-zero and positive everywhere and approaches the finite value a csc θ P csc θ(1 − cos θ P cos θ) (3.35) asR P tends to infinity (see (3.22)). To calculate the contribution of large |k| toward the values of the integrals in (3.31) and (3.33) we can therefore replace the functions I P and I Q that appear in the integrands of these integrals by the following leading terms in their asymptotic expansions for large |k| before performing the integrations. (Note that (3.37) is found more easily by a direct asymptotic evaluation of the integral over ξ in the first member of (3.30) than by the calculation of the limiting value of the second member of this equation.) To calculate the corresponding contribution from small values of |k| we should replace I P and I Q by the following leading terms in their expansions about k = 0 (see [9]). It turns out that the functions F e 1 and F e 2 , defined in (3.31)-(3.34), can be accurately approximated by and where C and S are the Fresnel cosine and sine integrals, respectively. The exact (F e 1 , F e 2 ) and approximate (F a 1 , F a 2 ) versions of the functions F 1 and F 2 are compared in figure 3.1 for k 1 = 3 5 a −3 and k 2 = 12 5 a −3 . The goodness of the fit in this figure and the shortness of the intervals 0 ≤ k ≤ k 1,2 relative to k 1,2 ≤ k < ∞ show that (i) the difference in values of the exact and approximate versions of these two functions is negligibly small in the intervals over which F 1 and F 2 make their main contributions toward the value of the integrand of the τ -integral in (3.28) and (ii) the contributions of large |k| toward the values of F e 1 and F e 2 is by far greater than those of small |k|. In order that we can derive an analytic expression for the radiation field E uc , we will accordingly base the analysis in the rest of this section on the following approximate values of I P and Replacing I P and I Q in (3.28) with their approximate values I a P and I a Q and interchanging the orders of integration with respect to k and τ , we write (3.28) as and proceed to evaluate the integral over τ first. That the unconventional radiation field E uc is given by the real part of the above expression (in which the integrals over −∞ < k < ∞ are written as twice the integrals over 0 ≤ k < ∞) is understood.  The relative positions of the critical angles θ P 1S and θ P 2S along the θ P -axis for α < π/3 and α > π/3. The ranges of values of θ P for which f 1C and f 2C have two turning points as functions of τ are also shown. Outside the shown intervals, both f 1C and f 2C vary monotonically with τ .

Critical points of the phase functions f lC andf lC
The first two derivatives with respect to τ of the functions f lC andf lC -defined by (3.21) and (3.43)-which appear in (3.42) are given in appendix A. It follows from (A.1)-(A.5) that the nature of the critical points of f lC is determined by the value of the coordinate θ P of the observation point: this function can have two turning points (a maximum and a minimum), can have a single inflection point, or can be monotonic. As indicated by (3.21), the changes θ → π − θ, θ P → π − θ P and ϕ P → ϕ P + π transform f 2C (τ ) into f 1C (τ ). This means that f 2C has the same kind of critical points as f 1C (τ ) but in a different hemisphere (in θ > π/2 and θ P > π/2 instead of θ < π/2 and θ P < π/2 and vice versa). Moreover, the function f lC for α > π/2 follows from that for α < π/2 if we replace θ by π − θ at the same time as replacing α by π − α (see (3.21)). It is sufficient, therefore, to consider only the cases in which α < π/2. Note that, owing to the presence of the factor w 1 in the expression for the density of the current sheet, the field E uc is zero for α = π/2 (see (1.7)).
In this paper, we denote the values of τ , θ P and ϕ P at which by τ lS , θ P lS , and ϕ P lS , respectively. It turns out that θ P 1S and θ P 2S always lie on opposite  for α = π/4 in frames (i) and (ii) and for α = 3π/8 in frames (iii) and (iv). In frames (i) and (ii), θ P > θ P 1S > θ P 2S for the blue curves a, θ P = θ P 1S > θ P 2S for the red curves b and θ P 2S < θ P < θ P 1S for the green curves c. In frames (iii) and (iv), θ P 1S < θ P < θ P 2S for the blue curves a, θ P = θ P 1S < θ P 2S for the red curves b and θ P < θ P 1S < θ P 2S for the green curves c. Note that a change in the coordinate ϕ P of the observation point shifts the above curves (which have here been plotted forR P = 10 6 ) vertically without changing their shapes. These curves illustrate that while in the case of π/3 < α < π/2 both f 1C and f 2C have maxima, minima or inflection points, in the case of 0 < α < π/3 either f 1C or f 2C is a monotonic function of τ .
Changes in values of the free parameter ϕ P simply shift the curve representing f lC versus τ up or down without altering its shape (see (3.21)). But changes in values of the remaining free parameterR P does alter the relative positions of the two turning points of f lC when θ P is close to θ P lS . The length of the interval separating the τ coordinates of the maximum and minimum of f lC decreases with increasingR P in a case where θ P is close to θ P lS and so this interval is small. If we denote the τ coordinates of the maximum and minimum of f lC by τ lmax and τ lmin , respectively, then it turns out that |τ lmax − τ lmin | ∝R −1/2 P forR P 1 when θ P has the value θ P lS |R P →∞ for any given α, i.e., when θ P is such that |τ lmax − τ lmin | shrinks to zero asR P tends to infinity. We shall see in § 4.4 that this property of f lC results in a decay of the radiation's intensity with distance in the direction of θ P lS that is slower than predicted by the inverse-square law.
Note that f 1C and f 2C for a given value of n are equal at τ = 0 and differ by 2π at τ = π (see (3.21) and (3.14)). In cases where the absolute value of f lC | τ =π − f lC | τ =0 for either l = 1 or l = 2 is greater than 2π (e.g., when α = 15 • and θ P = 65 • ) or f lC | τ lmin and f lC | τ lmax differ by more than 2π (e.g., when α = 85 • and θ P = 5 • ), ordinates of the points on f 1C and f 2C in figure 3.3 also span intervals whose lengths exceed 2π. In such cases, therefore, f 1C and f 2C for several values of n (i.e., several cycles of retarded time) simultaneously contribute toward the intensity of the pulse that is observed during a single period.
The above discussion applies also to the modified functionsf 1C andf 2C : their critical points differ from those of f 1C and f 2C only in their positions, not in their nature. The generic forms assumed byf 1C andf 2C are the same as those illustrated in figure 3.3 except that the role of α = π/3 in figure 3.2 is played by α = 0.8707129958.
Depending on relative positions of the coordinate θ P of the observation point and the inclination angle α, therefore, the four phase functions f 1C , f 2C ,f 1C andf 2C in the integrand of the τ -integral in (3.45) can jointly have a set of isolated stationary points with one to eight members (where ∂f lC /∂τ and/or ∂f lC /∂τ for l = 1 and/or l = 2 vanish) or can have one or two degenerate stationary points (where ∂ 2 f lC /∂τ 2 and/or ∂ 2f lC /∂τ 2 simultaneously vanish with ∂f lC /∂τ and/or ∂f 2C /∂τ ). The number of contributing stationary points is higher in cases where the ordinates of the curves depicted in figure 3.3 span intervals whose lengths exceed 2π and so the contributions from several cycles of retarded time are received during a single period of observation time.

A uniform asymptotic approximation to the integral over the colatitude of the source elements for large |k|
Since this paper is concerned with determining the radiation field E uc also at observation points for which the turning points of the phase functions f lC andf lC are close to one another or coalescent, we need to obtain an asymptotic approximation to the value of the τ -integral in (3.45) that is uniform with respect to the interval separating the nearby saddle points of these phase functions (see [5,8]).
In cases where f lC andf lC each have two turning points, one at their maxima τ lmax and τ lmax and one at their minima τ lmin andτ lmin , each of these functions can be tansformed into a cubic function via The transformation of the integration variable in (3.45) from τ to λ thus results in whose asymptotic value for large k can be written, as shown by [5], as (see also [8]). Evaluating the integrals over λ (cf. [9]), we obtain where Ai and Ai are the Airy function and the derivative of the Airy function with respect to its argument, respectively. Note that (P l , Q l ) are different from the vectors (P l , Q l ) defined in (3.26)-(3.27).
The indeterminate quantities (dτ /dµ) f lC | τ =τ lmin and (dτ /dµ) f lC | τ =τ lmax that appear in (3.61) have to be found by repeated differentiation of (3.48) with respect to λ and the evaluation of the resulting relations ∂f lC ∂τ and (3.63) at τ = τ lmin and τ = τ lmax . This procedure results in in which λ = ±σ l1 are the images of τ = τ lmin and τ = τ lmax , respectively. Likewise, a result that can be obtained by applying the procedure described above to the functionf lC . In this expression, too, λ = ±σ l1 are the images of τ =τ lmin and τ =τ lmax , respectively. Note that every one of the terms appearing in (3.61) would contribute toward the value of the radiation field only when the phase functions f 1C , f 2C ,f 1C andf 2C each have two turning points (see figure 3.2). If any one of these functions varies monotonically with τ , for the prescribed values of α and θ P , then the terms entailing the (non-existent) τ coordinates of its maximum and minimum should be omitted from (3.61).

The remaining integration with respect to k
To perform the k-integration in (3.61) we first need to render the k-dependence of the coefficients P l , Q l ,P l andQ l that appear in this equation explicit. This can be done by rewriting (3.55)-(3.60) as The quantities appearing in (3.68)-(3.74) are all independent of k, so that the remaining integrals in (3.61) are-according to (3.66)-of two types: those that extend over a finite interval 0 ≤ k ≤ k i and those that extend over a semi-infinite interval k i ≤ k < ∞, where i is either 1 or 2.
The rapid oscillations of the Airy functions for large k (cf. [9]) ensure that the integrals over 0 ≤ k ≤ k i in (3.61) receive their main contributions from the vicinity of k = 0. The ranges of these integrals can therefore be extended to 0 < k < ∞ without introducing an appreciable error. Moreover, each of the integrals over k i < k < ∞ is accurately approximated (according to numerical integrations) if it is written as the difference between two integrals with the same integrands but with the ranges 0 ≤ k < ∞ and 0 ≤ k ≤ k i and the Airy function in the integrand of the integral over 0 ≤ k ≤ k i is replaced by its value at k = 0, as in (3.80) below.
Once these approximations are applied to (3.61) and the Airy functions that appear in the resulting equation are expressed in terms of Bessel functions, the k-integrals in question assume the form of tabulated Fourier transforms (cf., [9,10]) and can be evaluated analytically to arrive at In the above expressions, 2 F 1 is the hypergeometric function, Γ(ν, x) is the incomplete gamma function and the integration variable κ is related to k via κ = 2 3 σ 3 l1 k. Note that G 4 is here found by performing the integrations in the first line of (3.81) for ν > 0 and evaluating the analytic continuation of the resulting expression at ν = −1/6. This yields a value for G 4 that exactly agrees with the outcome of the numerical evaluation of the integrals defining this function.
According to (3.82), the variable η l that appears in the arguments of the functions G 1 , · · · , G 4 equals 1 when f lC | τ =τ lmin = 0 and equals −1 when f lC | τ =τ lmax = 0. This holds true, as indicated by (3.83), also forη l whenf lC vanishes at its maximum or minimum. Moreover, η l assumes an infinitely large value at the point where the maximum and minimum of f lC coalesce and so σ l1 vanishes, an unbounded upper limit that is also approached byη l as the turning points off lC coalesce to form an inflection point.

The divergence that arises from the vanishing width of the current sheet and its regularization
The functions G 1 , · · · , G 4 that appear in (3.76) and (3.77) approach the following divergent values as their argument x tends to ±1 from inside and outside the interval −1 < x < 1 (see figure 3.4): ). The radiation field E uc correspondingly diverges when η l = ±1 (or η l = ±1), i.e., when f lC (orf lC ) vanishes at one of its turning points. The above singularities in the expression for the radiation field stem from assigning a zero width to the current sheet. Because its charge and current densities are proportional to a Dirac delta function, the current sheet described by (1 .18) and (1.19) has a vanishing thickness. The vanishing thickness of the current sheet in turn results in an infinitely wide range of values for the variable k that appears in its Fourier representation (see (3.11)). But, given that it is created by the coordinated motion of aggregates of subluminally moving particles, a superluminally moving source is necessarily volume-distributed: it can neither be point-like nor be distributed over a line or a surface [11]. In a physically more realistic model of the magnetosphere, where the processes that occur on plasma scales within the current sheet are taken into account, this sheet would have a non-zero thickness and the singularities in question would not occur.
To circumvent the divergence that arises from overlooking the finite width of the current sheet, we will here replace the integration domain 0 ≤ k < ∞ in (3.61) with 0 ≤ k ≤ k u and treat the upper limit k u ( 1) on the range of values of |k| as a free parameter. The thickness of the current sheet is dictated by microphysical processes that are not well understood: the standard Harris solution of the Vlasov-Maxwell equations that is commonly used in analysing a current sheet [12] is not applicable in the present case because the current sheet in the magnetosphere of a neutron star moves faster than light and so has no rest frame. Introducing the upper limit k u is tantamount to assuming that the (unknown) thickness of the current sheet is of the order of 1/k u in units of the light-cylinder radius c/ω.
The singularities that arise from k = ∞ can thus be regularized by (i) changing the ranges of those integrals in (3.78)-(3.82) that extend over 0 ≤ k < ∞ to 0 ≤ k ≤ k u , (ii) (iii) equating the integrals over (0, ∞) on the right-hand side of (3.88) to the expressions found in (3.78)-(3.82) and (iv) replacing the Airy functions in the integrands of the integrals over k u ≤ k < ∞ by the leading terms in their asymptotic expansions for large k (cf., [9]). The integrals over k u ≤ k < ∞ can then be performed analytically to arrive at the following regularized versions, G r 1 , · · · , G r 4 , of G 1 , · · · , G 4 : in which the subtracted contributions arising from k u ≤ k < ∞ are given by Here, κ u = 2 3 σ 3 l1 k u and Ci and Si are the cosine and sine integrals, respectively. It can be easily verified that the singularities of the functions G u i at x = ±1 are the same as those of the functions G i that were derived in (3.78)-(3.82) and so cancel out in G r i when G u i are subtracted from G i (see figure 3.5).
The following singularity-free expressions for the two parts E uc P and E uc P of the radiation field E uc are thus obtained by replacing the functions G i in (3.76) and (3.77) by their regularized versions G r i : Characteristics of the resulting radiation

Pulse profiles and polarization position angles
As a result of receiving contributions from multiple stationary points of the phase functions f lC andf lC , the unconventional component E uc of the radiation field E has an amplitude that exceeds that of the conventional component E v of this field by many orders of magnitude (see (2.40)-(2.42) and (3.28) and the figures below). We can therefore calculate the Stokes parameters of the present radiation from the expressions based on the spherical components E uc ϕ P and E uc θ P of E uc alone. (The superscript * in the above expressions denotes complex conjugation.) It can be seen from (3.75) and (3.94)-(3.95) that the unconventional radiation field E uc consists of the sum of two distinct parts: one part, E uc P , depending on the vectors P l . As suggested by the occurrence of the factor i in (3.95), these two parts turn out to be out of phase with one another by approximately π/2. By calculating the Stokes parameters for E uc P and E uc Q separately, we will show here that the polarization position angles associated with these two fields are approximately orthogonal to one another in general. Accordingly, the two distinct parts of the radiation field defined by E uc P and E uc Q are respectively referred to in this paper as the P and Q polarization modes.
In this section we begin with evaluating the Stokes parameters (I, V, L) of the radiation analysed in §4 in units of (see (1.22) and (3.94)-(3.95)) as functions of the longitude ϕ P atR P 1 for various values of the colatitude θ P , the inclination angle α and the lower bound k −1 u on the width (in units of the light-cylinder radius) of the current sheet. We will evaluate the exact expression for the radiated field at a suitably large value ofR P , rather than proceeding to the far-field limit R P → ∞, because as we have already pointed out in § 3.4 the relative positions of the stationary points of the phase functions f lC andf lC depend onR P at colatitudes for which the maxima and minima of these functions are close to one another. However, away from such colatitudes the shapes of the pulse profiles and position-angle distributions we will be discussing do not change perceptibly with distance once the value ofR P exceeds 10 6 . We will therefore useR P = 10 6 for plotting most of the figures in this section.
The number of components of a pulse profile is determined by the total number of stationary points of the four phase functions f 1C , f 2C ,f 1C andf 2C described in § 3.4 and the number of values of n (i.e., the number of cycles of retarded time) that contribute to the radiation received during a single period of observation time. The longitudinal interval occupied by various components of a pulse profile is determined by the separation between the τ coordinates of the stationary points of these phase functions. To be able to depict the pulse profiles over the entire longitudinal intervals occupied by their various components while displaying the finite width of each component, we will plot the profiles of high-intensity pulses for k u = 10 2 . It should be borne in mind, however, that intensity at the peaks of the main pulses in these figures is a linearly increasing function of k u (see the last paragraph of § 3.7). The higher the value of k u is, the narrower are the rapid low-amplitude modulations (microstructure) of the pulse profile. For k u = 10 4 and higher, these modulations are too sharp and dense to show up in most of the figures plotted here. We will also adopt those branches of the multi-valued function arctan appearing in (4.3) that yield continuous polarization position-angle distributions across various components of a given pulse.
In most of the examples given below, the inclination angle of the magnetic axis and the colatitude of the observation point are set in the upper hemisphere 0 < θ P < 90 • . At any given value of the inclination angle α, the pulse observed at 180 • − θ P differs from that observed at θ P only in that the intensity V of its circularly polarized part is replaced by −V and the longitude ϕ P is replaced by ϕ P + 180 • . Moreover, the results for α > 90 • follow from those for α < 90 • by replacing θ P , ϕ P and V by 180 • − θ P , ϕ P + 180 • and −V , respectively (see § 3.4). Note also that in these examples the range of values of the azimuthal coordinate ϕ P that is spanned by the observation point across the pulse window differs from that given by (3.94)-(3.95) for any n (see also (3.21)): the origin of this coordinate is shifted in each case to place the starting point of the plotted pulse profile at ϕ P = 0.
Some examples of the linearly and circularly polarized intensity and polarization position angle distributions of the pulses described by (3.75), (3.94) and (3.95) are shown in figures 4.1-4.8.
As suggested by figure 4.9, a radically different type of pulse is detected when the colatitude θ P of the observation point has a value close to (or equal to) θ P lS (R P , α) or θ P lS (R P , α) for which the extrema of one (or more) of the phase functions f lC andf lC coalesce into an inflection point (see § 3.4). In the example plotted in figure 4.9, the critical angle θ P 2S happens to lie within a distance of the order of 1/R P radians from the longitude π/2 of the observation point. When sampled over a wide longitudinal interval, the profile that is shown in part a of   In this case, the stationary points off C2 alone contribute toward the field. Although not easily discernible because of the low value of its intensity V , circular polarization reverses sense across the right-hand component of this pulse. In contrast to those shown in the previous figures, on the other hand, the position angles of the two polarization modes are essentially coincident in this case. This is the only pulse occurring in the entire pulse window.       of the azimuthal coordinate of this peak graphically: by plotting the distribution of the Stokes parameter I over successively shorter longitudinal intervals centred on the peak of the distribution until the maximum value of I stops growing. This zooming in procedure reveals not only the pulse shown in figure 4.9(b) but also another as narrow and as intense pulse at a longitudinal distance of 4 × 10 −18 deg from it. In addition, the coincidence of the limiting values of θ P 1S andθ P 2S for R P → ∞ results, in this case, in two identical pulses whose longitudes are separated by 180 • . The features exhibited by the pulse in figure 4.9 can be inferred from (3.75), (3.94) and (3.95) also analytically. The only variables in the expression for E uc that depend on the observer's longitude are σ l2 andσ l2 which vary linearly with ϕ P (see (3.51) and (3.52) and note that, according to (3.21), ϕ P drops out of the expressions for σ l1 andσ l1 ). The variables σ l2 andσ l2 , on the other hand, appear in (3.94)-(3.95) only in the combinations σ l2 /σ 3 l1 andσ l2 /σ 3 l1 . Hence, in cases where the turning points of the phase functions are sufficiently close to one another for σ l1 orσ 11 to be appreciably smaller than 1, as in figure 4.9, the arguments of the functions G r i that appear in the expression for E uc are highly sensitive functions of ϕ P . Not only the widths but, as indicated by (3.78)-(3.81) and (3.89)-(3.93), also the amplitudes of G r i vary sharply with ϕ P when σ l1 orσ 11 assume values that are close to zero.

Brightness temperature
The brightness temperature T b of the present radiation can be calculated by equating the magnitude of the Poynting flux of this radiation (c|E uc | 2 /4π) to the Rayleigh-Jeans law (2k B T b ν 2 ∆ν/c 2 ) for the energy that a black body of the same temperature would emit per unit time per unit area into the frequency band ∆ν centred on the frequency ν, where k B is the Boltzmann constant. The resulting equation can then be solved for T b to obtain This in conjunction with (4.1) shows that T b is related to the dimensionless Stokes param-eterÎ = I/I 0 by where r s0 and B 0 are the radius of the star and the magnitude of the star's dipolar field at its magnetic pole, respectively (see the first paragraph of § 1).
The limiting values of θ P 1S and θ P 2S in the second column are 48.533945294618400228 • and 33.932818533330613261 • , respectively.
we obtain in which the value ofÎ is specified, as in the case of figures 4.1-4.9, by the numerical evaluation of the Stokes parameter I in units of I 0 at the highest peak of the pulse detected atR P = 10 13 (i.e., atR P = 1.028 × 10 13 DP −1 when the factor 1.028DP −1 equals unity). The brightness temperature implied by (4.9) and k u = 10 7 is listed in Table 4.1 for the pulses depicted in figures 4.1, 4.4 and 4.9 and for a pair of examples of the pulses that are detected at the critical colatitudes lim R P →∞ θ P lS . Table 4.1 also shows the full width at half maximum δϕ P of the listed pulses (see figure 4.9(b)). Once resolved, the longitudinal distributions of the narrow pulses that stem from the focusing of the radiation when σ l1 and σ l1 are small all have the same shape as that of the pulse shown in figure 4.9(b). Note that, as indicated by the last column of Table 4.1, the pulse profiles depicted in figures 4.1-4.8 have to be plotted on considerably shorter longitudinal scales before their peaks assume the shape shown in figure 4.9(b) and the maximum values of their dimensionless intensityÎ can be discerned graphically (see § 4.1). In general, as one reduces the longitudinal interval over which I is plotted, the peak of the pulse splits in two before the finite widths of either of the partitioned pulses are visible.
Values of T b higher than those listed in Table 4.1 are predicted by (4.9) when the colatitude of the observation point lies closer to one of the critical angles θ P lS orθ P lS . In the case of α = 60 • , for example, the listed value (1.17 × 10 40 • K) of T b corresponds to θ P = lim R P →∞ θ P 1S = lim R P →∞ θ P 2S = π/2. For an observation point whose colatitude is closer to the critical angles in question than π/2 is, e.g., for θ P = θ P 1S + 10 −20 rad, T b and δϕ P have the values 1.17 × 10 54 • K and 2.26 × 10 −34 deg, respectively.

Frequency spectrum
Given that the radiation field E uc depends on the observation time t P only in the combination ϕ P − ωt P , the frequency spectrum of the present radiation is equally well described by the Fourier decomposition of the field E uc with respect to the azimuthal angle ϕ P of the observation point. In the present case, the content of this spectrum stems from two factors. One factor is the thickness of the current sheet ( c/(k u ω) which is manifested in the sharp small-amplitude modulations (microstructure) of the pulse profile (see figures 4.1 and 4.5): the wavelengths of such modulations are proportional to k −1 u . The other factor is the full width at half maximum (δϕ P ) of the pulse with the highest peak in the pulse profile (see figure 4.9(b)): the fraction δϕ P /2π of a period during which such narrow pulses propagate past a detector is by many orders of magnitude smaller than the fraction k −1 u of the light-cylinder-radius in cases where the colatitude of the observation point is close to or coincident with one of the critical angles described in § 3.4 (see Table 4.1). While the Fourier decomposition of the fluctuations associated with the first factor yields a frequency spectrum centred on radio waves when k u 10 5 , that of the fluctuations associated with the second factor yields a wide spectral distribution extending to gamma rays: the value δϕ P = 1.21 × 10 −26 radian appearing in the last column of Table 4.1 corresponds to a frequency spectrum that extends as far as ν ω/(2πδϕ P ) 1.31 × 10 27P −1 Hz.
Our replacing the Dirac delta function in (3.5) by its Fourier representation (3.11) is tantamount to Fourier analysing the fluctuations of the radiation field that arise from the short thickness of the current sheet with respect to ϕ P since the argument of that delta function depends on ϕ P linearly. The spectral distribution of the part of the radiation that stems from the thickness of the current sheet is therefore given by the k dependence of i.e., the k dependence of the square of the modulus of the integrand that appears in (3.61). The frequency ν of the radiation is related to k through ν = (2πk)/ω. At harmonic numbers k for which the arguments of the Airy functions in (4.10) are smaller or of the order of unity, these functions assume values that are independent of frequency. When their arguments are large, on the other hand, they reduce to Ai(−k 2/3 σ 2 l1 ) π −1/2 σ −1/2 l1 k −1/6 cos 2 3 kσ 3 l1 − 1 4 π , (4.11) and Ai (−k 2/3 σ 2 l1 ) π −1/2 σ 1/2 l1 k 1/6 sin 2 3 kσ 3 l1 − 1 4 π . |P l | k −1/2 |Q l | and |P l | k −1/3 |Q l | 2/3 1 |P l | k −1/2 |Q l | and |P l | k −1/3 |Q l | 4/3 1 |P l | k −1/2 |Q l | and |P l | k −1/3 |Q l | 5/3 2 |P l | k −1/2 |Q l | and |P l | k −1/3 |Q l | 7/3 2 Table 4.2: Values, in various regimes, of the spectral index β (defined in (4.13)) for the part of the radiation associated with the sharp small-amplitude modulations (microstructure) of the pulse profile.
Equation (4.10) and these limiting values of the Airy functions jointly yield the dependencẽ of the radiation intensityĨ on frequency and the values that the spectral index β can assume in various regimes. The variables σ l1 and σ l2 and the vector functions P l and Q l that appear in (4.10) are independent of k. When |P l | k −1/2 |Q l |, the vector K l and hence P l and Q l are also independent of k and are by a factor of the order of k 1/2 larger thanP l andQ l (see (3.55)-(3.60)). In this case, the possible values of the spectral index β are determined by the relative magnitudes of |P l | and |Q l | only. If |P l | k −1/3 |Q l |, then β = 2/3 when the Airy functions in the first square bracket in (4.10) are of the order of unity and β = 1 when these Airy functions have the limiting values given by (4.11) and (4.12) and so the first square bracket in (4.10) decays as k −1/6 . If |P l | k −1/3 |Q l |, then β = 4/3 when the arguments of the Airy functions in question are of the order of unity and β = 1 when the first square bracket in (4.10) decays as k −1/6 . When σ l1 is small and the second term of the first square bracket in (4.10) dominates, there is a short frequency interval in which the spectral intensityĨ increases with increasing k.
In the opposite regime |P l | k −1/2 |Q l |, the factor k −1/2 multiplying Q l in (3.59) and (3.60) reduces the value of the spectral index β by 1 everywhere (see Table 4.2).

Flux density and its rate of decay with distance
Flux density of a radiation is given, in general, by the magnitude of the Poynting vector, c|E uc | 2 /4π, which has the dimensions of erg/(cm 2 ×sec) in cgs units. In the present case, however, the linear extents in the azimuthal direction, R P δϕ P , of the focused radiation beams that embody the high-frequency (optical to gamma-ray) radiation are invariably smaller than 1 cm at R P = 1 kpc (see Table 4.1). The amount of energy that crosses a unit area per unit time is therefore given by S = c 4π |E uc | 2 δϕ P = 2.79 × 10 −3Î δϕ P B 0 d 2 P D 2 erg cm 2 × sec , (4.14) in which δϕ P is in radians (see (4.1), (4.4), (4.7) and (4.8)). Note that the linear extents in the latitudinal direction, R P δθ P , of the focused radiation beams that embody the highfrequency radiation are of the order of the light-cylinder radius, c/ω, in general: these beams remain fully in focus at all distancesR P over the latitudinal interval δθ P |θ P lS − lim R P →∞ θ P lS | which turns out to be of the order ofR −1 P independently of the values of the other parameters.
In the case of α = 60 • , θ P = 90 • , k u = 10 7 and D = 1 kpc depicted in figure 4.9, for example, the flux density S has the value 32.1 S 0 erg/(sec×cm 2 ), where (4.15) At latitudes closer to or further away from the critical angle for this example (θ P = 90 • ), the degree of focusing of the radiation beam and so the value of the flux density S is, respectively, higher or lower. As pointed out in § 3.4, the length of the interval |τ lmax − τ lmin | separating the τ coordinates of the maximum and minimum of the phase function f lC decreases asR −1/2 P with increasingR P in a case where this interval is small, i.e., when the colatitude θ P of the observation point has the critical value θ P lS (L, α), with L >R P . In particular, if the observation point has the colatitude lim R P →∞ θ P lS (or π − lim R P →∞ θ P lS ), then the maximum and minimum of f lC coalesce into an inflection point only atR P → ∞, rather than at a finite distance L. (These statements hold true also when f lC and θ P lS are replaced byf lC andθ P lS , respectively.) In the case illustrated in figure 4.9, for example, the colatitude of the observation point equals lim R P →∞ θ P 1S = 90 • so that at the finite distanceR P = 10 13 , the τ coordinates of maximum and minimum of f 1C are separated by the short interval 3.05 × 10 −5 degrees. It follows from the expression for ∂f 1C /∂τ in (A.1) that this separation has the value 3.05 × 10 −5 (R P /10 13 ) −1/2 for allR P .
The enhanced focusing of the radiation with distance that is caused by this shortening of the separation between the turning points of the phase functions results in a slower decay rate of the flux density with distance than that predicted by the inverse-sqaure law. Along colatitudes close to θ P lS orθ P lS , the flux density S of the radiation diminishes with increasing distance from its source asR −3/2 P instead ofR −2 P . This dependence of S on R P , or equivalently D, is illustrated in figure 4.10 in the case where α = 60 • , θ P = 90 • , k u = 10 7 , and D ranges from 0.1 to 10 5 kpc, i.e., from a galactic to a cosmological distance. The blue line with the slope −3/2 is the best fit to the red dots whose coordinates are determined by evaluating (4.14). The violation of the inverse-square law illustrated in this figure remains in force all the way to infinity whenever the colatitude of the observation point coincides with or is close to one of the eight angles given by lim R P →∞ θ P lS , lim R P →∞θP lS , π − lim R P →∞ θ P lS and π − lim R P →∞θP lS .
The violation of the inverse-square law encountered here (i.e., the fact that S is proportional to D −3/2 rather than being proportional to D −2 ) is not incompatible with the requirements of the conservation of energy because the radiation process discussed in this paper is intrinsically transitive. Temporal rate of change of the energy density of the radiation generated by this process has a time-averaged value that is negative (instead of being zero as in a conventional radiation) at points where the envelopes of of the wave fronts emanating from the constituent volume elements of the source distribution are cusped. The difference in the fluxes of power across any two spheres centred on the star is in this case balanced by the change with time of the energy contained inside the shell bounded by those spheres (see appendix C of [2] for a detailed discussion of this point).