# Atmospheric Effects on Neutron Star Parameter Constraints with NICER

Tuomo Salmi, Serena Vinciguerra, Devarshi Choudhury, Anna L Watts, Wynn C.G Ho, Sebastien Guillot, Yves Kini, Bas Dorsman, Sharon M Morsink, Slavko Bogdanov

► To cite this version:

Tuomo Salmi, Serena Vinciguerra, Devarshi Choudhury, Anna L Watts, Wynn C.G Ho, et al. Atmospheric Effects on Neutron Star Parameter Constraints with NICER. The Astrophysical Journal, 2023, 956 (2), pp.138. 10.3847/1538-4357/acf49d. hal-04196315

HAL Id: hal-04196315

https://hal.science/hal-04196315v1

Submitted on 9 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Atmospheric Effects on Neutron Star Parameter Constraints with NICER

Tuomo Salmi<fnref n="1" /> ID, Serena Vinciguerra<fnref n="1" /> ID, Devarshi Choudhury<fnref n="1" /> ID, Anna L. Watts<fnref n="1" /> ID, Wynn C. G. Ho<fnref n="2" /> ID, Sebastien Guillot<fnref n="3" /> ID, Yves Kini<fnref n="1" /> ID, Bas Dorsman<fnref n="1" /> ID, Sharon M. Morsink<fnref n="4" /> ID, and Slavko Bogdanov<fnref n="5" /> ID

Received 2023 March 23; revised 2023 July 14; accepted 2023 July 15; published 2023 October 17

## Abstract

We present an analysis of the effects of uncertainties in the atmosphere models on the radius, mass, and other neutron star parameter constraints for the NICER observations of rotation-powered millisecond pulsars. To date, NICER has applied the X-ray pulse profile modeling technique to two millisecond-period pulsars: PSR J0030+0451 and the high-mass pulsar PSR J0740+6620. These studies have commonly assumed a deep-heated, fully ionized hydrogen atmosphere model, although they have explored the effects of partial ionization and helium composition in some cases. Here, we extend that exploration and also include new models with partially ionized carbon composition, externally heated hydrogen, and an empirical atmospheric beaming parameterization to explore deviations in the expected anisotropy of the emitted radiation. None of the studied atmosphere cases have any significant influence on the inferred radius of PSR J0740+6620, possibly due to its X-ray faintness, tighter external constraints, and/or viewing geometry. In the case of PSR J0030+0451, both the composition and ionization state could significantly alter the inferred radius. However, based on the evidence (prior predictive probability of the data), partially ionized hydrogen and carbon atmospheres are disfavored. The difference in the evidence for ionized hydrogen and helium atmospheres is too small to be decisive for most cases, but the inferred radius for helium models trends to larger sizes around or above 14–15 km. External heating or deviations in the beaming that are less than 5% at emission angles smaller than 60°, on the other hand, have no significant effect on the inferred radius.

*Unified Astronomy Thesaurus concepts: Neutron stars (1108); X-ray astronomy (1810)*

*Supporting material: figure sets*

## 1. Introduction

Mass and radius measurements of neutron stars (NSs) can be used to probe the strongly degenerate matter at supranuclear densities in NS cores. This is because the NS mass–radius relation depends sensitively on the pressure–density relation, i.e., the equation of state of the interior matter (see, e.g., Lattimer & Prakash 2001; Özel & Psaltis 2009; Lattimer 2012; Hebeler et al. 2013; Baym et al. 2018). One useful approach to measure the masses and radii of different NSs is by modeling the observed X-ray pulses produced by hot spots on the surfaces of rotating NSs. This extensively studied technique exploits the general and special relativistic effects modifying the pulse shape depending on mass, radius, and other model parameters (see, e.g., Pechenick et al. 1983; Strohmayer 1992; Miller & Lamb 1998; Poutanen & Gierliński 2003; Morsink et al. 2007; Lo et al. 2013; AlGendy & Morsink 2014; Miller & Lamb 2015; Watts et al. 2016; Nättilä & Pihajoki 2018; Bogdanov et al. 2019b).

In recent years, the pulse profile modeling method has been applied by NASA’s Neutron Star Interior Composition Explorer (NICER; Gendreau et al. 2016), concentrating on rotation-powered millisecond pulsars (RMPs). In these pulsars,

the hot regions on the NS surface are heated by the flow-back of electrons or positrons that are accelerated in the magnetosphere (in contrast to, e.g., hot spots in accreting millisecond pulsar, which are heated due to the accretion from a companion star). The bombarding particles are produced in single photon magnetic pair cascades originating at the open field line region (see, e.g., Ruderman & Sutherland 1975; Arons 1981; Harding & Muslimov 2001). So far, results for two pulsars have been reported: PSR J0030+0451 (Miller et al. 2019; Riley et al. 2019; Vinciguerra et al. 2023b) and PSR J0740+6620 (Miller et al. 2021; Riley et al. 2021; Salmi et al. 2022). These studies have started to restrict dense matter models (see, e.g., Raaijmakers et al. 2021) and pointed to nonantipodal magnetic polar cap geometries (see, e.g., Bilous et al. 2019).

In the studies mentioned above, the energy and emission angle dependence (i.e., beaming pattern) of the radiation escaping from the NS surface is calculated from the input physics (for a given effective temperature<fnref n="6" /> and surface gravity) based mainly on the NSX deep-heated, nonmagnetic, and fully ionized hydrogen atmosphere model (Ho & Lai 2001). Sensitivity to chemical composition was tested using the NSX fully ionized helium atmosphere for PSR J0030+0451 in Miller et al. (2019) and for PSR J0740+6620 in Riley et al. (2021). In addition, the

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹ Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098XH Amsterdam, The Netherlands; t.h.j.salmi@uva.nl</footnote>

<footnote n="2">² Department of Physics and Astronomy, Haverford College, 370 Lancaster Avenue, Haverford, PA 19041, USA</footnote>

<footnote n="3">³ Institut de Recherche en Astrophysique et Planétologie, UPS-OMP, CNRS, CNES, 9 avenue du Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France</footnote>

<footnote n="4">⁴ Department of Physics, University of Alberta, 4-183 CCIS, Edmonton, AB, T6G 2E1, Canada</footnote>

<footnote n="5">⁵ Columbia Astrophysics Laboratory, Columbia University, 550 West 120th Street, New York, NY 10027, USA</footnote>

<footnote>Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>

<footnote n="6">⁶ Effective temperature $T_{\text{eff}}$ is the temperature of a blackbody producing a bolometric photon flux equal to the bolometric photon flux from the atmosphere model.</footnote>
sensitivity to NSX with partially ionized hydrogen was checked for PSR J0740+6620 in Miller et al. (2021). The alternative atmosphere cases resulted in consistent mass and radius values, except in the case of helium for PSR J0030+0451, which was discarded because of the very high inferred NS mass around 2.7 $M_⊙$, with most of the posterior probability above 2.0 $M_⊙$ (Miller et al. 2019). In this work, we explore additional cases with a helium atmosphere and show that for PSR J0030+0451 well-fitting solutions with helium atmospheres can be found also with more realistic NS masses (see discussion in Section 4.2).

The uppermost atmospheric layers, which determine the properties of escaping radiation, are usually expected to consist of hydrogen due to the rapid sinking of heavier elements via diffusive gravitational separation (Alcock & Illarionov 1980; Hameury et al. 1983; Brown et al. 2002; Zavlin & Pavlov 2002) on a timescale of ~1–100 s (Romani 1987). However, they could consist of helium if all the hydrogen originally present has been converted to helium via diffusive nuclear burning (Chang & Bildsten 2004; Wijngaarden et al. 2019), or if the accretion that recycled the RMP happened from a completely hydrogen-depleted companion star (see Bogdanov et al. 2021, and references therein). Less likely, but still possible, would be an atmosphere consisting of even heavier elements, e.g., due to nuclear burning, a wind exposing the underlying heavier elements, or a lack of light element accretion (Chang & Bildsten 2004; Bogdanov et al. 2021). We therefore also perform some new analysis using a partially ionized carbon NSX atmosphere model.

So far, the inferred effective temperatures in RMPs analyzed by NICER have been around $log_{10}(T_{eff}/K) \sim 6$ for which more than 99% of the hydrogen is completely ionized at all depths (see Figure 4 of Zavlin et al. 1996). However, neutral atoms could still exist in dense enough layers of the atmosphere, or there could be temperature layers cold enough to produce small bound-bound and bound-free opacity features in the escaping radiation (Zavlin et al. 1996). Therefore, we have also explored more partially ionized models. The caveat, however, is that the opacity tables needed for partially ionized models do not cover the entire range of energies, temperatures, and densities that are relevant in RMPs (Iglesias & Rogers 1996; Badnell et al. 2005; Colgan et al. 2016; Bogdanov et al. 2021).

Besides the composition and the ionization state, a few other possible sources of uncertainties in the atmosphere models were also discussed in Bogdanov et al. (2021; for a review see also, e.g., Potekhin 2014). These include the magnetic field strength, depth of energy deposition of the bombarding particles, and the accuracy of the computational methods used in atmosphere modeling. These effects were qualitatively predicted to produce deviations that are small for the NICER sources analyzed so far, but we quantify this for the cases of deposition depth and computational methods by performing pulse profile analyses using an independent numerical implementation of the atmosphere model, assuming an externally heated, fully ionized hydrogen atmosphere from Salmi et al. (2020; see also Bauböck et al. 2019). In Bogdanov et al. (2021), comparisons were previously done against another atmosphere code called McPHAC (Haakonsen et al. 2012), identifying a couple of issues in using the latter. Finally, in order to allow small deviations from the nominal beaming in any of the atmosphere models, we also study new pulse profile models with additional free parameters that can modify the beaming.

The remainder of this paper is structured as follows. In Section 2, we present the methods that we use to estimate the atmospheric uncertainties in pulse profile modeling and parameter constraints. We introduce the different model atmospheres and beaming correction formalism that are applied in our analysis. In addition, we compare the spectra and beaming patterns calculated using different atmosphere models. In Section 3, we present the parameter constraints obtained with different atmosphere model assumptions, considering first PSR J0740+6620 and then PSR J0030+0451. We discuss the implications of our results for future analysis in Section 4 and conclude in Section 5.

## 2. Modeling Procedure

In this Section, we present the methods we use to produce the X-ray pulses and fit the data. We focus first on the new features applied for the atmosphere modeling, and then summarize the other modeling aspects, which are mainly shared with those from Bogdanov et al. (2019b), Riley et al. (2019, 2021), Salmi et al. (2022), Vinciguerra et al. (2023a), and Vinciguerra et al. (2023b). For all posterior computations, we use the X-ray Pulse Simulation and Inference<fnref n="7" /> (X-PSI) code, with versions ranging from v0.7.3 to v1.2.1 (Riley et al. 2023).<fnref n="8" /> The exact version for each run, data products, posterior sample files, a complete set of output figures, and all of the analysis files in the Python language may be found in the persistent Zenodo repository of Salmi et al. (2023) at doi: 10.5281/zenodo.7449786. We use a new atmosphere extension that was developed for the free beaming runs and introduced in v1.2.0 of X-PSI; in Salmi et al. (2023), we also provide instructions on how to install that for the previous versions of X-PSI to allow reproducibility of this work.

### 2.1. Atmosphere Modeling

We employ atmosphere models to determine the specific intensity $I(E, \mu)$ emitted from the NS surface, where $E$ is the photon energy and $\mu$ is the cosine of the emission angle with respect to the surface normal. Since calculating this quantity is computationally expensive, we use pre-computed look-up tables of $I(E, \mu)$ when calculating the pulse profiles, similar to the previous NICER studies. The tables used in this study are presented next in Section 2.1.1, and the procedure for allowing deviation from the precalculated beaming pattern is described in Section 2.1.2.

#### 2.1.1. Numerical Model Atmospheres

We use five different look-up tables for the numerical atmosphere models. The first is the deep-heated and fully ionized hydrogen NSX model (nsxH; Ho & Lai 2001), used in the previous X-PSI analyses. This table provides the logarithm of the ratio of the specific intensity with respect to the cube of the effective temperature, $(I/T_{eff}^3)$, for a grid in the logarithm of the effective temperature, $log_{10}(T_{eff}/K)$, the logarithm of the surface gravity, $log_{10}(g_s/cm s^{-2})$, the cosine of the emission angle, $\mu$, and the logarithm of photon energy with respect to the temperature, $log_{10}(E/kT_{eff})$ (see Bogdanov et al. 2021, for more details). The ranges and the number of grid points for each

<footnote n="7">⁷ https://github.com/xpsi-group/xpsi</footnote>

<footnote n="8">⁸ The versions are practically identical for the considered models; the only actual difference is the fix of a numerical ray tracing issue since v0.7.12, affecting only a few parameter vectors.</footnote>
<table><caption>**Table 1**
Summary of the Numerical Atmosphere Tables</caption><thead><tr><th>Model</th><th>log<sub>10</sub>(T<sub>eff</sub>/K)</th><th>log<sub>10</sub>(g<sub>x</sub>/cm s<sup>-2</sup>)</th><th>log<sub>10</sub>(E/kT<sub>eff</sub>)</th><th>μ</th></tr></thead><tbody><tr><td>hatmH</td><td>[5.1, 6.8]: 18</td><td>[13.7, 14.6]: 10</td><td>[-1.3, 2.0]: 166<sup>a</sup></td><td>[0.0, 1.0]: 11</td></tr><tr><td>nsxCp</td><td>[5.6, 6.7]: 23</td><td>[13.6, 15.0]: 15</td><td>[-1.3, 2.0]: 166</td><td>[1 × 10<sup>-6</sup>, 1 − 1 × 10<sup>-6</sup>]: 67</td></tr><tr><td>nsxHp</td><td>[5.2, 6.7]: 31</td><td>[13.6, 15.0]: 15</td><td>[-1.3, 2.0]: 166</td><td>[1 × 10<sup>-6</sup>, 1 − 1 × 10<sup>-6</sup>]: 67</td></tr><tr><td>nsxHe</td><td>[5.1, 6.5]: 29</td><td>[13.7, 14.7]: 11</td><td>[-1.3, 2.0]: 166</td><td>[1 × 10<sup>-6</sup>, 1 − 1 × 10<sup>-6</sup>]: 67</td></tr><tr><td>nsxH</td><td>[5.1, 6.8]: 35</td><td>[13.7, 15.0]: 14</td><td>[-1.3, 2.0]: 166</td><td>[1 × 10<sup>-6</sup>, 1 − 1 × 10<sup>-6</sup>]: 67</td></tr></tbody></table>

**Notes.** The atmosphere grid ranges ([$\cdots$, $\cdots$]) and number of grid points ($\cdots$) are shown for fully ionized externally heated hydrogen (hatmH), partially ionized carbon (nsxCp), partially ionized hydrogen (nsxHp), fully ionized helium (nsxHe), and fully ionized hydrogen (nsxH). For NSX models, emission angles were placed at every 1°.5, with three additional values near $\mu = 0$, and $\mu = 1$. In the case of hatmH, the grid points were based on the sample points of Gauss-Legendre quadrature between 0 and 1.

<sup>a</sup> Originally a logarithmically spaced grid of E from 4.1 × 10⁻⁴ to 41 keV with 360 grid points was used for calculating this model, but the intensity values were then interpolated to the same log₁₀(E/kT<sub>eff</sub>) grid as in NSX, for simplicity.

atmosphere model parameter are shown in Table 1. The same information is also displayed for the other atmosphere tables used in this paper.

The second look-up table is the deep-heated and fully ionized helium NSX model (nsxHe; Ho & Lai 2001), included also, e.g., in the analysis of Riley et al. (2021). This table is otherwise similar to the hydrogen table, except for having slightly smaller upper limits for high temperature and surface gravity.

For the third and fourth look-up tables, we use the deep-heated partially ionized hydrogen NSX model (nsxHp), applied previously in Miller et al. (2021),<fnref n="9" /> and the partially ionized carbon NSX model, (nsxCp, which is essentially the same as in Ho & Heinke 2009, except for different grids used in the calculation). These tables have roughly similar ranges as the previous ones, as seen in Table 1, although the carbon model does not reach such cold temperatures as the others.

For the fifth look-up table, we use an externally heated and fully ionized hydrogen model (hatmH) from Salmi et al. (2020). This model is mostly similar to NSX, but it allows the bombarding return-current particles to release their kinetic energy as heat at different depths of the atmosphere instead of assuming that all the heat is released at very deep layers. For RMPs, deep heating is expected to be a reasonable assumption, as explained in Section 4.4, but we nevertheless check if relaxing this assumption (adopting an extreme case) might affect the results. In addition to external heating, the new model calculates the photon-electron scattering using the exact Compton scattering formalism and a temperature iteration scheme similar to those in Suleimanov et al. (2012), Salmi et al. (2019), and Suleimanov et al. (2020), which differ slightly from the Thomson scattering and iteration schemes used in NSX.

In this paper, the hatmH atmosphere was assumed to be heated by particles with a power-law energy distribution, characterized by a minimum Lorentz factor $\gamma_{min} = 10$, power-law slope $\delta = 2$, and a maximum Lorentz factor $\gamma_{max} = 2 \times 10^5$. These choices correspond to the fiducial parameter setup in Salmi et al. (2020), and they lead to an average Lorentz factor $\gamma_{avg} = 100$ of a bombarding particle. These values were selected to maximize the difference with the deep-heated atmosphere model, while not being incompatible with our

understanding of pulsar physics, given the uncertainties in simulations (see discussion in Section 4.4). Due to the computationally more expensive atmosphere model, fewer grid points were produced for the hatmH table in effective temperature, surface gravity, and emission angle, as seen in Table 1.

### 2.1.2. Beaming Correction Formalism and Priors

In some of our models (see Tables 2 and 3 for details), we introduce additional freedom in the beaming pattern predicted by the numerical atmosphere models. This is done by multiplying the intensity values $I(E, \mu)_{NUM}$ (interpolated from the atmosphere table) with a polynomial function modifying the beaming:

$$I(E, \mu, a, b, c, d) = I(E, \mu)_{NUM} f(E, \mu, a, b, c, d), \quad (1)$$

where

$$f(E, \mu, a, b, c, d) = C \left( 1 + a \left( \frac{E}{\text{keV}} \right)^c \mu + b \left( \frac{E}{\text{keV}} \right)^d \mu^2 \right), \quad (2)$$

$a, b, c$, and $d$ are beaming parameters. In the case of an energy-independent beaming correction, we fix $c = 0$, and $d = 0$, but in the case of an energy-dependent correction, we keep all of them free. The constant $C$ is determined so that $\int_0^1 \mu f(\mu, a, b, c, d) d\mu = 1/2$, which gives the following:

$$C = \frac{1}{1 + (2/3)a(E/\text{keV})^c + (1/2)b(E/\text{keV})^d} \quad (3)$$

This constant is chosen to minimize the change in the total flux spectrum $F(E) = \int_0^1 \mu I(E, \mu, a, b, c, d) d\mu$ for any given $a, b, c$, and $d$. This choice precisely conserves the emergent spectrum when the model intensity $I(E, \mu)_{NUM}$ does not depend on the emission angle. We found that this normalization factor already eliminates most of the difference in the spectrum caused by nonzero values of $a, b, c$, and $d$, leaving the spectrum unchanged with better than 1% accuracy for the beaming corrections we allow in our prior support. In addition, we performed a couple of low-resolution test runs either setting $C=1$ or finding the value of $C$ by numerical integration, requiring that $\int_0^1 \mu I(E, \mu)_{NUM} f(\mu, a, b, c, d) d\mu = \int_0^1 \mu I(E, \mu)_{NUM} d\mu$ separately for each emitted photon (increasing the computational

<footnote n="9">⁹ Note that the NICER analyses of Miller et al. (2019, 2021) do not use X-PSI but rather an independent ray-tracing and inference code.</footnote>
<table><caption>**Table 2**
Summary Table for PSR J0740+6620 Runs</caption><thead><tr><th>Model</th><th>Atmosphere</th><th>Free Beaming</th><th>Background</th><th>Settings</th><th>R<sub>eq</sub> (km)</th><th>M (M<sub>☉</sub>)</th><th>Evidence ln Z</th></tr></thead><tbody><tr><td>ST-U</td><td>nsxHe</td><td>No</td><td>SW</td><td>High</td><td>12.43<sup>+1.19</sup><sub>-1.00</sub> (13.23)</td><td>2.071<sup>+0.065</sup><sub>-0.065</sub> (1.949)</td><td>-16160.31 ± 0.07</td></tr><tr><td>ST-U</td><td>nsxHp</td><td>No</td><td>SW</td><td>High</td><td>12.40<sup>+1.14</sup><sub>-0.97</sub> (13.76)</td><td>2.074<sup>+0.065</sup><sub>-0.065</sub> (2.047)</td><td>-16160.83 ± 0.07</td></tr><tr><td>ST-U</td><td>nsxH</td><td>No</td><td>SW</td><td>High</td><td>12.30<sup>+1.16</sup><sub>-0.95</sub> (12.82)</td><td>2.072<sup>+0.064</sup><sub>-0.066</sub> (1.938)</td><td>-16162.06 ± 0.07</td></tr><tr><td>ST-U</td><td>hatmH</td><td>No</td><td>SW</td><td>High</td><td>12.41<sup>+1.21</sup><sub>-0.99</sub> (12.36)</td><td>2.069<sup>+0.065</sup><sub>-0.063</sub> (2.014)</td><td>-16162.07 ± 0.07</td></tr><tr><td>ST-U</td><td>nsxH</td><td>Max 10%</td><td>SW</td><td>High</td><td>12.18<sup>+1.02</sup><sub>-0.88</sub> (12.18)</td><td>2.073<sup>+0.064</sup><sub>-0.063</sub> (1.994)</td><td>-16162.27 ± 0.07</td></tr><tr><td>ST-U</td><td>nsxCp</td><td>No</td><td>SW</td><td>High</td><td>11.93<sup>+0.90</sup><sub>-0.72</sub> (12.42)</td><td>2.075<sup>+0.066</sup><sub>-0.067</sub> (2.145)</td><td>-16174.34 ± 0.08</td></tr></tbody></table>

**Note.** For explanations of the atmosphere tables, see Table 1. The “Free Beaming” column shows if a beaming correction is allowed, how large the correction can be, and whether it is energy-dependent (E-dep). The “Background” column shows whether background constraints have been applied using either XMM-Newton (XMM) or the NICER space weather estimate (SW) or if no constraints were applied (No). The “Settings” column shows whether high- or low-resolution settings are used for the X-PSI hot region cells, phase bins in the star frame, and energy bins for specific photon flux (see Vinciguerra et al. 2023a), and if the multimodal (MM) variant of MULTINEST was used. For all runs, we used $4 \times 10^3$ live points and 0.1 sampling efficiency. The $R_{\text{eq}}$ and $M$ columns show the radius and mass 68.3% credible intervals around the median values and the values corresponding to the maximum likelihood sample (in parentheses).

cost significantly). We found that the inferred radii do not vary by more than the 1 km fluctuation caused by the low-resolution sampling settings (see Section 2.4). This is probably a result of the inferred beaming parameters being close to zero, leading $C$ to be close to unity.

In the inference runs performed in this work, the prior support for the beaming parameters was based on a condition that the relative correction in the intensity (i.e., $|(f(E, \mu) - 1)|$) needs to remain below a predefined threshold (otherwise the sample is rejected). For PSR J0740+6620, we required the correction to be at most 10% at any energy and emission angle. For PSR J0030+0451, we required the correction to be at most 5% at any energy and emission angle below 60° (setting no limits for angles higher than that).$^{10,11}$

In both cases, the maximum allowed mean correction (averaged over all energies and angles) was around 6%–9%. These limits were chosen to avoid too much freedom in the parameter space, and they were roughly based on the differences between atmosphere models shown in Section 2.5.

In the case of PSR J0740+6620, the bounds of the beaming parameters $a$ and $b$ were set to $[-1.0, 1.0]$ while keeping $c$ and $d$ fixed to zero, and in the case of PSR J0030+0451, the bounds for all $a, b, c$, and $d$ were set to $[-0.7, 0.7]$. These ranges were found to always include at least 85% of the prior space of each parameter.

For simplicity, we also assumed the beaming parameters to be shared between all the emitting regions, even though the regions could have different temperatures, and the beaming correction might depend on the temperature.

## 2.2. Pulse Profile Modeling Using X-PSI

We calculate the energy-resolved X-ray pulses using the same pulse profile modeling technique, relying on the “Oblate Schwarzschild” approximation, as was done in the previous NICER analyses (see, e.g., Bogdanov et al. 2019b; and the references in Section 1). Our choices regarding the model parameters (besides the beaming parameters explained in Section 2.1.2), prior distributions, and likelihood calculation, follow those used in Vinciguerra et al. (2023a, 2023b) for PSR J0030+0451, and those used in Salmi et al. (2022) for PSR J0740+6620. For example, in both cases, we allow multiple imaging for the relativistic ray tracing (as explained in Riley et al. 2021). A few more features of the models are described next.

As in the previous works, for PSR J0740+6620, we use the ST-U (Single-Temperature-Unshared) hot region model, with two circular regions having uniform (but different) effective temperatures of the atmosphere. For PSR J0030+0451, the headline results from Riley et al. (2019) were produced using ST+PST (Single-Temperature+Protruding-Single-Temperature), in which one hot region is a circle, and the other is allowed to form ring and crescent-like structures (both with uniform but different temperatures). However, the ST+PST model is computationally expensive, so it is not feasible to fully explore the effects of atmosphere choices using ST+PST. Instead, we first performed an exhaustive exploration of the effects of different atmosphere models on PSR J0030+0451 using the much simpler ST-U spot model, which is the simplest tested model that does not show any structures in the residuals. This exploration allows us to find the cases where the ST-U results for PSR J0030+0451 are most sensitive to atmosphere choices. These most sensitive cases were selected for further analysis with the more complex ST+PST spot model.

When using ST+PST, we reduce the resolution in hot region cells, phase bins in the star frame, and energy bins for specific photon flux, as in the low-resolution runs in Vinciguerra et al. (2023a, see their Section 2.3.1 for details of the settings). This is done to mitigate the extra computational expense caused by the more complex model, while still allowing the use of high

<footnote>$^{10}$ These checks were done for 100 (or 50) emission angles linearly spaced in $\mu = [0.0, 1.0]$ (or $\mu = [0.5, 1.0]$) and for 100 energies linearly spaced in $\log_{10}(E/T_{\text{eff}}) = [-1.5, 2.2]$, using $T_{\text{eff}}$ of both the primary and secondary hot regions.</footnote>

<footnote>$^{11}$ These conditions produce an anticorrelation between $a$ and $b$, since a constant noncorrecting beaming function is found more easily if these factors have opposite signs so that the $\mu$ and $\mu^2$ terms cancel each other. For more efficient sampling, these priors can be mimicked by drawing $a$ and $b$ from normal distributions centered around 0 and $-a$, respectively. However, for simplicity, only initially uniform distributions were used for the final runs reported here.</footnote>
<table><caption>**Table 3**
Summary Table for PSR J0030+0451 Runs</caption><thead><tr><th>Model</th><th>Atmosphere</th><th>Free Beaming</th><th>Background</th><th>Settings</th><th>R<sub>eq</sub> (km)</th><th>M (M<sub>☉</sub>)</th><th>Evidence ln Ẑ</th></tr></thead><tbody><tr><td>ST+PST</td><td>nsxHe</td><td>No</td><td>No</td><td>Low-MM</td><td>15.39<sup>+0.43</sup><sub>-0.79</sub> (15.61)</td><td>1.83<sup>+0.13</sup><sub>-0.17</sub> (1.97)</td><td>-35778.17 ± 0.07</td></tr><tr><td>ST+PST</td><td>nsxH</td><td>No</td><td>No</td><td>Low-MM</td><td>13.11<sup>+1.30</sup><sub>-1.30</sub> (13.09)</td><td>1.37<sup>+0.17</sup><sub>-0.17</sub> (1.38)</td><td>-35778.42 ± 0.07</td></tr><tr><td>ST-U</td><td>nsxHe</td><td>E-dep, max 5%</td><td>No</td><td>High</td><td>15.44<sup>+0.41</sup><sub>-0.75</sub> (15.98)</td><td>1.78<sup>+0.16</sup><sub>-0.19</sub> (2.00)</td><td>-35783.17 ± 0.07</td></tr><tr><td>ST-U</td><td>nsxH</td><td>E-dep, max 5%</td><td>No</td><td>High</td><td>10.55<sup>+1.28</sup><sub>-0.96</sub> (10.02)</td><td>1.14<sup>+0.16</sup><sub>-0.10</sub> (1.18)</td><td>-35784.92 ± 0.07</td></tr><tr><td>ST-U</td><td>nsxHe</td><td>No</td><td>No</td><td>High</td><td>14.00<sup>+1.39</sup><sub>-1.88</sub> (14.89)</td><td>1.45<sup>+0.21</sup><sub>-0.22</sub> (1.70)</td><td>-35785.22 ± 0.07</td></tr><tr><td>ST-U</td><td>nsxH</td><td>No</td><td>No</td><td>High</td><td>10.59<sup>+1.29</sup><sub>-0.92</sub> (9.99)</td><td>1.12<sup>+0.15</sup><sub>-0.09</sub> (1.05)</td><td>-35786.21 ± 0.07</td></tr><tr><td>ST-U</td><td>hatmH</td><td>No</td><td>No</td><td>High</td><td>10.68<sup>+1.32</sup><sub>-1.00</sub> (11.08)</td><td>1.14<sup>+0.15</sup><sub>-0.09</sub> (1.18)</td><td>-35787.01 ± 0.07</td></tr><tr><td>ST-U</td><td>nsxHp</td><td>No</td><td>No</td><td>High</td><td>12.94<sup>+0.98</sup><sub>-0.81</sub> (12.19)</td><td>1.40<sup>+0.15</sup><sub>-0.14</sub> (1.40)</td><td>-35796.14 ± 0.07</td></tr><tr><td>ST+PST<sup>a</sup></td><td>nsxHp</td><td>No</td><td>No</td><td>Low-MM</td><td>11.18<sup>+0.88</sup><sub>-0.80</sub> (12.15)</td><td>1.22<sup>+0.14</sup><sub>-0.12</sub> (1.17)</td><td>-35821.04 ± 0.07</td></tr><tr><td>ST-U</td><td>nsxCp</td><td>No</td><td>No</td><td>High</td><td>10.95<sup>+0.77</sup><sub>-0.65</sub> (11.24)</td><td>1.03<sup>+0.04</sup><sub>-0.02</sub> (1.00)</td><td>-35963.05 ± 0.07</td></tr><tr><td>ST-U</td><td>nsxHe</td><td>No</td><td>XMM</td><td>High-MM</td><td>15.22<sup>+0.57</sup><sub>-1.10</sub> (15.82)</td><td>1.83<sup>+0.23</sup><sub>-0.26</sub> (2.11)</td><td>-42694.06 ± 0.07</td></tr><tr><td>ST-U</td><td>nsxH</td><td>No</td><td>XMM</td><td>High-MM</td><td>15.11<sup>+0.65</sup><sub>-1.32</sub> (15.19)</td><td>1.88<sup>+0.13</sup><sub>-0.20</sub> (1.81)</td><td>--42713.88 ± 0.07</td></tr><tr><td>ST-U</td><td>nsxHp</td><td>No</td><td>XMM</td><td>High-MM</td><td>11.02<sup>+1.21</sup><sub>-1.12</sub> (10.97)</td><td>1.27<sup>+0.16</sup><sub>-0.14</sub> (1.16)</td><td>-42728.79 ± 0.07</td></tr></tbody></table>

**Notes.** See the explanation for columns and their contents in Table 2. Unlike the previous table, here, all the runs were performed with $1 \times 10^4$ live points and 0.3 sampling efficiency. Note that the evidence for models with different background constraints is not comparable.

<sup>a</sup> In addition to the worse evidence, the maximum likelihood found for this model is significantly worse than the likelihood obtained by plugging the best-fit ST-U parameter vector (from the run with nsxHp atmosphere) into the ST+PST framework (−35769.98 instead of −35737.93 in ln units). This means that the sampler has not fully explored the parameter space.

enough sampling resolution of the parameter space (see Section 2.4).

For the instrument response model, we use only energy-independent effective area scaling factors for both PSR J0740+6620 and PSR J0030+0451. We combined the uncertainty in the distance and the scaling factor into a parameter called $\beta$ for the latter (as in Vinciguerra et al. 2023a), except when jointly fitting NICER and XMM-Newton data. Energy-dependent scaling factors were originally used in Riley et al. (2019) but were later replaced with energy-independent ones for the newer PSR J0030+0451 and PSR J0740+6620 analysis as explained in Vinciguerra et al. (2023a).

The calculation of the background-marginalized likelihood function is similar to that used in the previous works. For PSR J0030+0451, we apply background constraints using XMM-Newton (similarly to Vinciguerra et al. 2023b) for atmosphere cases that showed large differences in the radius when XMM-Newton was not included. For PSR J0030+0451, the background is typically inferred to include most of the unpulsed emission (Miller et al. 2019; Riley et al. 2019) if one does not use any constraints; with XMM-Newton, we can limit the NICER background to be smaller (by constraining the number of source photons). For all of our PSR J0740+6620 analyses, we apply the nonsmoothed NICER space weather background estimate as a lower limit for the background without using the XMM-Newton data (corresponding to the case labeled W21-0.9xSW in Salmi et al. 2022). This was shown to give similar results as constraining the background using XMM-Newton for PSR J0740+6620 in Salmi et al. (2022).

## 2.3. Data Sets

In this section, we summarize the event data sets used in the analysis of this paper.

### 2.3.1. PSR J0740+6620

For PSR J0740+6620 analysis, we used the same pre-processed NICER X-ray event data (described in Wolff et al. 2021) and instrument response files as in Miller et al. (2021), Riley et al. (2021), and we partly used the same data as in Salmi et al. (2022). Instead of including XMM-Newton data in the analysis, we used a lower limit for the NICER background based on the space weather background estimator (Gendreau 2020), as in some models of Salmi et al. (2022).

### 2.3.2. PSR J0030+0451

For PSR J0030+0451 analysis, we used the same pre-processed X-ray event data and instrument response files as in Vinciguerra et al. (2023b). The NICER data set is thus similar to that introduced in Bogdanov et al. (2019a), and used in Riley et al. (2019), Miller et al. (2019), but updated adopting the most recent NICER instrument response, as in Riley et al. (2021), Salmi et al. (2022). In addition, for some analyses, we also used the same XMM-Newton data as introduced in Bogdanov & Grindlay (2009), Vinciguerra et al. (2023b).

## 2.4. Posterior Computation

As in the previous works, we compute the posterior samples using MULTINEST (Feroz & Hobson 2008; Feroz et al. 2009; Buchner et al. 2014). For details of our nested sampling
protocol, we refer to Riley et al. (2019, especially the appendices), Riley (2019, especially chapter 3 and the associated appendix). The resolution settings for sampling were chosen to be adequate for an exploratory analysis, where the main focus is to compare the effects of different atmosphere model choices.

For the PSR J0740+6620 analysis, we employed the “low-resolution” settings from Riley et al. (2021): $4 \times 10^3$ live points; and 0.1 sampling efficiency<fnref n="12" />; and an estimated remaining log-evidence of $10^{-7}$. This was found to produce stable results that do not significantly change from run to run with different random seeds for sampling with MULTINEST. However, as noted in Riley et al. (2021), Salmi et al. (2022), using 10 times more live points would typically slightly broaden the credible intervals and shift the median radius up by 0.2–0.4 km.

For the PSR J0030+0451 analyses, we used slightly different resolution settings; $1 \times 10^4$ live points; and 0.3 sampling efficiency. These settings were chosen to make the computation reasonably fast (noting that PSR J0030+0451 analysis is more expensive than PSR J0740+6620 analysis with the same settings, likely due to the higher number of detected counts), while avoiding fluctuation in the results for runs with different MULTINEST seeds. The latter was found to be a problem, especially for our initial runs with additional beaming parameters, where $4 \times 10^3$ live points and 0.8 sampling efficiency (which seemed sufficient without the beaming parameters) resulted in the inferred NS radius varying by around 1 km between identical runs and always being biased toward higher radii compared to the run with higher sampling resolution. Therefore, the higher resolution was used in all the results reported in this paper. In most of our runs, the mode-separation variant of MULTINEST (known also as the multimodal option) was deactivated (i.e., isolated modes were not evolved independently), except for the joint NICER and XMM-Newton ST-U and NICER-only ST+PST fits of PSR J0030+0451. The settings used for each run are listed in Tables 2 and 3.

## 2.5. Model Comparison

We start by comparing emergent spectra and beaming patterns predicted by the different atmosphere models used in this paper. For the comparison, we select the effective temperature (of the hottest region) and surface gravity values that correspond to the maximum likelihood parameter vector inferred when applying the fully ionized hydrogen NSX ST-U model (model nsxH from Table 1) for either PSR J0030+0451 or PSR J0740+6620; these are $\log_{10}(T_{\text{eff}}/\text{K}) = 6.1059$, $\log_{10}(g_s/\text{cm s}^{-2}) = 14.143$, and $\log_{10}(T_{\text{eff}}/\text{K}) = 6.0409$, $\log_{10}(g_s/\text{cm s}^{-2}) = 14.194$ respectively. In addition, we also show one higher temperature example with $\log_{10}(T_{\text{eff}}/\text{K}) = 6.4$, and $\log_{10}(g_s/\text{cm s}^{-2}) = 14.15$. The different spectra are presented in Figure 1 for $\mu = 0.5$, i.e., for a 60° emission angle. This figure shows only small deviations between the different models when considering the relatively low temperatures typical for PSR J0030+0451 or PSR J0740+6620 (the difference in temperatures between the two is small). Only the partially ionized carbon model predicts a significantly different spectrum with several lines and edges (as in, e.g., Ho & Heinke 2009; Suleimanov et al. 2014). The impact of external heating starts to become major (larger than that between

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Upper panel: the emergent spectra for different atmosphere models at 60° emission angle. Temperature and surface gravity correspond to the best-fit results (for the hottest region) from PSR J0030+0451 analysis with the ST-U-nsxH model ($\log_{10}(T_{\text{eff}}/\text{K}) = 6.1059$, $\log_{10}(g_s/\text{cm s}^{-2}) = 14.143$). The different atmosphere models are explained in Table 1. Lower panel: the same as the upper panel, except with a higher temperature $\log_{10}(T_{\text{eff}}/\text{K}) = 6.4$, and surface gravity fixed to $\log_{10}(g_s/\text{cm s}^{-2}) = 14.15$. The spectra for best-fit parameters from PSR J0740+6620 analysis are available.</figcaption></figure>

(The complete figure set of 3 images is available.)

hydrogen and helium composition) only at higher temperatures, such as in the $\log_{10}(T_{\text{eff}}/\text{K}) = 6.4$ case shown in the lower panel of Figure 1. There, the partially ionized hydrogen spectrum also starts to deviate more from fully ionized hydrogen. Since features in opacity, such as resonances due to bound states, can cause absorption of photons and heat the outer atmospheric layers in partially ionized models, it is not too surprising that both partial ionization and external heating seem to have a similar effect on the spectrum, increasing the low-energy emission.

The different beaming patterns are presented in Figure 2. The beaming patterns for all of the atmosphere models except

<footnote n="12">¹² We note that this factor is not the number set as the MULTINEST sampling efficiency parameter, since it is first corrected by X-PSI to account for the nonuniform prior hypercube volume as explained in Riley (2019).</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Upper panel: the beaming patterns at 1.0 keV with temperature and surface gravity corresponding to the best-fit results from PSR J0030+0451 analysis with the ST-U-nsxH model. See Figure 1 and Table 1 for more details. Lower panel: the same as the upper panel, except with a higher temperature $log_{10}(T_{eff}/K) = 6.4$, and surface gravity fixed to $log_{10}(g_s/cm s^{-2}) = 14.15$. The beaming pattern at 0.5 keV and 2.0 keV with PSR J0030+0451 parameters and the pattern at 1.0 keV with PSR J0740+6620 parameters are available.

(The complete figure set of 5 images is available.)</figcaption></figure>

for the carbon model are very similar. We also see that the difference caused by composition or partial ionization is larger than the difference caused by a stopping layer for the heating particles, for the lower temperatures inferred for PSR J0030+0451 and PSR J0740+6620 (see the upper panel). As already presented in Bogdanov et al. (2021, see their Figure 1), the difference in the beaming pattern between hydrogen and helium (or between fully and partially ionized) atmosphere is typically smaller than 5% at 1 keV or smaller photon energies. However, the relative difference starts to increase when the

emission angle is larger than about 60°-70°. At 1 keV, both partially ionized hydrogen and fully ionized helium show a similar monotonic decrease in emission at higher angles; fully ionized hydrogen, however, decreases more slowly and turns to limb brightening at the highest angles. We note that this difference is expected to be important only if both hot regions are observed at high emission angles. Otherwise, the radiation is dominated by the small angle emission, which has higher intensity. The effect could be larger in the case of PSR J0030+0451 analysis without background estimates since, there, - only photons emitted with high emission angles can, typically, reach the observer, due to the inferred hot spot locations. For example, emission angles are always above 57° in the case of the maximum likelihood solution found for ST-U with a fully ionized hydrogen atmosphere.

### 3. Results

In this section, we show the inference results for both PSR J0740+6620 and PSR J0030+0451, when using different atmospheric choices. The main results, applying all the atmosphere tables, are shown in Figure 3 for PSR J0740+6620 and in Figure 4 for PSR J0030+0451. The cases impacting the results the most (for PSR J0030+0451) are shown in Figures 5–7. A summary of all the PSR J0740+6620 runs is presented in Table 2, and a summary of all the PSR J0030+0451 runs is presented in Table 3. In both cases, the rows are ordered based on decreasing evidence.

#### 3.1. Mass and Radius Constraints for PSR J0740+6620 with Different Atmospheres

We start by presenting the inference results for PSR J0740+6620 with different numerical atmospheres and using the ST-U hot region model in Figure 3 (see also the full list of inferred radius, mass, and evidence values in Table 2). As mentioned in Section 2.2, we also use the space weather background estimate as a lower limit for the NICER background. We find no significant differences in the inferred NS radii between any of the atmosphere models. The only model for which the posterior contours can be clearly distinguished from the rest is the partially ionized carbon atmosphere (especially in the inferred temperatures). However, as in the case of PSR J0030+0451 (see Section 3.2.1), the carbon model is significantly disfavored based on the evidence values (more than ~10 difference in natural logarithmic, i.e., ln units). The evidence differences between the other atmosphere cases are only marginal. The results also show that allowing 10% maximum deviation in the beaming pattern using the energy-independent beaming parameters (as explained in Section 2.1.2) for fully ionized hydrogen does not have any significant impact on the results.

#### 3.2. Mass and Radius Constraints for PSR J0030+0451 with Different Atmospheres

##### 3.2.1. ST-U NICER-only Fit

Next, we present the inference results for PSR J0030+0451 with different numerical atmospheres and using the ST-U hot region model in Figure 4 (see also Table 3). We see that the radius inferred using the externally heated, fully ionized hydrogen atmosphere model is consistent with the one from the corresponding deep-heated NSX model. There is only a
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Effect of the atmosphere model on the posterior distributions of radius, compactness, mass, and temperatures for both primary and secondary hot regions, using the PSR J0740+6620 NICER data set conditional on the ST-U model. Six types of posterior distribution are shown with five different atmosphere tables used in the analysis (see Table 1 for definitions), and one fully ionized hydrogen NSX model with a free beaming correction limited to 10% deviation at most (ST-U-nsxH-beam10). The credible intervals and the Kullback-Leibler (KL) divergence $D_{\text{KL}}$ are reported for the fully ionized hydrogen NSX case: ST-U-nsxH (see Table 2 for mass and radius values in other cases). The prior distributions are shown by the dashed-dotted functions. The 1D intervals (shaded in red) contain 68.3% of the posterior mass, and the contours in the 2D panels contain 68.3%, 95.4%, and 99.7% of the posterior mass. For more details of the figure elements, see Figure 5 of Riley et al. (2021) and Figure 5 of Salmi et al. (2022). Posterior distributions for the other parameters are shown in Figure 8 of the Appendix.</figcaption></figure>

small shift (smaller than ~0.03 in $\log_{10} T$ [K]) in the inferred temperature for both hot regions to higher values when using the externally heated model. This is expected because the spectrum peaks at slightly lower energy than the corresponding deep-heated model with the same model parameters (see Figure 1). However, at least in this temperature range and with the assumed return-current properties mentioned in Section 2.1.1, this effect is too small to cause any significant change in the inferred radius.

On the contrary, both the partial ionization and the chemical composition have significant effects on the inferred radius. The use of a partially ionized hydrogen model instead of a fully ionized atmosphere shifts the median radius and its 68.3% credible interval from $10.59_{-0.92}^{+1.29}$ to $12.94_{-0.81}^{+0.98}$ km, i.e., leaving no overlap in the 68.3% intervals. The inferred median mass also increases from around 1.1 to 1.4 $M_\odot$, and the temperatures of both hot regions become higher (again as expected based on the spectral shapes in Figure 1). As seen
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Effect of the atmosphere model on the posterior distributions of radius, compactness, mass, and temperatures for both primary and secondary hot regions, using the PSR J0030+0451 NICER data set conditional on the ST-U model. Five types of posterior distribution are shown with different atmosphere tables used in the analysis (see Table 1 for definitions). The 1D credible intervals and the KL-divergence estimates are reported for the fully ionized hydrogen NSX case: ST-U-nsxH (see Table 3 for mass and radius values in other cases). See the caption of Figure 3 for additional details about the figure elements. Posterior distributions for the other parameters are shown in Figure 9 of the Appendix.</figcaption></figure>

from Figure 9 of the Appendix, there are also small changes in the other model parameters. However, both hot regions are still inferred to be on the same hemisphere and highly inclined toward the observer. The evidence for the model with partially ionized hydrogen is worse than that for the model with fully ionized hydrogen, by 10 ln units (see Table 3).

Changes in the chemical composition of the atmosphere also cause the inferred radius to change, in the case of PSR J0030+0451. For partially ionized carbon, the inferred radius 68.3% credible interval is 10.95<sup>+0.77</sup><sub>-0.65</sub> km, which does not overlap with the corresponding interval of partially ionized hydrogen. In addition, major differences are seen in many of the other parameters. However, based on the evidence, the differences
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Effect of the atmosphere model on the posterior distributions of radius, mass, and beaming parameters *a*, *b*, *c*, and *d*, using the PSR J0030+0451 NICER data set conditional on the ST-U model and energy-dependent beaming parameterization. Priors (dashed-dotted curves) and posteriors for the four beaming parameters are shown only for the models where those parameters were kept free. The same two hydrogen and helium posteriors without beaming modifications are also shown, as in Figure 4 (ST-U-nsxH and ST-U-nsxHe), in the radius and mass panels. The posteriors for ST-U-nsxH and ST-U-nsxH-E-beam5 models are almost exactly overlapping. The credible intervals and KL-divergence estimates are reported for the latter (see Table 3 for mass and radius values in other cases). See the caption of Figure 3 for additional details about the figure elements. Posterior distributions for the other parameters are shown in Figure 10 of the Appendix.</figcaption></figure>

compared to the other models, which are more than ~200 in ln units (see Table 3), and the quality of the residuals between the data and best-fit model, the carbon model is clearly disfavored for PSR J0030+0451. For fully ionized helium, the evidence is similar (or even slightly higher) to that of hydrogen, and the effect on the shift in the inferred radius is the largest; from

10.59<sup>+1.29</sup><sub>-0.92</sub> (for hydrogen) to 14.00<sup>+1.39</sup><sub>-1.88</sub> km (for helium). The inferred mass also increases to around 1.5 M<sub>☉</sub>. Most of the other model parameters are not greatly affected, although there is more posterior mass for the hot regions and observer inclination to be closer to the equatorial plane compared to the hydrogen case.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Effect of the atmosphere model on the posterior distributions of radius, compactness, mass, radius, and temperatures for both primary and secondary hot regions, using the PSR J0030+0451 NICER data set conditional on the ST+PST model. Posterior distributions are shown for fully ionized hydrogen and helium, and for partially ionized hydrogen NSX models. The credible intervals and KL-divergence estimates are reported for the ST+PST-nsxH model (see Table 3 for mass and radius values of the other models). See the caption of Figure 3 for additional details about the figure elements. Posterior distributions for the other parameters are shown in Figure 11 of the Appendix.</figcaption></figure>

### 3.2.2. ST-U NICER-only Fit with Free Beaming

As explained in Section 2.1.2, we have also performed a few runs with additional beaming parameters introduced to add more freedom in the predicted beaming pattern. The main results of these runs, for fully ionized hydrogen and helium, are shown in Figure 5 and in Table 3. We see that limiting the beaming correction to 5% at maximum (at any energy or

emission angles below 60°) has no significant effects on the radius and mass constraints in the hydrogen case. For helium, the median radius shifts from 14.0 to 15.4 km, but the 68.3% credible intervals are still overlapping. For hydrogen, no significant information gain from prior to posterior distribution was found in any of the beaming parameters. The inferred values thus peak around zero in all the correction factors *a*, *b*, *c*,
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Effect of the atmosphere model on the posterior distributions of radius, compactness, mass, and temperatures for both primary and secondary hot regions, using the PSR J0030+0451 joint NICER and XMM-Newton (NxX) data sets conditional on the ST-U model. Posterior distributions are shown for fully ionized hydrogen and helium, and partially ionized hydrogen NSX models. The credible intervals and KL-divergence estimates are reported for the ST-U-NxX-nsxH model (see Table 3 for mass and radius values of the other models). See the caption of Figure 3 for additional details about the figure elements. Posterior distributions for the other parameters are shown in Figure 12 of the Appendix.</figcaption></figure>

and *d*, following closely the prior distribution. Also, the uncertainty in these parameters did not cause any noticeable broadening of the credible intervals of radius or mass. For helium, however, the inferred *a* and *b* differ significantly from zero (hitting the upper bound of 0.7 for *a*), and the radius and mass constraints become tighter most likely due to the cut in radius at 16 km.

### 3.2.3. ST+PST NICER-only Fit

We also performed a few runs with the more complex and computationally expensive ST+PST model for the hot region shapes, as mentioned in Section 2.2. This model corresponds to that used for the headline results of Riley et al. (2019). Since this model is computationally more expensive, runs were only done for the fully ionized hydrogen and helium, and partially
ionized hydrogen cases, in which the difference in radius was largest. The results are shown in Figure 6 and in Table 3. We see that both the inferred radius and mass increase when using a helium instead of a hydrogen atmosphere (both fully ionized). The radius changes from $13.11_{-1.30}^{+1.30}$ to $15.39_{-0.79}^{+0.43}$ km (the latter being limited by the 16 km prior upper limit), and mass from $1.37_{-0.17}^{+0.17}$ to $1.83_{-0.17}^{+0.13} M_{\odot}$ (which is still smaller than the 2.7 $M_{\odot}$ value found in Miller et al. 2019 when using a helium model; see discussion in Section 4.2). Again, there is no substantial difference in the evidence between helium and hydrogen.

In the case of partially ionized hydrogen, the inferred radius is $11.18_{-0.80}^{+0.88}$ km, which barely overlaps with the 68.3% interval of fully ionized hydrogen. Significant differences are found in many of the other parameters, for example, inferring more antipodal-like geometry with hot regions and observer inclination close to the equatorial plane (as seen from Figure 11 of the Appendix). However, the maximum likelihood found is more than 30 $\ln$ units worse than the likelihood for the corresponding best-fit ST-U solution in the ST+PST framework. This means that even the high-resolution settings ($1 \times 10^4$ live points; and 0.3 sampling efficiency) were not sufficient to fully explore the parameter space.

As mentioned in Section 2.4, the mode-separation variant was used for these runs to better identify and explore multiple modes found by the sampler. Using that, a secondary mode with worse evidence was found for all the models. The second mode for fully ionized hydrogen corresponds to a radius of around 11 km (instead of 13 km of the main mode), while the second modes of fully ionized helium and partially ionized hydrogen have similar radii to the corresponding main modes (i.e., 15 and 11 km, respectively).

### 3.2.4. ST-U NICER and XMM-Newton Fit

As in Vinciguerra et al. (2023b), we also performed a few runs where the NICER background signal is constrained by fitting jointly NICER and XMM-Newton data. The results are shown in Figure 7 and in Table 3. We see that the inferred radius is equally high for both the fully ionized hydrogen and helium models (unlike in the previous cases), i.e., close to the prior upper limit of 16 km. The mass is around $1.9 M_{\odot}$, and the hot region geometry is similar (more antipodal compared to the corresponding runs without XMM-Newton) in both the hydrogen and helium cases. However, a few other parameters, including the temperatures and interstellar absorption parameter $N_H$, are inferred to have notably different credible intervals. Interestingly, the evidence for the helium case is this time significantly higher than that for the hydrogen case (around 20 in $\ln$ units).

Hot regions close to the equatorial plane are also inferred using the partially ionized hydrogen atmosphere. In this case, however, the inferred radius is significantly smaller, around 11 km. Nevertheless, the evidence for the partially ionized hydrogen model is around 15 $\ln$ units worse than that for the fully ionized hydrogen.

The mode separation (see Section 2.4) allowed us to identify, in the case of the fully ionized hydrogen atmosphere, two distinct modes (similarly to Vinciguerra et al. 2023b). The secondary mode corresponds to a smaller radius around 11 km and a slightly different geometry, but with a significantly smaller local evidence. In the case of the fully ionized helium atmosphere, we found four distinct modes, which all correspond to radii above 15 km with slightly different

geometries and local evidence better than in the best mode found for the hydrogen run. In the case of the partially ionized hydrogen atmosphere, we found two modes, where the secondary mode corresponds to a radius of around 12 km instead of the 11 km of the primary mode.

## 4. Discussion

As shown in Section 3, we found that the radius inferred for PSR J0030+0451 is sensitive to some of the assumptions made for the atmosphere model. On the other hand, the results for PSR J0740+6620 were shown to be largely insensitive to any of the model assumptions that we considered here. In the following, we discuss possible reasons for this difference and what can be learned about the atmosphere properties of these sources based on our study.

### 4.1. Differences between J0030 and J0740

There may be several reasons why the PSR J0740+6620 analysis appears to be less affected by the atmosphere model choices than the PSR J0030+0451 analysis. One of them is the quality of the data; due to the larger number of total counts in the PSR J0030+0451 data set ($1.5 \times 10^6$ versus $0.6 \times 10^6$) and the higher background for PSR J0740+6620, the credible intervals are notably tighter for the former and more sensitive to systematic errors caused by unaccounted-for physics.

However, the sensitivity to atmosphere could also depend on the properties of the star itself. The radio observations of PSR J0740+6620 constrain the observer's inclination to be very close to the star's orbital plane (and thus the star's equatorial plane). The highest likelihood solutions have a hot spot that also tends to lie close to the equatorial plane. As a result of the constrained geometry, the hot spot on the equator is brightest when facing the observer, which is also when the light is emitted normal to the surface, minimizing the effect of beaming by the different atmosphere models. The inclination of PSR J0030+0451 is unconstrained by radio observations, and the highest likelihood solutions tend to favor, in most cases, the models where the spots are close to the limb of the star, as viewed by the observer. As a result, the emitted light that reaches the observer tends to be emitted close to the tangent to the surface, where differences in beaming are much more pronounced. This can result in a stronger dependence on the atmosphere model for PSR J0030+0451.

As presented in Section 3.2.4, the atmosphere effects, for hydrogen versus helium, on the inferred radii for PSR J0030+0451 are still small when we consider the joint NICER and XMM-Newton runs. The similarity could be caused either by the hard cut in the prior upper limit at 16 km or by the fact that the inferred results for both hydrogen (see also Vinciguerra et al. 2023b) and helium predict significantly different hot region geometries than those without background constraints. In this case, most of the detected photons are emitted at smaller emission angles (down to around 5° instead of 60°) where the relative difference between the atmosphere models is smaller. However, the inferred radius for the corresponding partially ionized hydrogen case is significantly smaller (around 11 km), even with a similar best-fit geometry where the detected photons are dominated by the small angle emission.

Besides the data quality and viewing geometry, it could also be that external parameter constraints help to reduce the free parameter space such that different results with different
atmosphere models are not found. Thus, the tight priors on mass, distance, and inclination for PSR J0740+6620 could also explain its insensitivity to the atmosphere model choices.

## 4.2. Implications for Chemical Composition

In the case of both PSR J0030+0451 and PSR J0740+6620, both hydrogen and helium compositions of the atmosphere were found to mostly fit the data equally well. However, the inferred radii of PSR J0030+0451 using helium tend to hit the 16 km upper limit of the prior, making these solutions less plausible. Differences in the evidence are only a few units in natural logarithm (with helium always being better), except in the case of the joint NICER and XMM-Newton analysis on PSR J0030+0451, where the evidence for helium is around 20 units larger than that for hydrogen. However, it could still be possible that small differences in the constrained background might lead to substantially different NS parameters and preferred atmosphere. It would therefore be beneficial to perform ancillary analysis using direct background estimates for NICER as in Salmi et al. (2022), which would give information about the background independently from XMM-Newton.

When comparing the evidence between hydrogen and helium, one could also weigh the evidence based on how much more likely the hydrogen composition is considered a priori. However, quantifying this is difficult. It depends on the probability that the NS has accreted from a hydrogen-poor companion and/or that the diffusive nuclear burning has converted all hydrogen to helium. In addition, if even a small amount of hydrogen (~$10^{-19}$ $M_\odot$, Romani 1987) has been accreted after the mass transfer from the companion ended, it would be enough to create an upper atmosphere of pure hydrogen (Blaes et al. 1992; Wijngaarden et al. 2019, 2020; Bogdanov et al. 2021).

As noted in Section 3.2.3, our helium results for PSR J0030+0451 differ from those in Miller et al. (2019); we infer a significantly lower NS mass, $1.8 M_\odot$ (for ST+PST), rather than $2.7 M_\odot$. This can likely be explained by the different prior assumptions, especially the more strict 16 km upper limit for the radius in the analysis of this paper, since this also limits how high a mass can be obtained for a given NS compactness (see Figure 6). In addition, Miller et al. (2019) use a completely independent analysis pipeline with the differences discussed, e.g., in Miller et al. (2021), Riley et al. (2021).

As mentioned in Section 1, a composition with elements heavier than helium might also be possible, although less likely. Given that we found significantly worse evidence for the tested carbon atmosphere cases (more than ~200 for PSR J0030+0451 and more than ~10 for PSR J0740+6620 in $\ln$ units), it seems clear that carbon atmospheres are very unlikely for these sources. Having any other heavy element composition or a mix of different elements was not considered in this study, as they are both considered unlikely. In particular, the latter seems improbable due to the rapid sinking of elements via diffusive gravitational separation as explained in Section 1.

## 4.3. Implications for Partial Ionization

Some of the atmosphere models we applied allowed the possibility of the atoms in the atmosphere to be partially ionized. The effect of the ionization state is usually expected to be highest for the colder atmospheres where the ionization

fraction is noticeably smaller than unity. However, since the atmosphere consists of layers with different temperatures, even an NS with a relatively high effective temperature might have shallow layers with low enough temperatures to produce bound-bound and bound-free opacity features in the escaping radiation. This could explain why the inferred radius results for PSR J0030+0451 (see Section 3.2.1) appear to be sensitive to the choice of ionization state, even if the inferred effective temperature is around $log_{10}(T_{\text{eff}}/\text{K}) \approx 6.1$ for both hot regions.

Nonetheless, at higher temperatures, the partially ionized atmosphere calculation becomes more inaccurate due to the opacity uncertainties for a range of photon energies (Iglesias & Rogers 1996; Badnell et al. 2005; Colgan et al. 2016; Bogdanov et al. 2021). Thus, interpreting the partially ionized results demands some extra caution. For $log_{10}(T_{\text{eff}}/\text{K}) \approx 6.1$, the opacity uncertainties should still be small, but larger at $log_{10}(T_{\text{eff}}/\text{K}) \approx 6.4$, which is inferred for both hot regions of PSR J0030+0451 with the partially ionized carbon case. Despite these uncertainties, we count the much worse evidence of the carbon model as a strong indication that the atmosphere of PSR J0030+0451 is not composed of carbon.

In the case of PSR J0740+6620, we found a good agreement between the fully and partially ionized hydrogen results, which is consistent with that from Miller et al. (2021), Riley et al. (2021), Salmi et al. (2022), where the code applying the partially ionized table (Miller et al. 2021) obtained consistent results with those using the fully ionized table (Riley et al. 2021). In our results, the evidence difference between the two models is not significant (partially ionized model performs better by less than 2 in $\ln$ units). In the case of PSR J0030+0451, however, we always find significantly worse evidence when using the partially ionized hydrogen model (at least 10 $\ln$ units in all cases). This suggests that the fully ionized models should be preferred over the partially ionized models for PSR J0030+0451; however, the latter predicts more reasonable NS radii when including XMM-Newton data (11 instead of 15 km) for the ST-U spot configuration.

## 4.4. Other Atmospheric Effects

In addition to the chemical composition and ionization state of the atmospheric elements, some other atmosphere model effects on the NS radius constraints were also considered in Section 3. One of them is the stopping depth of the return-current particles in the atmosphere. We studied this by using an externally heated, fully ionized hydrogen atmosphere model with certain assumptions in the properties of the return-current particle energy distribution (see Section 2.1.1). In this case, no significant changes in the inferred radius were observed for either PSR J0030+0451 or PSR J0740+6620 when compared to the corresponding deep-heated model. As mentioned in Section 3.2.1, this shows that even the slow bombarding particles (with average Lorentz factors of 100) do not stop at depths shallow enough to affect NS radius constraints. In principle, the effects could be larger for NSs with hotter spots (e.g., around $log_{10} T_{\text{eff}} = 6.4$), where the impact on the emergent intensity is larger, as seen in Section 2.5. However, even then, the effect is likely to be unimportant, since the return current particles are more likely notably faster than those assumed here, based on pulsar magnetosphere simulations (Harding & Muslimov 2011; Timokhin & Arons 2013; Chen & Beloborodov 2014; Cerutti & Beloborodov 2017; Philippov & Spitkovsky 2018). Most of them have focused on slowly
rotating and highly magnetic pulsars, but the particle energies
are expected to be larger for the millisecond pulsars with lower
magnetic fields (see Figure 10 of Harding & Muslimov 2011;
and the discussion in Bogdanov et al. 2021).

We also verified that the numerical implementation of the
atmosphere calculation is unlikely to have large effects on the
NS parameter constraints since the externally heated atmos-
phere table was produced with a completely different
algorithm. It seems that the smaller grid resolution of that
atmosphere table (see Table 1 for a summary) does not affect
the results either. However, there are still effects that we have
not considered in this study. One of them is the effect of the
magnetic field strength on the atmosphere. As discussed in
Bogdanov et al. (2021), such effects on the spectrum and
beaming pattern should be negligible if PSR J0030+0451 has a
purely dipolar magnetic field. In the presence of multipoles, the
field strength could get higher, but the NS atmosphere would
still likely be not significantly affected by the magnetic field
(Kalapotharakos et al. 2021). Another possibility with a high
magnetic field (unlikely for RMPs) would be the formation of a
bare condensed NS surface, instead of an atmosphere (Medin &
Lai 2007).

4.5. Beaming Uncertainty

Ideally, the uncertainty in the atmosphere model assump-
tions and numerical calculations could be parameterized (as
attempted in Section 2.1.2) instead of performing multiple runs
with different assumptions. Our empirical beaming correction
function allowed reasonable freedom in the predicted angular
department of radiation. However, the NS radius results were
shown to be largely not sensitive to the allowed beaming
correction in all the cases we considered (see Sections 3.2.2 and
3.1), which is in contrast to the much higher sensitivity found
when applying different atmosphere tables for PSR J0030
+0451. In the case of hydrogen, the values for all the beaming
correction factors a, b, c, and d were found to follow closely
their prior distributions, centered around zero. In the case of
helium, significantly different values were inferred for a and b,
but the nominal beaming was still changed only by 6% on
average at 1 keV, while the corresponding average difference
between the nsxH and nsxHe tables is about 2 times larger.

There are many reasons that could explain the smaller
sensitivity to the beaming correction compared to the
atmosphere table choice. One is that allowing a change in
only the beaming pattern may not be enough, and we should
also account for the change in the emergent spectrum, if we
want to allow similar parameter space to be explored when
using either a beaming-modified hydrogen or helium atmos-
phere. Another reason could be that we allowed no more than a
5% deviation in beaming in our PSR J0030+0451 runs at
emission angles below 60°, which can indirectly limit the
correction also at the highest angles (due to the functional form
of the beaming correction) to be less than the differences
between the atmosphere tables (which rapidly increase at the
highest angles for some cases).

An additional reason could be that the form of the correction
function is not general enough, to accurately capture the
differences between the different atmosphere models even at
smaller emission angles.

A more detailed analysis of this is left for future work.

4.6. Comparison to Other Modeling Uncertainties

Based on our results, and discussion in Section 4.1, it seems
clear that the sensitivity of the mass and radius constraints to
NS atmospheric effects can vary a lot from star to star. Due to
the complicated degeneracies between atmospheric beaming,
hot spot geometries, and background radiation, it is difficult to
determine which of these factors is most important. Different
assumptions in the atmosphere can lead to different modes
being preferred in the multimodal posterior surface, especially
when the solutions (e.g., in terms of geometry parameters)
favored in some atmosphere cases are pushed against the radius
or mass prior limits in other atmosphere cases. In this case,
there might be no clear trend between the two atmosphere
models for the same star, as seen in the fully versus partially
ionized hydrogen comparison for PSR J0030+0451, where
partial ionization results in significantly smaller or higher NS
radii than those of full ionization depending on the background
assumption.<fnref n="13" /> On the positive side, the sensitivity to the
atmosphere models allows us also to discriminate between the
models with lower evidence or with unexpected values for the
inferred parameters. However, the favored atmosphere models
and the inferred values may not necessarily remain the same
when external constraints or more complex spot patterns are
implemented.

5. Conclusions

We have explored the sensitivity of the constraints that we
obtain on NS parameters, especially the radius, to the
assumptions made for the NS atmosphere. We found that none
of our cases had any significant impact on the inferred radius
and other properties of PSR J0740+6620, possibly because of
its relatively low signal-to-noise ratio, its inferred viewing
geometry, or its external parameter constraints. However, the
results for PSR J0030+0451 were found to be potentially
sensitive to both the chemical composition of the atmosphere
and the choice of whether partial ionization in the atmosphere
was accounted for or not. The largest difference would appear
if the atmosphere were to consist of helium instead of
hydrogen, which would shift the radius from around 10.5 to
14 km in the case of the ST-U hot region model (circular hot
regions), or from 13 to 15 km in the case of ST+PST hot region
model (a circle and a crescent). In both cases, the difference in
the evidence is not enough to distinguish between the models,
although hydrogen composition is expected to be a priori more
likely based on evolutionary arguments.

When constraining the PSR J0030+0451 background using
XMM-Newton (as in Vinciguerra et al. 2023b) with the ST-U
model, the inferred radii using fully ionized hydrogen and
helium are consistent with each other but close to the prior
upper limit of 16 km. This consistency could be, at least partly,
caused by the different inferred hot region geometry, where
more radiation is emitted at smaller emission angles, for which
the fractional differences between the atmosphere models are
smaller. However, the corresponding partially ionized hydro-
gen case predicts a significantly smaller NS radius (around
11 km) even with a similar hot region geometry. Thus, more
studies with accurate background and other constraints may be

<footnote n="13">¹³ It should be also noted that the adequate sampling of the parameter space has not been rigorously proved for all the different atmospheres, hot spot shapes, and background cases.</footnote>
needed to properly resolve the typically multimodal structure of
the posterior surface for PSR J0030+0451.

In addition, we also found atmosphere cases that either did
not significantly impact the PSR J0030+0451 and PSR J0740
+6620 results or were shown to be unlikely based on the
inferred evidence. The former includes the runs with an
externally heated model with certain assumptions for the
return-current energy distribution, and runs with parameterized
uncertainty in atmospheric beaming pattern. The latter includes
the runs with carbon composition, which does not fit the data as
well as the other models.

Acknowledgments

This work was supported in part by NASA through the NICER mission and the Astrophysics Explorers Program. T.S., S.V., D.C., A.L.W., Y.K., and B.D. acknowledge support from European Research Council (ERC) Consolidator Grant (CoG) No. 865768 AEONS (PI: Watts). This work used the Dutch national e-infrastructure with the support of the SURF Cooperative using grant No. EINF-4664. Part of the work

was carried out on the HELIOS cluster including dedicated
nodes funded via the above-mentioned ERC CoG. The
computer resources of the Finnish IT Center for Science
(CSC) and the FGCI project (Finland) are also acknowledged
for part of the new atmosphere grid computation. W.C.G.H.
appreciates the use of computer facilities at the Kavli Institute
for Particle Astrophysics and Cosmology and acknowledges
support through grant 80NSSC22K0397 from NASA. S.G.
acknowledges the support of the Centre National d'Etudes
Spatiales (CNES). S.M.M. thanks NSERC for its support. T.S.
would like to thank Juri Poutanen, Thomas Riley, Nathalie
Degenaar, Cole Miller, Joonas Nättilä, and the anonymous
referee for useful discussions and comments.

**Appendix**
**Additional Figures**

In Figures 8–12, we show the posterior distributions for the
NS radii together with other model parameters, not including
those presented in Figures 3–7 of the main text.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Effect of the atmosphere model on the posterior distributions of radius and the other parameters not including those shown in Figure 3 (except the beaming parameters *a* and *b*) using the PSR J0740+6620 NICER data set conditional on the ST-U model. The additional parameters are cosine Earth inclination to spin axis (*cos* *i*), primary region initial phase (*φ*<sub>*p*</sub>), primary region center colatitude (*Θ*<sub>*p*</sub>), primary region angular radius (*ζ*<sub>*p*</sub>), secondary region initial phase (*φ*<sub>s</sub>), secondary region center colatitude (*Θ*<sub>s</sub>), secondary region angular radius (*ζ*<sub>s</sub>), distance (*D*), effective area scaling factor of NICER (*α*<sub>NICER</sub>), and interstellar neutral hydrogen column density (*N*<sub>H</sub>). See the caption of Figure 3 for more details.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Effect of the atmosphere model on the posterior distributions of radius and the other parameters not including those shown in Figure 4 using the PSR J0030+0451 NICER data set conditional on the ST-U model. The parameters are the same as described in the caption of Figure 8, except that *D* and *α<sub>NICER</sub>* are replaced with *β* = *α<sub>NICER</sub>*D<sup>-2</sup>. See the caption of Figure 4 for additional details.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Effect of the atmosphere model on the posterior distributions of radius and the other parameters not including those shown in Figure 5 using the PSR J0030+0451 NICER data set conditional on the ST-U model and energy-dependent beaming parameterization. The shown parameters are the same as in Figure 9. See the caption of Figure 5 for additional details.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Effect of the atmosphere model on the posterior distributions of radius and the other parameters not including those shown in Figure 6 using the PSR J0030+0451 NICER data set conditional on the ST+PST model. In addition to the parameters defined in Figures 8 and 9, we show the posteriors for azimuth offset between the emitting and the masking spherical caps of the secondary hot spot ($\chi_a$), colatitude of the masking region of the secondary hot spot ($\Theta_{o,x}$), and radius of the masking region of the secondary hot spot ($\zeta_{o,x}$). See the caption of Figure 6 for additional details.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Effect of the atmosphere model on the posterior distributions of radius and the other parameters not including those shown in Figure 7 using the PSR J0030+0451 joint NICER and XMM-Newton (NxX) data sets conditional on the ST-U model. In addition to the parameters defined in Figure 8, we show the posteriors for the effective area scaling factor of XMM-Newton ($\alpha_{XMM}$). See the caption of Figure 7 for additional details.</figcaption></figure>

## ORCID iDs

Tuomo Salmi https://orcid.org/0000-0001-6356-125X

Serena Vinciguerra https://orcid.org/0000-0003-3068-6974

Devarshi Choudhury https://orcid.org/0000-0002-2651-5286

Anna L. Watts https://orcid.org/0000-0002-1009-2354

Wynn C. G. Ho https://orcid.org/0000-0002-6089-6836

Sebastien Guillot https://orcid.org/0000-0002-6449-106X

Yves Kini https://orcid.org/0000-0002-0428-8430

Bas Dorsman https://orcid.org/0000-0002-9407-0733

Sharon M. Morsink https://orcid.org/0000-0003-4357-0575

Slavko Bogdanov https://orcid.org/0000-0002-9870-2742

## References

Alcock, C., & Illarionov, A. 1980, ApJ, 235, 534

AlGendy, M., & Morsink, S. M. 2014, ApJ, 791, 78

Arons, J. 1981, ApJ, 248, 1099

Badnell, N. R., Bautista, M. A., Butler, K., et al. 2005, MNRAS, 360, 458
Bauböck, M., Psaltis, D., & Özel, F. 2019, ApJ, 872, 162

Baym, G., Hatsuda, T., Kojo, T., et al. 2018, RPPh, 81, 056902

Bilous, A. V., Watts, A. L., Harding, A. K., et al. 2019, ApJL, 887, L23

Blaes, O. M., Blandford, R. D., Madau, P., & Yan, L. 1992, ApJ, 399, 634

Bogdanov, S., Dittmann, A. J., Ho, W. C. G., et al. 2021, ApJL, 914, L15

Bogdanov, S., & Grindlay, J. E. 2009, ApJ, 703, 1557

Bogdanov, S., Guillot, S., Ray, P. S., et al. 2019a, ApJL, 887, L25

Bogdanov, S., Lamb, F. K., Mahmoodifar, S., et al. 2019b, ApJL, 887, L26

Brown, E. F., Bildsten, L., & Chang, P. 2002, ApJ, 574, 920

Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125

Cerutti, B., & Beloborodov, A. M. 2017, SSRv, 207, 111

Chang, P., & Bildsten, L. 2004, ApJ, 605, 830

Chen, A. Y., & Beloborodov, A. M. 2014, ApJL, 795, L22

Colgan, J., Kilcrease, D. P., Magee, N. H., et al. 2016, ApJ, 817, 116

Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601

Gendreau, K., 2020 NICER Background Estimator Tools, https://heasarc.gsfc.nasa.gov/docs/nicer/tools/nicer_bkg_est_tools.html

Gendreau, K. C., Arzoumanian, Z., Adkins, P. W., et al. 2016, Proc. SPIE, 9905, 99051H

Haakonsen, C. B., Turner, M. L., Tacik, N. A., & Rutledge, R. E. 2012, ApJ, 749, 52

Hameury, J. M., Heyvaerts, J., & Bonazzola, S. 1983, A&A, 121, 259

Harding, A. K., & Muslimov, A. G. 2001, ApJ, 556, 987

Harding, A. K., & Muslimov, A. G. 2011, ApJ, 743, 181

Hebeler, K., Lattimer, J. M., Pethick, C. J., & Schwenk, A. 2013, ApJ, 773, 11

Ho, W. C. G., & Heinke, C. O. 2009, Natur, 462, 71

Ho, W. C. G., & Lai, D. 2001, MNRAS, 327, 1081

Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943

Kalapotharakos, C., Wadiasingh, Z., Harding, A. K., & Kazanas, D. 2021, ApJ, 907, 63

Lattimer, J. M. 2012, ARNPS, 62, 485

Lattimer, J. M., & Prakash, M. 2001, ApJ, 550, 426

Lo, K. H., Miller, M. C., Bhattacharyya, S., & Lamb, F. K. 2013, ApJ, 776, 19

Medin, Z., & Lai, D. 2007, MNRAS, 382, 1833

Miller, M. C., & Lamb, F. K. 1998, ApJL, 499, L37

Miller, M. C., & Lamb, F. K. 2015, ApJ, 808, 31

Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJL, 887, L24

Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2021, ApJL, 918, L28

Morsink, S. M., Leahy, D. A., Cadeau, C., & Braga, J. 2007, ApJ, 663, 1244

Nättilä, J., & Pihajoki, P. 2018, A&A, 615, A50

Özel, F., & Psaltis, D. 2009, PhRvD, 80, 103003

Pechenick, K. R., Ftaclas, C., & Cohen, J. M. 1983, ApJ, 274, 846

Philippov, A. A., & Spitkovsky, A. 2018, ApJ, 855, 94

Potekhin, A. Y. 2014, PhyU, 57, 735

Poutanen, J., & Gierliński, M. 2003, MNRAS, 343, 1301

Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021, ApJL, 918, L29

Riley, T. E. 2019, PhD thesis, Univ. of Amsterdam, https://hdl.handle.net/11245.1/aa86fcf3-2437-4bc2-810e-cf9f30a98f7a

Riley, T. E., Choudhury, D., Salmi, T., et al. 2023, JOSS, 8, 4977

Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJL, 887, L21

Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJL, 918, L27

Romani, R. W. 1987, ApJ, 313, 718

Ruderman, M. A., & Sutherland, P. G. 1975, ApJ, 196, 51

Salmi, T., Suleimanov, V. F., Nättilä, J., & Poutanen, J. 2020, A&A, 641, A15

Salmi, T., Suleimanov, V. F., & Poutanen, J. 2019, A&A, 627, A39

Salmi, T., Vinciguerra, S., Choudhury, D., et al. 2022, ApJ, 941, 150

Salmi, T., Vinciguerra, S., Choudhury, D., et al. 2023, Atmospheric Effects on Neutron Star Parameter Constraints with NICER, v1.0.0, Zenodo, doi:10.5281/zenodo.7449786

Strohmayer, T. E. 1992, ApJ, 388, 138

Suleimanov, V., Poutanen, J., & Werner, K. 2012, A&A, 545, A120

Suleimanov, V. F., Klochkov, D., Pavlov, G. G., & Werner, K. 2014, ApJS, 210, 13

Suleimanov, V. F., Poutanen, J., & Werner, K. 2020, A&A, 639, A33

Timokhin, A. N., & Arons, J. 2013, MNRAS, 429, 20

Vinciguerra, S., Salmi, T., Watts, A. L., et al. 2023a, arXiv:2308.08409

Vinciguerra, S., Salmi, T., Watts, A. L., et al. 2023b, arXiv:2308.09469

Watts, A. L., Andersson, N., Chakrabarty, D., et al. 2016, RvMP, 88, 021001

Wijngaarden, M. J. P., Ho, W. C. G., Chang, P., et al. 2019, MNRAS, 484, 974

Wijngaarden, M. J. P., Ho, W. C. G., Chang, P., et al. 2020, MNRAS, 493, 4936

Wolff, M. T., Guillot, S., Bogdanov, S., et al. 2021, ApJL, 918, L26

Zavlin, V. E., & Pavlov, G. G. 2002, in Proc. of the 270th Heraeus Seminar on Neutron Stars, Pulsars, and Supernova Remnants, ed. W. Becker, H. Lesch, & J. Trümper (Garching bei München: Max-Plank-Institut für extraterrestrische Physik), 263

Zavlin, V. E., Pavlov, G. G., & Shibanov, Y. A. 1996, A&A, 315, 141