  Magnetic neutron stars (NSs) often exhibit a cyclotron resonant scattering
feature (CRSF) in their X-ray spectra, but the site of the CRSF formation is
still an open puzzle. A promising candidate for high-luminosity sources has
always been the radiative shock (RS) in the accretion column. Yet, no
quantitative calculations of spectral formation at the RS have been performed
so far. Here we explore the scenario where the shock is the site of the CRSF
formation. We study spectral formation at the RS and the emergent spectral
shape across a wide range of the parameter space. We developed a Monte Carlo
code to conduct radiation transfer simulations at the RS, adopting a fully
relativistic scheme for the interaction between radiation and electrons. We
properly treated bulk-motion Comptonization in the pre-shock region, thermal
Comptonization in the post-shock region, and resonant Compton scattering in
both regions. We calculated the angle- and energy-dependent emergent X-ray
spectrum from the RS, focusing on both the CRSF and the X-ray continuum, under
diverse conditions. We find that a power law, hard X-ray continuum and a CRSF
are naturally produced by the first-order Fermi energization as the photons
criss-cross the shock. The CRSF shape depends mainly on the transverse optical
depth and the post-shock temperature. We show that the CRSF energy centroid is
shifted by ~(20-30)% to lower energies compared to the classical cyclotron
energy, due to the Doppler boosting between the shock frame and the bulk-motion
frame. We demonstrate that a "bump" feature arises in the right wing of the
CRSF due to upscattering of photons by the accreting plasma and extends to
higher energies for larger optical depths and post-shock temperatures. The
implications of the Doppler effect on the centroid of the emergent CRSF must be
considered if an accurate determination of the magnetic field strength is
desired.
