The recent discovery by LIGO/Virgo of a merging binary having a $\sim 23 M_\odot$ black hole and a $\sim 2.6 M_\odot$ compact companion has triggered a debate regarding the nature of the secondary, which falls into the so-called mass gap. Here we explore some consequences of the assumption that the secondary was a neutron star (NS). We show with concrete examples of heretofore viable equations of state (EOSs) that rapid uniform rotation is certainly adequate to explain the existence of a stable $\sim 2.6 M_\odot$ NS for moderately stiff EOS but may not be adequate for soft EOSs. Moreover, rotation may not even be necessary for sufficiently stiff ones. In particular, several soft EOSs favored by GW170817 and with maximum spherical masses of $\sim 2.1 M_\odot$ \textit{cannot} be invoked to explain this object as a uniformly rotating NS. Absolute upper mass limits on the maximum spherical NS derived from GW170817 immediately suggest that this unknown compact companion can be a slowly or even a nonrotating NS. Stiff EOSs that can achieve this are neither rejected nor favored by GW170817, and they are in accord with the results of NICER.

Recently, we have shown that dynamically stable ergostar solutions (equilibrium neutron stars that contain an ergoregion) with a compressible and causal equation of state exist [A. Tsokaros, M. Ruiz, L. Sun, S. L. Shapiro, and K. Ury\=u, Phys. Rev. Lett. 123, 231103 (2019)]. These stars are hypermassive, differentially rotating, and highly compact. In this work, we make a systematic study of equilibrium models in order to locate the position of ergostars in parameter space. We adopt four equations of state that differ in the matching density of a maximally stiff core. By constructing a large number of models both with uniform and differential rotation of different degrees, we identify the parameters for which ergostars appear. We find that the most favorable conditions for the appearance of dynamically stable ergostars are a significant finite density close to the surface of the star (i.e., similar to self-bound quark stars) and a small degree of differential rotation.

Binary neutron star (NSNS) mergers can be sources of gravitational waves coincident with electromagnetic counterpart emission. To solidify their role as multimessenger sources, we present fully 3D, general relativistic, magnetohydrodynamic simulations of spinning NSNSs initially on quasicircular orbits that merge and undergo delayed collapse to a black hole (BH). The NSNSs consist of two identical stars modeled as $\Gamma=2$ polytropes with spin $\chi_{NS}= 0.36$ aligned along the direction of the total orbital angular momentum $L$. Each star is initially threaded by a dynamical unimportant interior dipole B-field. The field is extended into the exterior where a nearly force-free magnetosphere resembles that of a pulsar. The magnetic dipole moment $\mu$ is either aligned or perpendicular to $L$ and has the same initial magnitude for each orientation. For comparison, we also impose symmetry across the orbital plane in one case where $\mu$ in both stars is aligned along $L$. We find that the lifetime of the transient hypermassive neutron star remnant, the jet launching time, and the ejecta are very sensitive to the B-field orientation. By contrast, the physical properties of the BH + disk remnant, such as the mass and spin of the BH, the accretion rate, and the electromagnetic luminosity, are roughly independent of the initial B-field orientation. In addition, we find imposing symmetry across the orbital plane does not play a significant role in the final outcome of the mergers. Our results show that an incipient jet emerges only when the seed B-field has a sufficiently large-scale poloidal component aligned to $L$. The lifetime [$\Delta t\gtrsim 140(M_{NS}/1.625M_\odot)\rm ms$] and Poynting luminosities [$L_{EM}\simeq 10^{52}$erg/s] of the jet, when it forms, are consistent with typical short gamma-ray bursts, as well as with the Blandford--Znajek mechanism for launching jets.

Can one distinguish a binary black hole undergoing a merger from a binary neutron star if the individual compact companions have masses that fall inside the so-called mass gap of $3-5\ M_\odot$? For neutron stars, achieving such masses typically requires extreme compactness and in this work we present initial data and evolutions of binary neutron stars initially in quasiequilibrium circular orbits having a compactness $C=0.336$. These are the most compact, nonvacuum, quasiequilibrium binary objects that have been constructed and evolved to date, including boson stars. The compactness achieved is only slightly smaller than the maximum possible imposed by causality, $C_{\rm max}=0.355$, which requires the sound speed to be less than the speed of light. By comparing the emitted gravitational waveforms from the late inspiral to merger and postmerger phases between such a binary neutron star vs a binary black hole of the same total mass we identify concrete measurements that serve to distinguish them. With that level of compactness, the binary neutron stars exhibit no tidal disruption up until merger, whereupon a prompt collapse is initiated even before a common core forms. Within the accuracy of our simulations the black hole remnants from both binaries exhibit ringdown radiation that is not distinguishable from a perturbed Kerr spacetime. However, their inspiral leads to phase differences of the order of $\sim 5$ rad over an $\sim 81$ km separation (1.7 orbits) while typical neutron stars exhibit phase differences of $\geq 20$ rad. Although a difference of $\sim 5$ rad can be measured by current gravitational wave laser interferometers (e.g., aLIGO/Virgo), uncertainties in the individual masses and spins will likely prevent distinguishing such compact, massive neutron stars from black holes.

Spin can have significant effects on the electromagnetic transients accompanying binary neutron star mergers. The measurement of spin can provide important information about binary formation channels. In the absence of a strong neutron star spin prior, the degeneracy of spin with other parameters leads to significant uncertainties in their estimation, in particular limiting the power of gravitational waves to place tight constraints on the nuclear equation of state. Thus detailed studies of highly spinning neutron star mergers are essential to understand all aspects of multimessenger observation of such events. We perform a systematic investigation of the impact of neutron star spin---considering dimensionless spin values up to $a_{\rm NS}=0.33$---on the merger of equal mass, quasicircular binary neutron stars using fully general-relativistic simulations. We find that the peak frequency of the post-merger gravitational wave signal is only weakly influenced by the neutron star spin, with cases where the spin is aligned (antialigned) with the orbital angular momentum giving slightly lower (higher) values compared to the irrotational case. We find that the one-arm instability arises in a number of cases, with some dependence on spin. Spin has a pronounced impact on the mass, velocity, and angular distribution of the dynamical ejecta, and the mass of the disk that remains outside the merger remnant. We discuss the implications of these findings on anticipated electromagnetic signals, and on constraints that have been placed on the equation of state based on multimessenger observations of GW170817.

Multi-messenger astrophysics is a fast-growing, interdisciplinary field that combines data, which vary in volume and speed of data processing, from many different instruments that probe the Universe using different cosmic messengers: electromagnetic waves, cosmic rays, gravitational waves and neutrinos. In this Expert Recommendation, we review the key challenges of real-time observations of gravitational wave sources and their electromagnetic and astroparticle counterparts, and make a number of recommendations to maximize their potential for scientific discovery. These recommendations refer to the design of scalable and computationally efficient machine learning algorithms; the cyber-infrastructure to numerically simulate astrophysical sources, and to process and interpret multi-messenger astrophysics data; the management of gravitational wave detections to trigger real-time alerts for electromagnetic and astroparticle follow-ups; a vision to harness future developments of machine learning and cyber-infrastructure resources to cope with the big-data requirements; and the need to build a community of experts to realize the goals of multi-messenger astrophysics.

We construct the first dynamically stable ergostars (equilibrium neutron stars that contain an ergoregion) for a compressible, causal equation of state. We demonstrate their stability by evolving both strict and perturbed equilibrium configurations in full general relativity for over a hundred dynamical timescales ($\gtrsim 30$ rotational periods) and observing their stationary behavior. This stability is in contrast to earlier models which prove radially unstable to collapse. Our solutions are highly differentially rotating hypermassive neutron stars with a corresponding spherical compaction of $C=0.3$. Such ergostars can provide new insights into the geometry of spacetimes around highly compact, rotating objects and on the equation of state at supranuclear densities. Ergostars may form as remnants of extreme binary neutron star mergers and possibly provide another mechanism for powering short gamma-ray bursts.

Following the detection of GW170817 and the accompanying kilonova AT2017gfo, it has become crucial to model and understand the various channels through which mass is ejected in neutron-star binary mergers. We discuss the impact that high stellar spins prior to merger have on the ejection of mass focussing, in particular, on the dynamically ejected mass by performing general-relativistic magnetohydrodynamic simulations employing finite-temperature equations of state and neutrino-cooling effects. Using eight different models with dimensionless spins ranging from $\chi\simeq-0.14$ to $\chi\simeq0.29$ we discuss how the presence of different spins affects the angular distribution and composition of the ejected matter. Most importantly, we find that the dynamical component of the ejected mass can be strongly suppressed in the case of high spins aligned with the orbital angular momentum. In this case, in fact, the merger remnant has an excess angular momentum yielding a more extended and "colder" object, with reduced ability to shed mass dynamically. We discuss how this result impacts the analysis of the recent merger event GW170817 and its kilonova afterglow.

A new code for computing fully general relativistic solutions of strongly magnetized rapidly rotating compact stars is developed as a part of the COCAL (Compact Object CALculator) code. The full set of Einstein's equations, Maxwell's equations and magnetohydrodynamic equations are consistently solved assuming perfect conductivity, stationarity, and axisymmetry, and strongly magnetized solutions associated with mixed poloidal and toroidal components of magnetic fields are successfully obtained in generic (non-circular) spacetimes. We introduce the formulation of the problem and the numerical method in detail, then present examples of extremely magnetized compact star solutions and their convergence tests. It is found that, in extremely magnetized stars, the stellar matter can be expelled from the region of strongest toroidal fields. Hence we conjecture that a toroidal electro-vacuum region may appear inside of the extremely magnetized compact stars, which may seem like the neutron star becoming the strongest toroidal solenoid coil in the universe.

We perform long-term simulations of spinning binary neutron stars, with our highest dimensionless spin being $\chi \sim 0.32$. To assess the importance of spin during the inspiral we vary the spin, and also use two equations of state, one that consists of plain nuclear matter and produces compact stars (SLy), and a hybrid one that contains both nuclear and quark matter and leads to larger stars (ALF2). Using high resolution that has grid spacing $\Delta x\sim 98$ m on the finest refinement level, we find that the effects of spin in the phase evolution of a binary system can be larger than the one that comes from tidal forces. Our calculations demonstrate explicitly that although tidal effects are dominant for small spins ($\lesssim 0.1$), this is no longer true when the spins are larger, but still much smaller than the Keplerian limit.

Events GW170817 and GRB 170817A provide the best confirmation so far that compact binary mergers where at least one of the companions is a neutron star (NS) can be the progenitors of short gamma-ray bursts (sGRBs). An open question for GW170817 remains the values and impact of the initial NS spins. The initial spins could possibly affect the remnant black hole (BH) mass and spin, the remnant disk and the formation and lifetime of a jet and its luminosity. Here we summarize our general relativistic magnetohydrodynamic simulations of spinning, NS binaries undergoing merger and delayed collapse to a BH. The binaries consist of two identical NSs, modeled as $\Gamma=2$ polytropes, in quasicircular orbit, each with spins $\chi_{\rm{NS}}=-0.053,\,0,\,0.24$, or $0.36$. The stars are endowed initially with a dipolar magnetic field extending from the interior into the exterior, as in a radio pulsar. Following merger, the redistribution of angular momentum by magnetic braking and magnetic turbulent viscosity in the hypermassive neutron star (HMNS) remnant, along with the loss of angular momentum due to gravitational radiation, induce the formation of a massive, nearly uniformly rotating inner core surrounded by a magnetized Keplerian disk-like envelope. The HMNS eventually collapses to a BH, with spin $a/M_{\rm BH} \simeq 0.78$ independent of the initial spin of the NSs, surrounded by a magnetized accretion disk. The larger the initial NS spin the heavier the disk. After $\Delta t\sim 3000-4000 M \sim 45-60(M_{\rm NS}/1.625M_\odot)\rm ms$ following merger, a mildly relativistic jet is launched. The lifetime of the jet [$\Delta t\sim 100-140(M_{\rm NS}/1.625M_\odot)\rm ms$] and its outgoing Poynting luminosity [$L_{\rm EM}\sim 10^{51.5\pm 1}\rm erg/s$] are consistent with typical sGRBs, as well as with the Blandford--Znajek mechanism for launching jets and their associated Poynting luminosities.

We present a new initial data formulation to solve the full set of Einstein equations for spacetimes that contain a black hole under general conditions. The method can be used to construct complete initial data for spacetimes (the full metric) that contain a black hole. Contrary to most current studies the formulation requires minimal assumptions. For example, rather than imposing the form of the spatial conformal metric we impose 3 gauge conditions adapted to the coordinates describing the system under consideration. For stationary, axisymmetric spacetimes our method yields Kerr-Schild black holes in vacuum and rotating equilibrium neutron stars. We demonstrate the power of our new method by solving for the first time the whole system of Einstein equations for a nonaxisymmetric, self-gravitating torus in the presence of a black hole. The black hole has dimensionless spin $J_{\rm bh}/M_{\rm bh}^2=0.9918$, a rotation axis tilted at a $30^\circ$ angle with respect to the angular momentum of the disk, and a mass of $\sim 1/5$ of the disk.

Rapidly and differentially rotating compact stars are believed to be formed in binary neutron star merger events, according to both numerical simulations and the multi-messenger observation of GW170817. The lifetime and evolution of such a differentially rotating star, is tightly related to the observations in the post-merger phase. Various studies on the maximum mass of differentially rotating neutron stars have been done in the past, most of which assume the so-called $j$-const law as the rotation profile inside the star and consider only neutron star equations of state. In this paper, we extend the studies to strange star models, as well as to a new rotation profile model. Significant differences are found between differentially rotating strange stars and neutron stars, with both differential rotation laws. A moderate differential rotation rate for neutron stars is found to be too large for strange stars, resulting in a rapid drop in the maximum mass as the differential rotation degree is increased further from $\hat{A}\sim2.0$, where $\hat{A}$ is a parameter characterizing the differential rotation rate for $j$-const law. As a result the maximum mass of a differentially rotating self-bound star drops below the uniformly rotating mass shedding limit for a reasonable degree of differential rotation. The continuous transition to the toroidal sequence is also found to happen at a much smaller differential rotation rate and angular momentum than for neutron stars. In spite of those differences, $\hat{A}$-insensitive relation between the maximum mass for a given angular momentum is still found to hold, even for the new differential rotation law. Astrophysical consequences of these differences and how to distinguish between strange star and neutron star models with future observations are also discussed.

This report provides an overview of recent work that harnesses the Big Data Revolution and Large Scale Computing to address grand computational challenges in Multi-Messenger Astrophysics, with a particular emphasis on real-time discovery campaigns. Acknowledging the transdisciplinary nature of Multi-Messenger Astrophysics, this document has been prepared by members of the physics, astronomy, computer science, data science, software and cyberinfrastructure communities who attended the NSF-, DOE- and NVIDIA-funded "Deep Learning for Multi-Messenger Astrophysics: Real-time Discovery at Scale" workshop, hosted at the National Center for Supercomputing Applications, October 17-19, 2018. Highlights of this report include unanimous agreement that it is critical to accelerate the development and deployment of novel, signal-processing algorithms that use the synergy between artificial intelligence (AI) and high performance computing to maximize the potential for scientific discovery with Multi-Messenger Astrophysics. We discuss key aspects to realize this endeavor, namely (i) the design and exploitation of scalable and computationally efficient AI algorithms for Multi-Messenger Astrophysics; (ii) cyberinfrastructure requirements to numerically simulate astrophysical sources, and to process and interpret Multi-Messenger Astrophysics data; (iii) management of gravitational wave detections and triggers to enable electromagnetic and astro-particle follow-ups; (iv) a vision to harness future developments of machine and deep learning and cyberinfrastructure resources to cope with the scale of discovery in the Big Data Era; (v) and the need to build a community that brings domain experts together with data scientists on equal footing to maximize and accelerate discovery in the nascent field of Multi-Messenger Astrophysics.

For isentropic fluids, dynamical evolution of a binary system conserves the baryonic mass and circulation; therefore, sequences of constant rest mass and constant circulation are of particular importance. In this work, we present the extension of our Compact Object CALculator (\cocal{}) code to compute such quasiequilibria and compare them with the well-known corotating and irrotational sequences, the latter being the simplest, zero-circulation case. The circulation as a measure of the spin for a neutron star in a binary system has the advantage of being exactly calculable since it is a local quantity. To assess the different measures of spin, such as the angular velocity of the star, the quasilocal, dimensionless spin parameter $J/M^2$, or the circulation $\mathcal{C}$, we first compute sequences of single, uniformly rotating stars and describe how the different spin diagnostics are related to each other. The connection to spinning binary systems is accomplished through the concept of circulation and the use of the constant rotational velocity formulation. Finally, we explore a modification of the latter formulation that naturally leads to differentially rotating binary systems.

Black hole-neutron star (BHNS) mergers are one of the most promising targets for multimessenger astronomy. Using general relativistic magnetohydrodynamic simulations of BHNS undergoing merger we showed that a magnetically--driven jet can be launched by the remnant if the NS is endowed with a dipole B field extending from the interior into the exterior as in a radio pulsar. These self-consistent studies considered a BHNS system with mass ratio $q=3:1$, BH spin $a/M_{BH}=0.75$ aligned with the total orbital angular momentum (OAM), and a NS that is irrotational, threaded by an aligned B field, and modeled by an $\Gamma$--law equation of state with $\Gamma=2$. Here, as a crucial step in establishing BHNS systems as viable progenitors of central engines that power short gamma--ray bursts (sGRBs) and thereby solidify their role as multimessenger sources, we survey different BHNS configurations that differ in BH spin ($a/M_{BH} =-0.5,\,0,\,0.5,\,0.75$), in the mass ratio ($q=3:1$ and $q=5:1$), and in the orientation of the B field (aligned and tilted by $90^\circ$ with respect to the OAM). We find that by $\Delta t\sim 3500M-4000M \sim 88(M_{NS}/1.4M_\odot){\rm ms}-100(M_{NS}/1.4M_\odot)\rm ms$ after the peak gravitational wave signal a jet is launched in the cases where the initial BH spin is $a/M_{BH}= 0.5$ or $0.75$. The lifetime of the jets[$\Delta t\sim 0.5(M_{NS}/1.4M_\odot){\rm s-0.7}(M_{NS}/1.4M_\odot)\rm s$] and their Poynting luminosities [$L_{jet}\sim 10^{51\pm 1}\rm erg/s$] are consistent with sGRBs, as well as with the Blandford--Znajek mechanism. By the time we terminate our simulations, we do not observe either an outflow or a large-scale B field collimation in the other configurations we simulate. These results suggest that future multimessenger detections from BHNSs are more likely produced by binaries with highly spinning BH companions and small tilt-angle B fields.

Recent numerical simulations in general relativistic magnetohydrodynamics (GRMHD) provide useful constraints for the interpretation of the GW170817 discovery. Combining the observed data with these simulations leads to a bound on the maximum mass of a cold, spherical neutron star (the TOV limit): ${M_{\rm max}^{\rm sph}}\lesssim 2.74/\beta$, where $\beta$ is the ratio of the maximum mass of a uniformly rotating neutron star (the supramassive limit) over the maximum mass of a nonrotating star. Causality arguments allow $\beta$ to be as high as $1.27$, while most realistic candidate equations of state predict $\beta$ to be closer to $1.2$, yielding ${M_{\rm max}^{\rm sph}}$ in the range $2.16-2.28 M_\odot$. A minimal set of assumptions based on these simulations distinguishes this analysis from previous ones, but leads to a similar estimate. There are caveats, however, and they are enumerated and discussed. The caveats can be removed by further simulations and analysis to firm up the basic argument.

Quasi-equilibrium models of uniformly rotating axisymmetric and triaxial quark stars are computed in general relativistic gravity scenario. The Isenberg-Wilson-Mathews (IWM) formulation is employed and the Compact Object CALculator (COCAL) code is extended to treat rotating stars with finite surface density and new equations of state (EOSs). Besides the MIT bag model for quark matter which is composed of de-confined quarks, we examine a new EOS proposed by Lai and Xu that is based on quark clustering and results in a stiff EOS that can support masses up to $3.3M_\odot$ in the case we considered. We perform convergence tests for our new code to evaluate the effect of finite surface density in the accuracy of our solutions and construct sequences of solutions for both small and high compactness. The onset of secular instability due to viscous dissipation is identified and possible implications are discussed. An estimate of the gravitational wave amplitude and luminosity based on quadrupole formulas is presented and comparison with neutron stars is discussed.

Outcomes of numerical relativity simulations of massive core collapses or binary neutron star mergers with moderate masses suggest formations of rapidly and differentially rotating neutron stars. Subsequent fall back accretion may also amplify the degree of differential rotations. We propose new formulations for modeling differential rotations of those compact stars, and present selected solutions of differentially rotating, stationary, and axisymmetric compact stars in equilibriums. For the cases when rotating stars reach break-up velocities, the maximum masses of such rotating models are obtained.

Targets for ground-based gravitational wave interferometers include continuous, quasiperiodic sources of gravitational radiation, such as isolated, spinning neutron stars. In this work we perform evolution simulations of uniformly rotating, triaxially deformed stars, the compressible analogues in general relativity of incompressible, Newtonian Jacobi ellipsoids. We investigate their stability and gravitational wave emission. We employ five models, both normal and supramassive, and track their evolution with different grid setups and resolutions, as well as with two different evolution codes. We find that all models are dynamically stable and produce a strain that is approximately one-tenth the average value of a merging binary system. We track their secular evolution and find that all our stars evolve towards axisymmetry, maintaining their uniform rotation, kinetic energy, and angular momentum profiles while losing their triaxiality.

We study quasiequilibrium solutions of triaxially deformed rotating compact stars -- a generalization of Jacobi ellipsoids under relativistic gravity and compressible equations of state (EOS). For relatively stiff (piecewise) polytropic EOSs, we find supramassive triaxial solutions whose masses exceed the maximum mass of the spherical solution, but are always lower than those of axisymmetric equilibriums. The difference in the maximum masses of triaxial and axisymmetric solutions depends sensitively on the EOS. If the difference turns out to be only about 10%, it will be strong evidence that the EOS of high density matter becomes substantially softer in the core of neutron stars. This finding opens a novel way to probe phase transitions of high density nuclear matter using detections of gravitational waves from new born neutron stars or magnetars under fallback accretion.

Stationary and axisymmetric solutions of relativistic rotating stars with strong mixed poloidal and toroidal magnetic fields are obtained numerically. Because of the mixed components of the magnetic field, the underlying stationary and axisymmetric spacetimes are no longer circular. These configurations are computed from the full set of the Einstein-Maxwell equations, Maxwell's equations and from first integrals and integrability conditions of the magnetohydrodynamic equilibrium equations. After a brief introduction of the formulation of the problem, we present the first results for highly deformed magnetized rotating compact stars.

We present a new code, named COCAL - Compact Object CALculator, for the computation of equilibriums and quasiequilibrium initial data sets of single or binary compact objects of all kinds. In the cocal code, those solutions are calculated on one or multiple spherical coordinate patches covering the initial hypersurface up to the asymptotic region. The numerical method used to solve field equations written in elliptic form is an adaptation of self-consistent field iterations in which Green's integral formula is computed using multipole expansions and standard finite difference schemes. We extended the method so that it can be used on a computational domain with excised regions for a black hole and a binary companion. Green's functions are constructed for various types of boundary conditions imposed at the surface of the excised regions for black holes. The numerical methods used in cocal are chosen to make the code simpler than any other recent initial data codes, accepting the second order accuracy for the finite difference schemes. We perform convergence tests for time symmetric single black hole data on a single coordinate patch, and binary black hole data on multiple patches. Then, we apply the code to obtain spatially conformally flat binary black hole initial data using boundary conditions including the one based on the existence of equilibrium apparent horizons.

We construct an asymptotic series for a general solution of the Einstein equations near a sudden singularity. The solution is quasi isotropic and contains nine independent arbitrary functions of the space coordinates as required by the structure of the initial value problem.