Middle aged γ-ray pulsar J1957+5033 in X-rays: pulsations, thermal emission and nebula
D. A. Zyuzin, A. V. Karpova, Y. A. Shibanov, A. Y. Potekhin, V. F. Suleimanov
MMNRAS , 1–14 (2020) Preprint 26 January 2021 Compiled using MNRAS L A TEX style file v3.0
Middle aged 𝛾 -ray pulsar J1957+5033 in X-rays: pulsations, thermalemission and nebula D. A. Zyuzin ★ , A. V. Karpova , Y. A. Shibanov , A. Y. Potekhin , V. F. Suleimanov , , Ioffe Institute, Politekhnicheskaya 26, St. Petersburg, 194021, Russia Institut für Astronomie und Astrophysik, Sand 1, D-72076 Tübingen, Germany Kazan (Volga region) Federal University, Kremlevskaja str., 18, Kazan 420008, Russia Space Research Institute of the Russian Academy of Sciences, Profsoyuznaya Str. 84/32, Moscow 117997, Russia
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We analyze new
XMM-Newton and archival
Chandra observations of the middle-aged 𝛾 -ray radio-quiet pulsar J1957+5033. Wedetect, for the first time, X-ray pulsations with the pulsar spin period of the point-like source coinciding by position with thepulsar. This confirms the pulsar nature of the source. In the 0.15–0.5 keV band, there is a single pulse per period and the pulsedfraction is ≈ ± (cid:38) . Γ ≈ .
6. We construct new hydrogen atmosphere models for neutron stars with dipole magnetic fields and non-uniformsurface temperature distributions with relatively low effective temperatures. We use them in the spectral analysis and derivethe pulsar average effective temperature of ≈ ( − ) × K. This makes J1957+5033 the coldest among all known thermallyemitting neutron stars with ages below 1 Myr. Using the interstellar extinction–distance relation, we constrain the distance tothe pulsar in the range of 0.1–1 kpc. We compare the obtained X-ray thermal luminosity with those for other neutron stars andvarious neutron star cooling models and set some constraints on latter. We observe a faint trail-like feature, elongated ∼ Γ = . ± . Key words: stars: neutron – pulsars: general – pulsars: individual: PSR J1957+5033
Neutron stars (NSs) are born in supernova explosions at very hightemperatures ∼ K (e.g. Müller 2020, and references therein).They lose the initial thermal energy via neutrino emission fromtheir interiors and then via photon emission from their surfaces (e.g.Yakovlev & Pethick 2004, and references therein). The cooling rateis sensitive to the physical properties of the superdense matter insideNSs, which are still poorly known (e.g. Yakovlev et al. 2005). Theequation of state (EoS) of such matter can be constrained by compar-ison of cooling theories with NSs surface temperatures derived fromobservational data.The middle-aged PSR J1957+5033 (hereafter J1957) is a radio-quiet 𝛾 -ray pulsar discovered with the Fermi
Large Area Telescope(LAT) (Saz Parkinson et al. 2010). Its spin period 𝑃 =
375 ms andperiod derivative (cid:164) 𝑃 = . × − s s − imply the characteristicage 𝑡 c ≡ 𝑃 / (cid:164) 𝑃 ≈
840 kyr, the spin-down luminosity (cid:164) 𝐸 = . × erg s − and the characteristic (spin-down) magnetic field 𝐵 = . × G . The distance to the pulsar is poorly known. Theonly available estimate, 0.8 kpc, is the so-called ‘pseudo-distance’ ★ E-mail: [email protected] Parameters are calculated using the timing solution for J1957based on five years of the
Fermi data obtained from https://confluence.slac.stanford.edu/display/GLAMCOG/LAT+ obtained using empirical correlation between the distance and the 𝛾 -ray flux above 100 MeV (Abdo et al. 2013). It is known to be uncertainby a factor of 2–3. Analysing the 𝛾 -ray pulse profile, Pierbattista et al.(2015) estimated the magnetic inclination and line of sight angles ofJ1957 for different 𝛾 -ray emission geometries.The pulsar X-ray counterpart was identified by position in the25-ks Chandra
Advanced CCD Imaging Spectrometer (ACIS-S) ob-servation (Marelli et al. 2015). Its spectrum in the 0.3–10 keV bandwas found to be well described by the absorbed single power-law(PL) with a photon index Γ = . ± .
3, an energy integrated unab-sorbed flux 𝐹 𝑋 = ( . ± . ) × − erg s − cm − and an absorptioncolumn density 𝑁 H < . × cm − (Marelli et al. 2015).We reanalyzed the Chandra data and confirmed the results byMarelli et al. (2015). However, we found an unexpectedly large countnumber in the 0.1–0.3 keV band – 30 against 90 counts detected in the0.3–10 keV band. This indicates the presence of a second soft com-ponent in the pulsar spectrum likely related to the thermal emissionfrom the surface of the NS with a very low effective temperature. Inthis case, J1957 becomes especially interesting for comparison withNS cooling theories according to which at its characteristic age the
Gamma-ray+Pulsar+Timing+Models . See also Kerr et al. (2015)for details. ObsID 14828, PI M. Marelli, observation date 2014-02-01. © a r X i v : . [ a s t r o - ph . H E ] J a n Zyuzin et al.
N E . - k e V A J e bu l a? C h a nd r a / A C I S - S : . : . : : . Galactic longitude G a l ac ti c l a tit ud e M O S + M O S +P N n e bu l a? X MM - N e w t on N E A J : . : . : . : : . : : . Galactic longitude G a l ac ti c l a tit ud e Figure 1.
Images of the J1957 field.
Left: Chandra
ACIS-S3 chip field-of-view in the 0.5–2 keV band.
Right:
18 arcmin ×
18 arcmin
XMM-Newton combinedMOS1+MOS2+PN image (red = 0.2–0.5 keV, green = 0.5–2 keV). J1957, the nearby star ‘A’ and the presumed trail-like nebula are marked. Compassescorrespond to the equatorial coordinate system. pulsar should have already passed from a relatively slow neutrinocooling stage to a significantly faster photon stage where observa-tional data on thermal emission from cooling NSs are particularlyscarce. Unfortunately, the soft component cannot be confirmed usingthe
Chandra data since the ACIS energy scale is not calibrated below0.3 keV. 3.2 s time resolution of the observations does not also allowto detect pulsations with the pulsar period. Therefore, we performeddedicated
XMM-Newton observations of J1957. Here we present theanalysis of these data. The X-ray data are described in Section 2.Timing and spectral analysis of J1957 are presented in Sections 3and 4. The results are discussed in Section 5 and summarized inSection 6. Some details of the analysis are given in the Appendices.
The J1957 field was observed with
XMM-Newton on 2019 October 5(ObsID 0844930101, PI D. Zyuzin). The total exposure was about 87ks. The European Photon Imaging Camera Metal Oxide Semiconduc-tor (EPIC-MOS) detectors were operated in the full-frame mode withthe imaging area of about 28 arcmin ×
28 arcmin and the mediumfilter and the EPIC-pn camera – in the large window mode with theimaging area of about 13.5 arcmin ×
26 arcmin and the thin filter.We also used the
Chandra /ACIS-S archival dataset (ObsID 14828)where the pulsar was exposed on the S3 chip. To analyze the data,we utilized the
XMM-Newton
Science Analysis Software (xmm-sas)v. 17.0.0 and
Chandra
Interactive Analysis of Observations (ciao)v. 4.12 packages.The
Chandra dataset was reprocessed using the chandra_reprotool. Applying the fluximage task, we created the exposure-corrected image of the ACIS-S3 chip which is presented in the leftpanel of Fig. 1 where the pulsar counterpart and the nearby star ‘A’ aremarked. The wavdetect command was used to obtain coordinates ofpoint-like sources. For J1957, we derived R.A. = 19 h m . s ◦ (cid:48) . (cid:48)(cid:48) 𝜎 purestatistical uncertainties).We combined the data from both MOS and PN detectors to obtaindeeper X-ray images using the ‘images’ script (Willatt & Ehle 2016).The resulting image is shown in the right panel of Fig. 1. Since XMM-Newton has lower spatial resolution than
Chandra , J1957 issomewhat blurred with the star ‘A’. One can see a faint thin featureprotruding from J1957 almost perpendicularly to the Galactic plane.It has a clumpy structure and is also visible in the
Chandra imagewhere it extends up to the edge of the ACIS-S detector ( ≈ XMM-Newton image its length seems to be longer, at least up to ∼ − for pn and 0.6counts s − for both MOS cameras (see Fig. 2). The resulting effectiveexposures are about 79.8, 79.8 and 48.7 ks for the MOS1, MOS2and pn detectors, respectively. We selected single to quadruple pixelevents (pattern ≤
12) for the MOS data and single and double pixelevents (pattern ≤
4) for the pn data.
The timing resolution of the EPIC-pn detector operating in the largewindow mode is ≈
48 ms which allows us to search for pulsations .MNRAS , 1–14 (2020) iddle aged 𝛾 -ray pulsar J1957+5033 in X-rays Time [ks] . . . . . . . R a t e [ c n t / s ] MOS1MOS2
Time [ks] R a t e [ c n t / s ] PN Figure 2.
High-energy light curves obtained from the FOVs of MOS (10–12 keV; left ) and pn (12–14 keV; right ) detectors. Dash-dotted lines indicate thresholdsused to filter out periods of background flares. . . . . . . . . . Frequency [Hz] Z Figure 3. 𝑍 -test periodogram for J1957. Dashed lines show confidencelevels. from J1957. We used the event list unfiltered from background flaressince multiple time gaps in the data may hamper signal detection.The barycenter correction was applied by the xmm-sas barycencommand using DE405 ephemeris and the pulsar coordinates de-rived from the Chandra data. We extracted events in the 0.15–0.5keV band using the 12.5-arcsec radius circle around the
Chandra position of J1957. Such a small aperture was chosen to eliminatethe contribution from the star ‘A’ located at about 20 arcsec fromthe pulsar. We searched for pulsations utilising 𝑍 -test (Buccheriet al. 1983) and 2.667–2.669 Hz frequency range, encapsulating thepredicted pulsar rotation frequency of about 2.6680281 Hz obtainedfrom extrapolation of the Fermi timing solution to the epoch of the
XMM-Newton observations (MJD 58761), and with a step of 0.1 𝜇 Hz.The resulting periodogram is shown in Fig. 3. The maximum 𝑍 is ≈ 𝐶 = [ − N exp (− Z , max / )] × = .
996 per cent (or ≈ 𝜎 ) where N = Δ 𝑓 × 𝑇 obs is the number of statistically independent trials, Δ 𝑓 is the frequency range and 𝑇 obs is the duration of the observation.The corresponding frequency is 2.6680249(12) Hz (the frequency1 𝜎 uncertainty was calculated using the formula from Chang et al. C o un t s / b i n . . . . . . . . . . . Phase
Figure 4.
Folded X-ray light curves for J1957 in different energy bandsindicated in the panels (periods of background flares are removed). Dottedlines in the middle and bottom panels indicate the background level. In the0.15–0.5 keV band the background contribution ( ≈
12 counts phase bin − ) isvery low in comparison with the pulsar count rate and thus it is not shown.The best-fitting sine curve is overlaid. 𝜎 with the predicted value from Fermi timing solution and firmly establishes the pulsar nature of theX-ray source.The J1957 X-ray pulse profile obtained using the derived frequencyis presented in Fig. 4. We see that pulsations are clearly detectedonly in the very soft band. Non-detection in harder bands is likelydue to the low number of counts from the pulsar. The background-corrected pulsed fraction in the 0.15–0.5 keV band PF = ( 𝑓 max − 𝑓 min )/( 𝑓 max + 𝑓 min ) , where 𝑓 max and 𝑓 min are the maximum andminimum intensity of the folded light curve, is ≈ ± ≈ ≈
64 per cent in the 1–10 keVband.
MNRAS , 1–14 (2020)
Zyuzin et al. C hand r a / AC I S - S . - ke VA J1957 : . : : . : . Right ascension D ec li n a ti on Figure 5.
Chandra
ACIS image of the J1957 field in the 0.5–2 keV band.J1957 and the nearby star ‘A’ also marked. The background for the pulsar wasextracted from the dashed circle. The solid rectangle shows the region used to extract the spectrum of the trail-like nebula while the dashed rectangleswere used for the background. Stars excluded from the regions are shown bycrossed circles.
We extracted the time integrated pulsar spectra from both the
XMM-Newton and
Chandra data. In the latter case we utilized the specex-tract routine and the 2.5-arcsec radius aperture. The
XMM-Newton spectra were extracted from the 12.5-arcsec radius circle using evse-lect task and redistribution matrix and ancillary response files werecreated by rmfgen and arfgen commands. The background spec-trum was obtained from the source-free region (see Fig. 5). For theinterstellar medium (ISM) absorption, we applied tbabs model withthe wilm abundances (Wilms, Allen & McCray 2000). The spec-tra were fitted simultaneously in the X-Ray Spectral Fitting Packagexspec v.12.10.1 (Arnaud, Gordon & Dorman 2018). We used thefollowing energy ranges: 0.3–10 keV for the Chandra , 0.2–10 keVfor the MOS and 0.15–10 keV for the pn spectra. The resulting num-ber of source counts after background subtraction is 232(MOS1) +254(MOS2) + 902(pn) + 88(ACIS).As a first step, to check how different models fit the data, weapplied the 𝜒 -statistic and grouped the data to ensure 25 countsper energy bin. The single absorbed PL (which describes the pulsarnon-thermal emission of magnetospheric origin) or blackbody (BB,which describes the thermal emission from the NS surface) modelsresulted in unacceptable fits with reduced 𝜒 𝜈 = 2.36 and 5.5 for 54degrees of freedom (dof), respectively. Then we tried the compositeBB + PL model and found that it fits the spectra well with 𝜒 𝜈 = 1.13(52 dof).For the thermal component, we also tried the NS magnetized atmo-sphere models nsmaxg (Ho, Potekhin & Chabrier 2008) presentedin the xspec package assuming an NS mass 𝑀 NS = 1.4M (cid:12) and radius 𝑅 =
13 km. However, in this case the temperature tends to the lowest https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/ value available for these models lg ( T / K ) = . ; here and hereafterlg ≡ log . Therefore, we calculated another grid of NS atmospheremodels, nsmdip, which include lower temperatures for the magneticfields that seem to be likely for this pulsar. We assumed a dipole mag-netic field (taking the effects of General Relativity into account) anda corresponding distribution of the local effective temperature overthe stellar surface. In these models, the angle 𝛼 between the rotationand the magnetic axes and the angle 𝜁 between the rotation axis andthe line of sight can be used as free parameters. The total thermalluminosity 𝐿 ∞ that would be measured by a distant observer is cal-culated by a proper surface integration and converted into the globaleffective temperature 𝑇 ∞ . The details of calculations are presentedin Appendix A.We have tested several nsmdip models (Table A1). In model ns-mdip1, we have assumed the canonical NS mass 𝑀 = . (cid:12) andmagnetic field at the pole 𝐵 p = × G, which corresponds tothe characteristic magnetic field 𝐵 ≈ . × G of J1957 at theequator derived from the spin-down. The radius 𝑅 = . 𝑧 g = .
22, and surface gravity 𝑔 s = . × cm s − . In order to test a higher redshift, we considera more compact NS model (nsmdip2) with 𝑀 = (cid:12) and 𝑅 = . 𝑧 g = .
44 and 𝑔 s = × cm s − ), which approximatelycorresponds to the EoS BSk26. Model nsmdip3 has the same 𝑧 g and 𝑅 as model 1 but a more consistent magnetic field estimate. Thecanonical characteristic magnetic field is the equatorial field of anorthogonal rotating dipole in vacuum with the given spin period andits derivative, assuming 𝑅 =
10 km and 𝐼 =
1, where 𝐼 is themoment of inertia in units of 10 g cm (e.g., Manchester & Tay-lor 1977). However, according to the EoS BSk24, for 𝑀 = . (cid:12) we have 𝑅 = . 𝐼 = .
51. Moreover, the spin-downof a pulsar is affected by its magnetosphere. Results of numericalsimulations of plasma behaviour in the pulsar magnetosphere sug-gest that the characteristic magnetic field should be multiplied bya factor of 0 . ( 𝑅 /
10 km ) − ( + sin 𝛼 ) − / √ 𝐼 , where 𝛼 is theangle between rotational and magnetic axes (Spitkovsky 2006). Forthe above-mentioned values of 𝑅 and 𝐼 this implies a 2–3 timesweaker field compared to the characteristic one. Thus we adopted 𝐵 p = . × G in model nsmdip3 which corresponds to the fieldstrength at the equator 𝐵 eq ≈ × G, that is 2.7 times smallerthan the pulsar characteristic (spin-down) field. We find that all threeabsorbed nsmdip + PL models describe the data equally well as theBB + PL model giving 𝜒 𝜈 ≈ 𝑊 -statistic (Arnaud et al.2018) which is 𝐶 -statistic (Cash 1979) suitable for Poisson data withPoisson background. We then performed the fitting using a Markovchain Monte-Carlo (MCMC) sampling procedure. We employed theaffine-invariant MCMC sampler developed by Goodman & Weare(2010) and implemented in a python package emcee by Foreman-Mackey et al. (2013). In addition, to estimate the distance to J1957,we used the interstellar absorption–distance relation towards the pul-sar as a prior (see Appendix B for details). About 100 walkers and13000 steps were typically enough to ensure fit convergences. Usingthe sampled posterior distribution, we obtained the best-fitting pa-rameters of the models with uncertainties, which are defined as theirmaximal-probability density values and respective credible intervals.In the case of nsmdip + PL models, we found that our data arerather insensitive to the 𝛼 and 𝜁 angles. From Fig. 4, one can see MNRAS000
10 km ) − ( + sin 𝛼 ) − / √ 𝐼 , where 𝛼 is theangle between rotational and magnetic axes (Spitkovsky 2006). Forthe above-mentioned values of 𝑅 and 𝐼 this implies a 2–3 timesweaker field compared to the characteristic one. Thus we adopted 𝐵 p = . × G in model nsmdip3 which corresponds to the fieldstrength at the equator 𝐵 eq ≈ × G, that is 2.7 times smallerthan the pulsar characteristic (spin-down) field. We find that all threeabsorbed nsmdip + PL models describe the data equally well as theBB + PL model giving 𝜒 𝜈 ≈ 𝑊 -statistic (Arnaud et al.2018) which is 𝐶 -statistic (Cash 1979) suitable for Poisson data withPoisson background. We then performed the fitting using a Markovchain Monte-Carlo (MCMC) sampling procedure. We employed theaffine-invariant MCMC sampler developed by Goodman & Weare(2010) and implemented in a python package emcee by Foreman-Mackey et al. (2013). In addition, to estimate the distance to J1957,we used the interstellar absorption–distance relation towards the pul-sar as a prior (see Appendix B for details). About 100 walkers and13000 steps were typically enough to ensure fit convergences. Usingthe sampled posterior distribution, we obtained the best-fitting pa-rameters of the models with uncertainties, which are defined as theirmaximal-probability density values and respective credible intervals.In the case of nsmdip + PL models, we found that our data arerather insensitive to the 𝛼 and 𝜁 angles. From Fig. 4, one can see MNRAS000 , 1–14 (2020) iddle aged 𝛾 -ray pulsar J1957+5033 in X-rays . . . . p . d . f . . . l g T ∞ [ K ] . . . p . d . f . − l g D [ k p c ] . . . p . d . f . α [ d e g ] . . . p . d . f . ζ [ d e g ] . . . p . d . f . . . Γ . . . p . d . f . .
025 0 . N H [10 cm − ] − . − . l g K [ ph / ( c m s k e V ) ] . . lg T ∞ [K] − lg D [kpc] α [deg] ζ [deg] . . Γ − . − . lg K [ph/(cm s keV)] . . . p . d . f . Figure 6.
1D and 2D marginal posterior distribution functions (p.d.f.) for parameters of the nsmdip1 + PL model (see Table 1). Vertical dashed lines in 1Ddistributions indicate the best-fitting values while light gray strips show 1 𝜎 credible intervals. In 2D distributions, 40, 68, 90 and 99 per cent confidence contoursare shown. that there is one pulse per period. This implies that 𝛼 + 𝜁 ≤ ◦ ,which was used as a prior in the spectral fitting. As an example,1D and 2D marginal posterior parameter distribution functions forthe nsmdip1 + PL model are shown in Fig. 6. One can see that theangles, in contrast to other parameters, cannot be well constrainedfrom the fit. Their formal best-fitting values are close to 0 ◦ . Ourfits show that the similar situation occurs for the nsmdip2 + PL andnsmdip3 + PL models. However, for any of the models the probabilitydensity function for angles remains relatively high for the wholemeaningful ranges of 0 ◦ – 90 ◦ . This results in a very large formaluncertainties of the angles, and in fact they can accept any value. Forthese reasons, angles are not included in the parameter list of Table 1 presenting the fit results. We found that other fit parameters dependon a specific value of angles only weakly. We also fitted spectra usingtwo models of the J1957 𝛾 -ray emission geometry, polar cap (PC)and slot-gap (SG), obtained from the Fermi data by Pierbattista et al.(2015) which satisfy the condition 𝛼 + 𝜁 ≤ ◦ . The results for the Pierbattista et al. (2015) used polar cap (PC), slot gap (SG), outer gap(OG) and one pole caustic (OPC) models which assume different regions ofthe pulsar magnetosphere where particles are accelerated and emit 𝛾 -rays.In the PC model this occurs at low altitudes near the magnetic poles andin the SG model – in a slot gap, which is a narrow gap extending from thepolar cap surface to the light cylinder. In the OG model the gap extendsMNRAS , 1–14 (2020) Zyuzin et al.
Table 1.
Best-fitting parameters for different models with thawed 𝛼 and 𝜁 angles and accounting for 𝛼 + 𝜁 <= 90 ◦ .Model 𝑁 H , 𝑘 B 𝑇 ∞ , 𝑅 ∞ , lg 𝐿 ∞ , Γ 𝐾 , 𝐷 , − ln L BIC10 cm − eV km erg s − ph cm − s − keV − pcBB + PL 2 . + . − . + − . + . − . . + . − . . + . − . . + . − . × − + − . + . − . . + . − . 𝑓 . + . − . . + . − . . + . − . × − + − . + . − . . + . − . 𝑓 . + . − . . + . − . . + . − . × − + − . + . − . . + . − . 𝑓 . + . − . . + . − . . + . − . × − + − † 𝑁 H is the absorbing column density, 𝑇 ∞ = 𝑇 /( + 𝑧 g ) is the effective temperature as measured by a distant observer, 𝑅 ∞ = 𝑅 ( + 𝑧 g ) is the radius of theequivalent emitting sphere as seen by a distant observer, 𝐿 ∞ = 𝐿 /( + 𝑧 g ) = 𝜋𝑅 𝜎 SB 𝑇 /( + 𝑧 g ) is the bolometric thermal luminosity as measuredby a distant observer ( 𝜎 SB is the Stefan-Boltzmann constant), Γ is the photon index, 𝐾 is the PL normalization and 𝐷 is the distance. For the gravitationalredshift 𝑧 g and unredshifted radius 𝑅 , see Table A1. The last two columns give the values of the maximum log-likelihood ln L and the Bayesian informationcriterion (BIC). All errors are at 1 𝜎 credible interval. 𝑓 Fixed parameters.
Table 2.
The same as in Table 1 but for fixed angles at 𝛼 = 66 ◦ and 𝜁 = 24 ◦ † .Model 𝑁 H , 𝑘 B 𝑇 ∞ , 𝑅 ∞ , lg 𝐿 ∞ , Γ 𝐾 , 𝐷 , − ln L BIC10 cm − eV km erg s − ph cm − s − keV − pcnsmdip1 + PL 3 . + . − . . + . − . 𝑓 . + . − . . + . − . . + . − . × − + − . + . − . . + . − . 𝑓 . + . − . . + . − . . + . − . × − + − . + . − . . + . − . 𝑓 . + . − . . + . − . . + . − . × − + − † Notations are the same as in Table 1.
SG with 𝛼 = ◦ and 𝜁 = ◦ are presented in Table 2 and Fig. 7.This model geometry is more preferable as it gives an acceptableX-ray pulse shape and pulsed fraction (see Section 5 for details).To understand which of the models is statistically more preferable,Tables 1 and 2 also present a Bayesian evidence L and informa-tion criteria BIC = 𝑛 ln 𝑁 b − L , where 𝑛 is the number of freemodel parameters and 𝑁 b is the number of spectral bins (which is373). When picking from several models, the one with the lowestBIC is preferred. As seen, the BB + PL model has the smallest BICand appears to be preferable. The strength of the evidence againstnsmdip + PL models with the higher BICs is defined by Δ BIC > ∼ verystrong evidence . At the same time, Δ BICs for any pair of nsmdip + PLmodels in Table 1 is < ∼
2, which is qualified either as weakly positive or not worth more than a bare mention . On the other hand, fixingof 𝛼 and 𝜁 almost does not change L while BICs values becomesmaller. In this case, the difference between nsmdip + PL and BB+PLmodels Δ BIC < ∼ verystrong to decisive . Moreover, the models nsmdip+PL are preferredfrom the physics point of view, because they assume the plausibleNS radii and temperature distributions, whereas the BB+PL modelresults in a best-fitting radius incompatible with thermal emissionfrom the entire NS surface (see the discussion in Section 5). from the null charge surface to high altitudes along the last-open-field lines.The OPC is a variation of the OG model which suggests different gap widthand energetics. The SG and OG can provide wide 𝛾 -ray beams and implyexponential spectral cut-off at high energies while the PC provides narrowbeams and predicts super-exponential spectral cut-off due to magnetic pairproduction process. (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) P ho t on s c m (cid:0) s (cid:0) k e V (cid:0) XMM/MOS1XMM/MOS2XMM/PNChandra/ACIS χ Energy [keV]
Figure 7.
The J1957 unfolded spectrum and the model nsmdip1 + PL withfixed angles at 𝛼 = 66 ◦ and 𝜁 = 24 ◦ (see the first line in Table 2). Dotted linesshow the model components. Data from different instruments are indicatedby different colours as indicated in the panel. Spectra were grouped to ensureat least 10 counts per energy bin for illustrative purposes. For the spectral analysis of the trail-like nebula, we used only the
Chandra data set since in the
XMM-Newton data it is somewhatblurred with several stars. We extracted spectrum of the trail usingthe 30 arcsec ×
300 arcsec rectangle region shown in Fig. 5 togetherwith the regions used for the background. It was binned to ensureat least 25 counts per energy bin and fitted in the 0.3–10 keV band.The total number of counts in the spectrum is 811 while only 105of them come from the source. We tried the PL model and found
MNRAS000
MNRAS000 , 1–14 (2020) iddle aged 𝛾 -ray pulsar J1957+5033 in X-rays (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) no r m a li ze d c oun t s s (cid:0) k e V (cid:0) χ Energy [keV]
Figure 8.
Chandra spectrum of the trail-like nebula, the best-fit PL modeland residuals. that the column density is highly uncertain. Thus, we fixed it at3 × cm − which is compatible with all pulsar models (seeTable 1). The resulting parameters Γ = . ± . 𝐾 = . + . − . × − ph cm − s − keV − , the unabsorbed flux in the 0.3–10 keV band 𝐹 𝑋 = . + . − . × − erg s − cm − and 𝜒 𝜈 = .
17 (dof = 29).Such spectrum can have a synchrotron nature. The spectrum andthe best-fit model are shown in Fig. 8. The spectrum can alsobe equally well described ( 𝜒 𝜈 = .
16, dof = 29) by the thermalbremsstrahlung model with a temperature 3 . + . − . keV and theunabsorbed flux in the 0.3–10 keV band 𝐹 𝑋 = . + . − . × − erg s − cm − . The time integrated X-ray spectrum of J1957 the in 0.15–10 keVrange can be well described by the composite model consisting ofthe thermal and PL components while the single PL model suggestedpreviously by Marelli et al. (2015) for a more narrow range of 0.3 – 10keV is statistically unacceptable in the extended range. In the case ofthe BB + PL model, the obtained effective temperature 𝑇 ∞ ≈
54 eV(6 . × K) is typical for the emission from the bulk of an NSsurface but the radius of the emitting area is smaller than an expectedNS radius of 10 – 15 km (see Table 1). The thermal emission could beproduced by hot polar caps of the NS heated by relativistic particlesfrom pulsar magnetosphere. For J1957, the ‘standard’ pulsar polarcap radius is about 0.3 km (Sturrock 1971) which is compatiblewith the lower bound of the derived radius. However, the derivedtemperature is too low for the polar cap emission (cf. Potekhin et al.2020), rejecting this possibility. Thus, the BB + PL model mightdescribe some hotter part of the pulsar surface while the other partis cooler and not observed in X-rays as takes place, e.g., for themiddle-aged PSR B1055 −
52 (Mignani, Pavlov & Kargaltsev 2010).Note, however, that such interpretation implies that the obtained BBtemperature cannot be used for a comparison with the predictions ofthe cooling theory; instead, the bolometric thermal luminosity shouldbe used for such a comparison (as discussed, e.g., in Potekhin et al.2020). Remarkably, the best-fitting thermal luminosities given by theBB+PL and nsmdip+PL models are compatible within uncertainties(see Table1).The nsmdip + PL models also give acceptable fits. Although they are possibly less preferable according to the Bayesian criteria thatdisregard prior theoretical constraints on NS radii, they are morephysically motivated. In particular, they are based on the magneticatmosphere models, computed for realistic NS parameters, and theyconsistently take into account the distributions of temperature andmagnetic field over the NS surface. Combining the results fromall these models (Table 1), the estimated redshifted NS effectivetemperature 𝑇 ∞ ≈ ± . ± .
05 MK). This makes the pulsarone of the coldest among all known NSs with estimated thermalluminosities, whose measured thermal emission from the surface ispowered by cooling (see Section 5.3). We cannot constrain the pulsarviewing geometry (angles 𝛼 and 𝜁 ) from the time integrated spectra.For the first time, we detected X-ray pulsations with the pulsar spinperiod. The pulsations are significant only in the soft band of 0.15–0.5 keV where the thermal emission component strongly dominatesin the spectrum of the pulsar. The pulse profile is a sine-like with asingle pulse per period and the pulsed fraction of ≈ ± 𝛼 and 𝜁 that the rotation axis makes with the magnetic axis andthe line of sight (magnetic obliquity and pulsar inclination), assumingthat the flux in the photon energy band 0.15–0.5 keV is dominatedby thermal emission. In the left panel of Fig. 9, the hatched regionscorrespond to the 𝛼 and 𝜁 values that provide the observed pulsedfractions of thermal radiation from 0.12 to 0.24 in the 0.15–0.5 keVenergy band for the atmosphere models nsmdip 1–3, according to thelegend. We see that for all the considered cases 𝛼 + 𝜁 > ∼ ◦ . Because ofa higher redshift in model nsmdip 2, it shows a stronger gravitationallight bending, which leads to a stronger smearing of the light curves,compared with the two other models. Therefore, higher 𝛼 and 𝜁 areneeded to obtain the same pulsed fraction. In particular, the upperlimit of 24 per cent is never reached (that is why the hatched area hasa shape of a bell rather than a horseshoe in this case). The absenceof an inter-pulse on the observed light curve in Fig. 4 suggests that 𝛼 + 𝜁 (cid:46) ◦ . The curves in Fig. 9, left, show the combinations of 𝛼 and 𝜁 that provide the pulsed fraction 18 per cent. For comparison, weshow an analogous line for the best-fitting BB model from Table 1( 𝑅 ∞ = . 𝑀 = . (cid:12) and 𝑅 = . 𝛼 + 𝜁 values. For nsmdip 2, the above-mentioned smearingof the light curves due to the light bending effect is so strong thatthe pulsed fraction does not exceed 5 per cent, meaning that the BBmodel is inapplicable in this case. We also show the tentative 𝛼 and 𝜁 values obtained by Pierbattista et al. (2015), who fitted geometricalmodels to the observed 𝛾 -ray pulse profile, using different theoreticalmodels of 𝛾 -ray pulse formation. We see that one of the models (theslot gap model) is compatible with nsmdip 1 and nsmdip 2. Examples MNRAS , 1–14 (2020)
Zyuzin et al.
Figure 9.
Left : Regions of magnetic obliquity 𝛼 and rotation axis inclination 𝜁 (hatched regions), compatible with the observed pulsed fraction 0 . ± .
06 ofemission in the soft X-ray energy band 0.15–0.5 keV, and the 𝛼 and 𝜁 combinations (curves) that provide the pulsed fraction 0.18, shown for the three nsmdipatmosphere models according to the legend. For comparison, a curve for a BB model (see text) is also plotted. The diagonal line separates the parts of thehatched regions in the lower triangle, which are compatible with the singe-pulsed light curves, as observed in Fig. 4. The violet diamonds show the tentative 𝛼 and 𝜁 values obtained by Pierbattista et al. (2015) using polar cap (PC), slot gap (SG), outer gap (OG) and one pole caustic (OPC) models of gamma-ray pulseformation. Right : Light curves in the redshifted energy band 0.15–0.5 keV, computed for models nsmdip 1, 2 and 3 at lg 𝑇 ∞ (K) = . 𝛼 and 𝜁 angles shown by the respective points A, B and C in the left panel (solid, dashed, and dot-dashed lines). For comparison, the bottom panelshows the light curves computed for the best-fit BB model. The violet dotted line in each panel shows the light curve for 𝛼 = ◦ and 𝜁 = ◦ given by the fitwith the SG model (Pierbattista et al. 2015). of theoretical light curves, calculated as described in Appendix A,are shown in the right panel of Fig. 9. For each model we show threecases, marked by letters A, B and C, corresponding to the magneticobliquity 𝛼 and inclination 𝜁 values shown in the left panel. Note,that case C does not agree with the observed profile since it predictsthe presence of an inter-pulse and thus can be excluded.The low number of counts does not allow us to perform the phase-resolved spectral analysis which could help to distinguish betweendifferent geometries and spectral models. Deeper X-ray observationsare necessary to do that. Ultraviolet (UV) observations could alsoclarify the situation whether the thermal emission is best descried bythe nsmdip or by the BB (solid state) spectral models as they predictdifferent fluxes in this range.Implementation of the extinction–distance relation in the fittingprocedure allowed us to estimate the distance to J1957. We note, thatdue to the large relative uncertainties in this relation especially atlow distances (see Appendix B), the constraints on 𝐷 are rather weakespecially for the BB + PL model where 𝑁 H uncertainties are largerthan for the atmosphere models. Thus, the BB + PL model gives thedistance up to ≈ 𝛾 -ray efficiency 𝜂 𝛾 = 𝐿 𝛾 / (cid:164) 𝐸 of ≈ 𝐿 𝛾 = 𝜋𝐷 𝐺 is the 𝛾 -ray luminosityand 𝐺 = . × − erg s − cm − (Marelli et al. 2015) is the 𝛾 -ray flux above 100 MeV.As for the PL spectral component, the derived photon index rangeof 1.5–1.9 is typical for pulsars (e.g. Kargaltsev & Pavlov 2008).For all the models, the unabsorbed non-thermal flux in the 2–10 keV band is 𝐹 𝑋 ≈ × − erg s − cm − . This corresponds toX-ray luminosity 𝐿 𝑋 of ≈ ( . − . ) × erg s − for thedistance range of 100–400 pc provided by the atmosphere modelsand up to 2 . × erg s − erg s − with the upper bound on thedistance of 1 kpc as provided by the BB model. The correspondingX-ray efficiency of 𝜂 𝑋 = 𝐿 𝑋 / (cid:164) 𝐸 = − . − − . and up to10 − . , respectively. These values are compatible with empiricaldependencies of pulsar X-ray nonthermal luminosity and efficiencyvs. characteristic age (see e.g. Zharikov, Shibanov & Komarova 2006;Zharikov & Mignani 2013). Upper limits on the pulsed fractionin harder bands of ≈ −
60 per cent do not give any additionalinformative constraints on the properties of the non-thermal X-rayemission from the pulsar magnetosphere.
As was noted in Section 2, there is a thin straight feature protrudingnorth-west from the pulsar (in equatorial coordinates; Fig. 5) likelyassociated with it. Its spectrum can be well described by the PL modelwith parameters typical for X-ray synchrotron nebulae powered bypulsars, known as PWNe. Its length is ≈ ≈ . − . − MNRAS , 1–14 (2020) iddle aged 𝛾 -ray pulsar J1957+5033 in X-rays et al. 2017b, and references therein). If so, we do not see any hintof a ‘normal’ tail-like PWN protruding behind the pulsar in X-rayimages. Torus or jet structures typical for PWNe of younger pulsarsare also not detected. However, the situation can be similar to theGuitar nebula powered by PSR B2224+65: the guitar-shaped bow-shock nebula was detected in H 𝛼 while no head-tail PWN was foundwithin the shock in X-rays (Reynolds et al. 2017). Instead, X-rayobservations revealed a jet-like feature inclined by ≈ ◦ to thep.m. direction of B2224+65. For our pulsar, no H 𝛼 emission wasdetected around it with the 3.6 m WIYN telescope at a 300 s exposure(Brownsberger & Romani 2014). However, the exposure may be tooshort to detect a fainter bow-shock at a high galactic latitude of ≈ ◦ where the density of the interstellar matter is small and deeperobservations are needed.The tail interpretation suggests that J1957 should move towardsthe Galactic plane (Fig. 1) raising a question about its birth sitesomewhere in the Galactic halo. For the misaligned outflow, thedirection of p.m. remains unclear. In the latter case, we can assumethat J1957 was born in the Galactic disk since this possibility ismore plausible than the birth in the halo. Due to the natal kick inthe supernova explosion it could achieve high velocities and moveaway from the disk. In such a case, we can estimate the J1957 p.m.perpendicular to the Galactic plane. At the pulsar galactic latitude of ≈ ◦ the p.m. 𝜇 ≈ 𝑡 − mas yr − , where 𝑡
840 kyr is the pulsartrue age normalized to its characteristic age 𝑡 c . This correspondsto the transverse velocity 𝑣 ≈ 𝑡 − 𝐷 kpc km s − which iscompatible with the pulsars velocity distribution (Hobbs et al. 2005)for the estimated distance range even if the true age is ∼ . + . − . keV.In this case the nebula emission comes from the shocked ISM andthe pulsar p.m. must be aligned with the nebula axis. Such situationappears for the tail of PSR J0357+3205 (Marelli et al. 2013). How-ever, due to low count statistics it is hard to distinguish between thePL and the thermal models in our case.Measurement of the pulsar’s p.m. is necessary to understand thenature of the feature and deeper X-ray observations are needed toconstrain its shape and spectral properties. As noted above, J1957 can be one of the coldest cooling NSs withmeasured surface temperatures (e.g., Potekhin et al. 2020 and ref-erences therein). In Fig. 10 we compare the estimated thermal lu-minosities of different isolated NSs with theoretical cooling curves(i.e., redshifted luminosities as functions of ages). For J1957 we usethe luminosity estimates reported in Table 1. The error bar unitesthe uncertainties for the models nsmdip 1–3. The model nsmdip 1is adopted as the best estimate, because it provides the lowest BICvalue among these three models. The observational estimates for theother cooling NSs are taken from Potekhin et al. (2020). For theNSs lacking timing-independent age estimates, including J1957, weplot 𝐿 ∞ against their characteristic ages 𝑡 c . In these cases the left-ward arrows indicate that 𝑡 c is likely to be larger than the true age,which is the common case, although there are exceptions where 𝑡 c is somewhat lower than the true age (see, e.g., examples, discussionand references in Potekhin et al. 2020). We see that J1957 has thelowest thermal luminosity among all cooling NSs with ages < Figure 10.
Cooling curves of NSs with different masses (coded with color),described by the EoS BSk24, for accreted and non-accreted heat blanketingenvelopes, compared with observations. The data are plotted as indicated inthe legend for NS classes following Potekhin et al. (2020): weakly magnetizedthermally emitting isolated NSs (TINS), X-ray isolated NSs (XINS), high-B pulsars (HB PSRs) and rotation powered (‘ordinary’) pulsars (RP PSR).Vertical error bars show the estimated uncertainties on bolometric thermalluminosities, as seen by a distant observer. Horizontal error bars show theestimated age intervals, whenever available; otherwise horizontal arrows markthe characteristic ages (which are usually, although not always, larger thanthe true ages). The position of J1957 in this figure corresponds to the best-fitting model nsmdip 1, while the error bar embraces the models nsmdip 1–3in Table 1.
The theoretical cooling curves in Fig. 10 are calculated using thenumerical code presented by Potekhin & Chabrier (2018). The BSk24model (Pearson et al. 2018) is used for the composition and EoSpof the NS matter. The NSs are supposed to have either non-accreted(ground state) heat blanketing envelopes or accreted envelopes com-posed of helium, carbon, and oxygen up to the densities and tempera-tures where these chemical elements can survive (see, e.g., Potekhinet al. 2003; Potekhin & Chabrier 2018). The accreted envelopes aremore heat-transparent than the ground-state ones. For this reason,the stars with the accreted envelopes are brighter at the early stage oftheir evolution (at 𝑡 (cid:46) yr), but they cool down faster and becomecolder at the late stage ( 𝑡 (cid:38) yr). An envelope may consist ofthe accreted material only partially. In such cases the cooling rate isintermediate between the non-accreted and fully accreted extremesshown in the figure.The critical temperatures for singlet neutron superfluidity in theinner crust and for proton and triplet neutron types of superfluidityin the core of a NS are evaluated, as functions of density, using theMSH, BS, and TTav parametrizations of Ho et al. (2015), which arebased on theoretical models computed, respectively, by Margueron,Sagawa & Hagino (2008), Baldo & Schulze (2007), and Takatsuka& Tamagaki (2004). As can be seen from Fig. 10, thermal luminosityof J1957 is higher than the predictions of the cooling model for theage 𝑡 = 𝑡 c , so that theoretical cooling curves pass to the left of theerror bar in this figure. However, as we mentioned above, it is likely MNRAS , 1–14 (2020) Zyuzin et al.
Figure 11.
Late-time cooling of a strongly magnetized NS, consistent with the spectral fit models nsmdip 1, 2 and 3 (the left, middle and right panels,respectively). The EoS model BSk24 with 𝑀 = . (cid:12) is used for nsmdip 1 and 3, and BSk26 with 𝑀 = (cid:12) is used for nsmdip 2. The polar magnetic field is 𝐵 p = × G (average (cid:104) 𝐵 (cid:105) = . × G) for nsmdip 1 and 2, while for nsmdip 3 we adopt 𝐵 p = . × G ( (cid:104) 𝐵 (cid:105) = . × G). Different line stylescorrespond to different models of baryon superfluidity in the stellar core, as shown in the legend: BS or EEHOr models for the proton pairing gap and TTavor D+av models for the neutron triplet pairing (see text). In each panel, the right/upper line of each type shows the cooling of a NS covered by a non-accretedmagnetized heat blanket made of iron, while the left/lower line of each type corresponds to the fully accreted heat-blanketing envelope composed of layers ofhydrogen, helium, carbon, and oxygen (Potekhin et al. 2003; Potekhin & Chabrier 2018). that the true age 𝑡 is smaller than 𝑡 c . If we treat 𝑡 c as an upper limitto the true age, than thermal luminosity of J1957 is compatible withthe considered theoretical model. The smallest discrepancy betweenthe best-fitting point and the theoretical cooling curves in Fig. 10 isobserved for the model of an NS with 𝑀 ≈ . (cid:12) , covered by thenon-accreted heat-blanketing envelope. In this case, an agreementbetween the model and observations is reached, if we assume that 𝑡 ∼ . ∼ . 𝑡 c .However, the cooling NS models shown in Fig. 10, being non-magnetic, are not fully consistent with the models of strongly mag-netized atmosphere spectra (Appendix A). To produce cooling curvesfully consistent with the spectral fitting, we employ a model of mag-netized heat-blanketing envelope with mass 10 − M (cid:12) . The bottomof such an envelope lies in a deep layer of the outer crust (at densities ∼ g cm − ), which is nearly isothermal at the considered NSages. The interior of a NS is treated as spherically symmetric, buttemperature distribution in the envelope is essentially anisotropic. Atthe magnetic pole, the mechanical and thermal structure of the enve-lope is computed numerically, following Potekhin et al. (2003) andusing updated microphysics input as per Potekhin & Chabrier (2018).Thermal conductivity treatment in the partially degenerate H and Helayers of the accreted envelope has been updated following Blouinet al. (2020). The distribution of the effective temperature over thesurface is then taken consistent with the dipole field, including the ef-fects of General Relativity, as described in Appendix A. The interiorof the star at densities (cid:38) g cm − is treated as spherically sym-metric, with the microphysics (in particular, thermal conductivities,synchrotron neutrino emission rates, etc.) pertinent to the average magnetic field strength for each model. In the cases nsmdip 1 and3 ( 𝑀 = . (cid:12) ), we use the composition and EoS model BSk24,while for nsmdip 2 we use BSk26 (Pearson et al. 2018). The latterEoS is softer and provides the core composition preventing neutronstar from the enhanced cooling at the assumed redshift 𝑧 g = . 𝑀 ≈ (cid:12) .The results are shown in Fig. 11. To test the effects of baryonsuperfluidity we show, in addition to the models with the BS and TTavsuperfluidity of protons and neutrons (as in Fig. 10), also the cooling curves computed using alternative proton and neutron superfluiditymodels in the core. For an alternative proton superfluidity, we usethe EEHOr parametrization by Ho et al. (2015) to the microscopiccalculations of proton critical temperature by Elgarøy et al. (1996),which predicted substantially stronger proton superfluidity than BS.For an alternative superfluidity of neutrons, we use the ‘Av18 SRC+P’model of Ding et al. (2016). which is marked ‘D+av’ in Fig. 11. In thelatter case, the triplet neutron superfluidity is suppressed due to theeffects of many-body correlations (although the maximum criticaltemperatures are similar in the TTav and D+av models, the lattermodel predicts a narrower range of densities where the neutrons aresuperfluid).We see that the obtained estimates of thermal luminosity of J1957can constrain the theoretical NS cooling models, constructed consis-tently with the fitting. Assuming that the true age of the pulsar doesnot exceed its characteristic age, we can conclude that the NS modelnsmdip 1 is compatible with all employed models of superfluidityand envelope composition, but the NS models nsmdip 2 and 3 arehardly compatible with the suppressed neutron superfluidity (D+av),if the envelope is non-accreted. If we assume that the true age isclose to 𝑡 c , then for each spectral model we can select the best-fittingsuperfluidity and heat-blanketing envelope models in Fig. 11, forwhich the cooling curves are in a good agreement with the spectralfitting results. This reinforces the suggestion that the thermal-likepart of the X-ray radiation of J1957 comes from the entire surfaceand is powered by passive cooling. Large uncertainties of the thermalluminosity provided by the BB model do not allow us to constraincooling models. Using the
XMM-Newton and
Chandra observations of the middle-aged 𝛾 -ray pulsar J1957, we detected, for the first time, the ther-mal spectral component and X-ray pulsations. We performed self-consistent modelling of thermal spectrum and cooling for models ofNSs with strong dipole magnetic field. We applied it to the obser- MNRAS000
Chandra observations of the middle-aged 𝛾 -ray pulsar J1957, we detected, for the first time, the ther-mal spectral component and X-ray pulsations. We performed self-consistent modelling of thermal spectrum and cooling for models ofNSs with strong dipole magnetic field. We applied it to the obser- MNRAS000 , 1–14 (2020) iddle aged 𝛾 -ray pulsar J1957+5033 in X-rays vations and estimated the J1957 thermal luminosity. It is consistentwith the NS cooling theory and provides certain constraints on the NSmodel parameters. We found that J1957 is one of the coldest middle-aged NSs with measured effective temperatures: its redshifted ther-mal luminosity is 𝐿 ∞ ≈ ( . − . ) × erg s − . This indicatesthat J1957 has already passed from a relatively slow neutrino coolingstage to a significantly faster photon stage. The X-ray pulse-profilewith a single pulse per period and with the pulse fraction of 18 ± ◦ (cid:46) 𝛼 + 𝜁 (cid:46) ◦ .Using the interstellar extinction-distance relation, we estimatedthe distance to the pulsar to be of ≈ ∼ ACKNOWLEDGEMENTS
We would like to thank the anonymous referee for useful com-ments and A. A. Danilenko for helpful discussion. The work waspartially supported by the Russian Foundation for Basic Research(RFBR) according to the project 19-52-12013. VFS thanks DeutscheForschungsgemeinschaft (DFG) for financial support (grant WE1312/53-1). His work was also partially funded by the subsidy 0671-2020-0052 allocated to Kazan Federal University for the state assign-ment in the sphere of scientific activities. DAZ thanks Pirinem Schoolof Theoretical Physics for hospitality. The scientific results reportedin this article are based on observations obtained with
XMM-Newton ,an ESA science mission with instruments and contributions directlyfunded by ESA Member States and NASA.
DATA AVAILABILITY
XMM-Newton data and https://cxc.harvard.edu/cda/ for
Chandra data.
REFERENCES
Abdo A. A., et al., 2013, ApJS, 208, 17Arnaud K., Gordon C., Dorman B., 2018, Xspec: an X-Ray spectral fittingpackage. Users’ Guide for version 12.10.1. https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XspecManual.html
Auchettl K., et al., 2015, ApJ, 802, 68Baldo M., Schulze H. J., 2007, Phys. Rev. C, 75, 025802Beloborodov A. M., 2002, ApJ, 566, L85Blouin S., Shaffer N. R., Saumon D., Starrett C. E., 2020, ApJ, 899, 46Brazier K. T. S., 1994, MNRAS, 268, 709Brownsberger S., Romani R. W., 2014, ApJ, 784, 154Buccheri R., et al., 1983, A&A, 128, 245Capitanio L., Lallement R., Vergely J. L., Elyajouri M., Monreal-Ibero A.,2017, A&A, 606, A65 Cash W., 1979, ApJ, 228, 939Chang C., Pavlov G. G., Kargaltsev O., Shibanov Y. A., 2012, ApJ, 744, 81Cutri R. M., et al. 2014, VizieR Online Data Catalog, p. II/328Ding D., Rios A., Dussan H., Dickhoff W. H., Witte S. J., Carbone A., PollsA., 2016, Phys. Rev. C, 94, 025802Elgarøy Ø., Engvik L., Hjorth-Jensen M., Osnes E., 1996, Phys. Rev. Lett.,77, 1428Flewelling H. A., et al., 2016, arXiv e-prints, p. arXiv:1612.05243Foight D. R., Güver T., Özel F., Slane P. O., 2016, ApJ, 826, 66Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125,306Ginzburg V. L., Ozernoy L. M., 1965, Sov. Phys. JETP, 20, 689Goodman J., Weare J., 2010, Communications in Applied Mathematics andComputational Science, 5, 65Greenstein G., Hartke G. J., 1983, ApJ, 271, 283Ho W. C. G., Potekhin A. Y., Chabrier G., 2008, ApJS, 178, 102Ho W. C. G., Elshamouty K. G., Heinke C. O., Potekhin A. Y., 2015, Phys.Rev. C, 91, 015806Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974Kargaltsev O., Pavlov G. G., 2008, in Bassa C., Wang Z., Cumming A., KaspiV. M., eds, American Institute of Physics Conference Series Vol. 983, 40Years of Pulsars: Millisecond Pulsars, Magnetars and More. pp 171–185( arXiv:0801.2602 ), doi:10.1063/1.2900138Kargaltsev O., Pavlov G. G., Klingler N., Rangelov B., 2017a, J. PlasmaPhys., 83, 635830501Kargaltsev O., Pavlov G. G., Klingler N., Rangelov B., 2017b, Journal ofPlasma Physics, 83, 635830501Kerr M., Ray P. S., Johnston S., Shannon R. M., Camilo F., 2015, ApJ, 814,128Lallement R., Vergely J. L., Valette B., Puspitarini L., Eyer L., CasagrandeL., 2014, A&A, 561, A91Lallement R., et al., 2018, A&A, 616, A132Manchester R. N., Taylor J. H. J., 1977, Pulsars. W. H. Freeman & Co., SanFranciscoMarelli M., et al., 2013, ApJ, 765, 36Marelli M., Mignani R. P., De Luca A., Saz Parkinson P. M., Salvetti D., DenHartog P. R., Wolff M. T., 2015, ApJ, 802, 78Margueron J., Sagawa H., Hagino K., 2008, Phys. Rev. C, 77, 054309Mignani R. P., Pavlov G. G., Kargaltsev O., 2010, ApJ, 720, 1635Misner C. W., Thorne K. S., Wheeler J. A., 1973, Gravitation. W. H. Freemanand Co., San FranciscoMüller B., 2020, Living Reviews in Computational Astrophysics, 6, 3Pavlov G. G., Zavlin V. E., 2000a, in Martens P. C. H., Tsuruta S., Weber M. A.,eds, IAU Symposium Vol. 195, Highly Energetic Physical Processes andMechanisms for Emission from Astrophysical Plasmas. p. 103Pavlov G. G., Zavlin V. E., 2000b, ApJ, 529, 1011Pearson J. M., Chamel N., Potekhin A. Y., Fantina A. F., Ducoin C., DuttaA. K., Goriely S., 2018, MNRAS, 481, 2994Pechenick K. R., Ftaclas C., Cohen J. M., 1983, ApJ, 274, 846Pierbattista M., Harding A. K., Grenier I. A., Johnson T. J., Caraveo P. A.,Kerr M., Gonthier P. L., 2015, A&A, 575, A3Potekhin A. Y., Chabrier G., 2003, ApJ, 585, 955Potekhin A. Y., Chabrier G., 2018, Astronomy and Astrophysics, 609, A74Potekhin A. Y., Yakovlev D. G., Chabrier G., Gnedin O. Y., 2003, ApJ, 594,404Potekhin A. Y., Lai D., Chabrier G., Ho W. C. G., 2004, ApJ, 612, 1034Potekhin A. Y., Chabrier G., Ho W. C. G., 2014, A&A, 572, A69Potekhin A. Y., Zyuzin D. A., Yakovlev D. G., Beznogov M. V., ShibanovY. A., 2020, MNRAS, 496, 5052Poutanen J., Beloborodov A. M., 2006, MNRAS, 373, 836Reynolds S. P., Pavlov G. G., Kargaltsev O., Klingler N., Renaud M.,Mereghetti S., 2017, Space Sci. Rev., 207, 175Saz Parkinson P. M., et al., 2010, ApJ, 725, 571Schlafly E. F., Meisner A. M., Green G. M., 2019, ApJS, 240, 30Spitkovsky A., 2006, ApJ, 648, L51Sturrock P. A., 1971, ApJ, 164, 529Suleimanov V., Potekhin A. Y., Werner K., 2009, A&A, 500, 891Takatsuka T., Tamagaki R., 2004, Prog. Theor. Phys., 112, 37MNRAS , 1–14 (2020) Zyuzin et al.
Taverna R., Turolla R., Suleimanov V., Potekhin A. Y., Zane S., 2020, MN-RAS, 492, 5057Thorne K. S., 1977, ApJ, 212, 825Willatt R., Ehle M., 2016, Guide for use of the images script.
Wilms J., Allen A., McCray R., 2000, ApJ, 542, 914Yakovlev D. G., Pethick C. J., 2004, ARA&A, 42, 169Yakovlev D. G., Gnedin O. Y., Gusakov M. E., Kaminker A. D., LevenfishK. P., Potekhin A. Y., 2005, Nuclear Phys. A, 752, 590Zavlin V. E., 2009, in Becker W., ed., Astrophysics and Space Sci. Library,Vol. 357, Neutron Stars and Pulsars. Springer, Berlin, p. 181Zharikov S., Mignani R. P., 2013, MNRAS, 435, 2227Zharikov S., Shibanov Y., Komarova V., 2006, Advances in Space Research,37, 1979
APPENDIX A: ATMOSPHERE MODELS FOR NEUTRONSTARS WITH DIPOLE MAGNETIC FIELDS
Magnetized plane-parallel NS atmosphere models are computedusing an advanced version of the code described in Suleimanov,Potekhin & Werner (2009). The code has been modified to accountfor different inclinations 𝜃 𝐵 of the magnetic field with respect to thelocal surface normal. Hydrogen composition is considered, takinginto account incomplete ionization at relatively low temperatures.The effects of the strong magnetic field and the atomic thermal mo-tion across the field on the plasma opacities are treated followingPotekhin & Chabrier (2003) with the improvements described inPotekhin, Chabrier & Ho (2014). Polarization vectors and opacitiesof normal electromagnetic modes are calculated as in Potekhin et al.(2004).Physical models of emission from NSs should take into accountmagnetic field and temperature distributions over the surface. Weassume a dipolar magnetic field, after accounting for the effect ofGeneral Relativity, according to Ginzburg & Ozernoy (1965) (seealso Pavlov & Zavlin 2000b): 𝐵 = 𝐵 p √︃ cos 𝛾 + 𝑓 sin 𝛾 / , cos 𝜃 𝐵 = ( 𝐵 p / 𝐵 ) cos 𝛾, (A1)where 𝐵 is the field strength, 𝜃 𝐵 is the field inclination to the surfacenormal at a magnetic colatitude 𝛾 , 𝑓 = √ − 𝑢 𝑢 − 𝑢 − ( − 𝑢 ) ln ( − 𝑢 ) 𝑢 + 𝑢 + ( − 𝑢 ) , (A2) 𝑢 = 𝑟 g / 𝑅 is the compactness parameter, 𝑟 g = 𝐺 𝑀 / 𝑐 is the gravi-tational radius, 𝐺 is the gravitational constant and 𝑐 is the speed oflight. The distribution of local effective temperature 𝑇 s over the stellarsurface is calculated using the results of Potekhin et al. (2003). In or-der to minimize model dependence, we assume the 𝑇 s -distribution ofan iron heat-blanketing envelope. This assumption does not changeour results since, for any chemical composition of the envelope, thedependence of 𝑇 s on 𝜃 𝐵 is similar to that given in Greenstein &Hartke (1983).Model atmospheres for an inclined magnetic field require solvingthe transfer problem in two dimensions. The optical properties of themagnetized plasma depend on the angle 𝜂 between the photon wavevector 𝒌 and the local magnetic field 𝑩 . On the other hand, under theplane-parallel approximation the radiation field naturally depends onthe angles ( 𝜃 𝑘 , 𝜙 𝑘 ) , where 𝜃 𝑘 is the angle between 𝒌 and the surfacenormal and 𝜙 𝑘 is the angle between the projections of 𝒌 and 𝑩 onthe surface. To avoid interpolation of the opacities over such a two-dimensional grid, the code solves the transfer problem over an ( 𝜂, 𝜓 ) angular grid, where 𝜓 is the azimuth associated to the polar angle 𝜂 , and then the transformation between the angular coordinates ( 𝜂, 𝜓 ) and ( 𝜃 𝑘 , 𝜙 𝑘 ) is used (see Taverna et al. 2020 for more details):cos 𝜂 = sin 𝜃 𝐵 sin 𝜃 𝑘 cos 𝜙 𝑘 + cos 𝜃 𝐵 cos 𝜙 𝑘 , (A3)cos 𝜓 = cos 𝜃 𝑘 − cos 𝜂 cos 𝜃 𝐵 | sin 𝜂 sin 𝜃 𝐵 | . (A4)The photon wave-vector at infinity 𝒌 (cid:48) differs from 𝒌 at the surfacedue to gravitational redshift and light-bending (Pechenick, Ftaclas &Cohen 1983; Pavlov & Zavlin 2000b). The photon energy at infinityis smaller than the photon energy at the surface by factor √ − 𝑢 = /( + 𝑧 g ) , where 𝑧 g is gravitational redshift. The period of PSRJ1957+5033 is sufficiently long so that we can assume that the stellarsurface is spherical and neglect the effects of rotation (see, e.g.,Poutanen & Beloborodov 2006 for description of these effects). Letus consider a surface element d 𝑆 = 𝑅 d cos 𝜗 d 𝜑 , where 𝜗 and 𝜑 arethe polar and azimuthal angles. For the relation between 𝜗 and 𝜃 𝑘 we use approximation (Beloborodov 2002)1 − cos 𝜗 = ( − cos 𝜃 𝑘 )/( − 𝑢 ) . (A5)The flux observed from this surface element is proportional to thespecific intensity 𝐼 ∞ 𝐸 and the solid angle d Ω occupied by this elementon the observer’s sky. This solid angle equals (Beloborodov 2002)d Ω = d 𝑆 cos 𝜃 𝑘 𝐷 − 𝑢 d cos 𝜃 𝑘 d cos 𝜃 (cid:48) 𝑘 , (A6)where 𝜃 (cid:48) 𝑘 is the angle between the local normal to the surface andthe photon momentum at infinity. Without loss of generality wecan choose the polar coordinate axis along the line of sight; then 𝜃 (cid:48) 𝑘 = 𝜗 . Since 𝐼 𝐸 / 𝐸 is invariant (Misner, Thorne & Wheeler 1973,Sect. 22.6), the observed specific intensity 𝐼 ∞ 𝐸 is related to the emittedintensity 𝐼 𝐸 by a constant redshift factor 𝐼 ∞ 𝐸 = 𝐼 𝐸 /( + 𝑧 g ) . Thusthe flux isd 𝐹 ∞ 𝐸 = 𝐼 ∞ 𝐸 d Ω = 𝑅 cos 𝜃 𝑘 𝐷 𝐼 𝐸 ( + 𝑧 g ) d 𝜑 d cos 𝜗. (A7)The monochromatic spectral flux density is then computed by in-tegrating the emission from different local patches over the stellarsurface. Making use of equation (A6), we obtain 𝐹 ∞ 𝐸 = 𝑅 𝐷 ( + 𝑧 g ) ∫ 𝜋 d 𝜑 ∫ 𝜋 / 𝐼 𝐸 ( 𝜃 𝑘 , 𝜙 𝑘 ) cos 𝜃 𝑘 sin 𝜃 𝑘 d 𝜃 𝑘 . (A8)For an axisymmetric magnetic field, the angle 𝜙 𝑘 and the magneticcolatitude 𝛾 are determined at every 𝜃 𝑘 and 𝜑 by the relations (cf.Ho et al. 2008)cos 𝛾 = cos 𝜑 sin 𝜗 sin Θ m + cos 𝜗 cos Θ m , (A9)cos 𝜙 𝑘 = cos 𝛾 cos 𝜗 − cos Θ m | sin 𝛾 sin 𝜗 | , (A10)where Θ m is the angle between the magnetic axis and line of sight.The integration in equation (A8) is restricted by those angles 𝜃 𝑘 thatcorrespond to real values of 𝜗 .For each neutron star parameter set, the local plane-parallel at-mosphere models are computed at the magnetic pole and at threemagnetic latitudes between the pole and the equator, according toTable A1. The atmosphere model at the equator needs not to becalculated, because temperature at low latitudes is so low that it This is equivalent to equation (8) of Ho et al. (2008), where we have restoredthe missed factor cos 𝜃 𝑘 .MNRAS000
Magnetized plane-parallel NS atmosphere models are computedusing an advanced version of the code described in Suleimanov,Potekhin & Werner (2009). The code has been modified to accountfor different inclinations 𝜃 𝐵 of the magnetic field with respect to thelocal surface normal. Hydrogen composition is considered, takinginto account incomplete ionization at relatively low temperatures.The effects of the strong magnetic field and the atomic thermal mo-tion across the field on the plasma opacities are treated followingPotekhin & Chabrier (2003) with the improvements described inPotekhin, Chabrier & Ho (2014). Polarization vectors and opacitiesof normal electromagnetic modes are calculated as in Potekhin et al.(2004).Physical models of emission from NSs should take into accountmagnetic field and temperature distributions over the surface. Weassume a dipolar magnetic field, after accounting for the effect ofGeneral Relativity, according to Ginzburg & Ozernoy (1965) (seealso Pavlov & Zavlin 2000b): 𝐵 = 𝐵 p √︃ cos 𝛾 + 𝑓 sin 𝛾 / , cos 𝜃 𝐵 = ( 𝐵 p / 𝐵 ) cos 𝛾, (A1)where 𝐵 is the field strength, 𝜃 𝐵 is the field inclination to the surfacenormal at a magnetic colatitude 𝛾 , 𝑓 = √ − 𝑢 𝑢 − 𝑢 − ( − 𝑢 ) ln ( − 𝑢 ) 𝑢 + 𝑢 + ( − 𝑢 ) , (A2) 𝑢 = 𝑟 g / 𝑅 is the compactness parameter, 𝑟 g = 𝐺 𝑀 / 𝑐 is the gravi-tational radius, 𝐺 is the gravitational constant and 𝑐 is the speed oflight. The distribution of local effective temperature 𝑇 s over the stellarsurface is calculated using the results of Potekhin et al. (2003). In or-der to minimize model dependence, we assume the 𝑇 s -distribution ofan iron heat-blanketing envelope. This assumption does not changeour results since, for any chemical composition of the envelope, thedependence of 𝑇 s on 𝜃 𝐵 is similar to that given in Greenstein &Hartke (1983).Model atmospheres for an inclined magnetic field require solvingthe transfer problem in two dimensions. The optical properties of themagnetized plasma depend on the angle 𝜂 between the photon wavevector 𝒌 and the local magnetic field 𝑩 . On the other hand, under theplane-parallel approximation the radiation field naturally depends onthe angles ( 𝜃 𝑘 , 𝜙 𝑘 ) , where 𝜃 𝑘 is the angle between 𝒌 and the surfacenormal and 𝜙 𝑘 is the angle between the projections of 𝒌 and 𝑩 onthe surface. To avoid interpolation of the opacities over such a two-dimensional grid, the code solves the transfer problem over an ( 𝜂, 𝜓 ) angular grid, where 𝜓 is the azimuth associated to the polar angle 𝜂 , and then the transformation between the angular coordinates ( 𝜂, 𝜓 ) and ( 𝜃 𝑘 , 𝜙 𝑘 ) is used (see Taverna et al. 2020 for more details):cos 𝜂 = sin 𝜃 𝐵 sin 𝜃 𝑘 cos 𝜙 𝑘 + cos 𝜃 𝐵 cos 𝜙 𝑘 , (A3)cos 𝜓 = cos 𝜃 𝑘 − cos 𝜂 cos 𝜃 𝐵 | sin 𝜂 sin 𝜃 𝐵 | . (A4)The photon wave-vector at infinity 𝒌 (cid:48) differs from 𝒌 at the surfacedue to gravitational redshift and light-bending (Pechenick, Ftaclas &Cohen 1983; Pavlov & Zavlin 2000b). The photon energy at infinityis smaller than the photon energy at the surface by factor √ − 𝑢 = /( + 𝑧 g ) , where 𝑧 g is gravitational redshift. The period of PSRJ1957+5033 is sufficiently long so that we can assume that the stellarsurface is spherical and neglect the effects of rotation (see, e.g.,Poutanen & Beloborodov 2006 for description of these effects). Letus consider a surface element d 𝑆 = 𝑅 d cos 𝜗 d 𝜑 , where 𝜗 and 𝜑 arethe polar and azimuthal angles. For the relation between 𝜗 and 𝜃 𝑘 we use approximation (Beloborodov 2002)1 − cos 𝜗 = ( − cos 𝜃 𝑘 )/( − 𝑢 ) . (A5)The flux observed from this surface element is proportional to thespecific intensity 𝐼 ∞ 𝐸 and the solid angle d Ω occupied by this elementon the observer’s sky. This solid angle equals (Beloborodov 2002)d Ω = d 𝑆 cos 𝜃 𝑘 𝐷 − 𝑢 d cos 𝜃 𝑘 d cos 𝜃 (cid:48) 𝑘 , (A6)where 𝜃 (cid:48) 𝑘 is the angle between the local normal to the surface andthe photon momentum at infinity. Without loss of generality wecan choose the polar coordinate axis along the line of sight; then 𝜃 (cid:48) 𝑘 = 𝜗 . Since 𝐼 𝐸 / 𝐸 is invariant (Misner, Thorne & Wheeler 1973,Sect. 22.6), the observed specific intensity 𝐼 ∞ 𝐸 is related to the emittedintensity 𝐼 𝐸 by a constant redshift factor 𝐼 ∞ 𝐸 = 𝐼 𝐸 /( + 𝑧 g ) . Thusthe flux isd 𝐹 ∞ 𝐸 = 𝐼 ∞ 𝐸 d Ω = 𝑅 cos 𝜃 𝑘 𝐷 𝐼 𝐸 ( + 𝑧 g ) d 𝜑 d cos 𝜗. (A7)The monochromatic spectral flux density is then computed by in-tegrating the emission from different local patches over the stellarsurface. Making use of equation (A6), we obtain 𝐹 ∞ 𝐸 = 𝑅 𝐷 ( + 𝑧 g ) ∫ 𝜋 d 𝜑 ∫ 𝜋 / 𝐼 𝐸 ( 𝜃 𝑘 , 𝜙 𝑘 ) cos 𝜃 𝑘 sin 𝜃 𝑘 d 𝜃 𝑘 . (A8)For an axisymmetric magnetic field, the angle 𝜙 𝑘 and the magneticcolatitude 𝛾 are determined at every 𝜃 𝑘 and 𝜑 by the relations (cf.Ho et al. 2008)cos 𝛾 = cos 𝜑 sin 𝜗 sin Θ m + cos 𝜗 cos Θ m , (A9)cos 𝜙 𝑘 = cos 𝛾 cos 𝜗 − cos Θ m | sin 𝛾 sin 𝜗 | , (A10)where Θ m is the angle between the magnetic axis and line of sight.The integration in equation (A8) is restricted by those angles 𝜃 𝑘 thatcorrespond to real values of 𝜗 .For each neutron star parameter set, the local plane-parallel at-mosphere models are computed at the magnetic pole and at threemagnetic latitudes between the pole and the equator, according toTable A1. The atmosphere model at the equator needs not to becalculated, because temperature at low latitudes is so low that it This is equivalent to equation (8) of Ho et al. (2008), where we have restoredthe missed factor cos 𝜃 𝑘 .MNRAS000 , 1–14 (2020) iddle aged 𝛾 -ray pulsar J1957+5033 in X-rays Figure A1.
Unabsorbed redshifted spectral flux density as function of redshifted energy, calculated according to equation (A8) for models 1, 2 and 3 listed inTable A1 (left, middle and right panels, respectively) for different inclinations Θ m of magnetic dipole axis (drawn with different line styles, according to thelegend) and redshifted effective temperatures (lg 𝑇 ∞ (K) = 5.1, 5.2, 5.3, 5.4, 5.5, 5.6 and 5.7, from bottom to top). For comparison, the BB spectra at the sametemperatures are shown by dotted lines. Table A1.
Parameter sets used in calculation of atmosphere models: cosinesof magnetic colatitude 𝛾 and magnetic field inclination 𝜃 𝐵 , logarithms ofmagnetic field strength 𝐵 and local effective temperatures 𝑇 s .lg 𝑇 ∞ (K) 5.1 5.2 5.3 5.4 5.5 5.6 5.7cos 𝛾 cos 𝜃 𝐵 lg 𝐵 (G) lg 𝑇 s (K)nsmdip 1: 𝑀 = . (cid:12) , 𝑅 = . 𝑧 g = . 𝐵 p = × G0.25 0.424 12.25 5.075 5.183 5.285 5.383 5.483 5.583 5.6840.50 0.724 12.32 5.188 5.292 5.393 5.492 5.591 5.691 5.7910.75 0.900 12.40 5.238 5.337 5.437 5.537 5.637 5.737 5.8371 1 12.48 5.266 5.361 5.459 5.560 5.661 5.762 5.861nsmdip 2: 𝑀 = (cid:12) , 𝑅 = . 𝑧 g = . 𝐵 p = × G0.25 0.398 12.28 5.139 5.246 5.348 5.447 5.546 5.647 5.7480.44 0.633 12.32 5.236 5.340 5.442 5.541 5.640 5.740 5.8400.73 0.875 12.40 5.308 5.408 5.508 5.608 5.708 5.807 5.9071 1 12.48 5.341 5.437 5.535 5.636 5.737 5.837 5.937nsmdip 3: 𝑀 = . (cid:12) , 𝑅 = . 𝑧 g = . 𝐵 p = . × G0.19 0.335 11.80 5.037 5.137 5.236 5.337 5.439 5.542 5.6450.56 0.775 11.90 5.206 5.306 5.405 5.505 5.605 5.705 5.8050.87 0.953 12.00 5.249 5.349 5.450 5.550 5.650 5.749 5.8481 1 12.04 5.260 5.359 5.461 5.561 5.661 5.760 5.859 cannot noticeably affect the observed flux (in practice, we use theblackbody spectrum at the equator, but we have checked that withalternative models 𝐹 ∞ 𝐸 remains the same within 3 per cents). At eachfixed magnetic colatitude 𝛾 , the opacities, polarizabilities, and EoSof the hydrogen plasma were computed on a fixed grid of plasmatemperature and density, from which the values required during theradiative-transfer calculation were obtained by interpolation.To calculate the integral (A8), the specific intensity 𝐼 𝐸 is evaluatedat arbitrary 𝛾 , 𝐸 , 𝜃 𝑘 and 𝜙 𝑘 from the computed values by interpola-tion. It should be noted that 𝐵 -dependent absorption features wouldproduce series of lines, if we kept 𝐸 fixed during this interpolation.In reality, such features are broadened due to the smooth variation of 𝐵 with 𝛾 . In order to reproduce this broadening and thus get rid ofthe non-physical series of lines, we first remap our calculated 𝐼 𝐸 as afunction of ratio 𝐸 / 𝐵 and then interpolate it in 𝜃 𝑘 , 𝜙 𝑘 , and 𝐵 ( 𝛾 ) for every fixed 𝐸 / 𝐵 (cf. Ho et al. 2008). The result of such integrationis shown in Fig. A1.The local effective temperature 𝑇 s and the global effective temper-ature 𝑇 eff are defined by the Stefan-Boltzmann law 𝜎 SB 𝑇 = 𝐹 𝑟 , 𝜋𝑅 𝜎 SB 𝑇 = 𝐿 𝑟 (A11)where 𝐹 𝑟 = ∫ ∞ d 𝐸 ∫ 𝜋 d 𝜙 𝑘 ∫ 𝜋 / 𝐼 𝐸 ( 𝜃 𝑘 , 𝜙 𝑘 ) cos 𝜃 𝑘 sin 𝜃 𝑘 d 𝜃 𝑘 (A12)is the local flux density, which depends on the magnetic colatitude 𝛾 , and 𝐿 𝑟 = 𝑅 ∫ 𝜋 d 𝜑 ∫ 𝜋 𝐹 𝑟 ( 𝛾 ) sin 𝜗 d 𝜗 (A13)is the local bolometric luminosity. The redshifted (‘apparent’) lu-minosity, effective temperature and radius as detected by a distantobserver are (e.g., Thorne 1977) 𝐿 ∞ = 𝐿 𝑟 /( + 𝑧 g ) = 𝜋𝜎 SB ( 𝑇 ∞ ) ( 𝑅 ∞ ) , (A14) 𝑇 ∞ = 𝑇 eff /( + 𝑧 g ) , 𝑅 ∞ = 𝑅 ( + 𝑧 g ) . (A15)For each model (nsmdip 1, 2 or 3), we fix radius 𝑅 , redshift 𝑧 g and 𝐵 p to the values indicated in Table A1 and treat the effectivetemperature 𝑇 ∞ , magnetic axis inclination Θ m and distance 𝐷 ascontinuous adjustable parameters to fit the observed spectral fluxesusing calculated grids of 𝐹 ∞ 𝐸 .In the axisymmetric model, the pulsar geometry is determined bythe angles 𝛼 and 𝜁 that the spin axis makes with the magnetic axisand with the line of sight, respectively (e.g., Pavlov & Zavlin 2000b).To produce phase-resolved spectra, it is sufficient to calculatecos Θ m = sin 𝜁 sin 𝛼 cos 𝜙 + cos 𝛼 cos 𝜁 (A16)for each rotation phase 𝜙 . Some light curves computed by integrationof such phase-resolved spectra are shown in Fig. 9. MNRAS , 1–14 (2020) Zyuzin et al. . . . . . . . . D [ kpc ] . . . . . . . N H [ c m − ] Figure B1.
The relation between the interstellar absorption 𝑁 H and the dis-tance 𝐷 (the dark gray line) with uncertainties shown in light gray in thedirection towards J1957. The dashed black line shows the maximum 𝑁 H de-rived from the spectral analysis of two AGNi in the pulsar field (see text fordetails). APPENDIX B: INTERSTELLARABSORPTION–DISTANCE RELATION
We estimated the distance to J1957 including the interstellarabsorption–distance relation shown in Fig. B1 as a prior in thefitting procedure. The relation was derived in the following way.We used the 3D map of the local interstellar medium presented in https://stilism.obspm.fr/ (see Lallement et al. 2014; Capi-tanio et al. 2017; Lallement et al. 2018 for details) to obtain therelation between the interstellar extinction 𝐸 ( 𝐵 − 𝑉 ) and the distance 𝐷 in the pulsar direction. Then 𝐸 ( 𝐵 − 𝑉 ) was converted to the X-rayabsorbing column density 𝑁 H utilising the relation by Foight et al.(2016). We used linear interpolation to derive the 𝑁 H – 𝐷 dependencebetween the obtained points.We also independently estimated the maximum absorptionin the pulsar direction using two brightest extra-galactic X-ray sources in the J1957 field with coordinates R.A., Dec. =(19 h m . s ◦ (cid:48) . (cid:48)(cid:48) h m . s ◦ (cid:48) . (cid:48)(cid:48) XMM-Newton . According to their spectralenergy distributions, the sources are active galactic nuclei (AGNi).We extracted their X-ray spectra, grouped them to ensure at least 25counts per energy bin and fitted with the absorbed model for AGNoptxagn. The resulting column density 𝑁 H shown by the dashedblack line in Fig. B1 is about 10 cm − which is in agreement withthe value obtained from the 𝐸 ( 𝐵 − 𝑉 ) – 𝐷 relation.The spectral models which we used to describe the pulsar thermalemission includes the ratio of the emitting area radius and the distanceas a parameter. Implementation of the 𝑁 H – 𝐷 relation allowed us toseparate it into two independent parameters. This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000