  We present a new theoretical model for the luminosity functions (LFs) of
Lyman alpha (Lya) emitting galaxies in the framework of hierarchical galaxy
formation. We extend a semi-analytic model of galaxy formation that reproduces
a number of observations for local and high-z galaxies, without changing the
original model parameters but introducing a physically-motivated modelling to
describe the escape fraction of Lya photons from host galaxies (f_esc). Though
a previous study using a hierarchical clustering model simply assumed a
constant and universal value of f_esc, we incorporate two new effects on f_esc:
extinction by interstellar dust and galaxy-scale outflow induced as a star
formation feedback. It is found that the new model nicely reproduces all the
observed Lya LFs of the Lya emitters (LAEs) at different redshifts in z ~ 3-6.
Especially, the rather surprisingly small evolution of the observed LAE Lya LFs
compared with the dark halo mass function is naturally reproduced. Our model
predicts that galaxies with strong outflows and f_esc ~ 1 are dominant in the
observed LFs. This is also consistent with available observations, while the
simple universal f_esc model requires f_esc << 1 not to overproduce the
brightest LAEs. On the other hand, we found that our model significantly
overpredicts LAEs at z > 6, and absorption of Lya photons by neutral hydrogen
in intergalactic medium (IGM) is a reasonable interpretation for the
discrepancy. This indicates that the IGM neutral fraction x_HI rapidly evolves
from x_HI << 1 at z < 6 to a value of order unity at z ~ 6-7, which is broadly
consistent with other observational constraints on the reionization history.
