  The cross section for e^+e^- \to \eta J/\psi between \sqrt{s}=3.8 GeV/c^2 and
5.3 GeV/c^2 is measured via initial state radiation using 980 fb^{-1} of data
on and around the \Upsilon(nS)(n=1,2,3,4,5) resonances collected with the Belle
detector at KEKB. Two resonant structures at the \psi(4040) and \psi(4160) are
observed in the \eta J/\psi invariant mass distribution. Fitting the mass
spectrum with the coherent sum of two Breit-Wigner functions, one obtains
BR(\psi(4040)\to\eta J/\psi)\cdot\Gamma_{ee}^{\psi(4040)} = (4.8\pm0.9\pm1.4)
eV and BR(\psi(4160)\to\eta J/\psi)\cdot\Gamma_{ee}^{\psi(4160)} =
(4.0\pm0.8\pm1.4) eV for one solution and BR(\psi(4040)\to\eta
J/\psi)\cdot\Gamma_{ee}^{\psi(4040)} = (11.2\pm1.3\pm1.9) eV and
BR(\psi(4160)\to\eta J/\psi)\cdot\Gamma_{ee}^{\psi(4160)} = (13.8\pm1.3\pm2.0)
eV for the other solution, where the first errors are statistical and the
second are systematic. This is the first measurement of this hadronic
transition mode of these two states, and the partial widths to \eta J/\psi are
found to be about 1 MeV. There is no evidence for the Y(4260), Y(4360),
\psi(4415), or Y(4660) in the \eta J/\psi final state, and upper limits of
their production rates in e^+e^- annihilation are determined.
