© ESO, 2017

1. Introduction

Accreting X-ray pulsars are rotating highly magnetised (B> 10¹² G) neutron stars (NSs) that accrete matter from their normal companion star. Close to the NS, within the so-called Alfvén radius, magnetic pressure dominates the ram pressure of the ionised accretion flow such that the infalling plasma travels along the magnetic field lines towards the NS magnetic poles, where it is abruptly decelerated from a free-fall velocity of ~0.6c and heated to ≳10⁸ K. At the NS polar regions, this hot plasma forms structures that are often called accretion columns, and it radiates in the hard X-ray domain. To date, more than 200 accreting pulsars are known¹.

X-ray spectra of some accreting pulsars exhibit absorption-line-like features called cyclotron resonance scattering features (CRSFs), or cyclotron lines. The first such line was discovered in 1976 during a balloon observation of Her X-1 (Trümper et al. 1978). The cyclotron feature permits a direct measurement of the neutron star magnetic field strength because the latter is directly related to the cyclotron line centroid energy E_CRSF as E_CRSF = 11.6 [keV] B₁₂n/ (z + 1), where B₁₂ is the magnetic field strength in units of 10¹² gauss, n is the Landau level number, and z is the gravitational redshift at the line-forming region. The fundamental cyclotron line corresponds to n = 1. Cyclotron lines corresponding to higher n are referred to as harmonics.

The cyclotron feature has been found to change with the X-ray luminosity L_X in the spectra of a few accreting pulsars. The correlation between the CRSF centroid energy E_CRSF and L_X has so far been firmly confirmed in six sources. Four of them exhibit a positive correlation for the fundamental line, one for the first harmonic, and the last pulsar exhibits a negative correlation for the fundamental line. Below we provide a brief summary of the E_CRSF/L_X correlations observed in these sources. The correlation in Cep X-4 is described separately.

• Her X-1 is the first source where a positiveE_CRSF/L_X correlation was reported by Staubert et al. (2007) using RXTE observations. Klochkov et al. (2011) confirmed the correlation on a short timescale using a pulse-to-pulse analysis technique, which is also used in our work. Her X-1 demonstrates another phenomenon, however. The energy of the CRSF in this source apparently changes with time (Staubert et al. 2016, and references therein). This is the first source for which a secular change of the CRSF was reported.

• GX 304-1. Yamamoto et al. (2011) discovered a CRSF during the 2010 outburst and reported a marginal indication of a positive correlation between its energy and the X-ray flux. The correlation was observed with high significance using INTEGRAL (Klochkov et al. 2012; Malacaria et al. 2015) and RXTE (Rothschild et al. 2017) data.

• A 0535+26. The source spectrum has two CRSFs: the fundamental and the first harmonic (Kendziorra et al. 1994; Caballero et al. 2008). A positive E_CRSF/L_X correlation for the fundamental line was first found on a short timescale, using the pulse-to-pulse technique applied to RXTE and INTEGRAL/ISGRI observations of a source outburst (Klochkov et al. 2011). The correlation was subsequently confirmed on longer timescales by Sartore et al. (2015) using INTEGRAL/SPI data.

• Vela X-1. The source spectrum also has two CRSFs interpreted as the fundamental line and the first harmonic, but the fundamental feature is considerably weaker than the first harmonic. A positive E_CRSF/L_X correlation for the harmonic line is observed in the source spectrum (Fürst et al. 2014; La Parola et al. 2016).

• V 0332+53. The fundamental cyclotron line in the source spectrum was discovered during the 1983 outburst with Tenma (Makishima et al. 1990). During subsequent outbursts, two harmonics were detected (Pottschmidt et al. 2005; Kreykenbohm et al. 2005). The energy of the fundamental CRSF is negatively correlated with the source X-ray flux (Mowlavi et al. 2006; Tsygankov et al. 2010; Doroshenko et al. 2017). This is the only source for which a negative E_CRSF/L_X correlation is confirmed so far. There are indications that the source, similarly to Her X-1, might exhibit another type of variability: Cusumano et al. (2016) reported a change in the energy of the fundamental CRSF at the same level of luminosity.

There are also X-ray pulsars in which an indication for a correlation between the CRSF energy and luminosity is reported but so far not unambiguously confirmed, such as 4U 0115+63 (Müller et al. 2013b; Iyer et al. 2015) and 4U 1538-522 (Hemphill et al. 2014).

Here, we focus on the accreting pulsar Cep X-4 discovered in 1972 with the Orbiting Solar Observatory-7 (Ulmer et al. 1973) and re-discovered in 1988 with Ginga as GS 2138+56 during a one-month-long X-ray outburst (Makino & GINGA Team 1988). Several other outbursts have occurred since the discovery. Pulsations with a period of 66.2 s were detected (Koyama et al. 1991) during the 1988 outburst, which established the nature of the source as an accreting magnetised NS. The CRSF was discovered at around 30 keV during the same outburst (Mihara et al. 1991). A Be star was proposed as an optical counterpart during observations of the 1993 outburst (Roche et al. 1997). This identification was confirmed later by Bonnet-Bidaud & Mouchet (1998), who also estimated the distance to the source to be 3.8 ± 0.6 kpc. During the 2002 outburst, the CRSF energy was re-measured with RXTE to be at $\mathrm{30.}{\mathrm{7}}_{-1.9}^{\mathrm{+}\mathrm{1.8}}$ keV without a considerable dependence on X-ray luminosity, although the continuum was reported to become harder with increasing luminosity (McBride et al. 2007). The source went into outburst again in 2014 and was observed with the Nuclear Spectroscopic Telescope Array (NuSTAR, Harrison et al. 2013) by the two detectors FPMA and FPMB. Two NuSTAR observations of Cep X-4 were performed during this outburst (Fig. 1). The averaged fluxes in these observations differ by a factor of about 3.5, providing an opportunity to investigate the spectrum-luminosity dependence. The CRSF energy was found to be different between the two observations: the lower cyclotron line energy corresponds to the fainter observation, indicating a possible positive correlation between the CRSF energy and luminosity, although with only two data points (Fürst et al. 2015).

In our work, we consider the behaviour of the CRSF energy and of the continuum hardness ratio as a function of X-ray luminosity within individual NuSTAR observations of the 2014 outburst using the pulse-to-pulse variability of the source. The structure of the paper is as follows. In Sect. 2 we present a summary of the observations and data reduction. In Sect. 3 we describe our analysis of the pulse-phase averaged spectra and of the pulse-amplitude resolved spectra. Theoretical implications of the results are discussed in Sect. 4. Section 5 provides a summary and conclusions.

2. Observations and data reduction

NuSTAR observed the source twice during the 2014 outburst, on June 18–19 and on July 1–2 (Fig. 1). The first observation is near the outburst maximum with a total on-source time of 47.9 ks (ObsID 80002016002, MJD 56 826.92–56 827.84). It is referred to as observation I in the following. The second observation, referred to as observation II, was performed at the decay of the outburst with a total on-source time of 45.2 ks (ObsID 80002016004, MJD 56 839.43–56 840.31). Figure 1 also represents the Swift/BAT observations of the outburst².

Fig. 1 Light curve of the 2014 outburst of Cep X-4 taken with Swift/BAT (black data points). The NuSTAR observations analysed here are indicated with red crosses.

In our pulse-to-pulse analysis described in Sect. 3.2, all source spectra of the two NuSTAR detectors FPMA and FPMB are extracted from a circular region with a radius of 80′′ on the detector plane. The background spectra are extracted from a circular region of the same radius at the periphery of the NuSTAR field of view, similarly to Fürst et al. (2015).

The same extraction regions were initially also used for our pulse-averaged spectra. For a detailed investigation of the possible cyclotron harmonic feature around 55 keV, we used modified extraction regions as described in Sect. 3.1, however.

The data processing was performed using the standard nupipeline and nuproducts utilities and XSPEC³ v12.9 as part of HEASOFT⁴ 6.18.

3. Spectral analysis and results

3.1. Pulse-averaged spectrum

Fig. 2 Pulse-averaged NuSTAR spectra of obs. I (a) and obs. II (b). The blue and green data points correspond to FPMA and FPMB, respectively. Panel B) shows residuals to the model without cyclotron lines (red line). Panel C) shows residuals to the model setting the strengths of the two cyclotron lines to zero, but leaving the other parameters at their best-fit values obtained using the model with the lines (in order to illustrate the contribution of the lines to the continuum). Residuals to the best-fit model (black line in panel A)) are shown in panel D). The best-fit parameters of the model with the lines are listed in Table 1.

The spectral continuum of highly magnetised accreting neutron stars is typically described with phenomenological power law models with an exponential cutoff at higher energies (e.g. Coburn et al. 2002). Some models, however, include a detailed description of the physics of the polar emitting region (Wang & Frank 1981; Becker & Wolff 2007; Farinelli et al. 2012; Ceccobello et al. 2014). These physical models have been successfully applied in some cases (Ferrigno et al. 2009; Farinelli et al. 2016; Wolff et al. 2016). However, assumptions and approaches used in these models are still being tested and revised. Since our work is focused on the behaviour of the CRSF centroid energy and of the general properties of the continuum, like its hardness, it is sufficient to describe the spectra using the traditional phenomenological models, which also allows us to compare our results with the previous works.

McBride et al. (2007) showed that the continuum of the Cep X-4 pulse-averaged spectrum between 3.5 and 70 keV is described well by a power-law function (powlaw in XSPEC) with a Fermi-Dirac cut-off (fdcut, Tanaka 1986) $\mathrm{fdcut}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{\propto }{\mathit{E}}^{\mathrm{-}\mathrm{\Gamma }}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\exp \mathrm{(}{\frac{\mathit{E}\mathrm{-}{\mathit{E}}_{\mathrm{cut}}}{{\mathit{E}}_{\mathrm{fold}}}}^{\mathrm{)}}}\mathit{,}$(1)with the inclusion of a wide additive Gaussian emission line (gauss in XSPEC) at around 14 keV required to account for a bump-like feature in the residuals between 8 and 20 keV.

A broad line-like feature like this in the residuals between these energies, known as the 10 keV feature, is reported for some other accreting pulsars observed with RXTE, Ginga, BeppoSAX, and other missions (Coburn et al. 2002). To flatten the 10 keV feature in the Cep X-4 spectrum, we can use not only a broad emission line, but also an absorption line or a combination of the two. The additional wide Gaussian emission line does not correspond to a separate physical component. The 10 keV feature thus apparently reflects limitations of the simple phenomenological function used to describe the continuum of the high-quality NuSTAR spectrum of Cep X-4.

We describe the fundamental CRSF, clearly seen in the NuSTAR spectra at around 30 keV (Fig. 2), as a multiplicative absorption line with a Gaussian optical depth profile (gabs in XSPEC). The Fe K_α fluorescent line at around 6.4 keV is modelled as an additive Gaussian emission line. We also take the photoelectric absorption into account in the soft part of the spectra using the phabs model in XSPEC with photoelectric cross-sections from Verner et al. (1996) and element abundances from Wilms et al. (2000).

Panels B of Fig. 2 show sharp residuals around 55 keV. These residuals are present in the spectra of the two NuSTAR observations in the data of FPMA and FPMB. To describe this line-like feature as well as possible and to clarify whether the feature could be a harmonic of the CRSF, we analysed the spectrum in the energy range between 50 and 70 keV. First, we used background spectra extracted from different areas on the detector plane and found that the feature does not disappear. The dependence of the effective area on energy does not contain any peculiarities around 55 keV either. Second, we performed an analysis searching for an optimal size of the source extraction region at which the signal-to-noise ratio has a maximal value in the energy range of 50–70 keV. We found that the optimal size lies between 50′′ and 70′′ and extracted the source spectrum from a 60′′ radius region. Background spectra were extracted from a circular region of 160′′ radius at the periphery of the NuSTAR field of view. This adjustment of the extraction region does not significantly change the best-fit parameters of the continuum and of the fundamental cyclotron line. We then built the source spectrum describing the feature with the gabs model (Fig. 2). A very similar feature is present in the Suzaku spectrum of Cep X-4 taken during observations of the same outburst. Jaisawal & Naik (2015) analysed the Suzaku data and interpreted the feature as the first harmonic. We infer from the above that the feature at around 55 keV in the Cep X-4 spectrum is indeed the first harmonic of the CRSF.

Panel C of Fig. 2 shows the contribution of the lines to the continuum. Both lines are well localised and do not dominate the continuum, in contrast to the case of 4U 0115+63, for example, where inadequately modelled cyclotron lines could fully distort the continuum, as illustrated in Fig. 8 of Müller et al. (2013b) for instance.

The applied model can therefore be written as $\mathrm{ABS}\mathrm{\times }\mathrm{(}\mathrm{FDCUT}\mathrm{\times }\mathrm{CRSF}\mathrm{\times }{\mathrm{CRSF}}_{\mathit{h}}\mathrm{+}\mathrm{ADD}\mathrm{+}\mathrm{Fe}\mathrm{)}\mathit{,}$(2)where ABS is the photoelectric absorption (phabs), FDCUT is the fdcut function, ADD is the additional emission Gaussian line in the continuum, the CRSF and CRSF_h are gabs models for the fundamental CRSF and the harmonic, respectively, and Fe is the iron K_α complex modelled as a Gaussian emission line. The best-fit parameters of this model for the pulse-phase averaged spectra for the two NuSTAR observations are listed in Table 1. Owing to lower statistics, the spectra of observation II do not allow us to constrain E_cut. Since we prefer to use the same spectral model in both observations, we fixed E_cut in observation II to the best-fit value we obtained in observation I.

3.2. Pulse-to-pulse analysis

To investigate the variation of spectral parameters as a function of X-ray luminosity within individual observations, we applied the pulse-to-pulse technique elaborated in Klochkov et al. (2011) and Müller et al. (2013a). This method is aimed at investigating the spectrum-luminosity dependence on short timescales using the pulse-to-pulse variability that is generally exhibited by accreting pulsars (Staubert et al. 1980; Frontera et al. 1985; Tsygankov et al. 2007).

Fig. 3 NuSTAR pulse profiles in the 3–79 keV energy range from observations I (top) and II (bottom). The shaded areas indicate the phase interval considered for the pulse-to-pulse analysis. The part of the profile within the interval is referred to as pulse (see text).

Fig. 4 Distributions of counts in individual pulse (i.e. of pulse amplitudes) for observations I (top) and II (bottom). The red dashed lines indicate the boundaries of the amplitude bins we used to extract amplitude-resolved spectra. In each distribution, the total number of counts is evenly distributed among the bins.

Our pulse-to-pulse analysis includes the following stages. First, all photon arrival times are barycentered using the DE-200 solar system ephemeris. We measured pulse periods of 66.3352 ± 0.0003 and of 66.3336 ± 0.0002 s for observation I and observation II, respectively. These periods are used to accumulate the pulse profiles shown in Fig. 3. We then specify a phase interval containing the brightest part of the pulse profile, which is referred to as pulse (shaded areas in Fig. 3). To characterise the brightness of the pulse, we define the pulse amplitude as the total number of counts in the entire NuSTAR energy band (3–79 keV) within the specified phase interval of an individual pulsation cycle divided by the width of this interval. The number of counts in a single pulse is insufficient to extract a meaningful spectrum, therefore we group pulses with similar amplitudes together, similarly to Klochkov et al. (2011). For this purpose, we construct the distribution of pulses over amplitude, which is shown in Fig. 4. Although the average flux within each observation does not change much, the distribution clearly covers a dynamical range of a factor of two in amplitude, indicating a strong intrinsic pulse-to-pulse variability of the source. We then split the distribution into four bins and keep approximately the same number of counts in each bin (the dashed vertical lines in Fig. 4). The width of the bins is wider than the statistical fluctuation (variance) of pulse amplitudes, which guarantees that the observed variability reflects the intrinsic flux changes. Finally, we extract a broadband spectrum for each pulse amplitude bin.

Since the number of photons in each pulse amplitude bin spectrum is lower than 25% of the pulse-averaged spectrum, the harmonic becomes insignificant, especially in observation II. Since our purpose in this part is to investigate the variability of the energy of the fundamental CRSF, we extracted spectra between 3 and 48 keV to eliminate a possible influence of the residuals around the harmonic energy on the fundamental line. Therefore we here use the same model as in the analysis of the pulse-averaged spectrum described in Sect. 3.1, except for the harmonic component, and we do not fix the cut-off E_cut energy in observation II.

As mentioned above, to describe the continuum successfully, we introduced an additional component to flatten the 10 keV feature. To ensure that this component does not influence the CRSF centroid energy (see e.g. Boldin et al. 2013), we constructed the contour plots for the CRSF energy versus the parameters of the 10 keV component (Fig. 5). The confidence regions in the plots do not indicate any significant correlation between the parameters. We therefore conclude that the reported variations in the line energy are not associated with possible variations in the 10 keV feature.

Fig. 5 Confidence regions to investigate possible intrinsic correlation between the CRSF centroid energy ECRSF and the parameters of the Gaussian line describing the 10 keV feature in the continuum. The left and right columns show results for observations I and II, respectively. The contours correspond to the 1σ, 2σ, and 3σ confidence levels for the two parameters of interest. The colour scale indicates Δχ2: darker areas correspond to lower Δχ2. Aadd is in photons s-1 cm-2 keV-1 at 1 keV.

From a circular region with a radius of 80′′ (see Sect. 2), we extracted four spectra per detector for each observation, corresponding to the four pulse amplitude bins. In our fits, we left the relative normalisation of the two detectors free because there is a slight difference in their calibration. The model parameters of the spectra of the amplitude bins were coupled except for the parameters whose behaviour as a function of pulse amplitude is investigated. The investigated parameters are the energy of the fundamental cyclotron line E_CRSF and the power-law index Γ. The continuum hardness HR considered below is defined as the ratio of the flux in the 8–12 keV range to the flux in the 4–6 keV range. All other parameters (except for the relative normalisation, as mentioned above) were coupled. We verified that the coupled parameters do not show any noticeable variability with pulse amplitude when they are left free.

The dependencies of the fundamental CRSF centroid energy and of the continuum hardness ratio on pulse amplitude are shown in Fig. 6. The bottom axis represents the bin-averaged count rate in the 3–79 keV range, while the top axis represents the corresponding X-ray luminosity in the same energy band for the assumed distance of 3.8 kpc. The vertical and horizontal error bars on the data points in Fig. 6 correspond to 1σ uncertainties.

Since the emission of the pulsar is not isotropic, the flux in the pulse cannot be converted into luminosity by simply multiplying by 4πD, where D is the source distance. Therefore we first calculated the pulse-averaged flux of each observation. This flux may to some extend be considered as the source flux averaged over the solid angle and can therefore be converted into luminosity in the usual way. Using the pulse averaged flux, we calculated the average luminosity for each of the two observations. We then assumed that the relative changes in the pulse amplitude with respect to the average pulse amplitude in an observation are equal to the relative change in the luminosity with respect to the average luminosity in this observation. Under this assumption, we calculated the luminosity of each pulse or the average luminosity in each pulse amplitude bin, which is represented by the top axis of Fig. 6.

Fig. 6 Fundamental CRSF centroid energy (top) and the spectral hardness ratio (bottom) of Cep X-4 as a function of pulse amplitude, which is converted to luminosity on the top axis. The orange solid line represents a common linear fit to the data points for the two observations. The dashed lines mark 1σ cones of the allowed slope ranges from the separate fits to the two observations.

As shown in the top panel of Fig. 6, the data reveal a strong positive correlation between the fundamental CRSF energy and the X-ray luminosity. Furthermore, the correlation is present in the data of observation I alone (data points in the top right), where the differences in luminosity between the data points are caused by the pulse-to-pulse variability. The data of observation II alone (data points on the bottom left) have insufficient statistics to claim a correlation.

The common linear fit to all the data points is shown with a solid line and yields a slope of 0.034 ± 0.002 keV (cts/s)^-1. This slope is formally consistent with separate linear fits to the data points of the two observations with slopes of 0.043 ± 0.009 and of 0.021 ± 0.014 keV (cts/s)^-1 for observations I and II, respectively. The dashed lines mark 1σ cones of the allowed slope ranges from the separate fits of the two observations. It is therefore very likely that the E_CRSF/L_X dependence is roughly the same in the two observations and between them.

The bottom panel of Fig. 6 shows the spectral hardness as a function of pulse amplitude or luminosity. The hardness ratio clearly rises with increasing luminosity from observation II to observation I as well as within each observation. It can also be seen that the dependence is substantially flatter in observation I. Linear fits have slopes of (2.0 ± 0.2) × 10^-3 and of (9.2 ± 3.3) × 10^-3 (cts/s)^-1 for observations I and II, respectively (dashed lines). The difference in slopes is thus statistically significant.

To ensure the reliability of the dependencies, we fitted the spectra with alternative phenomenological models. Specifically, we used a power law with an exponential cut-off $\mathrm{highecut}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{\propto }\{\begin{array}{c}\\ & \mathit{E}\mathrm{\le }{\mathit{E}}_{\mathrm{cut}}\\ & \mathit{E}\mathit{>}{\mathit{E}}_{\mathrm{cut}}\end{array}$(3)(highecut in XSPEC, White et al. 1983) and an additive emission Gaussian line to flatten the 10 keV feature. We also added a gabs component in the continuum to smooth the derivative discontinuity at E_cut in this case. The CRSF was also described as a multiplicative absorption line with a Lorentzian optical depth profile (lorentz in XSPEC). In addition to our model, we thus tested three alternative models: the highecut*gabs+gauss model for the continuum with either the gabs or the lorentz model for the cyclotron line, and the fdcut+gauss model for the continuum with the lorentz model for the cyclotron line. The iron line was described as an additive emission gauss in all models. All the dependencies described above are found to be present in all the alternative models, and all our conclusions remain valid independently of the choice of a specific spectral function.

4. Discussion

4.1. Timescales of the correlations

The aforementioned dependencies between the fundamental CRSF energy, the spectral hardness ratio, and the X-ray luminosity in separate observations are results from the pulse-to-pulse analysis, which makes use of an intrinsic source pulse-to-pulse variability. In other words, the dependencies for separate observations are present on a timescale corresponding to one cycle of X-ray pulsations, that is, on the order of a minute.

In the case of the E_CRSF/L_x correlation, however, the same linear dependence is likely to be present on both short (a pulse cycle) and long (outburst duration) timescales. We can therefore assume that the physical processes behind the positive E_CRSF/L_x dependence observed in Cep X-4 are linear, in a first approximation, on timescales from a minute to a dozen days, or in other words, the properties of these processes remain constant on these timescales.

In the case of the spectral hardness, on the other hand, such a linearity is not observed, since the linear dependencies in separate observations are not consistent with the common linear dependence. Thus, the shape of the HR/L_x dependence is likely to be more complex on long timescales of days, while on short timescales it seems to be linear.

4.2. Accretion regimes

As mentioned in Sect. 1, positive and negative correlations between the CRSF energy E_CRSF and the X-ray luminosity L_x are observed in some X-ray pulsars. A certain threshold or critical luminosity L_cr ~ 10³⁷ erg s^-1 that reflects the mass accretion rate Ṁ is believed to separate the accretion modes that are characterised by different signs of the correlation (Basko & Sunyaev 1976; Becker et al. 2012; Mushtukov et al. 2015b). A positive E_CRSF/L_x correlation is observed below L_cr, while a negative correlation is observed above it.

Above L_cr, the radiation pressure in the polar emitting region of the NS dominates gas pressure such that the infalling plasma is decelerated mainly through interactions of the plasma electrons with photons. This mode is often referred to as supercritical. The hight of the emitting region, which in this case is believed to have roughly a column-like shape (accretion column), is expected to grow with mass accretion rate (Davidson 1973; Basko & Sunyaev 1976; Wang & Frank 1981; Burnard et al. 1991; Postnov et al. 2015; Mushtukov et al. 2015a). If the CRSF is formed in a region associated with the accretion column as suggested in Becker et al. (2012), for example, the cyclotron line energy will be anti-correlated with the height of the column. Indeed, the strength of the NS dipole magnetic field decreases with increasing height above the stellar surface and so does the cyclotron line energy, which is proportional to the field strength.

Alternatively, the CRSF might form in the radiation reflected by the NS surface surrounding the accretion column, as was suggested by Poutanen et al. (2013). A growth of the column would widen the stellar surface that is illuminated by the radiation from the column, including the areas farther away from the magnetic pole where magnetic field is weaker. As a result, an anti-correlation between the line energy and the accretion rate should be observed.

In accreting pulsars with luminosities below L_cr, which are sometimes referred to as subcritical sources, the radiation pressure is unimportant compared to gas pressure. In this case, the infalling plasma has to be decelerated through Coulomb collisions, collective plasma effects, or through a collisionless shock that forms above the stellar surface. The following ideas have been proposed in the literature to explain the positive E_CRSF/L_x correlation in the subcritical case. In Staubert et al. (2007) and Becker et al. (2012), the line-forming region was again assumed to be associated with the stopping region such that variations in the height of this region above the stellar surface lead to corresponding changes in the CRSF energy. It was shown that contrary to the supercritical case, the stopping region approaches the surface with increasing mass accretion rate. Therefore an anti-correlation between the CRSF centroid energy and luminosity is expected. An alternative idea has been elaborated on in Mushtukov et al. (2015c), who suggested that the variations in the energy of the cyclotron line are due to the Doppler effect. The dependence of the mean velocity of the infalling matter in the line-forming region on the mass accretion rate is responsible for the positive E_CRSF/L_x correlation in this model.

4.3. Colissionless shock model

Here, we consider a so far poorly explored possibility that infalling matter decelerates in a collisionless shock forming above the NS surface (Langer & Rappaport 1982). This is likely to be the main stopping agent of the accretion flow when the photon luminosities are low enough for radiation pressure to be unimportant for the accretion flow breaking (but in the CRSF, where the optical depth is very high even at low luminosities, see Bykov & Krasilshchikov 2004). This model has recently been used to describe the observed non-linear dependencies of the CRSF parameters on luminosity in another subcritical X-ray pulsar GX 304−1 (Rothschild et al. 2017). We develop this model in more detail so that it can also be used to describe the evolution of the spectral hardness ratio (HR) with luminosity.

4.3.1. CRSF energy – luminosity dependence

The height of the collisionless shock forming in the infalling flow near the NS surface is expected to depend on electron density n_e as H_s ∝ 1 /n_e (Shapiro & Salpeter 1975). The electron density n_e ∝ Ṁ/A, where Ṁ is the mass accretion rate and A is the accretion area of the accretion flow on the stellar surface depending on the magnetospheric radius R_m as 1 /R_m in a dipole field approximation. The magnetospheric radius R_m ∝ Ṁ^{− x}, where x = 2/7 for disc accretion or Bondi quasi-spherical accretion, and x = 2/11 for quasi-spherical settling accretion (Shakura et al. 2012), so that A ∝ Ṁ^x and ${\mathit{H}}_{\mathrm{s}}\mathrm{\propto }\mathrm{1}\mathit{/}{\mathit{n}}_{\mathrm{e}}\mathrm{\propto }\mathit{A}\mathit{/}\mathit{Ṁ}\mathrm{\propto }\mathit{Ṁ}\mathit{x}\mathrm{-}\mathrm{1}\mathrm{\propto }\mathit{Ṁ}\mathrm{-}\mathit{\alpha }\mathit{,}$(4)where α = 1−x = 5/7 for disc accretion or Bondi quasi-spherical accretion, and α = 9/11 for quasi-spherical settling accretion. Assuming the disc accretion in Cep X-4, we use α = 5/7 in all formulas below.

With increasing mass accreting rate, the height of the collisionless shock thus decreases such that the line-forming region moves towards stronger magnetic fields and the cyclotron line shifts to higher energies. In a dipole field, E_CRSF(r) ∝ B(r) ∝ 1 /r³, and considering that r = R_NS + H_s, where R_NS is the NS radius and H_s is a function of accretion rate, we express the energy of the cyclotron line as a function of mass accretion rate: ${\mathit{E}}_{\mathrm{CRSF}}\mathrm{\left(}\mathit{Ṁ}\mathrm{\right)}\mathrm{=}{\mathit{E}}_{\mathrm{0}}{\left(\frac{{\mathit{R}}_{\mathrm{NS}}}{{\mathit{H}}_{\mathrm{s}}\mathrm{\left(}\mathit{Ṁ}\mathrm{\right)}\mathrm{+}{\mathit{R}}_{\mathrm{NS}}}\right)}^{\mathrm{3}}\mathit{,}$(5)where E₀ is the cyclotron line energy corresponding to the magnetic field at the NS surface on the magnetic poles. Assuming that the observed X-ray flux F_x is directly proportional to the mass accretion rate Ṁ and using Eq. (4), we obtain that ${\mathit{H}}_{\mathrm{s}}\mathrm{\propto }{\mathit{F}}_{\mathrm{x}}^{\mathrm{-}\mathit{\alpha }}$. Introducing the dimensionless quantity H_s/R_NS, we can write ${\mathit{H}}_{\mathrm{s}}\mathit{/}{\mathit{R}}_{\mathrm{NS}}\mathrm{=}\mathit{K}{\mathit{F}}_{\mathrm{x}}^{\mathrm{-}\mathit{\alpha }}\mathit{,}$(6)and finally, ${\mathit{E}}_{\mathrm{CRSF}}\mathrm{\left(}{\mathit{F}}_{\mathrm{x}}\mathrm{\right)}\mathrm{=}{\mathit{E}}_{\mathrm{0}}\mathrm{(}\mathit{K}{\mathit{F}}_{\mathrm{x}}^{\mathrm{-}\mathit{\alpha }}\mathrm{+}\mathrm{1}{\mathrm{)}}^{-3}\mathit{,}$(7)where the factor K determines the height of the collisionless shock. E₀ and K are found as fit parameters of observational data. In our case, E₀ = 31.4 ± 0.2 keV and K = 0.024 ± 0.003 (erg/s)^α.

As shown in Fig. 7, the collisionless shock model represented with the solid curve describes the observed positive E_CRSF/L_x correlation reasonably well. We note the apparent turn-up of CRSF energies at the high end of the explored X-ray luminosity range, however. The model presented above does not predict such a behaviour, but only yields a flattening at high luminosities. If this were real, the visible turn-up may indicate gradual deviations from the model being considered here. Indeed, with increasing mass accretion rate, the radiation pressure becomes more and more important, so that the simplest version of the collisionless shock may not be fully applicable for the brightest states. This intermediate case between the two regimes is the most difficult to treat, and apparently full-radiation hydro simulations are needed to investigate this interval of luminosities. We also note that the limited number of data points and their accuracy do not exclude the possibility that the indication of the second derivative is a pure statistical or systematical effect.

Substituting measured F_x, found from Eq. (7) K and assuming R_NS = 10 km in Eq. (6), we obtain averaged (over each observation) shock heights of ~200 and 100 m above the NS surface for observations II (in the outburst decay) and I (near the outburst maximum), respectively.

Fig. 7 Centroid energy of the fundamental cyclotron line in Cep X-4 as a function of X-ray luminosity (data points with 1σ error bars). The solid curve represents a fit using our model assuming a collisionless shock (see text), ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathit{/}\mathrm{d}\mathrm{.}\mathrm{o}\mathrm{.}\mathrm{f}\mathrm{.}\mathrm{=}\mathrm{1.75}\mathit{/}\mathrm{7}$.

4.3.2. Hardness ratio – luminosity dependence

The X-ray spectral analysis presented in Sect. 3 revealed a strong and apparently non-linear correlation between the spectral hardness ratio and X-ray luminosity in Cep X-4. Generally, an increase of spectral hardness with luminosity is expected to be due to the increase of the electron number density n_e with increasing leading to more effective Comptonisation of the emerging X-ray spectrum. In the model with a collisionless shock, most of the energy is released near the NS surface where the continuum spectrum Comptonised by electrons between the NS surface and the shock is produced. Numerical simulations performed by Bykov & Krasilshchikov (2004) suggest that at a given accretion rate neither T_e and n_e change significantly downstream of the shock. The shape of the emerging spectrum depends on the energy exchange between photons and electrons (the dimensionless Comptonisation y-parameter, Kompaneets 1956), which in the strong magnetic field is (Basko & Sunyaev 1975) $\mathit{y}\mathrm{=}\frac{\mathrm{2}}{\mathrm{15}}\frac{\mathit{k}{\mathit{T}}_{\mathrm{e}}}{{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}\mathrm{Max}\mathrm{\left\{}\mathit{\tau ,}{\mathit{\tau }}^{\mathrm{2}}\mathrm{\right\}}\mathit{,}$(8)where τ is the characteristic optical depth of the region and all other notations are standard. Consider a cylindrical geometry of the accretion region with the polar cap area A and height H_s. In real X-ray observations, the geometry is averaged over the viewing angle by the rotation of the NS, therefore the appropriate characteristic scale of the problem is ${\mathit{r}}_{\mathrm{eff}}\mathrm{=}\mathrm{(}\mathit{A}{\mathit{H}}_{\mathrm{s}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{.}$(9)Comptonisation is effective for y ~ 1, and the number of scatterings in this region is $\mathrm{~}{\mathit{\tau }}_{\mathrm{eff}}^{\mathrm{2}}\mathrm{\propto }{\mathit{n}}_{\mathrm{e}}^{\mathrm{2}}{\mathit{r}}_{\mathrm{eff}}^{\mathrm{2}}$. The effective y-parameter of the problem therefore scales as ${\mathit{y}}_{\mathrm{eff}}\mathrm{\propto }{\mathit{\tau }}_{\mathrm{eff}}^{\mathrm{2}}\mathrm{\propto }{\mathit{n}}_{\mathrm{e}}^{\mathrm{2}}{\mathit{r}}_{\mathrm{eff}}^{\mathrm{2}}\mathrm{\propto }{\left(\frac{\mathit{Ṁ}}{\mathit{A}}\right)}^{\mathrm{2}}\mathrm{(}\mathit{A}{\mathit{H}}_{\mathrm{s}}{\mathrm{)}}^{\frac{\mathrm{2}}{\mathrm{3}}}\mathrm{\propto }\mathit{Ṁ}\frac{\mathrm{10}}{\mathrm{7}}\mathit{Ṁ}\frac{\mathrm{4}}{\mathrm{21}}\mathit{Ṁ}\mathrm{-}\frac{\mathrm{10}}{\mathrm{21}}\mathrm{\propto }\mathit{Ṁ}\frac{\mathrm{24}}{\mathrm{21}}\mathit{,}$(10)where we have made use of the standard dependence A ∝ Ṁ^2/7 (a dipole magnetic field geometry). Next, we use the standard solution of the unsaturated Comptonisation problem (Pozdnyakov et al. 1983), where the emerging spectral flux at photon energies hν<kT_e can be written as ${\mathit{F}}_{\mathit{\nu }}\mathrm{=}\mathit{C}{\mathit{\nu }}^{\mathrm{-}\mathit{\alpha }}$(11)with the spectral index $\mathit{\alpha }\mathrm{=}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\mathrm{+}\sqrt{\frac{\mathrm{9}}{\mathrm{4}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{15}{\mathit{y}}_{\mathrm{eff}}}}\mathrm{·}$(12)(Here we made use of the decrease in the efficiency of the energy exchange between electrons and photons in a strong magnetic field, see Eq. (8) above). The spectral hardness ratio is determined in the usual way: $\mathrm{HR}\mathrm{=}\frac{\underset{{\mathit{\nu }}_{\mathrm{3}}}{\overset{{\mathit{\nu }}_{\mathrm{4}}}{\mathrm{\int }}}{\mathit{F}}_{\mathit{\nu }}\mathrm{d}\mathit{\nu }}{\underset{{\mathit{\nu }}_{\mathrm{1}}}{\overset{{\mathit{\nu }}_{\mathrm{2}}}{\mathrm{\int }}}{\mathit{F}}_{\mathit{\nu }}\mathrm{d}\mathit{\nu }}\mathit{,}$(13)where, in our case, the frequencies ν₁, ν₂, ν₃, and ν₄ correspond to energies of 4, 6, 8, and 12 keV, respectively. Plugging Eqs. (11) and (12) into Eq. (13) and taking Eq. (10) into account, we arrive at $\mathrm{HR}\mathrm{=}\frac{\underset{{\mathit{\nu }}_{\mathrm{3}}}{\overset{{\mathit{\nu }}_{\mathrm{4}}}{\mathrm{\int }}}{\mathit{\nu }}^{\mathrm{3}\mathit{/}\mathrm{2}\mathrm{-}\sqrt{\mathrm{9}\mathit{/}\mathrm{4}\mathrm{+}\mathrm{2}\mathit{/}\mathrm{\left(}\mathrm{15}{\mathit{y}}_{\mathrm{eff}}\mathrm{\right)}}}\mathrm{d}\mathit{\nu }}{\underset{{\mathit{\nu }}_{\mathrm{1}}}{\overset{{\mathit{\nu }}_{\mathrm{2}}}{\mathrm{\int }}}{\mathit{\nu }}^{\mathrm{3}\mathit{/}\mathrm{2}\mathrm{-}\sqrt{\mathrm{9}\mathit{/}\mathrm{4}\mathrm{+}\mathrm{2}\mathit{/}\mathrm{\left(}\mathrm{15}{\mathit{y}}_{\mathrm{eff}}\mathrm{\right)}}}\mathrm{d}\mathit{\nu }}\mathrm{=}\frac{\underset{{\mathit{\nu }}_{\mathrm{3}}}{\overset{{\mathit{\nu }}_{\mathrm{4}}}{\mathrm{\int }}}{\mathit{\nu }}^{\mathrm{3}\mathit{/}\mathrm{2}\mathrm{-}\sqrt{\mathrm{9}\mathit{/}\mathrm{4}\mathrm{+}\mathrm{2}\mathit{/}\mathrm{\left(}\mathrm{15}{\mathit{K}}_{\mathit{y}}\mathit{Ṁ}\mathrm{24}\mathit{/}\mathrm{21}\mathrm{\right)}}}\mathrm{d}\mathit{\nu }}{\underset{{\mathit{\nu }}_{\mathrm{1}}}{\overset{{\mathit{\nu }}_{\mathrm{2}}}{\mathrm{\int }}}{\mathit{\nu }}^{\mathrm{3}\mathit{/}\mathrm{2}\mathrm{-}\sqrt{\mathrm{9}\mathit{/}\mathrm{4}\mathrm{+}\mathrm{2}\mathit{/}\mathrm{\left(}\mathrm{15}{\mathit{K}}_{\mathit{y}}\mathit{Ṁ}\mathrm{24}\mathit{/}\mathrm{21}\mathrm{\right)}}}\mathrm{d}\mathit{\nu }}\mathit{,}$(14)where the parameter K_y is found from fitting of the observed HR change with X-ray flux.

Fig. 8 Similar to Fig. 7, but for the spectral hardness ratio. The solid curve represents a model expressed by Eq. (14), ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathit{/}\mathrm{d}\mathrm{.}\mathrm{o}\mathrm{.}\mathrm{f}\mathrm{.}\mathrm{=}\mathrm{1.36}\mathit{/}\mathrm{7}$.

The best-fit model is presented in Fig. 8. Clearly, the model of Eq. (14) describes the observed behaviour of hardness ratio with increasing X-ray luminosity of Cep X-4 well. The best-fit value of the sole parameter is K_y ≃ 0.14, which corresponds to an effective y-parameter range ~0.2−1.2 during the observed luminosity change in Cep X-4. These values of y_eff are consistent with the unsaturated Comptonisation regime, which justifies the use of Eqs. (11) and (12) for the model fitting.

To illustrate the qualitative analysis given above, we performed numerical simulations of radiative transfer in the scattering medium above the neutron star polar cap using a two-dimensional Feautrier method to obtain the hardness ratio. The details of the numerical calculations are presented in Appendix A.

Fig. 9 Model spectra with the observed fdcut continua for observations I and II. (I): Ṁ = 5 × 1016 g s-1, the polar cap radius r0 = 0.38 km, the height of the column z0 = 77 m, ne = 2.66 × 1021 cm-3. (II): Ṁ = 1.5 × 1016 g s-1, r0 = 0.32 km, z0 = 180 m, electron density ne = 1.12 × 1021 cm-3. The temperature of the scattering electrons is Te = 5 keV in both cases.

Figure 9 shows the emerging model spectra from the accretion column calculated as specific flux per unit surface integrated over the entire column for observations I (dashed line) and II (dash-dotted line). To compare them with the observed continua (solid and dotted lines), the modelled specific luminosity from two columns was divided by some factor to match the continuum given by Eq. (1). This factor is equal to 4πd², where d ~ 4 kpc is the distance to the source.

The hardness ratios of the modelled spectra are calculated to be 1.93 and 1.74 for observations I and II, respectively, which is very similar to the observed values (see Fig. 6). This shows that the model spectra are in good agreement with the observed continua, confirming the expected hardening of the continuum with increasing mass accretion rate. It should be noted that here we did not study the E_CRSF/L_x correlation obtained from the simulations because the accuracy of the calculations is not high enough to analyse the line centroid energy. The application of the Feautrier method for the cyclotron line energy variations will be explored in the future work.

We therefore conclude that changing the electron density downstream of the collisionless shock above the NS polar cap is consistent with the observed hardness ratio evolution with X-ray luminosity in Cep X-4.

5. Summary

We have investigated the spectral-luminosity variability in the X-ray pulsar Cep X-4 using NuSTAR observations of the 2014 outburst. The intrinsic pulse-to-pulse variability of Cep X-4 allowed us to reveal strong positive E_CRSF/L_x correlations on short timescales of one cycle of X-ray pulsations (minutes) as well as on long timescales of weeks (at least 12 days during which the source becomes fainter by a factor of about 3.5). Moreover, linear fits to the data points seem to have the same slope, which might suggest linearity of the processes behind the E_CRSF/L_x correlations.

We also confirm that the continuum becomes harder with X-ray luminosity, that is, the slope of the power law decreases, or in other words, the spectral hardness ratio (F_{8−12 keV}/F_{4−6 keV}) increases with X-ray flux. However, the slope of the linear fit to the points of the measured hardness ratio clearly changes from the first observation to the second, which implies a more complex form of the dependence than a linear one on long timescales (during at least one outburst).

The evolution of the spectral hardness ratio with X-ray luminosity is successfully described by a model in which a collisionless shock above the neutron star surface is formed. This model also describes the observed positive E_CRSF/L_x correlation reasonably well. Within this model, heights of cyclotron line emerging layers were estimated.

Furthermore, we found a harmonic of the cyclotron feature in both NuSTAR observations. We thus confirm the harmonic of the cyclotron line in the Cep X-4 spectrum, as previously found by Jaisawal & Naik (2015) in Suzaku data.

Acknowledgments

The research is supported by the joint DFG grant KL 2734/2-1 and Wi 1860 11-1 and the RFBR grants 14-02-91345 and 18-502-12025. We thank the anonymous referee for the useful comments and suggestions that substantially improved the manuscript. The research used the data obtained from HEASARC Online Service provided by the NASA/GSFC.

References

Appendix A: Numerical modelling

The Feautrier method is a powerful tool of radiation transfer calculations in regions with given structure (temperature and density distributions) and is widely used for spectrum calculations (Nagel 1981a,b; Nagel & Meszaros 1985; Alexander et al. 1989). In these works a one-dimensional model for two particular geometries was applied. In a collisionless shock model, where the radiation from the top of the column could be comparable to that of the sidewalls, it is necessary to take a two-dimensional radiative transfer into account.

For this reason, we used a generalised two-dimensional Feautrier method (Cannon 1970) to calculate the continuum and (with less accuracy) the fundamental CRSF feature. The calculation domain is taken to be a cylindrical column filled with plasma of constant temperature and density. The radiative transfer for two polarisation modes is used in the geometrical optics limit. We employ cylindrical coordinates, assuming axial symmetry; z = 0 at the NS surface and r = 0 at the centre of the column.

The radiative transfer equation for the scattering atmosphere reads $\frac{{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{i}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{i}}}\left(\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{u}}_{\mathit{i}}}{\mathit{\partial }{\mathit{r}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }{\mathit{u}}_{\mathit{i}}}{\mathit{\partial r}}\right)\mathrm{+}\frac{{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{i}}}{{\mathit{\sigma }}_{\mathit{i}}}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{u}}_{\mathit{i}}}{\mathit{\partial }{\mathit{z}}^{\mathrm{2}}}\mathrm{-}{\mathit{\sigma }}_{\mathit{i}}{\mathit{u}}_{\mathit{i}}\mathrm{+}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{S}}_{\mathit{ij}}{\mathit{u}}_{\mathit{j}}\mathrm{=}\mathrm{0.}$(A.1)Here θ is the angular coordinate, ${\mathit{u}}_{\mathit{i}}\mathrm{\left(}\mathit{\omega ,\theta ,r,z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left[{\mathit{I}}_{\mathit{i}}\mathrm{\left(}\mathit{\omega ,\theta ,r,z}\mathrm{\right)}\mathrm{-}{\mathit{I}}_{\mathit{i}}\mathrm{(}\mathit{\omega ,}\mathrm{-}\mathit{\theta ,r,z}\mathrm{)}\right]$, where I is the intensity; σ_i = σ_i(ω,θ,z) is the scattering coefficient independent of the r-coordinate because the magnetic field is assumed to be oriented along the z-axis.

The redistribution matrix elements S_ij are calculated using the equations for differential scattering cross sections from Bussard et al. (1986) derived in a fully relativistic approximation. Thus the second-order Doppler effect and the non-harmonicity of the electron Landau levels are included in our calculations. We restricted ourselves to considering the n = 1 Landau level for the intermediate excited state and did not take spin-flip transitions into account. The derivation of the fully relativistic normal modes including vacuum and plasma effects requires special treatment, and here we used the non-relativistic vacuum-plasma normal modes from Nagel & Meszaros (1985), for the sake of simplicity.

The calculations were performed assuming the mass accretion rates Ṁ = 5 × 10¹⁶ g s^-1 and Ṁ = 1.5 × 10¹⁶ g s^-1, corresponding to the averaged luminosities of observations I and II, respectively. The electron temperature of the column was set to 5.8 keV for both accretion rates, which is close to the electron temperature obtained in numerical modelling (Bykov & Krasilshchikov 2004). The electron number density is a function of Ṁ and the polar cap radius (Eq. (4)).

The bottom of the column (polar cap with area A where most of the accretion power is released) serves as the source of photons. Because scattering is important in a thin emission region above the polar cap, we assume a modified black-body spectrum of seed photons. Other boundary conditions are that there is no incident radiation onto the top and the sidewall of the column, ∂u/∂r = 0 along the z-axis.

The temperature of the seed radiation is a function of the accretion rate. Taking into account that L ∝ T⁴A ∝ Ṁ and using the standard dependence for the polar cap area, A ∝ Ṁ^2/7, the temperature ratio for different mass accretion rates is $\frac{{\mathit{T}}_{\mathrm{1}}}{{\mathit{T}}_{\mathrm{2}}}\mathrm{=}{\left(\frac{\dot{{\mathit{M}}_{\mathrm{1}}}}{\dot{{\mathit{M}}_{\mathrm{2}}}}\right)}^{\mathrm{5}\mathit{/}\mathrm{28}}\mathrm{·}$(A.2)The temperature of the radiation coming from the polar cap was fixed to 12 keV for observation II with the lower mass accretion rate.

The results of the simulation for the two accretion rates are shown in Fig. 9 in the main text.

All Tables

All Figures

	Fig. 1 Light curve of the 2014 outburst of Cep X-4 taken with Swift/BAT (black data points). The NuSTAR observations analysed here are indicated with red crosses.
In the text

Fig. 2 Pulse-averaged NuSTAR spectra of obs. I (a) and obs. II (b). The blue and green data points correspond to FPMA and FPMB, respectively. Panel B) shows residuals to the model without cyclotron lines (red line). Panel C) shows residuals to the model setting the strengths of the two cyclotron lines to zero, but leaving the other parameters at their best-fit values obtained using the model with the lines (in order to illustrate the contribution of the lines to the continuum). Residuals to the best-fit model (black line in panel A)) are shown in panel D). The best-fit parameters of the model with the lines are listed in Table 1.

In the text

	Fig. 3 NuSTAR pulse profiles in the 3–79 keV energy range from observations I (top) and II (bottom). The shaded areas indicate the phase interval considered for the pulse-to-pulse analysis. The part of the profile within the interval is referred to as pulse (see text).
In the text

	Fig. 4 Distributions of counts in individual pulse (i.e. of pulse amplitudes) for observations I (top) and II (bottom). The red dashed lines indicate the boundaries of the amplitude bins we used to extract amplitude-resolved spectra. In each distribution, the total number of counts is evenly distributed among the bins.
In the text

Fig. 5 Confidence regions to investigate possible intrinsic correlation between the CRSF centroid energy ECRSF and the parameters of the Gaussian line describing the 10 keV feature in the continuum. The left and right columns show results for observations I and II, respectively. The contours correspond to the 1σ, 2σ, and 3σ confidence levels for the two parameters of interest. The colour scale indicates Δχ2: darker areas correspond to lower Δχ2. Aadd is in photons s-1 cm-2 keV-1 at 1 keV.

In the text

	Fig. 6 Fundamental CRSF centroid energy (top) and the spectral hardness ratio (bottom) of Cep X-4 as a function of pulse amplitude, which is converted to luminosity on the top axis. The orange solid line represents a common linear fit to the data points for the two observations. The dashed lines mark 1σ cones of the allowed slope ranges from the separate fits to the two observations.
In the text

	Fig. 7 Centroid energy of the fundamental cyclotron line in Cep X-4 as a function of X-ray luminosity (data points with 1σ error bars). The solid curve represents a fit using our model assuming a collisionless shock (see text), ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathit{/}\mathrm{d}\mathrm{.}\mathrm{o}\mathrm{.}\mathrm{f}\mathrm{.}\mathrm{=}\mathrm{1.75}\mathit{/}\mathrm{7}$.
In the text

	Fig. 8 Similar to Fig. 7, but for the spectral hardness ratio. The solid curve represents a model expressed by Eq. (14), ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathit{/}\mathrm{d}\mathrm{.}\mathrm{o}\mathrm{.}\mathrm{f}\mathrm{.}\mathrm{=}\mathrm{1.36}\mathit{/}\mathrm{7}$.
In the text

	Fig. 9 Model spectra with the observed fdcut continua for observations I and II. (I): Ṁ = 5 × 1016 g s-1, the polar cap radius r0 = 0.38 km, the height of the column z0 = 77 m, ne = 2.66 × 1021 cm-3. (II): Ṁ = 1.5 × 1016 g s-1, r0 = 0.32 km, z0 = 180 m, electron density ne = 1.12 × 1021 cm-3. The temperature of the scattering electrons is Te = 5 keV in both cases.
In the text