Analytical expressions for pulse profile of neutron stars in plasma environments

We present an analytical study of light curves of slowly rotating radio pulsars with emphasis on the chromatic effects derived from the presence of a plasma environment; analyzing the effects of the compactness, the metric model, and the electronic plasma density profile. After doing a numerical integration of the trajectories and luminosity curves of pulsars for different spherically symmetric metrics representing the exterior region of the pulsar, we generalize the approximate Beloborodov formula in order to include plasma corrections, obtaining simple analytical expressions for the trajectories and the observed flux and significantly simplifying the calculation of the pulse profiles by a drastic reduction of their computational cost. We study the errors committed by our approximation, comparing the numerical and analytical procedures. We also show how to use the new formalism to model the flux coming from different emission caps, not necessarily circular or antipodal and including the case of ring-shaped hot spots. Finally, we extend the classification introduced by Beloborodov to the case of two distinguishable, non-antipodal, finite size emission caps, showing the respective classification maps and some of the characteristic pulse profiles.


I. INTRODUCTION
Neutron stars represent one of the most relativistic stellar objects accessible to observation [1,2].Due to the extreme conditions to which they are subjected, most of their properties remain uncertain, making the determination of their equation of state an open problem in nuclear astrophysics [3][4][5].However, this difficulty can be overcome by observing the neutron stars themselves and the phenomena associated with them, such as their luminosity or their lensing power.At the same time, neutron stars are one of the best natural candidates to probe physics in extreme environments and test the theory of General Relativity in the strong field regime [6].Indeed, the discoveries of neutron stars of around two solar masses [7], as well as the observation of gravitational waves resulting from mergers of binary systems [8], impose severe constraints on the equation of state of these objects, in the high-density regime.More recently, the mass and radius of neutron stars have been estimated by fitting the rotating hot spot patterns using data of Neutron Star Interior Composition Explorer (NICER) and X-ray Multi-Mirror (XMM-Newton) X-ray observations [9][10][11].
Most of the neutron stars emit highly collimated jets of electromagnetic radiation [12] coming from highly localized "hot spots" that emit significant quantities of observable radiation.Neutron stars exhibiting these "emission caps" are called pulsars, and because of its rotation, they produce observable luminosity curves, similar to a lighthouse.
These pulsar pulse profiles encode information about the physical properties of the emission caps, such as their size, geometry, temperature distribution and spectra, as well as containing information about the geometry of the spacetime itself around the star [13][14][15][16][17][18][19][20][21].
The observation of these periodic electromagnetic pulses offers another way to investigate the properties of these extremely compact objects.Due to their strong gravitational field, the radiation in the vicinity of neutron stars experiments a strong deflection, generating a gravitational lensing effect, whose magnitude is intimately related to the mass and radius of the star, so that their observation allows in principle to infer these parameters, or at least to constrain their values, providing information about the stellar structure that facilitates the determination of the equation of state [22,23].
The problem of modeling pulse profiles has been extensively studied in the literature.Studies have been done considering the Doppler effect, aberration; time delay [24] and fast rotation or stellar oblongation [25][26][27][28].In particular, In [13] it was studied gravitational lensing properties in slowly rotating neutron stars described by the Schwarzschild metric, consisting of one or two antipodal, uniform, circular polar caps on the surface that emit photons to the observer at infinity.In general, the treatment of photon propagation in a spherically symmetric and stationary spacetime leads to elliptic-type integrals or even more complicated expressions that need to be solved by numerical evaluation.The implementation of numerical methods is often unavoidable when considering a continuous surface temperature distribution, anisotropic emission or an arbitrary shape of the emission regions.While these models may be able to describe the pulse profile of stars of any compactness, it requires a complicated and computationally expensive numerical treatment.In this context, Beloborodov found in [29] a simple analytical expression for the deflection angle in a Schwarzschild metric background, that replaces the elliptic integral, approximating with high accuracy the pulse profiles of point emission caps, assuming that the compactness of the star is not too big.In turn, Turolla and Nobili [30] generalized this analysis to extended, uniform, and circular caps, giving simple analytical expressions that closely approximates the total observed flux.The Beloborodov formula not only has been used to study pulse profiles of NS under different situations [31][32][33][34][35][36][37] but also to study polarimetric images of accretion disks around black holes [38][39][40].
It is to be expected that neutron stars are immersed in a dense, plasma-rich magnetospheres [41,42].In the visible spectrum, the modification in photon trajectories and pulse profiles due to the presence of plasma are negligible, so generally gravitational lensing theory deals with the propagation of light rays in vacuum, where the trajectories and deflection angles are independent of the photon frequency, and thus the effects are achromatic.However, the effects of the optical medium cannot be safely neglected in the radio frequency range, where the refraction index of the plasma causes strong frequency-dependent modifications of the usual gravitational lensing behavior.In this range, photon propagation becomes frequency dependent, giving rise to chromatic phenomena.These effects have been studied in the context of black holes and its shadows [43][44][45][46][47][48], in weak fields under the weak and strong lensing regime [49][50][51][52][53][54][55][56][57][58], in studies of microlensing [59], and on its effects on the fast radio burst [60,61] to cite some examples.The influence of plasma in light propagation and associated light curves produced by radio pulsars have been numerically studied in the past by [62][63][64][65] assuming a Schwarzschild model for the exterior of the metric and a cold non-magnetized plasma (see also [66] where a Minkowski metric is assumed, but the magnetic field influence is taken into account).From the observational point of view, in [67] the pulsar emission regions were resolved and amplified by plasma lensing effect in an eclipse binary.Similar results have been carried out by other pulsars [68,69].
In this work, using numerical evaluations and introducing an analytical approach that generalize the Beloborodov's formula, we study the gravitational lensing phenomena in slowly rotating radio pulsars taking into account the dispersive effect of the pulsar's magnetosphere on the propagation of light rays.We will make use of alternative metrics to Schwarzschild to describe the exterior region in order to analyze the difference in their influence on the luminosity curves.More precisely, we will focus on the Schwarzschild (Sch) and the so called Reissner-Nordström like metrics consisting in solutions from alternatives gravitational theories which are functionally similar to the RN metric but with a parameter Q unrelated to the electrical charge; for example associated to tidal charges in brane theories [70], associated to parameters of the theory, as in Horndeski theories [71] or being proportional to the gravitational mass of the star in Modified Gravity [72].The particular choice of these metrics is only to demonstrate the differences that appear in the properties of compact objects in the strong field regime and that can potentially be observed.We will consider circular, homogeneous and isotropic emission caps (either thermal emission hot spots or electro-magnetic emission poles) of finite size, in highly compact slowly rotating radio pulsars, so we will not take into account Doppler effect, aberration, time delay, fast rotation or stellar oblongation.As in [62][63][64][65] the neutron star will be embedded in a pressureless and non-magnetized plasma environment, whose electron density will be given by a power law inversely proportional to the radial coordinate r, resulting in a chromatic analysis of the problem 1 .This paper is organized as follows.In Sec.(II) we introduce the family of metrics that we will consider in the rest of the work.In Sec.(III) we study numerically the photon trajectories considering the different metric models and plasma density and distribution, emphasizing the chromatic effects derived from the presence of the plasma environment, whose interaction with the light will depend on its wavelength.In Sec.(IV) we discuss the luminosity curve resulting from a single circular, homogeneous and isotropic emission cap on the stellar surface and how this is affected by the spacetime model, the distribution and density of the plasma, the angles between the cap and the pulsar rotation axis and between the rotation axis and the observer or the compactness and charge of the star.This analysis will be extended to the case of two identical and antipodal caps.In Sec.(V), motivated by Beloborodov's approximate analytical model, we modify the well known analytical formula in order to introduce corrections that take into account plasma environments, including now analytical expressions to describe uniform and circular extended caps immersed in a plasma environment.Then we apply the new formalism to the different metric models used above, obtaining for each of them simple analytical expressions for the photon trajectories (V A) and the observed flux (V B) that considerably simplify the calculation of the pulse profiles, drastically reducing their computational cost.We will compare the numerical results with the analytical approximations, studying the relative errors committed by the approach and its dependence on different magnitudes of interest according to the characteristics of the problem.Once the validity ranges of our model have been established, in Sec.(VI) we will show how to use the resulting approximations to easily model more realistic, non-antipodal, homogeneous or circular emission caps.In particular, we will focus on the case of ring-shaped caps.Under this new formalism it will be possible to examine in detail the effects of the shape and size of individual hot spots, without any great loss of accuracy or generality, thus facilitating the study and understanding of these objects.In Sec.(VII) we expand the classification system intro-duced by Beloborodov for the case of two non-antipodal, distinguishable caps, which can be also applied even in the case that the influence of plasma on the propagation of light rays can be neglected.Finally, in Appendix (A) (see supplementary material) we include some plots obtained for different plasma distributions in order to show how sensible are the photon trajectories and pulse profiles to the behavior of the electron density profiles, expanding the results of Secs.(III) to (V).

II. METRIC
Let us use normalized fundamental units, c = 1, G = 1 and = 1, being c the speed of light in vacuum, G the universal gravitational constant and Planck's constant, and the metric signature (−, +, +, +).For the exterior region of the star, we employ simple metric models, consisting of asymptotically flat and static spacetimes with spherical symmetry.These assumptions allow us to express the metric as follows 2 where the functions A(r), B(r) and C(r) depend only on the radial coordinate r, satisfying also the asymptotically flat conditions, i.e., ( The simplest metric in GR describing a static, spherically symmetric and asymptotically flat spacetime with a gravitating M mass is the Schwarzschild (Sch) metric.This solution is a useful approximation for describing slowly rotating astronomical objects.
A generalization of the Schwarzschild metric, that describes a spherically symmetric metric of a massive body with non-zero net electrical charge Q is the Reissner-Nordström (RN) metric, which is a solution of the Einstein-Maxwell equations [73,74].
As mentioned in the introduction there exist solutions of alternative gravity's theories which are functionally similar to the RN but with a charge parameter q * unrelated to the electrical charge (See Table I).In particular, in Brane world (BW) theories, a negative q * is theoretically preferred [75].Studies of NS in braneworld theory can be found in [76] (where it is shown that a uniform 2 Even when it is always possible to choose a suitable coordinate system for which the spacetime is characterized by only two metric functions A and B, in what follows we will prefer the form (1) because it allows to write general expressions that remain valid for a vast family of coordinate systems.density star can be matched to a RN-like solution for the exterior of the star, see also the discussion in [77]) and observational constraints on the value of q * in NS are presented.For a study of the plasma magnetosphere around NS in BW theories we refer to [78].
Given the similarities between the models for alternative theories when q * = 0, when q * < 0 we will refer to them as RN-like metrics, associating the parameter Q 2 = −q * , while when q * > 0 we will refer to them as RN metrics, associating the parameter Q 2 = q * , so that both models turn out structurally identical and can be described by a single parameter q * .Even though in the metrics used in this work C(r) = r 2 , in the following, we prefer to keep the coefficient C(r) in order to obtain general expressions that hold valid for alternative metrics to those used here.

III. RAY TRACING
The stellar surface will be placed at the radial coordinate r = R, the region r < R will be totally opaque since it corresponds to the stellar interior.The light emission will originate from this surface, with an intensity I dependent on the angle δ between the normal direction to the stellar surface and the photon emission direction (see Fig. (1)).The observer will be considered static at r = r O with r O → ∞.We will employ Greek indices to denote spacetime coordinates, while we will reserve Latin indices for spatial coordinates only.
In this section we will follow the work presented in [62], recovering his results for the Schwarzschild metric and extending its application to the different spacetimes mentioned in Table (I).Consider a spacetime of the form given by Eq. ( 1) where the metric functions A, B and C depend parametrically on the mass M and charge q * of the pulsar.In the geometrical optic limit, photon trajectories within a plasma medium are described by the Hamiltonian [79], being the refractive index of a general dispersive medium not necessarily coming from a plasma, while p α is the linear photon 4−momentum, and V α the plasma 4-velocity.
The light rays are obtained as solutions to the equations of motion We will assume a spherically symmetric plasma distribution around the pulsar.When light passes through it, the plasma will act as a dispersive medium, with a refractive index n = n(x α , ω) which will depend on the frequency ω = −p α V α of the photon.Assuming that the plasma is a static medium such that V t = −g tt and V i = 0, it follows that so that the Hamiltonian does not depend on t or φ, so p t and p φ are conserved.We set the timelike component of the momentum as the asymptotic frequency of the photon as measured by an asymptotic observer at rest From Eqs. ( 6) and ( 7) we can find the relation for the effective redshift, where ω(r) is the frequency of light measured by a static observer at the radial coordinate r.
In the rest of this work, we will assume a nonmagnetized pressureless plasma with the index of refraction taking the form where being ω e (r) the plasma frequency [80], e and m e the electron charge and mass respectively, 0 the vacuum permittivity and N (r) the number electron density of the plasma 3 .Together with Eqs. ( 3) and (9), the above expression allows us to obtain the dispersion relation for photons, Given the spherical symmetry of the problem, we can always choose a referential such that p φ = 0, corresponding to planar meridian trajectories (φ=constant) 4 .For these trajectories, it follows that p θ is also a constant of motion.
To construct the ray tracing, i.e., the spatial trajectory of photons, we can get rid of the parameter λ of Eq. ( 5) as follows, At the same time, from Eqs. (11) and ( 5) it can be deduced that Here the positive solution corresponds to the case where both r and θ are increasing, while the negative solution otherwise.We now introduce the impact parameter b and the asymptotic group velocity n 0 to express the constant p θ , which can be written as The group velocity usually tends asymptotically to 1, except in cases where one has a constant, or non-zero, plasma distribution around the observer, in which case we have with r O is the radial coordinate of the observer.Thus, the trajectory equation becomes [51] This equation allows us to describe the photon trajectory in a meridian plane in the general case of strong deflection, which is appropriate for orbits near the surface of the star.The refractive index of the plasma in this expression introduces the frequency dependence of the orbits.In order to facilitate the numerical integration, we introduce the variable u = 1/r, for which it turns out that where both the metric elements and the refractive index n, must be written in terms of u.We now would like to FIG. 1: A photon's trajectory, with impact parameter b, from the stellar surface at R to the observer at rO.The deflection angle β is calculated as θ − δ.
make use of Eq. ( 17) to study the outgoing ray tracing of a neutron star with a cold, nonmagnetic plasma atmosphere modeled by (9), where the photons move with φ = const.
For this, we will place the observer on the θ = 0 axis of the coordinate system, at a large distance from the star, at r = r O → ∞, such that the line of sight connecting the center of the star with the observer corresponds to the impact parameter b = 0. We are interested in describing rays that leave the surface of the star at u R = 1/R, θ R = θ(R) and reach the observer at u O = 0, θ O = 0.The angle θ(R) of the photon at the instant it leaves the pulsar can be calculated as Now, if the neutron star is sufficiently compact, it is possible that the value of θ given by Eq. ( 18) may be larger than π, so this angle will not be exactly the photon's colatitude coordinate at the moment of leaving the stellar surface, but rather the angle from the observer's position vector to the photon's position vector measured counterclockwise, as shown in Figure 1.From now on, when we refer to the θ angle of the photon we will mean the angle we have just defined, while we will use the symbol Θ to talk about its colatitude.
At this point it only remains for us to find an expression for the impact parameter of the rays whose trajectory connects the surface of the star with the observer.For this, we define δ as the angle between the radial and angular components of the photon momentum.At the surface of the star r = R, δ is the angle at which the ray is emitted with respect to the normal.
Under these assumptions, it is possible to express the impact parameter b in terms of δ as follows where we see that the maximum impact parameter b max occurs for δ = π/2.Thus, photons coming from the stellar surface will arrive at the observer with an impact parameter less than or equal to b max .The presence of the refractive index in Eq. (19) predicts that the apparent size of the pulsar for a distant observer vanishes when n(R) → 0.
From these expressions we can obtain the maximum angle θ max (b max ) of the stellar surface that is visible to the observer.Thus, the observer has access to all points on the surface with Θ ≤ θ max (b max ), corresponding to b ≤ b max , while the points with Θ ≥ θ max (b max ) will form an invisible region.As in the case of light propagation without plasma, the resulting fraction δf of the visible stellar surface is given by δf = (1 − cos[θ max (b max )])/2 (cf.[17], Eq. 26).
The expressions shown up to this point are valid for a generic plasma distribution, with a refractive index n = n(x α , ω) that depends on both photon frequency and position.For the case where ω e ω, n → 1, thus recovering the case of pure gravity.The propagation of a photon through this medium requires that ω e (r) < ω(r) at every point along the trajectory.The necessary condition for the rays to propagate from the star's surface to infinity then results for all r ≥ R.
Numerical ray tracing requires a specific plasma distribution.Following [62], we choose some number density profiles that contain as particular case the Goldreich-Julian density [41].Normally, the pulsar magnetosphere contains plasma with densities that are sustained by Goldreich-Julian currents, transferring charge carriers from the neutron star surface to the magnetosphere.Therefore, these currents provide a lower limit for the concentration of charge in a standard magnetosphere.Using the plasma frequency presented in Eq. (10), let us consider a radial power law with h ≥ 0. In order to simplify the numerical integration and to get rid of the largest number of trivial constants, we now introduce the parameter , defined as the ratio between the plasma frequency and the photon frequency over the stellar surface, so that the index of refraction can be simply rewritten as while the propagation condition is now expressed as for all r ≥ R. For the metric models and plasma distributions that we will use in this paper, the above condition can be reduced to in the cases where the plasma density decays to zero faster than A(r) grows as r goes from R to ∞, or to when we consider a constant plasma density that does not vanish around the observer.
Coming back to Eq. ( 19), if we assume a vanishing electronic density at the neighborhood of the asymptotic observer, it follows that independently of the plasma profile model under consideration, the maximum impact parameter b max is expressed in terms of as Hence, varying the observation frequency (and therefore the value of for a given system), we obtain information of the value of the metric quotient A(R)/C(R) at the surface of the star.Of course, b max is not a direct observable, however, due to the relation between b max and θ max we can infer that observation of pulse profiles of a given NS at different radiofrequencies codifies information on the geometric properties of the spacetime at the neighborhood of the star.

A. Summary
We already have all the necessary tools to perform ray tracing from a highly compact and relativistic object immersed in a pressureless non-magnetized plasma medium to an observer located at infinity.Below, for the reader's convenience, we list the steps to be followed to carry out the ray tracing.
1. Select a metric.Some examples are shown in the Table I.
3. Calculate the asymptotic group velocity n 0 (at infinity) from Eq. ( 15) and check that satisfies the pertinent propagation condition.
4. Obtain the maximum impact parameter b max from Eq. ( 19) and choose an impact parameter b in the range 0 < b < b max .
5. Calculate the angle θ(b) from which the ray will leave the stellar surface according to Eq. ( 18).

B. Results
The ray orbits obtained for the Schwarzschild, RN and RN like metrics are shown in Figures 2 and 3 below.A highly relativistic star with radius R/r S = 1.60, with r S = 2M the Schwarzschild radius, was considered.The traced rays count with positive impact parameters in the range from 0 to b max , whose values correspond to integer multiples of b max /5.The absolute value of the charge in each figure is We chose this value since it allows us to distinguish the different trajectories obtained for each metric model.In turn, each figure corresponds to a particular plasma distribution, and compares the results obtained for different ratios of frequencies.
For this, we wrote a code in Fortran 90 that uses the Runge-Kutta 4th order method to solve the differential Eq. ( 17) and a Gaussian quadrature rule to integrate Eq. (18).
We will start considering the plasma profile suggested in [62], where h = 3 was used as it considers the Goldreich-Julian density and its dependence on the polar magnetic field strength.Thus, the index of refraction is expressed in the whole space as The results obtained for this distribution of plasma are shown in Figure 2.Moreover, we shall consider a homogeneous plasma distribution with constant density, so that the refractive index is expressed as As it is well known, there is a correspondence between the dynamics of light rays of frequency ω in a homogeneous, non-magnetized, pressureless plasma, with frequency ω e , and the timelike geodesic motion of test massive particles of mass µ and energy E ∞ measured by an asymptotic observer in the same gravitational field [81].In particular, given the transformation ω e → µ and ω ∞ → E ∞ it is possible to use the same Hamiltonian to describe both phenomena.Thus, a constant plasma density is a model of great interest, where photons behave as if they have an effective inertial mass.In fact, it has been suggested in the literature [82,83] that neutrinos from magnetosphere caps in the PeV regime could be detected by various neutrino detectors such as the IceCube.Their dynamics can then be equivalently studied from the study of photons in a homogeneous plasma.The plots obtained for this distribution are shown in Fig. 3.In addition, in order to see how sensible are the photon orbits to the behavior of the electron density profiles, two others plasma densities are shown in the Appendix A, one of them of the form of Eq.( 21), with h = 2 and another one with an exponential decay profile given by n 2 e (r) = 1 The dotted line corresponds to Schwarzschild without plasma, green to Schwarzschild, blue to q * > 0 and red to q * < 0. As the ratio = ωe(R)/ω(R) increases, both the apparent size of the star bmax and the maximum angle of visibility θmax decrease.We see that bmax and θmax increases for q * < 0 and decreases for q * > 0, with Schwarzschild metric in between.
In general, as in [62], we observe that plasma concentrations produce a divergent lensing effect that competes with the convergent effect of gravitational lensing.Thus, while the gravitational potential tries to deflect the light rays towards the star, increasing the values of θ max and b max , the diffraction produced by the plasma deflects the rays in the opposite direction, moving them away from the star and reducing the aforementioned values.This produces the photon trajectories to have a S shape, where the convergent effect of the gravitational lens dominates in the vicinity of the surface, being overcome by the divergence of the plasma lens as we move away from it.These effects become more pronounced as the plasma density increases (or alternatively the observation frequency decreases), i.e., for → 1 the trajectories converge more and more at y = 0 axis in Fig. (2), that is, θ max → 0 and b max → 0, so that the star is no longer visible.
From these plots, and those shown in Appendix (A) we can see that the divergent lensing effects produced by the plasma distribution are magnified as the plasma concentration decay rate increases (i.e. higher h or exponential decay).Thus, the plasma profile with a decay given by a FIG.3: Same as in Fig. 2 but for h = 0.As the ratio = ωe(R)/ω(R) increases, both the apparent size of the star bmax and the maximum visibility angle θmax diverge.Note there is not plot at = 0.9, because it does not satisfy the propagation condition (26).In addition, for q * < 0, as the value of increases so does the radius of the photon sphere, eventually exceeding the stellar radius, so that trajectories with b ≤ bmax do not reach the observer.That is why we do not show the q * < 0 example in the bottom panel.
power law exhibits marked divergent lensing properties, while a homogeneous plasma distribution does not produce any divergent effect.This is in accordance with the basic notions of optics in refractive media, which indicate that the greater the variation of the refractive index at an interface, the greater its deflection power.
A curious result happens when considering a homogeneous plasma distribution.As we already said, the divergent lensing effect of the plasma is completely cancelled out in this profile.In Fig. 3 it is clear that both the impact parameter and the visible surface of the star increase at larger values of .If we analyze Eqs. ( 15) and ( 19), we see that this is because, given a constant plasma density, n 0 → 0 as → A(R) < 1, and the impact parameter b diverges as 1/n 0 .However, as we will shown in Eq. (39), the observed light intensity goes to zero as n 2 0 .We add that, the quotient = 0.90 does not satisfy the propagation conditions for this profile, so it was not plotted.
We note that for h > 0 increasing the charge q * of the star tends to decrease both the observable surface area by reducing the value of θ max , and the apparent size of the star by reducing the value of b max .The opposite effect results for decreasing q * .For the particular case of the RN metric (q * > 0) similar conclusions were observed by [18] and [15].

IV. PULSE PROFILE
Consider a radio pulsar with bright polar regions of emission on its surface that will produce periodic pulses in the luminosity curve.The observed flux F will be given by the integral of the intensity times the solid angle measured by the observer on his celestial sphere.In the case that the same emission region produces multiple images, all of them must be considered when integrating the flux.While this approach is slightly idealized, it allows a first approximation to describe qualitatively the effects associated with the various parameters of the star, its emission caps, the plasma environment and the spacetime itself.
The geometry of this model is detailed in [13] for the gravity pure situation.Again, in this section we will follow the work [62], reproducing its results for the Schwarzschild metric and expanding the procedure to other spacetimes.
The surface of the star, at r = R, is described by the angular spherical coordinates Θ and φ that we used above.In general, we will use θ instead of Θ.We will place the observer on the positive semi-axis z (Θ = 0), at r = r O → ∞, defining its coordinate system with respect to a plane perpendicular to the line of sight (the "detector").The detector surface registers the flux received from the object (the image) over the angles of its celestial sphere θ and φ .Since the star is spherically symmetric, the trajectories will remain in the same meridian plane, so we can assume φ = φ .
Hence, we can express the solid angle element in the observer's coordinate system as follows dΩ = sin θ dθ dφ, where the angle θ is assumed to be sufficiently small, so This allows us to express the latitudinal angle θ in terms of the impact parameter b, where b ≈ r O θ , resulting in where b is given by Eq. ( 19) and varies between 0 and b max for trajectories starting from the stellar surface.Thus, employing the approximation b ≈ r O θ and Eq. ( 18) we can map each observed image point in (θ (b), φ ) to its location on the surface of the star in (θ(b), φ), recalling that φ = φ .That is, we can map the observer's two-dimensional solid angle element onto the corresponding surface area element on the star.To find the pulse profile in question, we must find the flux of the emission region on the stellar surface and calculate the projection of this area on the observer's sky.Let us consider for now a single polar emission cap, also called hot spot, of circular shape and angular halfaperture θ c , whose position with respect to the center of the star is indicated by the versor Ĉ, being centered at θ = θ 0 .We further define the versor along the line of sight L in the outgoing direction from the pulsar to the observer, resulting in θ 0 = cos −1 Ĉ • L the angle between both versors.In general, θ 0 will be a time-dependent quantity since the star rotates about an R axis that need not coincide with the line of sight L. The angle between the rotation axis R and the line of sight is defined as ξ = cos −1 R • L, while the angle between the rotation axis and the center of the cap is given by χ = cos −1 R • Ĉ.Both ξ and χ are constants.These angles are represented in Fig. 4.
To find the angular position of the cap θ 0 as a function of time t, we must express the phase of the pulsar γ p (t) = Ωt, where Ω = 2π/P is the angular velocity, with P being the rotation period of the star.Thus, the latitudinal orientation of the cap can be expressed as where it can be seen that γ p (t) = 2nπ corresponds to the maximum value of θ 0 (t).For a single cap, is convention- In a spherical coordinate system fixed to the cap, whose pole is given by the versor Ĉ, an arbitrary point on the edge of the cap is expressed as (θ c , Φ), while the same point is expressed as (θ, φ b ) in the object coordinates, marking φ b the boundary of the cap in the parallel θ.Using the transformation between the Cartesian components of the two systems, we arrive at the expression where 2φ b (θ) should be understood as the angular length of the circle segment in the observer's sky with the where ∆ is the argument of the function cos −1 in Eq. (34).
Actually, the definition of φ b is more complicated than the one described above, since it is not well defined for any θ 0 orientation.This point is discussed in more detail in [13] and [18].However, for the purpose of this paper, it is sufficient for us with Eq. ( 35), which will allow us to simplify the solid angle element.
Following to [62], we introduce the normalized impact parameter x = b/M (M = 1).The solid angle segment results so that, to determine the observable flux of an emission cap, only a one-dimensional integral over the impact parameter x is required.In order to contemplate the occurrence of multiple images, we must allow θ to take all values in the interval [0, θ max ].For convenience, we rewrite Eq. ( 18) in terms of x as The emitted intensity I em of a bright spot on the stellar surface at r = R in a dispersive medium is related to the observed intensity I obs of that spot by a detector at r = r O by [84] The observed intensity is then We will now consider the possibility that the emitted intensity of the polar cap depends on the emission angle δ with respect to the stellar surface normal.Let us take for convenience where I 0 is a parametric constant and f B describes the emission anisotropy.The factor r 2 O has been included so that it can be cancelled out when joining this expression with Eq. ( 36), and has been divided by C(R) to normalize the flux by the emission area of the star.In this way, the expression for the observed flux can be rewritten as On the other hand, the emission angle at the surface is a function of x since, as can be deduced from Eq. ( 19), it turns out to be so that we can express the observed flux differential dF in the detector's solid angle differential dΩ as (43) To obtain the total observed flux F , it only remains to integrate the previous expression in x, obtaining This result, derived in [62] for the particular case of a Schwarzschild metric, is similar to that obtained in [13].
However, the values of θ and x max now depend on the frequencies ω and ω e , thus introducing the effects produced by the presence of plasma and its distribution.The existence of a second antipodal cap identical to the first simply requires the addition of a second component with θ 0,2 = π − θ 0,1 , where F T is the total observed flux, and F is the flux produced by each of the caps, given by Eq. (44).

A. Summary
For the reader's convenience, we list below the steps to be followed to obtain the pulse profile of a neutron star with circular and uniform, single or antipodal emission caps.
1. Choose the metric elements A(r), B(r) and C(r).
2. Set R, q * , and h, taking into account their respective limitations.
3. Select the angles between the axis of rotation and the line of sight ξ (0 ≤ ξ ≤ π/2), between the rotation axis and the center of the cap χ (0 ≤ χ ≤ π) and half-opening of the cap θ c (0 ≤ θ c ≤ π).

Choose a surface emission function
5. Calculate the maximum impact parameter x max [Eq.( 19)].
8. Repeat the procedure sweeping values of γ p between 0 and π (the single-cap profile is symmetric with respect to π) to obtain the single-cap pulse profile.
9. In the case of presence of a second cap identical and antipodal to the first one, add the flux of the latter according to Eq. ( 45) by sweeping the values of γ p between 0 and π/2 (the profile of the antipodal caps is symmetric with respect to π/2).In this case, we take χ 1 such that 0 ≤ χ ≤ π/2.

B. Results
We implemented in Fortran 90 a program that uses Simpson's rule to solve the integral in Eq. (44).
To obtain the luminosity curves we must plot the observed flux F as a function of the star period γ p , recalling the time dependence of θ 0 .This procedure was performed for a variety of frequency ratios = ω e (R)/ω(R), charges q * and radius ratios R/r S , with r S = 2M being the radius of the event horizon for the Schwarzschild metric.The most compact radius ratio used is such that it presents multiple images of the surface, which is now completely visible.The plasma density distributions used are the same as in the previous section, resulting in the refractive indices given by Eq.( 28) and Eq.(29).We consider polar caps with angular halfopening θ c = 5 o , single or antipodal, assuming isotropic emission (f B (δ) = 1).The specific value of I 0 only affects the scale of the pulse profile leaving its morphology unchanged, so I 0 = 1 was taken.
First, a single cap was considered in an orthogonal configuration (χ = ξ = π/2).The results are shown in Figures 6 and 7.
Highly relativistic stars with radii R/2M = 1.675 and 2 were considered.The values |q * | = (0.25M ) 2 and (0.75M ) 2 were taken.In turn, each figure corresponds to a particular plasma distribution, and compares the results obtained for different frequency ratios, with values = 0.00, 0.30, 0.60, 0.90 and 0.99.The first thing to note is that, for extremely compact stars and under certain conditions, the observed flux is much higher when the emitting cap is located at the back of the star (Ωt = π) than when it is located at the front (Ωt = 0).This highly counterintuitive phenomenon is due to the fact that, because the star is so compact, a spot at its back produces multiple images.Thus, when the cap is in this position, its image forms a ring, not too thick but with a large angular aperture, which occupies a bigger part of the observer's sky.As the cap enters into the multiple imaging zone (i.e., for sufficiently large θ 0 ), the luminosity is boosted resulting in the observed pulse profiles.In general, this effect increases with the compactness of the star, canceling out for stars of larger radius whose surface is no longer fully visible, so they do not form multiple images.We also see that the greater the radius of the star, the greater its relative luminosity at Ωt = 0.This is because, being less compact, the trajectories of the light rays are less disturbed by the gravitational field, so they are not so strongly deflected and more of them manage to reach the observer.At the same time, increasing the radius decreases the portion of the period during which the cap is visible.
As already discussed in [15] for the RN case, we see that, as increasing q * , the fraction of the period during which the cap is visible or, alternatively, the magnitude of the peak when multiple images occur decrease, while it increases as decreasing q * .This agrees with the results obtained in the previous section, where we saw that increasing q * decreases the observable surface of the star (period in which it is observable) and vice versa.At the same time, we see that the greater the value of the smaller the observed flux and the faster it decays to zero.Again, this is in agreement with what was observed in the previous section, where increasing the plasma density (or decreasing the observation frequency) reduced both the apparent size of the star and its observable surface.On the other hand, we note that the flux at Ωt = 0 is higher for q * > 0, while as Ωt increases, decreases more rapidly than for q * < 0, which presents a higher peak at Ωt = π but decaying to zero later.These differences are accentuated by increasing |q * |.
Contrary to what happens when considering plasma densities that decay asymptotically, for constant distribution, we see from Fig. (7) that higher result in higher flux.However, this is a mathematical artefact due to the The purple, blue, green, red, and brown lines correspond to = 0.00, 0.30, 0.60, 0.90 and 0.99 respectively.Solid lines correspond to q * > 0 while dashed lines correspond to q * < 0. In the left column we see the characteristic peak due to the cap in opposition generating multiple images.If the star is not compact enough, the flux goes to zero when the cap leaves the zone of visibility.In the top row, for small values of q * the curves obtained for q * > 0 and q * < 0 are almost identical.In the bottom row, we see that for q * < 0 the initial flux is smaller than for q * > 0 but decays more slowly.
Flux (•100) Ωt/π Ωt/π FIG.7: Same configuration as in Fig. 6 but for h = 0. We see that, for a constant plasma density, low frequencies do not satisfy the propagation conditions in the vicinity of the star.
normalization chosen for the emitted flux in Eq. ( 40), where I em was set to be inversely proportional to n 2 0 , and therefore divergent for large enough.Indeed, from Eq. ( 39) we can see that, for a physically realistic finite value of I em , I obs should vanish as grows.In addition, we see that for large enough frequency ratios that do not satisfy the propagation condition do not produce an observable pulse profile.
Let us repeat the previous procedure by adding an antipodal emission cap identical to the first one to study its impact on the pulse profile morphology.We keep the Ωt/π Ωt/π FIG.8: Pulse profile for antipodal caps.h = 3, ξ = χ = π/2.Purple, blue, green, red and brown lines correspond to = 0.00, 0.30, 0.60, 0.90 and 0.99 respectively.Solid lines correspond to q * > 0 while dashed lines correspond to q * < 0.
In the left column we see the characteristic peak due to the cap in opposition generating multiple images.When no cap produces multiple images or is out of the zone of visibility, the flux remains almost constant.In the top row, we see that for small values of |q * | the curves obtained for q * > 0 and q * < 0 are almost identical.In the bottom row, we see that the light curves associated to q * < 0 differ from those in q * > 0, as observed in Fig. (6).
Flux (•100) Ωt/π Ωt/π FIG.9: Same configuration as in Fig. 8 but for h = 0. We see that, for a constant plasma density, low frequencies do not satisfy the propagation conditions in the vicinity of the star.
orthogonal configuration of the angles (χ = ξ = π/2) and the same values for R, and q * and plasma profiles.The results are shown in Figs. 8 and 9.We see that now, given the antipodal cap, the increase in brightness due to the lensing effect also occurs at Ωt = 0, resulting in a flux peak at the initial instant, which will be more pronounced the higher the compactness of the star.In addition, we note that at later times, when none of the caps are either producing multiple images or hidden in the star's non-visibility zone, the resulting Purple, blue, green, red and brown lines correspond to = 0.00, 0.30, 0.60, 0.90 and 0.99 respectively.Solid lines correspond to q * > 0 while dashed lines to q * < 0. The rows, from top to bottom, correspond to Beloborodov classes I, II, III and IV.Both the q * and modify the morphology according to the metric.As q * increases, R decreases.
flux remains largely constant despite the change in the position of the caps.If the star is not compact enough to produce multiple images, there is a minimum with a small value in the flux when the two behaviors are spliced, but the magnitude does not vary so much from one regime to the other.On the other hand, when the surface of the star is fully visible, there is no minimum when switching from one behavior to the other, and the flux intensity is much higher when multiple images are produced.Beyond the addition of a second cap, the qualitative description of the results is otherwise identical to that for the single caps.
Up to this point, these results are consistent with those found in [13] and [18], considering the effects of compactness and charge on the pulse profile.
Finally, we vary the angles χ and ξ to study the frequency-dependent morphologies of the pulses for a variety of angles between the line of sight and the rotation axis and between the rotation axis and the position of the emitting cap.Employing the classification scheme proposed in [29], in Figs. 10 and  horizon, given by r S = 2M for the Sch metric, while for RN -like metrics, we have Note that both radius ratios, R/2M and R/r h , can be used to define the compactness of the star, depending on whether we want to express it in terms of its mass or the radius of its event horizon.In this work, we will use both alternatives for different purposes.
In the Beloborodov's scheme, class I has a primary cap that is always visible, while the antipodal is invisible at all times.In class II, the primary cap is always visible and the antipodal is only visible for part of the time.Class III are generated when each of the caps can be visible or invisible depending on the phase and class IV is when both caps are visible at all times, producing a relatively constant profile.The angular configurations (χ, ξ) = (π/9, π/6), (π/6, π/3), (π/3, 4π/9) and (π/9, 4π/9) were used, corresponding respectively to the four classes just introduced when = 0 and q * = 0. Note that, in order to represent all the classes, the star must have a low compactness, thus containing a region of its surface out of reach of the observer, where the caps are invisible.With this in mind, we use R/r h = 3.
Note that the considered stars do not have the same radius for q * < 0 and q * > 0. It can be seen from Eq. ( 46) that r h is smaller for q * > 0 than for q * < 0, so, by keeping the R/r h ratio identical, the surface radius R in q * > 0 will be smaller.
By observing all these results, it can be seen that both the metric model, the plasma frequency and the charge have influence on what will be the final classification of the pulse profile.We note that class I profiles present sinusoidal like oscillations around a positive mean value greater than the oscillation amplitude.Classes II and III combine periods of sinusoidal oscillations and periods of constant flux.Finally, class IV profiles are almost constant, presenting slight deformations near the extreme values of the phase.

V. ANALYTICAL APPROACH
In a static gravitational field with spherical symmetry, the exact deflection angle β = θ − δ is given by an elliptic or similar integral, corresponding to Eq. ( 18), which generally has no analytical solution.The formalism employed in the previous sections correctly describes the physics of light rays in the vicinity of neutron stars of any size, where the deflection angle β can become larger than π/2.
However, the integrals to be solved can be complicated and require a high computational cost.Expressions such as Eq. ( 37) or Eq. ( 44) may not represent a great challenge when one has a single spherical cap that emits isotropically, but the issue becomes more complicated when considering multiple sources of irregular shapes emitting anisotropically, such as the cases discussed in [13, 16,18,30].When working with a more realistic model, solving the flux integral to obtain the pulse profile can become a real inconvenience, and represent several hours of computation.
In this section we will focus on finding analytical approximations that allow us to generalize to plasma environment situations expressions already known for the case of pure gravity and thus simplifying the calculations of luminosity curves and reduce their computational cost, without great loss of accuracy or generality, facilitating a clear understanding of the deflection effects.
Other analytical developments neglecting plasma modifications to the orbits of light rays, (similar to the ones we developed here) can be found in [17] for scalar-tensor theories of gravity, or in [15] for a generic range of metrics.

A. Trajectory Equation
A simple formula that relates δ (the photon emission angle with respect to the stellar surface normal) to θ (the angle that gives the position on the stellar surface from which the emission occurs, given by Eq. ( 18)), and which replaces the elliptic or similar integral with high accuracy, is the so called cosine relation introduced in [29] 5 6 being A the metric element −g tt , which must be evaluated at the radial coordinate where the angles δ and θ are being measured, in this case R. Eq. ( 47), however, is valid for any point along the trajectory, being δ the angle between the position vector of the ray and the tangent to its trajectory at that point.A limitation of Beloborodov's approach is that it does not take into account the presence of plasma.For this reason, we have been forced to imitate its development, introducing now a spherically symmetric plasma distribution according to Eq. ( 23).The influence of the plasma has been retained up to order 2 , obtaining a correction for the Beloborodov cosine relation that can be expressed in general form as where P q * h (r) is a plasma correction factor which depends on the metric (through q * and M ) and the power h in the plasma density profile, Note that Eq. ( 48) tends to Eq. ( 47) in the low plasma density limit or higher observational frequencies ( → 0), so that both descriptions coincide in pure gravity.It can be seen that this factor vanishes (P q * h (r) → 0) as we move far enough away from the star (r → ∞), which corresponds to a plasma density that decays asymptotically to zero.At the same time, as q * goes to zero, P q * h (r) tends to its Schwarzschild value P q * =0 h (r) ≡ P Sch h (r) when q * → 0.
For the reader's convenience, we list below the steps to be followed for obtaining the plasma correction factors P q * h .
2. Extend the results in Taylor's series in terms of , keeping up to order 2 . 5Beloborodov finds this relation for the particular Schwarzschild case, with A(R) = 1−2M/R, however it is possible to verify that at least for the metrics worked here the relation is generalizable with their respective A(R).Note also, that there are several analytical approaches in the literature superior to the one found by Beloborodov, such as the one shown in [15].However, they all agree to first order with Eq. ( 47) 6 For R ≥ 2r S , this Eq.estimates the deflection angle β = θ − δ with high accuracy, with a maximum relative error of order 3% for R = 3r S .Standard pulsar models, on the other hand, predict R ≥ 2r S with a typical R ≈ 3r S .
3. Expand the above expression as a series around cos δ = 1 and express the result as a polynomial.

Evaluate cos θ and perform
Taylor series around cos δ = 1, retaining only until first order terms.
6. Expand the result by Taylor as a power series of ε, retaining up to second order terms.
Eq. ( 48) then gives us the relationship between the angles θ and δ for all r along the ray trajectory.On the other hand, Eq. ( 19) allows us to express δ in terms of the impact parameter b (or x) and the radius r.Therefore, by combining both expressions, we obtain an equation that directly relates θ to r for every point on the photon's trajectory, where the subscript a indicates that it is an analytical approximation.This is an approximate equation for the trajectory, which will allow us to find the set of coordinates (r, θ) that describes the photon trajectory.While we have previously used Eq. ( 18) and ( 19) to describe the parameters of the photon at the instant it leaves the stellar surface, we should clarify that these are also valid for the rest of its trajectory.Thus, θ and θ a correspond to the emission angle of the photon on the pulsar surface only when evaluated at r = R, while otherwise they describe the angular coordinate of the trajectory corresponding to that r.It can be seen that in [29] the relation found is the inverse, i.e., Beloborodov finds an expression for r as a function of θ.Since the plasma corrections introduce new dependencies on r, and since these are, along with the metric functions, dependent on the specific spacetime model being used, we find it convenient (and much easier) to use Eq.(50).
It is now necessary to verify that the new expressions generalize the approximation presented by Beloborodov, improving it substantially when considering a plasma environment.
With this purpose, we plot the θ and θ a curves as a function of δ, comparing the results obtained by numerical integration and by the analytical approximation, in order to get a general idea about the errors produced by our approximation.Since we want to know the "standard" errors introduced by our model, we take R = 4M (cf.footnote (6)) and q * = −(0.50M ) 2 .On the other hand, since the development of the analytical approximation assumes small , we only take values up to = 0.5.The results are shown in Fig. 12.
It can be seen that the deviations of the analytical trajectory are negligible when δ is small and that they 12: θ as a function of δ. q * = −(0.50M ) 2 , R = 4M .The purple, blue, green, red and brown lines correspond respectively to = 0.1, 0.2, 0.3, 0.4 and 0.5.Continuous lines were obtained numerically [Eq.( 37)] while dashed lines were obtained analytically from Eq. ( 48) (with plasma correction).Errors are minimum for small δ and grow with it, overestimating θmax.
increase as δ becomes larger, so studying the error produced for the marginal ray (δ = π/2 and b max ) will be of great interest for our analysis.Let us note that the analytical approximation tends to overestimate the value of θ max , resulting in a larger visible area of the stellar surface.Such a defect will have an impact on the pulse profiles, since a larger visible area implies a larger lensing contribution.Particularly for the constant plasma case (h = 0) there is a range of δ values for which θ is underestimated by the analytical approximation.However, we see that as we get closer to δ = π/2, θ becomes overestimated, reaching in some cases the maximum value allowed by the approach, θ a = π.Let's now analyze the dependence of the error in θ a with for δ = π/2 and q * = −(0.50M ) 2 , comparing the results obtained with and without the plasma correction.Also, we compare the errors in different R, to observe how the approach improves as the compactness decreases.The error curves obtained are shown in Fig. 13.From this figure we can conclude that the analytical approximation introduced in Eq. ( 48) describes more accurately the photon trajectory in plasma environments than the cosine relation proposed by Beloborodov (Eq.( 47)).Except for some specific behaviors obtained for h = 0, where anomalies generally occur, we see that the accuracy improves by up to an order of magnitude near = 0.5.For small values, as expected, both estimates approach until they coincide at = 0. On the other hand, for ≥ 0.5 the accuracy of our method is still superior, although the error committed starts to become excessive.
We see that, for small , the error in θ max strongly depends on the radius R, being considerably smaller for less compact stars.It is remarkable that, for smaller radius stars (R ≤ 4M ), the error decreases with increasing until it reaches its minimum value around ≈ 0.5, where it becomes very similar to the error obtained for larger radius stars.For larger , the error exhibits a weak dependence on R.
As always, the h = 0 case has its peculiarities.First, we see that there are certain peaks where the error decays drastically, which is explained by the fact that at those points the plasma correction goes from underestimating to overestimating the value of θ, causing the error to decay more by coincidence than by merit.We also see that the error stops being computed for ≈ 0.75, which happens because for larger values the propagation conditions expressed in Eq. ( 26) are no longer satisfied.On the other hand, we see that in the most compact stars the error is not computed for small values of , which arises because the argument within cos −1 in Eq. ( 50) is greater than one.
Finally, we study how the error in θ a changes as we move along the trajectory.For this, we take again R = 4M , q * = −(0.50M ) 2 and δ(R) = π/2, comparing the results obtained with and without consideration of plasma corrections for different values of .These results are shown in Fig. 14.Once again, we see that the error decreases significantly when applying the plasma correction.Note also that the accuracy increases progressively as r increases.This is because, in a sense, we are sending the rays from the observer to the star. 7.With this in mind, it is evident that the deviation increases as r decreases since we accumulate errors due to interaction with  48) (with plasma correction) while dashed lines from Eq. ( 47) (without plasma correction).At all times plasma corrections increase the accuracy of the approach.
plasma.We also see that the error is smaller for lower , which is based on the fact that plasma corrections in the Beloborodov formula were performed by conserving only terms of order 2 .These results are consistent with those shown in [29].

B. Pulse Profile
Making use of the analytical expressions introduced in the previous subsection, we will now focus on finding an analytical expression for the observed flux, which was previously calculated by numerical integration according to Eq. (44).For this, we follow the paper of [30], introducing plasma correction and using the specific intensity as was presented in SEq.IV.The procedure is also similar to the one found in [17].With these considerations, we can express the flux differential as8 being I 0 a constant proportional to the surface intensity emitted by the cap, whose value can be taken as required.For example, in [30] a Planckian emission at uniform T temperature is considered, although other choices for I 0 are equally valid.The total observed flux is obtained by integrating the above expression over the entire visible area of the emitting cap, which we denote S V .We remark that this approximation is valid only if I 0 = const over the entire cap surface S V .Thus, the flux results as follows This expression is consistent with those that can be found in [15] and [17] for point caps.
Using the plasma corrected cosine relation (Eq.( 48)) we can calculate d cos δ/d cos θ = A(R)/(1 − P q * h (R) 2 ).Then, expressing cos δ in terms of θ, we arrive at the following expression The flux is then expressed as the sum of two contributions, the first proportional to the surface area and the second to the projected area of the visible part of the emitting region, both modulated by the anisotropic emission function f B (δ).Let us note that, for A(R) → 1 and → 0, this expression is analogous to the Newtonian result, obtaining the projection of the observed surface.The problem of calculating the flux, once the geometry is fixed, is therefore reduced to determining S V and evaluating the following integrals Let us consider first a simple example, a uniform circular cap of half-opening θ c centered at θ 0 .For simplicity, we take θ c ≤ π/2 and 0 ≤ θ 0 ≤ π.The integral in φ in Eq. ( 54), which yields the angular length of the linear intersection between the parallel given by θ and the emission cap, was solved earlier, resulting in the function h(θ, θ c , θ 0 ) which is expressed in Eq. ( 35) [30].Thus, it only remains to find the θ min and θ max limits of the interval of integration in θ for Eq.(54).
To be visible to the observer, a point on the stellar surface must satisfy the condition δ ≤ π/2.The extreme case occurs when photons leave the star with a trajectory tangential to its surface.These trajectories correspond to the maximum observable angle of the star, θ F (we change the notation to avoid confusion with the upper limit of integration, which we will introduce later), which can be obtained from Eq. ( 48) by taking δ = π/2, resulting in from which the following propagation conditions can be deduced 1 − 2A(R) Note that this restricts the compactness of the star to cases where not every point on its surface is visible to the observer, being the limiting case precisely θ F = π.This is in agreement with the limitations we have been imposing, since in the analytical development of the cosine relation r h /R is assumed to be small.
Since any point on the surface of the star with θ > θ F will not be visible to the observer, and the cap extends from θ 0 − θ c to θ 0 + θ c , it can be seen that the limits of integration will be given by [30] being where the arbitrary constant has been taken zero.It can be seen that, if the cap is completely visible, results in which hugely simplifies the calculations.Thus, Eq. ( 61) is only required when a part of the cap escapes the visible area, in which case it must be evaluated in θ F .Introducing now the effective area we can express the observed net flux as 65) For more than one emitting cap, and considering that each one may have a different intensity, the total flux is obtained by simply adding the contributions of each cap, being able to approximate in this way non-spherical and non-homogeneous caps, obtained by combination.
For the reader's convenience, we list the steps to be followed to obtain the analytical pulse profile of a neutron star.
1. Choose the metric elements A(r), B(r) and C(r).
2. Set the parameters R, q * , and h, taking into account their respective constraints.
3. Select the angle between the rotation axis and the line of sight ξ (0 ≤ ξ ≤ π/2), 4. Set for each cap its intensity I 0,i , the angle between the axis of rotation and the center of the cap χ i (0 ≤ χ i ≤ π) and its angular half-opening θ c,i (0 ≤ θ c,i ≤ π).60), ( 61) and ( 63) must be modified, and if the choice of f B leads to the integrals no longer being analytic, this procedure is no longer valid.

Choose a surface emission function
6. Calculate the maximum visible angle on the surface of the star θ F from Eq. (55).
8. For each cap, calculate the limits θ min,i and θ max,i according to Eq. (57).9. Evaluate I 1,i and I 2,i on the limits θ min,i and θ max,i to obtain the integrals I p,i and I s,i as given in Eq. (59).
10. Calculate the effective area of each cap A ef f,i from Eq. ( 64).
11. Sum the contributions from each cap to obtain the total flux F tot at a given phase γ p according to Eq. (66).
12. Repeat the procedure for each value of the phase γ p ∈ [0, 2π] to be evaluated.
We now show the pulse profiles obtained analytically using the approach developed here and compare them with the corresponding numerical result.For this, we will consider pulsars of a typical size with R/r h = 3, charge q * = −(0.75M ) 2 and ≤ 0.5, with two antipodal and identical caps of angular half-aperture θ c = π/36, in the angular configurations (χ, ξ) = (π/9, π/6), (π/6, π/3), (π/3, 4π/9) and (π/9, 4π/9) corresponding (when = 0 and q * = 0) to classes I, II, III and IV defined in [29] (these and other new classes will be discussed in Sec.VII).This was performed for the Schwarzschild and RNlike metric models on the plasma distributions previously employed.The results are shown in Figs. 15 and 16.
We see that the analytically obtained pulse profile significantly approximates the numerical results for frequencies with ≤ 0.3, reaching a high accuracy for = 0.1, while for frequency ratios higher than 0.3 the approximation loses accuracy, although it is still qualitatively reasonable.
We note that the approximation tends to overestimate the flux generated by a single cap, which corresponds to the peaks observed in Figures 15 and 16 at Ωt = 0, π and 2π, or the profiles corresponding to (χ, ξ) = (π/9, π/6), in the top row of each figure.
At the same time, this effect is magnified when both caps are visible.This is due to the fact that the approximation overestimates the maximum angle θ F , resulting in a greater lensing effect and generating an increase in the flux observed for the caps close to "sunset".In this regime, moreover, the approximation presents a nearly constant profile, omitting the subtle ups and downs that can be observed in the numerical flux.It is in the transition between these behaviors (which is directly related to the location of θ F ) where the model is less accurate.Thus, we conclude that the accuracy increases for plasma distributions with lower density and decreases as we consider emission areas with higher θ, tending to overestimate the flux.The pulse profiles analytically obtained here are consistent with those shown in [29] and are morphologically very similar to those in [15].In summary, for class I profiles (top row of Figures ( 15)-( 16) remembering that at the moment we are using the classification system introduced in [29], which will be revised and generalized in Sec.(VII)) we see that the approximation qualitatively respects the morphology of the curves.The class II and III are in general accurate, apart from the already mentioned irregularities that occur mainly at the peaks and at the confluence of the flux of both caps.On the other hand, the class IV are the most irregular, given that as increases, they degenerate into class II or III depending on the value of θ F , which is different depending on the used method.
Regarding the plasma, we see that for h = 0 the approximation now tends to underestimate the flux.
This could be related to the fact that, in general, a constant plasma distribution inverts the lensing effect from divergent to convergent, thus reversing the previously mentioned effects.In any case, the curves are extremely similar to the numerical curves except for a vertical translation.On the other hand, the power-law decay profiles with h = 3 perform reasonably well, given the already mentioned peculiarities of peak accuracy and overestimation of the flux at the confluence of the two caps.

VI. RINGS
It has been mentioned in the literature that emission regions in neutron stars commonly have irregular shapes, sometimes resembling rings or crescent moons [85,86].Solving these cases numerically can be relatively complicated and involve high computational cost, so the analytical approximations introduced in Sect.V represent a helpful alternative that allows us to obtain the pulse profiles in a much shorter time (about an order of magnitude, depending on the optimization of the programs) with acceptable accuracy and without loss of generality.
In the papers by Sotani et.al., [15,16,87,88] analytical approaches are also developed to deal with this type of cases, with special emphasis on ring-shaped emission caps (for more information or guidance on this topic, check the cited papers), although the influence of a plasma environment is completely neglected.
We will demonstrate the utility and power of our model by using it in the resolution of pulse profiles generated by ring-shaped caps.Suppose we have a pulsar with a single ring-shaped emitting cap, whose outer edge has an angular half-aperture θ e while its inner edge has an angular half-aperture θ i .This can be easily obtained by calculating the flux produced by a circular cap of half-aperture θ e and subtracting the flux produced by a circular cap of half-aperture θ i , both caps with the same position vector Ĉ, In Figures 17 and 18 we show the pulse profiles generated by single ring-shaped emission caps with θ e = 35 o and 70 o respectively, for different values of θ i .For this, the parameters R = 4M , q * = −(0.75M ) 2 , = 0.15 and h = 3 (see Eq. ( 23)) were used on the RN-like metrics, for the χ = ξ = π/3 and χ = ξ = π/2 configurations.The fluxes obtained were normalized according to their maximum value.
Although the magnitude of the brightness decreases noticeably with increasing value of θ i (this is not visible in the plots given the normalization), the morphology of the profile minimally differs from that of the circular cap, at least qualitatively.The resulting flux presents the sinusoidal profile that we expect from a single cap.
By removing part of the cap (i.e., by increasing θ i ) the amplitude of the oscillation decreases, tending to flatten the curve.This is because the contribution from the center of the cap that we are subtracting has less weight the farther it is from the line of sight, so the minimum value of the flux will decrease less than the maximum value (this looks inverted in graphs due to normalization, where F max remains constant while F min increases).
On the other hand, we see that this effect on the morphology is only noticeable in relatively large rings (θ e ≈ 70 o ), while for small rings (θ e ≈ 35 o ) we see nearly no difference between the different θ i .This is because, as the ring is small, it occupies an environment where the contributions of each point to the flux are similar, so removing part of the cap proportionally decreases the net flux without affecting the profile morphology.By considering small rings approaching the point cap approximation, so that what we are doing could be compared to subtracting two point caps of different intensity.Moreover, we see that by changing the configuration (ξ, χ) the behavior is as expected.For (ξ, χ) = (π/3, π/3), the ring does not move too far away from the front of the star, keeping the value of θ 0 low, so the minimum flux is relatively high.In contrast, for the configuration (ξ, χ) = (π/2, π/2) the ring is opposite the observer at Ωt = π, so the observed flux will be much lower, since much of it is located in the invisible region of the star, or at too high θ values.
The analysis developed in this section is in full agreement with that found in [88], where similar results were obtained.We see that there is little difference between the profiles obtained by these metrics.
We emphasize that this is only one of the possible applications of our analytical approach.With it, the calculation of the pulse profile produced by caps of all kinds of irregular shapes or with different surface temperatures can be simplified, making the study of the properties of these systems more accessible and facilitating their understanding.

VII. CLASSIFICATION
When considering two emitting caps which are no longer identical and which are also not in an antipodal configuration, the Beloborodov system [29] is insufficient  to classify the different observable profiles.As the caps are distinguishable, the profile is no longer symmetric with respect to which cap is which, so by inverting their positions we should add a new class.On the other hand, by abandoning the antipodal configuration we give rise to profiles where, for example, neither of the caps is visible during the entire phase, or one is visible only occasionally while the other is always invisible.These phenomena can also be achieved by increasing the plasma density, resulting in new profile classes.We have then a total of nine profile classes for two distinguishable non-antipodal caps, the four classes introduced by Beloborodov, the inversions of his classes I and II, two classes with partial and null visibility caps and a null profile class.For convenience, we organize this new classification according to whether the caps have total (i.e., the cap is visible at all times for the observer), partial (it is visible at some times and not at others) or null visibility (it is not visible for the observer at any time).This classification is explained in Table II, where the correspondence with Beloborodov's classes is also indicated.Following, in Figures 19 to 21 are shown the location maps of the classes in the ξ −χ p plane, being χ p the angle between the rotation axis and the center of the primary cap, for different values of ∆χ = χ p − χ s , being χ s the angle between the rotation axis and the secondary cap.We take R = 6M , q * = −(0.5M ) 2 and = 0.3 for Sch and RN-like metrics, and also = 0 for Sch.The primary and secondary caps have angular half-apertures θ c,p = 3 o and θ c,s = 10 o respectively.
Analyzing the class location maps, we see that the different classes are separated by straight lines marked in red on the graphs, which correspond to the values of ξ and χ p for which one of the caps, at its farthest point from the line of sight (i.e., at its maximum elongation at θ), is partially in the zone of visibility and partially in the zone of invisibility, being in this way transition zones for caps from total to partial visibility, or from partial to null visibility.Thus, the thickness of these lines corresponds to the angular size of the caps (this issue is also discussed in [88] and [16], where extended caps are considered).Moreover, the location of these eight red lines that we see in each graph are given by the following equations being the maps symmetric with respect to the straight lines ξ = 0 and ξ = 180 o , that is, Class(π + δ, χ p ) = Class(π − δ, χ p ).On the other hand, we have that Class(ξ, π + δ) = Class(π − ξ, δ).Given these symmetries, we see that in each graph two slanted rectangles are formed, whose left vertex is on the χ p axis and whose upper left and lower right edges are shorter than the Given the differences between the numerical and analytical values of θ F , there is a small difference in the location of these lines which, in extreme cases, may result in the emergence of classes that should not be there, or in their omission.However, the morphology of the class maps obtained by the numerical and analytical methods are generally indistinguishable.
To complete the discussion, in Figure 22 we show some representative pulse profiles for each classes introduced in Table II.
However, when we look at Figure 23, where pulse profiles with caps in different configurations are included, we see that, contrary to what happened in Beloborodov's classification, by allowing the caps to be not only not antipodal but also not on the same meridian (i.e., ∆γ p = γ p,p − γ p,s = 0) we have that pulses of the same class present distinctly different profiles in terms of their morphology.
To solve this problem, we have elaborated the following subclassifications based on the scheme presented in Table II, taking into account only some of the most evident morphological differences between profiles of the same class.
Class I: • a: The caps are in phase opposition, so that the decay in brightness of one counteracts the increase in the other, resulting in a constant profile.• b: The caps are in phase conjunction or out of phase, so that their luminosity is not counterbalanced, resulting in a non-constant profile.
Class II and IV: • a: The caps are out of phase, so that the profile has two distinguishable peaks.
• b: The caps are out of phase, so that the profile has a single peak and periods of constant flux.
• c: The caps are in conjunction, so that the profile presents a single peak with two characteristic slopes.Class V: • a: The caps do not share a period of visibility, resulting in two peaks separated by two periods of zero flux.
• b: The caps share a visibility period interval and have intervals where they are the only one visible, resulting in two overlapping peaks and a single zero flux period.
• c: The visibility period of one of the caps is completely contained within the visibility period of the second, resulting in a single peak and a single period of zero flux.
To close this section, in Figure 24 we show some characteristic pulse profiles for the subclasses just mentioned.

VIII. CONCLUSIONS
A correction to the cosine relation presented in [29] [Eq.(47)] was found to now include the refraction effects produced by a plasma environment whose density is given by a power law, arriving at Eq. ( 48), whose exact form depends on the plasma profile h, the metric model and the quotient between frequencies .This formulation was applied to the study of pulsars, allowing a much simpler treatment of the problems seen above by replacing the complicated numerical integrations [30], giving simple expressions for the deflection angle β and for θ in terms of r and the impact parameter x.The study of the errors in θ showed that our approximation is effectively an accurate improvement of the Beloborodov's formula when considering plasma environments, presenting different levels of accuracy depending on the compactness of the star, the density and distribution of the plasma or the metric model.In general, the errors are acceptable for stars of low compactness (R ≥ 6) with low frequency ratios ( ≤ 0.3).We emphasize that the analytical approximations obtained can be used to considerably simplify the calculation of pulse profiles produced by circular and homogeneous caps, and even combinations of these, in environments where the plasma has a significant impact on light deflection, under realistic compactness and electric charge conditions, thus facilitating the study and understanding of these systems.In our case, we gave a demonstration of its use to solve the profile produced by ring-shaped caps [16,88], but this is only one of the wide range of models that can be treated with this approach.
Given these new facilities, we were able to consider neutron stars with two non-antipodal, non-identical caps, which gave rise to new pulse profiles that did not fit Beloborodov's classification scheme.For this reason we developed our own classification system, taking into account the new peculiarities of the model.The main effect was the enlargement of the non-specular area of the class map in the χ − ξ plane, i.e., the area that cannot be obtained as a reflection or translation of the others.At the same time, due to the flexibility in the model configuration, many profiles of the same class have distinctly different morphologies, so we have extended the classification into different subclasses.
We also mention that one or more of the assumptions made in our model are likely to be violated in nature, so it is important to know their range of validity.It should be noted that the formalism presented here can only be applied to slowly rotating neutron stars.Otherwise, metrics such as Kerr's or similar ones must be considered.Moreover, the caps of a realistic pulsar have no need to emit isotropically according to Lambert's law (f B (δ) = 1), either because the star is covered by an atmosphere, or because the emissivity is strongly constrained to energies below the electron plasma frequency.Realistic emission models predict an angular dependence in the emitted intensity.An interesting possibility is that the emission originates above the neutron star surface, generating emission with δ > π/2 [13].A more realistic model would have to take into account the effect of the strong magnetic fields in the vicinity of the star on the propagation of light rays [89][90][91].A study of this nature has been made in the past to study light rays in plasma environment around black holes [90].
Even though our formulas and plots are expressed in terms of dimensionless coefficients to maintain as much universality as possible, regardless of the specific details of the model, in order to have an rough estimation of the order of magnitude in the involved parameters, let us consider a NS rotating at the angular velocity Ω and with an exterior dipole magnetic field B. We also assume that the particle number density N (R) of the magnetospheric plasma at the surface of the star is proportional to the Goldriech-Julian density N GJ (R), i.e., with κ the so called multiplication parameter, taking values in the range 10 3 −10 5 [92,93] and B the characteristic value of magnetic field strength at the surface R [42], with I = 2/5M R 2 the Newtonian moment of inertia of the star (assumed as a rigid sphere).Moreover, due to the relativistic motion of electrons and positrons along the magnetic field lines, the angular plasma frequency is decreased by a factor γ−3/2 [94], where γ is the average of the Lorentz factor γ = E/(m e c 2 ) of the mentioned particles, which takes a minimum value γ ≈ 10 2 , ( [93][94][95]).Note also that the frequency ν ∞ of the electromagnetic wave as measured for an asymptotic observer is related to by: with K = 1 − 2 M R + q * R 2 the gravitational redshift factor.In that setting, a NS of mass M = 1.8M with a rotation period P = 1s varying as Ṗ = 10 −15 has the associated parameters shown in Table (III).Corrections to the Goldriech-Julian density when one consider a RNlike metric are given by [78] (see their Eq.(9)).It can be checked that the correction terms are of order O(1), therefore the estimated values for the parameters in RNlike metrics remain of the same order of magnitude as those presented in Table(III).
On the other hand, a strong magnetic field leads to an anisotropic emissivity with a preferential direction along the field [29].In such a situation, one surely should abandon semi-analytic models and perform a full numerical integration of the ray path equations taking into account a dielectric permittivity tensor.Despite all these restrictions on the range of validity and being limited to extreme cases, the present work constitutes a simple model that allows us to describe and approximate to a large extent the photon dynamics around compact objects embedded in plasma media, achieving a better understanding and comprehension of the different morphologies in the light curves and their dependence on different physical and geometric parameters involved.
Finally, we would like to emphasize that the generalized Beloborodov formula that we have obtained can also be applied to the study of polarization of light rays in plasma environments, as well as to determine luminosity profiles associated with accretion disks around black holes and other possible compact objects.It is worth mentioning that in the last years alternative formulas to Beloborodov's (many of them empirical) that improve it have also been proposed in the literature [96][97][98][99][100][101][102].It would be desirable to be able to obtain its analogs in situations where a plasma medium (not necessarily a cold plasma) is present.Purple, blue, green, green, red and brown lines correspond to = 0.00, 0.30, 0.60, 0.90 and 0.99 respectively.The solid lines correspond to q * > 0 while the dashed lines correspond to q * < 0. Ωt/π Ωt/π FIG.28: Pulse profile for antipodal caps.Exponential decay, ξ = χ = π/2.Purple, blue, green, green, red and brown lines correspond to = 0.00, 0.30, 0.60, 0.90 and 0.99 respectively.The solid lines correspond to q * > 0 while the dashed lines correspond to q * < 0. Purple, blue, green, green, red and brown lines correspond to = 0.00, 0.30, 0.60, 0.90 and 0.99 respectively.The solid lines correspond to q * > 0 while the dashed lines correspond to q * < 0. The solid lines correspond to q * > 0 while the dashed lines correspond to q * < 0.

FIG. 4 :
FIG.4: Angular configuration between the observer L, the rotation axis of the star R and the center of the emitting cap Ĉ, of angular half-aperture θc.

FIG. 5 :
FIG. 5: The figure on the left shows the dependence of the observed position of the cap θ0 on the phase of the pulsar Ωt.In the figure on the right we see a graphical representation of the φ integral, φ b , at a given θ.

FIG. 6 :
FIG.6: Pulse profile for single caps.h = 3, ξ = χ = π/2.The purple, blue, green, red, and brown lines correspond to = 0.00, 0.30, 0.60, 0.90 and 0.99 respectively.Solid lines correspond to q * > 0 while dashed lines correspond to q * < 0. In the left column we see the characteristic peak due to the cap in opposition generating multiple images.If the star is not compact enough, the flux goes to zero when the cap leaves the zone of visibility.In the top row, for small values of q * the curves obtained for q * > 0 and q * < 0 are almost identical.In the bottom row, we see that for q * < 0 the initial flux is smaller than for q * > 0 but decays more slowly.

2 Flux 6 , ξ = π 3 FluxFIG. 11 :
FIG.11: Same configuration as in Fig.10but for h = 0. We see that, for a constant plasma density, low frequencies do not satisfy the propagation conditions in the vicinity of the star.

6 , ξ = π 3 FluxFIG. 15 :
FIG.15: Analytical pulse profile for two identical, antipodal caps in Schwarzschild metric, R/r h = 3, θc = 5 o .Purple, blue, green, red and brown lines correspond to = 0.1, 0.2, 0.3, 0.4 and 0.5 respectively.Solid lines correspond to the numerical result, dashed lines to analytical.The rows, from top to bottom, correspond to Beloborodov classes I, II, III and IV.For h = 3, the flux is overestimated at maximums and underestimates at minimums.For h = 0 is in general overestimated.

Class
Main c. vis.Secondary c. vis.

TABLE I :
Family of metrics around compact objects with a gravitating mass M and a charge q * , in static, spherically symmetric and asymptotically flat spacetimes.

TABLE II :
New classification system for pulsars with two distinguishable non-antipodal caps, compared to Beloborodov's.

TABLE III :
Characteristic values of the parameters involved in the description of the plasma magnetosphere for a assumed mass of 1.8 solar masses of the NS.