We present a new method for the numerical solution of the radiative-transfer equation (RTE) in multidimensional scenarios commonly encountered in computational astrophysics. The method is based on the direct solution of the Boltzmann equation via an extension of the Lattice Boltzmann (LB) methods and allows to model the evolution of the radiation field as it interacts with a background fluid, via absorption, emission, and scattering. As a first application of this method, we restrict our attention to a frequency independent (``grey'') formulation within a special-relativistic framework, which can be employed also for classical computational astrophysics. For a number of standard tests that consider the performance of the method in optically thin, optically thick and intermediate regimes with a static fluid, we show the ability of the LB method to produce accurate and convergent results matching the analytic solutions. We also contrast the LB method with commonly employed moment-based schemes for the solution of the RTE, such as the M1 scheme. In this way, we are able to highlight that the LB method provides the correct solution for both non-trivial free-streaming scenarios and the intermediate optical-depth regime, for which the M1 method either fails or provides inaccurate solutions. When coupling to a dynamical fluid, on the other hand, we present the first self-consistent solution of the RTE with LB methods within a relativistic-hydrodynamic scenario. Finally, we show that besides providing more accurate results in all regimes, the LB method features smaller or comparable computational costs compared to the M1 scheme. We conclude that LB methods represent a competitive and promising avenue to the solution of radiative transport, one of the most common and yet important problems in computational astrophysics.

The recent detection of GW190814 featured the merger of a binary with a primary having a mass of $\sim 23\,M_{\odot}$ and a secondary with a mass of $\sim 2.6\,M_{\odot}$. While the primary was most likely a black hole, the secondary could be interpreted as either the lightest black hole or the most massive neutron star ever observed, but also as the indication of a novel class of exotic compact objects. We here argue that the secondary in GW190814 needs not be an ab-initio black hole nor an exotic object; rather, based on our current understanding of the nuclear-matter equation of state, it can be a rapidly rotating neutron star that collapsed to a rotating black hole at some point before merger. Using universal relations connecting the masses and spins of uniformly rotating neutron stars, we estimate the spin, $0.49 \lesssim \chi \lesssim 0.68$, of the secondary -- a quantity not constrained so far by the detection -- and a novel strict lower bound on the maximum mass, $M_{\rm TOV} > 2.08^{+0.04}_{-0.04}\, \,M_{\odot}$, of nonrotating neutron stars, consistent with recent observations of a very massive pulsar. The new lower bound also remains valid even in the less likely scenario in which the secondary neutron star never collapsed to a black hole.

We provide a systematic description of the steps necessary -- and of the potential pitfalls to be encountered -- when implementing a two-moment scheme within an Implicit-Explicit (IMEX) scheme to include radiative-transfer contributions in numerical simulations of general-relativistic (magneto-)hydrodynamics. We make use of the M1 closure, which provides an exact solution for the optically thin and thick limit, and an interpolation between these limits. Special attention is paid to the efficient solution of the emerging set of implicit conservation equations. In particular, we present an efficient method for solving these equations via the inversion of a $4\times 4$-matrix within an IMEX scheme. While this method relies on a few approximations, it offers a very good compromise between accuracy and computational efficiency. After a large number of tests in special relativity, we couple our new radiation code, \texttt{FRAC}, with the general-relativistic magnetohydrodynamics code \texttt{BHAC} to investigate the radiative Michel solution, namely, the problem of spherical accretion onto a black hole in the presence of a radiative field. By performing the most extensive exploration of the parameter space for this problem, we find that the accretion's efficiency can be expressed in terms of physical quantities such as temperature, $T$, luminosity, $L$, and black-hole mass, $M$, via the expression $\varepsilon=(L/L_{\rm Edd})/(\dot{M}/\dot{M}_{\rm Edd})= 7.41\times 10^{-7}\left(T/10^6\,\mathrm{K}\right)^{0.22} \left(L/L_\odot\right)^{0.48} \left(M/M_\odot\right)^{0.48}$, where $L_{\mathrm{Edd}}$ and $\dot{M}_{\mathrm{Edd}}$ are the Eddington luminosity and accretion rate, respectively. Finally, we also consider the accretion problem away from spherical symmetry, finding that the solution is stable under perturbations in the radiation field.

The first binary neutron-star merger event, GW170817, and its bright electromagnetic counterpart have provided a remarkable amount of information. By contrast, the second event, GW190425, with $M_{\rm tot}=3.4^{+0.3}_{-0.1}\,M_{\odot}$ and the lack of an electromagnetic counterpart, has hardly improved our understanding of neutron-star physics. While GW190425 is compatible with a scenario in which the merger has lead to a prompt collapse to a black hole and little ejected matter to power a counterpart, determining the mass ratio and the effective spin $\tilde{\chi}$ of the binary remains difficult. This is because gravitational waveforms cannot yet well constrain the component spins of the binary. However, since the mass of GW190425 is significantly larger than the maximum mass for nonrotating neutron stars, $M_{_{\rm TOV}}$, the mass ratio $q$ cannot be too small, as the heavier star would not be gravitationally stable. Making use of universal relations and a large number of equations of state, we provide limits in the $(\tilde{\chi},q)$ plane for GW190425, namely: $q_{\rm min} \geq 0.38$ and $ \tilde{\chi}_{\rm max}\leq 0.20$, assuming $M_\mathrm{tot} \simeq 3.4\, M_\odot$. Finally, we show how future observations of high-mass binaries can provide a lower bound on $M_{_{\rm TOV}}$.

With the first detection of gravitational waves from a binary system of neutron stars, GW170817, a new window was opened to study the properties of matter at and above nuclear-saturation density. Reaching densities a few times that of nuclear matter and temperatures up to $100\,\rm{MeV}$, such mergers also represent potential sites for a phase transition (PT) from confined hadronic matter to deconfined quark matter. While the lack of a postmerger signal in GW170817 has prevented us from assessing experimentally this scenario, two theoretical studies have explored the postmerger gravitational-wave signatures of PTs in mergers of binary systems of neutron stars. We here extend and complete the picture by presenting a novel signature of the occurrence of a PT. More specifically, using fully general-relativistic hydrodynamic simulations and employing a suitably constructed equation of state that includes a PT, we present the occurrence of a "delayed PT", i.e. a PT that develops only some time after the merger and produces a metastable object with a quark-matter core, i.e. a hypermassive hybrid star. Because in this scenario, the postmerger signal exhibits two distinct fundamental gravitational-wave frequencies -- before and after the PT -- the associated signature promises to be the strongest and cleanest among those considered so far, and one of the best signatures of the production of quark matter in the present Universe.

Exploiting a very large library of physically plausible equations of state (EOSs) containing more than $10^{7}$ members and yielding more than $10^{9}$ stellar models, we conduct a survey of the impact that a neutron-star radius measurement via electromagnetic observations can have on the EOS of nuclear matter. Such measurements are soon to be expected from the ongoing \textit{Neutron Star Interior Composition Explorer} mission and will complement the constraints on the EOS from gravitational-wave detections. Thanks to the large statistical range of our EOS library, we can obtain a first quantitative estimate of the commonly made assumption that the high-density part of the EOS is best constrained when measuring the radius of the most massive, albeit rare, neutron stars with masses $M\gtrsim2.1\,M_\odot$. At the same time, we find that radius measurements of neutron stars with masses $M\simeq1.7-1.85\,M_\odot$ can provide the strongest constraints on the low-density part of the EOS. Finally, we quantify how radius measurements by future missions can further improve our understanding of the EOS of matter at nuclear densities.

We explore in a parameterized manner a very large range of physically plausible equations of state (EOSs) for compact stars for matter that is either purely hadronic or that exhibits a phase transition. In particular, we produce two classes of EOSs with and without phase transitions, each containing one million EOSs. We then impose constraints on the maximum mass, ($M < 2.16 M_{\odot}$), and on the dimensionless tidal deformability ($\tilde{\Lambda} <800$) deduced from GW170817, together with recent suggestions of lower limits on $\tilde{\Lambda}$. Exploiting more than $10^9$ equilibrium models for each class of EOSs, we produce distribution functions of all the stellar properties and determine, among other quantities, the radius that is statistically most probable for any value of the stellar mass. In this way, we deduce that the radius of a purely hadronic neutron star with a representative mass of $1.4\,M_{\odot}$ is constrained to be $12.00\!<\!R_{1.4}/{\rm km}\!<\!13.45$ at a $2$-$\sigma$ confidence level, with a most likely value of $\bar{R}_{1.4}=12.39\,{\rm km}$; similarly, the smallest dimensionless tidal deformability is $\tilde{\Lambda}_{1.4}\!>\!375$, again at a $2$-$\sigma$ level. On the other hand, because EOSs with a phase transition allow for very compact stars on the so-called `twin-star' branch, small radii are possible with such EOSs although not probable, i.e. $8.53\!<\!R_{1.4}/{\rm km}\!<\!13.74$ and $\bar{R}_{1.4}=13.06\,{\rm km}$ at a $2$-$\sigma$ level, with $\tilde{\Lambda}_{1.4}\!>\!35.5$ at a $3$-$\sigma$ level. Finally, since these EOSs exhibit upper limits on $\tilde{\Lambda}$, the detection of a binary with total mass of $3.4\,M_{\odot}$ and $\tilde{\Lambda}_{1.7}\!>\!461$ can rule out twin-star solutions.

Combining the GW observations of merging systems of binary neutron stars and quasi-universal relations, we set constraints on the maximum mass that can be attained by nonrotating stellar models of neutron stars. More specifically, exploiting the recent observation of the GW event GW 170817 and drawing from basic arguments on kilonova modeling of GRB 170817A, together with the quasi-universal relation between the maximum mass of nonrotating stellar models $M_{\rm TOV}$ and the maximum mass supported through uniform rotation $M_{\rm max}=\left(1.20^{+0.02}_{-0.05}\right) M_{\rm TOV}$ we set limits for the maximum mass to be $ 2.01^{+0.04}_{-0.04}\leq M_{\rm TOV}/M_{\odot}\lesssim 2.16^{+0.17}_{-0.15}$, where the lower limit in this range comes from pulsar observations. Our estimate, which follows a very simple line of arguments and does not rely on the modeling of the electromagnetic signal in terms of numerical simulations, can be further refined as new detections become available. We briefly discuss the impact that our conclusions have on the equation of state of nuclear matter.

The stability properties of rotating relativistic stars against prompt gravitational collapse to a black hole are rather well understood for uniformly rotating models. This is not the case for differentially rotating neutron stars, which are expected to be produced in catastrophic events such as the merger of binary system of neutron stars or the collapse of a massive stellar core. We consider sequences of differentially rotating equilibrium models using the $j$-constant law and by combining them with their dynamical evolution, we show that a sufficient stability criterion for differentially rotating neutron stars exists similar to the one of their uniformly rotating counterparts. Namely: along a sequence of constant angular momentum, a dynamical instability sets in for central rest-mass densities slightly below the one of the equilibrium solution at the turning point. In addition, following Breu & Rezzolla (2016), we show that "quasi-universal" relations can be found when calculating the turning-point mass. In turn, this allows us to compute the maximum mass allowed by differential rotation, $M_{\rm max,dr}$, in terms of the maximum mass of the nonrotating configuration, $M_{_{\rm TOV}}$, finding that $M_{\rm max, dr} \simeq \left(1.54 \pm 0.05\right) M_{_{\rm TOV}}$ for all the equations of state we have considered.