  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.

General Relativity and Quantum Cosmology;Astrophysics