Houshang Ardavan, Radiation by the superluminally moving current sheet in the magnetosphere of a neutron star, Monthly Notices of the Royal Astronomical Society, Volume 507, Issue 3, November 2021, Pages 4530–4563, https://doi.org/10.1093/mnras/stab2354

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 firmly 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. However, the role played by the superluminal motion of this current sheet in generating the multiwavelength, focused pulses of radiation that we receive from neutron stars is unknown. 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 salient features (brightness temperature, polarization, spectrum and profile with microstructure and with a phase lag between the radio and gamma-ray peaks) are strikingly similar to those of the emission received from pulsars. In addition, the flux density of this radiation diminishes with the distance D from the star as D−3/2 (rather than D−2) in certain latitudinal directions: a result that suggests that the high energetic requirements normally attributed to magnetars and the sources of fast radio bursts and gamma-ray bursts could be artefacts of the assumption that the radiation fields of all sources necessarily decay as predicted by the inverse-square law.

From the results obtained by numerical simulations of the magnetospheric structure of an obliquely rotating neutron star (Spitkovsky 2006; Kalapotharakos, Contopoulos & Kazanas 2012; Tchekhovskoy, Spitkovsky & Li 2013; Philippov & Spitkovsky 2018), Tchekhovskoy, Philippov & Spitkovsky ( 2016) have derived a semi-analytic description of the distributions of the electric and magnetic fields that permeate the plasma surrounding these objects (see also Bogovalov 1999). 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 (Section  2). 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 (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 (Bolotovskii & Ginzburg 1972; Ginzburg 1972; Bolotovskii & Bykov 1990).

(a) Snapshot of a single turn of the current sheet about the light cylinder. 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.

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 multiwavelength, 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 (Philippov & Spitkovsky 2018; Philippov et al. 2019).

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 (Section  2). 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 (Figs  1b– d). 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 (Section  3).

Constructive interference of the emitted waves and formation of caustics thus play a crucial role in determining the radiation field of the current sheet (Section  4). Its double-peaked pulse profiles and S-shaped polarization position angle distributions stem from the caustics associated with the nearby stationary points of certain phase functions (Section  5.1). What underpin the high brightness temperatures and the broad frequency spectra of this radiation are the extraordinary values of the amplitude and width of the pulses that are generated when the maxima and minima of the phase functions in question coalesce into inflection points (Sections  5.2 and  5.3). There is always a latitudinal direction along which the separation between any pair of these nearby maxima and minima decreases with increasing distance from the source (Section  4.4). The enhanced focusing of the emitted waves that takes place as a result of the diminishing separation between the nearly coincident stationary points of a phase function gives rise, in turn, to a lower rate of decay of the flux density of the radiation with distance (Section  5.3). That the decay of the present radiation with distance disobeys the inverse-square law is not incompatible with the requirements of the conservation of energy because the radiation process discussed here is intrinsically transient: the difference in the fluxes of power across any two spheres centred on the star is balanced by the change with time of the energy contained inside the shell bounded by those spheres (see Ardavan 2019, Appendix C). In addition to providing an explanation for the salient features of the multiwavelength emission received from pulsars (Lyne & Graham-Smith 2012; Abdo et al. 2013), the results reported in this paper suggest, therefore, that the high energetic requirements normally attributed to magnetars and the sources of fast radio bursts and gamma-ray bursts (Kumar & Zhang 2015; Kaspi & Beloborodov 2017; Petroff, Hessels & Lorimer 2019) could be artefacts of the assumption that the radiation fields of all sources necessarily decay as predicted by the inverse-square law (Section  6).

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 Ardavan 2019, Section  3).

The function |$G_j({\hat{r}},{\hat{\varphi }},{\hat{z}};{\hat{r}}_P,{\hat{\varphi }}_P,{\hat{z}}_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 jr = jz = 0 and |$j_\varphi =r^\prime \omega q\delta (r-r^\prime)\delta ({\hat{\varphi }})\delta (z)$| with a constant r′, then equation ( 43) 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/ω.

Generic forms of the function g(φ) for source points whose |$({\hat{r}},{\hat{z}})$| coordinates lie across the boundary Δ = 0 delineating the projection of the cusp curve of the bifurcation surface on to the |$({\hat{r}},{\hat{z}})$| plane. Depending on whether ϕ lies outside or inside the interval (ϕ−, ϕ+), contributions are made towards the observed field 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 for |${\hat{r}}=3$| and shows only two rotation cycles.

The source points inside the bifurcation surface, close to its cusp, make their contributions towards the observed value of the field at three distinct retarded positions in their trajectory (where a horizontal line g = ϕ in Fig.  2 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 Fig.  2). For the source points on the bifurcation surface (i.e. those for which g = ϕ± in Fig.  2), two of the contributing retarded positions coalesce at the extrema of the curve Δ > 0 in Fig.  2 giving rise to a divergent value of the Green’s function at P (see Ardavan 2019, figs 9 and 10). For the source points located on the cusp locus C of the bifurcation surface (i.e. those for which Δ = 0 in Fig.  2), all three of the contributing retarded positions coalesce at the inflection point of the curve Δ = 0 in Fig.  2 giving rise to a higher order singularity in Gj.

It follows from equations ( 55) and ( 56) that the factor |$\partial G_j/\partial {\hat{\varphi }}$| in the integrand of the integral in equation (42) 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 (31) under the integral sign when calculating the field. Had we evaluated the integral in equation ( 31) 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 Hadamard ( 2003), Hoskins ( 2009) and the illustrative example in Appendix A of Ardavan ( 2019).

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 |$[\Phi ({\textbf {x}}_P,t_P),{\bf A}({\textbf {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 |${\hat{\varphi }}$| coordinates |${\hat{\varphi }}_\pm$| of the two sheets of the bifurcation surface depend on the observation time tP (see equations 61a,b and  39), so that these two sheets move across the |${\hat{\varphi }}$| extent of the source distribution as tP 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 |${\hat{\varphi }}$| designating the initial azimuthal positions of the constituent volume elements of a source that fully occupies an annular region range over the interval |$0\le {\hat{\varphi }}\lt 2\pi$|⁠ . The |$({\hat{r}},{\hat{z}})$| coordinates of these source elements fall either in Δ ≥ 0 or in Δ < 0. The elements in Δ ≥ 0 are always divided into two sets: a set inside the bifurcation surface for which |${\hat{\varphi }}_-\le {\hat{\varphi }}\le {\hat{\varphi }}_+$| and so the Green’s function Gj has the form |$G_j^{\rm in}$| and a set outside for which |${\hat{\varphi }}$| lies either in |$(0,{\hat{\varphi }}_-)$| or in |$({\hat{\varphi }}_+,2\pi)$| and so Gj has the form |$G_j^{\rm out}$| (see equation 55). On the other hand, if the position of the observation point is such that Δ < 0 for all values of |$({\hat{r}},{\hat{z}})$| within the magnetosphere, then the source lies entirely outside the bifurcation surface and Gj has the form |$G_j^{\rm sub}$|⁠ . Note that, for certain space–time coordinates of the observation point P, the values of |${\hat{\varphi }}_-$| and |${\hat{\varphi }}_+$| that lie in the interval (0, 2π) could correspond to different rotation periods, i.e. to different values of m (see equations 49, 52, and  53). To simplify the notation, here we adopt an observation time tP at which the values of |${\hat{\varphi }}_-$| and |${\hat{\varphi }}_+$| that lie in the interval (0, 2π) correspond to the same rotation period m.

The task of the rest of this section is to evaluate the right-hand side of equation ( 90) by treating it as a repeated integral.

In this section, we assess the expectation that the main contribution towards the value of the radiation field should come from large values of |k|. We will do this by evaluating the k-integral in equation ( 100) 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 function |${\cal F}_1^{\rm e}$| diverges as (1 − x2)−1/2 at |x| = 1 −. Moreover, the generalized hypergeometric functions appearing in equation ( 106) are singular in the limit where their arguments approach unity. While the singularities of the individual terms of equation ( 106) at x = 2, x = −2 −, and x = 0 cancel out when added together, the function |${\cal F}_2^{\rm e}$| diverges as (x + 2)−1/2 at x = −2 + (see Fig.  3).

The solid (red) lines in this figure depict the exact values (⁠|${\cal F}_1^{\rm e}$| and |${\cal F}_2^{\rm e}$|⁠ ) of the functions derived in equations ( 104) and ( 106). The dashed (blue) lines depict the approximate versions (⁠|${\cal F}_1^{\rm a}$| and |${\cal F}_2^{\rm a}$|⁠ ) of these functions given by equations ( 113) and ( 114) for k1 = 3/(5a3) and k2 = 12/(5a3).

The first two derivatives with respect to τ of the functions flC and |${\bar{f}}_{lC}$| – defined by equations ( 93) and ( 118) – which appear as phases of the exponential factors in the integrand of equation ( 117) are given in Appendix  A.

It follows from equations ( A1)–( A5) that the nature of the critical points of flC 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 equation ( 93), the changes θ → π − θ, θP → π − θP, and φP → φP + π transform f2C(τ) into f1C(τ). This means that f2C has the same kind of critical points as f1C(τ) but in a different hemisphere (in θ > π/2 and θP > π/2 instead of θ < π/2 and θP < π/2 and vice versa). Moreover, the function flC for α > π/2 follows from that for α < π/2 if we replace θ by π − θ at the same time as replacing α by π − α (see equation 93). It is sufficient, therefore, to consider only the cases in which α < π/2. Note that, owing to the presence of the factor w1 in the expression for the density of the current sheet, the field Euc is zero for α = π/2 (see equation 7).

Relative positions of the critical angles θP1S and θP2S along the θP-axis for α < π/3 and α > π/3. The ranges of values of θP for which f1C and f2C have two turning points as functions of τ are also shown. Outside the shown intervals, both f1C and f2C vary monotonically with τ.

Plots of flC versus τ for various values of the free parameters θP and α display the following features. When α < π/3, (i) f1C is monotonic for 0 < θP < θP1S and has a maximum and a minimum for θP1S < θP < π, (ii) f2C has a maximum and a minimum for 0 < θP < θP2S and is monotonic for θP2S < θP < π, and (iii) the extrema of f1C and f2C coincide at an inflection point for θP = θP1S and θP = θPS2, respectively. Hence, both f1C and f2C are monotonic for θP2S < θP < θP1S when 0 < α < π/3. On the other hand, when π/3 < α < π/2, (i) f1C is monotonic in 0 < θP < θP1S and has two turning points in θP1S < θP < π, (ii) f2C has a maximum and a minimum in 0 < θP < θP2S and is monotonic in θP2S < θP < π, and (iii) once again the extrema of f1C and f2C coincide at an inflection point for θP = θP1S and θP = θP2S, respectively. Hence, f1C and f2C each have a maximum and a minimum for θP2S < θP < θP1S when π/3 < α < π/2. Fig.  5 illustrates some of the generic forms assumed by f1C and f2C as functions of τ inside the integration domain 0 < τ < π.

Dependence of the phase functions f1C and f2C on the integration variable τ for α = π/4 in frames (i) and (ii) and for α = 3π/8 in frames (iii) and (iv). In frames (i) and (ii), θP > θP1S > θP2S for the blue curves a, θP = θP1S > θP2S for the red curves b, and θP2S < θP < θP1S for the green curves c. In frames (iii) and (iv), θP1S < θP < θP2S for the blue curves a, θP = θP1S < θP2S for the red curves b, and θP < θP1S < θP2S 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 for |${\hat{R}}_P=10^6$|⁠ ) vertically without changing their shapes. These curves illustrate that while in the case of π/3 < α < π/2 both f1C and f2C have maxima, minima, or inflection points, in the case of 0 < α < π/3 either f1C or f2C is a monotonic function of τ.

Changes in values of the free parameter φP simply shift the curve representing flC versus τ up or down without altering its shape (see equation 93). But changes in values of the remaining free parameter |${\hat{R}}_P$| do alter the relative positions of the two turning points of flC when θP is close to θPlS. The length of the interval separating the τ coordinates of the maximum and minimum of flC decreases with increasing |${\hat{R}}_P$| in a case where θP is close to θPlS and so this interval is small. If we denote the τ coordinates of the maximum and minimum of flC by τlmax and τlmin, respectively, then it turns out that |$\vert \tau _{l{\rm max}}-\tau _{l{\rm min}}\vert \propto {\hat{R}}_P^{-1/2}$| for |${\hat{R}}_P\gg 1$| when θP has the value |$\theta _{PlS}\vert _{{\hat{R}}_P\rightarrow \infty }$| for any given α, i.e. when θP is such that |τlmax − τlmin| shrinks to zero as |${\hat{R}}_P$| tends to infinity. We shall see in Section  5.5 that this property of flC results in a decay of the radiation’s flux density with distance in the direction of θPlS that is slower than predicted by the inverse-square law.

Note that f1C and f2C for a given value of n are equal at τ = 0 and differ by 2π at τ = π (see equations 93 and  86). In cases where the absolute value of flC|τ = π − flC|τ = 0 for either l = 1 or l = 2 is greater than 2π (e.g. when α = 15° and θP = 65°) or |$f_{lC}\vert _{\tau _{l{\rm min}}}$| and |$f_{lC}\vert _{\tau _{l{\rm max}}}$| differ by more than 2π (e.g. when α = 85° and θP = 5°), ordinates of the points on f1C and f2C in Fig.  5 also span intervals whose lengths exceed 2π. In such cases, therefore, f1C and f2C for several values of n (i.e. several cycles of retarded time) simultaneously contribute towards the intensity of the pulse that is observed during a single rotation period.

The above discussion applies also to the modified functions |${\bar{f}}_{1C}$| and |${\bar{f}}_{2C}$|⁠ : their critical points differ from those of f1C and f2C only in their positions, not in their nature. The generic forms assumed by |${\bar{f}}_{1C}$| and |${\bar{f}}_{2C}$| are the same as those illustrated in Fig.  5 except that the role of α = π/3 in Fig.  4 is played by α = 0.8707129958 rad.

Depending on relative positions of the coordinate θP of the observation point and the inclination angle α, therefore, the four phase functions f1C, f2C, |${\bar{f}}_{1C}$|⁠ , and |${\bar{f}}_{2C}$| in the integrand of the τ-integral in equation ( 117) can jointly have a set of isolated stationary points with one to eight members (where ∂flC/∂τ and/or |$\partial {\bar{f}}_{lC}/\partial \tau$| for l = 1 and/or l = 2 vanish) or can have one or two degenerate stationary points (where ∂2flC/∂τ2 and/or |$\partial ^2{\bar{f}}_{lC}/\partial \tau ^2$| simultaneously vanish with ∂flC/∂τ and/or |$\partial {\bar{f}}_{2C}/\partial \tau$|⁠ ). The number of contributing stationary points is higher in cases where the ordinates of the curves depicted in Fig.  5 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.

Since this paper is concerned with determining the radiation field Euc also at observation points for which the turning points of the phase functions flC and |${\bar{f}}_{lC}$| are close to one another or coalescent, we need to obtain an asymptotic approximation to the value of the τ-integral in equation ( 117) that is uniform with respect to the interval separating the nearby stationary points of these phase functions (see Chester et al. 1957; Bleistein & Handelsman 1986).

Note that every one of the terms appearing in equation ( 133) would contribute towards the value of the radiation field only when the phase functions f1C, f2C, |${\bar{f}}_{1C}$|⁠ , and |${\bar{f}}_{2C}$| each have two turning points (see Fig.  4). 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 equation ( 133).

The rapid oscillations of the Airy functions for large k (cf. Olver et al. 2010) ensure that the integrals over 0 ≤ k ≤ ki in equation ( 133) 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 ki < 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 ≤ ki and the Airy function in the integrand of the integral over 0 ≤ k ≤ ki is replaced by its value at k = 0, as in equation ( 152) below.

According to equation ( 154), the variable ηl that appears in the arguments of the functions |${\cal G}_1,\cdots ,{\cal G}_4$| equals 1 when |$f_{lC}\vert _{\tau =\tau _{l{\rm min}}}=0$| and equals −1 when |$f_{lC}\vert _{\tau =\tau _{l{\rm max}}}=0$|⁠ . This holds true, as indicated by equation ( 155), also for |${\bar{\eta }}_l$| when |${\bar{f}}_{lC}$| vanishes at its maximum or minimum. Moreover, ηl assumes an infinitely large value at the point where the maximum and minimum of flC coalesce and so σl1 vanishes, an unbounded upper limit that is also approached by |${\bar{\eta }}_l$| as the turning points of |${\bar{f}}_{lC}$| coalesce to form an inflection point.

Real part (in red) and imaginary part (in blue) of the functions |${\cal G}_1,\cdots ,{\cal G}_4$| given by equations ( 150)–( 153) for κ1 = 10−2 and κ2 = 4 × 10−2. The limiting values of these functions across their discontinuities at x = ±1 are given by equations ( 156)–( 159).

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 equations ( 18) and ( 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 equation 83). 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 (Bolotovskii & Bykov 1990). 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 remove the singularity that arises from overlooking the finite width of the current sheet, we will here replace the integration domain 0 ≤ κ < ∞ in equations ( 150)–( 153) with 0 ≤ κ ≤ κu and treat the upper limit κu (≫ 1) on the range of values of κ as a free parameter. This is tantamount to setting a lower limit on the scale of fluctuations of the pulse profile that are manifested as its microstructure: it can be seen from equations ( 161)–( 165) and Fig.  7 that |$\kappa _{\rm u}^{-1}$| is proportional to the wavelength of modulations of the functions |${\cal G}^{\rm u}_i$|⁠ , and hence, of the microstructure of the pulse profile (see e.g. Fig.  8). The free parameter |$\kappa _{\rm u}^{-1}$| cannot of course be smaller than the (unknown) thickness of the current sheet in units of the light-cylinder radius.

Real part (in red) and imaginary part (in blue) of the regularized versions (⁠|${\cal G}^{\rm r}_1,\cdots ,{\cal G}^{\rm r}_4$|⁠ ) of the functions |${\cal G}_1,\cdots ,{\cal G}_4$| for κ1 = 10−2, κ2 = 4 × 10−2, and κu = 102. The higher the value of κu, the larger are the absolute values of the maxima and minima of these functions and the shorter are the wavelengths of their modulations (microstructure). In the limit κu → ∞, these functions approach those depicted in Fig.  6.

The Stokes parameters I, V, and the position angles ψ of the polarization modes P and Q at an observation point with the colatitude θP = 5° for the inclination angle α = 5° and κu = 102. Distribution of the Stokes parameter L is essentially coincident with that of I. In this case, only the stationary points of f2C and |${\bar{f}}_{2C}$| contribute towards the field. The large value of the intensity and the short duration of this pulse stem from the proximity of the maximum and minimum of f2C. At its peak, the right-hand component of this pulse has the intensity I = 1.05 × 1012I0 and the longitudinal width 6.76 × 10−9 deg when κu is 107. Note the reversal of the sense of circular polarization and the swing through 180° of the position angle of the Q polarization mode across the left-hand component of the pulse profile and the orthogonality of the position angles of the two polarization modes over the right-hand component of the pulse profile. This is the only pulse occurring in the entire pulse window.

It can be seen from equations ( 147) and ( 166)–( 167) that the unconventional radiation field Euc consists of the sum of two distinct parts: one part, |${\bf E}^{\rm uc}_{\cal P}$|⁠ , depending on the vectors |$\boldsymbol {{\cal P}}^{(j)}_l$| and |${{\boldsymbol {\bar{\cal P}}}}^{(j)}_l$| and another part, |${\bf E}^{\rm uc}_{\cal Q}$|⁠ , depending on the vectors |$\boldsymbol {{\cal Q}}^{(j)}_l$| and |${{\boldsymbol {\bar{\cal Q}}}}^{(j)}_l$|⁠ . As suggested by the occurrence of the factor i in equation ( 167), these two parts turn out to be out of phase with one another by approximately π/2. By calculating the Stokes parameters for |${\bf E}^{\rm uc}_{\cal P}$| and |${\bf E}^{\rm uc}_{\cal Q}$| separately, we will show that, in general, the polarization position angles associated with these two fields are approximately orthogonal to one another. Accordingly, the two distinct parts of the radiation field defined by |${\bf E}^{\rm uc}_{\cal P}$| and |${\bf E}^{\rm uc}_{\cal Q}$| are, respectively, referred to in this paper as the P and Q polarization modes.

The number of components of a pulse profile is determined by the number of stationary points of the four phase functions f1C, f2C, |${\bar{f}}_{1C}$|⁠ , and |${\bar{f}}_{2C}$| described in Section  4.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 positions of the peaks of various components of a pulse profile are determined by the φP coordinates of the points at which these phase functions and their derivatives with respect to τ vanish simultaneously. 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 κu = 102. It should be borne in mind, however, that the value of intensity at the peak of each of the shown pulse components is a linearly increasing function of κu (see the last paragraph of Section  4.7). The higher the value of κu, the narrower and closer to one another are also the rapid low-amplitude modulations of the pulse profile. For values of κu higher than 104, these modulations are too sharp and dense to show up in most of the figures plotted here. Furthermore, we will adopt those branches of the multivalued function arctan appearing in equation ( 170) 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 Section  4.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 equations ( 166)–( 167) for any n (see also equation 93): 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 full intensity, the linearly and circularly polarized intensity, and polarization position–angle distributions of the pulses described by equations ( 147), ( 166), and ( 167) are shown in Figs  8– 15.

As suggested by Figs  16 and  17(a), a radically different type of pulse is detected when the colatitude θP of the observation point has a value close to (or equal to) |$\theta _{PlS}({\hat{R}}_P,\alpha)$| or |${\bar{\theta }}_{PlS}({\hat{R}}_P,\alpha)$| for which the extrema of one (or more) of the phase functions flC and |${\bar{f}}_{lC}$| coalesce into an inflection point (see Section  4.4). In the example plotted in Figs  16 and  17, the critical angle θP2S happens to lie within a distance of the order of |$1/{\hat{R}}_P$| radians from the longitude π/2 of the observation point. When sampled over a wide longitudinal interval, the profile that is shown in Fig.  16 does not substantially differ from those of the pulses shown in other figures. The exceptionally high peak intensity and narrow width of the pulse at longitude 210° of this profile, i.e. what is depicted in part a of Fig.  17, show up only when the exact position of this peak is resolved. Here, we have set |${\hat{R}}_P= 10^{13}$| and have determined the value φP = 134.80887902374020168766219° (where n = 2 in equation 93) 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 Fig.  17(a) but also another as narrow and as intense pulse at a longitudinal distance of 4 × 10−18 deg from it. The 180° separation of the left- and right-hand components of the pulse profile shown in Fig.  16 stems from the coincidence in this example of the limiting values of θP1S and θP2S for RP → ∞.

The features exhibited by the pulse in Figs  16 and  17 can be inferred from equations ( 147), ( 166), and ( 167) also analytically. The only variables in the expression for Euc that depend on the observer’s longitude are σl2 and |${\bar{\sigma }}_{l2}$|⁠ , which vary linearly with φP (see equations 123 and  124 and note that, according to equation 93, φP drops out of the expressions for σl1 and |${\bar{\sigma }}_{l1}$|⁠ ). The variables σl2 and |${\bar{\sigma }}_{l2}$|⁠ , on the other hand, appear in equations ( 166)–( 167) only in the combinations |$\sigma _{l2}/\sigma _{l1}^3$| and |${\bar{\sigma }}_{l2}/{{\bar{\sigma }}_{l1}}^3$|⁠ . Hence, in cases where the turning points of the phase functions are sufficiently close to one another for σl1 or |${\bar{\sigma }}_{11}$| to be appreciably smaller than 1, the arguments of the functions |${\cal G}^{\rm r}_i$| that appear in the expression for Euc are highly sensitive functions of φP and so the pulse component that arises from two nearly coalescent critical points of flC or |${\bar{f}}_{lC}$| is exceedingly narrow (see Fig.  17a). Not only the widths but, as indicated by equations ( 150)–( 153) and ( 161)–( 165), also the amplitudes of |${\cal G}^{\rm r}_i$|⁠ , and hence of corresponding pulse components, vary sharply with φP when σl1 or |${\bar{\sigma }}_{11}$| assumes values that are close to zero.

The brightness temperatures implied by equation ( 176) and κu = 107 are listed in Table  1 for the pulses depicted in Figs  8, 11, and  16 and for a pair of examples of the pulses that are detected at the critical colatitudes |$\lim _{R_P\rightarrow \infty }\theta _{PlS}$|⁠ . Table  1 also shows the full width at half-maximum δφP of the components of the listed pulses that have the highest peaks (see Fig.  17a). Note that, as indicated by the last column of Table  1, the pulse profiles depicted in Figs  8– 16 have to be plotted on considerably shorter longitudinal scales before their peaks assume the shape shown in Fig.  17(a) and the maximum values of their dimensionless intensity |${\hat{I}}$| can be discerned graphically (see Section  5.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. Once resolved, the longitudinal distributions of the narrow pulses that are observed at critical colatitudes all have the same shape as that of the pulse shown in Fig.  17(a).

Brightness temperature Tb and full width at half-maximum δφP of the pulse detected at colatitude θP for the inclination angle α and |${\hat{R}}_P=10^{13}$|⁠ , κu = 107. The dimensionless factor |${\hat{T}}_{\rm b}$|⁠ , defined by equations ( 174)–( 176), is of the order of unity in the case of most radio pulsars. The limiting values of the critical colatitudes θP1S and θP2S in the second column are 48.533945294618400228° and 33.932818533330613261°, respectively.

Brightness temperature Tb and full width at half-maximum δφP of the pulse detected at colatitude θP for the inclination angle α and |${\hat{R}}_P=10^{13}$|⁠ , κu = 107. The dimensionless factor |${\hat{T}}_{\rm b}$|⁠ , defined by equations ( 174)–( 176), is of the order of unity in the case of most radio pulsars. The limiting values of the critical colatitudes θP1S and θP2S in the second column are 48.533945294618400228° and 33.932818533330613261°, respectively.

It should be noted that evaluating the critical colatitude |$\theta _{PlS}(\alpha ,{\hat{R}}_P)$| at |${\hat{R}}_P=\infty$| is not the same as placing the observation point at infinity. When one moves away from the source along a critical colatitude that is evaluated at |${\hat{R}}_P=\infty$|⁠ , the decreasing separation between the nearby stationary points of the corresponding phase function reduces to zero at infinity. Moreover, the separation between the stationary points in question turns out to have a very small value at all distances |${\hat{R}}_P$| along the latitudinal direction |$\lim _{R_P\rightarrow \infty }\theta _{PlS}$|⁠ . Since the coherence of the present radiation stems from the closeness of two stationary points of a phase function, the brightness temperature Tb in equation ( 173) thus has a high value in the latitudinal direction |$\lim _{R_P\rightarrow \infty }\theta _{PlS}$| at any distance from the source.

Values of Tb higher than those listed in Table  1 are predicted by equation ( 176) when the colatitude of the observation point lies closer to one of the critical angles θPlS or |${\bar{\theta }}_{PlS}$|⁠ . In the case of α = 60°, for example, the listed value of Tb (⁠|$8.67\times 10^{40}{\hat{T}}_{\rm b}\, {}^\circ$| K) corresponds to |$\theta _P=\lim _{R_P\rightarrow \infty }\theta _{P1S}=\lim _{R_P\rightarrow \infty }\theta _{P2S}=\pi /2$|⁠ . For an observation point whose colatitude is closer to the critical angle |$\theta _{PlS}(\alpha ,{\hat{R}}_P)$| for α = 60°, |${\hat{R}}_P=10^{13}$|⁠ , e.g. for θP = θP1S + 10−20 rad, Tb and δφP have the values |$8.67\times 10^{54}{\hat{T}}_{\rm b}\, {}^\circ$| K and 2.26 × 10−34 deg, respectively.

Given that the radiation field Euc depends on the observation time tP only in the combination φP − ωtP, the frequency spectrum of the present radiation is equally well described by the Fourier decomposition of the field Euc with respect to the azimuthal coordinate φP of the observation point. On the other hand, since the argument of the Dirac delta function in equation ( 77) depends on φP linearly, our replacing that delta function by its Fourier representation in equation ( 83) is tantamount to Fourier analyzing the fluctuations of the radiation field with respect to φP. The spectral distribution |${\bf E}_\nu ^{\rm uc}$| of the radiation field Euc can be inferred, therefore, from the k dependence of the integrand that appears in equation ( 133) and the relationship ν = kω/(2π) between k and the frequency of the radiation ν.

The variables σl1 and σl2 and the vector functions Pl and Ql that appear in equation ( 177) are independent of k. The step function H(k2 − k) in equation ( 131), moreover, can be set equal to 1 as pointed out in Section  4.6. So, when |Pl| ≫ k−1/2|Ql|, the vector |${{\boldsymbol {\cal K}}}_l$| and hence |$\boldsymbol {{\cal P}}_l$| and |$\boldsymbol {{\cal Q}}_l$| are also independent of k and the contributions of the terms containing them are by a factor of the order of k1/2 larger than those of the terms containing |${{\boldsymbol {\bar{\cal P}}}}_l$| and |${{\boldsymbol {\bar{\cal Q}}}}_l$| (see equations 127– 132). In this case, the possible values of the spectral index β are determined by the relative magnitudes of |$\vert \boldsymbol {{\cal P}}_l\vert$| and |$\vert \boldsymbol {{\cal Q}}_l\vert$| only. If |$\vert \boldsymbol {{\cal P}}_l\vert \gg k^{-1/3}\vert \boldsymbol {{\cal Q}}_l\vert$|⁠ , then β = 2/3 when the Airy functions in the first square bracket in equation ( 177) are of the order of unity and β = 1 when these Airy functions have the limiting values given by equations ( 179) and ( 180), and so the first square bracket in equation ( 177) decays as k−1/6. If |$\vert \boldsymbol {{\cal P}}_l\vert \ll k^{-1/3}\vert \boldsymbol {{\cal Q}}_l\vert$|⁠ , 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 equation ( 177) decays as k−1/6. In the opposite regime |Pl| ≪ k−1/2|Ql|, the factor k−1/2 multiplying Ql in equations ( 131) and ( 132) increases the value of the spectral index β by 1 everywhere (see Table  2).

Values, in various regimes, of the spectral index β defined in equation ( 181).

Values, in various regimes, of the spectral index β defined in equation ( 181).

The content of the spectrum in equation ( 177) stems from two types of fluctuations distinguished by the value of σl1 (or |${\bar{\sigma }_{l1}}$|⁠ ). One are the fluctuations arising from the isolated critical points of the phase function flC (or |${\bar{f}}_{lC}$|⁠ ) for which σl1 (or |${\bar{\sigma }}_{l1}$|⁠ ) is of the order of, or greater than, unity. In the time domain, these are manifested not only in the longitudinal width of a pulse but also in the sharp small-amplitude modulations (microstructure) of its profile (see Figs  8, 9, and  12): both the widths of the pulse components and the wavelengths of such modulations are proportional to |$\kappa _{\rm u}^{-1}$|⁠ . The other type of fluctuation is that arising from the imminent coalescence of two stationary points of a phase function and is manifested in the narrow width at half-maximum (δφP) of the pulse that is detected at a critical colatitude θlS (or |${\bar{\theta }}_{lS}$|⁠ ). As indicated by Fig.  17(a), the fraction δφP/2π of a rotation period during which such narrow pulses propagate past a detector is by many orders of magnitude smaller than unity when σl1 (or |${\bar{\sigma }}_{l1}$|⁠ ) is small (see Table  1). While the Fourier decomposition of the first type of fluctuations yields a spectrum centred on radio frequencies in cases where κu ≳ 107, that of the second type of fluctuation yields a wide spectral distribution extending from radio to gamma-ray frequencies: the width δφP = 1.21 × 10−26 rad of the pulse depicted in Fig.  17(a) translates into a frequency spectrum that extends as far as |$\nu \simeq \omega /(2\pi \delta \varphi _P)\simeq 1.31\times 10^{27} {\hat{P}}^{-1}$| Hz (see Fig.  18).

The Stokes parameters I, V and the position angles ψ of the polarization modes P and Q at an observation point with the colatitude θP = 2.5° for the inclination angle α = 5° and κu = 104. In this case, only the stationary points of f2C and |${\bar{f}}_{2C}$| contribute towards the field. Note the S-shaped swing through 180° of the position angle of the Q polarization mode across the pulse profile. The pulse window encompasses another similar pulse at a longitudinal distance of about 140° from this one.

Fig.  18 exemplifies the on-pulse spectrum of the radiation detected close to a critical colatitude. The peak amplitudes of the oscillations of this spectrum, which decrease as ν−2/3 with increasing frequency, are manifested as a set of discrete frequency bands.

Given that the sensitivity of any detector with which they are observed is limited to a finite bandwidth, the time-domain pulse profiles depicted in the above figures do not have the same shapes in the frequency domain. The range of frequencies over which a given component of a time-domain pulse profile is observable depends on the full width at half-maximum of that pulse component: a width illustrated in Fig.  17 that can be determined by the procedure described in Section  5.1.

The full width at half-maximum of the two higher peaks of the pulse depicted in Fig.  16 (those at longitudes 30° and 210°) is considerably shorter than that of the two lower peaks of this pulse (at longitudes 15° and 195°): while δφP = 1.19 × 10−29 rad for the higher peaks when κu = 1010, the corresponding width of the lower peaks is δφP = 6 × 10−8 rad. Thus, the two higher peaks of this time-domain pulse (which are detectable at frequencies around |$\omega /(2\pi \delta \varphi _P)=1.33\times 10^{30}{\hat{P}}^{-1}$| Hz) lag its two lower peaks (which are detectable at frequencies around |$2.65\times 10^8{\hat{P}}^{-1}$| Hz) by 15°.

Fig.  19 (in which |${\hat{R}}_P=10^{13}$|⁠ , |$\theta _P=\lim _{R_P\rightarrow \infty }\theta _{P2S}$|⁠ , α = 50°, κu = 1010, and only the stationary points of f2C and |${\bar{f}}_{2C}$| contribute towards the field) shows another example of such a phase lag, one in which the radio pulse is double-peaked but the gamma-ray pulse has a single component. In Fig.  19, we have plotted the contributions towards the Stokes parameter I that arise from the two isolated stationary points of the function |${\bar{f}}_{2C}$| and from the two nearly coalescent stationary points of f2C separately: the former in part (a) and the latter in part (b) of the figure. Once resolved, the two radio peaks at longitudes 30° and 95° (in Fig.  19a) each have the width δφP = 3.36 × 10−9 rad while the single gamma-ray peak at longitude 112° of this profile (in Fig.  19b) has the width δφP = 1.83 × 10−29 rad. The phases of the adjacent radio and gamma-ray peaks are thus separated by 17° in this case.

A phase lag between the peak of a narrow high-frequency pulse and that of an associated lower frequency pulse occurs whenever the radiation is observed at or near one of the critical latitudes. This can be seen from the relations |$k={\textstyle \frac{3}{2}}\sigma _{l1}^{-3}\kappa$| and κ ≤ κu (encountered in Sections  4.6 and  4.7) by noting that the frequency ν = ωk/(2π) associated with a given peak of the pulse profile sensitively depends on the value of σl1, i.e. on the nature of the critical point of the phase function that gives rise to the constructive interference of the waves embodying that pulse peak. For a pulse peak that arises from the imminent coalescence of two stationary points of the phase function flC, the variable σl1 is vanishingly small and so the value of the frequency ν is correspondingly high, irrespective of the value of κu. On the other hand, for a pulse peak that arises from an isolated stationary point of flC the variable σl1 is of the order of or greater than unity and so the maximum value of ν is determined by the upper limit κu on the value of κ, i.e. by the scale of fluctuations of the microstructure of the pulse profile. In the case depicted in Fig.  19, for example, the high-frequency peak at longitude 112° arises from the coalescing maximum and minimum of the phase function f2C while the low-frequency peaks at longitudes 30° and 95°, respectively, arise from the maximum and minimum of |${\bar{f}}_{2C}$|⁠ . The frequency ν at which the low-frequency peaks are detectable is determined not by the value of |${\bar{\sigma }}_{21}$| but by the value of κu.

The low-frequency and high-frequency peaks of a pulse profile would coincide only if κuω/(2π) has a value comparable to that of the higher frequency. Both types of peaks could then arise from the same isolated critical point of a phase function and share a broad spectrum.

Note that the linear extents in the latitudinal direction, RPδθP, of the focused radiation beams that embody the high-frequency radiation are of the order of the light-cylinder radius, c/ω, in general: these beams remain fully in focus at all distances |${\hat{R}}_P$| over the latitudinal interval δθP ≃ |θ − θPlS|, which turns out to be of the order of |${\hat{R}}_P^{-1}$| independently of the values of the other parameters (see Fig.  20). In the case of α = 60°, θP = 90°, κu = 107, and D = 1 kpc depicted in Figs  16 and  17, for example, the flux density S has the value |$32.1\, {\hat{S}}$| erg cm2 s−1. At latitudes closer to, or further away from, the critical angle for this example (⁠|$\lim _{R_P\rightarrow \infty }\theta _{PlS}=90^\circ$|⁠ ), the degree of focusing of the radiation beam and so the value of the flux density S are, respectively, higher or lower.

As pointed out in Section  4.4, the length of the interval |τlmax − τlmin| separating the τ coordinates of the maximum and minimum of the phase function flC decreases as |${\hat{R}}_P^{-1/2}$| with increasing |${\hat{R}}_P$| when this interval is small, i.e. when the colatitude θP of the observation point lies within an interval of length |${\hat{R}}_P^{-1}$| of the critical angle θPlS(α, L) with |$L \gt {\hat{R}}_P$| (see Fig.  20). In particular, if the observation point has the colatitude |$\lim _{R_P\rightarrow \infty }\theta _{PlS}$| (or |$\pi -\lim _{R_P\rightarrow \infty }\theta _{PlS}$|⁠ ), then the maximum and minimum of flC coalesce into an inflection point only at |${\hat{R}}_P\rightarrow \infty$|⁠ , rather than at a finite distance L. (These statements hold true also when flC and θPlS are replaced by |${\bar{f}}_{lC}$| and |${\bar{\theta }}_{PlS}$|⁠ , respectively.) In the case illustrated in Figs  16 and  17, for example, the colatitude of the observation point equals |$\lim _{R_P\rightarrow \infty }\theta _{P1S}=90^\circ$|⁠ , so that at the finite distance |${\hat{R}}_P=10^{13}$|⁠ , the τ coordinates of maximum and minimum of f1C are separated by the short interval 3.05 × 10−5 deg. It follows from the expression for ∂f1C/∂τ in equation ( A1) that this separation has the value |$3.05\times 10^{-5}({\hat{R}}_P/10^{13})^{-1/2}$|  deg for all |${\hat{R}}_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 flC and |${\bar{f}}_{lC}$| results in a decay rate of the flux density with distance that is slower than that predicted by the inverse-square law. Along colatitudes close to θPlS or |${\bar{\theta }}_{PlS}$|⁠ , the flux density S of the radiation diminishes with increasing distance from its source as |${\hat{R}}_P^{-3/2}$| instead of |${\hat{R}}_P^{-2}$|⁠ . This dependence of S on |${\hat{R}}_P$|⁠ , or equivalently D, is illustrated in Fig.  21 in the case where α = 60°, θP = 90°, κu = 107, and D ranges from 1 to 107 kpc, i.e. from a galactic to a cosmological distance.

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 transient. 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 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 Ardavan 2019, Appendix C).

The results reported in this paper are consequences mainly of the shape and motion of the current sheet: two features of the magnetosphere, which are the same not only for a wide class of magnetic-field configurations and plasma types (Bogovalov 1999; Beskin 2018) but also for both close to and far from the light cylinder. These results are found by inserting the charge and current densities (equations 18 and  19) that follow from the semi-analytic description of the magnetosphere (Tchekhovskoy et al. 2016) in the classical expression for the retarded potential (equation 31), thereby evaluating the radiation field (equation 32) in the time domain analytically (Section  4). Characteristics of the resulting radiation, i.e. its pulse profile and polarization position-angle distribution (Section  5.1), its brightness temperature (Section  5.2), its frequency spectrum (Section  5.3), the phase lag between its low-frequency and high-frequency peaks (Section  5.4), and its flux density (Section  5.5), were then obtained with the aid of its Stokes parameters and Fourier representation.

The examples in Figs  8– 16 illustrate the longitudinal distributions of the obtained pulse profiles for a range of values of the inclination angle α and observer’s colatitude θP. The electric field of the radiation generated by the current sheet turns out to be the sum of two distinct parts with differing polarization position angles (Section  4.6). Figs  8– 16 also show 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 multivalued function ψ in equation 170 are adopted that yield continuous position-angle distributions across various components of a given pulse.) Depending on the values of α and θP, which determine the number and positions of the stationary points of various phase functions, the pulse profiles can be peaked singly (Fig.  9), doubly (Figs  8, 10, 11, and  13), or multiply (Figs  14, 15, and  16), can have both narrow (Figs  8, 9, 13, and  16) and wide (Figs  12, 14, and  15) widths, and can comprise components that lie in either one (Figs  8, 10, and  12) or two widely separated (Figs  9, 11, and  16) longitudinal intervals.

The Stokes parameters I, V and the position angles ψ of the polarization modes P and Q at an observation point with the colatitude θP = 20° for the inclination angle α = 5° and κu = 102. In this case, the stationary points of |${\bar{f}}_{2C}$| alone contribute towards the field. This is the only pulse occurring in the entire pulse window.

The Stokes parameters (I, V, L) and the position angles ψ of the polarization modes P and Q at an observation point with the colatitude θP = 55° for the inclination angle α = 45° and κu = 102. In this case, only the stationary points of f2C and |${\bar{f}}_{2C}$| contribute towards the field. Note that the polarization of this pulse changes from linear to circular across it and the position angle of the Q mode swings through 180° across each one of its two components. The pulse window encompasses another pulse at a longitudinal distance of about 115° from this one.

The Stokes parameters (I, V, L) and the position angles ψ of the polarization modes P and Q at an observation point with the colatitude θP = 35° for the inclination angle α = 45° and κu = 102. In this case, only the stationary points of f2C and |${\bar{f}}_{2C}$| contribute towards the field. Position angle of the Q mode turns through 360° across the narrow gap at longitude 70° in the curve depicting its distribution. These are the only pulses occurring in the entire pulse window.

The degree of linear polarization is considerably higher than that of circular polarization for most of these examples. There are, however, cases such as that depicted in Fig.  11 for which the Stokes parameter V is comparable to, or greater than, the Stokes parameter L over parts of the pulse window. As can be seen from Figs  8– 10 and  14– 16, the sense of circular polarization is often reversed across certain components of the pulse profile.

S-shaped swings through 180° are frequently exhibited by the longitudinal distributions of the polarization position angles of either one (Figs  8 and  9) or both (Figs  11– 15) polarization modes. But there are cases in which the position angles of either one or both modes remain constant across the pulse profile (Figs  8– 10). Position angles of the two concurrent modes can be approximately orthogonal over a section of the pulse profile (Figs  8, 9, 12, 14, and  15) or can be coincident (Figs  10 and  13).

The Stokes parameters (I, V, L) and the position angles ψ of the polarization modes P and Q at an observation point with the colatitude θP = 30° for the inclination angle α = 65° and κu = 102. In this case, only the stationary points of f2C and |${\bar{f}}_{2C}$| contribute towards the field. Note the simultaneous swing through 180° of the position angles of both modes across the right-hand component of the pulse. Note also the high degree of circular polarization of the pulse throughout the depicted longitudes. The pulse window encompasses in addition a weaker pulse at a longitudinal distance of about 20° from these ones.

The Stokes parameters (I, V, L) and the position angles ψ of the polarization modes P and Q at an observation point with the colatitude θP = 77.5° for the inclination angle α = 65° and κu = 107. In this example of a multipeak pulse profile the stationary points of all of the phase functions (f1C, f2C, |${\bar{f}}_{1C}$|⁠ , and |${\bar{f}}_{2C}$|⁠ ) contribute towards the field. Position angle of the Q mode turns through 360° across the narrow gap at longitude 100° in the curve depicting its distribution.

The Stokes parameters (I, V, L) and the position angles ψ of the polarization modes P and Q at an observation point with the colatitude θP = 100° for the inclination angle α = 80°, κu = 104 (pulse profile), and κu = 107 (position angle). In this case, the stationary points of all of the phase functions (f1C, f2C, |${\bar{f}}_{1C}$|⁠ , and |${\bar{f}}_{2C}$|⁠ ) contribute towards the field. This is an example of a multicomponent pulse the polarization position angles of whose different components have differing longitudinal variations.

Though occurring in every case, the rapid small-amplitude modulations of the pulse profiles and position-angle distributions are visible only in Figs  8– 10 and  12 for which κu ≤ 104. The wavelength of such modulations is proportional to |$\kappa _{\rm u}^{-1}$|⁠ , so that the microstructure they superpose on pulse profiles and position-angle distributions is too sharp and dense to show up in the rest of the plotted figures. The free parameter |$\kappa _{\rm u}^{-1}$| cannot of course be smaller than the (unknown) thickness of the current sheet in units of the light-cylinder radius (see Section  4.7).

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.

Figs  16 and  17 illustrate an example of a radically different type of pulse: one detectable near those observation points for which two nearby stationary points of one of the phase functions in the expression for the radiation field coalesce, thus giving rise to a much tighter focusing of the emitted waves. For any given values of the inclination angle α and the observer’s distance RP, there are at least two, and at most eight, critical colatitudes (θPlS and/or |${\bar{\theta }}_{PlS}$|⁠ ) at which two nearby stationary points of one of the phase functions (flC and/or |${\bar{f}}_{lC}$|⁠ ) coalesce to give rise to a higher-order focusing of this kind. Though their profiles over the entire pulse window look similar to those of other pulses (Fig.  16), such pulses display extraordinarily large amplitudes and short widths once they are inspected over sufficiently short longitudinal intervals to resolve their peaks (Fig.  17a). As pointed out in Section  1, the extraordinary values of the amplitudes and widths of such pulses, illustrated by the example in Figs  16 and  17(a), are what underpin the high brightness temperatures and broad frequency spectra of the radiation generated by the current sheet.

The Stokes parameters (I, V, L) and the position angles ψ of the polarization modes P and Q at an observation point with the coordinates |${\hat{R}}_P=10^{13}$| and θP = 90° for the inclination angle α = 60° and κu = 107. In this case, the stationary points of all of the phase functions (f1C, f2C, |${\bar{f}}_{1C}$|⁠ , and |${\bar{f}}_{2C}$|⁠ ) contribute towards the field. While the two higher peaks of this profile (at longitudes 30° and 210°) arise from nearly coalescent stationary points of f2C, the lower peaks (at longitudes 15° and 195°) arise from two isolated stationary points of |${\bar{f}}_{2C}$|⁠ . The full widths at half-maxima of the higher and lower peaks of this pulse profile can be inferred from Figs  17(a) and (b), respectively.

In part (a) of this figure, the width of the component at longitude 210° of the distribution of the Stokes parameter I that is depicted in Fig.  16 is resolved by plotting this distribution over successively shorter longitudinal intervals centred at its peak until the value of I at its maximum stops growing. In part (b), the same procedure is applied to the component at longitude 195° of the distribution of the Stokes parameter I depicted in Fig.  16. Besides the difference in width, there is a generic difference between the shape of a resolved pulse that arises from two nearly coalescent stationary points of a phase function (as in a) and that of a resolved pulse that arises from an isolated stationary point of a phase function (as in b).

By equating the magnitude of the Poynting flux of the radiation to the Rayleigh–Jeans law for the energy that a blackbody of the same temperature would emit per unit time per unit area into a given frequency band, it was shown in Section  5.2 that the brightness temperature of the pulsed radiation depicted in Figs  16 and  17(a) has the value |$8.67\times 10^{40}{\hat{T}}_{\rm b}~{}^\circ$| K, where the dimensionless factor |${\hat{T}}_{\rm b}$| is of the order of unity in the case of most radio pulsars. In contrast, the brightness temperatures of the pulses depicted in Figs  11 and  8 have the values |$1.73\times 10^{18}{\hat{T}}_{\rm b}~{}^\circ$| K and |$9.56\times 10^{22}{\hat{T}}_{\rm b}~{}^\circ$| K, respectively (Table  1). There are also conditions under which the pulse generated by the current sheet is so focused as to attain a brightness temperature of the order of |$10^{55}{\hat{T}}_{\rm b}~{}^\circ$| K (see Section  5.2).

Given that the radiation field of the current sheet depends on the observation time tP and the azimuthal coordinate φP of the observation point only in the combination φP − ωtP, the frequency spectrum of this radiation is equally well described by the Fourier decomposition of its longitudinal distribution. Depending on the colatitude of the observation point, the content of this spectrum stems either from the fluctuations manifested in the microstructure of the pulse profile or from the much sharper variations characterized by the full width at half-maximum (δφP) of the pulse that arises from the imminent coalescence of two stationary points of a phase function (Table  1). In the high-frequency limit, the spectral flux density Sν of the radiation 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, and 7/3 (Table  2). On the other hand, at observation points for which two stationary points of a phase function in the expression for the radiation field are nearly coincident, the narrow width of the peak of the pulse profile is reflected in a broad spectrum that extends from radio waves to gamma rays: the width δφP = 1.21 × 10−26 rad of the pulse depicted in Fig.  17(a), for example, implies a frequency spectrum that extends as far as ν ≃ ω/(2πδφP) ≃ 1.31 × 1027 Hz when ω = 102 rad/s. Fig.  18 shows the on-pulse spectral distribution of this high-frequency radiation. The peak amplitudes of the oscillations of this spectrum, which decrease as ν−2/3 with increasing frequency, are manifested in a dynamic spectrum as a set of discrete frequency bands.

Distribution of the spectral flux density of the narrow pulse shown in Fig.  17(a). Thus, the dynamic spectrum of the radiation consists of a discrete set of bands whose amplitudes decrease with increasing frequency.

A phase lag between the peak of a high-frequency pulse and that of an associated lower-frequency pulse occurs whenever the radiation is observed at or near one of the critical latitudes (Section  5.4). The range of frequencies over which a given component of a time-domain pulse profile is observable depends on the width of that pulse component. The full width at half-maximum of each component of a time-domain pulse profile is in turn determined by the nature of the stationary point of the phase function that gives rise to the constructive interference of the waves embodying that component: the widths of the pulse components that arise from two nearly coalescent stationary points are by many orders of magnitude shorter than those of the components that arise from isolated stationary points. Figs  16, 17, and  19 illustrate this feature both for singly and doubly peaked pulse profiles and for differing values of the phase lag.

The Stokes parameters I and V for a three-component pulse profile in the time-domain. In this figure, |${\hat{R}}_P=10^{13}$|⁠ , |$\theta _P=\lim _{R_P\rightarrow \infty }\theta _{P2S}$|⁠ , α = 50°, and κu = 1010, and the origin of the horizontal axis in (b) is shifted by 111.9695512° relative to that of the horizontal axis in (a). While the two components shown in (a) can be detected only in the radio band, the single component shown in (b) is observable as a much shorter gamma-ray pulse (see Section  5.4). The radio and gamma-ray peaks of this pulse profile are thus separated by about 17°.

The separation between two nearly coincident stationary points of the integral over the latitudinal distribution of the current sheet in the expression for the radiation field decreases as |$R_P^{-1/2}$| with increasing distance RP of the observation point from the star (Fig.  20). The enhanced focusing of the radiation that is caused by this shortening of the separation between the stationary 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 change differently with RP, this effect results in a flux density for the radiation that diminishes with increasing distance as |$R_P^{-3/2}$| rather than |$R_P^{-2}$| (Fig.  21). As indicated by Fig.  20, the latitudinal width δθP of the radiation beam whose flux density decreases as |$R_P^{-3/2}$| with increasing |${\hat{R}}_P$| is of the order of (RPω/c)−1 rad. However, the gradual change in the rate of decay of the flux density with distance, from |${\hat{R}}_P^{-3/2}$| to |${\hat{R}}_P^{-2}$| away from a critical latitude, takes place over a latitudinal interval of the order of a radian (see equation 184). In as much as there are several critical latitudes for any given values of α and |${\hat{R}}_P$| (Section  4.4), the solid angle over which the non-spherically decaying pulses are observable, though limited, is finite.

Dependence of the separation |τlmax − τlmin| of the maximum and minimum of the phase function flC on |${\hat{R}}_P$| and on θP − θPlS for α = 5° and l = 2. The critical colatitude θPlS has the value 5.0190029119° in this case. This figure shows that |τlmax − τlmin| decreases as |${\hat{R}}_P^{-1/2}$| with increasing distance when |$(\theta _P-\theta _{PlS})\lesssim {\hat{R}}_P^{-1}$|⁠ .

|$log(S/{\hat{S}})$| versus log(D) for α = 60°, θP = 90°, and κu = 107. The blue line with the slope −3/2 is the best fit to the red dots whose ordinates are determined by evaluating equation ( 182). 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\rightarrow \infty }\theta _{PlS}$|⁠ , |$\lim _{R_P\rightarrow \infty }{\bar{\theta }}_{PlS}$|⁠ , |$\pi -\lim _{R_P\rightarrow \infty }\theta _{PlS}$|⁠ , or |$\pi -\lim _{R_P\rightarrow \infty }{\bar{\theta }}_{PlS}$|⁠ .

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 transient. 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 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 therefore balanced by the change with time of the energy contained inside the shell bounded by those spheres (see Ardavan 2019, Appendix C).

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 movement (e.g. precession or a glitch) of the star’s rotation or magnetic axes, the radiation beams embodying such pulses sweep past the Earth. Using the decay rate |$R_P^{-2}$|⁠ , instead of |$R_P^{-3/2}$|⁠ , we would overestimate the power emitted by the sources of the bursts of radiation we would receive in this way by as large a factor as 1015 if the neutron stars that generate the bursts lie at cosmological distances. It is widely maintained that the powers emitted by distant sources of gamma-ray and fast radio bursts are by many orders of magnitude greater than that emitted by a pulsar (Kumar & Zhang 2015; Petroff et al. 2019). The unquestioned assumption on which this consensus is based is that the radiation fields of all sources necessarily decay as predicted by the inverse-square law. This assumption is brought into question by the results of the present analysis, however. The putative difference between the energetic requirements of what are regarded to be different types of sources could arise from the difference in the latitudinal direction along which an obliquely rotating neutron star is observed.

It should be emphasized that in contrast to every one of the models currently considered in the literature on the subject (Beskin 2018; Melrose, Rafat & Masterano 2021), which are at least in part phenomenological, the all-encompassing body of explanations provided by the present analysis for the salient features of the observed emission from pulsars (Lyne & Graham-Smith 2012; Abdo et al. 2013) and for the putative energetic requirements of magnetars (Kaspi & Beloborodov 2017) and sources of fast radio bursts and gamma-ray bursts (Kumar & Zhang 2015; Petroff et al. 2019) is a consequence purely of the basic equations that govern the magnetospheric structure of a neutron star and the emission of electromagnetic waves.

A few final remarks concerning methodology are in order:

The thickness assigned to the current sheet by the description in Tchekhovskoy et al. ( 2016) is zero (Section  2). However, a superluminally moving source is necessarily volume-distributed (Bolotovskii & Bykov 1990): it would give rise to a divergent field if it has no thickness (Section  4.7). We have removed the singularity that arises from overlooking the finite width of the current sheet by setting a lower limit, |$\kappa _{\rm u}^{-1}$|⁠ , on the scale of the fluctuations associated with the microstructure of the pulse profile and treating κu as a free parameter (Section  4.7). The thickness of the present 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 analyzing a current sheet (Harris 1962) 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.

Certain characteristics of the pulses described here (e.g. their widths) would be modified in the course of their propagation through the plasma between the current sheet and the interstellar medium. However, the effect of this plasma on the present radiation is no different from its effect on any other coherent radiation. The mechanisms by which the emissions commonly considered in the literature on pulsars achieve their coherence are different from that described in the present paper (Beskin 2018; Melrose et al. 2021) but, once generated, every one of these emissions is subjected to the same propagation effects before it reaches the interstellar medium. The non-spherically decaying 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 (Philippov & Spitkovsky 2018; Philippov et al. 2019).

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 Ardavan 2019, Section  3).

Since the problem discussed in this paper entails the formation of caustics, we cannot proceed to the far-field limit |xP| → ∞ before evaluating the radiation field, as is customarily done in radiation theory. The far-field approximation of the argument of the delta function in equation ( 31) 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. Achkasov & Zhuravlev 2020). The stringent requirements set by constraint (35) in the superluminal regime and further technical reasons why a conventional approach to the present problem does not work are discussed in Appendix B of Ardavan ( 2019).

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 (Spitkovsky 2006; Kalapotharakos et al. 2012). Irrespective of the formalism on which they are based (whether force-free, magnetohydrodynamic, 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 Section 3 of Ardavan ( 2019), 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 type of 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 Ardavan 2019, Section  3).

No new data were generated or analyzed in support of this research.

Abdo A. A. et al.  , 2013, ApJS, 208, 17 10.1088/0067-0049/208/2/17

Achkasov V. V. , Zhuravlev M. Y. , 2020, Rep. Math. Phys., 85, 375 10.1016/S0034-4877(20)30042-2

Ardavan H. , 2019, J. Plasma Phys., 85, 905850304 10.1017/S0022377819000382

Beskin V. S. , 2018, Phys.-Usp., 61, 353 10.3367/UFNe.2017.10.038216

Bleistein N. , Handelsman R. A. , 1986, Asymptotic Expansions of Integrals. Dover, New York

Bogovalov S. V. , 1999, A&A, 349, 1017

Bolotovskii B. M. , Bykov V. P. , 1990, Sov. Phys.-Usp., 33, 477 10.1070/PU1990v033n06ABEH002601

Bolotovskii B. M. , Ginzburg V. L. , 1972, Sov. Phys.-Usp., 15, 184 10.1070/PU1972v015n02ABEH004962

Burridge R. , 1995, SIAM J. Appl. Math., 55, 390

Chester C. , Friedman B. , Ursell F. , 1957, Proc. Cambridge Philos. Soc., 53, 599 10.1017/S0305004100032655

Erdélyi A. , ed., 1954, Tables of Integral Transforms. Vol. 1. McGraw-Hill, New York

Ginzburg V. L. , 1972, Sov. Phys.-JETP, 35, 92

Hadamard J. , 2003, Lectures on Cauchy’s Problem in Linear Partial Differential Equations. Dover, New York

Harris E. G. , 1962, Nuovo Cim., 23, 115 10.1007/BF02733547

Hoskins R. F. , 2009, Delta Functions: An Introduction to Generalised Functions, 2nd edn. Woodhead, Oxford

Jackson J. D. , 1999, Classical Electrodynamics, 3rd edn. Wiley, New York

Kalapotharakos C. , Contopoulos I. , Kazanas D. , 2012, MNRAS, 420, 2793 10.1111/j.1365-2966.2011.19884.x

Kaspi V. M. , Beloborodov A. M. , 2017, ARA&A, 55, 261 10.1146/annurev-astro-081915-023329

Kumar P. , Zhang B. , 2015, Phys. Rep., 561, 1 10.1016/j.physrep.2014.09.008

Lyne A. , Graham-Smith F. , 2012, Pulsar Astronomy, 4th edn. Cambridge Univ. Press, Cambridge

Melrose D. B. , Rafat M. Z. , Masterano A. , 2021, MNRAS, 500, 4530 10.1093/mnras/staa3324

Olver F. W. J. , Lozier D. W. , Boisvert R. F. , Clark C. W. , eds, 2010, NIST Handbook of Mathematical Functions. Cambridge Univ. Press, Cambridge

Petroff E. , Hessels J. W. T. , Lorimer D. R. , 2019, A&A Rev., 27, 4 10.1007/s00159-019-0116-6

Philippov A. A. , Spitkovsky A. , 2018, ApJ, 855, 94 10.3847/1538-4357/aaabbc

Philippov A. , Uzdensky D. A. , Spitkovsky A. , Cerutti B. , 2019, ApJ, 876, L6 10.3847/2041-8213/ab1590

Tchekhovskoy A. , Spitkovsky A. , Li J. C. , 2013, MNRAS, 435, L1 10.1093/mnrasl/slt076

Tchekhovskoy A. , Philippov A. , Spitkovsky A. , 2016, MNRAS, 457, 3384 10.1093/mnras/stv2869

Oxford University Press is a department of the University of Oxford. It furthers the University's objective of excellence in research, scholarship, and education by publishing worldwide

Sign In or Create an Account

This PDF is available to Subscribers Only

For full access to this pdf, sign in to an existing account, or purchase an annual subscription.