  We develop a method to compute low-eccentricity initial data of binary
neutron stars required to perform realistic simulations in numerical
relativity. The orbital eccentricity is controlled by adjusting the orbital
angular velocity of a binary and incorporating an approaching relative velocity
of the neutron stars. These modifications improve the solution primarily
through the hydrostatic equilibrium equation for the binary initial data. The
orbital angular velocity and approaching velocity of initial data are updated
iteratively by performing time evolutions over ~3 orbits. We find that the
eccentricity can be reduced by an order of magnitude compared to standard
quasicircular initial data, specifically from ~0.01 to <~0.001, by three
successive iterations for equal-mass binaries leaving ~10 orbits before the
merger.
