Multiwavelength Pulsations and Surface Temperature Distribution in the Middle-aged Pulsar B1055–52

We present a detailed study of the X-ray emission from PSR B1055–52 using XMM-Newton observations from 2019 and 2000. The phase-integrated X-ray emission from this pulsar is poorly described by existing models of neutron star atmospheres. Instead, we confirm that, similar to other middle-aged pulsars, the best-fitting spectral model consists of two blackbody components, with substantially different temperatures and emitting areas, and a nonthermal component characterized by a power law. Our phase-resolved X-ray spectral analysis using this three-component model reveals variations in the thermal emission parameters with the pulsar’s rotational phase. These variations suggest a nonuniform temperature distribution across the neutron star’s surface, including the cold thermal component and probable hot spot(s). Such a temperature distribution can be caused by external and internal heating processes, likely a combination thereof. We observe very high pulse fractions, 60%–80% in the 0.7–1.5 keV range, dominated by the hot blackbody component. This could be related to temperature nonuniformity and potential beaming effects in an atmosphere. We find indication of a second hot spot that appears at lower energies (0.15–0.3 keV) than the first hot spot (0.5–1.5 keV) in the X-ray light curves and is offset by about half a rotation period. This finding aligns with the nearly orthogonal rotator geometry suggested by radio observations of this interpulse pulsar. If the hot spots are associated with polar caps, a possible explanation for their temperature asymmetry could be an offset magnetic dipole and/or an additional toroidal magnetic field component in the neutron star crust.

1. INTRODUCTION PSR B1055-52 (B1055 hereafter, also known as PSR J1057-5226) is a middle-aged isolated neutron star (NS) sharing common observational properties with Geminga and PSR B0656+14.These three pulsars were dubbed the "Three Musketeers" by Becker & Truemper (1997) due to having similar spin-down energy loss rates, inferred surface dipole magnetic fields, distances, γ-ray emission and X-ray spectral properties.The X-ray emission of these pulsars includes thermal and nonthermal Corresponding author: Armin Vahdat mv.armin@gmail.comcomponents.It is commonly accepted that most of the thermal radiation emanating from the bulk surface of these pulsars results from the heat transfer from the NS interior, hot spots are caused either by anisotropic heat transfer or returning currents from the magnetosphere or both, while the nonthermal emission is due to synchrotron radiation in the magnetosphere (e.g., Pavlov et al. 2002;Harding & Muslimov 1998).
The source was discovered as a radio pulsar by Vaughan & Large (1972).The X-ray emission was discovered with the Einstein observatory (Cheng & Helfand 1983).The X-ray pulsations were detected with ROSAT by Oegelman & Finley (1993), who suggested that the X-ray spectrum consists of two components, including soft thermal emission from the NS surface.The results were confirmed with ASCA observations (Becker & Truemper 1997, and references therein).
Chandra observations of B1055 revealed that its Xray spectrum, like those of the other two Musketeers, is best described by a three-component model, comprising cold and hot blackbodies as well as a power-law component (Pavlov et al. 2002).XMM-Newton observations of the Three Musketeers allowed for a phase-resolved, lowresolution X-ray spectroscopy which highlighted some differences between the three Musketeers in surface temperature distribution, and orientations of the rotation and magnetic axes (De Luca et al. 2005).There are also differences regarding the X-ray detection of pulsar wind nebulae (PWNe) around these close (≲ 500 pc) NSs.B1055 and PSR B0656+14 have surprisingly faint PWNe while the three-tail Geminga's PWN is relatively bright (e.g., Caraveo et al. 2003;De Luca et al. 2006;Posselt et al. 2015;Bîrzan et al. 2016;Posselt et al. 2017).
In radio, B1055 exhibits a unique emission profile, featuring a main pulse (MP) and an interpulse (IP) separated by approximately 160 • in pulse longitude (Weltevrede & Wright 2009).These components likely originate from distinct magnetic poles.The MP-IP separation remains consistent across different radio frequencies, suggesting a highly inclined magnetic axis relative to the pulsar's rotation axis (Weltevrede & Wright 2009, and references therein).
In the γ-ray regime, B1055 exhibits an unusual profile, as emphasized in the Third Fermi-LAT pulsar catalog (3PC; Smith et al. 2023).It has three overlapping principal peaks with indications of two small peaks preceding and following the outer steep wing of the principal γ-ray emission components.Only two of the 150 young pulsars in the 3PC show such a complex profile, and the origin of these profiles is not understood.
The distance to B1055 is still a matter of debate.The ATNF Pulsar Catalogue (ver.1.68) reports 0.093 kpc, based on the DM-based value (DM = 29.69±0.01,Petroff et al. 2013) and the electron density model of Yao et al. (2017), while a distance of 0.714 kpc is derived using the NE2001 model (Cordes & Lazio 2002).Mignani et al. (2010) suggested a much smaller distance of 0.2-0.5 kpc by investigating the contributions of the individual thermal components to the multiwavelength spectrum.This distance can also explain the unusually high γ-ray efficiency of the pulsar that has been previously obtained assuming a distance of 750 pc (Abdo et al. 2010).Throughout this paper, we will scale the distance to 0.35 kpc to calculate the luminosity and other related quantities.
In this paper we report the results of two new XMM-Newton observations of B1055 carried out in 2019 (PI B. Posselt).Posselt et al. (2023a) analyzed the phaseintegrated data for a long-term flux variability, but did not find significant differences between these new and the 2000 XMM-Newton data.However, the new data allow for a better characterization of the spectral and timing properties of the source and tighter constraints on the geometry of the pulsar.
The paper is structured as follows.In Section 2, we describe the data and our data reduction procedure.Section 3 is dedicated to phase-integrated spectroscopy, for which we combine the previous 2000 XMM-Newton data with our new observations to better constrain temperatures and test various atmosphere models.Timing analysis results are reported in Section 4, where we discuss the pulse profile and its properties such as pulsed fraction and phase-energy map.In Section 5 we focus on phase-resolved spectroscopy of B1055.Lastly, in Section 6 we discuss our results.

X-ray data
The XMM-Newton observatory observed B1055 on 2019 June 21 for 85 ks and 2019 July 9 for 81 ks (Table 1).During both observations, the European Photon Imaging Camera (EPIC) (Strüder et al. 2001) was operated in small-window (SW) mode for EPIC-pn (time resolution 5.7 ms), while EPIC-MOS1 and EPIC-MOS2 cameras (Turner et al. 2001) were operated in full-frame (FF) mode (time resolution 2.6 s).To better constrain the phase-integrated spectral fit parameters (Section 3), we also reanalyzed the previous XMM-Newton observations of B1055 taken on on 2000 December 14 and 15 for 24 ks and 57 ks, respectively (De Luca et al. 2005).
Unfortunately, due to the the collapse of data along the CCD columns, the EPIC-pn data of 2000 obtained in Timing mode, are faced with significantly higher background contamination resulting in a lower signalto-noise ratio (S/N) above a few keV compared to the 2019 observations.
The data processing of the four EPIC observations was done with the XMM-Newton Science Analysis System (SAS) ver.20.0.0 (Gabriel et al. 2004) applying standard tasks.The dead-time corrected net exposure times, corresponding camera modes and filters, and count rates are given in Table 1.

Radio data
We obtained radio pulsar data simultaneously with the XMM-Newton observations.Two 34 ks observations (program ID P1010) were carried out on MJD 58 655 and MJD 58 673 with the CSIRO Parkes 64 m radio telescope (also known as Murriyang).We used the UWL receiver (0.7-4.0 GHz) and the MEDUSA backend (e.g., Hobbs et al. 2020).The observations were excised of radio frequency interference (RFI) and calibrated using standard PSRCHIVE tools (Hotan et al. 2004;van Straten et al. 2012).From the averaged pulse profile, we produced a smoothed, high S/N template.By crosscorrelating this template in the Fourier-domain (Taylor 1992), we generated pulse times of arrival (ToAs) for 57 subintegrations for each Parkes observation.
We also obtained simultaneous and contemporaneous data on B1055 from the MeerTime Thousand Pulsar Array (TPA) programme (Johnston et al. 2020).
The observations were carried out with the 64-dish SARAO MeerKAT radio telescope.The TPA used the L-band receiver (centred at a frequency of 1284 MHz) and a total bandwidth of 775 MHz.Details on the data reduction and RFI removal can be found in Parthasarathy et al. (2021); Lazarus et al. (2016), details on the TPA procedures (such as template production) are reported by Posselt et al. (2023b).We used three observing epochs: MJD 58 655.65 (simultaneous with XMM-Newton , 5.4 ks), MJD 58 678.68 (202 s), MJD 58 689.72 (202 s) to generate ToAs with the appropriate MeerKAT-based template.

γ-ray data
We used the Fermi Large Area Telescope (LAT) γray data set for this pulsar that was provided with the 3rd Fermi-LAT Pulsar Catalog (hereafter 3PC; Smith et al. 2023).The data set consists of photons detected between MJD 54682 and 58791, including the epochs of our X-ray and radio observations, with energies between 50 MeV and 300 GeV.

PHASE-INTEGRATED SPECTRAL ANALYSIS
3.1.Spectral extraction XMM-Newton observations of 2000 December 14-15 were carried out in Timing mode for pn and Full Frame (imaging) mode for MOS1 and MOS2.Following the works of Posselt et al. (2015Posselt et al. ( , 2023a)), we extracted pn spectra covering energies higher than 0.4 keV since soft background noise can impact the extracted spectra below this energy.We used 33 ⩽ RAWX pixel ⩽ 39 to extract the source spectra of both observations (101 and 201 hereafter) whereas 5 ⩽ RAWX ⩽ 7 and 4 ⩽ RAWX ⩽ 6 were used for background regions of 101 and 201, respectively.
The 2000 EPIC-MOS1/2 data were extracted in the same way as described in Posselt et al. (2015) (e.g., circular apertures of 45 ′′ and 60 ′′ radii for source and background regions) The EPIC data of 2019 June 20 and July 19 were collected in imaging modes (see Table 1).For EPIC-pn the source region is a circle with a radius of 30 ′′ .We utilized the SAS tool eregionanalyse to fine-tune the position and the aperture that maximizes the source-tobackground count ratio.For EPIC-MOS1/2 we used the circular aperture of the 45 ′′ radius, similar to the data of 2000, to extract the phase-integrated spectra.As shown in Figure 1, we have selected the background regions with a circular aperture of 60 ′′ radius, sufficiently close to the pulsar to get an accurate estimate of the nearby background contribution but far enough to avoid contamination from the point-spread function of the source.

Spectral fitting
We started our analysis by comparing the two observations of 2019 (101N and 201N hereafter) with each other.We have checked and confirmed the consistency of the two data sets in terms of source and background count ratio and distribution and spectral fit parameters.
We restricted our energy range to 0.3 − 8 keV.Above 8 keV, not only the background contribution exceeds the source flux but also the accuracy of our fit parameters did not improve further with increasing the upper energy bound.
We used the Tübingen-Boulder model through its XSPEC implementation tbabs to fit the interstellar absorption by setting the abundance table to wilm (Wilms et al. 2000) and photoelectric cross-section table to vern (Verner et al. 1996).
We first fitted the high energy part (2.3 − 8 keV) of the spectra with an absorbed power-law (PL) model (powerlaw in XSPEC).Although the absorbing hydrogen column density, N H , generally does not affect the spectrum at high energies, we fixed N H and varied N H at different values and checked the confidence contours kT1, K1, R1 = 153.9+3.7 3.8 eV, 37.7    2).Hot BB, cold BB, and PL components are displayed with red, green, and yellow colors.The fit parameters with a margin of 1σ error are also shown at top of the figure . of PL-index (Γ)-normalization for the two epochs and found an agreement within 1σ in each case.
In the next step, we fit the 0.3 − 8 keV broad-band spectra using a three-component model.This model is comprised of two blackbody (BB) models (bbodyrad in XSPEC) and a PL model with a best-fit value of Γ = 1.8, determined from the PL-only fit.We also allowed N H to vary freely.For the two observations of 2019, we compared the confidence contours of temperaturenormalization for the hot and cold BB components and found that their values agreed within a 3σ range.
After establishing consistency, we combined the two data sets of 2019 and compared the combined data with the data of 2000 by repeating the same checks.While we obtained a similar consistency level between the 2000 and 2019 data as before, the count statistic and background noise level were different, which was most noticeable at high energies.
Finally we combined the 12 spectra, and carried out a 2BB+PL fit.We obtained the best fit (χ 2 υ = 1.09;where υ=822 is the number of degrees of freedom) with BB temperatures kT 1 = 153.9±3.7 eV, kT 2 = 68.6±0.8 eV, and N H = (1.9 ± 0.3) × 10 20 cm −2 (Figure 2 and Table 2).We also tested several NS atmosphere models available in XSEPC and applicable to atmospheres with magnetic fields of about 10 12 G: NSATMOS (Heinke et al. 2006), NSA (Pavlov et al. 1995;Zavlin et al. 1996), and NSMAXG (Mori & Ho 2007;Ho et al. 2008).Fitting with single-component and two-component atmosphere models resulted in large systematic residuals at higher energies (χ 2 υ ≳ 8) whereas two-component and threecomponent models (by adding a PL or BB) also produced unacceptable fits as a consequence of large systematic residuals at low energies, ≲ 1.5 keV.

Radio timing solutions
We applied the pulsar-timing code, tempo2 (Hobbs et al. 2006) for our ToAs to determine the timing solutions for either the two Parkes or the three MeerKAT observations.For this, we used reference epoch MJD 58 664.376 (near  tion from Posselt et al. (2023a).From the tempo2 fits, we obtain the frequency, ν, and its time derivative, ν, see Table 3.The difference in ν is 0.13 nHz between the Parkes and MeerKAT timing solutions.This difference is 1.9σ of the MeerKAT timing solution, the difference ν has a similarly low significance level (1.3σ).Hence, the two independent radio timing solutions agree with each other.For comparison with the X-ray data, we will use in the following the Parkes radio timing solution.

Relative γ-ray timing solution
We folded this Fermi-LAT data for this pulsar using the 3PC ephemeris (Smith et al. 2023), which was developed from 12 years of monitoring with Parkes.We replaced the "TZRMJD" and "TZRFRQ" parameters with those from our own radio timing ephemeris to  phase-align the γ-ray pulse profile with our own radio profile, and found a radio/γ-ray alignment very similar to that found in 3PC.

X-ray timing solutions
Since the time resolution of the MOS detectors in FF mode is as large as 2.6 seconds, they are not suitable for the timing analysis of B1055.Therefore, we only use the 101N and 201N pn-observations which have a nominal frame time of 5.7 ms in SW mode.We corrected all event times in our data sets to the solar system barycenter using the standard SAS task barycen.
The first event in 101N observation was detected at MJD 58654.849534barycentric dynamical time; TDB (or t 1 = 677449468.950103s in XMM-Newton mission reference time; MRT).For the second observation, the first event was detected at MJD 58673.018533TBD (or t 2 = 679019270.421713s in MRT).
The times elapsed between the first and last detected events for 101N and 201N data sets were T span1 = 79988.797120s (≈ 0.9258 days) and T span2 = 79135.1205240s (≈ 0.9159 days), respectively.
We used the Z 2 n statistic (Buccheri et al. 1983) to study X-ray pulsations in the vicinity of the frequency expected from the radio timing solution.We optimized GTI screenings, energy regions, and extraction apertures to ensure high signal-to-noise (S/N ) extraction and large Z 2 n,max values.Prior to the ν and ν search in our X-ray data, we subtracted (t 1 + T span2 + t 2 )/2 from the event times of 101N and 201N to minimize the correlation between ν and ν, so that the measured ephemeris corresponds to the middle of the two observations (i.e., T ref = MJD 58664.376681).
We first performed the Z 2 n test separately for the two data sets at fixed ν = −1.478×10−13 Hz s −1 (the Parkes frequency derivative; see Table 3), using the stingray python package (Huppenkothen et al. 2019).For 101N (201N), we extracted 44113 (41984) events from GTIfiltered 0.3 − 7 keV data using the extraction region of a 30 ′′ radius, which provided the largest Z 2 n,max values.From the H-test (de Jager et al. 1989) we found that the H-statistic reaches a maximum for three harmonics (n = 3) for each of the data sets.We obtained Z 2  1,max = 1925  (1802) and Z 2 3,max = 2123 (2029) for 101N (201N) data at the frequency ν = 5.07317468( 13) Hz ( 5.07317498( 14) Hz).The numbers in brackets indicate the 1σ uncertainties calculated as: The frequencies obtained from the two individual Xray observations agree with each other within 2σ, and within 5σ with the frequency measured by the Parkes telescope.
In the next step, we combined the two EPIC-pn data sets (86,097 events in the 0.3 − 7 keV band) and chose the middle of our two observations as the reference time (MJD 58664.376681).For the total time span T span(1+2) ≈ 19.085 d, the combined observation should allow one to measure not only frequency, but also frequency derivative.Therefore, we calculated Z 2 n on a ν-ν grid and found Z 2 1,max = 3739 (Z 2 3,max = 4155) at ν(Z 2 1max ) = 5.073174881(4) Hz and ν(Z 2 1,max ) = +7 × 10 −14 Hz s −1 (see Figure 4).
The Z 2 1,max and Z 2 3,max values in the summed data are very close to the sums of Z 2 1,max and Z 2 3,max in separate data sets (see Table 4), indicating that we reached a good phase connection between the two data sets and obtained a coherent timing solution.However, the difference of 72 nHz between the frequencies measured in the X-ray and Parkes radio data exceeds the formal statistical uncertainty of 4 nHz by a factor of 18.The positive value of ν(Z 2 1,max ) is in conflict with the well constrained, negative radio timing value, likely indicating an issue with the accuracy of the X-ray timing solution.Since the two independent radio timing solutions agree with each other, we explored technical reasons for the deviating XMM-Newton results and possibilities for correction.Firstly, we checked possible time jumps in our two X-ray time-series.Time jumps can occur due to the known temperature effect on the frame time of the EPIC-pn detector as well as 'reset' errors of the counter clock (Kirsch et al. 2004).In 101N, no real-time jumps have been recognized by SAS and in 201N they were successfully corrected (we set the SAS environment parameter SAS JUMP TOLERANCE to a value of 44). 1  The second observation lacks some internal housekeeping data that would allow a detailed monitoring of the detector temperature.Without these additional data, SAS relies on approximations for correcting the temperature-based clock drifts, which, however, usually work very well.Because of the lack of these additional data, we have investigated in particular whether the second observation may cause faulty results.Unrecognized time jumps could in principle occur in ground station switches.Gotthelf & Halpern (2020) had an XMM-Newton observation in the same instrument mode immediately after our second observation.Their results do not indicate any unusual timing offset or timing problems.Instead of the three ground stations in our observation, the data by Gotthelf & Halpern (2020) uses only two.We excluded in our observation those times of the additional ground station, suspecting that it may introduce some problem.However, this did not change the X-ray determined ν (which was different by only 2 nHz from the full X-ray data solution).Thus, ground station switches in the second observation are unlikely to have introduced any timing error.
Comparing radio and X-ray determined pulsation periods of different pulsars in different XMM-Newton instrument modes, Martin-Carrillo et al. ( 2012) reported a relative timing accuracy, △P/P radio of better than 10 −8 , where △P is the difference between the measured period.For our B1055 observations, the relative error is △P/P radio = 1.4 × 10 −8 .This value seems to be reasonably close to expectation.Together with our checks for any possible time jumps, we therefore conclude that we have to add a systematic uncertainty on the order of 60 nHz to our statistical uncertainty of the determined X-ray frequency.
Similarly, the νX has a larger uncertainty.Although the timing solution derived from X-rays is less accurate than the one obtained from radio, using this specific Xray timing solution maximises the pulsed X-ray signal.Since the resulting errors in the photon phases are small (see Section 4.4), we employed the derived values of ν X and νX in the subsequent X-ray analyses to maximise statistical robustness and consistency in X-rays.

Timing properties
We studied the harmonic behavior of the X-ray pulse profile using the Fourier coefficients.For comparison, Table 4. Maximum values of Z 2 n (ν, ν) for different energy ranges and numbers of harmonics (n).N and %bgd are the number of total counts and percentage of the average background counts in the source aperture of the combined and individual (101N & 201N) data.χ 2 υ characterizes the agreement between the binned histograms and the harmonic description of the phase-folded light curve in each band.we also obtained the folded light curve in the form of a histogram in 0.3 − 7 keV with the timing solution determined above.The phase-folded light curve and contributions from each harmonic are displayed in Figure 5.
The appropriate number of harmonics can be different in different energy bands.We used a chi-square test from the Python library SciPy (Virtanen et al. 2020) to compare the sum of 3 harmonics with the 20-bin histogram of the folded light curve.The two descriptions of the pulse profile are in very good agreement.This result implies that binning is not necessary to examine the pulse profile in such a case of smooth pulsations, and additional uncertainties associated with the binning procedure can be avoided using the Fourier analysis.
The normalized X-ray pulse profile F(ϕ) (see Equation (A1) in Hare et al. 2021 for its definition), in the 0.3−7 keV range, is plotted over the Parkes radio profile at 2.4 GHz, and the Fermi-LAT γ-ray light curve from 50 MeV to 300 GeV is shown in Figure 6. Figure 7 highlights the pulse profiles and their 3σ uncertainties in four energy bands together with Parkes radio and Fermi-LAT γ-ray light curves.In order to estimate the intrinsic uncertainty in the alignment of the radio and X-ray profiles, we calculated the maximum phase shift, ∆ϕ max , between the X-ray timing solution (ν X , νX ) and Parkes radio timing solution (ν R , νR ) as follows.
For a given timing solution, the change of phase between times t 1 and t 2 (the first event in the first observation and the last event in the second observation) is: where ∆t = t 2 − t 1 and T ref is the reference time.The difference of phase changes between the radio solution and X-ray solution is: We estimated ∆ϕ max = 0.12, considering the midpoint of each observation as the reference times for the respective solutions, and ∆ϕ max = 0.13, considering the beginning of observation 1 and the end of observation 2 as the reference times.
The pulsed fraction (PF) is one of the important properties of a folded light curve which primarily depends on the phases and amplitudes of Fourier harmonics and their dependence on energy.We estimated PFs for three different PF definitions, namely, the amplitude PF (p amp , sometimes referred to as "peak-to-peak" or "max-to-min" PF), the area PF (p area ), and the rootmean-square PF (p rms ; see Appendix C in Hare et al. 2021 for exact definitions).The energy dependence of the PFs is shown in Figure 8.Both the amplitude and area PFs increase with the photon energy, saturate at ∼ 70%-75% above 1 keV and show a hint of a decrease towards 5 keV.Similarly, p rms increases with energy in the 0.15 to ∼ 1.2 keV range, and then decreases at higher energies.
While the pulse profiles presented in Figure 7 do not explicitly show all the energy ranges discussed below, it is important to note that the analysis includes a broader range of energy bands to explore the pulsar's emission behavior comprehensively.In accordance with Figure 7, the pulse profile at the lower energies, E ≲ 0.4 keV, displays the lowest maxima and minima.
The energy range of 0.6 − 0.8 keV, where the transition between the cold and hot BB components takes place (see Figure 2), exhibits the highest maxima and deepest minima.Similarly, the 1.2-2 keV corresponds to the transition region between the hot BB and PL components, with the maxima aligned in phase, but the minimum at 1.2 − 2 keV occurring approximately 0.2 phases earlier.
A significant phase difference of about 0.5 is observed between the maxima and minima in the 0.45 − 0.6 keV (cold BB), 0.6 − 0.8 keV (transition from cold to hot BB), and 0.8 − 1.2 keV (hot BB) ranges.This suggests the presence of a hotter region on the NS's surface once per period, with the closest approach to the line of sight occurring at phases of 0.15 − 0.2.
Lastly, in the 2 − 5 keV range (PL), the light curve seemingly exhibits two maxima and two minima per period.These maxima occur at phases around 0.2 and 0.7.
There is a noticeable difference in the shape and amplitude of the light curve between 0.15 − 0.3 keV (red color with a maximum amplitude at ∼ F(ϕ) ∼ 1.15) and other bands.It is noteworthy that about half (N = 81103) of the total events (N = 167600) are confined in 0.15 − 0.3 keV.
We also investigated the pulsation behavior through two different representations of the 2D phase-energy space of the X-ray events.Following the work of Arumugasamy et al. ( 2018), we calculated the deviations of the counts, N ij in each phase-energy bin from the phase-averaged value.These deviations from the phaseaveraged values are then represented by the significance map with pixel values: where i and j enumerate the energy and phase intervals, respectively, and N i = J −1 J j=1 N i,j is the phaseaveraged counts in the i th energy bin, and J is the number of phase bins.Here we have binned the events in the phase-energy space, using J = 10 equal-sized phase bins and choosing variable-size energy bins to maintain ≳ 40 counts per phase-energy bin.Note that this phaseenergy deviation map is restricted to the 0.15 − 7.0 keV energy range, where we have enough source counts (Figure 9 left).
Using alternatively a normalization method by Tiengo et al. (2013), the numbers of counts in phase-energy bins, N i,j , are divided first by the phase-averaged counts in the i th energy bin, N i , and then by the energy-averaged (for the 0.15 − 7.0 keV energy band) counts in the j th phase bin, N j = I −1 I i=1 N ij , where I is the number of energy channels.The resulting normalized phase-energy map has the following pixel values: These results are displayed in Figure 9 (right).
Upon examining the right panel (the normalized phase-energy map according to Tiengo et al. 2013), we can classify the counts in the map into four distinct regions based on energy.In the first region (0.15−0.3 keV), a minimum pulsation between phases 0.1 − 0.4 is observed, followed by a faint peak between phases 0.7−0.9.
The second region (0.3 − 0.5 keV) shows a minimal to negligible pulsation across all phases.Interestingly, this contrasts with both the left panel (the phase-energy deviation map) and the pulse profile in Figure 7, where continuous pulsation is more evident throughout the soft to mid X-ray range.
In the fourth region (0.5 − 1.5 keV), a pronounced pulsation maximum is present between phases ϕ = 0.1−0.4,and a pulse minimum between phases ϕ = 0.6 − 1, both of which are corroborated by Figure 7.
The count distribution in the last region (2 − 5 keV) appears dispersed in both the right and left panels, posing a challenge in determining the nature of the observed pulsation.However, it is worth noting that a pulsation minimum and maximum are observed between 2−3 keV.

PHASE-RESOLVED SPECTRAL ANALYSIS
We defined 5 equal-sized phase bins and extracted spectra in each bin using the same extraction radius, filtering, and binning criteria as those in the phaseintegrated analysis (Section 3).We verified that these bins contained enough (at least 15,000) counts to sufficiently constrain the spectrum.Building on the phaseintegrated best 2BB+PL fit, we carry out spectral fits, freezing the absorbing hydrogen column and photon index to their best-fit values, N H = 2 × 10 20 cm −2 and Γ = 1.8.Each spectrum, its best-fit model parameters, and χ 2 υ values are shown in Figure 10.In contrast to the approach by De Luca et al. (2005), where the temperatures were fixed due to a low S/N ratio, our analysis benefits from a higher S/N ratio, which allowed us to vary both the temperatures and their corresponding normalizations in our spectral fit.For each bin, a good fit was obtained, constraining the parameters of the three spectral components throughout all phases.However, in the fourth bin (ϕ = 0.6-0.8), the ftest indicated that the inclusion of the Hot BB component is statistically kT1, K1, R1 = 122.7 +6.1 5.6 eV, 364.4 +165.7 117.1 , 0.67 +0.14 0.12 km kT2, K2, R2 = 62.7 +2.2 2.4 eV, 39.8 +9.6 6.8 E3, 7.0 +0.8 0.6 km , KPL = 1.8, 14.9 +1.not required.
Variations of hot and cold thermal components with rotational phase are demonstrated by the temperaturenormalization confidence contours for 5 phase bins (Figure 11), as well as the phase dependencies of temperature, radius and unabsorbed flux (Figure 12).For the Cold BB component, the phase-integrated spectrum yields the best-fit temperature kT BB,c = 69 eV and an equivalent emitting sphere radius of R BB,c = 0.035K 1/2 BB,c d 350 ∼ 5 km.We see that the best-fit tem-perature varies from 75 eV in Bin 1 to approximately 65 eV in Bins 2 and 3, subsequently returning to a value of 70 eV in Bin 5.The apparent radius exhibits variations anti-correlated with temperature, increasing from 4.4 km in Bin 1 to around 7 km in Bin 2 and then decreasing to 4.8 km in Bin 5.Although the temperaturenormalization contours move almost along the direction of their maximal stretch, caused by the anticorrelation between the temperature and radius, the variation is likely real as demonstrated by the parameter uncertain-   ties in Figure 12.In addition, the 99% confidence contours of Bin 1 and Bin 5 are clearly separated from the respective contours of the other bins in Figure 11.
The variability in the variations of unabsorbed 0.3 − 8 keV flux F unab BB,c exhibits a pattern akin to that of varies approximately in phase with kT BB,c (perhaps with a slight phase shift) because the flux in the Wien tail of a thermal spectrum is particularly sensitive to temperature variations.However, the relative variations of the flux are small, within ∼ 10% of its average value, because the effect of temperature variations is partly compensated by antiphase variations of emitting area, ∝ R 2 .
The hot BB component exhibits more pronounced phase variations.The hot temperature oscillates with phase, with minimum in Bin 2 and maximum in Bin 5, in concert with the cold BB temperature.The variations of radius (and emitting area) are anti-correlated with the temperature variations, similar to the cold BB.In contrast to the cold BB component, the unabsorbed flux associated with the hot BB component (F unab BB,h ) varies in phase with R BB,h (hence with the projected area of the hot region), and in antiphase with hot and cold temperatures and the cold BB flux.Such behavior is consistent with the phase shift between the low-energy and mid-energy light curves (Figure 7).The hot BB flux is substantially lower (by a factor of about 6, on average) than the cold BB flux, and it shows much stronger variations, with a minimum value close to zero, and relative amplitude of about 50%.
We also checked whether allowing Γ to vary would change the variations thermal parameters.The results are also plotted in Figure 12, showing no significant differences between the cases of free and fixed Γ.For both fixed and free Γ, the PL normalization varies in antiphase with the emitting area of the hot BB component (Figure 13).

Discrepancy between X-ray and Radio Timing
Solutions of B1055-52 As described in Section 4.3, we find a relative difference of the measured X-ray to the radio period of △P/P radio = 1.4 × 10 −8 .Although this is in agreement with the relative timing accuracy in different XMM-Newton instrument modes (Martin-Carrillo et al. 2012), the high statistical significance of this difference prompted us to also consider a pulsar glitch as a potential non-technical explanation.
However, the simultaneous radio observations from Parkes, which are sensitive to frequency variations, show no indication of a glitch during the observed period.There is also no indication of a glitch from the three considered epochs of MeerKAT data.Moreover B1055 is not known to exhibit glitching behavior.From a monthly monitoring program with Parkes since the beginning of 2007, Lower et al. (2021) reported only a 90% upper limit on the size of undetected glitches of △ν 90% g /ν < 3.4 × 10 −9 for B1055.This glitch limit is below our X-ray/radio discrepancy.For these reasons, we regard it as highly unlikely that B1055 experienced an undetected glitch between the X-ray observations or after the simultaneous Parkes coverage ended in the second XMM-Newton observation.

The phase-resolved thermal X-ray spectra
One or two BB components are commonly used to describe the thermal spectra of middle-aged pulsars (e.g., Caraveo et al. 2004;De Luca et al. 2005 for the Three Musketeers, Schwope et al. 2022 for PSR B0656+14;Rigoselli et al. 2022 for PSR 1740+1000).However, our phase-resolved spectral analysis shows that such a twotemperature model is not a fully consistent description of the data because both the hot and cold BB component show significant variations with phase.Hence, the temperature is not uniformly distributed, neither in the hot spot nor in the bulk NS surface.Arumugasamy et al. (2018) showed similar phase variations of the cold and possibly the hot BB temperatures for PSR B0656+14 (their Figure 13), although without providing 2D confidence contours for normalization and temperature.
Based on the phase-resolved spectra of B1055 (Figure 10) and the respective parameter evolution over phase (left panels of Figures 11 and 12), the behavior of the hot spot BB ("S1") can be summarised as follows: Starting from Bin 1, which includes the maximum of the 0.5-1.5 keV light curve (where the hot BB component dominates; see Figure 7), the temperature of the hot BB component decreases with increasing phase and reaches a minimum near the boundary between Bin 2 and Bin 3, i.e., on the descending part of the light curve.While the light curve reaches its minimum near the boundary between Bin 3 and Bin 4, the temperature grows up to its maximum at Bin 5, at the ascending part of the light curve.The phase shift between the temperature and the 0.5 − 1.5 keV flux oscillations is caused by the normalization (projected area) varying in opposite direction with respect to the temperature variation.In particular, the temperature increase between Bin 3 and Bin 5 is partly compensated by the normalization decrease between Bin 2 and Bin 5, so that the increase of the flux is substantially smaller than it would be at a constant normalization.
The cold BB component exhibits similar behavior to the hot BB component between Bin 1 and Bin 2, where the temperature starts to decrease and the normalization (radius) begins to increase (right panels of Figures 11 and 12).However, during the transition from Bin 2 to Bin 4, both the kT and the normalization remain nearly constant, within a 1σ range.The cold BB component shows antiphase variations of the temperature and radius, similar to the hot BB component.
The observed behavior of the temperature and visible area is incompatible with the assumption of a uniform temperature in the hot and cold BB component.At this assumption the kT -normalization confidence contours would be shifting along the normalization axis, in contrast to the actual observations (Figure 11,left).
Evidence for at least one other thermal component is also provided by the phase-energy maps.They indicate the presence of another, secondary spot ("S2"), best visible in the right panel of Figure 9, at phases ϕ ∼ 0.2-0.5 and energies E ≲ 0.3 keV.For a centered dipole, one would expect two similar hot spots.For a nearly-orthogonal rotator, these two hot spots should correspond to two maxima in the light curves.In principle, the location of the second hot spot could be such that only a portion of it is visible during a rotational cycle.B1055 is likely a nearly-orthogonal rotator, as indicated by the radio pulse profile with two pulses per period (Figure 7).Therefore, second hot spots would be the natural explanation for the pulsar's S2 emissions.However, the S2's apparent colder temperature (lower energy in the phase-energy map) is puzzling.The S2 is also noticeable in the light curves.In the softer energy band, ∼ 0.3 − 0.5 keV, a second smaller hump, which might be associated with S2, appears.It is aligned (perhaps coincidentally) with the radio inter-pulse.In addition, the phase-dependent cold BB temperature (Figure 12) shows a flattening around the minimum (ϕ = 0.2 − 0.8), possibly due to the contribution from the warm second spot S2.
Taking into account bending of photon trajectories by the NS gravitational field and assuming locally isotropic BB emission emerging from the surface, Turolla & Nobili (2013) calculated light curves and the pulsed fraction, PF amp , as a function of the viewing angle and hot spot position for one and two spots with various angular radii.In particular, they found that for two antipodal spots of equal temperatures and sizes on the surface of a NS with mass M = 1.4M ⊙ and radius R = 15 km, the maximum pulsed fraction, PF amp ≈ 29% for energyintegrated flux, is reached for an orthogonal rotator.
The measured PF in the 0.8 − 1.2 keV range for B1055, corresponding to the hot spot emission, is significantly higher, reaching approximately 65%-70%.But Turolla & Nobili (2013) also showed that much higher PFs (up to 100%) can be obtained if there is only one spot, or two rather different non-antipodal spots.However, this finding alone does not allow us to conclude that B1055 has very different or non-antipodal hotspots because another reason for high PFs can be the presence of a light-element (H or He) NS atmosphere.Emission from such an atmosphere is beamed along the magnetic field (Pavlov et al. 1994), and we can expect a higher PF if the hot spot is associated with the magnetic pole, where the magnetic field is normal to the emitting surface.We note that the (badly fitting) atmosphere models mentioned in Section 3.2 are not applicable to phase-resolved spectra.
Hot spots could be formed by two different mechanisms.Firstly, they could be due to anisotropic internal heating, i.e., anisotropic heat transfer from the very hot NS interiors to colder surface layers in the presence of a nonuniform magnetic field.Even a moderately strong dipole magnetic field, such as B1055's B ∼ 10 12 G, can lead to noticeable temperature nonuniformity, with a colder equatorial belt and a hotter area of the rest of the surface, slowly increasing towards the magnetic poles of a centered dipole (e.g., Yakovlev 2021).However, although the spectrum of such a NS can be fitted with a BB+BB model, the ratios of hot and cold temperatures and radii are very different from the observed ones, for B1055 and other middle-aged pulsars (Figure 14).
If, in addition to the dipolar magnetic field, there is a toroidal field component in the NS crust, then the two hot spots around the poles of the dipolar component, can have different temperature and sizes (Geppert et al. 2006).The toroidal component leads to the blanketing effect, channeling heat flow along the polar axis and creating an extended cold equatorial belt with an asymmetric temperature distribution that covers most of the stellar surface.
Any scenario that involves hot spot(s) should explain the observed high PFs.Geppert et al. (2006) estimated PFs for various surface temperature distributions involving two warm spots, assuming semi-isotropic BB emission.They obtained a maximum pulsed fraction of PF amp = 33% for an orthogonal rotator, which remains much below the observed values for B1055 (65%-70% in the 0.8−1.2keV energy range and even 38% for energies ≲ 0.5 keV).Perna et al. (2013) considered coupled thermal and magnetic evolution of a NS and computed surface temperature and magnetic field distributions at different ages for a variety of initial magnetic field configurations.Assuming local BB spectra, they considered not only the (semi-)isotropic angular distribution of emitted radiation but also radiation moderately beamed along the local magnetic field direction, mimicking emission from a light-element atmosphere.Interestingly, addition of a large-scale toroidal magnetic component disrupts the north-south symmetry of the temperature distribution due to the Hall effect, resulting in single-peaked light curves and notably heightened PFs.For specific models, these PFs exceed 50%, within the 0.5 − 2 keV energy band.Overall, the presence of a dominant large-scale toroidal magnetic component and/or atmospheric beaming effects could explain high PFs.Conversely, multipolar configurations yield intricate temperature profiles with relatively smaller PFs (Perna et al. 2013).
For B1055, a 50% PF is still too low.Only in combination with atmospheric effects could internal heating (in the presence of a toroidal magnetic field component) explain the large observed PFs.
Hot spots could be also due to external heating by relativistic particles accelerated in the pulsar's magnetosphere, which subsequently precipitate onto the NS polar caps (PCs) and heat them (Ruderman & Sutherland 1975).This mechanism is responsible for thermal X-ray emission from hot spots in old pulsars, including millisecond pulsars, whose bulk surface is too cold to emit X-rays, but it can also contribute to the hot thermal components in middle-aged pulsars (e.g., Pavlov et al. 2002).The range of bolometric luminosities, L BB,h , associated with B1055's hot BB confidence contours is ∼ (3-5)×10 30 erg s −1 (Figure 3).Harding & Muslimov (2001), who considered the PC heating in the space charge limited flow acceleration model, predicted a polar cap luminosity L PC ∼ 2 × 10 −4 Ė = 6 × 10 30 erg s −1 for a pulsar with a characteristic age τ ∼ 500 kyr and a spin period P ∼ 0.2 s, consistent with our estimates of the hot BB luminosity.This suggests that the hot spots, associated with the hot thermal component in the B1055's X-ray emission, could be polar caps heated by relativistic particles precipitating from the pulsar's magnetosphere.
Furthermore, the conventional polar cap radius of B1055, calculated based on a simple "centered" dipole magnetic field geometry, R pc = (2πR 3 NS /cP ) 1/2 ≈ 900 m (assuming a standard NS radius of 13 km), is consistent with the maximum radius from the phase-resolved analysis (Figure 12).Yet, in this scenario, we still lack a clear explanation for the differences in heating between the two magnetic poles.One plausible explanation for the anisotropy in external heating could be associated with off-centered dipole magnetic field configurations (e.g., Harrison & Tademaru 1975).A displacement of the dipole center from the NS center may lead to differences for internal and external heating, thus different pole temperatures.
The effect of an off-centered dipole magnetic field in the case of internal heating was recently studied by Igoshev et al. (2023).They investigated how these configurations affect the temperature patterns on the surface, the light curves, and the spectra of middle-aged pulsars.
For an off-centered dipole field with B = 10 13 G, Igoshev et al. (2023) reported kT BB,h /kT BB,c ≈ 1.1 and R BB,h /R BB,c ≈ 2. These values are notably different from the typical phase-integrated values of kT BB,h /kT BB,c ≈ 2 and R BB,h /R BB,c ≈ 0.1 for middle aged pulsars (Figure 14).In summary, both the internal and external heating mechanisms can contribute to the formation of the observed hot spots on B1055's surface.Detailed atmosphere modeling, with account for beamed emissions, is needed for full understanding of the emission properties of middle-aged pulsars like B1055.However, such modeling is beyond the scope of this paper and is planned for future research.

Nonthermal emission and multiwavelength pulse profile of B1055
The phase shifts between the peaks in multiwavelength pulse profiles (see Figures 6 and 7) can be attributed to different emission mechanisms operating at distinct emission regions that are located at different heights and azimuthal angles.For example, thermal Xray emission comes directly from the NS surface while radio emission is thought to come from different heights above the magnetic poles.The relative shifts and shapes of the multiwavelength pulse profiles are governed by the location of the magnetic poles with respect to the rotation axis and the line of sight, in short the pulsar geometry.

Nonthermal X-ray emission
In our spectral fits with the BB c +BB h +PL model, the nonthermal (PL) component dominates at E ≳ 2 keV (Figures 2 and 10).Its phase-integrated photon index, Γ X ≈ 1.8, substantially exceeds Γ γ ≈ 0.9 in the GeV γ-ray range (Posselt et al. 2023a).The X-ray photon index might show some variations with phase (see Figure 13), but they are not statistically significant even in our relatively deep observation.
The 2 − 5 keV light curve shows one clear peak with a maximum at ϕ ≈ 0.1, i.e., within the broad γ-ray pulse (Figure 7).The maximum of the PL flux phase dependence, F PL (ϕ), is seen at about the same phase (Figure 13).The X-ray maximum slightly lags the leading edge of the γ-ray pulse, but it is ahead of the mean phase of the γ-ray pulse as well as of the radio MP.This Xray peak is also seen in the energy-phase map (Figure 9, right panel).
The 2 − 5 keV light curve also shows a hint of second peak at ϕ ≈ 0.8, near the radio IP.This peak, however is not seen in the F PL (ϕ) curve, which suggests that it is due to an admixture of thermal hot component rather that to the nonthermal emission.
Thus, we cannot exclude the possibility that the nonthermal X-ray emission is connected with the γ-ray emission and perhaps with the radio MP.Unfortunately, the small number of detected nonthermal X-ray photons does not allow us to investigate this connection in more detail.

Constraints on pulsar geometry from the radio and γ-ray light curves
Similar to other radio pulsars with an IP, separated from the MP by about half a period, B1055 is usually interpreted as a nearly-orthogonal rotator, i.e., the magnetic inclination α (the angle between the rotation and magnetic axes) and the viewing angle ζ (between the rotation axis and the line of sight) are not strongly different from 90 • .Although Manchester & Lyne (1977) noted that such pulsars could also be interpreted as a nearly aligned hollow-cone rotators, B1055 is typically interpreted as a nearly-orthogonal rotator (see, e.g., Weltevrede & Wright 2009, and references therein), even though the lag of the IP behind the MP is ∆ϕ = 0.44 (peak-to-peak) rather than 0.5.Based on analysis of swings of linear polarization angle in both the MP and IP and assuming a dipole magnetic field, Weltevrede & Wright (2009) estimate α ≈ 75 • , ζ ≈ 111 • for the MP.They also suggest that the IP arises from emission formed on open field lines close to the magnetic axis at a height ∼ 700 km above the magnetic pole, while the MP originates from field lines lying well outside the polar cap boundary beyond the null surface, and farther away from the magnetic axis, at about the same height.
Additional constraints on the origin of nonthermal emission can be obtained from the analysis of γ-ray light light curves.Pierbattista et al. (2015) computed emission patterns for Fermi-LAT pulsars, including B1055, for four γ-ray emission models: Polar Cap, Slot Gap, Outer Gap and One Pole Caustic, for the core plus cone radio emission model.For each of these models, they generated γ-ray and radio light curves across a parameter space defined by α and ζ.However, none of these models provided a satisfactory explanation for both the γ-ray and radio light curves of B1055.Exceptionally for γ-ray pulsars, B1055's γ-ray pulse can be explained by the Polar Cap model but this model does not work for the radio light curve.The Outer Gap model offered the most reasonable, although far from being perfect, description of the general characteristics of both the γray and radio light curves, with values of α ∼ 77 • and ζ ∼ 87 • .Overall, it seems that those models cannot provide an accurate description of both the γ-ray and radio pulsations.
These models, as well as new models explaining pulsar γ-ray emission, are currently further developed and light curve fitting is used to constrain these models, for instance for the Vela pulsar (e.g., Barnard et al. 2022;Venter et al. 2017, and references therein).However, we are not aware of a recent work focusing on the unusual radio and γ-ray light curves of B1055.

Pulsar geometry with account for X-ray pulsations
Over the past two decades, X-ray observations of B1055 have consistently shown a single peak at energies corresponding to the hot BB component, i.e., one hot spot.However, if this pulsar is indeed an orthogo-nal rotator, as suggested by the radio data, one might expect to observe two peaks in the folded light curve, corresponding to two magnetic poles of a centered magnetic dipole.
In our XMM-Newton data we see not only the peak from the dominating hot spot S1, centered at ϕ ≈ 0.3 − 0.4, but also a hint of a secondary thermal X-ray peak, possibly associated with another hot spot (S2), which lags the main peak by approximately ∆ϕ ≈ 0.4-0.5 (Figure 9, right panel).The S2 hot spot also appears as a hump at ϕ ≈ 0.8 (near the phase of the radio interpulse) in the 0.3 − 0.5 keV light curve (Figure 7).
Thus, the thermal X-ray light curves of B1055 are consistent with two observed magnetic poles, hence a orthogonal rotator geometry, but the poles show an asymmetry in their temperatures.
The radio IP and MP are separated by a distance of approximately ∼ 159 • (∆ϕ = 0.44, Weltevrede & Wright 2009).As illustrated in Figure 5, the peak of the MP exhibits a phase lag of 0.222 with respect to the X-ray peak in the 0.3 − 7 keV, while the IP aligns with the smaller hump (ϕ = 0.7 − 0.8) of the X-ray profile.
Based on the the observed radio pulse profiles and phase dependence of the polarization angles, Weltevrede & Wright (2009) estimated emission heights and plotted a polar cap map showing sites of radio emission production projected onto the NS surface (see their Figure 7).In their cartographic representation, the IP emerges from directly above the polar region while the site of MP generation is offset from the central polar region.This spatial arrangement could potentially explain the time delay between the observed thermal X-ray pulse and the MP, as well as the alignment between the IP and S2.Qualitatively, the lags between the X-ray and radio pulses are also consistent with the two-pole interpretation and further supports the notion of B1055 being an orthogonal rotator.
The 2−5 keV light curve maximum closely aligns with that of the 0.5-1.5 keV light curve (see Figure 7), a phenomenon that may be attributed to a peculiar coincidence or an emission region of the non-thermal Xrays being close to the hot polar region.Conversely, the 2-5 keV maximum falls within the broader γ-ray pulse, which may indicate similar origin sites for the high-energy non-thermal emission.
The phase-energy map above 1.5 keV seems to show a slight trend of the maxima/minima smoothly shifting to smaller phases with increasing energy.For the respective magnetospheric X-ray emitting particles, this trend may indicate different emission heights or a different angular separation from the magnetic pole.Figure 7 also shows a statistically insignificant indication of a second hump in proximity to the location of S2 and the IP.This observation is in agreement with the interpretation of S2 as a second magnetic pole.Puzzlingly, the γ-ray light curve shows a minimum at the respective rotation phase.

Absorption lines in middle-aged pulsars
Absorption features or hints for them were reported in the phase-resolved X-ray spectra of the two other Musketeers, B0656+14 (Arumugasamy et al. 2018;Schwope et al. 2022) and Geminga (Jackson & Halpern 2005), as well as for PSR J1740+1000, another interesting middleaged pulsar (Kargaltsev et al. 2008).However, in later XMM-Newton observations, Rigoselli et al. (2022) found no evidence of spectral lines in the phase-averaged and phase-resolved spectra of this pulsar suggesting a time variation in the spectrum.
An explanation for phase-dependent absorption features could be provided by multipolar magnetic field arcs that trap electrons in close proximity to the NS surface.Through the process of cyclotron resonance scattering, this may lead to the formation of phase-dependent absorption lines in the X-ray spectra (e.g, Arumugasamy et al. 2018).Such absorption lines can be alternatively explained by proton cyclotron lines if the local magnetic fields are high (e.g., Tiengo et al. 2013).Viganò et al. (2014) showed that small-scale temperature variations can also mimic absorption lines.
Our analysis of the phase-resolved spectrum of B1055 did not show any evidence for absorption features, similar to the earlier study by De Luca et al. (2005).However, B1055 is notably fainter than the two other Musketeers.Considering this faintness and the variability seen for J1740+1000, the current non-detection for B1055 does not exclude the possibility of detecting phase-dependent spectral features, deeper X-ray observations achievable with more sensitive telescopes such as the upcoming Athena X-ray Observatory.

CONCLUSION
Using the most comprehensive X-ray data set to date, we find that, in agreement with the findings by De Luca et al. (2005), the spectrum of B1055 is best described by a tripartite model consisting of two blackbody components, presumably emitted from the bulk of the NS surface and small hot spots, and a power-law component, emitted from the pulsar's magnetosphere.
A phase-resolved spectral analysis shows periodic variations in the temperature of both the hot and cold blackbody components.These variations, in tandem with periodic changes in the projected emission areas, suggest non-uniform temperature distributions both over the bulk NS surface and within the alleged hot spot(s).Previously published results on X-ray pulsations of the thermal components showed only one peak per period, likely associated with one visible hot spot, which did not agree with the nearly-orthogonal rotator geometry of B1055, implied by the detailed radio studies of B1055's main pulse and interpulse.In the new X-ray data, we find indications for a second hot spot.However, it appears to be cooler than the already known hot spot.
We explore two potential mechanisms to explain the thermal X-ray emission patterns of B1055: external heating by relativistic particles accelerated in the pulsar's magnetosphere, and internal heating resulting from anisotropic heat transfer due to an offset of the dipole magnetic field and/or the presence of an additional toroidal magnetic field component within the NS's crust.Both mechanisms can, in principle, explain the observed high pulse fractions and phase dependence of the spectral parameters, particularly if beaming effects due to an atmosphere above the hot spots are additionally considered.The faintness and lower temperature of the putative second hot spot could perhaps be explained by an offset dipole, but more detailed modeling is needed.
The complexity of B1055's spectral and temporal characteristics underscores the need for further investigations, including detailed atmosphere modeling and considerations of beamed emissions.
We express our gratitude to Colin J. Clark for his help with the Fermi-LAT data, Michael Freyberg for his assistance with EPIC-pn calibration issues, Prakash Arumugasamy for his support in resolving software-related challenges, as well as to Andrei Igoshev for useful discussions.This work was supported by the Bundesministerium für Wirtschaft und Energie through Deutsches Zentrum für Luft-und Raumfahrt (DLR) under the grant number 50 OR 1917.Support for this work was also provided by the National Aeronautics and Space Administration through the XMM-Newton award No. 80NSSC20K0806.BP acknowledges funding from the STFC consolidated grant to Oxford Astrophysics, code ST/000488.AV thanks the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant #1829740, the Brinson Foundation, and the Moore Foundation; his participation in the program has benefited this work.

Figure 1 .
Figure 1.Count map of the field in the direction of the B1055 for the first observation.The source extraction regions of 30 ′′ (EPIC-pn, top) and 45 ′′ (EPIC-MOS2, bottom) radii are displayed in blue whereas the background region of 60 ′′ radii is shown with magenta color.The images are smoothed with a Gaussian kernel with σ = 1.′′ 25..

Figure 2 .
Figure 2. 0.3-8 keV phase-integrated spectral fit to the combined (2000 and 2019) data set of 12 spectra with 2BB+PL model (Table2).Hot BB, cold BB, and PL components are displayed with red, green, and yellow colors.The fit parameters with a margin of 1σ error are also shown at top of the figure.

Figure 5 .
Figure5.The 0.3 − 7 keV phase-folded light curve (X-ray pulse profile) plotted as a histogram with 20 phase bins and as the sum of 3 harmonics, (N/20)F(ϕ), where F(ϕ) is the normalized pulse profile and N = 86097 is the total number of events in the chosen energy range.The ±3σ uncertainty of the X-ray pulse profile is shown with transparent blue color.The orange, green, and red sine waves correspond to the first, second, and third harmonics, respectively.Area, amplitude and rms pulse fractions in percents, with their 1σ uncertainties, are displayed at the top right corner.
Figure 8.Energy dependence of three backgroundcorrected pulse fractions plotted with 90% confidence level errors for the energy ranges defined in Table4.The orange, blue and green dots correspond to amplitude, area and rootmean-square pulse fractions.

Figure 9 .
Figure 9. Different representations of the EPIC-pn X-ray events of B1055 in phase-energy space.The right panel shows a phase-energy map where only the normalization is applied, while the left panel shows deviations from the phase-averaged values in each phase-energy bin -see Equations (3) and (4).

Figure 10 .
Figure 10.Fits to phase-resolved spectra for 5 equal-sized phase bins of the 2019 EPIC-pn data.Hot BB, cold BB, and PL components are displayed with red, green, and yellow colors.The fit parameters with a margin of 1σ error are shown at the top of each panel.

Figure 11 .
Figure11.Variations of hot BB (left) and cold BB (right) with phase in the temperature-normalization plane.Red and green contours correspond to confidences levels of 68.3% and and 99% respectively.Additionally, we depict the confidence contours of phase-integrated spectra using brown (68.3%) and black (99%) colors.Errors are estimated for two parameters of interest.Small dots in some contour plots are small contours themselves representing local minima, indicating a complicated shape of the χ 2 surface.

Figure 12 .
Figure 12.Phase dependencies of the temperatures (top), radii of equivalent emitting sphere (middle), and unabsorbed fluxes in 0.3 − 8.0 keV band (bottom) for hot BB (left) and cold BB (right) components.The vertical error bars show 1σ uncertainties of the fitting parameters.Blue lines represent sine function fits to the data for visualization purposes.

Figure 13 .
Figure 13.Phase dependence of the photon index Γ (top), normalization KPL (middle), and unabsorbed fluxes in 0.3 − 8.0 keV band (bottom) for the PL component.The vertical error bars show 1σ uncertainties of the fitting parameters.Blue lines represent sine function fits to the data for visualization purposes.

Table 1 .
XMM-Newton EPIC observations of B1055 used in this studyNote-Ti-ME: timing mode with medium filter, SW-TN: small window with thin filter, FF-ME: full frame with medium filter.Total count rate (cts ks −1 ) in 0.3 − 10 keV for 2019 observations and 0.4 − 10 keV for 2000 observations.

Table 3 .
The radio timing properties of B1055.