Efficiency of Nonthermal Pulsed Emission from Eight MeV Pulsars

We report on the properties of pulsed X-ray emission from eight MeV pulsars using XMM-Newton, NICER, NuSTAR, and HXMT data. For five of the eight MeV pulsars, the X-ray spectra can be fit by a broken power-law model with a break energy of ∼5–10 keV. The photon indices below and above the break energy are ∼1 and ∼1.5, respectively. In comparison with the X-ray emission of the Fermi-LAT pulsars, the MeV pulsars have a harder spectrum and a higher radiation efficiency in the 0.3–10 keV energy bands. When isotropic emission is assumed, the emission efficiency in the keV–MeV bands is estimated to be η MeV ∼ 0.01–0.1, and this is similar to the efficiency of the GeV emission of the Fermi-LAT pulsars with a similar spin-down power. To explain the observed efficiency of the MeV pulsars, we estimate the required pair multiplicity as 104–7, which depends on the emission process (curvature radiation or synchrotron radiation) and on the location in the magnetosphere. The high multiplicity indicates that the secondary pairs that are created by a pair-creation process of the GeV photons produce the X-ray/soft gamma-ray emission of the MeV pulsars. We speculate that the difference between MeV pulsars and Fermi-LAT pulsars can be attributed to the difference in viewing angle measured from the spin axis if the emission originates from a region inside the light cylinder (canonical gap model) or to the difference in the inclination angle of the magnetic axis if the emission is produced in the equatorial current sheet outside the light cylinder.


INTRODUCTION
Since the launch of F ermi-Large Area Telescope (hereafter F ermi-LAT) in 2008, more than 300 gammaray emitting rotation-powered pulsars have been identified (Smith et al. 2023).The F ermi-LAT observations have revealed that the typical spectral energy distribution of the energetic pulsars shows a peak in the GeV bands.Moreover, the GeV/TeV observations also show that the observed pulsed gamma-rays from the pulsars are produced in the outer magnetosphere, because the gamma-ray photons produced near the stellar surface cannot avoid the pair-creation process with the strong magnetic field (Aliu et al. 2008;MAGIC Collaboration et al. 2020).The discovery of ∼ 20 TeV pulsed emission from the Vela pulsar also supports this hypothesis (H.E. S. S. Collaboration et al. 2023).Either takata@hust.edu.cncanonical gap model (Cheng et al. 1986;Arons 1981) or the equatorial current sheet model (Lyubarskii 1996;Cerutti et al. 2015;Chernoglazov et al. 2023) has been considered as the possible acceleration region in the outer magnetosphere.
The so-called MeV pulsars form a distinct group from the F ermi-LAT pulsars, based on their observed emission properties (Wang et al. 2014;Kuiper & Hermsen 2015;Harding & Kalapotharakos 2017;Kuiper et al. 2018).PSR B1509-58 located in the supernova remnant MSH 15-52 has been known as an energetic young pulsar since the high-energy observation in the Energetic Gamma-Ray Experiment Telescope era (Thompson et al. 1999).The F ermi-LAT observation confirms the pulsed emission from PSR B1509-58 in 0.1-1 GeV bands, but the spectrum from X-ray to GeV bands has a peaking at ∼ 1 − 10 MeV, which is about three orders of magnitude smaller than that of the typical F ermi-LAT pulsars (Kuiper et al. 1999;Abdo et al. 2010).Kuiper et al. (2018) report the F ermi-LAT detection of PSR J1846-0258 and reveal that the spectral shape in the X-ray to GeV bands is similar to that of PSR B1509-58.Based on the spectral properties reported in Wang et al. (2014) and Kuiper & Hermsen (2015), we select the eight sources as the MeV pulsars.Table 1 summaries the spin-period (P s ), dipole magnetic field (B s ), the spin-down power (L sd ), and the distance (d) of the eight pulsars.
It has been revealed that the X-ray emission of the MeV pulsars has a distinct property from that of the F ermi-LAT pulsars.The X-ray emission of the F ermi-LAT pulsars is in general described by the thermal component from the neutron star surface and/or the non-thermal component that has a photon index of Γ 1 ∼ 1.5 − 2 (Abdo et al. 2013;Hsiang & Chang 2021).It has been discussed that the non-thermal Xray emission is produced by the synchrotron radiation from the secondary electron/positron pairs that are produced by the pair-creation process of the GeV gammarays (Kisaka & Tanaka 2017).For the MeV pulsars, on the other hand, the previous X-ray studies have collected evidences that the broadband X-ray spectrum is described by a single power law with a photon index of Γ 1 < 1.5 or a power law with a turnover in 1 − 10 keV bands [see e.g.Chen et al. (2016), Hu et al. (2017) and references therein for PSR B1509-58, Kargaltsev et al. (2009) for PSR J1617-5055, Madsen et al. (2020) for PSR J1811-1925, Kuiper & Hermsen (2015) for PSR J1813-1749, Lin et al. (2009) and references therein for PSR J1838-0655, Gotthelf et al. (2021) for PSR J1846-0258, Gotthelf et al. (2011) and Kim et al. (2024) for PSR J1849-0001, and Lu et al. (2007) for PSR J1930+1852].We also refer Kuiper & Hermsen (2015) for a comprehensive observational study of the hard X-ray/soft gamma-ray emissions of the young pulsars.
Although several emission models have been proposed (Wang et al. 2013(Wang et al. , 2014;;Harding & Kalapotharakos 2017), the emission process and even the emission region of the MeV pulsars are still not clearly understood.
In this paper, we revisit the X-ray analysis for the eight MeV pulsars using the data taken by XMM-Newton telescope, the Neutron Star Interior Composition Explorer (NICER), Nuclear Spectroscopic Telescope Array (NuSTAR), and Hard X-ray Modulation Telescope (HXMT).Main purpose of this study is (i) to estimate the efficiency of the non-thermal emission of the MeV pulsars (η MeV ), which is the ratio of the luminosity in X-ray/soft gamma-ray bands to the spin-down power, and (ii) to compare it with the efficiency of the "GeV emission" (η GeV ) of the F ermi-LAT pulsars.The emission efficiency can be used to constrain the particle acceleration and emission processes.For the F ermi- Values of the parameters are obtained from ATNF Pulsar Catalog (Manchester et al. 2005).
LAT pulsars, for example, the efficiency of the GeV emission is correlated with the spin-down power (L sd ) as , and the origin of the correlation has been investigated with the emission model (Takata et al. 2010;Kalapotharakos et al. 2019).The emission efficiency can be linked with so-called multiplicity (κ) that is the ratio of the number density of the emitting electrons/positrons over the Goldreich-Julian value, n GJ = Ω s B/(2πce), where Ω s is the spin angular frequency, B is the local magnetic field strength, c is the speed of light and e is the elementary charge (Goldreich & Julian 1969).The required multiplicity to explain the observed efficiency can be used to investigate the emission process, too.
The paper is organized as follows.In section 2, we describe the data reduction of each telescope.In section 3, we estimate the efficiency of the non-thermal emission of the MeV pulsars by fitting the observed spectra.In section 4, we compare the efficiency of the MeV pulsars with that of the F ermi-LAT pulsars and we discuss the required multiplicity for different emission processes.Then we summarize our results in section 5.In this paper, we apply Gaussian-cgs units system to express the physical quantities.

XMM-Newton
Because the MeV pulsars are usually located in their pulsar wind nebulae, it is required to carry out a spinphase resolved spectroscopy to extract the pulsed spectrum.For the XMM-Newton, therefore, we analyze the data taken by PN-CCD (hereafter PN), because the timing resolution data taken by MOS-CCD is insufficient for the MeV pulsars.We use the SAS task4 epproc to obtain the cleaned event file and then we apply the task barycen to carry out barycentric correction for the arrival time of the each event (ephemeris DE 405).We remove epochs with background flaring by limiting the count rate above 10 keV to < 0.4 count s −1 .
For each event file analyzed in this study, we create Z 2 1 -periodogram (Buccheri et al. 1983) and identify the periodic signal at the spin-frequency of the pulsar.We fold each event file using the best frequency in each periodogram, and create the light curve to determine the onpulse and off-pulse phases.For all pulsars in this study, the folded light curve is described by a single broad peak, such as the pulse profile of PSR B1509-58 in Figure 1 (Figure B1 for eight MeV pulsars).We therefore define the on-pulse phase covering 60% of the single spin phase and the off-pulse phase with other 40%, as indicated in Figures 1 and B1.We create the spectra of the on-pulse phase and the off-pulse phase, and obtained the spectrum of the pulsed component by subtracting the spectrum of the off-pulse phase from the spectrum of the on-pulse phase.For XMM-Newton data, we group the source spectrum such that each spectral bin has a minimum signal-to-noise ration of 3. We use the standard tasks of the SAS to generate the spectral file, response matrix file (rmf) and auxiliary response file (arf).We confirm with SAS command epatplot5 that the effect of the pile-up on the spectrum is negligible for the eight MeV pulsars.

NICER
NICER observed our targets except for PSRs J1617-5055 and J1930+1852.We apply nicerl2 to each data set to obtain the cleaned event file and we apply a barycentric time correction to the created event file using the barycorr task of HEASoft.NICER repeatedly observes one target and the exposure of each observation is usually several kilo-seconds.For each pulsar, we create a local ephemeris by combining several data sets; all ephemerides and obs.IDs of the NICER data used in this study are summarized in Table A3.We also refer the ephemerides reported in Ho et al. (2020) for J1831-1749 and Hu et al. (2023) for J1846-0258 We merge the multiple observations using the task niobsmerge6 .We create the good time interval files for on-pulse and off-pulse phases, and extract the events using extractor.We use niextspect to extract the spectrum from the event files, and create the rmf and arf from the merged event file using the task nicerl3-spect.We group the events so that each spectral bin contains a count of > 1000, since each channel before binning records a few hundred events.We checked the different threshold of the photon counts for the binning and confirmed that results of the fitting is insensitive to how we bin the data.

NuSTAR
First, we carry out the tasks nupipeline and nuproducts of HEASoft to obtain the clean event list after barycentric time correction.In Z 2 1 -periodogram, we can confirm the strong spin signal in the data extracted from the source region.We fold each data set with the best frequency in Z 2 1 -periodogram and create the good time interval files for on-pulse phase and off-pulse phase.We rerun to the task nuproducts with the option usrgtifile and obtain the spectra of the on-pulse/offpulse phases and response files (rmf and arf).Using the pipeline nuproducts, we group the on-pulsed spectra to ensure at least 100 counts per spectral bin, since each channel before binning records a few ten events.

HXMT
HXMT has three payloads [high-energy (20-250 keV), medium-energy (8-35 keV) and low-energy (1-12 keV) X-ray telescopes] and can cover a wide energy range in X-ray bands.In this study, we use the data of the highenergy X-ray telescope of PSR B1509-58 taken in 2017 July (Obs.ID.P0101324001); we cannot find a significant pulsation in the data taken by medium-and lowenergy telescopes.We apply the task hepical to remove spike events caused by electronic systems and hegtigen to select good time interval file with the conditions that (i) the pointing offset angle is smaller than < 0.04 the pointing direction is > 1 • above Earth, (iii) the geomagnetic cut-off rigidity is bigger than > 8 GeV, and (iv) the South Atlantic Anomaly does not occur within 300 s.We use hxbary for the barycentric correction of the HXMT data, and we confirm a significant periodic signal in the data of the high-energy X-ray telescopes (Figure 1).We use the hspec-merge to merge the spectra and the response files extracted from the different data sets.For HXMT data, we group the on-pulsed spectra to ensure at least 1000 counts per spectral bin.

RESULTS
In this section, we present the results of the spectral fitting for the pulsed components of the eight MeV pulsars.We carry out a joint fit of the spectra taken by the four telescopes.Using Xspec, we fit the spectra with a single power law model constant*tbabs*powerlaw (hereafter s-pl model) or a broken power law model constant*tbabs*bknpower (hereafter bk-pl model), where constant accounts for cross-calibration mismatch among the telescopes.
In previous studies (Cusumano et al. 2001;Chen et al. 2016;Hu et al. 2017), a log-parabola (logpar) model, which is described by a function of F (E) ∝ E −[a+blog(E/E0)] with a and b being fitting parameters and E 0 the energy scale, is used to describe the curved X-ray spectrum of the MeV pulsars.In this study, the bk-pl model is preferred over the logpar model, because it is more intuitive to physically interpret.As we discussed in the following section 4.2.2, for example, we will use the fitting break energy of bk-pl model to constrain the multiplicity of the pairs.
In our analysis, we fix the constant factor for the XMM-Newton observation at unity.Figures 2 and B2-B4 show the pulsed spectra for the eight pulsars.Although we choose only one spectrum of each detector in the figure for the representation, we use more available data (as shown in Table A1-A3) for the spectral fitting and fit all spectra simultaneously; the Figures 2 and B2-B4 represent the spectra taken by the earlier observation if there are multiple observations.Table 2 summarizes the results of the fitting spectral parameter, the estimated isotropic luminosity, and its emission efficiency in 0.3-150 keV bands.

PSR B1509-58
For PSR B1509-58, there is one archival XMM-Newton data taken with the PN camera under the small widow mode, whose timing resolution (∼ 5.7 ms7 ) is sufficient enough to carry out the phase-resolved spectroscopy.There are four data sets obtained by NuSTAR observations.NICER observed the target in 2017 July and August, for which we create local ephemerides to fold the data.The left panel of Figure 2 shows the extracted spectra of the pulsed component.
First we fit the spectrum with a s-pl model and obtain the hydrogen column density of N H = 1.32(5) × 10 22 cm −2 and the photon index of Γ 1 = 1.41(1) with χ 2 /D.O.F.= 2023.99/2023(Table 2).We find that the predicted N H is higher than N H ∼ 0.9 × 10 22 cm −2 reported in the previous studies using Bepposax observation (Cusumano et al. 2001;Romani et al. 2023) (the previous studies used Phabs for the absorption model in Xspec, but we confirm that the fitting results are not sensitive to the absorption model and obtain a similar to current value).
Second, we fit the spectra with a bk-pl model.As Table 2 indicates, the improvement of the fitting from the s-pl model is significant with an F static value of 32 (p = 1.5 × 10 −14 ).The obtained hydrogen column density of N H = 1.0(2) × 10 22 cm −2 is more consistent with the previous value.Our fitting indicates a spectral break at E b = 5(2) keV, and the photon index below and above the break are Γ 1 = 1.1(2) and Γ 2 = 1.43(3), respectively.The photon index below the break energy is consistent with the previous Bepposax observation (Cusumano et al. 2001) and Chandra observation (Hu et al. 2017), in which a single power law model is applied to fit the observation data.The photon Histogram show the pulse profiles measured by XMM-Newton PN.Hydrogen column density is fixed to NH = 0.9 × 10 22 cm −2 for PSR B1509-58 and 6.0 × 10 22 cm −2 for PSR J1838-0655.The break energy at the wings of the pulse (∼ 0.25 and ∼ 0.75 spin-phase) of PSR J1838-0655 cannot be constrained.
Because of enough count rate, we investigate the evolution of the break energy (E b ) in the bk-pl model with the spin-phase.We divide the on-pulsed phase (Figure 1) into six segments, and we fit each spectrum with a bk-pl model; during the fitting, we fix the hydrogen column density to N H = 0.9×10 22 cm −2 .The left panel of Figure 3 shows the evolution of the break energy with the spin phase for PSR B1509-58.We can see that the break energy is almost constant with E b = 4 − 8 keV, although a small variation may exist.
We use cflux in Xspec to estimate the unabsorbed flux in 0.3-10 keV bands for each fitting model and calculate the isotropic luminosity, L X = 4πd 2 f X , where f X is the estimated flux and d is the distance obtained from ATNF Pulsar Catalog (Manchester et al. 2005).Figure 4 compares the radiation efficiency (≡ η 0.3−10 ) in 0.3-10 keV bands of the MeV pulsar wind those of the F ermi-LAT pulsar (see section 4.1.1 for a detailed discussion).To evaluate the radiation efficiency in the wide energy bands, we also estimated the efficiency in 0.3-150 keV using Xspec.We choose the maximum energy of 150 keV because there is no HXMT observations for other MeV pulsars and because Xspec extrapolates the flux up to 150 keV using the NuSTAR data and the assumed spectral shape.The estimated luminosity and its efficiency (≡ η X ) in 0.3-150 keV bands are summarized in Table 2.For PSR B1509-58, we obtain the efficiency of η X ∼ 0.03 − 0.04 in 0.3-150 keV bands.

PSR J1657-5055
XMM-Newton observed the target with the timingmode of the PN camera.The spin signal P s ∼ 69.4 ms is confirmed in the Z 2 1 -periodogram.Because of the onedimensional imaging capability of the data taken by the timing mode, we select the pixels so that the extracted events maximize the power of the periodic signal.There is one archival data of the NuSTAR observation, but no observation has been carried out by the NICER.The right panel of Figure 2 shows the extracted spectrum of the pulsed component.The spectrum is well fitted by a s-pl model with N H = 4.0(7) × 10 22 cm −2 and Γ 1 ∼ 1.52(5).We find that the best fitting hydrogen column density is slightly smaller than N H = 6.5(7) × 10 22 cm −2 reported in Hare et al. (2021), in which only NuSTAR data is used.The photon index Γ 1 = 1.52(5), on the other hand, is consistent with ∼ 1.51(2) of the previous study.As the right panel of Figure 2 indicates, the quality of the data is insufficient to constrain the parameters of the bk-pl model.
The spectrum of the NuSTAR data covers the energy range up to ∼ 70 keV, and no spectrum cut-off has been observed in 10-70 keV bands.By assuming that the spectrum above 70 keV extends with a single power law, we use Xspec to extrapolate the flux up to the energy 150 keV and obtain the efficiency of η X = 0.0050(5) in 0.3-150 keV bands.
In the spectrum taken by NuSTAR in Figure 2, one may notice an absorption feature at ∼ 50 keV, where the data points of both NuSTAR's FPMA and FPMB detectors drop from the average of the residual.We will also observe a similar feature in the spectrum of PSR J1849-0001 with a higher significance level (section 3.7).In this study, we apply the Bayesian method to evaluate the absorption-line feature and utilize a python tool, Bayesian X-ray Analysis (BXA, Buchner et al. 2014).By grouping the NuSTAR's spectra to ensure at least one count per spectral bin, a spl with Gaussian absorption model is compared with a pure s-pl model.For PSR J1617-5055, we find that the inclusion of the absorption feature does not improve the fitting results and could not give an evidence of the absorption feature.

PSR J1811-1925
For PSR J1811-1925, there are one archival data taken by the XMM-Newton PN camera under the small window mode and one archival data by NuSTAR.We also use the NICER observations carried out in 2022 March to September.We fit the pulsed spectrum (Figure B2 and Table 2) with a s-pl model and obtains N H = 4.0(8) × 10 22 cm −2 and Γ 1 = 1.40(8), which is consistent with the previous study using Chandra and NuS-TAR data (Madsen et al. 2020).Similar to PSR J1657-5055, the quality of the data is insufficient to constrain the parameters of the bk-pl model.The isotropic emission efficiency with the s-pl model is estimated to be η X = 0.013(2) in 0.3-150 keV bands.

PSR J1813-1749
We use two archival data taken by the XMM-Newton PN camera under the small window mode and one archival data of NuSTAR.We also analyze the data of the NICER observation in 2019 June-August, and refer the ephemeris provided by Ho et al. (2020).The observed spectrum (Figure B2) is fitted by a s-pl model with a large column density of N H = 13(4) × 10 22 cm −2 and photon index of Γ 1 = 1.8(3).Because of the strong absorption, the observed spectrum cannot constrain the parameters of the bk-pl model.In order to compare with previous results, we fix the hydrogen column density to N H = 9.8 × 10 22 cm −2 .The photon index of the s-pl model is Γ 1 = 1.5(2), which is consistent with ∼ 1.48 of the INTEGRAL observation reported in Kuiper & Hermsen (2015).The efficiency in 0.3-150 keV bands is estimated to be η X = 0.0008(4).

PSR J1838-0655
Both XMM-Newton and NuSTAR observed the source twice and the observations of the XMM-Newton PN were carried out under the small window mode.We also analyze the data of the NICER observations carried in 2019, 2020 and 2021.The observed spectrum (right panel in Figure B2) can be fitted by a s-pl model with hydrogen column density of N H = 6.9(2) × 10 22 cm −2 and photon index of Γ 1 = 1.40(2).We find in Table 2 that the bk-pl model significantly improves the fitting with an F static value of 24 (p = 5 × 10 −11 ).In the bk-pl model, we obtain a break energy of E b = 8.1(9)keV.The photon index below and above the energy break are Γ 1 = 1.04(9) and Γ 2 = 1.49(4), respectively, which are consistent with the values reported in Lin et al. (2009), who analyze the Chandra, RXTE and the SUZAKU data.The current hydrogen column density, N H ∼ 6.0 × 10 22 cm −2 , is larger than N H ∼ 4.5(3.7 − 5.2) × 10 22 cm −2 obtained with Chandra observation (Gotthelf & Halpern 2008), a Photon index for the s-pl fitting or index below the break energy for the bk-pl fitting.
b Break energy for the bk-pl fitting.
c Photon index above break energy for the bk-pl fitting.
d Luminosity in 0.3-150 keV energy bands, e Efficiency of the emission in 0.3-150 keV energy bands.
f Distance is assumed to be d = 7 kpc (Gotthelf et al. 2011).
g Fixed to the value reported in Temim et al. (2010).
We obtain the efficiency of η X ∼ 0.06 − 0.07 in 0.3-150 keV bands, which is the highest among our targets.We perform the phase-resolved spectroscopy using the XMM-Newton and NuSTAR data.We extract the spectrum of 0.1 phase interval in the on-pulse phase and fit each spectrum with the bk-pl model by fixing the hydrogen column density to N H = 6.0 × 10 22 cm −2 .The right panel of Figure 3 shows the evolution of the break energy with the spin phase.The break energy (E b ) in the wing of the pulsed shape (∼ 0.25 and ∼ 0.75 spin phase) cannot be constrained as the figure shows.With the large uncertainty, we cannot discuss the evolution of the break energy with the spin-phase.However, the break energy may be smaller at ∼ 0.6 − 0.7 spin-phase of the right tip on the peak.

PSR J1846-0258
PSR J1846-0258 is known as a high-magnetic field pulsar and it experienced a magnetar-like X-ray outburst in 2006 (Gavriil et al. 2008) and in 2020 (Blumer et al. 2021b;Hu et al. 2023).The X-ray emission after the outburst is fitted by the thermal emission from the surface plus single power law tail in hard X-ray bands, and it shows the temporal evolution (Hu et al. 2023).Since we are interested in the emission in a quiescent state, we analyze the data taken before 2020 outburst.We use the data of the XMM-Newton observation (small window mode for PN) and NuSTAR observation carried out in 2017 September.We also generate the spectrum with the NICER data taken in 2018, 2019 and 2020, for which we refer to the ephemerides reported in Hu et al. (2023).

PSR J1849-0001
We use one archival data taken by the XMM-Newton PN camera under the small-window mode and one data taken by NuSTAR.For the NICER observation, we construct a local ephemeris for the 2018 data.Although Bogdanov et al. (2019) provide the ephemeris covering 2018 February to September, we extend the ephemeris to cover 2018 February to November.The right panel of Figure B4 shows the extracted spectrum of PSR J1849-0001.
The s-pl model provides an acceptable fitting result (Table 2), and its parameters are N H = 5.1(3)×10 22 cm 2 and Γ 1 = 1.48(5).Although the obtained N H is smaller than N H = 8.1(2) × 10 22 cm 2 reported in Kim et al. (2024), who use the same XMM-Newton and NuSTAR data as our study, the photon index is consistent with their Γ 1 = 1.42(3) within the error.The bk-pl model improves the fitting results with an F -static value of 8.2 (p = 0.00031), and the photon index (Γ 1 = 0.8(4)) below the break energy (E b = 4.9(8) keV) is consistent with ∼ 1.1 in the previous studies (Gotthelf et al. 2011;Kuiper & Hermsen 2015), in which XMM-Newton data is fitted by a s-pw model.The obtained N H = 4.0(9) × 10 22 cm 2 is also consistent with N H = 4.3(6) × 10 22 cm 2 of the previous study.Assuming the distance 7 kpc to the source (Gotthelf et al. 2011), the efficiency in 0.3-150 keV bands is estimated to be η X ∼ 0.025.Kim et al. (2024) suggest the existence of a peaking of the spectrum at ∼ 60 keV, which is lower than ∼ 1 − 5 MeV of PSRs B1509-58 and J1846-0258.As the left panel in Figure B4 shows, on the other hand, our extracted spectrum may indicate the feature of an absorption at around 40 keV rather than the peaking of the spectrum.In our study, we group the spectrum of the NuSTAR data at least 100 counts per spectral bin.The absorption-like feature in the spectrum becomes more obvious if we group the spectrum with smaller photon counts per spectral bin.Similar to the descriptions in section 3.2, we evaluate the absorption feature with BXA.We fit the NuSTAR spectra with a s-pl plus Gaussian absorption model, and obtain the energy (E ab ) and width (σ ab ) of the Gaussian absorption line as E ab = 36 +3 −4 keV and σ ab = 0.6 +1 −0.5 keV, respectively.We compare the s-pl with a Gaussian absorption line model with a pure s-pl model and obtain a Bayesian factor of ∼ 2.5, which is not large enough to claim the strong evidence of the absorption-line feature; a Bayesian factor of >10 is usually desired to confirm the absorption feature (Buchner et al. 2014).If the absorption future were real, the non-thermal emissions from PSR J1849-0001 would be likely produced in the region located near the neutron star surface.This is because such an absorption feature is possibly associated with a cyclotron absorption in a highly magnetized and dense region around the neutron star.Since the current data is not enough to draw any conclusion, it would be desired to investigate the absorption as well as peaking feature in the spectrum by a deeper observation.

PSR J1930+1852
XMM-Newton observed the target once under the large window mode and once under the full frame mode, for which the timing resolutions are still sufficient enough to carry out the phase-resolved spectroscopy for PSR J1930+1852.NuSTAR observed the target once, but NICER has not observed the target.As Table 2 shows, the observed spectrum (Figure B4) is fitted by a s-pl model with a photon index of Γ 1 ∼ 1.53(8).For the bk-pl model, the current data cannot constrain the hydrogen column density, and therefore we fix the value to N H = 1.95 × 10 22 cm −2 determined by the Chandra observation (Temim et al. 2010).We obtain the parameters of Γ 1 = 0.8( 4), E b = 5(2) keV and Γ 2 = 1.5(1).With the current large uncertainty, the photon index below the break energy is consistent with Γ 1 = 1.2 − 1.3 obtained by the Chandra observation (Lu et al. 2002;Temim et al. 2010).The estimated efficiency in 0.3-150 keV bands is η X = 0.007 − 0.008.Figure 4 compares the observed efficiency of the isotropic emission in 0.3-10 keV energy bands (η 0.3−10 ) of the MeV pulsars with those of the F ermi-LAT pulsars, for which we refer the X-ray flux of the pulsars reported in the F ermi-LAT second pulsar catalog (Abdo et al. 2013); in our study, PSR B1509-58 is treated as a MeV pulsar, although it is listed in the LAT pulsar catalog.Since the F ermi-LAT catalog reports the phase-averaged value, we scale the pulsed flux of the MeV pulsars to the phase-average value by multiplying 0.6.
From Figure 4, we can see the tendency that currently known MeV pulsars are distributed in a narrower range of the spin-down power than the F ermi-LAT pulsars, and no MeV pulsars with L sd < 2 × 10 36 erg s −1 is found.We find that the MeV pulsars have a larger efficiency than the averaged value (∼ 10 −4 ) of the F ermi-LAT pulsars.There are six F ermi-LAT pulsars (Crab, PSRs J1119-6127, J1741-2054, J1747-2958, J1801-2451 and J2021+3651), who efficiency is greater than η 0.3−10 > 5 × 10 −4 and has an efficiency similar to  the MeV pulsars.We investigate the spectral properties of those F ermi-LAT pulsars in the literatures.Among the F ermi-LAT pulsars, the Crab pulsar (L sd = 4.36 × 10 38 erg s −1 ) has a unique broadband spectrum of the non-thermal emission.The spectrum has a photon index of ∼ 1.5 in 0.1-10 keV bands, ∼ 2 in 10-100 keV bands and > 2 in > 100 keV bands (Kuiper et al. 2001).The X-ray and soft γ-ray emissions of the Crab pulsars have been interpreted as a result of the synchrotron emission of the secondary pairs around the light cylinder (Takata & Chang 2007;Harding & Kalapotharakos 2015).The spectral evolution from X-ray to MeV bands of the Crab pulsar is different from those of the MeV pulsars (e.g.PSRs B1509-58.J1838-0655 and J1849-0001), in which the photon index below E b ∼ 5 − 10 keV is Γ 1 ∼ 1 and ∼ 1.5 above E b (Table 2).Moreover, F ermi-LAT confirmed that the spectra of the two MeV pulsars, PSRs B1509-58 and J1846-0258, do not extend in the GeV energy band.These differences may suggest that the emission process and/or emission region of the MeV pulsars is different from that of the Crab pulsar.
For four F ermi-LAT pulsars (PSRs J1119-6127, J1741-2054, J1801-2451 and J2021+3651) with η 0.3−10 > 5 × 10 −4 in 0.3-10 keV bands, the nonthermal X-ray spectrum is described with a photon index of ∼ 1.6 − 2.0 (Van Etten et al. 2008;Ng et al. 2012;Karpova et al. 2014;Marelli et al. 2014;Coti Zelati et al. 2020;Blumer et al. 2021a;Woo et al. 2023), which is greater than ∼ 1 − 1.5 of MeV pulsars.It may be worth to note that PSRs J1119-6127 of the F ermi-LAT pulsar and J1846-0258 of the MeV pulsar have been known as high-magnetic field pulsars, and they have experienced the magnetar-like Xray outbursts (Gavriil et al. 2008;Gögüş et al. 2016;Archibald et al. 2017;Wang et al. 2020;Blumer et al. 2021b).Their spin-down properties are also resemble to each other; PSR J1119-6127 has P s ∼ 407 ms, B s ∼ 4.1 × 10 13 G and L sd ∼ 2.3 × 10 36 erg s −1 .However, the X-ray emission properties in a quiescent state are different.The X-ray emission of PSR J1119-6127 is dominated by thermal component, of which the flux in 0.5-8 keV bands is a several factor more than that of the non-thermal component of Γ 1 ∼ 2.1 (Ng et al. 2012).For PSR J1846-0258, on the other hand, the observed X-ray spectrum is described by a non-thermal emission, as we described in 3.6.The investigations of PSRs J1119-6127 and J1846-0258 will suggest that the pulsar spin-down parameter is not key parameters to differentiate between the F ermi-LAT pulsars and the MeV pulsars.
Finally, PSR J1747-2958 is the F ermi-LAT pulsar having the highest X-ray efficiency η 0.3−10 ∼ 5 × 10 −3 in Figure 4, and its pulsed X-ray emission is reported by Li et al. (2018).The X-ray emission from the pulsar region taken by the Chandra observation can be described by a non-thermal emission with a photon index of 1.25 − 1.55 (Klingler et al. 2018;Li et al. 2018).The relatively hard X-ray spectrum with a high efficiency of PSR J1747-2958 is similar to the X-ray properties of the MeV pulsars.However, because PSR J1747-2958 is associated with bright pulsar wind nebula (so called Mouse), it cannot exclude the possibility that a phase-averaged spectrum reported in the previous studies is dominated by the nebula emission (Klingler et al. 2018); the phaseresolved spectroscopy is not carried out in the previous studies.We check the archival XMM-Newton PN data and NuSTAR data, but we could not confirm the pulsed emission in the periodogram.A deeper observation will be required to understand the X-ray emission properties of PSR J1747-2958.

Efficiency of total non-thermal pulsed emission
Multi-wavelength observations have revealed that the non-thermal pulsed emissions from the MeV pulsars and F ermi-LAT pulsars have in general the most radiation efficiency in 1 − 10 MeV and 1-10 GeV bands, respectively.For PSR B1509-58 of the MeV pulsar, the observed flux in 100 keV-10 MeV bands (Kuiper et al. 1999) indicates that the isotropic efficiency increases by about factor of two to three from the value (η X ) reported in Table 2.By assuming that other MeV pulsars have a similar spectral behavior to PSR B1509-58, the typical efficiency of the MeV pulsar would be estimated to be η MeV ∼ 2 − 3 × η X ∼ 0.01 − 0.1; the isotropic efficiency of PSR J1813-1749 could be still smaller than 0.01.Since the efficiency (η X ) in the 0.3-150 keV bands in Table 2 represents, in order of magnitude, the total efficiency of the MeV pulsars, we compare the efficiencies of the eight MeV pulsars with the efficiencies in 0.1-300 GeV bands of the F ermi-LAT pulsars in Figure 5.We plot the efficiency mainly referred to the bk-pl model.If the spectra cannot be fitted by the bk-pl, we then derive the value from the s-pl model instead.We note that few F ermi-LAT pulsars appears with an efficiency greater than unity.This is likely due to over-estimated distance and/or the assumption of the isotropic radiation (Abdo et al. 2013).
Based on Figure 5, we may conclude that the efficiency of the total non-thermal pulsed emission of the MeV pulsars is similar to those of the F ermi-LAT pulsars having similar spin-down powers.Even though the total efficiency of the MeV pulsar could increase by a factor of two to three from the plotted values, it is still consistent with that of the F ermi-LAT pulsars.This similarity of the efficiency may suggest that the emission of the MeV pulsars is either (i) a primary emission but the maximum energy of the electrons/positrons is smaller than that of F ermi-LAT pulsars or (ii) a sec-ondary emission from the electron/positron pairs, into which most of the primary GeV emission are converted.

Emission process
In the last section, we estimate that if the emission of the MeV pulsars are isotropic, the typical efficiency is η MeV ∼ 10 −2 −10 −1 .The actual radiation geometry will depend on various factors such as the emission process, inclination angle and the size of the emission region, and it is not trivial to determine the true efficiency from the observed flux.Since the observed pulsed width covers ∼ 60 % of the one rotation, the solid angle of the radiation may not be much smaller than the isotropic value of 4π.Hence we assume the efficiency η MeV ∼ 10 −2 − 10 −3 as a more conservative value.The actual efficiency could take an even smaller value, if the line of sight, the axis of the spin and axis of the radiation beam were almost aligned in the same direction.In this paper, however, we do not consider such a peculiar situation.
In the followings, we discuss either the possibility of the curvature radiation process or the synchrotron radiation process.The inverse-Compton process (IC process) can also produce the high-energy photons in the pulsar magnetosphere.For the young pulsars, however, it has been predicted that the energy of the scattered photons are in the TeV energy bands, which can explain the recent observation of the Vela pulsar (H.E. S. S. Collaboration et al. 2023).Ng et al. (2014) study non-thermal X-ray emission from the millisecond pulsars whose X-ray pulse is in phase with the gamma-ray pulse.They suggest that the IC process of the accelerated electrons (or positrons) off the radio wave produces the non-thermal X-ray emission.However, this model predicts an accompany of the GeV emission and is incompatible with the non-detection of the GeV emission of the MeV pulsars.We, therefore, will not discuss the IC process as the origin of the MeV pulsar's emission in the following discussion.

Curvature radiation process
It is widely accepted that the GeV emission is produced in the region located near or beyond the light cylinder, because the GeV photons produced near the stellar surface cannot avoid the pair-creation process with the magnetic field.For canonical "gap" models of the pulsar magnetosphere (eg.Cheng et al. 1986;Arons 1981), the GeV emission is produced by the curvature radiation process, due to the bending of the magnetic field lines, of the primary electrons and/or positrons accelerated by the electric field along the magnetic field line.For this model, the number density in the acceleration region cannot be much bigger than the Goldreich-Julian value, n GJ .The characteristic energy and radiation power of the curvature radiation process from single electron (or positron) are provided by E curv = 3hcγ 3 /(4πR c ) and P curv = 2e 2 cγ 4 /3R 2 c , where γ is the Lorentz factor of the particles, h is the Planck constant and R c is the curvature radius of the magnetic field line.
One can show that when the dynamical timescale τ dyn ∼ R c /c is equal to the cooling timescale of the curvature radiation process, τ curv ∼ 3m e cR 2 c /(2e 2 γ 3 ), the characteristic energy of the radiation is E curv ∼ 9m e c 2 /(4α) ∼ 158 MeV, where m e the mass of an electron and α is the fine structure constant.Hence, if the curvature radiation is the origin of the non-thermal emission in X-ray to soft γ-ray bands of the MeV pulsars, it will be taken place in a slow cooling regime.
In the slow cooling regime, the fraction of the energy of the particles that is released as the curvature radiation process during the dynamical timescale will be τ curv /τ dyn ∼ E curv /158 MeV, which is ∼ 0.01 for E curv ∼ 1 MeV of the MeV pulsars.If the emission efficiency of the non-thermal emission of the MeV pulsars is η MeV ∼ 10 −3 − 10 −2 , it is necessary that 10 ∼ 100 % of the spin-down energy is converted into the particle energy within or near the light cylinder.Moreover, we can also show that the required multiplicity (κ) and Lorentz factor (γ) to explain the observations (E curv ∼ 1 MeV and η MeV ∼ 10 −3 − 10 −2 ) are κ ∼ 10 2−5 and γ ∼ 10 5−6 , respectively.Such a high-multiplicity indicates that the emission originates from so-called "secondary" electrons/positrons that are created by the pair-creation process.To produce the population of the pairs with γ ∼ 10 5−6 and κ ∼ 10 2−5 , it will be required (i) the pair-creation process of most of an 0.1-1 TeV gammarays that are produced with an efficiency of 0.1−1 or (ii) an acceleration of the secondary pairs if they are created by the GeV gamma-rays.

Synchrotron radiation
The GeV gamma-rays emitted in the magnetosphere could be converted into the electron and positron pairs by the photon-photon pair-creation process or magnetic pair-creation process.The created pairs can lose their energy via the synchrotron radiation, which is mainly taken place in the fast-cooling regime in the pulsar magnetosphere.In this section, we consider the possibility that the MeV emission originates from the synchrotron radiation from the created pairs.We will not discuss a detailed model of the pair-creation process, while we will estimate the required Lorentz factor (γ), pitch angle (θ p ) and the multiplicity (κ) of the pairs to explain the spectral properties if the magnetic field strength of   the emission region is specified.Because of insignificant evolution of the pulse shape with the energy bands, we may make the following assumptions: 1. Emission in X-ray to soft gamma-ray bands is produced by the same population of the pairs under the synchrotron cooling.
2. The magnetic field strength and the pitch angle of the pairs do not vary in the emission region.
3. The spectral break at E b ∼ 5 keV bands is caused by either (3-i) the cooling timescale, τ syn = 3m 3 e c 5 /(2e 4 B 2 γ min sin 2 θ p ), equaling to the dynamical timescale of τ dy = r/c or (3-ii) the pairs are already cooled down to the lowest Landau level before running τ dy = r/c.
For the case of (3-i), if the emission region is located at distance r from the center of the neutron star, the required Lorentz factor (γ min ) and pitch angle are , and sin θ p = 4πm e cE b 3γ 2 min ehB , (1) respectively, where B is the magnetic field strength in the emission region.For the case of (3-ii), they are , and sin θ p = 1/γ min , respectively.The maximum Lorentz factor of the pairs, from which the synchrotron cooling starts, are estimated from where E peak is the peak energy in the spectral energy distribution, and E peak = 1 MeV is assumed in this study.Because of the fast cooling regime, the luminosity of the synchrotron emission will be in the order of where Ṅe is the number of the pairs created per unit time.The "required" particle injection rate to explain the observed luminosity is Ṅe = η MeV L sd γ max sin θ p m e c 2 . (5) In this section, we calculate the required multiplicity, κ, from the ratio of Ṅe to a particle injection rate from the polar cap of the Goldreich-Julian value, ṄGJ = 2π 2 B s R 3 N S /(ecP 2 s ), where R N S is the radius of the neutron star.By combining the equations (1) [or (2)], (3) and (5) and using the observed parameters of E b ∼ 5 keV, E peak ∼ 1 MeV and η MeV ∼ 0.01 − 0.001, we can estimate γ min , sin θ p , γ max and Ṅe as a function of the magnetic field strength.We note that the required multiplicity is roughly proportional to (E peak /E b ) 1/2 if we apply the different values for E peak and/or E b [see equation ( 6)].
Figure 6 represents the expected γ min (solid line) and γ max (dashed line) for PSR B1509-58 as a function of the radial distance from the center of the neutron star, where we assume the dipole magnetic field of B(r) = B s (R N S /r) 3 .Figure 7 shows the required multiplicity by assuming the emission efficiency of η MeV = 10 −2 (solid line) and 10 −3 (dashed line), respectively; we do not plot for r/R N S < 4, since γ min calculated with the dipole magnetic field for PSR B1509-58 becomes smaller than unity.In both figures, we can see a change in the slopes of the lines at about r/R N S = 80.This is because that the parameters (γ min , sin θ p ) are determined by the equation (1) for r/R N S > 80, while it is determined by the equation (2) for r/R N S < 80.As Figure 7 shows, we can find that the required multiplicity in r/R N S < 80 for the case (3-ii) is independent on the magnetic field strength of the emission region, and it is expressed by We can see in Figure 7 that the required multiplicity when the emission region is located at r/R N S < 80 is κ > 10 6 .Such a large pair-multiplicity will be realized only by pair-creation cascade developed in a region near the stellar surface.For example, Timokhin & Harding (2019) theoretically investigate the cascade process of the polar cap accelerator and estimate the maximum possible multiplicity as κ max ∼ 3 × 10 6 , although the predicted multiplicity depends on the electric current structure of the acceleration region; for example, the polar cap model penetrated by a return electric current produces a multiplicity of one or two orders of magnitude smaller than κ max (Timokhin & Harding 2015).As indicated by Figure 7, therefore, if the efficiency of PSR B1509-58 is η > 0.001 and the emission originates from r/R N S < 80, the multiplicity being comparable to or more than the predicted maximum possible value will be required.We note that the difference in the estimated multiplicities of the different MeV pulsars is by several factors, assuming the efficiency is a constant.
The pair-creation cascade near the stellar surface can be also initiated if the inward GeV emission produced in the outer magnetosphere irradiates the strong magnetic field region near the stellar surface (Zhang & Cheng 1997;Lin et al. 2009).Such a cascade model can create the pairs with a rate of Ṅe ∼ η GeV L sd /2E es , where η GeV is the efficiency of the inward GeV emission and E es is the maximum photon energy that can escape from the cascade region.Assuming η GeV ∼ 10 −2 and E es ∼ 5 MeV, the multiplicity can be of the order of κ ∼ 10 6 .Uncertainty of the cascade model of the inward emission is the efficiency of the GeV emission and how fraction of the inward emission irradiates near the stellar surface; a detailed modeling will be necessary (Wang et al. 2014).This model predicts that the viewing angle for the MeV pulsar is closer to the spin axis to avoid the GeV emission from the outer magnetosphere.
For PSR B1509-58, the energy corresponding to the cyclotron frequency with B s ∼ 1.5 × 10 13 G is heB s /(2πm e c) ∼ 160 keV.The current synchrotron emission model, therefore, expects that the emission region is located at r/R N S > 4, as indicated in Figure 6.It may be difficult to realize the multiplicity κ > 10 6 in the pair-creation cascade developed at r/R N S > 4 due to the weaker magnetic field strength (Timokhin & Harding 2019).Hence if the observed emission originates from r/R N S > 4 with a multiplicity κ > 10 6 (r/R N S < 80), the emitting pairs must be created in the region located near the stellar surface.Moreover, the pairs created near the surface immediately lose the perpendicular momentum by the synchrotron radiation.Hence, an additional energy injection processed at the vicinity of the emission region is required to gain the perpendicular momentum.For example, Harding & Kalapotharakos (2015) argue that the pairs streaming from the polar cap region gains the pitch angles through the resonant cyclotron absorption of radio waves.
As Figures 6 and 7 show, the required particle multiplicity decreases as the emission region moves farther away from the neutron star, while the required Lorentz factor increases.For example, if the emission is produced near the light cylinder, the maximum Lorentz factor of the pairs is of the order of γ max ∼ 10 5 .To explain the MeV pulsar's emission with the synchrotron emission of the pairs created outer magnetosphere, E GeV ∼ 100 GeV of the photon energy with an emission efficiency 10 −3 − 10 −2 will be required.Such an emission component has not been observationally and theoretically predicted.Hence, we still need to invoke an energy injection to the pairs that happen in the outer magnetosphere.
As we have seen that if the emission of the MeV pulsars originates from the synchrotron radiation from the region inside the light cylinder, the primary GeV emission, which is produced in either the polar cap region or outer magnetosphere, is necessary to create the electron/positron pairs.Within this framework of the scenario, therefore, it is not unexpected that the MeV pulsars also produce a GeV emission in the outer mag-netosphere, which does not direct toward the Earth.Considering that the GeV emission in the outer magnetosphere will be more concentrated to the equatorial plane (Watters & Romani 2011;Kalapotharakos et al. 2018), we speculate that the MeV pulsars have in general a smaller viewing angle than that of the F ermi-LAT pulsars, as suggested by Wang et al. (2014).
We note that the pulse profile of the non-thermal X-ray emission from the F ermi-LAT pulsars (e.g. the Vela pulsar and Crab pulsar) show, in general, narrow and multiple peaks (Kuiper & Hermsen 2015;Kargaltsev et al. 2023).This feature of the pulse profile is different from the signal broad peak of the MeV pulsars.The difference in the pulse profile may also indicate that the viewing geometry and/or the emission region of the MeV pulsars is different from those of the F ermi-LAT pulsars.

Emission from the current sheet
It has been discussed that the current sheet in equatorial plane beyond the light cylinder is the origin of the high-energy emission from the pulsars (Cerutti et al. 2016;Hakobyan et al. 2023;Chernoglazov et al. 2023).In this section, therefore, we estimate the required multiplicity (κ) of the MeV pulsars with a current sheet model for the MeV pulsars.For the particle acceleration due to the magnetic reconnection at the current sheet, the typical energy with the Lorentz factor of the particles may be obtained from the relation that (Lyubarskii 1996;Chernoglazov et al. 2023) where Γ bulk is the Lorentz factor of the the bulk flow of the electron/positron pairs, γ ′ p is the Lorentz factor of the pairs in the flow rest frame and B lc is the magnetic field strength at the light cylinder.The characteristic energy of the synchrotron radiation from the current sheet is written as Using the equations ( 7) and ( 8), the multiplicity is estimated from where we scale the value using the typical bulk Lorentz factor of Γ bulk = 10 (Lyubarskii 1996) and the pulsar parameters of PSR B1509-58.For MeV pulsars, the multiplicity with κ ∼ 10 6 will be required if the emission originates from the current sheet.
The rate of the energy injection to the particles from the current sheet may be written as (Lyubarskii 1996) where x is the radial extent of the energy dissipation region and a represents the sheet width that may be described by a ∼ (3/ √ r e )(c/ω B ) 3/2 , where r e = e 2 /(m e c 2 ) and ω B = eB/m e c is the classical electron radius and gyrofrequency, respectively.By anticipating that the energy injection rate of the equation ( 10) corresponds to the radiation luminosity, L = η MeV L sd , the required radial size of the current sheet is estimated to be which will provide a reasonable size of x ∼ R lc .
For the current sheet model, the difference between F ermi-LAT and MeV pulsars will be attributed to the difference in the multiplicities of the pairs injected into the current sheet.The equation (9) implies that the multiplicity for the F ermi-LAT pulsar is of the order of κ ∼ 10 5 with a photon energy of E peak ∼ 1 GeV, and it is more than one order of magnitude smaller than that of the MeV pulsar.A young pulsar may be observed as a MeV pulsar if the multiplicity of the current sheet is κ ∼ 10 6 , while it becomes a GeV pulsar if κ ∼ 10 5 .One possibility of such a difference in the injected multiplicity is the difference in the inclination angle of the magnetic axis from the spin axis.As we described in section 4.2.2, the large multiplicity of κ ∼ 10 6 can be realized in the cascade process in the polar cap region.Although a detailed study is required, (i) the number of the created pairs in the polar cap region will depend on the magnetic field lines and (ii) how fraction of the created pairs injected into the current sheet depends on the magnetic inclination angle.Intuitively, since a pulsar with a larger magnetic inclination has a smaller angle between the magnetic axis and equatorial plane, it may inject more pairs to the equatorial current sheet.The typical inclination angle of the MeV pulsar, therefore, may be greater than that of the F ermi-LAT pulsars.It will be worth to carry out the study of the polar cap cascade process with the effect of the magnetic inclination angle to investigate this possibility.

SUMMARY
We reported on the pulsed X-ray emission properties of the MeV pulsars, as measured by XMM-Newton, NICER, NuSTAR and HXMT.We fit the broadband X-ray spectra with the s-pl model and bk-pl model.For bk-pl model, the spectrum below the break energy (E b ∼ 5 − 10 keV) shows a hard spectrum with a photon index of Γ 1 , while it above the break energy is ∼ 1.5 (Table 2).In comparison with the X-ray emission of the F ermi-LAT pulsars, the MeV pulsars have in general a harder spectrum and a higher efficiency in 0.3-10 keV energy bands.By assuming the isotropic emission, we found that total efficiency of the non-thermal emission of the MeV pulsars can reach to η MeV ∼ 0.01−0.1,which is similar to the GeV efficiency of the F ermi-LAT pulsars having a similar spin-down power.
We discussed the curvature radiation and synchrotron radiation as a possible emission process of the MeV pulsars.For the curvature radiation process, it was shown that ∼ 10 − 100% of the energy conversion efficiency from the spin-down power to the secondary pairs are required if the emission efficiency is η MeV = 0.001 − 0.01.For the synchrotron radiation inside the right cylinder, if the emission region is r/R lc < 80, the required multiplicity is κ ∼ 10 6 for η MeV ∼ 0.001, while κ ∼ 10 7 for η MeV ∼ 0.01.Such a large pair-multiplicity will be realized only by the pair-creation cascade developed near the stellar surface.If the synchrotron emission in the outer magnetosphere is the origin of the emission of the MeV pulsars, an energy injection to the secondary pairs is necessary.We speculate that if the emission of the MeV pulsars originates from the region inside the magnetosphere, the observer's viewing angle of the MeV pulsars is closer to the spin axis than that of the F ermi-LAT pulsars.
We also discussed the synchrotron emission of the pairs from the current sheet and estimated the required multiplicity as κ ∼ 10 6 for the MeV pulsars.We argued that based on the current sheet model, the difference between the MeV pulsars and F ermi-LAT pulsars is attributed to the difference in the multiplicity of the pairs injected to the current sheet.We speculated that the MeV pulsars have a larger magnetic inclination angle, and inject more pairs to the current sheet than the F ermi-LAT pulsars.

Figure 4 .
Figure 4. Efficiency (η0.3−10) in 0.3-10 keV bands for the MeV pulsars (blue square symbols) and F ermi-LAT pulsars (black circle symbols).The efficiency of the MeV pulsars is rescaled to the phase-averaged value by multiplying the factor of 0.6 to the pulsed flux obtained in this study.The F ermi data is taken from Abdo et al. (2013).

Figure 5 .
Figure 5.Comparison between the emission efficiencies of MeV pulsars (ηX , blue square symbols) in 0.3-150 keV bands and F ermi-LAT young pulsars (black circle symbols) in 0.1-300 GeV bands.The efficiency of the MeV pulsars is rescaled to the phase-averaged value by multiplying the factor of 0.6 to the pulsed flux obtained in this study.The F ermi data is taken from Abdo et al. (2013).

Figure 6 .
Figure 6.Parameters of the synchrotron emission model for the MeV pulsars.The solid line and dashed line show the minimum and the maximum Lorentz factors in the emission region as a function of the assumed radial distance to the emission region from the center of the neutron star.A magnetic dipole field, B(r) = Bs(Rs/r) 3 is assumed to calculate the magnetic field strength at the distance r.The results are for the spin-down parameters of PSR B1509-58.

Figure 7 .
Figure 7.The required multiplicity of the pairs for the synchrotron emission model of PSR B1509-58.The solid line and the dashed line show the cases for the radiation efficiency of ηMeV = 10 −2 and 10 −3 , respectively.

Table 2 .
Spectral parameters for s-pl fitting or bk-pl fitting.

Table A1 .
Journal of the XMM-Newton PN observation.