#  XMM-Newton observations of PSR J0554+3107: pulsing thermal emission from a
cooling high-mass neutron star

A. S. Tanashkin, A. V. Karpova, A. Y. Potekhin, Y. A. Shibanov and D. A.
Zyuzin
Ioffe Institute, Politekhnicheskaya 26, St. Petersburg, 194021, Russia E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

XMM-Newton observations of the middle-aged radio-quiet $\gamma$-ray pulsar
J0554+3107 allowed us, for the first time, firmly identify it in X-rays by
detection of pulsations with the pulsar period. In the 0.2–2 keV band, the
pulse profile shows two peaks separated by about a half of the rotation phase
with the pulsed fraction of $25\pm 6$ per cent. The profile and spectrum in
this band can be mainly described by thermal emission from the neutron star
with the hydrogen atmosphere, dipole magnetic field of $\sim 10^{13}$ G and
non-uniform surface temperature. Non-thermal emission from the pulsar
magnetosphere is marginally detected at higher photon energies. The spectral
fit with the atmosphere$+$power law model implies that J0554+3107 is a rather
heavy and cool neutron star with the mass of 1.6–2.1 M⊙, the radius of
$\approx 13$ km and the redshifted effective temperature of $\approx 50$ eV.
The spectrum shows an absorption line of unknown nature at $\approx 350$ eV.
Given the extinction–distance relation, the pulsar is located at $\approx 2$
kpc and has the redshifted bolometric thermal luminosity of $\approx 2\times
10^{32}$ erg s-1. We discuss cooling scenarios for J0554+3107 considering
plausible equations of state of super-dense matter inside the star, different
compositions of the heat-blanketing envelope and various ages.

###### keywords:

stars: neutron – pulsars: general – pulsars: individual: PSR J0554+3107

††pubyear: 2022††pagerange: XMM-Newton observations of PSR J0554+3107: pulsing
thermal emission from a cooling high-mass neutron star–B

##  1 Introduction

Studies of thermal emission from neutron stars (NSs) is one of the ways to
investigate the properties of super-dense nuclear matter in their interiors
(e.g. Yakovlev & Pethick, 2004). If thermal emission originates from the
entire NS surface, the comparison of the measured thermal luminosity (or,
equivalently, mean effective temperature) with predictions of NS cooling
theories can set constraints on the equation of state (EoS) of such matter.
Middle-aged (104–106 year-old) stars are the best targets for such analysis
since thermal components usually dominate over non-thermal ones in their X-ray
spectra. However, effective temperatures and ages are estimated only for a few
dozens of NSs, which is not enough to make definite conclusions (Potekhin et
al., 2020).

Table 1: J0554 parameters obtained from Pletsch et al. (2013). Numbers in parentheses denote 1$\sigma$ uncertainties relating to the last significant digit quoted. R.A. (J2000) | 05h54m05$\aas@@fstack{s}$01(3)
---|---
Dec. (J2000) | +31∘07′41′′(4)
Galactic longitude $l$, deg | 179.058
Galactic latitude $b$, deg | 2.697
Spin period $P$, ms | 465
Spin frequency $f$, Hz | 2.15071817570(7)
Frequency derivative $\dot{f}$, Hz s-1 |  $-$0.659622(5) $\times$ 10-12
Frequency second derivative $\ddot{f}$, Hz s-2 | 0.18(2) $\times$ 10-23
Epoch of frequency determination, MJD | 55214
Data time span, MJD | 54702–56383
Solar system ephemeris model | DE405
Characteristic age $t_{\mathrm{c}}$, kyr | 51.7
Spin-down luminosity $\dot{E}$, erg s-1 | 5.6 $\times$ 1034
Characteristic magnetic field $B_{\mathrm{c}}$, G | 8.2 $\times$ 1012

The middle-aged radio-quiet PSR J0554+3107 (hereafter J0554) was discovered by
Einstein@Home111http://einstein.phys.uwm.edu in a blind search of Fermi Large
Area Telescope (LAT) $\gamma$-ray data (Pletsch et al., 2013). Its parameters
are presented in Table 1. Upper limits on the integral flux density at 111,
150 and 1400 MHz are 0.5, 1.2 and 0.066 mJy, respectively (Tyul’bashev et al.,
2021; Grießmeier et al., 2021; Pletsch et al., 2013).

The pulsar is projected onto the supernova remnant (SNR) G179.0+2.6 and likely
associated with it (Pletsch et al., 2013). This is a shell-type oxygen-rich
remnant (How et al., 2018) with the diameter of about 70 arcmin. Its age is
unclear: the large diameter and the low surface brightness imply the age of
$\sim$ 10–100 kyr, but the radial configuration of the magnetic field is
typical for young SNRs (Fuerst & Reich, 1986; Gao et al., 2011). There are
different estimates of the distance to the remnant, from about 3 to 6 kpc
(Case & Bhattacharya, 1998; Guseinov et al., 2003; Pavlovic et al., 2014),
obtained through empirical correlations between the radio surface brightness
and the diameter of the SNR. However, recent studies show that G179.0+2.6 can
be much closer, at $\approx 0.9$ kpc (Zhao et al., 2020). This result is based
on an apparent interstellar extinction jump at this distance along the remnant
line of sight and on the assumption that it is associated with the dust formed
by the SNR. However, the detected extinction jump may correspond to a
foreground molecular cloud and not the remnant itself. If so, this value
should be considered as the lower limit to the distance.

Since J0554 has not been detected in radio, it is not possible to obtain the
distance basing on the dispersion measure. The only available estimate is the
so-called ‘pseudo-distance’ of $\approx 1.9$ kpc calculated from the empirical
relation between the distance and the pulsar $\gamma$-ray flux (Saz Parkinson
et al., 2010). Though such an estimate is rather uncertain (within a factor of
2–3), it is consistent with the association of J0554 with G179.0+2.6.

Zyuzin, Karpova & Shibanov (2018) found the likely X-ray counterpart of J0554
in the 1st Swift-XRT Point Source Catalogue (1SXPS Evans et al., 2014).
Seventeen counts were detected from the source in the 11 ks exposure, and
thirteen of them are in the 0.3–1 keV band, which indicates that the source is
soft. Indeed, fitting the X-ray spectrum with a power law model, Zyuzin et al.
(2018) obtained the photon index $\Gamma>4$, which is too high for pulsars and
may indicate the presence of a thermal-like spectral component. Pure thermal
emission models – the black body and the NS atmosphere – resulted in NS
surface temperatures of $\approx$ 50–100 eV. In any case, the obtained model
parameters remained very uncertain because of the low Swift count statistics.
We therefore performed deeper observations with XMM-Newton to confirm the
counterpart of J0554 and clarify its X-ray properties.

Here we present results of these observations. The data and the reduction
procedure are presented in Sec. 2. Imaging is described in Sec. 3. Timing and
spectral analyses are presented correspondingly in Sec. 4 and Sec. 5. We
discuss the results in Sec. 6 and make conclusions in Sec. 7.

##  2 Observations and data reduction

Figure 1: High-energy light curves extracted from the FoVs of the MOS ($>10$
keV; top) and PN (10–12 keV; bottom) cameras. The dotted lines show the
thresholds applied to filter the MOS1 and PN data for periods of background
flaring activity (no filtering was applied to the MOS2 data).

45-ks XMM-Newton observations of the J0554 field were carried out on 2021
October 7 (ObsID 0883760101, PI A. Karpova). The European Photon Imaging
Camera Metal Oxide Semiconductor (EPIC-MOS) detectors were operated in the
full frame mode and the EPIC-pn (PN hereafter) detector – in the large window
mode. The corresponding imaging areas are 28 arcmin $\times$ 28 arcmin and
13.5 arcmin $\times$ 26 arcmin. For all instruments, the thin filter was
chosen.

We used the XMM-Newton Science Analysis Software (xmm-sas) v.19.1.0 for
analysis. The emproc and epproc routines were utilised to reprocess the data.
To filter out the periods of background flares, we extracted high-energy light
curves from the fields of view (FoVs) of all EPIC detectors. Only one short
flare is present (see Fig. 1), mostly affecting the PN detector, but also
having a slight impact on the MOS1 light curve, so we cleaned only the PN and
MOS1 data, leaving the MOS2 data unfiltered. As a result, the effective
exposures are 44.3, 44.5 and 38.8 ks for MOS1, MOS2 and PN cameras,
respectively. Single, double, triple and quadruple pixel events were selected
for MOS (pattern $\leq$ 12) and single and double pixel events – for PN
(pattern $\leq$ 4).

##  3 Imaging

Figure 2: XMM-Newton combined (MOS1+MOS2+PN) image of the J0554 field obtained
in soft (0.2–1 keV, red), medium (1–2 keV, green) and hard (2–12 keV, blue)
energy bands. The bright soft point source in the centre of the image is the
X-ray counterpart of the pulsar, whose $\gamma$-ray position is marked by the
cross. The position uncertainty is significantly smaller than the marker size
and is not shown. The dashed 38.5-arcsec radius circle shows the region used
for the background extraction in the timing and spectral analyses.

Figure 3: Comparison of the XMM-Newton PSF and the observed radial brightness
profile obtained from the PN data in the 0.2–10 keV energy band.

Using the ‘images’ script222https://www.cosmos.esa.int/web/xmm-newton/images
(Willatt & Ehle, 2016), we created combined MOS+PN exposure-corrected images
in different energy bands. The resulting image is shown in Fig. 2. The J0554
candidate counterpart is seen as a bright source in the centre. Its position
R.A. = 05h54m05$.\\!\\!^{\mathrm{s}}$067(10) and Dec. =
+31∘07′41$\aas@@fstack{\prime\prime}$40(13) was derived by the edetect_chain
task using the data from all EPIC detectors (numbers in parentheses are
1$\sigma$ pure statistical uncertainties). Taking into account the XMM-Newton
absolute pointing accuracy of 1.2
arcsec333https://xmmweb.esac.esa.int/docs/documents/CAL-TN-0018.pdf, the
coordinates are in agreement within 1$\sigma$ with the J0554 ones obtained
from the Fermi data (Table 1).

From Fig. 2 one can see that the J0554 candidate counterpart is a soft source.
No extended emission like a pulsar wind nebula (PWN) or misaligned outflows
are seen in the pulsar vicinity. However, the XMM-Newton point spread function
(PSF) is rather wide and the compact PWN if exists can be blurred with the
pulsar. To treat this possibility more carefully we used the eradial tool,
which extracts the source radial brightness profile and fits the PSF to it.
The result is shown in Fig. 3, and we can conclude that the radial profile is
consistent with the PSF and no compact PWN is resolved. We also do not resolve
any extended emission in the XMM-Newton FoV that could be identified with the
SNR G179.0+2.6, whose radio shell is outside the FoV.

##  4 Timing analysis

The PN detector operating in the large window mode provides the time
resolution of $\approx 48$ ms, which is sufficient to search for regular X-ray
pulsations with the 465 ms spin period of J0554. We used the barycen task and
DE405 ephemeris to apply the barycentric correction and then extracted
filtered for background flaring activity events in the 0.2–2 keV band from the
22-arcsec radius circle aperture centred at the pulsar position. The aperture
was selected using the eregionanalyse task which produces the optimum
extraction radius basing on signal-to-noise ratio. This results in 1026 source
counts or $\approx 98$ per cent of the total number of source counts detected
in the whole PN energy band. To search for the pulsations, we performed the
$Z^{2}_{n}$-test (Buccheri et al., 1983) for the number of harmonics $n$ from
1 to 5 using the 0.8 mHz window around the spin frequency of 2.150474376(14)
Hz expected at the epoch of the XMM-Newton observation (MJD 59494) according
to the Fermi timing results (Table 1). We detected pulsations at the frequency
of 2.150493(2) Hz and found statistically significant contributions from two
leading harmonics. The resulting periodogram is shown in Fig. 4. The maximum
$Z^{2}_{2}=42.7$, which implies the detection confidence level of $\approx
4.7\sigma$444See Appendix A for details..

Figure 4: $Z^{2}_{2}$-test periodogram for J0554 in the 0.2–2 keV energy band.

Figure 5: Pulse profiles of J0554 in the 0.2–2 keV and 2–4 keV energy bands.
Two phase cycles are plotted for clarity. The dashed red line in the bottom
panel indicates the background level. In the 0.2–2 keV band the background
level ($\approx 19$ counts per bin) is significantly below the minimum of the
pulse profile and the bottom boundary of the plot.

We note that the frequency detected in X-rays is somewhat higher than the
predicted one. This can be due to the pulsar timing noise and/or glitches,
which might have occurred during the $\sim 8$ yr period between the last
$\gamma$-ray observations of 2013, included by Pletsch et al. (2013) into the
Fermi timing solution, and the XMM-Newton observations of 2021. For instance,
the relative difference $\log(\Delta f/f)\approx-5.1$ is consistent with the
distribution of the relative glitch sizes $\log(\Delta f_{\mathrm{g}}/f)$
observed for pulsars (Lower et al., 2021).

Table 2: Best-fitting parameters for different spectral models†.

Model | 2BB + PL | 2BB + PL‡ | (2BB + PL) $\times$ gabs | $(\mathrm{NS130+PL})\times\textsc{gabs}$ | $(\mathrm{NS133+PL})\times\textsc{gabs}$
---|---|---|---|---|---
$N_{\mathrm{H}}$, $10^{21}$ cm-2 | $3.1_{-0.2}^{+0.2}$ | $5.2_{-0.8}^{+0.4}$ | $2.3_{-0.4}^{+0.5}$ | $1.62_{-0.06}^{+0.08}$ | $1.66_{-0.07}^{+0.08}$
$D$, kpc | $5.5_{-0.5}^{+0.3}$ |  | $4.1_{-1.4}^{+1.4}$ | $2.0_{-0.4}^{+0.2}$ | $2.0_{-0.2}^{+0.4}$
$\Gamma$ | $2.0^{\mathrm{fixed}}$ | $2.0^{\mathrm{fixed}}$ | $2.0^{\mathrm{fixed}}$ | $2.2_{-0.4}^{+0.6}$ | $2.3_{-0.4}^{+0.5}$
$K$, $10^{-6}$ ph cm-2 s-1 keV-1 | $2.2_{-0.6}^{+0.6}$ | $1.9_{-0.7}^{+0.6}$ | $2.1_{-0.5}^{+0.6}$ | $1.6_{-0.7}^{+1.3}$ | $2.0_{-0.7}^{+1.5}$
log $L_{\mathrm{X}}$, erg s-1 | $31.27_{-0.15}^{+0.13}$ | $29.62_{-0.12}^{+0.20}+\mathrm{log}D_{0.9}^{2}$ | $31.08_{-0.28}^{+0.20}$ | $30.21_{-0.30}^{+0.15}$ | $30.27_{-0.19}^{+0.19}$
$\alpha$, deg |  |  |  | $70_{-20}^{+20}$ | $60_{-20}^{+20}$
$\zeta$, deg |  |  |  | $60_{-10}^{+20}$ | $80_{-30}^{+10}$
$M$, M⊙ |  |  |  | $1.9_{-0.2}^{+0.2}$ | $1.8_{-0.2}^{+0.3}$
$R$, km |  |  |  | $13.5_{-1.7}^{+1.2}$ | $13.0_{-1.3}^{+1.5}$
$R^{\infty}$, km |  |  |  | $16.9_{-1.3}^{+1.5}$ | $17.1_{-1.3}^{+1.3}$
$k_{\mathrm{B}}T^{\infty}$, eV |  |  |  | $47_{-2}^{+2}$ | $49_{-2}^{+2}$
$R^{\infty}_{\mathrm{cold}}$, km | $19_{-5}^{+1}$ | $14_{-6}^{+5}\ D_{0.9}$ | $19_{-9}^{+1}$ |  |
$R^{\infty}_{\mathrm{hot}}$, km | $1.12_{-0.11}^{+1.05}$ | $0.33_{-0.20}^{+0.23}\ D_{0.9}$ | $1.19_{-0.17}^{+0.87}$ |  |
$k_{\mathrm{B}}T^{\infty}_{\mathrm{cold}}$, eV | $86_{-4}^{+5}$ | $65_{-3}^{+7}$ | $84_{-10}^{+6}$ |  |
$k_{\mathrm{B}}T^{\infty}_{\mathrm{hot}}$, eV | $156_{-25}^{+12}$ | $140_{-17}^{+23}$ | $135_{-13}^{+19}$ |  |
log $L^{\infty}$, erg s-1 | $33.32_{-0.19}^{+0.10}$ | $32.79_{-0.39}^{+0.11}+\mathrm{log}D_{0.9}^{2}$ | $33.23_{-0.61}^{+0.12}$ | $32.25_{-0.11}^{+0.12}$ | $32.33_{-0.10}^{+0.10}$
$E_{0}$, eV |  |  | $370_{-70}^{+30}$ | $340_{-40}^{+40}$ | $350_{-50}^{+30}$
$\sigma$, eV |  |  | $24_{-14}^{+18}$ | $25_{-9}^{+13}$ | $24_{-8}^{+14}$
$\tau$, eV |  |  | $>40$ | $>680$ | $>890$
EW, eV |  |  | $<430$ | $150_{-40}^{+120}$ | $180_{-70}^{+70}$
$W$/d.o.f. | $240/216$ | $231/216$ | $228/213$ | $228/211$ | $228/211$
$\chi^{2}_{\nu}$/d.o.f. | $1.46/40$ | $1.16/40$ | $1.14/37$ | $1.27/34$ | $1.27/34$

  * •

† $N_{\mathrm{H}}$ is the absorbing column density, $D$ is the distance,
$\Gamma$ is the photon index, $K$ is the PL normalisation, $L_{\mathrm{X}}$ is
the non-thermal luminosity in the 2–10 keV band, $\alpha$ is the angle between
the axis of rotation and the magnetic axis, $\zeta$ is the angle between the
rotation axis and the line of sight, $M$ is the NS mass, $R$ and
$R^{\infty}=R(1+z_{\mathrm{g}})$ are the intrinsic and apparent radii of the
NS, $T^{\infty}=T/(1+z_{\mathrm{g}})$ is the atmosphere redshifted effective
temperature, $R^{\infty}_{\mathrm{cold}}$ and $R^{\infty}_{\mathrm{hot}}$ are
the apparent radii of the equivalent emitting spheres,
$T^{\infty}_{\mathrm{cold}}$ and $T^{\infty}_{\mathrm{hot}}$ are the
redshifted effective temperatures of the cold and hot BB components,
$L^{\infty}=L/(1+z_{\mathrm{g}})^{2}$ is the apparent bolometric thermal
luminosity, $E_{0}$, $\sigma$, $\tau$ and EW are the absorption line centre,
width, depth and equivalent width, $z_{\mathrm{g}}$ is the gravitational
redshift. All errors are at 1$\sigma$ credible intervals, while the lower and
upper limits are set at 98 per cent. The last two rows provide the minimum
values of $W$-statistics, which was used as a log-likelihood in the MCMC
procedure, and $\chi^{2}_{\nu}$ calculated for the spectra which were grouped
to ensure at least 25 counts per bin.

  * •

‡ In this case, in contrast to all other models, the $N_{\mathrm{H}}$–$D$
relation was not used in the fitting procedure (see text for details).
$D_{0.9}\equiv D/0.9$ kpc, where 0.9 kpc is the lower limit of the distance to
J0554 obtained from its association with the SNR G179.0+2.6 (Zhao et al.,
2020).

Folded X-ray pulse profile in the 0.2–2 keV energy band is presented in the
upper panel of Fig. 5. Following the approach suggested by Becker & Trümper
(1999) we calculated the optimal number of phase bins to be 11. However, this
value exceeds the maximum of 9 bins set by the time resolution of the PN
camera and the pulsar period, so we used the latter value to construct the
folded light curve. Two peaks separated by about a half of the pulsar rotation
phase are clearly resolved. The pulsed fraction (PF) of the emission was
calculated from the photon phases using the bootstrap method outlined in
Swanepoel, de Beer & Loots (1996). The result obtained with this technique is
not affected by binning effects, which is particularly important in our case
since the number of counts is low and the phase bins are wide. The resulted
background-corrected PF in the 0.2–2 keV band is $25\pm 6$ per cent. In the
lower panel of Fig. 5 we also show the pulse profile in the 2–4 keV energy
band. The pulsar possibly demonstrates single-peaked pulsations, but the count
statistics in this band is too low to make definite conclusions.

##  5 Spectral analysis

We extracted the pulsar spectra from the 19-arcsec radius aperture using
evselect task. This radius was chosen using the eregionanalyse task as the
optimum value in terms of signal-to-noise ratio for the 0.2–10 keV band. The
redistribution matrix and the ancillary response files were generated by
rmfgen and arfgen tools. For the background, we chose the 38.5-arcsec radius
circular source-free region, located at the same CCD chip and approximately
the same CCD RAWY pixel position as the source (see Fig. 2), in accordance
with general recommendations by the EPIC Consortium555See, e.g., footnote 3,
pages 28–29. As a result, we obtained 221, 246 and 1001 net counts in the
0.2–10 keV band from the MOS1, MOS2 and PN data, respectively. We fitted the
spectra simultaneously with the X-Ray Spectral Fitting Package (xspec)
v.12.11.1 (Arnaud, 1996) using the tbabs model with the wilm abundances
(Wilms, Allen & McCray, 2000) to take into account the interstellar medium
(ISM) absorption.

At first, we grouped the spectra to ensure at least 25 counts per energy bin
and used $\chi^{2}$-statistics to preliminarily check, which of the models,
typically used to describe pulsars X-ray emission, can fit the data. The
single powerlaw (PL) model, corresponding to the NS magnetosphere emission,
though statistically acceptable ($\chi^{2}_{\nu}$/d.o.f. = 1.03/41, d.o.f.
$\equiv$ degrees of freedom) resulted in a photon index $\Gamma\approx 7$,
which is too high to be consistent with typical slopes of non-thermal X-ray
spectra of pulsars (e.g. Kargaltsev & Pavlov, 2008). Addition of the blackbody
(BB) component (model bbodyrad in xspec) still gives a too high photon index,
while the single BB model fits the spectrum worse ($\chi^{2}_{\nu}$/d.o.f. =
2.02/43). In contrast, the BB + BB model, composed of two BB components with
different temperatures and radii of equivalent emitting spheres, is
statistically acceptable with $\chi^{2}_{\nu}$/d.o.f. = 1.09/41. This model is
usually assumed to describe thermal emission from some colder and hotter areas
of the NS surface. However, there is some flux excess above the model at
energies $\gtrsim 2$ keV. Although its significance is low, this may indicate
the presence of the non-thermal emission. Thus, we added the PL component to
the BB + BB model, which resulted in $\chi^{2}_{\nu}$/d.o.f. = 1.04/39 and a
reasonable photon index $\Gamma\lesssim 3$.

We also constructed the hydrogen atmosphere models nsmdintb for the NS with a
dipole magnetic field to fit the thermal spectral component. In these models,
which are described in Appendix B, free parameters are: the NS mass $M$ and
radius $R$, the distance $D$, the angle $\alpha$ between the axis of rotation
and the magnetic axis, the angle $\zeta$ between the rotation axis and the
line of sight, and the redshifted effective temperature
$T^{\infty}=T/(1+z_{\mathrm{g}})$, where $z_{\mathrm{g}}$ is the gravitational
redshift at the NS surface. The effective surface temperature $T$ is defined
by the total thermal luminosity $L$ according to the relation $L\equiv
4\pi\sigma_{\mathrm{SB}}R^{2}T^{4}$, where $\sigma_{\mathrm{SB}}$ is the
Stefan-Boltzmann constant. We used two values of the magnetic field at the
pole close to the estimates based on the dipole spin-down formula (see
Appendix B), $B_{\mathrm{p}}=10^{13}$ G and $2\times 10^{13}$ G (hereafter
NS130 and NS133, respectively). We obtained rather poor fits with
$\chi^{2}_{\nu}$/d.o.f. = 1.60/39 for NS130 and $\chi^{2}_{\nu}$/d.o.f. =
1.55/39 for NS133. However, we found that addition of an absorption Gaussian
line (model gabs in xspec) at $\approx 0.35$ keV significantly improves the
fit resulting in $\chi^{2}_{\nu}$/d.o.f. = 1.15/36 and 1.17/36. Addition of
the PL component to describe the flux excess at high energies slightly
improves the fits leading to $\chi^{2}_{\nu}$/d.o.f. = 1.15/34 and 1.07/34.
Basing on these preliminary tests, we further focused on the models 2BB + PL
and $(\textsc{nsmdintb}+\mathrm{PL})\times\textsc{gabs}$ as the ones providing
the best fit statistics.

Figure 6: Interstellar extinction–distance relation in the J0554 direction
obtained from the dust map of Green et al. (2019).

The number of collected source counts is not large and for more rigorous
analysis we regrouped the spectra to ensure at least 1 count per energy bin.
This allows us to obtain the most robust estimates of the J0554 spectral
parameters. We applied $W$-statistics (Wachter, Leach & Kellogg, 1979)
appropriate for Poisson data with Poisson background666See
https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html.
and performed the fitting using the Markov chain Monte Carlo (MCMC) technique.
We utilised the Bayesian parameter estimation procedure using pyxspec
interface and the python package emcee (Foreman-Mackey et al., 2013), which
implements the affine-invariant MCMC sampler developed by Goodman & Weare
(2010). The best-fitting model parameters, defined as the ones corresponding
to the maximum values of probability density, were derived from the sampled
posterior distributions together with their 1$\sigma$ credible intervals.

The Bayesian inference also allows us to include some additional information
in the fitting procedure. We used the 3D map777http://argonaut.skymaps.info/
of the dust distribution in the Galaxy based on Gaia, Pan-STARRS 1 and 2MASS
data (Green et al., 2019) to obtain the extinction–distance relation in the
J0554 direction. In Fig. 6 we show five representative samples of this
relation drawn from the posterior distribution of the distance–reddening
profiles (see Green et al. 2019 for details). The following procedure was
used. At each step of the MCMC fitting, we randomly take one of the five
samples and use the relation $N_{\mathrm{H}}$ = $a\times 10^{21}E(B-V)$ to
convert the selective reddening $E(B-V)$ into the equivalent hydrogen column
density $N_{\mathrm{H}}$, which is responsible for the ISM absorption in
X-rays. The conversion factor $a$ is drawn from the normal distribution with
the mean of 8.9 and the standard deviation of 0.4 according to the empirical
relation $N_{\mathrm{H}}$ = $(2.87\pm 0.12)\times 10^{21}A_{V}$ cm-2 (Foight
et al., 2016) and implying the standard reddening law $A_{V}=3.1E(B-V)$, where
$A_{V}$ is the optical extinction. Then we compute the distance $D$ using
linear interpolation between the closest sample values
$\underline{N_{\mathrm{H}}}$ and $\overline{N_{\mathrm{H}}}$ such that
$\underline{N_{\mathrm{H}}}\leq$ $N_{\mathrm{H}}$
$\leq\overline{N_{\mathrm{H}}}$.

In contrast to nsmdintb models, where both the NS radius and the distance are
free parameters, the BB model has a normalisation
$N=(R_{\mathrm{BB}}^{\infty}/D)^{2}$ as a free parameter, where
$R_{\mathrm{BB}}^{\infty}$ is the apparent radius of the equivalent emitting
sphere. Thus, in the case of the 2BB + PL model at each step we independently
sampled the radii $R^{\infty}_{\mathrm{cold}}$ and $R^{\infty}_{\mathrm{hot}}$
for colder and hotter emitting areas and then calculated normalisations using
these values together with the computed distance. We constrained both radii to
be $\leq 20$ km as a reasonable value of an NS radius measured by a distant
observer.

For the $(\textsc{nsmdintb}+\mathrm{PL})\times\textsc{gabs}$ models, we used a
prior $\alpha+\zeta~{}\geq 90^{\circ}$, which follows from the shape of the
folded light curve and the PF value. The PL photon indices $\Gamma$ were
constrained in the range of 0.5–3 typical for pulsars (e.g. Kargaltsev &
Pavlov, 2008). However, in the case of the 2BB + PL model, $\Gamma$ tends to
the upper bound, while the temperature of the hotter BB component in general
takes slightly lower values than that in the case of the pure BB + BB model.
This behaviour is typical for situations when the PL component competes with
the hotter BB component at lower energies, while the number of counts above
$\sim 2$ keV is very low and the PL slope there is poorly constrained. This is
exactly the case of J0554. Thus, we applied the fixed photon index $\Gamma=2$
for the 2BB + PL model.

The best-fitting parameters are presented in Table 2. The pulsar spectrum with
the best-fitting model $(\mathrm{NS133+PL})\times\textsc{gabs}$ is shown in
Fig. 7. In Fig. 8 we present a corner plot of posterior distribution functions
for some parameters of this model. To check the fit quality, we calculated
values of the $\chi^{2}$-statistics using the spectra grouped to ensure at
least 25 counts per bin. They are provided in the last row of Table 2. These
$\chi^{2}$ values differ from the preliminary ones partially due to inclusion
of parameter constraints in the fitting procedure.

The Gaussian line model gabs describing the spectral feature is represented as

$\mathfrak{G}(E)=\mathrm{exp}\left(-\frac{\tau}{\sqrt{2\pi}\sigma}e^{-\frac{(E-E_{0})^{2}}{2\sigma^{2}}}\right),$
(1)

where $E_{0}$, $\sigma$ and $\tau$ are the line centre, width and depth.
$\tau$ is poorly defined from the spectral fits and we provide only its lower
limits. We also calculated the line equivalent width (EW) defined as:

$\mathrm{EW}=\int(1-\mathfrak{G}(E))\mathrm{d}E,$ (2)

which is better constrained by the fits.

Figure 7: J0554 spectrum, best-fitting model
$(\mathrm{NS133+PL})\times\textsc{gabs}$ and residuals. The dashed lines show
the model components. Data from different detectors are shown by different
colours. For illustrative purposes, spectra were regrouped using xspec command
rebin to reach at least 2$\sigma$ detection significance in each new bin (but
no more than 12 adjacent bins were allowed to be combined).

From the 2nd column of Table 2, one can see that the 2BB + PL model fails to
fit the data when the $N_{\mathrm{H}}$–$D$ relation is used. We have found
that it describes the spectrum poorly at energies $\lesssim 0.4$ keV, where
the model flux becomes systematically higher than the observed one. Exclusion
of the $N_{\mathrm{H}}$–$D$ prior removes the discrepancy in the soft band and
results in the statistically acceptable fit (the 3rd column of Table 2).
However, for reasonable BB normalisations, which imply the apparent radii of
$\lesssim 20$ km and the distances of $\gtrsim 0.9$ kpc, the obtained column
densities are about 1.5 times higher than the maximum value provided by the
dust map of Green et al. (2019, see Fig. 6). Recalling that the addition of
the Gaussian absorption line significantly improves the fit in the case of the
atmosphere models, we tried to add the line to the 2BB + PL model as well.
This leads to the good fit without $N_{\mathrm{H}}$ anomalies even when the
$N_{\mathrm{H}}$–$D$ relation is implemented (see the 4th column in Table 2).

Figure 8: 1D and 2D marginal posterior distribution functions for some
parameters of the $(\mathrm{NS133+PL})\times\textsc{gabs}$ model. Contours in
2D distributions correspond to 39, 68 and 86 per cent confidence levels. In 1D
distributions, vertical dashed lines indicate 1$\sigma$ credible intervals.

##  6 Discussion

###  6.1 J0554 spectral and timing properties

The time-integrated X-ray spectrum of J0554 can be well described by the
composite emission models containing the thermal and non-thermal spectral
components. For the former, we tried the BB model, which imitates the emission
from the solid or liquid state NS surface, and the hydrogen atmosphere models
with dipole magnetic fields.

In the case of the 2BB + PL model without the $N_{\mathrm{H}}$–$D$ prior, we
assume that the cold BB component originates from the bulk of the stellar
surface while the hot component describes emission from a hot spot (e.g. polar
caps). Then the distance to J0554 should be $\lesssim 2$ kpc, otherwise the
radius of the cold component will be inadequately large (see the 3rd column of
Table 2). However, the model requires a much higher column density than
expected for such distance according to the dust map by Green et al. (2019).
This means that either the model, being formally statistically acceptable,
results in a physically inappropriate $N_{\mathrm{H}}$ and hence wrong NS
parameters, or the considered $N_{\mathrm{H}}$–$D$ relation is incorrect in
the J0554 direction. Since the pulsar is projected onto the SNR G179.0+2.6,
the derived absorption excess could be provided by some filament of the
remnant even if it is not visible in X-rays. Unfortunately, there are no X-ray
extragalactic objects in the pulsar field which are bright enough to make an
independent estimate of the maximum absorption in the pulsar direction.
However, we checked some other extinction maps by Drimmel, Cabrera-Lavers &
López-Corredoira (2003), Chen et al. (2014), Chen et al. (2019), Sale et al.
(2014) and Lallement et al. (2014, 2018); Capitanio et al. (2017). All of them
provide $N_{\mathrm{H}}\lesssim 3.5\times 10^{21}$ cm-2 at 2 kpc which is
significantly smaller than the best-fitting value. Moreover, basing on the [O
iii]/H$\alpha$ emission ratio measured for the SNR G179.0+2.6 associated with
J0554, How et al. (2018) argued that $E(B-V)$ should not be much greater than
0.3 (i.e. $N_{\mathrm{H}}$ $\lesssim 2.7\times 10^{21}$ cm-2). This makes the
2BB + PL model hardly acceptable without additional constraints and model
components.

Implementation of the $N_{\mathrm{H}}$–$D$ prior along with the addition of
the low energy absorption feature solves the 2BB + PL problem, resulting in
about twice lower column density and acceptable fit statistics. If we assume
that the cold BB component describes the emission from the bulk of the stellar
surface while the hot component corresponds to a hot spot, then the size of
the latter is larger than the ‘standard’ polar cap of a radius of a few
hundred meters, expected according to the model of Sturrock (1971). The
obtained effective temperatures are in agreement with results for other NSs of
similar age (Potekhin et al., 2020).

The $(\textsc{nsmdintb}+\mathrm{PL})\times\textsc{gabs}$ models are also
plausible. Parameters obtained for two different magnetic fields are very
similar (see the last two columns of Table 2).888We also tried analogous
models with other field strengths and found that a decrease of
$B_{\mathrm{p}}$ to $\sim 2\times 10^{12}$ G worsens the fitting statistics,
while an increase to $5\times 10^{13}$ G does not noticeably change the
results. These models indicate that J0554 should be a rather heavy NS, with
the mass of 1.6–2.1 M⊙. Its redshifted effective temperature $T^{\infty}=48\pm
3$ eV ($0.56\pm 0.3$ MK), which is about twice lower than the temperature of
the cold BB component from the bulk surface of the NS in the 2BB + PL model.
This is a typical situation when the thermal NS component is equally well
described by the BB and atmosphere models, since the latter spectra are harder
(Potekhin, 2014).

It is important, that all three statistically acceptable models implementing
the $N_{\mathrm{H}}$–$D$ prior require the absorption feature at $\approx
0.35$ keV regardless the local continuum shape. This supports the presence of
the feature in the data. In the 2BB + PL model without the
$N_{\mathrm{H}}$–$D$ relation, which we consider as hardly acceptable, the
absence of the absorption line is compensated by the implausibly high
$N_{\mathrm{H}}$ value. We note that there are only a few rotation powered
pulsars for which absorption lines have been reported: PSR J1740$+$1000
(Kargaltsev et al., 2012), PSR J0659$+$1414 (Arumugasamy et al., 2018;
Zharikov et al., 2021; Schwope et al., 2022), PSR J0726$-$2612 (Rigoselli et
al., 2019), PSR J1819$-$1458 (McLaughlin et al., 2007; Miller et al., 2013)
and Calvera (Shevchuk et al., 2009; Shibanov et al., 2016; Mereghetti et al.,
2021). For none of these objects the nature of lines has been unambiguously
established.

The nature of the J0554 feature is also unclear. One of possible explanations
is a cyclotron absorption line. Such line position measured by a distant
observer is given by:

$E_{\mathrm{cyc}}^{\infty}=11.577\ (1+z_{\mathrm{g}})^{-1}Z\
\frac{m_{\mathrm{e}}}{m}\frac{B}{10^{12}\ \mathrm{G}}\ \mathrm{keV},$ (3)

where $m_{\mathrm{e}}$ is the electron mass, $Z$ and $m$ are the charge number
and mass of the particle that is responsible for the cyclotron absorption.
Hence, the surface magnetic field is $\approx 4\times 10^{10}$ G if the line
is produced by electrons, and $\approx 7\times 10^{13}$ G if it is produced by
protons or $\approx(3-4)\times 10^{13}$ G if it is produced by heavier ions.
The first value is $\sim 200$ times lower and the second is about an order of
magnitude higher than the estimated spin-down magnetic field (see Table 1 and
Appendix B). Thus, the electron cyclotron line can be created at a few stellar
radii above the NS surface, e.g. in a radiation belt, where the magnetic field
is much weaker (see e.g. Luo & Melrose, 2007). Otherwise, if it is an ion
cyclotron line, then the magnetic field at the surface is considerably
stronger than the $B$ values estimated in Appendix B and used in our
atmosphere modelling. Such a cyclotron line might indicate the presence of
strong multipole field components (cf. Bilous et al., 2019; Lockhart et al.,
2019) or magnetic loops (Tiengo et al.,, 2013; Mereghetti, Pons & Melatos,
2015; Rodríguez Castillo et al.,, 2016), which are not considered in our
models. We also note that magnetized atmosphere models predict too low
equivalent width of the proton cyclotron line in comparison with the observed
one (cf. the discussion of a similar case by Hambaryan et al. 2017). In
particular, our dipolar models show that the proton cyclotron line is damped
by the smearing due to magnetic field variations over the surface (see Figs.
13).

Alternatively, the feature may be formed by atomic transitions in a non-
hydrogen atmosphere (e.g., Mori & Ho 2007), the ISM or a cloud near the outer
part of the magnetosphere (cf. Hambaryan et al., 2009; Pires et al., 2019).

Finally, we cannot exclude the possibility that the absorption feature is an
instrumental artefact. The addition of the line improves significantly the fit
of the PN spectrum but does not influence the statistics for the MOS data.
However, this is not surprising since the former spectrum contains much more
counts. To check whether it is the artefact, it would be useful to examine
whether the spectra of other sources in the PN FoV show similar features.
Unfortunately, all other sources are not bright enough to perform such
analysis. Phase-resolved spectral analysis of J0554 could also help to clarify
the nature of the feature in its emission. However, this is impossible because
of the low count statistics.

The non-thermal X-ray luminosity of J0554 in the 2–10 keV band
$L_{\mathrm{X}}\approx(1.6-1.9)\times 10^{30}$ erg s-1 for the
$(\textsc{nsmdintb}+\mathrm{PL})\times\textsc{gabs}$ models and $1.2\times
10^{31}$ erg s-1 for the (2BB + PL) $\times$ gabs model. The corresponding
X-ray efficiencies $\eta_{\mathrm{X}}=L_{\mathrm{X}}/\dot{E}\approx
10^{-4.8}$–$10^{-4.9}$ and $10^{-3.7}$, respectively999For the
$(\textsc{nsmdintb}+\mathrm{PL})\times\textsc{gabs}$ models, we recalculated
the spin-down luminosity using the best-fitting parameters and the formula for
the moment of inertia by Ravenhall & Pethick (1994). The resulting
$\dot{E}=(1.2-1.4)\times 10^{35}$ erg s-1. For the (2BB + PL) $\times$ gabs
model, we used $\dot{E}$ from Table 1.. For a 52 kyr old pulsar these values
are compatible with the dependencies $L_{X}(t_{\mathrm{c}})$ and
$\eta_{\mathrm{X}}(t_{\mathrm{c}})$ based on observations of other X-ray
emitting pulsars (see e.g. Zharikov & Mignani, 2013).

Using the $N_{\mathrm{H}}$–$D$ relation, we found the distance to J0554 to be
1.6–2.4 kpc in the case of the atmosphere models. This is compatible with the
‘pseudo-distance’ of 1.9 kpc based on the $\gamma$-ray data. If this estimate
is correct, the SNR G179.0+2.6 is located somewhat closer than relations
between the radio surface brightness and the diameter of the SNR predict. On
the other hand, the (2BB + PL) $\times$ gabs model resulted in the larger
distance of 2.7–5.5 kpc. This agrees with the upper limit on the distance to
G179.0+2.6 of about 5 kpc provided by How et al. (2018).

Figure 9: J0554 pulse profile in the 0.2–2 keV band (grey dotted histogram
with error bars) and theoretical light curves computed for the model NS133
with $\alpha=\zeta=60^{\circ}$, $k_{\mathrm{B}}T=49$ eV, and for three
different combinations of the mass $M$ and radius $R$ of the NS, drawn using
different line styles as specified in the legend, which cover the range of
pulsed fractions from $\approx 15$ (black solid line) to $\approx 19.5$ per
cent (red dashed line). The inset shows the contours at the 90 per cent
confidence level on the ($\alpha$,$\zeta$) plane for the same three NS models.

We detected, for the first time, X-ray pulsations with the J0554 spin period.
The pulse profile shows two peaks per period. The pulsed fraction in the 0.2–2
keV band is $25\pm 6$ per cent. This is a typical value for the thermal
emission originated from the bulk of an NS surface (e.g. Pavlov & Zavlin,
2000a). Using the pulse profile, we can set additional constraints on the
angles $\alpha$ and $\zeta$ in the case of the atmosphere models. The J0554
observed and theoretical light curves are shown in Fig. 9. The maximum PF
provided by the models nsmdintb in the 0.2–2 keV band is $\approx 20$ per
cent. This is somewhat lower than the measured value but is compatible with it
within uncertainties. The corresponding $\alpha$ and $\zeta$ both lie in the
range of 50∘–70∘ which is compatible with the results of the spectral analysis
(Table 2). Despite the low number of counts, there is a marginal peak in the
2–4 keV pulse profile which coincides with the smaller peak in the 0.2–2 keV
band (see Fig. 5). This may indicate the pulsations of the non-thermal
component. As can be seen from Fig. 7, its flux in the 0.2–2 keV band is
significantly lower than that of the thermal component. Nevertheless, it may
contribute to the pulsed flux and thus somewhat increase the model predicted
PF. It may also be responsible for some asymmetry of the pulse profile. The
inset in Fig. 9 shows the 90 per cent confidence contours of the 2D
distribution of the angles $\alpha$ and $\zeta$ at fixed values of $M$, $R$,
$B$ and $T^{\infty}$. These contours are obtained assuming that the flux
values, derived from the number of counts in the 9 phase bins, are normally
distributed and uncorrelated.

In contrast to the atmosphere models, the light curves for the 2BB model are
almost flat. This is because the BB components assume uniform temperature
distribution over the emitting areas and isotropic radiation (unlike the
peaked radiation from the magnetized atmospheres), but mainly because of the
smallness of the ratio of the fluxes from the hot and cold BB components.
Taking the most probable values for the model (2BB + PL) $\times$ gabs from
Table 2, we obtained the strict upper limit $\mbox{PF}<3.2$ per cent, which is
incompatible with the observed one.

###  6.2 J0554 and NS cooling theories

Figure 10: Theoretical cooling curves, calculated according to the EoS models
BSk24 (panel a) and APR∗ (panel b), compared with observational estimates of
the thermal luminosity and age of J0554 (the error bar). The cooling curves
for NS models with masses $M=1.4$, 1.6, 1.8 and $2.0~{}\mbox{M${}_{\odot}$}$
(coded by color), endowed with a dipole magnetic field having the strength
$B_{\mathrm{p}}=2\times 10^{13}$ G at the pole, are shown by solid and dashed
lines, respectively, for non-accreted and accreted heat blanketing envelopes
(see text for details). For comparison, cooling curves for analogous non-
magnetic NS models are shown, respectively, by dot-dashed and dotted lines.
The inset to panel b illustrates the role of the enhancement of the modified
Urca processes described by Shternin, Baldo & Haensel (2018, SBH). Here, the
solid and long-dashed cooling curves are computed, respectively, with and
without taking the SBH effect into account.

The comparison of the measured thermal luminosity with predictions of NS
cooling theories can provide constraints on the EoS and other properties of
the NS matter and tests for theoretical models of such matter. In Fig. 10 we
demonstrate examples of such comparison. Here, the cooling curves, which show
the time-dependence of the bolometric thermal luminosity of the star in the
reference frame of a distant observer, $L^{\infty}(t^{\infty})$, have been
produced using the computer code presented by Potekhin & Chabrier (2018) with
essentially the same microphysics input. The most uncertain microphysics
ingredients that can significantly affect the cooling are the EoS and
composition of the neutron star core, composition of the heat-blanketing outer
envelopes, and density-dependences of the critical temperatures of the baryon
superfluidity in the core (see, e.g., Potekhin, Pons & Page 2015, for review).
Here we show the cooling curves for two EoSs and for two models of the heat-
blanketing envelopes composition. For the EoS and composition of the core we
employ two widely used models: the model BSk24 (Pearson et al., 2018) and the
model A18+$\delta$v+UIX∗ (Akmal, Pandharipande & Ravenhall, 1998), which is
hereafter named APR∗ and parametrized according to Potekhin & Chabrier (2018).
Both EoSs describe the $npe\mu$ matter, which consists of the leptons and
nucleons without allowance for other baryons or free quarks. For the
composition of the non-accretred crust and envelopes, we use the BSk24 model
of the ground-state matter. The alternative model is the accreted envelope
with helium filling the heat-blanketing layers up to the density $\rho=10^{9}$
g cm-3 (cf. Beznogov, Potekhin & Yakovlev, 2021) and with the ground-state
composition at higher densities. We use the model of a magnetized envelope
with the surface distribution of magnetic field and local effective
temperature consistent with the atmosphere model described in Appendix B, with
the polar magnetic strength $B_{\mathrm{p}}=2\times 10^{13}$ G. We assume that
the field is a relativistic dipole not only in the atmosphere but also in the
stellar core. Accordingly, we take into account the effects of Landau
quantization on the EoS and thermal conductivities in the outer crust and
envelopes and the heat loss due to the synchrotron and fluxoid-scattering
mechanisms of neutrino emission, as well as magnetically induced modifications
of other neutrino emission mechanisms (see Potekhin et al. 2015 and references
therein). For comparison, for each EoS and each envelope composition, we also
plot a cooling curve for the same NS without magnetic field.

The baryon superfluidity in the $npe\mu$ matter of an NS can be of three
types, characterized by the neutron singlet or triplet and proton singlet
pairing. The theory of the neutron singlet superfluidity, which occurs in the
inner crust of an NS, is sufficiently robust (see Ding et al., 2016);
specifically, here we use the MSH model (Margueron et al., 2008) in the
parametrized form by Ho et al. (2015). In contrast, there are several
substantially different theoretical models for the other two types of
superfluidity that operate in the NS core and considerably affect the cooling
(see, e.g., Page et al. 2014; Sedrakian & Clark 2019, for review). For
illustration, we use the parametrizations from Ho et al. (2015) for the BS
(Baldo & Schulze, 2007) and TTav (Takatsuka & Tamagaki, 2004) models of the
proton singlet and neutron triplet superfluidity types, respectively.

The error bar in Fig. 10 embraces the 1$\sigma$ limits to the measured thermal
luminosity $L^{\infty}$ for the atmosphere models (the last two columns of
Table 2). It is placed at $t^{\infty}=t_{\mathrm{c}}$, and the leftward arrow
indicates that the true age of the pulsar is likely somewhat smaller. This
expectation is based on the fact that usually (although not always) true ages
of pulsars are smaller than their characteristic ages (see, e.g., Potekhin et
al. 2020).

Figure 10a shows the cooling curves for the BSk24 model. These cooling curves
pass above the observational error bar for the NS model with
$M=1.4\mbox{M${}_{\odot}$}$, but below it for the models with $M\geq
1.6\mbox{M${}_{\odot}$}$. The fast cooling of the massive NSs is mainly due to
the powerful direct Urca processes of neutrino emission, which operate at
sufficiently high densities in the cores of these stars and overpower the more
common modified Urca processes. This cooling enhancement occurs if $M$ exceeds
a certain threshold value (e.g., Haensel, 1995), which is slightly below
$1.6\mbox{M${}_{\odot}$}$ for the BSk24 model (Pearson et al., 2018). The
observations can be made compatible with the enhanced cooling models, if we
assume that the true age of J0554 is smaller than $t_{\mathrm{c}}$. Among the
cooling curves shown in Fig. 10a, the smallest discrepancy between the true
and characteristic ages is required for the NS with
$M=1.6\mbox{M${}_{\odot}$}$ and accreted envelope. In this case $M$ is only
slightly larger than $M_{\mathrm{DU}}$, so that the direct Urca processes
operate only in a small central part of the NS core.

Figure 10b shows the cooling curves for the APR∗ model. For this EoS, the
direct Urca processes cannot occur for $M\lesssim 2\mbox{M${}_{\odot}$}$.
However, the non-enhanced (so-called minimal) cooling can be compatible with
the observations of J0554 in this case, if the blanketing envelope is non-
accreted. It is worthwhile to note that the latter compatibility is achieved
due to the enhancement of the modified Urca processes, recently discovered by
Shternin, Baldo & Haensel (2018). This enhancement becomes very strong near
the threshold density for opening the direct Urca process. For the APR∗ model,
the threshold mass only slightly exceeds $2\mbox{M${}_{\odot}$}$, therefore
the effect of Shternin et al. (2018) significantly enhances the total neutrino
luminosity of the NS with $M=2\mbox{M${}_{\odot}$}$ and thus decreases its
temperature. In the absence of such enhancement, the luminosity would be
higher for NS models with higher masses in the minimal cooling scenario, as
illustrated by the long-dashed curves in the inset to Fig. 10b, because more
massive stars have larger heat capacities.

The presented analysis is self-consistent in the sense that cooling scenarios
for BSk24 and APR∗ models are considered taking into account the NS parameters
(mass, radius, magnetic field) obtained from the spectral fit. In order to
give preference to one of the EoSs, better constraints on the NS mass and true
age are required.

##  7 Summary

Using 45-ks XMM-Newton observations we detected a soft point-like X-ray source
within 1$\sigma$ area of the Fermi position of the $\gamma$-ray pulsar
J0554+3107, confirming the earlier Swift detection at much higher
significance. We firmly established the pulsar nature of the source by
detecting coherent pulsations with the J0554 frequency in the 0.2–2 keV band.
The pulse profile demonstrates two peaks separated by about a half of the
rotation phase. The background-corrected pulsed fraction is $25\pm 6$ per
cent. Marginal single-peaked pulsations are seen in the hard band above 2 keV,
but low number of counts precludes definite conclusions.

The spectral analysis shows that the thermal emission from the surface of the
NS dominates at energies below $\approx 2$ keV, while a weak non-thermal
magnetospheric component may be present in the hard band. To describe the
former, we constructed a set of the NS hydrogen atmosphere models with dipole
magnetic field. In order to fit the data they require addition of the
absorption line at $\approx 0.35$ keV. The nature of the feature is unclear.
Among the possibilities are cyclotron absorption, atomic transitions in the
ISM or NS atmosphere, or an instrumental artefact. Implementing the absorption
column density–distance relation, we estimated the distance to the pulsar to
be $\approx 2$ kpc. We note that the combination of two blackbody components,
corresponding to the cold NS surface with a hot spot on it, and the absorption
line also provides statistically acceptable fit, resulting in about twice as
large distance. However, this model is less physically realistic and it fails
to reproduce the observed pulse profile with the pulsed fraction $\gtrsim 20$
per cent. Therefore, we consider the atmosphere models as more appropriate
ones.

The best fitting parameters obtained for the atmosphere models suggest that
J0554 is a rather heavy NS with the mass in the range of 1.6–2.1 M⊙ and the
radius of about 13 km. The redshifted effective temperature of $\approx 50$ eV
corresponds to the bolometric luminosity of $\approx 2\times 10^{32}$ erg s-1
as seen by a distant observer. Utilising this value, together with the pulsar
characteristic age of 50 kyr, we investigated cooling scenarios for J0554 in
the frame of two popular EoSs. For the BSk24 model, the observed bolometric
luminosity is consistent with the predictions if the mass of the NS is close
to the lower limit of the range obtained from the spectral fit, and if the
true age of J0554 is substantially smaller than the characteristic one. This
model also favours the presence of the accreted heat-blanketing envelope. On
the other hand, the APR∗ model requires non-accreted envelope and the mass
close to 2.0M⊙, which leads to effective cooling of the NS through modified
Urca processes and compatibility of the model cooling curve with observations
if the true age of J0554 is close to the characteristic one.

Deeper X-ray observations are necessary to better constrain the J0554
parameters. Phase-resolved spectral analysis would allow one to unveil the
nature of the absorption feature and to better constrain the pulsar geometry.
Measurement of the pulsar proper motion could confirm the association of J0554
and SNR G179.0+2.6 and provide an independent estimate of their age.

## Acknowledgements

We are grateful to V. Suleimanov for providing calculations of specific
spectral fluxes from hydrogen atmospheres of NSs with strong magnetic fields.
We also thank the anonymous referee for valuable suggestions that improved the
quality of the paper. DAZ thanks Pirinem School of Theoretical Physics for
hospitality. The work of AST and AYP was supported by the Russian Science
Foundation grant 19-12-00133-P. Scientific results reported in this paper are
based on observations obtained with XMM-Newton, en ESA science mission with
instruments and contributions directly funded by ESA Member States and NASA.

## Data Availability

The XMM-Newton data are available through the data archive
https://www.cosmos.esa.int/web/xmm-newton/xsa.

## References

  * Akmal et al. (1998) Akmal A., Pandharipande V. R., Ravenhall D. G., 1998, Phys. Rev. C, 58, 1804
  * Arnaud (1996) Arnaud K. A., 1996, in Jacoby G. H., Barnes J., eds, Astronomical Society of the Pacific Conference Series Vol. 101, Astronomical Data Analysis Software and Systems V. p. 17
  * Arumugasamy et al. (2018) Arumugasamy P., Kargaltsev O., Posselt B., Pavlov G. G., Hare J., 2018, ApJ, 869, 97
  * Baldo & Schulze (2007) Baldo M., Schulze H. J., 2007, Phys. Rev. C, 75, 025802
  * Becker & Trümper (1999) Becker W., Trümper J., 1999, A&A, 341, 803
  * Bejger & Haensel (2002) Bejger M., Haensel P., 2002, A&A, 396, 917
  * Beznogov et al. (2021) Beznogov M. V., Potekhin A. Y., Yakovlev D. G., 2021, Phys. Rep., 919, 1
  * Bilous et al. (2019) Bilous A. V., et al., 2019, ApJ, 887, L23
  * Buccheri et al. (1983) Buccheri R., et al., 1983, A&A, 128, 245
  * Capitanio et al. (2017) Capitanio L., Lallement R., Vergely J. L., Elyajouri M., Monreal-Ibero A., 2017, A&A, 606, A65
  * Case & Bhattacharya (1998) Case G. L., Bhattacharya D., 1998, ApJ, 504, 761
  * Chen et al. (2014) Chen B. Q., et al., 2014, MNRAS, 443, 1192
  * Chen et al. (2019) Chen B. Q., et al., 2019, MNRAS, 483, 4277
  * Deutsch (1955) Deutsch A. J., 1955, Annales d’Astrophysique, 18, 1
  * Ding et al. (2016) Ding D., Rios A., Dussan H., Dickhoff W. H., Witte S. J., Carbone A., Polls A., 2016, Phys. Rev. C, 94, 025802
  * Drimmel et al. (2003) Drimmel R., Cabrera-Lavers A., López-Corredoira M., 2003, A&A, 409, 205
  * Evans et al. (2014) Evans P. A., et al., 2014, ApJS, 210, 8
  * Foight et al. (2016) Foight D. R., Güver T., Özel F., Slane P. O., 2016, ApJ, 826, 66
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Fuerst & Reich (1986) Fuerst E., Reich W., 1986, A&A, 154, 303
  * Gao et al. (2011) Gao X. Y., Han J. L., Reich W., Reich P., Sun X. H., Xiao L., 2011, A&A, 529, A159
  * Ginzburg & Ozernoy (1965) Ginzburg V. L., Ozernoy L. M., 1965, Sov. Phys. JETP, 20, 689
  * Goodman & Weare (2010) Goodman J., Weare J., 2010, Communications in Applied Mathematics and Computational Science, 5, 65
  * Green et al. (2019) Green G. M., Schlafly E., Zucker C., Speagle J. S., Finkbeiner D., 2019, ApJ, 887, 93
  * Grießmeier et al. (2021) Grießmeier J. M., Smith D. A., Theureau G., Johnson T. J., Kerr M., Bondonneau L., Cognard I., Serylak M., 2021, A&A, 654, A43
  * Guseinov et al. (2003) Guseinov O. H., Ankay A., Sezer A., Tagieva S. O., 2003, Astronomical and Astrophysical Transactions, 22, 273
  * Haensel (1995) Haensel P., 1995, Space Sci. Rev., 74, 427
  * Hambaryan et al. (2009) Hambaryan V., Neuhäuser R., Haberl F., Hohle M. M., Schwope A. D., 2009, A&A, 497, L9
  * Hambaryan et al. (2017) Hambaryan V., Suleimanov V., Haberl F., Schwope A. D., Neuhäuser R., Hohle M., Werner K., 2017, A&A, 601, A108
  * Ho et al. (2015) Ho W. C. G., Elshamouty K. G., Heinke C. O., Potekhin A. Y., 2015, Phys. Rev. C, 91, 015806
  * How et al. (2018) How T. G., Fesen R. A., Neustadt J. M. M., Black C. S., Outters N., 2018, MNRAS, 478, 1987
  * Kargaltsev & Pavlov (2008) Kargaltsev O., Pavlov G. G., 2008, in Bassa C., Wang Z., Cumming A., Kaspi V. M., eds, American Institute of Physics Conference Series Vol. 983, 40 Years of Pulsars: Millisecond Pulsars, Magnetars and More. pp 171–185 (arXiv:0801.2602), doi:10.1063/1.2900138
  * Kargaltsev et al. (2012) Kargaltsev O., Durant M., Misanovic Z., Pavlov G. G., 2012, Science, 337, 946
  * Lallement et al. (2014) Lallement R., Vergely J. L., Valette B., Puspitarini L., Eyer L., Casagrande L., 2014, A&A, 561, A91
  * Lallement et al. (2018) Lallement R., et al., 2018, A&A, 616, A132
  * Lockhart et al. (2019) Lockhart W., Gralla S. E., Özel F., Psaltis D., 2019, MNRAS, 490, 1774
  * Lower et al. (2021) Lower M. E., et al., 2021, MNRAS, 508, 3251
  * Luo & Melrose (2007) Luo Q., Melrose D., 2007, MNRAS, 378, 1481
  * Manchester et al. (2005) Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993
  * Margueron et al. (2008) Margueron J., Sagawa H., Hagino K., 2008, Phys. Rev. C, 77, 054309
  * McLaughlin et al. (2007) McLaughlin M. A., et al., 2007, ApJ, 670, 1307
  * Mereghetti et al. (2015) Mereghetti S., Pons J. A., Melatos A., 2015, Space Sci. Rev., 191, 315
  * Mereghetti et al. (2021) Mereghetti S., Rigoselli M., Taverna R., Baldeschi L., Crestan S., Turolla R., Zane S., 2021, ApJ, 922, 253
  * Miller et al. (2013) Miller J. J., McLaughlin M. A., Rea N., Lazaridis K., Keane E. F., Kramer M., Lyne A., 2013, ApJ, 776, 104
  * Mori & Ho (2007) Mori K., Ho W. C. G., 2007, MNRAS, 377, 905
  * Page et al. (2014) Page D., Lattimer J. M., Prakash M., Steiner A. W., 2014, in Bennemann K.-H., Ketterson J. B., eds, International Series of Monographs on Physics Vol. 2, Novel Superfluids. pp 505–579
  * Pavlov & Zavlin (2000a) Pavlov G. G., Zavlin V. E., 2000a, in Martens P. C. H., Tsuruta S., Weber M. A., eds, Proceedings of IAU Symposium Vol. 195, Highly Energetic Physical Processes and Mechanisms for Emission from Astrophysical Plasmas. p. 103
  * Pavlov & Zavlin (2000b) Pavlov G. G., Zavlin V. E., 2000b, ApJ, 529, 1011
  * Pavlovic et al. (2014) Pavlovic M. Z., Dobardzic A., Vukotic B., Urosevic D., 2014, Serbian Astronomical Journal, 189, 25
  * Pearson et al. (2018) Pearson J. M., Chamel N., Potekhin A. Y., Fantina A. F., Ducoin C., Dutta A. K., Goriely S., 2018, MNRAS, 481, 2994
  * Pétri (2019) Pétri J., 2019, MNRAS, 485, 4573
  * Pires et al. (2019) Pires A. M., Schwope A. D., Haberl F., Zavlin V. E., Motch C., Zane S., 2019, A&A, 623, A73
  * Pletsch et al. (2013) Pletsch H. J., et al., 2013, ApJ, 779, L11
  * Potekhin (2014) Potekhin A. Y., 2014, Physics Uspekhi, 57, 735
  * Potekhin & Chabrier (2003) Potekhin A. Y., Chabrier G., 2003, ApJ, 585, 955
  * Potekhin & Chabrier (2018) Potekhin A. Y., Chabrier G., 2018, A&A, 609, A74
  * Potekhin et al. (2003) Potekhin A. Y., Yakovlev D. G., Chabrier G., Gnedin O. Y., 2003, ApJ, 594, 404
  * Potekhin et al. (2004) Potekhin A. Y., Lai D., Chabrier G., Ho W. C. G., 2004, ApJ, 612, 1034
  * Potekhin et al. (2014) Potekhin A. Y., Chabrier G., Ho W. C. G., 2014, A&A, 572, A69
  * Potekhin et al. (2015) Potekhin A. Y., Pons J. A., Page D., 2015, Space Sci. Rev., 191, 239
  * Potekhin et al. (2020) Potekhin A. Y., Zyuzin D. A., Yakovlev D. G., Beznogov M. V., Shibanov Y. A., 2020, MNRAS, 496, 5052
  * Ravenhall & Pethick (1994) Ravenhall D. G., Pethick C. J., 1994, ApJ, 424, 846
  * Rigoselli et al. (2019) Rigoselli M., Mereghetti S., Suleimanov V., Potekhin A. Y., Turolla R., Taverna R., Pintore F., 2019, A&A, 627, A69
  * Rodríguez Castillo et al. (2016) Rodríguez Castillo G. A., et al., 2016, MNRAS, 456, 4145
  * Sale et al. (2014) Sale S. E., et al., 2014, MNRAS, 443, 2907
  * Saz Parkinson et al. (2010) Saz Parkinson P. M., et al., 2010, ApJ, 725, 571
  * Schwope et al. (2022) Schwope A., et al., 2022, A&A, 661, A41
  * Sedrakian & Clark (2019) Sedrakian A., Clark J. W., 2019, European Physical Journal A, 55, 167
  * Shevchuk et al. (2009) Shevchuk A. S. H., Fox D. B., Rutledge R. E., 2009, ApJ, 705, 391
  * Shibanov et al. (2016) Shibanov Y., Danilenko A., Zharikov S., Shternin P., Zyuzin D., 2016, ApJ, 831, 112
  * Shternin et al. (2018) Shternin P. S., Baldo M., Haensel P., 2018, Physics Letters B, 786, 28
  * Spitkovsky (2006) Spitkovsky A., 2006, ApJ, 648, L51
  * Sturrock (1971) Sturrock P. A., 1971, ApJ, 164, 529
  * Suleimanov et al. (2009) Suleimanov V., Potekhin A. Y., Werner K., 2009, A&A, 500, 891
  * Swanepoel et al. (1996) Swanepoel J. W. H., de Beer C. F., Loots H., 1996, ApJ, 467, 261
  * Takatsuka & Tamagaki (2004) Takatsuka T., Tamagaki R., 2004, Progress of Theoretical Physics, 112, 37
  * Tiengo et al. (2013) Tiengo A., et al., 2013, Nature, 500, 312
  * Tyul’bashev et al. (2021) Tyul’bashev S. A., Kitaeva M. A., Tyulbasheva G. E., 2021, Astronomy Reports, 65, 819
  * Wachter et al. (1979) Wachter K., Leach R., Kellogg E., 1979, ApJ, 230, 274
  * Willatt & Ehle (2016) Willatt R., Ehle M., 2016, Guide for use of the images script. https://www.cosmos.esa.int/documents/332006/641121/README.pdf
  * Wilms et al. (2000) Wilms J., Allen A., McCray R., 2000, ApJ, 542, 914
  * Yakovlev & Pethick (2004) Yakovlev D. G., Pethick C. J., 2004, ARA&A, 42, 169
  * Zhao et al. (2020) Zhao H., Jiang B., Li J., Chen B., Yu B., Wang Y., 2020, ApJ, 891, 137
  * Zharikov & Mignani (2013) Zharikov S., Mignani R. P., 2013, MNRAS, 435, 2227
  * Zharikov et al. (2021) Zharikov S., Zyuzin D., Shibanov Y., Kirichenko A., Mennickent R. E., Geier S., Cabrera-Lavers A., 2021, MNRAS, 502, 2005
  * Zyuzin et al. (2018) Zyuzin D. A., Karpova A. V., Shibanov Y. A., 2018, MNRAS, 476, 2177
  * Zyuzin et al. (2021) Zyuzin D. A., Karpova A. V., Shibanov Y. A., Potekhin A. Y., Suleimanov V. F., 2021, MNRAS, 501, 4998

##  Appendix A Pulsations detection significance and spin frequency
uncertainty

To estimate the confidence level of the detection of the periodic signal from
J0554, we generated 1 million synthetic light curves with the length and the
mean count rate equal to the ones observed from J0554, but consisting of pure
Poisson noise without any periodic component. For each of the light curves, we
performed $Z_{2}^{2}$-test using the same frequency window and the same number
of trial frequencies as for $Z_{2}^{2}$-test on the observed light curve. The
obtained highest $Z_{2}^{2}$ values were used to construct the cumulative
distribution function (CDF) of $Z_{2}^{2}$ in the absence of the periodic
signal (see Fig. 11). We found that the probability to get $Z_{2}^{2}=42.7$
from pure noise is about $2\times 10^{-6}$, which corresponds to the detection
confidence level of $\approx 4.7\sigma$.

The frequency uncertainty was estimated in a similar fashion. We fitted the
pulse profile with the sum of the first two harmonics. We simulated 1000 event
lists of the periodic signal with the measured frequency, amplitudes and
relative phases of two harmonics, keeping the mean count rate fixed and
varying the number of photons and their times of arrival according to the
Poisson statistics. For each event list, we found the frequency of the signal
performing $Z_{2}^{2}$-test identical to the one applied to the real data. The
nearly symmetrical 68 per cent confidence interval of the resulting frequency
distribution (see Fig. 12) was taken as the desired uncertainty.

Figure 11: Cumulative distribution function of $Z_{2}^{2}$ for a pure noise
signal with the mean count rate equal to the observed from J0554. The vertical
dashed line corresponds to $Z_{2}^{2}=42.7$.

Figure 12: Distribution of the best frequencies measured for 1000 simulated
event lists with the parameters fixed at the values obtained for J0554. The
solid vertical line corresponds to the most probable value, while the shaded
area indicates the 68 per cent credible interval which was taken as the error
of the frequency determination.

##  Appendix B Atmosphere models

Figure 13: Top panel: Thermal spectra of an NS model with
$M=1.4\mbox{M${}_{\odot}$}$, $R=14$ km and $B_{\mathrm{p}}=10^{13}$ G, as seen
by a distant observer, for $\log T_{\mathrm{eff}}^{\infty}\mbox{\,[K]}=5.4$,
5.5, 5.6, 5.7 and 5.8 (the five bunches of lines from bottom to top), and
seven values of the angle between the magnetic axis and the line of sight,
$\Theta_{\mathrm{m}}$ (different line styles and colors according to the
legend). The arrows indicate approximate positions of the features due to the
proton cyclotron absorption (pc), the principal bound-bound transition (bb)
and the principal photoionization threshold (bf), smeared by the field
distribution over the surface and/or by atomic motion effects. For comparison,
the dotted lines show the blackbody spectra for the same NS parameters and
redshifted effective temperatures $T^{\infty}$. Bottom panel: The same as in
the top panel but for an NS model with $M=2\mbox{M${}_{\odot}$}$ and $R=10$
km.

The computation of the spectrum that can be measured by a distant observer is
patterned after Zyuzin et al. (2021). We construct the integral spectrum by
assembling local spectra at different patches on the surface. We assume a
dipolar magnetic field, modified by the effects of General Relativity
(Ginzburg & Ozernoy, 1965; Pavlov & Zavlin, 2000b). The temperature
distribution, which is associated with this magnetic field, is calculated
following Potekhin et al. (2003).

For every selected field strength at the magnetic pole $B_{\mathrm{p}}$ and
selected NS mass and radius, the local radiative flux density was computed at
three magnetic latitudes, including the pole, for a set of 480 directions of
the photon wave vector and for 150–200 photon energies in the X-ray band,
using an advanced version of the code described in Suleimanov, Potekhin &
Werner (2009). The fourth latitude is the equator, which is too cold to allow
construction of an atmosphere model with the currently available opacities in
strong magnetic fields for the selected range of effective temperatures
log$T^{\infty}$ [K] = 5.4–5.8. However, its contribution to the total flux is
small, so we replace it by the blackbody model (we have checked that using
alternative models does not lead to a noticeable change in the total
spectrum).

The code of Suleimanov et al. (2009) has been modified to allow for different
angles $\theta_{B}$ between the magnetic field and the normal to the surface.
Hydrogen composition is considered, taking into account incomplete ionization.
The effects of the strong magnetic field and the atomic thermal motion across
the field on the plasma opacities are treated following Potekhin & Chabrier
(2003) with the improvements described in Potekhin, Chabrier & Ho (2014).
Polarization vectors and opacities of normal electromagnetic modes are
calculated as in Potekhin et al. (2004). Then flux values at arbitrary
latitudes, energies and directions are obtained by interpolation (or
extrapolation, whenever needed).

The monochromatic spectral flux density measured by a distant observer is
computed by integrating the emission from different local patches over the
stellar surface for any selected angle $\Theta_{\mathrm{m}}$ between the
magnetic dipole axis and the line of sight (see Appendix A of Zyuzin et al.
2021 for details). In the axisymmetric model, the pulsar geometry is
determined by the angles $\alpha$ and $\zeta$ that the spin axis makes with
the magnetic axis and with the line of sight, respectively (e.g., Pavlov &
Zavlin, 2000b). To produce phase-resolved spectra, it is sufficient to
calculate
$\cos\Theta_{\mathrm{m}}=\sin\zeta\sin\alpha\cos\phi+\cos\alpha\cos\zeta$ for
each rotation phase $\phi$. Then the light curve and the phase-integrated
spectrum are given by integration of the phase-resolved spectrum over the
energy or over the phase $\phi$, respectively.

For isolated pulsars, the widely used estimate of the magnetic field strength
is based on the expression

$B\approx 3.2\times 10^{19}\,C\,\sqrt{P\dot{P}}\textrm{~{}~{}G},$ (4)

where $P$ is the period in seconds, $\dot{P}$ is the period time derivative,
and $C$ is a coefficient, which depends on stellar parameters. For the non-
relativistic rotating magnetic dipole in vacuo (Deutsch, 1955), the magnetic
field strength at the equator $B_{\mathrm{eq}}$ is given by setting

$C=R_{10}^{-3}\,(\sin\alpha)^{-1}\,\sqrt{I_{45}},$ (5)

where $R_{10}\equiv R/(10\mbox{~{}km})$ and $I_{45}$ is the moment of inertia
in units of $10^{45}$ g cm2. The latter depends on the EoS, but in most
plausible settings it can be estimated with an accuracy within 10 per cent by
the approximation of Ravenhall & Pethick (1994), which can be written in the
form

$I_{45}\approx 0.42(M/\mbox{M${}_{\odot}$})(R^{\infty}/10\mbox{~{}km})^{2}$
(6)

(see Bejger & Haensel 2002 for a more general fitting formula). The
characteristic field $B_{\mathrm{c}}$ (Table 1) is defined by equation (4)
with $C=1$ (e.g., Manchester et al., 2005). For the likely values of
$\alpha\gtrsim 50^{\circ}$, $M\approx(1.6$–2.1) M⊙ and $R\approx 11.7$–14.7
km, implied by the fitting results in Table 2, equations (4)–(6) give
$B_{\mathrm{eq}}\sim(4$–$11)\times 10^{12}$ G, which implies, for the non-
relativistic dipole field, $B_{\mathrm{p}}\sim(0.8$–$2.2)\times 10^{13}$ G.
These values are consistent with the values $B_{\mathrm{p}}=10^{13}$ G and
$B_{\mathrm{p}}=2\times 10^{13}$ G that we used to construct the atmosphere
models. We have also tried models with lower and higher field strengths, but
found that they do not provide a better fit.

A real pulsar differs from a rotating magnetic dipole, because its
magnetosphere is filled with plasma, carrying electric charges and currents.
According to the results of numerical simulations of plasma behavior in the
pulsar magnetosphere (Spitkovsky, 2006), the equatorial magnetic field can be
approximately described by equation (4) with

$C\approx(0.8\pm 0.1)R_{10}^{-3}\,(1+\sin^{2}\alpha)^{-1/2}\,\sqrt{I_{45}},$
(7)

which gives estimates in the range of $B_{\mathrm{p}}\sim(3$–$9)\times
10^{12}$ G. Additional uncertainties arise from the effects of General
Relativity, pulsar wind and deviations from the pure dipole geometry (see
Pétri 2019 for discussion and references).

For each of the selected $B_{\mathrm{p}}$ values, we have considered $M$, $R$,
$T^{\infty}$, $\alpha$ and $\zeta$ as fitting parameters, using interpolation
and extrapolation based on the computed spectra for
$M=1.4\mbox{M${}_{\odot}$}$ and 2.0M⊙, $R=10$ km, 12 km and 14 km, $\log
T^{\infty}\mbox{\,[K]}=5.4$, 5.5, 5.6, 5.7 and 5.8, and various
$\Theta_{\mathrm{m}}$. Examples of the computed spectra are shown in the top
and bottom panels of Fig. 13 for less and more compact NS models,
respectively.