Mass Distribution and Maximum Mass of Neutron Stars: Effects of Orbital Inclination Angle

Matter at ultra-high densities finds a physical realization inside neutron stars. One key property is their maximum mass, which has far-reaching implications for astrophysics and the equation of state of ultra dense matter. In this work, we employ Bayesian analysis to scrutinize the mass distribution and maximum mass threshold of galactic neutron stars. We compare two distinct models to assess the impact of assuming a uniform distribution for the most important quantity, the cosine of orbital inclination angles ($i$), which has been a common practice in previous analyses. This prevailing assumption yields a maximum mass of $2.25$~$M_\odot$ (2.15--3.32~$M_\odot$ within $90\%$ confidence), with a strong peak around the maximum value. However, in the second model, which indirectly includes observational constraints of $i$, the analysis supports a mass limit of $2.56^{+0.87}_{-0.58}~M_\odot$ ($2\sigma$ uncertainty), a result that points in the same direction as some recent results gathered from gravitational wave observations, although their statistics are still limited. This work stresses the importance of an accurate treatment of orbital inclination angles, and contributes to the ongoing debate about the maximum neutron star mass, further emphasizing the critical role of uncertainties in the individual neutron star mass determinations.


Introduction
More than fifty years after the detection of the first pulsar [1], these compact objects continue to challenge our understanding.The study of neutron stars (hereafter NSs), of which pulsars are a subset, provides insights into fundamental aspects of the Universe, ranging from the state of matter at ultra-high densities to interstellar enrichment with heavy elements.Additionally, they offer the most precise tests for General Relativity (GR) [2] and contribute to the understanding of stellar evolution and binary interactions.The importance of studying such stars is clear, and the last 10-15 years of research have brought about paradigm shifts, correcting previously unquestioned beliefs in the field.
The proposal that NSs are formed at the end of the life of massive progenitor stars, initially suggested by Baade and Zwicky [3], coupled with their known associations with supernova remnants [4], helped construct an evolutionary scenario for their formation.Despite an initial variety of ideas [5], the prevailing notion was that NSs are formed from the collapse of an Iron ( 56 Fe and neighbour isotopes) core, developed at the center of progenitor stars with initial masses ranging from 8 to 25 M ⊙ .A misleading association of the iron core before the collapse with the Chandrasekhar limit, together with early mass measurements, contributed to establishing a paradigm of a unique formation channel for NSs, with a mass scale around 1.4 M ⊙ and minimal dispersion [6].
Simultaneously, the existence of a mass limit predicted by the Tolman-Oppenheimer-Volkoff (TOV) equation, derived from the theory of GR and not from fermion degeneracy, sparked extensive discussions.Since the determination of this limit relies on the equation of state (EoS) describing matter at extreme conditions, and this EoS cannot be securely determined from first principles or terrestrial experiments, an accurate maximum mass could not be determined.In 1974, Rhoades and Ruffini [7] proposed an "absolute" limit under the assumption of a static and isotropic metric.Assuming the equation of state above a certain density to be the stiffest possible 1 , they established an upper threshold for NS masses at 3.2 M ⊙ , which has since been adopted to distinguish NSs from black holes (BHs).While uncertainties in early X-ray mass measurements did not forbid the existence of high masses, and some theoretical studies suggested the possibility of a few physical EoSs leading to a maximum mass around 2 M ⊙ , it became a consensus in the scientific community that, for evolutionary reasons, NS masses should not exceed the "canonical" value of 1.4 M ⊙ , in agreement with the first precise mass measurements [6,8].Or, in other words, that the 1.4 M ⊙ was the value imprinted at birth by collapse physics.
However, observational efforts have steadily increased the number of measured masses over the years.For over a decade now, it has been known that the mass range covered by NSs is much larger than previously expected, with the current interval spanning from 1.17 M ⊙ to values exceeding 2.0 M ⊙ , a much broader range than was previously thought possible.The first pulsar discovered with a mass that deviates significantly from the "typical" value was Vela X-1, with 1.86 ± 0.16 M ⊙ [9].In the following years, a few other potentially massive NSs started being discovered, as is the case of PSR J0751+1807 (m = 2.1 ± 0.2 M ⊙ ) [10], PSR B1516+02B (m = 2.08 ± 0.19 M ⊙ ) [11], PSR J1748-2021B (m = 2.74 ± 0.21 M ⊙ ) [12], and PSR J1614-2230, with m = 1.97 ± 0.04 M ⊙ [13], which was considered the most accurate inference among the massive ones at that time.The discovery of PSR J0740+6620 [14], with m = 2.14 +0.10  −0.09 M ⊙ , proved once and for all the existence of NS masses >2.0 M ⊙ , highlighting the issue of the maximum mass.
From an evolutionary standpoint, the presence of a broad mass range was soon associated with different formation mechanisms.An analysis by Schwab et al. [15] with a selected sample of well-constrained masses (uncertainties less than ≲0.025M ⊙ ) revealed a double-peaked distribution, where a group of NSs centered at ∼1.35 M ⊙ was associated with the "standard" Fe core-collapse supernova, while the second group clustered around ∼1.25 M ⊙ , linked to the electron-capture supernova scenario, was expected to occur in degenerate cores of O − Ne − Mg [5,16].Although subsequent analyses with the complete sample of NS masses did not detected this lower peak, it is highly probable that it occurs for progenitor stars with initial masses between 8 and 10 M ⊙ [17], exploding via electron capture onto a degenerate O − Ne − Mg core.Instead, these analyses found a second peak around 1.75-1.80M ⊙ [18][19][20][21], in addition to the dominant one, at ∼1.4 M ⊙ .This massive group was rapidly associated with accretion onto NSs in binary systems, as expected from the existence of millisecond pulsars and other systems containing an NS.Even though accretion is still a prime candidate for the origin of massive NSs, recent works have also demonstrated the possibility of forming massive NSs directly from heavier iron cores [22], while others have shown that smaller Fe cores can collapse, leaving behind a 1.17 M ⊙ NS [23].
One particular class of binary systems, the "spiders", are prime candidates to populate the high-mass interval."Spiders" are close binary systems (P b < 1 day), where a millisecond pulsar orbits a low-mass companion star that is in the process of having its envelope ablated away by the pulsar's wind 2 .If the donor companion has a mass of 0.1-0.5 M ⊙ , it is classified as a redback, while those with companion masses ≲ 0.05 M ⊙ -and where the accretion has stopped-are named black widows.These systems are expected to experience a large accretion phase in the early-stage phases, accumulating a great amount (>0.8 M ⊙ ) of mass in the most extreme cases [24].Evidence suggests that these systems can host the most massive NSs in the Universe [25,26].The most recent massive spider discovered is PSR J0952-0607, with m = 2.35 ± 0.17 M ⊙ [27], placing the maximum NS mass as >2.19 M ⊙ (1σ confidence).
Finally, NSs may also be formed through the (single-degenerate) accretion-induced collapse (AIC) scenario, where a massive white dwarf (WD) exceeds its mass limit and collapses without igniting carbon, or the double-degenerate AIC, where two WDs merge.Although these events have never been positively identified, population synthesis has yielded the expectation of approximately 10 7 pulsars formed by AIC in the single-degenerate channel and a few times this figure comes from the double-degenerate channel [28].Additionally, the long gamma-ray burst GRB 211211A fed the possibility of an NS-WD merger [29], leaving behind a magnetar, another feasible scenario to populate the second peak of masses around ∼1.8 M ⊙ .
The detection of gravitational waves (GW) from the merger of two NSs in 2017, GW170817 [30], opened a new era in multimessenger and multiwavelength astronomy, providing a new tool to set constraints on NS physics.Together with galactic determinations, constraints placed by GW170817 contributed to placing the limit of masses below 2.2-2.3M ⊙ [31][32][33][34][35]. Later on, the detection of GW190814 [36], where one of the components that merged with a 23 M ⊙ BH had a mass of 2.59 +0.08 −0.09 M ⊙ , raised a tension about whether the threshold of the maximum mass should be above or below 2.2-2.3M ⊙ .The work of Nathanail et al. [37] concluded that the secondary of GW190814 needs to be a BH in order to avoid contradictions with the post-merger observations of GW170817.At the same time, an analysis made by the Ligo-Virgo (LV) Collaboration [38] of the binary BH (BBH) merging population in the second catalog found GW190814 to be an outlier, i.e., to be likely originated from an NS-BH merger.In addition, from an analysis of the available LV NS-BH events-and assuming GW190814 to be one of them-the maximum mass of a non-spinning NS was found to be 2.7 +0.5 −0.4 M ⊙ [39].A recent result combined constraints from all possible astronomical approaches, and derived a maximum mass in the range of 2.49-2.52M ⊙ [40], providing independent evidence supporting the existence of extremely heavy NS masses.
In this article, we focus on analyzing the up-to-date sample of NS masses found in binary systems in the Galaxy and in Globular Clusters, and the impact the orbital inclination angle has on the conclusions derived from this kind of analysis, especially in the estimation of the maximum mass (m max ).We start in Section 2 with an overview of all available measurement methods and a discussion of the detailed features leading to individual masses and uncertainties.In Section 3, we describe the statistical method and the "accuracy-dependent" and "accuracy-independent" models used to analyze the whole distribution, now featuring 125 members (listed in Appendix A).In Section 4, we present our results and compare the posterior distributions of each model.Finally, we draw some general conclusions in the last section.It is important to emphasize the difference between the TOV mass (M TOV ), supported by non-rotating pulsars, and the m max , which is derived from our analysis and is the maximum supported by pulsars with maximum uniform rotation 3 .Furthermore, according to the quasi-universal relation derived in Most et al. [41], m max can exceed by a factor of up to 1.2 the M TOV value.

Neutron Star Mass Measurements
The majority of NSs are observed as pulsars, rapidly-rotating and highly magnetized NSs emitting beams of radiation along their magnetic axes.This emission is observed on Earth as a pulse due to the lighthouse effect.The remarkable regularity of these pulses, with pulsars being one of the most stable clocks in the observable universe, makes pulsar timing the most accurate method for determining their masses and testing fundamental physics.This procedure involves monitoring the times of arrival (ToAs) of pulses over an extended period of time to determine the pulsar's rotation period.Thanks to the regularity, small deviations in ToAs are detectable with precision and are indicative of the presence of a companion.The greater the number of collected ToAs, the greater will be the precision achieved.Hence, several years of observations are necessary to achieve a high precision [42].
Currently, more than 3300 radio pulsars have been identified (see an online catalogue at ATNF [43]), yet only a few of their features can be directly inferred from observations, and mass measurements are possible for only a small fraction of the total sample.The Neutron star Interior Composition ExploreR (NICER), a telescope placed on board the International Space Station in 2017, facilitates timing-and rotation-resolved spectroscopy of thermal and non-thermal X-ray emissions from NSs.Recently, NICER enabled precise measurements of radii and masses for two pulsars, namely PSR J0740+6620 and PSR J0030+0451 [44,45].Despite this being a promising method, especially for isolated pulsars, the dominant means for inferring NS masses continues to be the study of orbital motion in binary systems determined through pulsar timing of radio sources.
The ToAs reveal the orbital properties of the system expressed in terms of Keplerian parameters: orbital period (P B ), eccentricity (e), semi-major axis projection onto the line of sight (x = a sin i), and the time (T 0 ) and longitude (ω) of periastron.From Kepler's third law, a mass function containing the pulsar mass (m p ), the companion mass (m c ), and the inclination angle (i) of the system can be written down as: where T ⊙ ≡ GM ⊙ /c 3 = 4.925490947 µs (G is the gravitational constant, c is the speed of light, and M ⊙ is the solar mass).If the mass function of the pulsar and the companion are measured ( f p , f c ), along with the mass ratio (q = m p /m c ), the individual masses of the system can be determined, provided the inclination angle i is known.However, until now we were only able to set precise constraints for both masses in a few particular cases, because of the difficult estimation of the inclination angle i, which we now address.

Orbital Inclination Angle
Since binary stars have small angular sizes, it is challenging to fully resolve the orbit.Consequently, if the radial and transverse velocities of the system cannot be measured, the orbital inclination of the orbit with respect to the plane of the sky (i) cannot be directly determined.
However, these systems often display variability of the companion light curves, from which an inclination-dependent radial velocity can be estimated.Pulsar mass estimates via radial velocity measurements depend on the inclination angle i of the orbit with respect to the line of sight, through the relation: Therefore, a high uncertainty in the inclination angle will result in a high uncertainty in the mass measurement.Although it is difficult to precisely determine i, constraints can be set if, for example, an eclipse is observed through spectroscopy of the companion.

Relativistic Binaries
In the special case of compact binaries, where the companion is a WD or an NS, relativistic effects may influence the orbit and can sometimes be measured.These effects are described in terms of the so-called post-Keplerian (pK) parameters, defined as: 1.
Orbital period decay, Ṗb : Range of Shapiro delay, r: Shape of Shapiro delay, s: . "Einstein delay", γ: Advance of periastron, ω: If at least two pK parameters are measured, the component masses can be individually determined.When more pK parameters are measured, it is also possible to test GR with very high precision, as demonstrated in Kramet et al. [2].Accretion torques in binary systems can circularize the orbits, resulting in many NS binaries with extremely low eccentricities, hampering the measurement of ω and γ.On the other hand, the Shapiro delay of the pulses, caused by the gravitational field of the companion, depends on the orbital inclination and is typically relevant for systems with high inclinations.Lastly, orbital decay due to gravitational wave radiation ( Ṗb ) is only measurable for very tight orbits.All of these conditions make pulsar mass measurements a challenging task.If relativistic effects are too small, they can go undetected even after years of pulsar timing.

Shapiro Delay
The Shapiro delay is the increase in light travel time through the curved space-time near a massive body.In binary pulsar systems that have highly inclined (nearly edge-on) orbits, excess delay in the pulse ToA can be observed when the pulsar is situated almost behind the companion during orbital conjunction.In combination with the mass function, the Shapiro delay offers one of the most precise methods with which to directly infer the mass of NSs.

Optical Spectroscopy
On the other hand, when the pulsar has an optically bright low-mass companion, such as a Main Sequence or post-Main Sequence star or a WD, phase-resolved spectroscopy of the companion can yield the orbital radial velocity amplitude (K c ).When combined with x and P b , this provides the binary mass ratio q = (m p /m c ) = (K c /K p ).In the case of WD companions, their radii can be estimated when the distance (d) to Earth is known, along with the optical flux (F O ) and effective temperature (T e f f ) measurements, as: where σ is the Stefan-Boltzmann constant.Their masses can thus be estimated by combining the effective temperature with the surface gravity obtained from an atmosphere model, which provides a model-dependent method.Combining m c with q, the pulsar mass is addressed.
Neutron star masses in spider systems are typically inferred through spectrophotometric methods.The system is filled with intra-binary material, causing the radio pulsation to be scattered and absorbed [46].As a consequence, their optical light curves are sensitive not only to the orbital inclination, but also to the heating models of the companion's surface, which can be challenging to predict.A large systematic error in the inclination angle estimate can result in a significant bias in the mass estimate for this class of pulsars.

Gamma-Ray Pulsars
Millisecond pulsars also present gamma-ray pulsations [47].In contrast to other wavelengths, it seems unlikely that gamma-rays are absorbed in the diffuse intra-binary material of spider systems.Consequently, the observed gamma-ray eclipses are potentially associated only with the occultation of the pulsar by the companion, providing a more robust determination of the inclination angle.The work of Clark et al. [48] conducted a search for gamma-ray eclipses in 49 confirmed and candidate spider systems, with significant detections in five of them, from which mass determinations were obtained (Table 1).In only one of these five systems, the inclination angle was found to be inconsistent with previous optical modeling.
For PSR B1957+20, photometric observations provided an estimate of 63 • ≲ i ≲ 67 • , resulting in a best fit with a high mass of m = 2.4 ± 0.1 M ⊙ [49].Observations in gamma-rays, on the other hand, require an inclination angle i > 84.1 • , corresponding to a significantly lower mass of m = 1.81 ± 0.07 M ⊙ .If this discrepancy is found to be consistent for the most massive objects measured through X-ray/Optical modeling, as well as confirmed for this particular system, it could have a substantial impact on the determination of the maximum mass of NSs and, consequently, on our understanding of matter at ultra-high densities.
However, it is essential to note that the results from Clark et al. [48] assume that gamma-ray photons detected at Fermi-LAT energies are not absorbed by the diffuse material and that the wind is isotropic (a questionable simplified assumption).Further studies on intra-binary shocks in spider systems are still necessary to firmly constrain their geometry and confirm such results.

Analysis of the Mass Distribution
As mentioned in the Introduction, studying the mass distribution of NSs can provide valuable insights into the evolutionary mechanisms leading to their formation and help constrain their maximum masses.The very first of these analyses was conducted on a small sample of eight NS masses from four double neutron star (DNS) systems, and employed Bayesian inference [6].Those results indicated that NS masses should predominantly fall within the range of 1.3 < m/M ⊙ < 1.6.Subsequent analyses with a larger sample of 19 NS masses arrived at a similar conclusion [8], with no evidence of a significant dispersion around the mean value of 1.35 ± 0.04 M ⊙ .
As the search for NSs continued and surveys yielded an increasing number of discovered objects, the landscape began to change.Over the years, different research groups employed frequentist and Bayesian inference techniques to extract information from the mass distribution of NSs [18][19][20][21]32,34,50,51].All these analyses consistently revealed the presence of a double-peaked distribution.These findings directly challenged the old idea of a single formation channel, since such a scenario could hardly account for the entire range of observed masses, suggesting that different processes are at work.
In this work, we employ a Bayesian analysis of an updated sample of NS masses using the advanced technique of Markov Chain Monte Carlo (MCMC) simulations.This technique allows us to determine the posterior distribution of a set of unknown parameters (θ) based on the a priori information available for each of these parameters, combined with information from observed data.According to Bayes' theorem, the posterior distribution is expressed as P(θ|d) = P(θ)L(d|θ), where P(θ) represents the a priori distribution, and L(d|θ) denotes the likelihood of the model.Assuming a parameterized model, the likelihood can be defined as L(d|θ) = P(d i |m i p )P(m i p |θ), where d i is the data.Consequently, the posterior distribution marginalized over pulsar masses (m i p ) is given by: To account for the double-peaked distribution, recent works have adopted a Gaussian mixture model for the model likelihood, P(m i p |θ), with n components.This model has been compared with other distribution families but, so far, no evidence has emerged to reject the preference for a Gaussian mixture model.Furthermore, if we wish to introduce a free parameter to estimate a cutoff in the distribution (motivated by the expectation of a maximum mass), we can employ a truncated Gaussian mixture model: with ∑ n j r j = 1 to ensure normalization.The set of model parameters we aim to infer includes the mean (µ j ) and standard deviation (σ j ) of each Gaussian component, along with their respective amplitudes (r j ) and the maximum mass (m max ) set as a free parameter, θ = {µ j , σ j , r j , m max }, with j = {1, n}.In the following subsections, we will delve into the results of our analysis assuming two different individual pulsar distributions, P(d i |m i p ).To sample our posterior distributions, we employ an MCMC algorithm with STAN [52].

"Accuracy-Dependent" Model
As discussed in Section 2, the most accurately determined pulsar masses are those where relativistic effects are revealed, and two or more pK parameters are observed.However, in general, satisfying these requirements is not straightforward, and in many instances only mass limits can be established.In this work, we employed the models used in the works of Alsing et al. [32], Shao et al. [34], Antoniadis et al. [50], which we refer to as accuracy-dependent models, where the pulsar mass likelihood depends on whether or not two or more pK parameters are measured, as we describe below.
For systems where the pulsar mass is accurately determined from observations, we assume the likelihood to be a normal distribution, N (m i , σ m i ), with mean and standarddeviation values provided in the sixth column of Table A1.In cases where, despite the mass function f , only the ω is constrained, it is possible to determine the total mass (m t ) of the system, given in the fourth column of the same table.In such cases, the pulsar's mass is marginalized as where mt , σ m t , and f are the measured total mass, its uncertainty, and the measured mass function, respectively.We assume that the total mass and the mass function are independent, and that m t has Gaussian uncertainties.Integrating the last line of the above equation over i with a flat a priori, the individual pulsar mass can be marginalized from: Note that the above equation is slightly different from Equation (3) in Alsing et al. [32], as suggested by Farr and Chatziioannou [53].Finally, if only the mass ratio (q) can be determined through phase-resolved optical spectroscopy, in addition to f , the likelihood of pulsar mass is given as: where q and σ q are the measured mass ratio and its uncertainty.Here, we assume that the mass ratio and the mass function are independent and that q has Gaussian uncertainties.Integrating the last line over i and assuming again a flat a priori, we are led to: The key assumption for deriving Equations ( 11) and ( 12) is that cos i follows an isotropic distribution, which means it is subject to a uniform a priori assumption, allowing i to span any value between 0 • and 90 • .

"Accuracy-Independent" Model
Given the requirement for a uniform distribution for cos i in the previous model, and considering the substantial impact that the orbital inclination angle has on the pulsar mass estimation (as discussed in Section 2.1 and shown later in Section 4), we looked for an alternative model that could avoid such an assumption.In the "accuracy-independent" model, we assume all pulsar masses listed in Table A1 to be modeled with a normal distribution, as performed in previous works [19,54,55].To illustrate the procedure, observations of PSR 2S 0921-630 (the first in our table), for example, found it to have a mass of 1.44 ± 0.10 M ⊙ [56].We then modeled this system with a normal distribution, with a mean value of 1.44 and a standard-deviation of 0.1, N (1.44, 0.10), and so on.Systems without any constraints on the individual mass were naturally excluded from this analysis.
Although to model all masses with Gaussians may appear less robust, it is worth noting that the pulsar mass values reported in Table A1 were calculated while taking into account all observational measurements and constraints of Keplerian and post-Keplerian parameters, which are tightly tied, available for each particular system, including i. Furthermore, Bayesian statistics is known to weight the sampling from data uncertainty, i.e., data with large uncertainties will have a lower weight in posterior distributions than those with small uncertainties.In this sense, mass determinations derived from spectrophotometry will have a lower weight in the posterior results than those derived from Shapiro delay, for example.

Results
The summaries of marginalized distributions for the "accuracy-dependent" and "accuracy-independent" models are presented in Tables 2 and 3, respectively.With the exception of m max in the last row, for which we report the mode value, the second column displays the mean value for each parameter, followed by its respective standard deviation.The third column shows the "highest posterior density interval" (HPDI) with 94% probability, indicating the shortest range of values that encompasses the given probability.
Table 2. Summary of marginal posterior distribution for the "accuracy-dependent" model.With the exception of m max , the second column displays the mean value, followed by respective standard deviation and the 94% highest posterior density interval, which defines the lowest interval that comprises 94% of the probability.For m max , the value shown in column 2 is the mode.As is apparent in both Tables, and in line with previous analyses available in the literature, both samplings yield a bimodal distribution featuring two distinct "groups" of NSs.The first group is centered around ∼1.35 M ⊙ with a small standard deviation of roughly 0.09 M ⊙ .In the second group, the NSs cluster at approximately 1.8 ± 0.26 M ⊙ .The predominance of objects in the first group is discernible from the amplitude values r.The real divergence between both models becomes apparent in the marginalized distribution of m max , which we will explore further below.

Mean
Figure 1 provides a visual representation of the uncertainties in each parameter, summarized in Table 2 (left panel) and Table 3 (right panel).In the figure, light grey lines depict 1000 posterior samples of the sampled pulsar mass distribution.The Maximum Posterior Probability, depicted in black, corresponds to the estimate that aligns with the mode of the posterior distribution.The blue line represents the average mass distribution, which is constructed assuming the mean values of each parameter.In contrast to the findings in Alsing et al. [32] and Shao et al. [34], our analysis did not reveal a sharp cut-off in the posterior distributions for either of the models.2 (left panel) and Table 3 (right panel).The blue curve is the posterior mean distribution and the black line is the maximum posterior distribution.

Effects on Individual Masses
To investigate the differences between the two models defined above, we analyzed the impact of a uniform distribution for cos i on the sampling of individual pulsar masses.In the "accuracy-dependent" model, if only f and m t or q are known, the individual masses of pulsars are marginalized from Equations ( 11) and (12).For the sake of simplicity, we present in the second and third columns of Table 4 the 95% HPDI of masses of four selected systems (namely the pulsars B1957+20, J1311-3430, B1516+02B and J1748-2021B) derived in the "accuracy-dependent" model, using Equations ( 11) and ( 12).These systems were known for a long time to be potentially extremely massive (≥2.0 M ⊙ ), and are important-although not exclusively-for the inference of m max .We now focus our attention on them to illustrate the impact of a flat distribution over cos i on posterior estimations.11) and (12), as functions of the total mass (second column) or mass ratio (third column).They are particularly interesting since their masses derived from observations are ≥2.0M ⊙ (fourth column).Below, we provide comments on the inferences for each of the four pulsars listed in Table 4.We compare sampled masses with the results obtained from spectrophotometry modeling, summarized in the fourth column.It is worth noting that, for all four pulsars, the masses sampled from the "accuracy-dependent" model are inconsistent with values derived from the original observations (last column).

PSR B1957+20
A black-widow system.Analysis of the light curve of the companion star yields a radial-velocity amplitude of K 2 = 353 ± 4 km s −1 .When combined with the pulsar's mass function, this measurement provides a minimum companion mass of m c,min = 0.022 M ⊙ .The best-fit values for the mass ratio and inclination angle are q = 69.2± 0.8 and i = 65 • ± 2 • , respectively.These parameters, when combined, result in a best-fit pulsar mass of m = 2.40 ± 0.12 M ⊙ [49], and a lower limit on the pulsar mass at 1.66 M ⊙ .As seen in the third column of Table 4, the "accuracy-dependent" model leads us to a marginalized mass between 1.16-2.14M ⊙ , inconsistent with spectroscopic determination.As we commented in Section 2.4, the work in [48] places a constraint of 50 • < i < 85 • , resulting in a mean mass at 1.81 ± 0.07 M ⊙ .Further investigation into the geometry and emissions of spider systems is necessary to pinpoint the correct mass value.

PSR J1311-3430
Until recently, the mass of this black-widow system was determined using a lightcurve analysis, resulting in a mass of m = 2.63 +0.3  −0.2 M ⊙ [57].However, constraints on the inclination angle i were poor.More recently, Kandel and Romani [24] conducted an analysis of heating models for the light curve, leading to improved determinations of the NS masses.In their preferred model, the inclination angle is estimated to be i = 68.7 • ± 2.1 • .With a radial velocity amplitude of K 2 = 641.2± 3.6 and a companion mass of m c = 0.012 ± 0.006, the pulsar's mass is inferred to be m = 2.22 ± 0.10 M ⊙ , above the 95%HPDI derived in the "accuracy-dependent" model.

PSR B1516+02B
This pulsar is located in the globular cluster NGC 5904 and is part of a binary system with a companion that is either a WD or a low-mass main sequence (MS) star.The binary system has a total mass of 2.29 ± 0.17 M ⊙ , leading to a best-fit pulsar mass estimate of m = 2.08 ± 0.19 M ⊙ [11].There is a 90% probability that the pulsar's mass is greater than 1.82 M ⊙ , and a 0.77% probability that the inclination angle is low enough for the pulsar's mass to fall within the range of 1.20-1.44M ⊙ .

PSR J1748-2021B
This pulsar is part of a massive binary system with a total mass of 2.92 ± 0.20 M ⊙ , which was obtained from precise measurements of ω under the assumption of a fully relativistic system [12].The probability that the pulsar's mass falls within the range of 1.20-1.44M ⊙ is only 0.10%, requiring an extremely low orbital inclination.The median mass of the companion star is estimated to be 0.142 with lower and upper 1σ limits at 0.124 and 0.228 M ⊙ , respectively.This range of companion masses suggests that the companion star could be a WD or a non-evolved MS star, which implies a pulsar mass of 2.74 M ⊙ .

Effects on the Maximum Mass
We now proceed to investigate the impact of a uniform distribution for cos i on the posterior distribution of m max .Following the approach of Alsing et al. [32] and Shao et al. [34], we obtained the marginalized posterior distribution shown by the solid line in the left panel of Figure 2, with a maximum of ∼2.2 M ⊙ .On the other hand, the marginalized posterior distribution for the model where all pulsars are assumed to be normally distributed ("accuracy-independent" model) is shown by the solid line in the right panel of Figure 2, with a maximum of around 2.6 M ⊙ .In addition to the maximum values that differ considerably between the two treatments, the shape of the m max distribution changes significantly too, as we can see in Figure 2.
For both models, we verified the impact of changing the masses of PSR B1957+20, PSR J1048+2339, PSR J1555-2908, and PSR J2129-0429 by the values derived from γ-ray observations, listed in Table 1.Since in the analysis of Clark et al. [48] the only value that changed considerably compared to previous estimates was that of PSR B1957+20, we did not expect a significant change in the overall results, as confirmed by the dashed line in both panels of Figure 2.Only the amplitude is mildly affected, but the point estimates for the distributions remain the same.
In the subsequent analysis, we examined the impact of changing the treatment of individual pulsar masses in the "accuracy-dependent" model for the four systems mentioned in the previous section (PSR B1957+20, PSR J1311-3430, PSR B1516+02B, PSR J1748-2021B).In order to do that, we kept the treatment described in Section 3.1 for all systems, with the exception of those mentioned above, which are now also assumed to follow a normal distribution.The result is shown by the dot-dashed line in the left panel of Figure 2. By changing the likelihood treatment for only these particular systems, the marginalized dis-tribution of m max shifts to the right, and the probability of higher m max values is increased.The dot-dashed line in the left panel tends to the solid line in the right panel as we assume more systems to be modeled with Gaussians, and equalizes when all systems are treated as Gaussian.As we can note, the lower individual mass estimates obtained when a flat distribution over i is invoked (and discussed in Section 4.1) are reflected, consequently, in a lower estimation of m max .Marginal posterior distribution of m max sampled from the "accuracy-dependent" model (left panel) and the "accuracy-independent" model (right panel).The solid line represents the analysis assuming the individual pulsars to be sampled from 3 types of likelihoods (left panel), and to be sampled from normal distributions only (right panel).In the dashed curve, we analyze the effect of assuming the masses found from γ-ray observations, listed in Table 1.
As demonstrated, the prior assumption regarding the orbital inclination angle significantly influences the determination of the mass limits, emphasizing the importance of a cautious approach.The high amplitude of values close to 2 M ⊙ in the left panel of Figure 2 is a consequence of lower estimates of individual masses found in this particular approach, leading to a density accumulation around 2 M ⊙ and a significant decrease in the probability for values ≥ 2.4 M ⊙ .
In this scenario, the "accuracy-independent" model is preferred, since the individual sampled masses are consistent with values derived from observations.The ideal scenario, however, would involve the adoption of the "accuracy-dependent" model, with consideration of observational constraints on the inclination angle (i) for each binary system.This meticulous treatment remains a subject for a future work.

Posterior Predictive Check
In Bayesian analysis, a crucial step to validate a model is to assess whether predictive simulated data resemble the observed data.In essence, we aim to understand whether new observations drawn from the posterior distribution would be consistent with the existing observed sample.Any significant disparity might indicate a misfit.This analytical process is known as a posterior predictive check (PPC).One approach to conducting a PPC is to visually compare summaries of real data with those of simulated data.In addition, it can be valuable to quantitatively assess the level of discrepancy.This can be achieved by defining a "test quantity" (T), often chosen as the mean, and calculating a Bayesian p-value.This p-value indicates the probability that the test quantity of simulated data (T sim ) exceeds the observed test quantity (T): Following the detection of GW170817, several studies [31,33,[58][59][60] derived mass limits for NSs based on observational constraints set by the GW event.All these analyses resulted in maximum masses that are significantly lower than the m max we found for the "accuracy-independent" model.
Our goal is to investigate the maximum value of the distribution and, faced with the substantial differences we have identified for this quantity between the two models we treated here and with the values derived from GW data, we sought further investigation to determine which value of m max complies with the available sample of NS masses.Consider that we are analyzing a distribution without specifying that it has to be truncated, the question we want to address is: what would be the value consistent with an upper threshold for this distribution?
For this purpose, we conducted a new MCMC sampling from a non-truncated Gaussian mixture, with a normal likelihood for all pulsars.The only difference to the previous model is the absence of the denominator in Equation (10).In Table 5, we summarize the mean and standard deviation of each model parameter.Using the mean values reported in Table 5, we then generated 10,000 posterior predictive distributions using the software MATHEMATICA [Wolfram Research Inc.].Subsequently, we defined the test quantity, T, as the number of elements in the observed sample with masses exceeding a specific value, referred to as "m max ".This approach enables us to examine the upper limit for the distribution.The p-value is derived from the number of times T sim > T. To illustrate, the top left panel in Figure 3 offers insight into the probability of detecting new objects with m > 2.09 M ⊙ in future observations, which is 47.8%, a high value.A p-value close to zero is consistent with a mass limit, i.e., a low p-value implies a reduced probability of detecting new objects with masses exceeding the associated value.Through this alternative analysis, we confirm that the distribution of galactic NSs supports high m max values, such as the 2.6 M ⊙ we found from the "accuracy-independent" model.
It is important to note that in this analysis all masses were modeled with Gaussians.An inconsistency between results found through the PPC and through the "accuracyindependent" model would also reveal a flaw in the "accuracy-independent" model, which is not the case.This analysis is not intended to be conclusive, but rather complementary.−0.08 M ⊙ .The last panel represents our result for the "accuracy-independent" model.We used the mean value of each referenced work, since they cover the whole range of high masses reasonably well.

Discussions and Conclusions
The old idea of a unique formation channel for NSs and the existence of a canonical mass have long been invalidated, possibly pointing towards a variety of formation mechanisms, including electron capture supernova and accretion-induced collapses.The current debate revolves around whether the mass limit supported by NSs is above or below the range of 2.2-2.3M ⊙ .
Recent analyses of the galactic sample of NSs have found evidence for a sharp cut-off at m max < 2.3 M ⊙ [32,34], even though the number of high-mass estimates in observed NSs continues to grow.The most recent addition to this sample is PSR J0952-0607, with a reported mass of 2.35 ± 0.17 M ⊙ [27].On the other hand, observations from the gravitational wave event, GW170817, have also played a crucial role in constraining the mass limit, with results falling in the interval between 2.1 and 2.3 M ⊙ [31,33,35,[58][59][60], in line with previous galactic analyses.However, the detection of the GW190814 signal, originating from the merger of a 23 M ⊙ BH with an unidentified m = 2.59 +0.08 −0.09 M ⊙ companion, raised a new tension.An analysis of the BBH merging population, performed by the LV team [38], found the GW190814 to be an outlier, i.e., it is potentially associated with an NS-BH merger.In addition, a maximum mass constraint for non-spinning NSs was also derived from the available LV sample of NS-BH events, with a result of 2.7 +0.5 −0.4 M ⊙ [39].More recently, an analysis combining as many astronomical observations as possible led to a preferred threshold value of 2.49-2.52M ⊙ [40], while Fan et al. [61] placed a minimum of 2.25 +0.08 −0.07 M ⊙ to M TOV , two completely independent pieces of evidence, to solve the puzzle of the maximum mass value.
To address this ongoing tension, we performed in this work an analysis of the mass distribution of galactic NSs to investigate the impact of the assumption of a uniform distribution between 0 and 1, associated with the orbital inclination angle i.As previously mentioned, this parameter is responsible for the largest uncertainties in mass measurements and requires careful consideration.Our findings reveal that this uniform assumption results in individual sampled masses that fall below the constraints set from observational modeling (see Section 4.1).As a consequence, the distribution threshold m max also shows a preference for lower values, at ∼2.2 M ⊙ .
In the "accuracy-independent" model, all pulsars are modeled as normal distributions.The mass values reported in Table A1 were derived while taking into account observational constraints at i, which in turn can be model-dependent.This approach is consistent with a scenario favoring the existence of extremely massive NSs (∼2.6 M ⊙ ).A complementary analysis (PPC in Section 3) strengthens the possibility of a high m max for the galactic population, in line with a few conclusions derived from GW analyses, mentioned throughout this work.The discrepancy between the two approaches we treat in this work revealed that it is necessary to be careful when making assumptions over i, since different hypotheses over this quantity can lead to drastically different pulsar masses, consequently influencing the m max determination.An ideal approach would involve implementing the "accuracydependent" model while accounting for the distribution of i for each particular system, constrained from observational data, as performed for the case of black hole masses in [62], where each object was analyzed individually and it was checked that the different Keplerian parameter values employed in each mass estimate did not rely on inconsistent assumptions.A similar treatment for NS masses remains a subject for future work.
Although determining the orbital inclination of a pulsar binary is a challenging task, we showed in this work that the constraints obtained from observations are crucial information about the system and must be treated accordingly.Whether the adopted spectrophotometric model adopted to set these constraints is the ideal or not is another question and needs to be investigated further.The work of Clark et al. [48], if confirmed by future research, can help solve the problem of NS mass accuracy, providing a solution to the NS maximum mass puzzle, which in turn has strong consequences for the supranuclear equation of state.But for now, based on the galactic sample, we cannot rule out the possibility that extremely massive NSs (m > 2.3 M ⊙ ) exists in nature.

5 PDFFigure 1 .
Figure 1.Grey lines represent 1000 posterior samples drawn from the truncated models summarized in Table2(left panel) and Table3(right panel).The blue curve is the posterior mean distribution and the black line is the maximum posterior distribution.

Figure 2 .
Figure 2.Marginal posterior distribution of m max sampled from the "accuracy-dependent" model (left panel) and the "accuracy-independent" model (right panel).The solid line represents the analysis assuming the individual pulsars to be sampled from 3 types of likelihoods (left panel), and to be sampled from normal distributions only (right panel).In the dashed curve, we analyze the effect of assuming the masses found from γ-ray observations, listed in Table1.

Figure 3 .
Figure 3.Posterior predictive check on the two-Gaussian model without truncation.The purpose is to investigate the tail of distributions.High p-values indicate that values higher than the one specified in the label are very common, thus they cannot be pointed to as valid thresholds.The adopted m max from NS-NS mergers are, from left to right and top to bottom, Ai et al.[60] with 2.09+0.11−0.09 M ⊙ ; Shao et al.[35] with 2.13 +0.08 −0.07 M ⊙ ; Rezzola et al.[58] with 2.16 +0.17 −0.15 M ⊙ ; Margalit and Metzger[31] with 2.17 M ⊙ ; Ruiz et al.[59] with 2.16-2.28M ⊙ ; Shibata et al.[33] with 2.3 M ⊙ ; Ai et al.[60] with 2.43 +0.10 −0.08 M ⊙ .The last panel represents our result for the "accuracy-independent" model.We used the mean value of each referenced work, since they cover the whole range of high masses reasonably well.

Table 3 .
Summary of marginal posterior distribution for the "accuracy-independent" model.With the exception of m max , the second column displays the mean value, followed by the respective standard deviation and the 94% highest posterior density interval, which defines the lowest interval that comprises 94% of the probability.For m max , the value shown at column 2 is the mode.

Table 4 .
Individual pulsar mass of systems listed in this table are sampled from Equations (

Table A1 .
Neutron star measurements for 125 binary systems, with 1σ uncertainties.