All the neutron star ( NS ) atmosphere models published so far have been calculated in the “ cold plasma approximation ” , which neglects the relativistic effects in the radiative processes , such as cyclotron emission/absorption at harmonics of cyclotron frequency . Here we present new NS atmosphere models which include such effects . We calculate a set of models for effective temperatures T _ { eff } = 1 –3 MK and magnetic fields B \sim 10 ^ { 10 } – 10 ^ { 11 } G , typical for the so-called central compact objects ( CCOs ) in supernova remnants , for which the electron cyclotron energy E _ { c,e } and its first harmonics are in the observable soft X-ray range . Although the relativistic parameters , such as kT _ { eff } / m _ { e } c ^ { 2 } and E _ { c,e } / m _ { e } c ^ { 2 } , are very small for CCOs , the relativistic effects substantially change the emergent spectra at the cyclotron resonances , E \approx sE _ { c,e } ( s = 1 , 2 , \ldots ) . Although the cyclotron absorption features can form in a cold plasma due to the quantum oscillations of the free-free opacity , the shape and depth of these features change substantially if the relativistic effects are included . In particular , the features acquire deep Doppler cores , in which the angular distribution of the emergent intensity is quite different from that in the cold plasma approximation . The relative contributions of the Doppler cores to the equivalent widths of the features grow with increasing the quantization parameter b _ { eff } \equiv E _ { c,e } / kT _ { eff } and harmonic number s . The total equivalent widths of the features can reach \sim 150–250 eV ; they increase with growing b _ { eff } and are smaller for higher harmonics .