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.

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.

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.

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.

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 report work towards a relativistic formulation for modeling strongly magnetized neutron stars, rotating or in a close circular orbit around another neutron star or black hole, under the approximations of helical symmetry and ideal MHD. The quasi-stationary evolution is governed by the first law of thermodynamics for helically symmetric systems, which is generalized to include magnetic fields. The formulation involves an iterative scheme for solving the Einstein-Maxwell and relativistic MHD-Euler equations numerically. The resulting configurations for binary systems could be used as self-consistent initial data for studying their inspiral and merger.

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.

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.

A fully geometrical treatment of general relativistic magnetohydrodynamics (GRMHD) is developed under the hypotheses of perfect conductivity, stationarity and axisymmetry. The spacetime is not assumed to be circular, which allows for greater generality than the Kerr-type spacetimes usually considered in GRMHD. Expressing the electromagnetic field tensor solely in terms of three scalar fields related to the spacetime symmetries, we generalize previously obtained results in various directions. In particular, we present the first relativistic version of the Soloviev transfield equation, subcases of which lead to fully covariant versions of the Grad-Shafranov equation and of the Stokes equation in the hydrodynamical limit. We have also derived, as another subcase of the relativistic Soloviev equation, the equation governing magnetohydrodynamical equilibria with purely toroidal magnetic fields in stationary and axisymmetric spacetimes.

Binary systems of compact objects with electromagnetic field are modeled by helically symmetric Einstein-Maxwell spacetimes with charged and magnetized perfect fluids. Previously derived thermodynamic laws for helically-symmetric perfect-fluid spacetimes are extended to include the electromagnetic fields, and electric currents and charges; the first law is written as a relation between the change in the asymptotic Noether charge $\dl Q$ and the changes in the area and electric charge of black holes, and in the vorticity, baryon rest mass, entropy, charge and magnetic flux of the magnetized fluid. Using the conservation laws of the circulation of magnetized flow found by Bekenstein and Oron for the ideal magnetohydrodynamic (MHD) fluid, and also for the flow with zero conducting current, we show that, for nearby equilibria that conserve the quantities mentioned above, the relation $\dl Q=0$ is satisfied. We also discuss a formulation for computing numerical solutions of magnetized binary compact objects in equilibrium with emphasis on a first integral of the ideal MHD-Euler equation.

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.

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.

Quasi-equilibrium models of rapidly rotating triaxially deformed stars are computed in general relativistic gravity, assuming a conformally flat spatial geometry (Isenberg-Wilson-Mathews formulation) and a polytropic equation of state. Highly deformed solutions are calculated on the initial slice covered by spherical coordinate grids, centered at the source, in all angular directions up to a large truncation radius. Constant rest mass sequences are calculated from nearly axisymmetric to maximally deformed triaxial configurations. Selected parameters are to model (proto-) neutron stars; the compactness is M/R = 0.001, 0.1, 0.14, 0.2 for polytropic index n = 0.3 and M/R = 0.001, 0.1, 0.12, 0.14 for n = 0.5. We confirmed that the triaxial solutions exist for these parameters as in the case of Newtonian polytropes. However, it is also found that the triaxial sequences become shorter for higher compactness, and those may disappear at a certain large compactness for the n = 0.5 case. In the scenario of the contraction of proto-neutron stars being subject to strong viscosity and rapid cooling, it is plausible that, once the viscosity driven secular instability sets in during the contraction, the proto-neutron stars are always maximally deformed triaxial configurations, as long as the compactness and the equation of state parameters allow such triaxial sequences. Detection of gravitational waves from such sources may be used as another probe for the nuclear equation of state.

We present our latest results for simulation for merger of black hole (BH)-neutron star (NS) binaries in full general relativity which is performed preparing a quasicircular state as initial condition. The BH is modeled by a moving puncture with no spin and the NS by the $\Gamma$-law equation of state with $\Gamma=2$ and corotating velocity field as a first step. The mass of the BH is chosen to be $\approx 3.2 M_{\odot}$ or $4.0M_{\odot}$, and the rest-mass of the NS $\approx 1.4 M_{\odot}$ with relatively large radius of the NS $\approx 13$--14 km. The NS is tidally disrupted near the innermost stable orbit but $\sim 80$--90% of the material is swallowed into the BH and resulting disk mass is not very large as $\sim 0.3M_{\odot}$ even for small BH mass $\sim 3.2M_{\odot}$. The result indicates that the system of a BH and a massive disk of $\sim M_{\odot}$ is not formed from nonspinning BH-NS binaries irrespective of BH mass, although a disk of mass $\sim 0.1M_{\odot}$ is a possible outcome for this relatively small BH mass range as $\sim 3$--4$M_{\odot}$. Our results indicate that the merger of low-mass BH and NS may form a central engine of short-gamma-ray bursts.

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 present numerical results of three-dimensional simulations for the merger of binary neutron stars (BNSs) in full general relativity. Hybrid equations of state (EOSs) are adopted to mimic realistic nuclear EOSs. In this approach, we divide the EOSs into two parts, i.e., the thermal part and the cold part. For the cold part, we assign a fitting formula for realistic EOSs of cold nuclear matter slightly modifying the formula developed by Haensel and Potekhin. We adopt the SLy and FPS EOSs for which the maximum allowed ADM mass of cold and spherical neutron stars (NSs) is ~ 2.04Mo and 1.80Mo, respectively. Simulations are performed for BNSs of the total ADM mass in the range between 2.4Mo and 2.8Mo with the rest-mass ratio Q_M to be in the range 0.9 < Q_M < 1. It is found that if the total ADM mass of the system is larger than a threshold M_{thr}, a black hole (BH) is promptly formed in the merger irrespective of the mass ratios. In the other case, the outcome is a hypermassive NS of a large ellipticity, which results from the large adiabatic index of the realistic EOSs adopted. The value of M_{thr} depends on the EOS: ~ 2.7Mo and ~ 2.5Mo for the SLy and FPS EOSs, respectively. Gravitational waves are computed in terms of a gauge-invariant wave extraction technique. In the formation of the hypermassive NS, quasiperiodic gravitational waves of a large amplitude and of frequency between 3 and 4 kHz are emitted. The estimated emission time scale is < 100 ms, after which the hypermassive NS collapses to a BH. Because of the long emission time, the effective amplitude may be large enough to be detected by advanced laser interferometric gravitational wave detectors if the distance to the source is smaller than ~ 100 Mpc.

We present results of three dimensional numerical simulations of the merger of unequal-mass binary neutron stars in full general relativity. A $\Gamma$-law equation of state $P=(\Gamma-1)\rho\epsilon$ is adopted, where $P$, $\rho$, $\varep$, and $\Gamma$ are the pressure, rest mass density, specific internal energy, and the adiabatic constant, respectively. We take $\Gamma=2$ and the baryon rest-mass ratio $Q_M$ to be in the range 0.85--1. The typical grid size is $(633,633,317)$ for $(x,y,z)$ . We improve several implementations since the latest work. In the present code, the radiation reaction of gravitational waves is taken into account with a good accuracy. This fact enables us to follow the coalescence all the way from the late inspiral phase through the merger phase for which the transition is triggered by the radiation reaction. It is found that if the total rest-mass of the system is more than $\sim 1.7$ times of the maximum allowed rest-mass of spherical neutron stars, a black hole is formed after the merger irrespective of the mass ratios. The gravitational waveforms and outcomes in the merger of unequal-mass binaries are compared with those in equal-mass binaries. It is found that the disk mass around the so formed black holes increases with decreasing rest-mass ratios and decreases with increasing compactness of neutron stars. The merger process and the gravitational waveforms also depend strongly on the rest-mass ratios even for the range $Q_M= 0.85$--1.

Binary neutron star mergers are expected to be one of the most promising source of gravitational waves (GW) for the network of laser interferometric and bar detectors becoming operational in the next few years. The merger wave signal is expected to be sensitive to the interior structure of the neutron star (NS). The structure of high density phases of matter is under current experimental investigation in heavy-ion collisions. We investigate the dependence of the merger process and its GW signal on the presence of quarks in these phases by performing numerical simulations, where the smoothed particle hydrodynamics (SPH) method and the conformally flat approximation for the 3-geometry in general relativistic gravity are implemented.

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 model the adiabatic inspiral of relativistic binary neutron stars in a quasi-equilibrium (QE) approximation, and compute the gravitational wavetrain from the late phase of the inspiral. We compare corotational and irrotational sequences and find a significant difference in the inspiral rate, which is almost entirely caused by differences in the binding energy. We also compare our results with those of a point-mass post-Newtonian calculation. We illustrate how the late inspiral wavetrain computed with our QE numerical scheme can be matched to the subsequent plunge and merger waveform calculated with a fully relativistic hydrodynamics code.

Gravitational wave emission is considered to be the driving force for the evolution of short-period cataclysmic binary stars, making them a potential test for the validity of General Relativity. In spite of continuous refinements of the physical description, a 10% mismatch exists between the theoretical minimum period ($P_{\rm turn} \simeq 70$ min) and the short-period cut-off ($P_{\rm min} \simeq 80$ min) observed in the period distribution for cataclysmic variable binaries. A possible explanation for this mismatch was associated with the use of the Roche model. We here present a systematic comparison between self-consistent, numerically constructed sequences of hydrostatic models of binary stars and Roche models of semi-detached binaries. On the basis of our approach, we also derive a value for the minimum period of cataclysmic variable binaries. The results obtained through the comparison indicate that the Roche model is indeed very good, with deviations from the numerical solution which are of a few percent at most. Our results therefore suggest that additional sources of angular momentum loss or alternative explanations need to be considered in order to justify the mismatch.

We investigate the equilibrium and stability of supermassive stars of mass $M \agt 10^5M_{\odot}$ in binary systems. We find that corotating binaries are secularly unstable for close, circular orbits with $r \alt 4R(M/10^6M_{\odot})^{1/6}$ where $r$ is the orbital separation and $R$ the stellar radius. We also show that corotation cannot be achieved for distant orbits with $r \agt 12 R (M/10^6M_{\odot})^{-11/24}$, since the timescale for viscous angular momentum transfer associated with tidal torques is longer than the evolution timescale due to emission of thermal radiation. These facts suggest that the allowed mass range and orbital separation for corotating supermassive binary stars is severely restricted. In particular, for supermassive binary stars of large mass $M \agt 6\times 10^6M_{\odot}$, corotation cannot be achieved, as viscosity is not adequate to mediate the transfer between orbital and spin angular momentum. One possible outcome for binary supermassive stars is the onset of quasi-radial, relativistic instability which drives each star to collapse prior to merger: We discuss alternative outcomes of collapse and possible spin states of the resulting black holes. We estimate the frequency and amplitude of gravitational waves emitted during several inspiral and collapse scenarios.

We perform 3D numerical simulations for merger of equal mass binary neutron stars in full general relativity preparing irrotational binary neutron stars in a quasiequilibrium state as initial conditions. Simulations have been carried out for a wide range of stiffness of equations of state and compactness of neutron stars, paying particular attention to the final products and gravitational waves. We take a fixed uniform grid in Cartesian coordinates with typical grid size $(293,293,147)$ in $(x,y,z)$ assuming a plane symmetry with respect to the equatorial plane. A result of one new large-scale simulation performed with grid size $(505,505,253)$ is also presented. We find that the final product depends sensitively on the initial compactness of the neutron stars: In a merger between sufficiently compact neutron stars, a black hole is formed in a dynamical timescale. As the compactness is decreased, the formation timescale becomes longer and longer. For less compact cases, a differentially rotating massive neutron star is formed instead of a black hole. In the case of black hole formation, the disk mass around the black hole appears to be smaller than 1% of the total rest mass. It is also indicated that waveforms of high-frequency gravitational waves after merger depend strongly on the compactness of neutron stars.

We investigate close binary neutron stars in quasiequilibrium states in a general relativistic framework. We assume conformal flatness for the spatial metric and irrotational velocity field for the neutron stars. We adopt the polytropic equation of state. The computation is performed for the polytropic index n(=0.5, 0.66667, 0.8, 1, 1.25), and compactness of neutron stars M/R(=0.03 - 0.3). Results of this paper are as follows. (i) The sequences of the irrotational binary are always terminated at an innermost orbit where a cusp (inner Lagrange point)appears at the inner edges of the stellar surface. The binaries with cusps are found to be dynamically unstable for n=0.5 and stable for n > 0.8 irrespective of M/R < 0.2. For n=0.66667, the stability changes depending on M/R. (ii) The gravitational wave frequency at the innermost orbit turns out to be between 800 and 1500 Hz for realistic compactness 0.14 < M/R < 0.2. (iii) The ISCO for n=0.5 appears to be determined by a hydrodynamic instability for M/R < 0.2. We derive fitting formulae for the relation between the orbital angular velocity at the ISCO and the compactness to clarify it. (iv) The maximum density of neutron stars in binary systems slightly decreases with decreasing the orbital separation and hence they are stable against individual radial collapse during the inspiral. (v) q = J_tot/M_ADM^2 at the innermost orbits is always less than unity for M/R > 0.13 irrespective of n, which indicates that the realistic binary neutron stars satisfy a necessary condition (q<1) for formation of a black hole before the merger. (vi) The specific angular momentum of any mass element in irrotational binary neutron stars at the innermost orbit appears to be too small to form a disk around black holes formed after the merger.

We propose a new numerical method to compute quasi-equilibrium sequences of general relativistic irrotational binary neutron star systems. It is a good approximation to assume that (1) the binary star system is irrotational, i.e. the vorticity of the flow field inside component stars vanishes everywhere (irrotational flow), and (2) the binary star system is in quasi-equilibrium, for an inspiraling binary neutron star system just before the coalescence as a result of gravitational wave emission. We can introduce the velocity potential for such an irrotational flow field, which satisfies an elliptic partial differential equation (PDE) with a Neumann type boundary condition at the stellar surface. For a treatment of general relativistic gravity, we use the Wilson--Mathews formulation, which assumes conformal flatness for spatial components of metric. In this formulation, the basic equations are expressed by a system of elliptic PDEs. We have developed a method to solve these PDEs with appropriate boundary conditions. The method is based on the established prescription for computing equilibrium states of rapidly rotating axisymmetric neutron stars or Newtonian binary systems. We have checked the reliability of our new code by comparing our results with those of other computations available. We have also performed several convergence tests. By using this code, we have obtained quasi-equilibrium sequences of irrotational binary star systems with strong gravity as models for final states of real evolution of binary neutron star systems just before coalescence. Analysis of our quasi-equilibrium sequences of binary star systems shows that the systems may not suffer from dynamical instability of the orbital motion and that the maximum density does not increase as the binary separation decreases.

We have performed 3D numerical simulations for merger of equal mass binary neutron stars in full general relativity. We adopt a $\Gamma$-law equation of state in the form $P=(\Gamma-1)\rho\epsilon$ where P, $\rho$, $\varep$ and $\Gamma$ are the pressure, rest mass density, specific internal energy, and the adiabatic constant with $\Gamma=2$. As initial conditions, we adopt models of corotational and irrotational binary neutron stars in a quasi-equilibrium state which are obtained using the conformal flatness approximation for the three geometry as well as an assumption that a helicoidal Killing vector exists. In this paper, we pay particular attention to the final product of the coalescence. We find that the final product depends sensitively on the initial compactness parameter of the neutron stars : In a merger between sufficiently compact neutron stars, a black hole is formed in a dynamical timescale. As the compactness is decreased, the formation timescale becomes longer and longer. It is also found that a differentially rotating massive neutron star is formed instead of a black hole for less compact binary cases, in which the rest mass of each star is less than 70-80% of the maximum allowed mass of a spherical star. In the case of black hole formation, we roughly evaluate the mass of the disk around the black hole. For the merger of corotational binaries, a disk of mass $\sim 0.05-0.1M_*$ may be formed, where M_* is the total rest mass of the system. On the other hand, for the merger of irrotational binaries, the disk mass appears to be very small : < 0.01M_*.

Circularly orbiting black hole-gaseous star close binary systems are examined by using numerically exact stationary configurations in the framework of Newtonian gravity. We have chosen a polytropic star for the fluid component of the binary system and considered two ideal situations: 1) a synchronously rotating star and 2) an irrotationally rotating star. They correspond to a rotating star under the influence of viscosity and to that in the inviscid limit, respectively. By analyzing the stationary sequences of binary systems with small separations, we can discuss the final stages of black hole-gaseous star close binary systems. Our computational results show that the binary systems reach the Roche(-Riemann) limit states or the Roche lobe filling states without suffering from hydrodynamical instability due to tidal force for a large parameter range of the mass ratio and the polytropic index. It is very likely that such stable Roche(-Riemann) limits or Roche lobe filling states survive even under the general relativistic effect. Therefore, at the final stage of the evolution which is caused by the emission of gravitational waves, the Roche overflow will occur instead of merging of a black hole and a star.

We propose a new numerical method to calculate irrotational binary systems composed of compressible gaseous stars in Newtonian gravity. Assuming irrotationality, i.e. vanishing of the vorticity vector everywhere in the star in the inertial frame, we can introduce the velocity potential for the flow field. Using this velocity potential we can derive a set of basic equations for stationary states which consist of (i) the generalized Bernoulli equation, (ii) the Poisson equation for the Newtonian gravitational potential and (iii) the equation for the velocity potential with the Neumann type boundary condition. We succeeded in developing a new code to compute numerically exact solutions to these equations for the first time. Such irrotational configurations of binary systems are appropriate models for realistic neutron star binaries composed of inviscid gases, just prior to coalescence of two stars caused by emission of gravitational waves. Accuracies of our numerical solutions are so high that we can compute reliable models for fully deformed final stationary configurations and hence determine the inner most stable circular orbit of binary neutron star systems under the approximations of weak gravity and inviscid limit.

In this paper we have examined numerically exact configurations of close binary systems composed of incompressible fluids with internal flows. Component stars of binary systems are assumed to be circularly orbiting each other but rotating nonsynchronously with the orbital motion, i.e. stars in binary systems have steady motions seen from a rotational frame of reference. We have computed several equilibrium sequences by taking fully into account the tidal effect of Newtonian gravity without approximation. We consider two binary systems consisting of either 1) a point mass and a fluid star or 2) a fluid star and a fluid star. Each of them corresponds to generalization of the Roche--Riemann or the Darwin--Riemann problem, respectively. Our code can be applied to various types of incompressible binary systems with various mass ratios and spin as long as the vorticity is constant. We compare these equilibrium sequences of binaries with approximate solutions which were studied extensively by Lai, Rasio and Shapiro (LRS) as models for neutron star--neutron star (NS--NS) binary systems or black hole--neutron star (BH--NS) binary systems. Our results coincide qualitatively with those of LRS but are different from theirs for configurations with small separations. For these binary systems, our sequences show that dynamical or secular instability sets in as the separation decreases. The quantitative errors of the ellipsoidal approximation amount to 2\sim 10% for configurations near the instability point. Compared to the results of LRS, the separation of the stars at the point where the instability sets in is larger and correspondingly the orbital frequency at the critical state is smaller for our models.

We have succeeded in obtaining exact configurations of irrotational binary systems for compressible (polytropic) equations of state. Our models correspond to binaries of equal mass neutron star systems in the viscosity free limit. By using the obtained sequences of stationary states, the evolution of binary systems of irrotational neutron stars due to gravitational wave emission has been examined. For inviscid binary systems the spin angular velocity of each component in a detached phase is smaller than the orbital angular velocity at a contact phase. Thus the irrotational approximation during evolution of binary neutron stars due to gravitational wave emission can be justified. Our computational results show that the binary will never reach a dynamically unstable state before a contact phase even for rather stiff polytropes with the index N > 0.7, as the separation of two components decreases due to gravitational wave emission. This conclusion is quantitatively different from that of Lai, Rasio & Shapiro who employed approximate solutions for polytropic binary systems.