  Numerical simulations of radiative two-temperature hot accretion flows (HAFs)
around Neutron stars (NSs) are performed. We assume that all of the energy
carried by the HAF around a NS will be thermalized and radiated out at the
surface of the NS. The thermal photons will propagate outwards radially and
cool the HAF vis Comptonization. We define $\dot m$ as the mass accretion rate
at the surface of the central object in unit of Eddington accretion rate ($\dot
M_{\rm Edd}=10L_{\rm Edd}/c^2$, with $L_{\rm Edd}$ and $c$ being Eddington
luminosity and speed of light, respectively). When $\dot m$ is lower than $\sim
10^{-4}$, the cooling of the HAF is not important and outflows are very strong.
When $\dot m > \sim 10^{-3}$, cooling becomes important and outflows are
significantly weak. In the range $10^{-4} < \dot m < 10^{-3}$, the HAFs
transients from a strong outflow phase to a very weak outflow phase with
increase of $\dot m$. The properties of the HAF around a NS are also compared
to those of the HAF around a BH. We find that with a similar $\dot m$, the
dynamical properties of the HAF around a NS are quite similar as those of the
HAF around a BH. However, the emitted spectrum of a HAF around a NS can be
quite different from that of a HAF around a BH due to the presence of a thermal
soft X-ray component coming from the surface of the NS.
