The first inspiral of two neutron stars observed in gravitational waves was remarkably close, allowing the kind of simultaneous gravitational wave and electromagnetic observation that had not been expected for several years. Their merger, followed by a gamma-ray burst and a kilonova, was observed across the spectral bands of electromagnetic telescopes. These GW and electromagnetic observations have led to dramatic advances in understanding short gamma-ray bursts; determining the origin of the heaviest elements; and determining the maximum mass of neutron stars. From the imprint of tides on the gravitational waveforms and from observations of X-ray binaries, one can extract the radius and deformability of inspiraling neutron stars. Together, the radius, maximum mass, and causality constrain the neutron-star equation of state, and future constraints can come from observations of post-merger oscillations. We selectively review these results, filling in some of the physics with derivations and estimates.

At second order in perturbation theory, the $r$-modes of uniformly rotating stars include an axisymmetric part that can be identified with differential rotation of the background star. If one does not include radiation-reaction, the differential rotation is constant in time and has been computed by S\'a. It has a gauge dependence associated with the family of time-independent perturbations that add differential rotation to the unperturbed equilibrium star: For stars with a barotropic equation of state, one can add to the time-independent second-order solution arbitrary differential rotation that is stratified on cylinders (that is a function of distance $\varpi$ to the axis of rotation). We show here that the gravitational radiation-reaction force that drives the $r$-mode instability removes this gauge freedom: The expontially growing differential rotation of the unstable second-order $r$-mode is unique. We derive a general expression for this rotation law for Newtonian models and evaluate it explicitly for slowly rotating models with polytropic equations of state.

At second order in perturbation theory, the unstable r-mode of a rotating star includes growing differential rotation whose form and growth rate are determined by gravitational-radiation reaction. With no magnetic field, the angular velocity of a fluid element grows exponentially until the mode reaches its nonlinear saturation amplitude and remains nonzero after saturation. With a background magnetic field, the differential rotation winds up and amplifies the field, and previous work where large mode amplitudes were considered suggests that the amplification may damp out the instability. A background magnetic field, however, turns the saturated time-independent perturbations corresponding to adding differential rotation into perturbations whose characteristic frequencies are of order the Alfv\'en frequency. As found in previous studies, we argue that magnetic- field growth is sharply limited by the saturation amplitude of an unstable mode. In contrast to previous work, however, we show that if the amplitude is small, i.e., of order 10^(-4), then the limit on the magnetic-field growth is stringent enough to prevent the loss of energy to the magnetic field from damping or significantly altering an unstable r-mode in nascent neutron stars with normal interiors and in cold stars whose interiors are type II superconductors. We show this result first for a toy model, and we then obtain an analogous upper limit on magnetic field growth using a more realistic model of a rotating neutron star. Our analysis depends on the assumption that there are no marginally unstable perturbations, and this may not hold when differential rotation leads to a magnetorotational instability.

We discuss magnetic field enhancement by unstable r-modes (driven by the gravitational radiation reaction force) in rotating stars. In the absence of a magnetic field, gravitational radiation exponentially increases the r-mode amplitude $\alpha$, and accelerates differential rotation (secular motion of fluid elements). For a magnetized star, differential rotation enhances the magnetic field energy. Rezzolla et al. (2000--2001) argued that if the magnetic energy grows faster than the gravitational radiation reaction force pumps energy into the r-modes, then the r-mode instability is suppressed. Chugunov (2015) demonstrated that without gravitational radiation, differential rotation can be treated as a degree of freedom decoupled from the r-modes and controlled by the back reaction of the magnetic field. In particular, the magnetic field windup does not damp r-modes. Here we discuss the effect of the back reaction of the magnetic field on differential rotation of unstable r-modes, and show that it limits the generated magnetic field and the magnetic energy growth rate preventing suppression of the r-mode instability by magnetic windup at low saturation amplitudes, $\alpha \ll 1$, predicted by current models.

A principal goal of gravitational-wave astronomy is to constrain the neutron star equation of state (EOS) by measuring the tidal deformability of neutron stars. The tidally induced departure of the waveform from that of point-particle (or spinless binary black hole (BBH)) increases with the stiffness of the EOS. We show that causality (the requirement that the speed of sound is less than the speed of light for a perfect fluid satisfying a one-parameter equation of state) places an upper bound on tidal deformability as a function of mass. Like the upper mass limit, the limit on deformabity is obtained by using an EOS with $v_{sound} = c$ for high densities and matching to a low density (candidate) EOS at a matching density of order nuclear saturation density. We use these results and those of [B.D. Lackey et al., Phys. Rev. D 89, 043009 (2014)] to estimate the resulting upper limit on the gravitational-wave phase shift of a black hole-neutron star (BHNS) binary relative to a BBH. Even for assumptions weak enough to allow an maximum mass of $4\ M_\odot$ (a match at nuclear saturation density to an unusually stiff low-density candidate EOS), the upper limit on dimensionless tidal deformability is stringent. It leads to a still more stringent estimated upper limit on the maximum tidally induced phase shift prior to merger. We comment in an appendix on the relation between causality, the condition $v_{sound} < c$, and the condition $dp/d\epsilon < 1$ for the effective EOS governing the equilibrium star.

Using an extended set of equations of state and a multiple-group multiple-code collaborative effort to generate waveforms, we improve numerical-relativity-based data-analysis estimates of the measurability of matter effects in neutron-star binaries. We vary two parameters of a parameterized piecewise-polytropic equation of state (EOS) to analyze the measurability of EOS properties, via a parameter {\Lambda} that characterizes the quadrupole deformability of an isolated neutron star. We find that, to within the accuracy of the simulations, the departure of the waveform from point-particle (or spinless double black-hole binary) inspiral increases monotonically with {\Lambda}, and changes in the EOS that did not change {\Lambda} are not measurable. We estimate with two methods the minimal and expected measurability of {\Lambda} in second- and third- generation gravitational-wave detectors. The first estimate, using numerical waveforms alone, shows two EOS which vary in radius by 1.3km are distinguishable in mergers at 100Mpc. The second estimate relies on the construction of hybrid waveforms by matching to post-Newtonian inspiral, and estimates that the same EOS are distinguishable in mergers at 300Mpc. We calculate systematic errors arising from numerical uncertainties and hybrid construction, and we estimate the frequency at which such effects would interfere with template-based searches.

Gravitational wave observations can potentially measure properties of neutron star equations of state by measuring departures from the point-particle limit of the gravitational waveform produced in the late inspiral of a neutron star binary. Numerical simulations of inspiraling neutron star binaries computed for equations of state with varying stiffness are compared. As the stars approach their final plunge and merger, the gravitational wave phase accumulates more rapidly if the neutron stars are more compact. This suggests that gravitational wave observations at frequencies around 1 kHz will be able to measure a compactness parameter and place stringent bounds on possible neutron star equations of state. Advanced laser interferometric gravitational wave observatories will be able to tune their frequency band to optimize sensitivity in the required frequency range to make sensitive measures of the late-inspiral phase of the coalescence.

The late inspiral, merger, and ringdown of a black hole-neutron star (BHNS) system can provide information about the neutron-star equation of state (EOS). Candidate EOSs can be approximated by a parametrized piecewise-polytropic EOS above nuclear density, matched to a fixed low-density EOS; and we report results from a large set of BHNS inspiral simulations that systematically vary two parameters. To within the accuracy of the simulations, we find that, apart from the neutron-star mass, a single physical parameter Lambda, describing its deformability, can be extracted from the late inspiral, merger, and ringdown waveform. This parameter is related to the radius, mass, and l=2 Love number, k_2, of the neutron star by Lambda = 2k_2 R^5/3M_{NS}^5, and it is the same parameter that determines the departure from point-particle dynamics during the early inspiral. Observations of gravitational waves from BHNS inspiral thus restrict the EOS to a surface of constant Lambda in the parameter space, thickened by the measurement error. Using various configurations of a single Advanced LIGO detector, we find that neutron stars are distinguishable from black holes of the same mass and that Lambda^{1/5} or equivalently R can be extracted to 10-40% accuracy from single events for mass ratios of Q=2 and 3 at a distance of 100 Mpc, while with the proposed Einstein Telescope, EOS parameters can be extracted to accuracy an order of magnitude better.

The properties of neutron star matter above nuclear density are not precisely known. Gravitational waves emitted from binary neutron stars during their late stages of inspiral and merger contain imprints of the neutron-star equation of state. Measuring departures from the point-particle limit of the late inspiral waveform allows one to measure properties of the equation of state via gravitational wave observations. This and a companion talk by J. S. Read reports a comparison of numerical waveforms from simulations of inspiraling neutron-star binaries, computed for equations of state with varying stiffness. We calculate the signal strength of the difference between waveforms for various commissioned and proposed interferometric gravitational wave detectors and show that observations at frequencies around 1 kHz will be able to measure a compactness parameter and constrain the possible neutron-star equations of state.

A new method is described for constructing initial data for a binary neutron-star (BNS) system in quasi-equilibrium circular orbit. Two formulations for non-conformally flat data, waveless (WL) and near-zone helically symmetric (NHS), are introduced; in each formulation, the Einstein-Euler system, written in 3+1 form on an asymptotically flat spacelike hypersurface, is exactly solved for all metric components, including the spatially non-conformally flat potentials, and for irrotational flow. A numerical method applicable to both formulations is explained with an emphasis on the imposition of a spatial gauge condition. Results are shown for solution sequences of irrotational BNS with matter approximated by parametrized equations of state that use a few segments of polytropic equations of state. The binding energy and total angular momentum of solution sequences computed within the conformally flat -- Isenberg-Wilson-Mathews (IWM) -- formulation are closer to those of the third post-Newtonian (3PN) two point particles up to the closest orbits, for the more compact stars, whereas sequences resulting from the WL/NHS formulations deviate from the 3PN curve even more for the sequences with larger compactness. We think it likely that this correction reflects an overestimation in the IWM formulation as well as in the 3PN formula, by $\sim 1$ cycle in the gravitational wave phase during the last several orbits. The work suggests that imposing spatial conformal flatness results in an underestimate of the quadrupole deformation of the components of binary neutron-star systems in the last few orbits prior to merger.

For an equation of state in which pressure is a function only of density, the analysis of Newtonian stellar structure is simple in principle if the system is axisymmetric, or consists of a corotating binary. It is then required only to solve two equations: one stating that the "injection energy", $\kappa$, a potential, is constant throughout the stellar fluid, and the other being the integral over the stellar fluid to give the gravitational potential. An iterative solution of these equations generally diverges if $\kappa$ is held fixed, but converges with other choices. We investigate the mathematical reason for this convergence/divergence by starting the iteration from an approximation that is perturbatively different from the actual solution. A cycle of iteration is then treated as a linear "updating" operator, and the properties of the linear operator, especially its spectrum, determine the convergence properties. For simplicity, we confine ourselves to spherically symmetric models in which we analyze updating operators both in the finite dimensional space corresponding to a finite difference representation of the problem, and in the continuum, and we find that the fixed-$\kappa$ operator is self-adjoint and generally has an eigenvalue greater than unity; in the particularly important case of a polytropic equation of state with index greater than unity, we prove that there must be such an eigenvalue. For fixed central density, on the other hand, we find that the updating operator has only a single eigenvector, with zero eigenvalue, and is nilpotent in finite dimension, thereby giving a convergent solution.

We report the results of a first study that uses numerical simulations to estimate the accuracy with which one can use gravitational wave observations of double neutron star inspiral to measure parameters of the neutron-star equation of state. The simulations use the evolution and initial-data codes of Shibata and Uryu to compute the last several orbits and the merger of neutron stars, with matter described by a parametrized equation of state. Previous work suggested the use of an effective cutoff frequency to place constraints on the equation of state. We find, however, that greater accuracy is obtained by measuring departures from the point-particle limit of the gravitational waveform produced during the late inspiral. As the stars approach their final plunge and merger, the gravitational wave phase accumulates more rapidly for smaller values of the neutron star compactness (the ratio of the mass of the neutron star to its radius). We estimate that realistic equations of state will lead to gravitational waveforms that are distinguishable from point particle inspirals at an effective distance (the distance to an optimally oriented and located system that would produce an equivalent waveform amplitude) of 100 Mpc or less. As Lattimer and Prakash observed, neutron-star radius is closely tied to the pressure at density not far above nuclear. Our results suggest that broadband gravitational wave observations at frequencies between 500 and 1000 Hz will constrain this pressure, and we estimate the accuracy with which it can be measured. Related first estimates of radius measurability show that the radius can be determined to an accuracy of ~1 km at 100 Mpc.

We introduce a parameterized high-density equation of state (EOS) in order to systematize the study of constraints placed by astrophysical observations on the nature of neutron-star matter. To obtain useful constraints, the number of parameters should be smaller than the number of neutron-star properties that have been measured or will have been measured in the next several years. And the set must be large enough to accurately approximate the large set of candidate EOSs. We find that a parameterized EOS based on piecewise polytropes with 3 free parameters matches to about 4% rms error an extensive set of candidate EOSs at densities below the central density of 1.4 solar mass stars. Adding observations of more massive stars constrains the higher density part of the EOS and requires an additional parameter. We obtain constraints on the allowed parameter space set by causality and by present and near-future astronomical observations. In particular, we emphasize potentially stringent constraints on the EOS parameter space associated with two measured properties of a single star; and we find that a measurement of the moment of inertia of PSR J0737-3039A can strongly constrain the maximum neutron-star mass. We also present in an appendix a more efficient algorithm than has previously been used for finding points of marginal stability and the maximum angular velocity of stable stars.

Equilibria of binary neutron stars in close circular orbits are computed numerically in a waveless formulation: The full Einstein-relativistic-Euler system is solved on an initial hypersurface to obtain an asymptotically flat form of the 4-metric and an extrinsic curvature whose time derivative vanishes in a comoving frame. Two independent numerical codes are developed, and solution sequences that model inspiraling binary neutron stars during the final several orbits are successfully computed. The binding energy of the system near its final orbit deviates from earlier results of third post-Newtonian and of spatially conformally flat calculations. The new solutions may serve as initial data for merger simulations and as members of quasiequilibrium sequences to generate gravitational wave templates, and may improve estimates of the gravitational-wave cutoff frequency set by the last inspiral orbit.

Two Fokker actions and corresponding equations of motion are obtained for two point particles in a post-Minkowski framework, in which the field of each particle is given by the half-retarded + half-advanced solution to the linearized Einstein equations. The first action is parametrization invariant, the second a generalization of the affinely parametrized quadratic action for a relativistic particle. Expressions for a conserved 4-momentum and angular momentum tensor are obtained in terms of the particles' trajectories in this post-Minkowski approximation. A formal solution to the equations of motion is found for a binary system with circular orbits. For a bound system of this kind, the post-Minkowski solution is a toy model that omits nonlinear terms of relevant post-Newtonian order; and we also obtain a Fokker action that is accurate to first post-Newtonian order, by adding to the post-Minkowski action a term cubic in the particle masses. Curiously, the conserved energy and angular momentum associated with the Fokker action are each finite and well-defined for this bound 2-particle system despite the fact that the total energy and angular momentum of the radiation field diverge. Corresponding solutions and conserved quantities are found for two scalar charges (for electromagnetic charges we exhibit the solution found by Schild). For a broad class of parametrization-invariant Fokker actions and for the affinely parametrized action, binary systems with circular orbits satisfy the relation $dE = \Omega dL$ (a form of the first law of thermodynamics), relating the energy, angular velocity and angular momentum of nearby equilibrium configurations.

We study the inertial modes of slowly rotating, fully relativistic compact stars. The equations that govern perturbations of both barotropic and non-barotropic models are discussed, but we present numerical results only for the barotropic case. For barotropic stars all inertial modes are a hybrid mixture of axial and polar perturbations. We use a spectral method to solve for such modes of various polytropic models. Our main attention is on modes that can be driven unstable by the emission of gravitational waves. Hence, we calculate the gravitational-wave growth timescale for these unstable modes and compare the results to previous estimates obtained in Newtonian gravity (i.e. using post-Newtonian radiation formulas). We find that the inertial modes are slightly stabilized by relativistic effects, but that previous conclusions concerning eg. the unstable r-modes remain essentially unaltered when the problem is studied in full general relativity.

We consider compact binary systems, modeled in general relativity as vacuum or perfect-fluid spacetimes with a helical Killing vector k^\alpha, heuristically, the generator of time-translations in a corotating frame. Systems that are stationary in this sense are not asymptotically flat, but have asymptotic behavior corresponding to equal amounts of ingoing and outgoing radiation. For black-hole binaries, a rigidity theorem implies that the Killing vector lies along the horizon's generators, and from this one can deduce the zeroth law (constant surface gravity of the horizon). Remarkably, although the mass and angular momentum of such a system are not defined, there is an exact first law, relating the change in the asymptotic Noether charge to the changes in the vorticity, baryon mass, and entropy of the fluid, and in the area of black holes. Binary systems with M\Omega small have an approximate asymptopia in which one can write the first law in terms of the asymptotic mass and angular momentum. Asymptotic flatness is precise in two classes of solutions used to model binary systems: spacetimes satisfying the post-Newtonian equations, and solutions to a modified set of field equations that have a spatially conformally flat metric. (The spatial conformal flatness formalism with helical symmetry, however, is consistent with maximal slicing only if replaces the extrinsic curvature in the field equations by an artificially tracefree expression in terms of the shift vector.) For these spacetimes, nearby equilibria whose stars have the same vorticity obey the relation \delta M = \Omega \delta J, from which one can obtain a turning point criterion that governs the stability of orbits.

We present calculations of the spin-down of a neutron star atmosphere due to hydrostatic expansion during a Type I X-ray burst. We show that (i) Cumming and Bildsten overestimated the spin-down of rigidly-rotating atmospheres by a factor of two, and (ii) general relativity has a small (5-10%) effect on the angular momentum conservation law. We rescale our results to different neutron star masses, rotation rates and equations of state, and present some detailed rotational profiles. Comparing with recent observations of large frequency shifts in MXB 1658-298 and 4U 1916-053, we find that the spin-down expected if the atmosphere rotates rigidly is a factor of two to three less than the observed values. If differential rotation is allowed to persist, we find that the upper layers of the atmosphere spin down by an amount comparable to the observed values; however, there is no compelling reason to expect the observed spin frequency to be that of only the outermost layers. We conclude that hydrostatic expansion and angular momentum conservation alone cannot account for the largest frequency shifts observed during Type I bursts.

Several recent surprises appear dramatically to have improved the likelihood that the spin of rapidly rotating, newly formed neutron stars (and, possibly, of old stars spun up by accretion) is limited by a nonaxisymmetric instability driven by gravitational waves. Except for the earliest part of the spin-down, the axial l=m=2 mode (an r-mode) dominates the instability, and the emitted waves may be observable by detectors with the sensitivity of LIGO II. A review of these hopeful results is followed by a discussion of constraints on the instability set by dissipative mechanisms, including viscosity, nonlinear saturation, and energy loss to a magnetic field driven by differential rotation.

We study the r-modes and rotational ``hybrid'' modes of relativistic stars. As in Newtonian gravity, the spectrum of low-frequency rotational modes is highly sensitive to the stellar equation of state. If the star and its perturbations obey the same one-parameter equation of state (as with isentropic stars), there exist {\it no pure r-modes at all} - no modes whose limit, for a star with zero angular velocity, is an axial-parity perturbation. Rotating stars of this kind similarly have no pure g-modes, no modes whose spherical limit is a perturbation with polar parity and vanishing perturbed pressure and density. We compute the post-Newtonian corrections to the $l=m$ r-modes of isentropic and non-isentropic uniform density stars.

A brief review of the stability of rotating relativistic stars is followed by a more detailed discussion of recent work on an instability of r-modes, modes of rotating stars that have axial parity in the slow-rotation limit. These modes may dominate the spin-down of neutron stars that are rapidly rotating at birth, and the gravitational waves they emit may be detectable.

Almost none of the r-modes ordinarily found in rotating stars exist, if the star and its perturbations obey the same one-parameter equation of state; and rotating relativistic stars with one-parameter equations of state have no pure r-modes at all, no modes whose limit, for a star with zero angular velocity, is a perturbation with axial parity. Similarly (as we show here) rotating stars of this kind have no pure g-modes, no modes whose spherical limit is a perturbation with polar parity and vanishing perturbed pressure and density. Where have these modes gone? In spherical stars of this kind, r-modes and g-modes form a degenerate zero-frequency subspace. We find that rotation splits the degeneracy to zeroth order in the star's angular velocity $\Omega$, and the resulting modes are generically hybrids, whose limit as $\Omega\to 0$ is a stationary current with axial and polar parts. Because each mode has definite parity, its axial and polar parts have alternating values of $l$. We show that each mode belongs to one of two classes, axial-led or polar-led, depending on whether the spherical harmonic with lowest value of $l$ that contributes to its velocity field is axial or polar. We numerically compute these modes for slowly rotating polytropes and for Maclaurin spheroids, using a straightforward method that appears to be novel and robust. Timescales for the gravitational-wave driven instability and for viscous damping are computed using assumptions appropriate to neutron stars.

Perturbations of rotating relativistic stars can be classified by their behavior under parity. For axial perturbations (r-modes), initial data with negative canonical energy is found with angular dependence $e^{im\phi}$ for all values of $m\geq 2$ and for arbitrarily slow rotation. This implies instability (or marginal stability) of such perturbations for rotating perfect fluids. This low $m$-instability is strikingly different from the instability to polar perturbations, which sets in first for large values of $m$. The timescale for the axial instability appears, for small angular velocity $\Omega$, to be proportional to a high power of $\Omega$. As in the case of polar modes, viscosity will again presumably enforce stability except for hot, rapidly rotating neutron stars. This work complements Andersson's numerical investigation of axial modes in slowly rotating stars.

We study nonaxisymmetric perturbations of rotating relativistic stars. modeled as perfect-fluid equilibria. Instability to a mode with angular dependence $\exp(im\phi)$ sets in when the frequency of the mode vanishes. The locations of these zero-frequency modes along sequences of rotating stars are computed in the framework of general relativity. We consider models of uniformly rotating stars with polytropic equations of state, finding that the relativistic models are unstable to nonaxisymmetric modes at significantly smaller values of rotation than in the Newtonian limit. Most strikingly, the m=2 bar mode can become unstable even for soft polytropes of index $N \leq 1.3$, while in Newtonian theory it becomes unstable only for stiff polytropes of index $N \leq 0.808$. If rapidly rotating neutron stars are formed by the accretion-induced collapse of white dwarfs, instability associated with these nonaxisymmetric, gravitational-wave driven modes may set an upper limit on neutron-star rotation. Consideration is restricted to perturbations that correspond to polar perturbations of a spherical star. A study of axial perturbations is in progress.

Causality alone suffices to set a lower bound on the period of rotation of relativistic stars as a function of their maximum observed mass. That is, by assuming a one-parameter equation of state (EOS) that satisfies v_sound < c and that allows stars with masses as large as the largest observed neutron-star mass, M_sph^max, we find P[ms] > 0.282 + 0.196 ( M_sph^max/M_odot-1.442). The limit does not assume that the EOS agrees with a known low-density form for ordinary matter, but if one adds that assumption, the minimum period is raised by a few percent. Thus the current minimum period of uniformly rotating stars, set by causality, is 0.28ms (0.29ms for stars with normal crust). The minimizing EOS yields models with a maximally soft exterior supported by a maximally stiff core. An analogous upper limit set by causality on the maximum mass of rotating neutron stars requires a low-density match and the limit depends on the matching density, epsilon_m. We recompute it, obtaining a slightly revised value, M_rot^max \sim 6.1( 2 * 10^14 g/cm^3 epsilon_m )^1/2 M_odot.

We present the first direct comparison of codes based on two different numerical methods for constructing rapidly rotating relativistic stars. A code based on the Komatsu-Eriguchi-Hachisu (KEH) method (Komatsu et al. 1989), written by Stergioulas, is compared to the Butterworth-Ipser code (BI), as modified by Friedman, Ipser and Parker. We compare models obtained by each method and evaluate the accuracy and efficiency of the two codes. The agreement is surprisingly good. A relatively large discrepancy recently reported (Eriguchi et al. 1994) is found to arise from the use of two different versions of the equation of state. We find, for a given equation of state, that equilibrium models with maximum values of mass, baryon mass, and angular momentum are (generically) all distinct and either all unstable to collapse or are all stable. Our implementation of the KEH method will be available as a public domain program for interested users.