  We present two observations of the high-mass X-ray binary GX 301-2 with
NuSTAR, taken at different orbital phases and different luminosities. We find
that the continuum is well described by typical phenomenological models, like a
very strongly absorbed NPEX model. However, for a statistically acceptable
description of the hard X-ray spectrum we require two cyclotron resonant
scattering features (CRSF), one at ~35 keV and the other at ~50 keV. Even
though both features strongly overlap, the good resolution and sensitivity of
NuSTAR allows us to disentangle them at >=99.9% significance. This is the first
time that two CRSFs are seen in GX 301-2. We find that the CRSFs are very
likely independently formed, as their energies are not harmonically related
and, if it were a single line, the deviation from a Gaussian shape would be
very large. We compare our results to archival Suzaku data and find that our
model also provides a good fit to those data. We study the behavior of the
continuum as well as the CRSF parameters as function of pulse phase in seven
phase bins. We find that the energy of the 35 keV CRSF varies smoothly as
function of phase, between 30-38 keV. To explain this variation, we apply a
simple model of the accretion column, taking the altitude of the line-forming
region, the velocity of the in-falling material, and the resulting relativistic
effects into account. We find that in this model the observed energy variation
can be explained simply due to a variation of the projected velocity and
beaming factor of the line forming region towards us.
