Constraining millisecond pulsar geometry using time-aligned radio and gamma-ray pulse profile

Context. Since the launch of the Fermi Gamma-Ray Space Telescope, several hundred gamma-ray pulsars have been discovered, some being radio-loud and some radio-quiet with time-aligned radio and gamma-ray light curves. In the second Fermi Pulsar Catalogue, 117 new gamma-ray pulsars have been reported based on three years of data collected by the Large Area Telescope on the Fermi satellite, providing a wealth of information such as the peak separation ∆ of the gamma-ray pulsations and the radio lag δ between the gamma-ray and radio pulses. Aim. We selected several radio-loud millisecond gamma-ray pulsars with period P in the range 2–6 ms and showing a double peak in their gamma-ray profiles. We attempted to constrain the geometry of their magnetosphere, namely the magnetic axis and line-of-sight inclination angles for each of these systems. Method. We applied a force-free dipole magnetosphere from the stellar surface up to the striped wind region – well outside the light cylinder – to fit the observed pulse profiles in gamma-rays, consistently with their phase alignment with the radio profile. In deciding whether a fitted curve is reasonable or not, we employed a least-square method to compare the observed gamma-ray intensity with that found from our model, emphasising the amplitude of the gamma-ray peaks, their separation, and the phase lag between radio and gamma-ray peaks. Results. We obtained the best fits and reasonable parameters in agreement with observations for ten millisecond pulsars. Eventually, we constrained the geometry of each pulsar described by the magnetic inclination α and the light-of-sight inclination ζ. We found that both angles are larger than approximately 45◦.


Introduction
The launch of the Fermi Gamma-Ray Space Telescope in June 2008 revolutionised our knowledge of the physics of pulsars when it discovered hundreds of new gamma-ray pulsars, among them many millisecond ones (Abdo et al. 2013), a surprising result at that time. The second gamma-ray pulsar catalogue (2PC) contains 117 pulsars. The most recent catalogue, 3PC, will soon be reported, now containing 300 pulsars. These gamma-ray pulsars are either radio-loud or radio-quiet. If seen in both wavelengths, their radio and gamma-ray light curves provide a wealth of information about their magnetospheric topology and emission sites close to the neutron star surface and in regions near the light-cylinder radius (r LC ), respectively. The particles are likely to corotate with the neutron star within r LC and flow beyond r LC as a pulsar wind, that is, a relativistic magnetised outflow of pair plasmas, forming a current sheet known as the striped wind (Coroniti 1990;Michel 1994).
Particle trajectories and photon-production mechanisms within the magnetosphere and the electromagnetic field topology itself are all intertwined with each other. Realistic modelling of these quantities therefore necessitates the computation of the plasma-radiation field interaction in a self-consistent manner. Several promising approaches in this direction have been undertaken by many authors using either a fluid descriptionsee for example Spitkovsky (2006), Bai & Spitkovsky (2010), Kalapotharakos et al. (2012), Pétri (2011Pétri ( , 2012aPétri ( , 2018) -or a particle description; see for example (Cerutti et al. 2016). The vacuum-retarded dipole solution (Deutsch 1955) entails vanishing particle density (no particle to accelerate) while possessing a huge electric field component parallel to the magnetic field inside the magnetosphere. Force-free electrodynamics (FFE) solutions, on the contrary, must have large values of the charge density, which shorts out the parallel component of the electric field and precludes the particle acceleration. Some resistivity or dissipation is required to avoid a full screening of the electric field, allowing for efficient particle acceleration to ultrarelativistic speeds.
Another complication while modelling pulsar electrodynamics stems from the magnetic obliquity with respect to the spin axis of the pulsar. Initial attempts to understand the physics underlying the pulsar magnetosphere commonly focused on the aligned rotator assumption. The reason was that finding a solution for the aligned case is much easier than for the more realistic oblique case. Contopoulos et al. (1999) produced successful numerical simulations to encounter FFE for the first time for the aligned rotator. Several authors constructed a force-free pulsar magnetosphere by considering time-dependent simulations for the oblique case (Spitkovsky 2006;Kalapotharakos & Contopoulos 2009;Pétri 2012b). Because the simulations would otherwise be very time consuming, in these previous works, a large ratio of R/r LC = 0.2-0.3 was used corresponding to an unrealistically fast pulsar with periods of ∼1 ms. Some authors proposed models including an 'electrosphere', in which a differentially rotating plasma flows around the equatorial plane, separated from the oppositely charged domes by huge gaps (Krause-Polstorff & Michel 1985;Pétri et al. 2002;McDonald & Shearer 2009).
For the first time, Watters et al. (2009) performed a detailed study of pulsar light curves with an individual assessment of polar cap, outer gap, and slot gap models. Since then, several authors have contributed to the goal of producing a complete atlas of pulsar light curves for different obliquities α and inclination angles ζ. Venter et al. (2009), Johnson et al. (2014), and Pierbattista et al. (2015) adopted models based on the vacuumretarded dipole solution, first given by Deutsch (1955), and later extended to include rotational sweep-back magnetic field lines (e.g., Dyks et al. 2004). Contopoulos & Kalapotharakos (2010) studied the emission of curvature radiation from the current sheet that develops around the equatorial plane in the force-free regime, which they called pulsar synchrotron, and were able to reproduce the high-energy light curves emitted by pulsars.
Estimation of the radio and gamma-ray time lag δ and gamma-ray peak separation ∆ for those pulsars showing doublepeaked light curves and application to the observed light curves of several gamma-ray pulsars provided in the first Fermi Catalogue Abdo et al. (2010a) was carried out by Pétri (2011) who used an oblique split monopole solution (Bogovalov 1999) for the field structure and a simple polar cap geometry for the radio counterpart. With the study presented here, our aim is to understand the characteristics of the gamma-ray and radio light curves by employing the FFE prescription where we also take into account the real effect of particle flow within the dipolar magnetosphere on the field structure. Contopoulos & Kalapotharakos (2010) also compared the predicted relation δ − ∆, assuming emission near the current sheet, to the calculated relations based on the first pulsar catalogue. Kalapotharakos et al. (2012) focused on the orthogonal rotator case, taking into account a large range of conductivity in dissipative pulsar magnetospheres. In Kalapotharakos et al. (2014), the orthogonal condition was relaxed to the whole range of magnetic obliquity to produce δ and ∆ measured by the observations. By comparing the observed ∆ − δ diagram with their estimations from the model, these authors favoured the cases with high conductivities; as the conductivity increases, the system converges to the FFE limit. Chen et al. (2020) recently suggested a force-free magnetospheric model and presented radio, X-ray, and gamma-ray light curves in agreement with the observations of the millisecond pulsar (MSP) PSR J0030+0451, providing an estimation of the magnetic inclination of the pulsar of ∼80 • . More recently, Kalapotharakos et al. (2020) proposed the FFE model by taking into account an off-centred dipole field and quadrupole components, the Kerr-metric for the ray-tracing of photons from a distant observer to the hot spots on the surface of the star.
Fitting the radio-loud gamma-ray pulsar high-energy light curve heavily relies on a good time alignment constraint between radio and gamma-rays. Radio profiles are often fitted by a single or multiple Gaussian components. In our approach, outlined below, the radio beam emitted from the polar cap is symmetric, and the intensity profile is defined by a Gaussian centred at the magnetic axis; the fiducial phase for the radio light curve of a pulsar corresponds to the phase at which the magnetic axis lies in the plane defined by the rotation axis and observer's line of sight. For normal (young) pulsars with periods above about 100 ms, the rotating vector model (Radhakrishnan & Cooke 1969) is usually exploited to constrain the fiducial phase and emission height.
However, this method is inappropriate for MSPs with periods of about a few milliseconds because of the complex structure of their pulse shapes. This is mainly due to the fact that radio emission emanates from altitudes of a substantial fraction of the light-cylinder radius where relativistic effects are important (aberration, magnetic sweep-back, plasma currents), and very close to the stellar surface where multipolar magnetic field components are expected to be important. Therefore, the phase-zero determination of the radio light curve for MSPs is problematic, not to mention that the radio duty cycle can also be very large. Therefore, for the definition of the phase zero of ten pulsars investigated in this work, two distinct methods were used in the second Fermi Pulsar Catalogue: one taking into account the fiducial point at peak intensity; and the other using the fiducial point placed at the point of symmetry of the prominent radio peak. PSR J0030+0451, PSR J0102+4839, PSR J0437-4715, PSR J1614-2230, PSR J2017+0603, and PSR J2043+1711 belong to the former category while PSR J0614-3329, PSR J1124-3653, PSR J1514-4946, and PSR J2302+4442 belong to the latter category (see Table 8 in Abdo et al. 2013).
In this study, we adopted the same fiducial phases as those given in the catalogue. Independently of the method used in the catalogue, all zero phases correspond to the alignment of the magnetic axis with the line of sight in our model.
All these recent studies are very promising for the constraint of the geometry of the pulsar and are the subject of this paper. The paper is organised as follows. In Sect. 2 we briefly mention the radio emission phenomenology of pulsars and draw attention to the complications in the radio pulse profiles of MSPs. In Sect. 3, we summarise the gamma-ray observations of Fermi/LAT presented in 2PC. In Sect. 4, we describe the forcefree model for the magnetosphere and the striped wind where gamma-ray pulses are assumed to be produced in our model. We explain our fitting method in Sect. 5. Remaining consistent with radio and gamma-ray time-alignment, we reproduced the observed gamma-ray pulse profiles for the sample of MSPs and attempted to determine the geometry of each pulsar. We present our results in Sect. 6. A broader discussion of our method is provided in Sect. 7. We draw conclusions in Sect. 8.

Radio observations and phenomenology
We aim to use radio and gamma-ray observations to find the dipolar magnetic field emission geometry of the pulsar with respect to the observer. The role of the radio observations enters the problem in the following manner. Firstly, one needs to assume that the radio emission arises deeper in the magnetosphere closer to the surface of the neutron star and secondly the emission needs to arise from regions where the magnetic field geometry can be closely approximated with a static dipole. These two assumptions can then be used to treat the radio emission as a proxy for identifying the fiducial dipolar plane of the magnetic field lines; the geometry of the gamma-ray emission can then be solved with respect to this plane. The assumptions hold true for normal-period pulses (i.e. pulsars with periods longer than about 100 ms), where detailed phenomenological pulse profile and polarisation studies suggest that the radio pulsar emission beam comprises a central core emission surrounded by two nested conal emissions (see Rankin 1983Rankin , 1990Rankin , 1993Mitra & Deshpande 1999). The emission arises from regions of open dipolar magnetic field lines below 10% of the light-cylinder radius; see Mitra & Li (2004) and Mitra (2017). These conclusions were based on demonstrations that the observed opening angle of normal pulsars follows the P −0.5 (where P is the pulsar period) behaviour (see, e.g., Skrzypczak et al. 2018) which is as expected from emission arising from the open dipolar magnetic field line region. Secondly, the polarisation position angle (PPA) across the pulse profile follows the rotating vector model (RVM), which is an indication of emission arising from regions of diverging magnetic field line geometry. Further, the effect of aberration and retardation in the linear approximation has been observed in normal radio pulsars -see Blaskiewicz et al. (1991) Pétri & Mitra (2020) -, whereby a shift between the centre of the total intensity profile and the steepest gradient point of the position angle curve is observed. This shift is proportional to r em /r L and is used to obtain the radio emission heights.
However, in MSPs, the normal period pulsar phenomenology cannot be applied. A radio pulse profile does not have any orderly behaviour, and the PPA traverse is significantly more complex and cannot be modelled using the RVM. Only recently, Rankin et al. (2017) tried to model a few MSPs in terms of the core cone model and suggested that the radio emission might arise from a few tens of kilometres above the neutron star surface. However, these estimates are highly uncertain, because in none of the cases can the RVM be fitted to the PPA traverse. Therefore, in observations of MSPs, the location of the radio emission region is not constrained. In this study, we acknowledge this drawback and assume that the emission arises from regions anchored on open dipolar magnetic field lines of about 10% of light-cylinder radius r L .

Gamma-ray observations
The Large Area Telescope (LAT) instrument onboard the Fermi Satellite, launched in June 2008, detects gamma-rays with energies between 20 MeV and ∼300 GeV. The full three years of data used to compile the 2PC in Abdo et al. (2013) was collected via the observations from 2008 August 4 to 2011 August 4. The gamma-ray pulsar search using the known rotation ephemerides of radio or X-ray pulsars led to the discovery of 61 gamma-ray pulsars in the 2PC. Phase alignment of gamma-ray pulses with radio pulses provides vital information about the geometry of the different emission regions. Another method used in the 2PC for discovering gamma-ray pulsars is the blind periodicity search which can provide ephemeris of pulsars relying only on the LAT gamma-ray data. In the catalogue, 36 new gamma-ray pulsars were discovered with this method. In addition, radio searches of several hundred LAT sources by the Fermi Pulsar Search Consortium provided 47 new pulsar discoveries, including 43 MSPs and 4 young or middle-aged pulsars (Ray et al. 2012). One of the most prominent differences of the 2PC from the 1PC is the much higher ratio of MSPs to the normal pulsar population.
Spectra of gamma-ray pulsars are relatively insensitive to the neutron star period, peaking around a cut-off energy of several GeV for all of them. Moreover, the subexponential decrease of the spectra above this cut-off rules out the polar cap scenario for gamma-ray emission because of the strong magnetic opacity at these energies (Abdo et al. 2009). Slot gaps and outer gaps have since been the favoured site of high-energy photon production. While there is no independent observational constraint on the location of the gamma-ray emission, because the radio time lag δ and peak separation ∆ values found in 2PC are supported by extensive force-free and dissipative magnetosphere simulations, these regions are shifted more and more towards the light cylinder and even outside in the striped wind (Pétri 2011). Therefore, in the present study, we use a magnetosphere model (see Sect. 4) where the gamma-ray emission starts from the light cylinder and continues into the striped wind.

Magnetosphere and emission model
As the neutron star is surrounded by plasma made of electronpositron pairs producing the observed radiation, we need an accurate solution for the neutron star magnetosphere. Therefore, the radio and gamma-ray light-curve computations rely on the electromagnetic field structure extracted from the force-free pulsar magnetosphere obtained by time-dependent numerical simulations. Such calculations have been performed by Pétri (2020) using his time-dependent pseudo-spectral code detailed in Pétri (2012b). Specifically, the neutron star radius is set to R/r L = 0.2 corresponding to a 1.2 ms period and the artificial outer boundary is set at r = 7 r L where the light-cylinder radius is r L = c/Ω and Ω = 2 π/P. The numerical solution of the electromagnetic field is therefore expressed as a series in Fourier-Chebyshev polynomials and can be accurately evaluated at any arbitrary point within the simulation box. The striped wind starting from the light cylinder is also self-consistently computed on almost one wavelength λ L = 2 π r L .
The two emission sites, radio and high-energy, are described as follows. Radio photons are expected to emanate from the open field-line regions above the polar caps whereas the gamma-ray photons are produced in the current sheet of the striped wind, outside the light cylinder at r ≥ r L . The emissivity of the striped wind drops sharply with distance, so we only considered emission in the shell r ∈ Abdo et al. (2010a,b) r LC . We extended our previous analysis made by Pétri (2011) where a split monopole wind was assumed to be directly connected to the stellar surface. Now this approximation, which in reality does not hold inside the light cylinder, is replaced by a realistic force-free dipole magnetosphere smoothly linking the quasi-static zone inside the light cylinder to the wind zone outside the light cylinder. Figure 1 shows the sky-maps produced by our striped wind model for six different magnetic inclination angles α. In each panel, the phase of prominent radio peak is set to zero, and gamma-ray light curves are plotted according to a time-alignment shift of the phase for synchronisation with the radio pulse profile.
The gamma-ray pulse width is a combination between the relativistic beaming effect and the depth of the line of sight crossing the current sheet. For infinitely thin current sheet, this width W γ scales as the inverse of the wind Lorentz factor Γ and W γ ∝ 1/Γ as shown by Kirk et al. (2002). If the thickness is finite, it will reach a lower limit reflecting the current sheet thickness scaled to the wind wavelength λ L for ultra-relativistic speeds Γ 1. This has been observed in the Geminga pulsar as reported by Abdo et al. (2010b). The gamma-ray peak intensity for each pulse also depends on the length of the line of sight crossing the current sheet in each pulse. Contrary to the split monopole case, the striped wind geometry is less symmetric because of the distorted dipole due to rotation and plasma effects. This loss of symmetry imprints the associated pulsed profile, which is also less symmetric, showing discrepancies in the peak intensities compared to the split monopole emission.
The sky maps are constructed as follows. First, we determine the polar cap rims by identifying the last closed field lines grazing the light cylinder. These field lines also draw the base of the current sheet in the striped wind. More precisely, radio emission is assumed to occur over the whole polar cap area at r = R (therefore at zero altitudes in our simulation box). As we are not interested in the exact radio pulse profiles, which often obtain a complex shape with multiple components, but only in the location of the maximum intensity taken as phase zero for time-alignment purposes do we modulate the polar cap emission with a Gaussian profile w peaking at the magnetic axis such that A101, page 3 of 10 where µ is the magnetic moment vector and n the unit vector pointing in the direction of the line of sight. This profile is not intended to reproduce exact pulse profiles but simply to easily locate the maximum of the radio peak that sets the phase zero for the radio and gamma-ray light curves. Emission in the striped wind is located around the current which is identified by the location outside the light cylinder where the radial magnetic field component reverses sign. In an ideal picture, this current sheet is infinitely thin, but due to dissipation, it will spread to a certain thickness depending on the plasma condition. The variation of the Gaussian profile does not affect the phase of prominent radio or gamma-ray peaks in our model. It only impacts on the width of the radio pulses. However, an exact fitting of the radio pulse profiles is out of the scope of this work; we only use its peak location to determine the phase zero for time alignment purposes. Moreover, by construction in our model, the radio peak intensity corresponds to the phase of the magnetic axis (north pole or south pole). We stress that only the radio pulse phase is needed, regardless of its shape, to constrain gamma-ray peak phases and thus the geometry of pulsars. The dipolar polar cap size has a half opening angle of θ pc ≈ sin θ pc = √ R/r LC . This corresponds to a radio pulse profile width of W ≈ 3/2 θ pc , which also takes into account the magnetic field line curvature, assuming a static dipole which is a good approximation at the surface as long as R r LC . For our 1.2 ms pulsar computed from the simulations, this corresponds to a width of W ≈ 38 • , but for the sample of MSPs we have fitted below, with a 3 ms period on average, this width shrinks to W ≈ 25 • . Therefore, when the observer line of sight ζ lies within the radio emission beam (α − W ≤ ζ ≤ α + W), a radio pulse will still be detected. We added this constraint in the fitting procedure explained in the following section.

Fitting method
The radio and gamma-ray light curves of the pulsars are produced by numerical calculations assuming a force-free magnetosphere. To produce a complete sky-map, we trace the inclination angle between the magnetic dipole moment and the rotation axis from α = 0 • to 90 • with 5 • resolution, for the inclination of the line of sight with respect to the rotational axis (ζ) varying from 0 • to 180 • with 2 • resolution. The maximum intensities for radio and gamma-rays are independently normalised to unity. We do not model the particle distribution within the magnetosphere with a self-consistent physical approach. In that sense, our model is purely geometric. The gamma-ray intensity is a free parameter of the model while the form of the light curves is obtained by solving the equations numerically for the force-free electrodynamics (see Sect. 4). We are mainly interested in the time lag δ between the radio and the first gamma-ray peak, the gamma-ray peak separation ∆ and the ratio of the peak amplitude. These are geometric constraints, and hence we cannot model the exact pulse profile shapes.
In our analyses, we take into account the radio pulse profiles provided in the 2PC to reproduce gamma-ray pulse profiles calculated by our model with the consistent time lags between gamma-ray and radio peaks. The radio observations of PSR J0030+0451, PSR J0614-3329, PSR J1614-2230, PSR J2017+0603 and PSR J2302+4442 by the Nançay Radio Observatory (1.4 GHz band), of PSR J0102+4839 (1.5 GHz band) and PSR J1124-3653 (0.8 GHz band) by the Green Bank Telescope, of PSR J0437-4715 and PSR J1514 by the Parkes Radio Telescope (1.4 GHz band), and of PSR J2043+1711 by the Arecibo Observatory (1.4 GHz band) were used in the 2PC; see Abdo et al. (2013) and the ATNF Pulsar Catalogue 1 (Manchester et al. 2005).
We set the time of maximum radio intensity to phase zero and shifted the gamma-ray profile accordingly. In our model, there is a lag δ between the gamma and radio peaks up to 0.5 in phase. For ζ and α values producing double-peaked gammaray light-curves in the direction of the observer, we measure the phase difference of these peaks (peak separation, ∆). The expected peak separation measured from our model for a given set ζ − α can be seen in Fig. 2. The maximum peak separation is 0.5 in units of the pulsar period P. We note that the distribution of the peak separation in the force-free model is compatible with that found from the asymptotic split-monopole solutions (Pétri 2011).
For each source, in order to determine the goodness of fits to the observed gamma-ray light curve, we use a simple leastsquare procedure, as follows, and try to constrain the reasonable α and ζ values that produce the best-fitting model curves. The reduced χ 2 is defined by where ν is the degree of freedom, which is equals to the number of data points minus the number of fitted parameters, I obs i and σ 2 i are the observed intensity and associated error of gammarays for the ith phase bin, and I model i is the model prediction at the observational phase bin. As the phase space in the model outputs is evenly distributed from 0 to 1 with intervals of 0.01, we interpolate the model light curve and find the intensity at the phase of observation so that we can apply the least-square analysis to the intensities. Regarding the observation data, we take the lower bounds for each bin in the light curve as the time of observation of the associated intensity. As the difference of lower and upper bounds for each bin is sufficiently marginal, e.g. 0.01 for gamma-ray and smaller for the radio light curves, using the minimum phase of each bin instead of the mean phase does not change the fit quality significantly. Also, we subtract the nominal value for the background (reported in the 100 MeV-100 GeV energy band for each source) from the total intensity in order to remove background noise effects.
We mainly focus on the sources with clearly measured gamma-ray and radio light curves with high signal-to-noise ratios. Of the numerous sources observed both in gamma-rays and radio bands and reported in the secong Fermi Catalogue (Abdo et al. 2013), we pick up pulsars with (1) double-peaks observed in gamma-rays for the complete phase of the pulsar, (2) a well-described prominent peak in the radio profile and clear peak locations, and (3) a gamma-ray peak detection time which is not very close to that of the prominent radio peak in order to obtain their gamma-ray peak phases (localised according to their radio peaks). We should note that the gamma-ray peak properties of pulsars, which do not suit item (3)  not be reproduced by our model. Some freedom in the emission height could add retardation effects, shifting the gamma-ray light curves in the right direction. However, this is at the expense of adding a new free parameter in our model, which we want to avoid in a first attempt to constrain the geometry.
As mentioned earlier, to obtain the geometry for our sample pulsars, we fit the observed gamma-ray light curves with the model light curves using the χ 2 ν minimisation process given by Eq.
(2). We note that in our model there are three factors that influence the geometry: (1) the height ratio of the gamma-ray peaks (amplitude), (2) the gamma-ray peak separation, ∆, and (3) the phase lag between radio and associated gamma-ray peaks, δ. However, bridge emission between the gamma-ray peaks, as often seen in so many observed gamma-ray light curves, is not predicted by the model. Therefore, while finding a best fit and calculating a χ 2 ν value, although we take into account off-peak gamma-ray data together with bridge emission between gammaray peaks, we aim to retrieve the light curves best suited to the three above conditions, ∆, δ. and peak amplitude ratio. Therefore, we still regard χ 2 ν values for the best fit of some pulsars as adequate, even though they are not close to unity. We also examined a weighted fit for the peaks, such as calculating χ 2 ν only around peaks or giving more weight to the data around peaks which do not lead to the estimation of different geometry. We therefore apply the unweighted fitting method for all pulsars here with one exception (see Table 1).

Millisecond pulsars with time-aligned radio and gamma-ray pulse profiles
In the framework of force-free electrodynamics, we aimed to fit the observed gamma-ray light-curves (peak properties) and tried to reasonably constrain the values of the angles α and ζ for MSPs with high signal-to-noise ratios. The best fits of ten pulsars were retrieved by minimising the square of the differences between the gamma-ray intensity observed and obtained from the model (2), integrated with all phase bins for each pulsar. We show and discuss our results below for each source, separately.
In our model, we assume radio emission heights to be close to the stellar surface of the simulated magnetosphere and the gamma-ray pulsations to start at r LC up to 2 r LC for the sake of simplicity. Radio emission height and the exact region where gamma-ray pulsations are produced are subject to a certain A&A 647, A101 (2021) Table 1. Estimated ranges of reasonable angles, α, ζ, by fits to the gamma-ray profiles with the artificial offsets, φ, for the best fits.

Source Name
Period (ms) Notes. χ 2 ν is the reduced chi-squared values for the best fits. (a) For PSR J1514-4946 the best fits were retrieved by accounting for gamma-ray intensities at phases 0.2-0.3 and 0.55-0.65, not the entire phase. amount of uncertainty, giving us some freedom to shift the gamma-ray light curves with respect to the radio peak, denoted by an additional phase shift φ accounting for a time-of-flight effect reflecting the variation in emission height of less than approximately 10% in the leading or trailing direction. For each pulsar, we set the offset φ to the value that gives the best fit to its gamma-ray light curve and plotted reduced chi-squared statistics for different angles of the line of sight and obliquity by adopting this φ. Although the least χ 2 ν can be found with any α − ζ, some angles cannot be excepted as good solutions. Possible solutions are restricted by the geometrical conditions: (1) α − W ≤ ζ ≤ α + W (we adopt W ∼ 20 • for all pulsars in this work) to be able to observe radio pulsation and (2) ζ α ≥ π/4 to see simultaneous radio and gamma-ray pulsations. We assume θ pc = 20 • for all MSPs investigated in this paper. Also, the retrieved gamma-ray light curves show symmetry in ζ around 90 • .

PSR J0030+0451
PSR J0030+0451 is a MSP with P = 4.87 ms which shows double pulses in gamma-rays and has a relatively clear light curve profile in gamma-ray energy bands. The primary radio peak of the source shows a multi-peaked structure rather than a smooth single peak. The interpulse observed with much smaller amplitude compared to the main pulse implies a geometry close to an orthogonal rotator while the proximity of α to 90 • depends on the details of the radio emission procedure. As can be seen in Fig. 3a, the best fit to the gamma-ray light curve of PSR J0030+0451 was obtained by fixing α = 70 • , ζ = 60 • and φ = −0.07 giving a reduced chi-squared of χ 2 ν = 7.2. We find that α = (55 • −75 • ) and ζ = (52 • −74 • ) gives fits with 7 χ 2 ν 12 (see Fig. 3b).

PSR J0102+4839
PSR J0102+4839 is another MSP with P = 2.96 ms which has a primary peak and a second peak with lower amplitude in gamma-rays. Its radio profile has a broad single peak featured with broader trailing edge compared to the leading edge of the peak. We obtained the best fit to the gamma-ray profile with α = 55 • , ζ = 70 • , φ = 0.03 and χ 2 ν 2.5; see Fig. 3c. By taking φ = 0.11, we calculated χ 2 ν for different α and ζ values which indicates that reasonable fits can be obtained with 2.5 χ 2 ν 4 for α = (50 • −70 • ) and ζ = (50 • −76 • ) (see Fig. 3b).

PSR J0614-3329
PSR J0614-3329, with P = 3.15 ms, has clear double peaks in the gamma-ray light curve. The amplitude of the leading peak is slightly lower than the trailing peak. The best fit to the gammarays with χ 2 ν = 10.4 is obtained with the parameters α = 75 • , ζ = 56 • and φ = 0.05; see Fig. 3g. The fits with 10 χ 2 ν 15 implies α = (70 • −75 • ) and ζ = (60 • −70 • ) (see Fig. 3h). In spite of high χ 2 ν which is mostly due to the off-pulse emission, there is a strong match between the peaks produced by the model and the observations which is what we are looking for as our main goal.

PSR J1124-3653
The gamma-ray profile of PSR J1124-3653 (P = 2.41 ms) can be interpreted as either (1) having a single peak feature with a trailing edge gradually decreasing in intensity or (2) having a less prominent second peak through the end of the trailing edge. We obtained the best fit to the gamma-ray profile with α = 65 • , ζ = 46 • and χ 2 ν = 2.4 (see Fig. 4a) while reasonable fits can be obtained with 2.4 χ 2 ν 4 for α = (50 • −65 • ) and ζ = (46 • −66 • ); see Fig. 4b. We do not need to shift the model curve along the phase to be able to find a plausible fit. We note that the amplitude of the second peak is overestimated in the best fit. 6.6. PSR J1514-4946 PSR J1514-4946 (P = 3.59 ms) shows double peaks in gammarays with bridge emission in between these peaks. Our model does not address this inter-pulse emission. Trying to find best fits for the entire light curve at any phase did not work for this pulsar because of bridge emission. We therefore focused on the two main peaks. The best fit in Fig. 4c was retrieved by accounting for the gamma-ray intensities at phases 0.2−0.3 and 0.55−0.65 which yielded a best fit with α = 55 • , ζ = 62 • , χ 2 ν 1 and φ = −0.05 (Fig. 4c). We found that reasonable light curves with 1 χ 2 ν 2 can be retrieved for α = (45 • −65 • ) and ζ = (46 • −66 • ) (Fig. 4d).

PSR J1614-2230
The best fit to the gamma-ray light curve of PSR J1614-2230 (P = 3.15 ms) is shown in Fig. 4e (Fig. 4f). However, solutions with α ∼ 90 • cannot represent the true geometry of the system because interpulse in radio was not observed for PSR J1614-2230.

PSR J2017+0603
PSR J2017+0603 is a MSP with P = 2.90 ms which has a similar gamma-ray profile to that of PSR J1514-4946, but with a narrower peak separation of ∆ ∼ 0.25. The best fit is obtained  (Fig. 4h).

PSR J2043+1711
Of the sources examined in this study, PSR J2043+1711 is the fastest spinning pulsar with P = 2.38 ms. The amplitudes and locations of both pulsations in the gamma-ray are well fitted with our model. We obtained the best fit by setting α = 55 • , ζ = 74 • and φ = −0.01 with χ 2 ν = 5.3 while the reasonable model curves could be found by the angles in the range α = (50 • −75 • ) and ζ = (46 • −80 • ) for 5.3 χ 2 ν 8; see Fig. 5b. We want to point out that we do not show the observed radio light curve of the source because we have not found any radio data in the 2PC files corresponding to this pulsar. However, it is not directly related to our analysis here as we only deal with the time-alignment of gamma-ray and radio not the exact shape of the radio profile.

PSR J2302+4442
PSR J2302+4442, with a period of P = 5.19 ms, shows double peaks in the gamma-rays. The best fit is produced with α = 50 • , ζ = 60 • and an offset φ = −0.03 with χ 2 ν = 11.1; see Fig. 5c. The ranges of the angles, provided that 11.1 χ 2 ν 17, are found to be α = (45 • −60 • ) and ζ = (48 • −62 • ); see Fig. 5d. In the best fit curve, the second peak perfectly matches the profile of the observed second peak. Although we overestimate the amplitude of the first peak, the peak locations are in line with the observations.

Discussion
In this study, we focused on the relationship between the radio time lag, the gamma-ray peak separation, the ratio of the peak amplitude to obtain the geometric parameters of pulsar magnetospheres, the obliquity, and the inclination of the line of sight. However, there is some uncertainty as to the precise location of the radio and high-energy emission sites. This leads to possible variations in the time lag between radio and gamma-rays. For instance, retardation effects within the light cylinder due to the finite propagation speed of light amount to an additional phase shift of φ ret = ∆r 2 π r L 1 2 π ≈ 0.16.
This shift has been implemented in our approach by adding a phase φ that also corrects for the period of MSP being longer than the 1.2 ms used for the force-free simulations and the emission sky maps. An additional unknown comes from the observation side. Indeed, there is no hint for the maximum of the radio peak to represent the centre of the pulse profile. It could well be that one or several cone emission patterns produce the pulses with sub-pulses of different intensities. In our simulations, we exactly know the location of the magnetic axis and the associated polar cap centre. Therefore, we used the radio peak as phase zero for our simulations. Nevertheless, for the observational data, there is no hint of a one-to-one correspondence between peak intensity and the middle of the pulse profile. This effect adds another unconstrained phase-shift between radio and gamma-rays. Further, while estimating the radio time lag δ, we assumed the radio emission to originate from regions of open dipolar magnetic field lines. However, estimates of δ will be affected if the radio emission arises from regions where there is an influence of the strong non-dipolar magnetic field.
A&A 647, A101 (2021) From a more fundamental physical perspective, the microphysics of pair creation, acceleration, and their outflow into the striped wind is still poorly understood. The particle density number, its energy distribution function, and their radiation spectra are not known with any accuracy. Nevertheless, the wind Lorentz factor, the current sheet thickness, and the flow velocity pattern all affect the gamma-ray pulse profiles. It is difficult to estimate the emissivity of the wind without the full description of the particle and radiation dynamics from the stellar surface to the wind. Clearly, further physical insight is required in order to precisely fit the pulse shape for individual pulsars. Such attempts have been avoided in this study because of the number of uncertainties in the model. Our approach, although simple, deals with the least number of parameters but already satisfactorily fits existing light curves. There is no doubt that more extensive works, including a kinetic description of the magnetosphere and its associated radiation mechanisms, will unveil many important fundamental issues in pulsar physics.

Conclusion
We numerically solved the equations for force-free pulsar magnetospheres for many obliquities α and computed the corresponding gamma-ray pulsed emission emanating from the striped wind. Through individual analyses of ten MSPs with spin periods in the range P ∼ (2−6) ms, we showed that their gammaray pulse profiles could be faithfully reproduced by assuming that radio pulses escape from open field-line regions close to the polar caps. In contrast, gamma-ray pulses are produced in the current sheet of the striped wind, outside the light cylinder. In this study, we constrained the range of parameters describing the geometry of each pulsar, namely their magnetic and light of sight inclination angles with respect to the spin axis of the star.
In order to retrieve good fits to the observed gamma-ray light curves of each source, we applied the least-square method to compare the reported intensities with those found from the model in each bin of the observation throughout one complete rotation of the star. We present our estimations together with their statistics in Table 1. We were able to fit all pulsars correctly.
The up-to-date measurement of radio pulse profiles and PPA traverse of MSPs indicates extremely irregular behaviour. Therefore, we were not able to use such techniques to deduce independent knowledge of their geometry. The distortion of the PPA can result from altitude-dependent aberration retardation effects (Mitra & Seiradakis 2004) or the presence of non-dipolar surface magnetic fields, which renders the PPA analysis useless for constraining α and ζ.
We plan to apply the same technique to young pulsars which possess much more well-defined radio pulse profiles and PPA. Indeed, thanks to the RVM model, the radio emission height has been located at about 5% of the light-cylinder radius (Mitra 2017). Even though such heights are much smaller than r L , they are large enough compared to the stellar radius, thus the dipole field is dominant. The force-free dipolar magnetosphere is therefore a good approximation for the computation of the observational signature in radio and gamma-rays. There is no more freedom for the radio altitude, and because of the narrow radio peaks, the uncertainties on the radio zero phase are no longer critical for fitting the δ − ∆ relation with high confidence. However, as young pulsars do not fall into the same category as MSPs, we will show the results in a separate forthcoming work.
The most recent and best-quality published gamma-ray data of pulsars were reported in the second Fermi Catalogue in 2013. Since then, almost one decade has passed and more observations have accumulated, increasing the signal-to-noise ratio of gamma-ray light curves and with improved intensity resolution. These updates are expected to be published as the third Fermi Pulsar Catalogue (3PC) in the future, which will undoubtedly increase the number of gamma-ray pulsars and provide even more stringent constraints on gamma-ray data to be tested with our model.