Abstract We report on the energy-resolved timing and phase-resolved spectral analysis of X-ray emission from PSR J0659+1414 observed with XMM-Newton and NuSTAR. We find that the new data rule out the previously suggested model of the phase-dependent spectrum as a three-component (two blackbodies + power law) continuum, which shows large residuals between 0.3 and 0.7 keV. Fitting neutron star atmosphere models or several blackbodies to the spectrum does not provide a better description of the spectrum and requires spectral model components with unrealistically large emission region sizes. The fits improve significantly if we add a phase-dependent absorption feature with central energy 0.5–0.6 keV and equivalent width up to ≈50 eV. We detected the feature for about half of the pulse cycle. Energy-resolved pulse profiles support the description of the spectrum with a three-component continuum and an absorption component. The absorption feature could be interpreted as an electron cyclotron line originating in the pulsar magnetosphere and broadened by the nonuniformity of the magnetic field along the line of sight. The significant phase variability in the thermal emission from the entire stellar surface may indicate multipolar magnetic fields and a nonuniform temperature distribution. The strongly pulsed nonthermal spectral component detected with NuSTAR in the 3–20 keV range is well fit by a power-law model with a photon index Γ = 1.5 ± 0.2.