  Prompted by the recent discovery of pulsed emission from an ultra-luminous
X-ray source, M82 X-2 ("ULX-pulsar"), we perform a two-dimensional
radiation-hydrodynamic simulation of a super-critical accretion flow onto a
neutron star through a narrow accretion column. We set an accretion column with
a cone shape filled with tenuous gas with density of $10^{-4} {\rm g}~ {\rm
cm}^{-3}$ above a neutron star and solve the two dimensional gas motion and
radiative transfer within the column. The side boundaries are set such that
radiation can freely escape, while gas cannot. Since the initial gas layer is
not in a hydrostatic balance, the column gas falls onto the neutron-star
surface, thereby a shock being generated. As a result, the accretion column is
composed of two regions: an upper, nearly free-fall region and a lower settling
region, as was noted by Basko \& Sunyaev (1976). The average accretion rate is
very high; ${\dot M}\sim 10^{2-3} L_{\rm E}/c^2$ (with $L_{\rm E}$ being the
Eddington luminosity), and so radiation energy dominates over gas internal
energy entirely within the column. Despite the high accretion rate, the
radiation flux in the laboratory frame is kept barely below $L_{\rm E}/(4\pi
r^2)$ at a distance $r$ in the settling region so that matter can slowly
accrete. This adjustment is made possible, since large amount of photons
produced via dissipation of kinetic energy of matter can escape through the
side boundaries. The total luminosity can greatly exceed $L_{\rm E}$ by several
orders of magnitude, whereas the apparent luminosity observed from the top of
the column is much less. Due to such highly anisotropic radiation fields,
observed flux should exhibit periodic variations with the rotation period,
provided that the rotation and magnetic axes are misaligned.
