  The survey of Lyman $\alpha$ emitters (LAEs) with Subaru Hyper Suprime-Cam,
called SILVERRUSH (Ouchi et al.), is producing massive data of LAEs at
$z\gtrsim6$. Here we present LAE simulations to compare the SILVERRUSH data. In
162$^3$ comoving Mpc$^3$ boxes, where numerical radiative transfer calculations
of reionization were performed, LAEs have been modeled with physically
motivated analytic recipes as a function of halo mass. We have examined $2^3$
models depending on the presence or absence of dispersion of halo Ly$\alpha$
emissivity, dispersion of the halo Ly$\alpha$ optical depth, $\tau_\alpha$, and
halo mass dependence of $\tau_\alpha$. The unique free parameter in our model,
a pivot value of $\tau_\alpha$, is calibrated so as to reproduce the $z=5.7$
Ly$\alpha$ luminosity function (LF). We compare our model predictions with
Ly$\alpha$ LFs at $z=6.6$ and $7.3$, LAE angular auto-correlation functions
(ACFs) at $z=5.7$ and $6.6$, and LAE fractions in Lyman break galaxies at
$5<z<7$. The Ly$\alpha$ LFs and ACFs are reproduced by multiple models, but the
LAE fraction turns out to be the most critical test. The dispersion of
$\tau_\alpha$ and the halo mass dependence of $\tau_\alpha$ are essential to
explain all observations reasonably. Therefore, a simple model of one-to-one
correspondence between halo mass and Ly$\alpha$ luminosity with a constant
Ly$\alpha$ escape fraction has been ruled out. Based on our best model, we
present a formula to estimate the intergalactic neutral hydrogen fraction,
$x_{\rm HI}$, from the observed Ly$\alpha$ luminosity density at $z\gtrsim6$.
We finally obtain $x_{\rm HI}=0.5_{-0.3}^{+0.1}$ as a volume-average at
$z=7.3$.
