Frequency resolved radio and high-energy emission of pulsars

Context. Pulsars are detected as broadband electromagnetic emitters from the radio wavelength up to high and very high energy in the MeV/GeV and sometimes even in the TeV range. Multi-wavelength phase-resolved spectra and light curves offer an unrivalled opportunity to understand their underlying radiation mechanisms and to localize their emission sites and therefore the particle acceleration regions. Aims. In this paper, we compute pulsars multi-wavelength phase-resolved light-curves and spectra assuming that curvature radiation operates from inside the magnetosphere of a rotating vacuum magnet. Radio emission arises from the polar caps whereas gamma-ray energy emanates from the slot gaps in the vicinity of the separatrix between closed and open field lines. Methods. By integrating particle trajectories within the slot gaps, we compute energy dependent photon sky maps in the radio band (MHz-GHz) and in the gamma-ray band (MeV/GeV) for mono-energetic distribution functions of leptons. Results. We obtained many details of the energy dependent light curves and phase-resolved spectra from the radio wavelength up to the gamma-ray energies. Choosing Lorentz factors of {\gamma} = 30 for radio emitting particles and {\gamma} = 10^7 for gamma-ray emitting particles, we found realistic spectra accounting for the wealth of multi-wavelength pulsar observations.


Introduction
Neutron stars are mainly observed as pulsars, emitting regular pulses in the whole electromagnetic spectrum from radio, through optical, X-rays up to gamma-rays. Although it is well established celerating electric field. They used test particle integration in the radiation reaction limit regime in the global force-free dissipative magnetosphere which is basically a fluid description avoiding the stringent strong field constrain faced by PIC codes. In such a way, they were able to deduce spectra for realistic pulsar field strengths and periods. Nevertheless, starting from PIC simulations, Kalapotharakos et al. (2018) recently found a relation between the particle injection rate and the spin down luminosity. This work shows the fruitful feedback between simulations and observations to extract useful information about the nature of particle acceleration and dissipation of the relativistic magnetized flow. Harding (2016) and Venter et al. (2018) gave recent reviews of the successful interplay between magnetospheric modelling and gamma-ray observations. Further attempts to fit particular pulsars were carried out by other groups. For instance Takata et al. (2007) and Hirotani (2008) used the vacuum retarded dipole to model the outer gap of the Crab pulsar. Du et al. (2011) performed computation in the annular gap context for the Vela pulsar whereas Du et al. (2012) did it for the Crab pulsar. Several millisecond pulsars were also fitted by Du et al. (2013) using a static dipole. Recently the impact of general-relativistic effects on pulsar light-curves have been investigated by Pétri (2018) followed by some improvement made by Giraud & Pétri (2020), including the Shapiro delay.
In these aforementioned works, our attention was focused on the shape of the pulses emitted in radio and gamma-ray, explicitly showing the impact of gravitation on the geometry of the profiles and on the phase shift between the main radio pulse and the first gamma-ray pulse. We did not take into account the frequency characteristics of this emission, assuming it to be curvature radiation, which is a function of the point of emission and more particularly of the curvature of the field lines. But what about their frequency distribution, i.e. the radio and gamma-ray spectrum? How does the shape of the pulses evolve as a function of frequency or energy? This paper answers these questions by providing precise quantitative results if the emission comes from curvature radiation all along the magnetic field lines in a flat spacetime geometry. The outline of our paper is as follows. In section 2 we expose our model magnetosphere model and the associated radiation mechanism. Next, in section 3, we show detailed results of high-energy spectra and light-curves from the slot gap models, localising the sites of efficient gamma-ray photon production depending on the magnetic inclination angle and line of sight. In section 4, the same approach is used to study the radio spectra and light-curves depending on the geometry. Conclusions and possible extensions are summarised in section 5.

A coherent multi-wavelength model
Multi-wavelength modelling of radio-loud gamma-ray pulsars is usually limited to the joint computation of only one radio and one gamma-ray light curve. Although this strategy already makes it possible to constrain the geometrical parameters of the magnetosphere, this is hardly satisfactory Article number, page 3 of 39 because the pulse profiles evolve as a function of frequency or energy. A detailed multi-frequency investigation, including coherently the calculation of radio and gamma-ray spectra from, for example, the curvature radiation, proves to be much more constraining and will allow to physically and quantitatively link the geometry of the field lines, the altitude of the photon production sites, the distribution functions of the emitting particles as well as the associated light curves.
In this paper, this detailed view of broadband electromagnetic emission processes is discussed in depth. To do so, we start with a reminder of the method of calculation of the properties of the curvature emission, its power and its characteristic frequency, applied to the pulsar magnetosphere.
Then we will determine the high energy spectrum around the GeV, thus in the Fermi/LAT band, as well as some characteristic light curves to conclude on the spectrum and radio pulse profiles.

Magnetosphere structure
We use the analytical expression of the electromagnetic field of a rotating dipole in vacuum as first given by Deutsch (1955). Charged particles are accelerated in gaps where the plasma present in the pulsar's magnetosphere is not dense enough to screen the component of the electric field parallel to the magnetic field lines. In our model, those gaps are either located along the last closed magnetic field lines from the star surface up to the light cylinder (slot gaps, (Arons 1983)) or above the polar caps delimited by the last closed magnetic field lines (polar gaps, (Ruderman & Sutherland 1975)).
We also assume that radio and gamma-ray emission are produced by the acceleration of charged particles in those gaps through curvature radiation. Let us therefore remind the essential feature of this emission mechanism.

Curvature radiation
The calculation of the frequency of curvature radiation is deduced from the radius of curvature ρ of the particle trajectory along the magnetic field lines, which is equal to the radius of curvature of these field lines in the corotating frame. The characteristic curvature frequency ω curv is (Jackson 2001) beyond which the particle no longer radiates significantly. Specifically, the shape of the curvature radiation spectrum for a particle with Lorentz factor γ is given in Jackson (2001) by where K 5/3 is the modified Bessel function of order 5/3 (Arfken & Weber 2005).
A trajectory or line in space has a curvature as well as a twist. Locally, at a point of this curvature, we associate a trihedron (T, N, B) also called Frenet's frame. The curvature κ indicates the change in the tangent vector T and the torsion τ indicates the variation in the osculating plane Article number, page 4 of 39 at the same point. In summary the Frenet formulas giving the variations of the basis vectors of the trihedron are ds being the curvilinear abscissa along the trajectory. Torsion does not intervene in the calculation of the particle curvature radiation. The curvature itself is thus deduced from the derivative of the tangent T = B/B to the field lines as a function of the curvilinear abscissa s along these field lines.
The curvature radius ρ is the inverse of the curvature κ of the magnetic field which is the derivative of the tangent T to the magnetic field line with respect to the curvilinear abscissa ds, so we have We recall that the variation of the curvilinear abscissa ds is related to the variations of the coordinates of the trajectory dx i by ds 2 = γ ik dx i dx k in any spatial geometry given by the spatial metric γ ik which is the spatial projection of the space-time metric g ik with respect to a given observer.
In our model, we assume that the particles follow the field lines in the reference frame corotating with the star, the frequency of the curvature radiation ω c is therefore that of the rotating frame of reference (1). In order to find the frequency of the radiation emitted by our rotating pulsar and measured by an distant inertial observer at rest ω c , we must take into account the Doppler effect and perform a Lorentz transformation from the corotation frame to the inertial frame of the observer.
This is equivalent to multiplying the equation (1) by a Doppler factor η due to the neutron star rotation. We thus have ω c = η ω c with the Doppler factor generally put in the form n being the unit vector giving the initial direction of propagation of the photon (tangent to the magnetic field line) as seen in the observer's frame of reference and taking into account the aberration and β the corotation velocity. So n is given by the equation γ is the Lorentz factor of the particle accelerated in the electromagnetic field, primary or secondary beam. The power radiated by this particle is given by Article number, page 5 of 39 with e the elementary charge and 0 the vacuum permittivity.
For curves of known geometry, this curvature can be determined analytically, as for example for a static dipole for which the equation of the field lines is known, the lines being themselves each contained in a plane. For a rotating dipole such as the Deutsch dipole or in general relativity, it is not possible to obtain an analytical expression of this curve. It must be estimated numerically by approximating the derivative (4).
Concretely, this curvature is calculated by measuring the variation ∆T of the vector tangent T to the magnetic field lines as one moves along them by a value ∆s of the curvilinear abscissa s small enough that the finite difference represents a good approximation of the derivative. Typically a fraction of the radius of the star ∆s R * is sufficient. We then replace κ in the equation (4)  In order to minimize discretization errors, we chose a method with centered finite differences.
Knowing the local radius of curvature and the value of the magnetic field at the point of emission in the rotating frame of reference, we switch back to the spectrum of the curvature radiation in the inertial frame of reference by imposing an energy distribution function for the photon emitting leptons. The simplest situation consists in taking a mono-energy distribution of particles although a general particle population is possible.

Particle population
The shape of the light curves and spectra strongly depend on the energy distribution of the radiating particles. Current models of acceleration of primary particles and creation of secondary and higher generation of pairs (tertiary, quaternary, ...) predict a peak of the Lorentz factor around γ ≈ 10 7 for primaries and around γ ≈ 10 2 for others (Beskin et al. 1993). Although our approach is easily adapted to any particle distribution function, in this work we only consider mono-energetic particles either in the primary beam or in the secondary beam. The primary particles are responsible for the high-energy emission in the GeV band in the slot gap while the secondary particles from the cascades produce the radio emission from the polar caps.
In a more realistic picture, the particle distribution function significantly deviates from a monoenergetic one, the radiation being linearly additive, the calculation of the spectra and light curves will follow the same linear addition scheme of the individual mono-energetic contributions with respective weights corresponding to the weights imposed by the particle distribution function. For example, for a population containing N 1 particles of energy γ 1 and N 2 particles of energy γ 2 , the total emission will be N 1 times that of the mono-energy spectrum at γ 1 to which we superimpose N 2 times that of the mono-energy spectrum at γ 2 . A continuous distribution is then sub-divided into several intervals of mono-energetic particles and the entire spectrum can be reconstructed. Our Article number, page 6 of 39 results therefore serve as a building block for the most general distribution function obtained by a post-treatment of mono-energetic distribution functions.
Let us now review the results for the high-energy mono-energetic and radio spectra: for each of the two cases, we will assume that at each emission point, located along the last magnetic field lines with a curvilinear spacing ∆ between each point, there is emission by curvature radiation considered to emanate from particles with Lorentz factors γ fixed for each wavelength. We assume that there is no emission for points located more than 95% of R cyl from the centre of the star.
The classical estimation of the maximum Lorentz factor attainable by particles in the electromagnetic field takes into account their braking by radiation reaction in order to compensate the acceleration by the parallel electric field E present in the gaps. The precise value of this field is strongly dependent on the dynamics within the gaps. By equating the radiated power (7) and the power provided by this electric field e E c, the Lorentz factor becomes In our vacuum model, E ≈ 10 12 V/m corresponds to the upper limit, that of the Deutsch field where the electric field is not screened at all. In reality, this value is much lower but the maximum Lorentz factor only varies with E 1/4 , so it is not very sensitive to a large variation of this field.
It is currently well accepted that the pulsar magnetosphere contains at least two different populations of particles flowing out along open magnetic field lines. The first component is made of primary particles accelerated by the strong parallel electric field in vacuum gaps up to Lorentz factors about γ ≈ 10 7 and with a particle density number close to the corotation density of n p ≈ 2 ε 0 Ω · B/e. These primary particles emit gamma-ray photons decaying into a secondary flow made of electron-positron pairs (Erber 1966) with Lorentz factor in the range γ ≈ 10 2 − 10 3 and a multiplicity factor of κ ≈ 10 3 −10 5 (Daugherty & Harding 1982). The typical energy of curvature photons produced by each beam will therefore vary by several orders of magnitude. The radio emission is generated by the dense secondary beam whereas the high energy MeV/GeV emission is generated by the ultra-relativistic but much less dense primary beam.
In our model, we assumed a simple mono-energetic population, we typically set γ g = 10 7 for the gamma-rays emission coming from the slot gaps whereas for the radio emission from the polar caps we take γ r = 30.

Emissivity
First of all, the calculation of the spectra we present in the following does not take into account the range of frequencies at which photons are emitted. We approximate the spectrum by a Dirac distribution centred on the characteristic frequency ω c and whose intensity corresponds to the total power of the integrated radiation on all frequencies of the true continuous spectrum (7). In other Article number, page 7 of 39 words, the spectrum is approximated in the corotating frame by Each point on a field line emits photons at a single frequency given by the local characteristic frequency ω c (ρ) (1) associated with the local curvature of this same field line. The number of photons to consider is proportional to the total emitted power ω curv (ρ). A Lorentz transform bring all these quantities into the observer inertial frame. This contrasts sharply with our previous work (Giraud & Pétri 2020) where each position on a given field line produced only one photon of indeterminate energy. In this paper, we remove the ambiguity on frequency by incorporating details of the emission mechanism, frequency and power. Emissivity will no longer be constant along a field line but will vary according to the local curvature.
In a forthcoming step, we will consider using the full expression of the curvature radiation spectrum given by (2). Finally, in another possible extension, the most realistic situation takes into account a particle distribution function depicted by a power law such that the number of particles of Lorentz factor between γ and γ + dγ is given by where p is the spectral index of the power law, the minimum Lorentz factor being γ min and a the maximum value being γ max . This particle distribution function produces another power law for the curvature emissivity j curv such that There is thus a simple relationship between the spectral index of the power distribution of particles and photons. Beyond the cut-off frequency associated with γ max the emissivity falls exponentially according to e −ω/ω max c . Below the cut-off frequency associated with γ min the emissivity decreases according to another power law independent of the particle distribution and scaling as (ω/ω min c ) 1/3 . In between these two limits, the photon spectral index is related to particle power-law distribution p.

Aberration, retardation and seep-back effects
The radio and high-energy emission of pulsars originate at high altitude in the magnetosphere, and well above the surface for young pulsars of period P 100 ms (Mitra 2017). The corotation of these emission sites at relativistic velocities imparts unique characteristics to the pulsar profiles by impacting the relationship between the geometry of the emitting zones and their observational signature. Phillips (1992) has detailed the impact of several important effects on the arrival time of pulses of all kinds. Let us consider two emission zones located at an altitude respectively r 1 and r 2 > r 1 .
Article number, page 8 of 39 Q. Giraud and J. Pétri: Frequency resolved emission of pulsars For photons moving radially away from the neutron star, the delay introduced by the difference in path to be travelled becomes a phase shift for the observer The dragging of the photon production areas by the rotation of the pulsar causes a projection of the direction of propagation of the photons in the direction of rotation due to aberration by a value of For both emission sites, this introduces an additional time delay of Finally, the rotation of the dipole curves the field lines in the opposite direction of rotation. Shitov (1983) gave a simple expression for this sweep-back effect, such that the delay induced amounts to The total delay ∆t Σ is then the sum of each of these contributions ∆t Σ = ∆t a + ∆t r + ∆t B .
We now discuss in some details the high-energy and radio light-curves and spectra for monoenergetic particle distribution functions applied to the polar cap and slot gap models of a vacuum rotating dipole.

High-energy emission
We are interested in the high-energy emission as measured in the Fermi/LAT (Large Area Telescope) band, between 100 MeV and 100 GeV. This Fermi/LAT telescope has detected nearly 300 gamma-ray pulsars all showing very similar spectra peaking around several GeV (Abdo et al. 2013). Our long lasting aim is to compare our predicted light curves and spectral modelling with a whole population of young and millisecond radio-loud gamma-ray pulsars. This section details the bottom line to construct general light-curves and spectra from the slot gaps knowing the particle distribution function.

Thin slot gap
The locus of the last open field lines taking their root at the polar caps generates a two-dimensional surface known as the separatrix. It is sustained by a current sheet separating the dead magnetosphere with closed field lines from the active zone supported by open field lines. In an ideal plasma picture, this current sheet is infinitely thin, representing a discontinuity in the solution. In reality, due to some microphysics dynamics within the plasma, we expect this current sheet to spread Article number, page 9 of 39 around this separatrix with a depth depending on the plasma condition. Therefore, in a first stage, we assumed an ideally infinitely thin gap and in a second stage we release this assumption by using a thick slot gap.

Sky maps and light curves
Let us start with the sky map in Fig. 1 showing an example of the photon energy distribution received by a distant observer from the slot gaps, in logarithmic scale, the energy E c = ω c being obtained by multiplying the frequency of the observed radiation ω c (including Doppler and aberration effects) by Planck reduced constant . The parameters used assume an inclination angle of χ = 60 • and a Lorentz factor γ = 10 7 for the particles accelerated in these cavities, which are typical values quoted for example in Becker et al. (2009) and in Beskin et al. (1993). For each point on the map marked by the phase and inclination of the line of sight (φ, ζ), in fact an area of size 0.5 • × 0.5 • , only the energy of the most energetic photons at that point are displayed. This map therefore reveals the GeV photon production efficiency for each couple of phase-inclination (φ, ζ) parameters, but without giving any indication of the actual shape of the spectra for each phase and each angle of the line of sight. We will detail these characteristics later.
Knowing the power of curvature radiation P curv in (7)

Spectra
Looking at the spectra offers another perspective of pulsar multi-wavelength emission properties.
So let us begin by studying the spectra of a thin slot gap without taking into account the weight assigned to each photon as a function of the radiated power (7).
The Fig. 9 represents the phase and line of sight averaged spectra from these gaps for different obliquities χ = {30 • , 60 • , 90 • } respectively in blue, green and red solid lines. We find a distribution consistent with what we saw in Fig. 1 for a obliquity χ = 60 • (green curve) with very few photons on the highest energies side, above 50 GeV and a high concentration of photons in the centre of the spectrum, around a few GeV. The range of the spectrum also varies significantly with the inclination of the magnetic dipole.
Nevertheless, from an observational as well as from a more realistic point of view, it is preferable not to average on the line of sight inclination angles since a particular observer only sees the part of the sky map in Fig. 1 corresponding to a fixed angle ζ (i.e. a horizontal line in this plot). An example of the variation of the phase-integrated spectrum (i.e. over the whole period of the pulsar) for χ = 60 • is shown in Fig. 10 for different values of the line of sight inclination ζ = {30 • , 60 • , 90 • } again respectively in blue, green and red solid lines. We notice a great disparity in the shape and the limits of the spectra depending on the obliquity χ although the peak in the spectra remains at Article number, page 12 of 39 approximately the same energy, around 4 GeV. The smallest obliquity χ = 30 • shows the tightest spectrum with an energy interval restricted to the band [3, 7] GeV. A higher obliquity with χ = 60 • moves this lower limit to about 1 GeV while the upper limit increases to 8-9 GeV. In addition, a double-peak spectrum appears, with one peak at low energy around 1.5 GeV and the other remaining at 4 GeV. We show later that taking into account the curvature radiation power (7) eliminates this second low energy peak. An even higher obliquity of χ = 90 • shifts the spectrum in the opposite direction, towards the higher energies with a lower limit of 4 GeV and an upper limit of 30 GeV. As a conclusion, the observer line of sight strongly impacts the received phase-averaged spectrum. The behaviour of these spectra is a direct consequence of the geometry of the field lines, their visibility by the observer, and the aberration and Doppler effects.
The sky map in Fig. 1 and the spectra in Fig. 9 and 10 have been made on the assumption that for each emission point along the last closed magnetic field line a single photon is emitted, i.e. by applying the same procedure as in Giraud & Pétri (2020), the successive points being separated by a certain length ∆s along the field line and arbitrarily fixed by the numerical integration code. A more realistic view must take into account the effectiveness of this curvature radiation, i.e. its power spectrum. Indeed, since we know the radius of curvature of the magnetic field lines at each point, we can deduce the power radiated by a particle accelerated along them using the equation (7).  in energy of photons produced by the curvature radiation. The units remain arbitrary, we will write them down in UA.
Taking the same approach as above, examples of spectra are shown in Fig. 11. As only one value of the Lorentz factor is used, these spectra show rather abrupt decreases and increases in intensity as a function of frequency. As in the spectrum of Fig. 9, a weak flux of photons is produced with an Article number, page 14 of 39 Q. Giraud and J. Pétri: Frequency resolved emission of pulsars   Fig. 9 and Fig. 11. Indeed, the radiated power changes the number of photons produced but not their energy.
As in Fig.9, these spectra are averaged over the phase and line of sight inclination. A real observer will only measure the spectrum associated with a fixed value of the line of sight ζ. Following the same procedure as for Fig.11, another example of a phase-integrated spectrum taking into account the power of curvature radiation is shown in Fig 12 for   It is not enough to simply count the photons emitted in isolation at each point.
Applying these methods for curvature radiation, we are able to compute phase-resolved spectra for a given magnetic axis inclination χ and line of sight ζ depicting a particular observer. Performing this calculation, the phase-integrated spectra have been divided into 10 regular phase inter-  Note that the maximum intensity of the spectrum is not very sensitive to the phase interval, whatever the inclination and phase considered, the maximum of the spectra is between 1 and 10 in our arbitrary units. Only the inclination of ζ = 60 • seems to show a strong amplitude or even a disappearance of the signal for some given phase intervals.

Thick slot gap
Thin slot gaps are idealized regions of infinitely small thickness producing very sharp gamma-ray peaks due to caustics effects (Morini 1983). Until now we only considered emission along the last closed field lines, neglecting a possible thickness of the current sheet. However, because of the plasma and radiation dynamics within these gaps, we expect an emission layer of finite thickness Article number, page 18 of 39 to form. Therefore, in this paragraph we consider slot gaps relying not only on last closed field lines but also on field lines in their vicinity, some crossing the light cylinder, being open lines, and some staying within the light-cylinder, being closed lines. Therefore, in order to tend to a more realistic description of the pulsar high-energy emission from slot gaps, we added in our simulations emission coming from magnetic field lines slightly above or slightly below the last closed field lines.
Sky maps and spectra shown in this paragraph have been realised assuming a thick layer of emission along the last closed magnetic field lines. The power radiated by particles is multiply by a Gaussian weight w g so the emission is maximum along the last closed field lines and decrease when there is an increase in distance between the line where the emission take place and the last closed field line. The colatitude of the last closed field line foot point on the stellar surface being marked by θ pc and the colatitude of any field line by θ, in the coordinate system aligned with the magnetic field axis, the Gaussian weight is written as We take δ = ∆θ/5 with an angular distance between two successive magnetic foot points at the same longitude ϕ equal to ∆θ = π/100.

Sky maps and light curves
The sky map of the highest photon energy is shown in fig. 14  as already stated before, these sky maps and light curves samples serve as building blocks for the elaboration of more realistic and complete radiation patterns in pulsar magnetospheres. The particle distribution function and the full curvature spectrum must be taken into account. This is why we postponed a meaningful comparison with radio and gamma-ray observations for future works.
Article number, page 20 of 39

Spectra
The emission spectra and maps computed in the previous section have been made with an infinitely thin emission area along the last open field lines. Fig. 19 shows the spectrum of the high energy emission for an emission area of a given thickness, using the same method as that used to draw the sky maps in Fig. 11. According to this spectrum, when a thickness transverse to the emission area is included, additional radiation is received including lower and higher photon energies than those of a thin gap. To the spectrum of the infinitely thin zone are thus added other spectra from neighbouring field lines but whose curvature varies slightly, increasing or decreasing according to the nature of the field line: open or closed, which will obviously lead to variations in the radiated power and energy of the radiation. The differences in the shape of the spectra in Figs. 11 and 19 can also be partly explained by the modulation of intensity under the influence of the Gaussian factor which we have introduced by the weight (17) and which has also been used here. Thus for each point of impact on the celestial sphere, the measured intensity is the power radiated by the particle multiplied by this Gaussian function.
The Fig 20 represents the evolution of the high energy spectra over different phase intervals for an inclination χ = 60 • and several fixed values of the observation angle ζ. It is the same as in  We will now determine from which altitude in the magnetosphere of the pulsar this high-energy pulsed emission originates.

Altitude of emission
Since the frequency of the curvature radiation emitted within the slot gaps has been calculated and since the position of the emission points along the last field lines is known, we can deduce from which altitude in the magnetosphere the energy of the radiation originates. [ field lines are more curved at larger distances from the neutron star. It shows that these points are located in the most remote parts of the magnetosphere and that their high energy may be due to the magnetic field being dragged by the rotation of the neutron star which generates a large curvature of the magnetic field lines. The large curvature of the magnetic field lines responsible for the highenergy photons is shown in Fig. 6 where the radius of curvature of the field lines calculated at each emission point is shown (the smaller the radius of curvature, the greater the curvature). However, only the high-energy emission points have a small radius of curvature, if the higher energy emission areas are located in the upper part of the magnetosphere, this may also be due to the Doppler effect caused by the rotation of the pulsar. The Doppler factor η given by equation (5) is indeed larger for emission points located at a high altitude from the neutron star because their linear velocity β is larger. Fig. 7 gives a map of the emission points with the Doppler factor calculated for each of them, this factor is more important near the high energy emission points. High energy emission near the light cylinder therefore requires a favourable combination of several effects: curvature of the magnetic field lines and blue shifting Doppler factor.
Article number, page 23 of 39

Radio emission
As is well known, the radio emission in the polar caps is caused by the acceleration of secondary pairs, that is particles generated by magnetic photon absorption in the strong magnetic field of the pulsar (Erber 1966), producing an avalanche of leptons. The Lorentz factor of these secondary particles is of the order of γ = 10 2 , as stated in Beskin et al. (1993). In what follows, we choose a mono-energetic distribution function with Lorentz factor γ = 30 in order to obtain frequencies calculated with equation (7) in the observational window of radio pulsars.
In order to avoid a sharp on/off boundary within the polar cap, not reflecting the real radio profiles, a weight is assigned to each photon emitted from the polar cap according to the angular distance θ of these points from the magnetic pole, so that the emission profile resembles a Gaussian function. The technique is similar to the introduction of a thickness in the slot gaps described in the previous section. Thus, when plotting the emission maps and the associated light curves, the intensity received over a 0.5 × 0.5 degree zone of the celestial sphere will not increase by one unit for each photon received. Instead it increases by a weighted intensity I taking into account the Gaussian shape for photons. Photons coming from the north polar cap are weighted according to axis. The Gaussian factor seen in equation (18) is used in order to keep pulses of Gaussian shapes and maximum emission at the centre of the polar cap.
We generated a random but homogeneous distribution of 10 8 points inside both polar caps as depicted in Fig. 25. To obtain such a distribution, the colatitude θ of the emission points in the spherical coordinate system oriented along the magnetic axis (where the colatitude θ mp of the north and south pole are respectively 0 • and 180 • ) is determined through the formula with θ pc the latitude of the point where the last closed field line crosses the stellar surface and X is a random number in [0, 1]. The longitude ϕ of the emission point is determined by introducing another random number Y in [0, 1] such that where ∆ϕ = 2 π/N ϕ is the polar angle between two neighbouring points of a polar cap rim and N ϕ the number of point samples in the longitudinal direction.

Sky maps
Let us now determine the radio emission maps as a function of frequency in order to study the evolution of the pulse shape as measured by a distant observer. of sky maps and a corresponding sample of light curves obtained for an emission volume in the range r ∈ [2, 3] R above the polar caps for different frequency intervals. Note that the frequency interval breakdown is not the same for the different obliquities χ because the spectra do not have the same range or the same upper and lower limits. The choice of frequencies has been optimized to make the intensity maps stand out best.
These maps reveal a narrowing of the pulse profile towards the high frequencies since these photons come from the deepest regions of the magnetosphere. There is also an east-west asymmetry in these profiles between the rising and falling ramps. Sometimes two pulses are observed in the profile, sometimes only one. It should be stressed that these results are only preliminary and that any comparison with observations will have to take into account a particle energy distribution in the form of a power law and not simply a mono-energetic distribution as done in the present work, see the particle distribution function discussion in section 2. The characteristic curvature frequency must also be replaced by the continuous curvature spectrum centred on this typical frequency. This will make the frequency variation of the pulse profile smoother and the transition between the frequency bands more continuous. After this further step, however left as a future work, detailed comparison with existing multi-wavelength observations can be performed.
According to these maps, especially those in Fig. 28, low frequencies are emitted rather at the edge of the area while the higher frequency components are emitted towards the centre of the area.
This may seem counter-intuitive since for a static dipole, the centre of the pulse, corresponding to the magnetic axis, is a straight line of zero curvature and therefore of zero frequency and power. It increases when the magnetic axis joins the rotation axis, i.e. when χ decreases. For example, for χ = 30 • this minimum radius is 0.05 R cyl and moves outside the centre of the cap. At the same time the maximum curvature radius increases significantly when the rotator is close to the aligned configuration. It is only 0.048 R cyl for χ = 90 • and reaches up to 0.35 R cyl for χ = 30 • which is an increase of almost one order of magnitude. The photons produced will be much less energetic in this latter case. Within the limit of a perfectly aligned rotator χ = 0 • , the radius of curvature becomes infinite at the centre of the cap, on the magnetic axis, and the curvature radiation is extinguished at its centre.
The altitude at which photons are produced and their detachment from the magnetosphere in the direction of the observer is not precisely constrained by the observations. The radio polarisation Article number, page 28 of 39 Fig. 22. Sky maps of photon energy for an inclination χ = 30 • , each map represents the contribution from a piece of the magnetosphere : those located in a spherical shell of thickness R whose altitude is increased by R to the next map.
data provide some indication for young pulsars, but the error bar remains appreciable. There is therefore some freedom in the choice of the height of the emission sites above the polar ice caps.
For this reason we have also drawn emission maps and light curves for an altitude in the range [3,4] R as shown in Figs. 32, 33 and 34 for the same frequency intervals above the polar caps.
Pushing the emission altitude further away from the star causes a spreading of the radio pulse profile since the field lines are divergent. This is reflected in the fact that high altitudes produce the lowest frequencies in relation to wider pulses while low altitudes produce the highest frequencies in relation to narrower pulses, in accordance with the radius-to-frequency mapping model. But in our case, we dispense with the static dipole to take into account all the effects due to the rotation of the magnetic field, the effects of signal aberration and delay.
In order to facilitate the visualisation of the evolution of the radio pulses as a function of frequency, the shape of a light curve has been plotted in Fig. 35  curved compared to the low-frequency photons produced at higher altitudes where the field lines are more curved, due to the sweep back effect, eq.(16). This effect depends on the ratio (R * /R cyl ), it is negligible for young pulsars but its observational signature should be perceptible for the fastest millisecond pulsars.
The phase shift can be estimated from eq.(13), (15) and (16) in paragraph 2.6. The maximum offset thus calculated (for a photon emitted at 1 R above the surface and another emitted at 2 R and always with χ = 60 • ) is of the order of 4% of the phase. These effects are therefore not sufficient to explain the shift of the radio emission peak in Fig. 35 but they can certainly play a role in the radio emission observed in these emission areas.

Spectra
As for the high energy radiation, fig. 36 shows the radio spectrum for different values of the obliquity χ of the magnetic field. The points of emission for this spectrum are located on the stellar surface and with the same spatial distribution as shown in fig. 25. in the opposite direction of the rotation as shown by Shitov (1983), delaying the radio signal by an amount ∆t B given by (16). The Doppler effect also contributes partly since the corotation speed v cp = r Ω sin χ is more significant for the perpendicular rotator because of the longer lever arm of length r sin χ . The amplitude of this Doppler effect is related to β cp = (r/R cyl ) sin χ . Generally speaking, the characteristic radio frequency ν radio depends on the ratio (r/R cyl ) and on sin χ . Therefore the range of radio waves around ν radio behaves as an increasing function of sin χ . On the other hand, this shift of the radio spectrum towards higher frequencies is attenuated for rotators slower than those studied within this paper. Indeed, for young pulsars, of period P greater than about P > 100 ms, the ratio between the neutron star radius and the light cylinder radius is much smaller than one, (r/R cyl ) 1, hence producing a less perceptible variation in frequency as a function of sin χ .
In the above proposed image, the photons escape directly from the star's surface. In our simulations, the period of rotation of the pulsar is limited to 2 ms because of the prohibitive computing time for a star of larger period. The neutron star surface is fixed at a radius 0.1 R cyl . Nevertheless, for slow pulsars, this radius of 0.1 R cyl corresponds to the average real altitude of radio emission (it is in fact a little lower, of the order of 0.05 R cyl , see for instance Mitra (2017)). The width of the pulses as well as the delay between the radio peak and the first gamma-ray peak thus remains realistic in spite of the too high period. However, the assumption that radio emission is produced at a fixed altitude is not entirely correct. Indeed, it is known that the highest frequencies are generated at lower altitudes because of the radio-to-frequency mapping. It will therefore be necessary to include an additional degree of freedom for the true position of the emission sites by introducing, for example, an interval of variable altitudes as was done in the high energy band.
The figures in 37 show the radio spectrum for different values of the obliquity χ , when the emission area is located above the polar caps, respectively in the range r ∈ [1, 2] R on the left part and in the range r ∈ [3, 4] R on the right part. Compared to Fig. 36, the emission no longer occurs in a given radius r but for an entire range of radii r ∈ [r 1 , r 2 ]. The emission is no longer integrated on a surface (the polar cap) but in a whole three-dimensional volume. The result is a wider range of values of radii of curvature and consequently a wider radio emission spectrum, regardless of the value of χ . This spread is maximum for χ = 30 • and less for other values of the inclination of the magnetic axis such as 60 • and 90 • . Shifting the emission volume to a higher altitude produces an additional spread of the spectrum, which can be seen by comparing the left and right plots in fig. 37. We conclude from these spectra that the farther the radio emission zone from the surface of the neutron star, the wider the frequency range of the received radio radiation. This phenomenon is very important for the smallest inclination investigated, namely for χ = 30 • .
This broadening of the spectra is due to the fact that our emission zones above the polar caps possess a certain thickness and is therefore the result of the variation in the curvature of magnetic field lines. However, other phenomena may affect the shape of these spectra, such as the geometry of the field lines, which is certainly varying with altitude, or also the Doppler factor, which depends on the instantaneous rotation speed v = Ω ∧ r = r Ω sin θ e ϕ and increases with altitude due to the corotation of the magnetosphere with the neutron star.
Article number, page 32 of 39

Conclusion
We studied one of the main emission mechanisms taking place within the pulsar magnetosphere, namely the curvature radiation from the slot gap and polar cap inside the light cylinder. We showed detailed radio and high-energy light-curves and spectra in the Fermi/LAT band. Using realistic values of the Lorentz factors for primary and secondary beams, our spectra fall into the right windows. We computed phase-resolved spectra and found that the spectral energy distribution depends on the phase interval. Also, the averaged radio spectra significantly depend on the pulsar obliquity whereas the high-energy part seemed much less sensitive to the geometry.  In order to get a more complete view of this emission, it might be useful to add synchrotron radiation and inverse Compton scattering to future simulations, as well as other models with different emission zones, such as outer gaps or the striped wind. The calculation of the spectra will also gain in realism if the mono-energetic particle distribution is replaced by a power law distribution and the δ approximation of the curvature spectrum by its spectral energy distribution.
The polarisation of the pulsed emission could also be added to our simulations of radio and high energy emission. In the radio domain, this polarisation already makes it possible to shift the emission sites at high altitude for slow pulsars with periods above 100 ms. The imminent observation of the polarisation in X-rays will bring a new strong constraint on the high energy part of the spectrum, determining its location either inside the light cylinder or outside within the wind.
Finally, a discussion of these results in relation to pulsar multi-wavelength observations should constrain many dynamical parameters of the relativistic plasma by comparing, for example, the light curves as well as the radio and gamma-ray spectra for pulsars detected simultaneously in these two energy bands.  were funded by the Equipex Equip@Meso project (Programme Investissements d'Avenir) and the CPER Alsacalcul/Big Data.