What constraints can one pose on the maximum mass of neutron stars from multi-messenger observations?

Abstract: The maximum mass of neutron stars ($M_{\rm TOV}$) plays a crucial role in understanding their equation of state (EoS). Previous studies have used the measurements for the compactness of massive pulsars and the tidal deformability of neutron stars in binary neutron star (BNS) mergers to constrain the EoS and thus the $M_{\rm TOV}$. The discovery of the most massive pulsar, PSR J0952-0607, with a mass $\sim 2.35M_{\odot}$, has provided a valuable lower limit for $M_{\rm TOV}$. Another efficient method to constrain $M_{\rm TOV}$ is by examining the type of central remnant formed after a BNS merger. Gravitational wave (GW) data can provide the total mass of the system, while accompanying electromagnetic signals can help infer the remnant type. In this study, we combine all the previous constraints and utilize the observational facts that about $24\%$ of the short gamma-ray bursts are followed by an X-ray internal plateau, which indicate that roughly this fraction of BNS mergers yield supermassive neutron stars, to perform (Markov Chain) Monte Carlo simulations. These simulations allow us to explore the probability density distribution of $M_{\rm TOV}$ and other parameters related to BNS mergers. Our findings suggest that $M_{\rm TOV}$ is likely around $2.49M_{\odot} - 2.52M_{\odot}$, with an uncertainty range of approximately $-0.16M_{\odot}$, $0.15M_{\odot}$ at $1\sigma$ ($2\sigma$) confidence level. Furthermore, we examine the type of merger remnants in specific events like GW170817 and GW190425 to further constrain $M_{\rm TOV}$ and other relevant parameters, which can help to understand the physical processes involved in BNS mergers.