Atmospheric Effects on Neutron Star Parameter Constraints with NICER

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 degrees, on the other hand, have no significant effect on the inferred radius.


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 temperature6 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 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.
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 e , with most of the posterior probability above 2.0 M e (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 T log K 6 10 eff ( )~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.

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. (2019Riley et al. ( , 2021)), Salmi et al. (2022), Vinciguerra et al. (2023a), andVinciguerra et al. (2023b).For all posterior computations, we use the X-ray Pulse Simulation and Inference7 (X-PSI) code, with versions ranging from v0.7.3 to v1.2.1 (Riley et al. 2023). 8The 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.

Atmosphere Modeling
We employ atmosphere models to determine the specific intensity I(E, μ) emitted from the NS surface, where E is the photon energy and μ 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, μ) 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.

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, T log K 10 eff ( ), the logarithm of the surface gravity, g log cm s 10 s 2 ( ) -, the cosine of the emission angle, μ, and the logarithm of photon energy with respect to the temperature, E kT log 10 eff ( ) (see Bogdanov et al. 2021, for more details).The ranges and the number of grid points for each 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 deepheated partially ionized hydrogen NSX model (nsxHp), applied previously in Miller et al. (2021),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), andSuleimanov 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 10 min g = , powerlaw slope δ = 2, and a maximum Lorentz factor max g = 2 10 5 ´.These choices correspond to the fiducial parameter setup in Salmi et al. (2020), and they lead to an average Lorentz factor γ 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.

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, μ) NUM (interpolated from the atmosphere table) with a polynomial function modifying the beaming: a, b, c, and d are beaming parameters.In the case of an energyindependent 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 ò m m m = , which gives the following: This constant is chosen to minimize the change in the total flux spectrum for any given a, b, c, and d.This choice precisely conserves the emergent spectrum when the model intensity I(E, μ) 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 separately for each emitted photon (increasing the computational Notes.The atmosphere grid ranges ([L, L ]) and number of grid points (:L) 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 μ = 0, and μ = 1.In the case of hatmH, the grid points were based on the sample points of Gauss-Legendre quadrature between 0 and 1. a Originally a logarithmically spaced grid of E from 4.1 × 10 −4 to 41 keV with 360 grid points was used for calculating this model, but the intensity values were then interpolated to the same E kT log 10 eff ( ) grid as in NSX, for simplicity.
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, μ) − 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.

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 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 × 10 3 live points and 0.1 sampling efficiency.The R 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).
10 These checks were done for 100 (or 50) emission angles linearly spaced in μ = [0.0,1.0] (or μ = [0.5, 1.0]) and for 100 energies linearly spaced in , using T eff of both the primary and secondary hot regions. 11These 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 μ and μ 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.
enough sampling resolution of the parameter space (see Section 2.4).
For the instrument response model, we use only energyindependent 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 β 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).

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

PSR J0740+6620
For PSR J0740+6620 analysis, we used the same preprocessed 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).

PSR J0030+0451
For PSR J0030+0451 analysis, we used the same preprocessed 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).

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 Notes.See the explanation for columns and their contents in Table 2. Unlike the previous table, here, all the runs were performed with 1 × 10 4 live points and 0.3 sampling efficiency.Note that the evidencefor models with different background constraints is not comparable.a 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.
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 "lowresolution" settings from Riley et al. (2021): 4 × 10 3 live points; and 0.1 sampling efficiency12 ; and an estimated remaining logevidence of 10 −1 .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.4km.
For the PSR J0030+0451 analyses, we used slightly different resolution settings; 1 × 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 × 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 modeseparation 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.

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 T log K 6.1059 .The different spectra are presented in Figure 1 for μ = 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 hydrogen and helium composition) only at higher temperatures, such as in the T log K 6.4 10 eff ( )= 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 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.

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.

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 energyindependent beaming parameters (as explained in Section 2.1.2) for fully ionized hydrogen does not have any significant impact on the results.

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

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 small shift (smaller than ∼0.03 in T log K 10 [ ]) 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 e , and the temperatures of both hot regions become higher (again as expected based on the spectral shapes in Figure 1).As seen 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 0.65 0.77 -+ 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  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  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 0.92 1.29 -+ (for hydrogen) to 14.00 1.88 1.39 -+ km (for helium).The inferred mass also increases to around 1.5 M e .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.

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, 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.

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 e (which is still smaller than the 2.7 M e 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 × 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).

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 e , 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.

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.

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 × 10 6 versus 0.6 × 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.

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 e , 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. 2019Wijngaarden et al. , 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 e (for ST+PST), rather than 2.7 M e .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.

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 T log K 6.1 10 eff ( )» 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 T log K 6.1 10 eff ( )» , the opacity uncertainties should still be small, but larger at T log K 6.4 10 eff ( )» , 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.

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 returncurrent 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 T log 6.4 10 eff = ), 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 atmosphere 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).

Beaming Uncertainty
Ideally, the uncertainty in the atmosphere model assumptions 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 dependency 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 atmosphere.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.

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. 13On 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.

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 hydrogen 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 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.

Figure 1 .
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 ( T log K 6.1059 10 eff ( )= , g log cm s 14.143 10 s 2 ( )= -).The different atmosphere models are explained in Table 1.Lower panel: the same as the upper panel, except with a higher temperature T log K 6.4 10 eff ( )= , and surface gravity fixed to g log cm s 14.15 10 s 2 ( )= -.The spectra for best-fit parameters from PSR J0740 + 6620 analysis are available.(The complete figure set of 3 images is available.)

Figure 2 .
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 T log K 6.4 10 eff ( )= , and surface gravity fixed to g log cm s 14.15 10 s 2 ( )= -.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.)

Figure 3 .
Figure3.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 Table1for 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 KL are reported for the fully ionized hydrogen NSX case: ST-U-nsxH (see Table2for 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 Figure5of 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.

Figure 4 .
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 Table1for definitions).The 1D credible intervals and the KL-divergence estimates are reported for the fully ionized hydrogen NSX case: ST-U-nsxH (see Table3for mass and radius values in other cases).See the caption of Figure3for additional details about the figure elements.Posterior distributions for the other parameters are shown in Figure9of theAppendix.
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 Table1for definitions).The 1D credible intervals and the KL-divergence estimates are reported for the fully ionized hydrogen NSX case: ST-U-nsxH (see Table3for mass and radius values in other cases).See the caption of Figure3for additional details about the figure elements.Posterior distributions for the other parameters are shown in Figure9of theAppendix.

Figure 5 .
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.

Figure 6 .
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.

Figure 7 .
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 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.

Figure 8 .
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 ( i cos ), primary region initial phase (f p ), primary region center colatitude (Θ p ), primary region angular radius (ζ p ), secondary region initial phase (f s ), secondary region center colatitude (Θ s ), secondary region angular radius (ζ s ), distance (D), effective area scaling factor of NICER (α NICER ), and interstellar neutral hydrogen column density (N H ). See the caption of Figure 3 for more details.

Figure 9 .
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 α NICER are replaced with β = α NICER D −2 .See the caption of Figure 4 for additional details.

Figure 10 .Figure 11 .
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.

Table 1
Summary of the Numerical Atmosphere Tables

Table 2
Summary Table for PSR J0740+6620 Runs

Table 3
Summary Table for PSR J0030+0451 Runs