Neutron star interiors provide the opportunity to probe properties of cold dense matter in the QCD phase diagram. Utilizing models of dense matter in accord with nuclear systematics at nuclear densities, we investigate the compatibility of deconfined quark cores with current observational constraints on the maximum mass and tidal deformability of neutron stars. We explore various methods of implementing the hadron-to-quark phase transition, specifically, first-order transitions with sharp (Maxwell construction) and soft (Gibbs construction) interfaces, and smooth crossover transitions. We find that within the models we apply, hadronic matter has to be stiff for a first-order phase transition and soft for a crossover transition. In both scenarios and for the equations of state we employed, quarks appear at the center of pre-merger neutron stars in the mass range $\approx 1.0-1.6\,{\rm M}_{\odot}$, with a squared speed of sound $c^2_{\rm QM}\gtrsim 0.4$ characteristic of strong repulsive interactions required to support the recently discovered neutron star masses $\geq 2\,{\rm M}_{\odot}$. We also identify equations of state and phase transition scenarios that are consistent with the bounds placed on tidal deformations of neutron stars in the recent binary merger event GW170817. We emphasize that distinguishing hybrid stars with quark cores from normal hadronic stars is very difficult from the knowledge of masses and radii alone, unless drastic sharp transitions induce distinctive disconnected hybrid branches in the mass-radius relation.

In simulations of binary neutron star mergers, the dense matter equation of state (EOS) is required over wide ranges of density and temperature as well as under conditions in which neutrinos are trapped, and the effects of magnetic fields and rotation prevail. Here we assess the status of dense matter theory and point out the successes and limitations of approaches currently in use. A comparative study of the excluded volume (EV) and virial approaches for the $np\alpha$ system using the equation of state of Akmal, Pandharipande and Ravenhall for interacting nucleons is presented in the sub-nuclear density regime. Owing to the excluded volume of the $\alpha$-particles, their mass fraction vanishes in the EV approach below the baryon density 0.1 fm$^{-3}$, whereas it continues to rise due to the predominantly attractive interactions in the virial approach. The EV approach of Lattimer et al. is extended here to include clusters of light nuclei such as d, $^3$H and $^3$He in addition to $\alpha$-particles. Results of the relevant state variables from this development are presented and enable comparisons with related but slightly different approaches in the literature. We also comment on some of the sweet and sour aspects of the supra-nuclear EOS. The extent to which the neutron star gravitational and baryon masses vary due to thermal effects, neutrino trapping, magnetic fields and rotation are summarized from earlier studies in which the effects from each of these sources were considered separately. Increases of about $20\% (\gtrsim 50\%)$ occur for rigid (differential) rotation with comparable increases occurring in the presence of magnetic fields only for fields in excess of $10^{18}$ Gauss. Comparatively smaller changes occur due to thermal effects and neutrino trapping. Some future studies to gain further insight into the outcome of dynamical simulations are suggested.