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.

32 pages, 25 bitmapped ps files, to appear in ApJ supplement