© ESO, 2014

1. Introduction

The accreting binary X-ray pulsar Her X-1 shows strong variability on several different timescales: the 1.24 s spin period of the neutron star, the 1.7-day binary period, the 35-day period, and the 1.65-day period of the pre-eclipse dips. The 35-day ON-OFF modulation can be understood as due to the precession of a warped accretion disk. Because of the high inclination (i> 80°) of the binary we see the disk nearly edge-on. The precessing warped disk therefore covers the central X-ray source during a substantial portion of the 35-day period. Furthermore, a hot X-ray heated accretion disk corona reduces the X-ray signal (energy independently) by Compton scattering whenever it intercepts our line of sight to the neutron star. As a result the X-ray source is covered twice during a 35-day cycle. Another 35 d modulation is present in the systematic variation of the shape of the 1.24 s pulse profile. It has been suggested (Trümper et al. 1986) that the reason for this is free precession of the neutron star, leading to a systematic change in our viewing angle to the X-ray emitting regions. By comparing the variation in flux (the turn-ons) and the variations in pulse shape, Staubert et al. (2009) had concluded, that if precession of the NS exists a strong interaction between the NS and the accretion disk is required, allowing a synchronization of the two clocks to nearly the same frequency through a closed loop physical feed-back (for which there is independent evidence). Further analysis of the variations in pulse profiles (Staubert et al. 2013), however, has shown that the histories of the turn-ons and of the variations in pulse shape are identical, with correlated variations even on short timescales (~300 d).

It is believed that the X-ray spectrum emerges from the hot regions around the magnetic poles where the accreted material is channeled by the ~10¹² G magnetic field down to the surface of the NS. The height of the accretion mound is thought to be a few hundred meters. If the magnetic and spin axes of the neutron star are not aligned, the view of a terrestrial observer is modulated at the rotation frequency of the star.

The X-ray spectrum of Her X-1 is characterized by a power law continuum with exponential cut-off and an apparent line-like feature. The continuum is believed to be due to thermal bremsstrahlung radiation from the ~10⁸ K hot plasma modified by Comptonization (Becker & Wolff 2007; Becker et al. 2012). The line feature was discovered in 1976 in a balloon observation (Trümper et al. 1978). This feature is now generally accepted as an absorption feature around 40 keV due to resonant scattering of photons off electrons on quantized energy levels (Landau levels) in the Terragauss magnetic field at the polar cap of the neutron star. The feature is therefore often referred to as a cyclotron resonant scattering feature (CRSF). The energy spacing between the Landau levels is given by E_cyc = ħeB/m_ec = 11.6 keV B₁₂, where B₁₂ = B/ 10¹² G, providing a direct method of measuring the magnetic field strength at the site of the emission of the X-ray spectrum. The observed line energy is subject to gravitational redshift z at the location where the line is formed, such that the magnetic field may be estimated by B₁₂ = (1 + z) E_obs/11.6 keV, with E_obs being the observed cyclotron line energy. The discovery of the cyclotron feature in the spectrum of Her X-1 provided the first ever direct measurement of the magnetic field strength of a neutron star, in the sense that no other model assumptions are needed. Originally considered an exception, cyclotron features are now known to be rather common in accreting X-ray pulsars, with ~20 binary pulsars now being confirmed cyclotron line sources, with several objects showing multiple lines (up to four harmonics in 4U 0115+63). Reviews are given by e.g., Coburn et al. (2002), Staubert (2003), Heindl et al. (2004), Terada et al. (2007), Wilms (2012), Caballero & Wilms (2012). Theoretical calculations of cyclotron line spectra have been performed either analytically (Ventura et al. 1979; Nagel 1981; Nishimura 2008) or making use of Monte Carlo techniques (Araya & Harding 1999; Araya-Góchez & Harding 2000; Schönherr et al. 2007).

Here we present new results (from the last five years) on the energy of the cyclotron resonance scattering feature E_cyc in the pulse averaged X-ray spectrum of Her X-1, combined with the historical long-term evolution. We present the first statistically significant evidence of a true long-term decrease in E_cyc and first evidence of a weak dependence of E_cyc on phase of the 35 d cycle. We speculate about the physics behind both the long-term decrease and the previously observed fast upward jump as being connected to changes in the configuration of the magnetic field structure at the responsible polar cap of the accreting neutron star. A preliminary report was published earlier by Staubert (2013).

2. Data base and method of analysis

Her X-1 is probably the best observed accreting binary X-ray pulsar. Its X-ray spectrum, including the CRSF, has been measured by many instruments since the discovery of the CRSF in 1976 (Trümper et al. 1978). The data base behind previously reported results are summarized in corresponding tables of the following publications: Gruber et al. (2001), Coburn et al. (2002), Staubert et al. (2007, 2009, 2013), Klochkov et al. (2008a, 2011), Vasco et al. (2011, 2013). Details about more recent observations (proposed by our group during the last five years) by RXTE, INTEGRAL, Suzaku, and NuSTAR are given in Tables 1 and 2. Of particular importance is the observation by NuSTAR in September 2012 which provided the most accurate value for the CRSF energy measured to date (Fürst et al. 2013). For the investigation of the long-term evolution of the cyclotron line energy, Main-On state observations at 35-day-phases <0.20 were used, in order to avoid interference with a dependence on 35 d phase (see Sect. 5). The spectral analysis was performed using the standard software appropriate for the respective satellites (see the publications cited above), with the addition of non-standard corrections based on our deeper analysis of the calibration of Imager on Board the INTEGRAL Satellite (IBIS, Ubertini et al. 2003), and the use of RECORN models (Rothschild et al. 2011) for RXTE. For the spectral model we have chosen the highecut¹ model which is based on a power law continuum with exponential cut-off, and the CRSF is modeled by a multiplicative absorption line with a Gaussian optical depth profile. Details of the fitting procedure can be found in the papers cited above.

3. Variation of the cyclotron line energy E_cyc

Variability in the energy of the CRSF in Her X-1 is found with respect to the following variables:

• Variation with X-ray luminosity (both on long and on shorttimescales).

• Variation with phase of the 1.24 s pulsation.

• Variation with phase of the 35 d precessional period.

• Variation with time, that is a true long-term decay.

The dependence on pulse phase, which is described in detail by Vasco et al. (2013), will only be discussed here to the extent necessary to understand the connection to the 35 d phase dependence.

Fig. 1 Positive correlation between the cyclotron line energy and the maximum X-ray flux of the corresponding 35 day cycle as measured by RXTE/ASM (see Fig. 2 of Staubert et al. 2007) with eight added points: INTEGRAL 2005 (Klochkov et al. 2008a), Suzaku of 2005 and 2006 (Enoto et al. 2008), RXTE 2009, INTEGRAL 2010, Suzaku 2010 and 2012, and NuSTAR 2012. The Suzaku points of 2005/2006 have been corrected upward by 2.8 keV, to account for the difference arising because the Lorentzian profile was used in the analysis by Enoto et al. (2008), while for all others the Gaussian profile was used. The blue rhombs are values observed until 2006, the red dots are from after 2006. The solid red line is a linear fit to data until 2006 with the original slope of 0.66 keV/(ASM cts/s), as found by Staubert et al. (2007). The dotted red line is the best fit to the data after 2006 with the slope fixed to the same value.

4. Variation of E_cyc with luminosity

For Her X-1, the dependence of the centroid energy of the phase averaged cyclotron line on X-ray flux was discovered by Staubert et al. (2007) while analyzing a uniform set of observations from RXTE. The original aim of the analysis at that time had been to investigate a possible decrease in the phase averaged cyclotron line energy with time during the first decade of RXTE observations. Instead, the dependence on X-ray flux was discovered and shown that the apparent decrease in the measured values of the line energy (see Sect. 6) was largely an artifact due to this flux dependence. The correlation was found to be positive, that is the cyclotron line energy E_cyc increases with increasing X-ray luminosity L_x.

Figure 1 reproduces the original correlation graph of Staubert et al. (2007) with new data points added (see Table 2). The first three new data points (INTEGRAL 05 and Suzaku 05/06) fit very well into the previous data set (and do not change the formal correlation – see the solid red line), but most of the values from 2006–2012 are significantly lower. As we will show below, it is these data which clearly establish a decrease in the cyclotron line energy with time. After 2006 the flux dependence is less obvious. However, the data points (except the one from INTEGRAL 2012) are consistent with the originally measured slope (0.66 keV/ASM-cts/s) with generally lower E_cyc values. The dotted red line is a fit through the data after 2006 with the same slope as the solid red line. We note that flux refers to the maximum Main-On flux as determined using the RXTE/ASM and/or the Swift/BAT monitoring data (since 2012 from BAT only); the conversion is: (2–10 keV ASM-cts/s) = 89 × (15–50 keV BAT-cts cm^-2 s^-1). The INTEGRAL 2012 point does clearly not follow this behavior, as will be more obvious below. We have invested a considerable effort to check the calibration of the INTEGRAL/ISGRI detector (INTEGRAL Soft Gamma-Ray Imager, Lebrun et al. 2003) for the time of observation and the data analysis procedure. The ISGRI response was closely examined by us for each of our Her X-1 observations. When necessary, the ARFs (Auxilliary Response Files) were checked (using the nearest Crab observations) and the energy scale was individually controlled by making use of observed instrumental background lines with known energy. Finally, spectra were generated using data from SPI (Spectrometer onboard INTEGRAL, Vedrenne et al. 2003): the resulting E_cyc values were always consistent with those of the ISGRI analysis. Since we have found no errors, we keep this point in our data base, but will exclude it from some of the analysis discussed below.

5. Variation of E_cyc with precessional phase

In order to investigate whether the cyclotron line energy has any dependence on phase of the 35 d precession, we had successfully scheduled a few Main-On observations at late 35 d phases. In addition to the full coverage of the Main-On of cycle No. 232 (a singular event!), we so far have four more measurements: from INTEGRAL in July/August 2007 and July 2011, as well as from Suzaku in September 2011 (all at 35 d phases 0.25, see Table 2) and from INTEGRAL in July 2005 at phase 0.24 (E_cyc = (37.3± 1.2) keV and (5.12 ± 0.37) ASM-cts/s). As we will show below, the cyclotron line energy changes with time. So, for a comparison of values measured at different times, they must be normalized to a common reference time. We are therefore making use of the results presented in Sect. 6 and compare E_cyc values which are normalized to the reference time of MJD 53500 (using the slope of fit 4 of Table 4): (37.11 ± 0.61) keV, (40.15 ± 0.81) keV, (36.92 ± 1.09) keV, and (37.36 ± 1.20) keV, respectively. Three out of these four values are indeed quite low in comparison to all other time normalized values, and gives an indication that E_cyc may indeed decrease at late 35 d phases. However, because of the lower fluxes at late 35 d phases, the uncertainties are fairly large for all of these measurements.

A more indirect, but perhaps more reliable method is the following: There are two well established observational facts with regard to pulse profiles and cyclotron line energies in Her X-1, both demonstrated in Fig. 2:

• 1.

  E_cyc varies strongly with pulse phase (by up to ~25%) (Voges et al. 1982; Soong et al. 1990; Vasco et al. 2013). The shape of the (E_cyc vs. pulse phase)-profile, is not dependent on 35 d phase (Vasco et al. 2013).

• 2.

  The main peak of the pulse profile moves to later pulse phases with increasing 35 d phase (Staubert et al. 2009). This is also true for the 30–45 keV profiles – the energy range which includes the CRSF (these profiles are shown in Fig. 2).

The combination of these two observational facts inevitably leads to a modulation of E_cyc with 35 d phase: with progressing 35 d phase, more and more photons are found at later pulse phases (in the main peak of the pulse) where the cyclotron line energy is decreasing. This means that the phase averaged cyclotron line energy must decrease with progressing 35 d phase (the above mentioned effect was first considered by Klochkov et al. 2008b). In order to quantitatively test this, we have performed a formal folding of 30–45 keV pulse profiles (for 10 different 35 d phases) with template (E_cyc vs. pulse phase) profiles for these 35 d phases. The template profiles were constructed by inter- and extrapolation of the four individual profiles at different 35 d phases as given in Figs. 3 and 4 of Vasco et al. (2013) (taking into account that both the maximum value and the peak-to-peak amplitude are slightly 35 d phase dependent). By folding the pulse profile with the corresponding E_cyc profile, the expected pulse phase averaged E_cyc value can be calculated.

Fig. 2 Mean dependence of cyclotron line energy on pulse phase for the Main-On of 35 d cycle 323, as observed by RXTE/PCA in 2002 November. The solid red line represents a best fit function (a combination of two cosine components). Normalized pulse profiles of the 30–45 keV range are shown for five different 35 d phases: 0.048 (black), 0.116 (red), 0.166 (green), 0.21 (blue) and 0.24 (purple). The main pulse is progressively moving to the right. The right hand scale is normalized flux (0–100) for the pulse profiles.

Fig. 3 Centroid pulse phase averaged cyclotron line energies at different 35 d phases of Main-On cycle 323. The data points with uncertainties are direct measurements for ten small integration intervals. The blue triangles connected by the solid blue line are values which are calculated by folding observed pulse profiles in the 30–45 keV range with template (Ecyc vs. pulse phase)-profiles for the same 35 d phases (see text for a detailed description).

In Fig. 3 these calculated values are shown as blue triangles (connected by the solid blue line): a slow increase up to phase ~0.19 is followed by a somewhat sharper decay. For comparison, we show the directly measured phase averaged E_cyc values (data points with uncertainties) for 13 small integration intervals covering the Main-On of cycle 323 (MJD 52 599/ Nov 2002). The directly measured values have relatively large uncertainties, but are overall consistent with the calculated modulation, and also with regard to the mean absolute value.

6. Variation of E_cyc with time – long-term variation

In Fig. 4 (an update of Fig. 1 of Staubert et al. 2007) we display observed values of the pulse phase averaged centroid cyclotron line energy as a function of time, covering the complete history of observations since the discovery of the line in 1976. We combine historical data, as taken from the compilation by Gruber et al. (2001, their Tables 2 and 3) for the time before the RXTE era, published values from observations with RXTE and INTEGRAL (Klochkov et al. 2006; Staubert et al. 2007; Klochkov et al. 2008a) and with Suzaku (Enoto et al. 2008), as well as recent values as given in Table 2 (see also Staubert 2013; Fürst et al. 2013). For the analysis of the long-term variation of E_cyc we exclude values with 35 d phases >0.20 in order to avoid contamination due to a possible third variable, the 35 d phase (a dependence, if any, is very weak for small phases; see Fig. 3).

Fig. 4 Centroid energy of the phase averaged cyclotron resonance line feature in Her X-1 since its discovery. Data from before 1997 were originally compiled by Gruber et al. (2001), where the original references can be found. The data after 1997 are from observations by RXTE, INTEGRAL (Klochkov et al. 2006, 2008a; Staubert et al. 2007) and Suzaku (Enoto et al. 2008), plus recent values as given in Table 2. Here only values measured at 35d phases <0.20 are shown.

Two features are apparent from Fig. 4: firstly, we confirm the apparent difference in the mean cyclotron line energy before and after 1991, first pointed out by Gruber et al. (2001). Taking the measured values of E_cyc and their stated uncertainties at face value, the mean cyclotron line energies ⟨ E_c ⟩ from all measurements before 1991 is 34.9 ± 0.3 keV, the corresponding value for all measurements between 1991 and 2006 is 40.3 ± 0.1 keV (40.2 ± 0.1 keV for RXTE results only, showing that the very high value measured by BATSE is not decisive). However, a comparison of measurements from different instruments is difficult because of systematic uncertainties due to calibration and analysis techniques. Nevertheless, we believe that the large difference of ~5 keV between the mean values and the good internal consistency within the two groups (5 different instruments before 1991 and four after 1991) most likely indicate real physics.

As already mentioned in Sect. 4, the first observations with RXTE in 1996 and 1997 showed lower E_cyc values than those found from CGRO/BATSE and Beppo/SAX, leading to the idea of a possible decay with time. This idea had then served successfully as an important argument to ask for more observations of Her X-1. In a series of RXTE observations until 2005 the apparent decrease seemed to continue until this date. At this time we were determined to publish a paper claiming evidence of a decay of the phase averaged cyclotron line energy E_cyc. However, working with a uniform set of RXTE data between 1996 and 2005, we discovered that there was a dependence of E_cyc on X-ray flux (Staubert et al. 2007), degrading the apparent decrease with time largely to an artifact: nature seemed to have conspired such that later measurements were (on average) taken when the flux happened to be low (Her X-1 is known for varying its flux within a factor of two, on timescales of a few 35 d cycles). When the cyclotron line energy was normalized to a common flux value, the time dependence vanished.

In the following we will demonstrate in a systematic way that today we have clear evidence for a reduction in the phase averaged cyclotron line energy E_cyc with time over the last 20 years. Both dependencies – on flux and on time – seem to be always present (they may, however, change their relative importance with time). Using a procedure of fitting with two variables simultaneously, the two dependencies can be separated and the formal correlation minimized. With the inclusion of new measurements (2005–2012), we are now able to present the first statistically significant evidence of a true long-term decay of the phase averaged cyclotron line energy.

Fig. 5 Upper panel, left scale: Her X-1 pulse phase averaged cyclotron line energies Ecyc normalized to a reference ASM count rate of 6.8 cts/s using a flux dependence of 0.54 keV/ASM-cts/s. A break in mean Ecyc after MJD 54 000 (>2006) is apparent. Lower panel, right scale: 2–10 keV X-ray flux (of 35 d maximum) from monitoring by RXTE/ASM (from Swift/BAT after MJD 56 000). Blue data points with uncertainties are those fluxes which are used to correlate with Ecyc, the green curve connects measurements of each 35 d cycle.

6.1. Normalizing E_cyc using the originally discovered flux dependence

Before introducing fits with two variables (flux and time), we repeat here the procedure applied in Staubert (2013), that neglects any time dependence and normalizes E_cyc by flux only, using the originally determined linear dependence with a slope of 0.66 keV/(ASM-cts/s) (Staubert et al. 2007). The normalization is done to the common reference flux of 6.80 ASM-cts/s. The results were shown in Fig. 2 (right) of Staubert (2013) (and, slightly updated, in Fig. 13 of Fürst et al. 2013), showing the flux normalized cyclotron line energy as a function of time (1996–2012). These figures demonstrate that there is a significant decrease in E_cyc at least after 2006 (>MJD 54 500). Based on a repeated analysis we give the following additional information.

• (1)

  The normalized E_cyc values between 1996 and 2006 (MJD 50 000–54 000) are consistent with a constant (based on χ²/d.o.f. = 0.77), supporting the neglect of a time dependence for this decade. The weighted mean is (40.06 ± 0.09) keV. However, even during this time period there is a slight downward trend with a slope of (− 2.5 ± 1.4) × 10^-4 keV d^-1. In order to test whether this downward trend is really a true time dependence and not an artifact due to the neglect of e.g., a more complicated flux dependence, a quadratic term was added when doing the flux normalization: this does not improve the fit and does not remove the downward trend.

• (2)

  The weighted mean of the normalized E_cyc values for 2006–2012 is (37.69 ± 0.10) keV (with a high χ²/d.o.f. of 6.6, mainly due to the INTEGRAL 2012 point).

• (3)

  The difference between these two mean values is highly significant (>17 standard deviations), demonstrating the decrease in E_cyc with time.

6.2. Normalizing E_cyc using fits to the 1996–2006 data with two variables

From Fig. 1, it is already evident that most of the E_cyc values measured after 2006 are significantly lower than those from before. In order to separate the dependence on time and the dependence on X-ray flux, we have performed fits to the 1996–2006 data with two variables – X-ray flux and time. We use the function (1)with F being the X-ray flux (the maximum flux of the respective 35 d cycle) in units of ASM-cts/s, as observed by RXTE/ASM (and/or Swift/BAT), with F₀ = 6.80 ASM-cts/s, and T being time in MJD with T₀ = 53 000 – the relationship between ASM and BAT is the following: (2–10 keV ASM-cts/s) = 89 * (15–50 keV BAT-cts cm^-2 s^-1).

In Table 3 we summarize results of fits to the 1996–2006 data set with successive numbers of free parameters. Assuming no dependence at all (a = b = 0.0), as well as dependence only on time (a = 0) leads to unacceptable fits with χ²/d.o.f. of 3.4 and 2.9, respectively. Allowing simultaneous dependence on flux and on time leads to a good fit (χ²/d.o.f. = 0.6), representing a description of these data with a satisfactory separation of flux and time dependence. It is again verified, that the flux dependence is the dominating effect, as clearly seen when the F-test is applied².

Allowing a simultaneous flux and time dependence, reduces the flux dependence slightly (as compared to neglecting the time dependence). If we now use the flux dependence of 0.54 keV/ (ASM-cts/s) for the normalization of the complete data set (1996–2012), we arrive at Fig. 5. We note the following features:

• 1)

  The mean cyclotron line energies before andafter 54 500 (~2007) are 40.1 ± 0.1 and 37.6 ± 0.1, respectively. The difference is significant to >17 standard deviations (similar to the result of Sect. 6.1, where the flux normalization was done with the original slope of 0.66 keV/(ASM-cts/s)).

• 2)

  As noted before, there is a small downward trend between 1996 and 2006. Again, an added quadratic term for the flux normalization is not significant.

• 3)

  Because of the lack of measurements in 2007 and 2008 (at 35 d phases <0.20), we cannot distinguish between a fairly abrupt drop within this period and a smooth change from the 2004–2006 period to the 2009–2012 period (possibly with a somewhat stronger decay than between 1996 and 2006).

• 4)

  The lower panel of Fig. 5 shows the 2–10 keV X-ray flux (the maximum flux of each 35 d cycle) from monitoring by RXTE/ASM and Swift/BAT. The typical variation of the flux by a factor of ~2 (as also evident from Fig. 1) is apparent. The mean flux, averaged over several 35 d cycles, is constant.

Fig. 6 Her X-1 pulse phase averaged cyclotron line energies Ecyc normalized to a reference ASM count rate of 6.8 cts/s using a flux dependence of 0.44 keV/ASM-cts/s. The data are now consistent with a linear decline of Ecyc with time with a slope of −7.22 × 10-4 keV d-1.

6.3. Normalizing E_cyc using fits to the 1996–2012 data with two variables

We finally turn to fits with two variables (flux and time) to the data set of 1996–2012. For these fits we exclude the data point measured in April 2012 by INTEGRAL since, as noted before, this value is not consistent with measurements by both Suzaku and NuSTAR about six months later.

As discussed in Sect. 4, we have invested a strong effort to check the calibration of the INTEGRAL/IBIS detector for the time of observation and the subsequent data analysis procedure. We have not found any errors, so we keep this point in our data base (and show it in all plots), but exclude it from the fits to be discussed now (but even if this data point is included, the general conclusion about the long-term decrease in the cyclotron line energy is not changed). As before, the bilinear function (formula (1)) is applied. The results of these fits with increasing numbers of free parameters are summarized in Table 4, the fits are numbered 1 through 4. Using the flux dependence found in the final (simultaneous) fit (fit 4) for normalizing the observed E_cyc values to the reference flux of 6.8 (ASM-cts/s), the remaining linear time dependence is shown in Fig. 6. The bilinear fit number 4 is acceptable with a χ² = 20.4 for 18 d.o.f.

From Fig. 6 we see that the INTEGRAL12 point is far (~6σ) from the best fit (the dashed line). If we now repeat the bilinear fit including this point, the linear time dependence is not changed significantly (from (−7.22 ± 0.39)×10^-4 to (−6.91 ± 0.39)×10^-4 keV d^-1), which means that the conclusions about the decay of E_cyc with time are not changed.

We note that a further improvement of the fit can be achieved by introducing a quadratic term in the time dependence (again for INTEGRAL12 excluded): the flux dependence is unchanged, the linear time term is now (− 6.59 ± 0.49) × 10^-4 keV d^-1 and the quadratic time term is (− 6.88 ± 3.18) × 10^-8 keV d^-2 (both with a reference time MJD 53 500), and a χ² = 15.7 for 17 d.o.f. However, the improvement is only marginally significant with an F-value of 5.1, corresponding to a probability of 3.7% for an improvement by chance.

Referring to Table 4 and Fig. 6 we note the following results:

• 1)

  Clearly, fit 1 is not acceptable (χ² = 430 for 20 d.o.f.): E_cyc is not constant.

• 2)

  Introducing a linear flux dependence (neglecting any time dependence (fit 2)), improves the χ² significantly and finds a flux dependence of 0.74 ± 0.09 keV/(ASM-cts/s), which is (within uncertainties) consistent with the value found for the previous fit over the shorter time range (Table 3). However, a flux dependence alone is not sufficient. In addition, when a quadratic term for the flux normalization is introduced, assuming no time dependence, the fit is not acceptable (χ² = 313, compare the corresponding values of Table 4).

• 3)

  Introducing a linear time dependence (neglecting any flux dependence (fit 3)), improves the χ² dramatically. The F-value corresponds to a formal probability of <10^-10 for the improvement to be just by chance.

• 4)

  Adding now the linear flux dependence (fit 4), χ² is further reduced significantly, meaning that the flux dependence is definitely needed in addition to the time dependence (the F-value corresponds to a chance probability of ~5 × 10^-5). The slope describing the flux dependence is somewhat reduced as compared to fit 2, but is now very close to the corresponding dependence found in fitting the <2006 data (see Table 3).

• 5)

  We note that now the time dependence is the more dominant variation, while for the <2006 data set it was the flux dependence. This is consistent with a ~4 keV reduction in E_cyc over the covered time range (1996–2012), while the increase is only ~3 keV over the flux range provided by nature (a factor of ~2 in flux). A reduction in importance of the flux dependence, when adding the data after 2006, can already be expected from Fig. 1: the >2006 data are consistent with no flux dependence at all (even though they are formally consistent with the same slope as the <2006 data).

• 6)

  An even better fit can be achieved by introducing a quadratic term in the time dependence (E_cyc/dt² = (−6.88 ± 3.18) × 10^-8 keV d^-2), albeit marginally significant (3.7% chance probability). The negative second derivative means that the decrease in E_cyc accelerates with time. We note, that one might rather expect the opposite, e.g., a kind of exponential decay.

• 7)

  It is not understood why the measurement of INTEGRAL in 2012 is so different from the nearby (<six months) observations by Suzaku and NuSTAR. We find no errors in our analysis. However, this one data point does not in any way change the general conclusion.

7. Summary of observational results

Here we summarize our observational results:

• 1.

  The main result of this research is that we finally can establish along-term decay of the pulse phase averaged cyclotron lineenergy E_cyc over time. The time range covered here is 16 years (1996 to 2012). The reduction is highly significant and can be well described by a linear decay with a change by (7.2 ± 0.4) × 10^-4 keV d^-1 (0.26 keV yr^-1 or 4.2 keV over 16 yrs). The data are, however, also consistent with two other models: First (see Fig. 5), a somewhat slower decay (~3 × 10^-4 keV d^-1) until 2006 and then a more sudden drop between 2006 and 2009, with a possible constant value (~37.7 keV) thereafter. Alternatively, even an acceleration of the decay of E_cyc over time is possible, since the fit with a quadratic term in the time dependence is formally the best one.

• 2.

  The flux dependence of E_cyc is confirmed with a value between 0.44 and 0.54 keV/ASM-cts/s (with a typical uncertainty of 0.1), corresponding to a ~5% increase for a factor of two increase in flux. This value is slightly lower than the value (0.66±0.10) from the original discovery (Staubert et al. 2007) (at that time, a time dependence was neglected). It is not excluded that there is a variation with time of this flux dependence, but this is difficult to judge since the fluxes observed after 2006 only cover a very small range.

• 3.

  One observation, the one from INTEGRAL in April 2012, does not fit into the overall picture and is inconsistent with values measured by Suzaku and NuSTAR six months later. Since we find no errors in our analysis, we can only speculate that it is due to real physics – that is a fluctuation on a timescale of a few months.

• 4.

  From the analysis of pulse profiles and pulse phase resolved spectroscopy (see Sect. 5), both as a function of 35 d phase, we find evidence of a variation of E_cyc with this precessional phase: while there is very little (if any) variation until 35 d phase 0.2, there may be a decrease at later phases. If it is indeed so, that the pulse phase dependence of E_cyc does not change with 35 d phase (Vasco et al. 2013), then this effect is inevitable. The fully covered Main-On of cycle 232 (Fig. 3) and the few dedicated observations at late 35 d phases support this finding. This new result adds another piece to the puzzle on the question about the physical nature of the 35 d modulation – precession of the accretion disk (plus?) precession of the neutron star? – and the mechanism of generating the varying pulse profiles and the varying spectra – continua and CRSF (see Sect. 8 and discussions in Staubert et al. 2009, 2013; Vasco et al. 2013).

8. Discussion

8.1. Dependence of E_cyc on luminosity

A negative correlation between the pulse phase averaged cyclotron line energy and the X-ray luminosity (a decrease in E_cyc with increasing L_x), had first been noted by Mihara (1995) in observations of a few high luminosity transient sources (4U 0115+63, Cep X-4, and V 0332+53) by Ginga. This negative correlation was associated with the high accretion rate during the X-ray outbursts, as due to a change in height of the shock (and emission) region above the surface of the neutron star with changing mass accretion rate, Ṁ. In the model of Burnard et al. (1991), the height of the polar accretion structure is tied to Ṁ. From this model one expects that an increase in accretion rate leads to an increase in the height of the scattering region above the neutron star surface, and therefore to a decrease in magnetic field strength and hence a decrease in E_cyc. During the 2004/2005 outburst of V 0332+53 a clear anti-correlation of the line position with X-ray flux was observed (Tsygankov et al. 2006). A similar behavior was observed in outbursts of 4U 0115+63 in March/April 1999 and Sep/Oct 2004: both Nakajima et al. (2006) and Tsygankov et al. (2007) had found a general anti-correlation between E_cyc and luminosity. However, Müller et al. (2013b), analyzing data of a different outburst of this source in March/April 2008, observed by RXTE and INTEGRAL, have found that the negative correlation for the fundamental cyclotron line is likely an artifact due to correlations between continuum and line parameters when using the NPEX continuum model.

The first positive correlation was discovered by Staubert et al. (2007) in Her X-1, and secured by a reanalysis of the same RXTE data by Vasco et al. (2011), using the bolometric X-ray flux as reference. This analysis confirmed that the originally used 2–10 keV flux is a good measure of the bolometric luminosity. While the above discussed analysis tests the correlated variability of E_cyc and L_x on long timescales (35 d and longer), the pulse-amplitude resolved analysis of Klochkov et al. (2011) does so on short timescales (down to the pulse period of 1.24 s). Selecting pulses with amplitudes in certain ranges and producing mean spectra for each pulse amplitude range, showed that the cyclotron line energy scales positively with the mean pulse amplitude. In addition, it was found that the photon index Γ of the underlying power law continuum scales negatively with the pulse amplitude (the absolute value of Γ gets smaller, that is the spectrum flattens). The same behavior was seen in data of the transient A 0535+26. A recent pulse phase resolved analysis of A 0535+26 observations by RXTE and INTEGRAL showed that data of one of the two peaks (of the double peak pulse profile) displays the same trend while data of the other peak do not (Müller et al. 2012, 2013a). Applying the same pulse-amplitude resolved technique to data of V 0332+53 and 4U 0115+63, Klochkov et al. (2011) found the same behavior as originally detected in data sets that were selected on much longer timescales: E_cyc decreases and Γ increases with increasing L_x. Finally, we mention that a positive correlation of E_cyc with L_x was also found in two more X-ray binary pulsars: in GX 304-1 (Yamamoto et al. 2011; Klochkov et al. 2012) and in NuSTAR observations of Vela X-1 (Fürst et al. 2014). We note, that the still small group of four objects with a positive E_cyc/L_x correlation now outnumbers the group of secure sources with the originally discovered opposite behavior.

Our current understanding of the physics behind these correlations assumes that we can distinguish between two accretion regimes in the accretion column above the polar cap of the neutron star: super- and sub-Eddington accretion. The former is responsible for the first detected negative correlation in high luminosity outbursts of transient X-ray sources (the reference source being V 0332+53): in this case the deceleration of the accreted material is provided by radiation pressure, such that with increasing accretion rate Ṁ, the shock and the scattering region move to larger height above the surface of the neutron star and consequently to weaker B-field (Burnard et al. 1991). Sub-Eddington accretion, on the other hand, leads to the opposite behavior. In this regime the deceleration of accreted material is predominantly through Coulomb interactions and an increase in Ṁ leads to an increase in electron density (due to an increase of the combined hydrostatic and dynamical pressure) resulting in a squeezing of the decelerating plasma layer to smaller height and stronger B-field (Staubert et al. 2007). More detailed physical considerations have recently been presented by Becker et al. (2012). The persistent sources Her X-1 and Vela X-1 are clearly sub-Eddington sources.

Despite the above discussed doubts about the reality of the negative correlation of the energy E_cyc of the fundamental CRSF with L_x in 4U 0115+63 (Müller et al. 2013b), and keeping in mind that Klochkov et al. (2011) had confirmed the correlation using the pulse-amplitude resolved technique (not using the NPEX function), we would like to note here that we are intrigued by the following plots about 4U 0115+63: Fig. 8 of Nakajima et al. (2006) and Figs. 11 and 12 of Tsygankov et al. (2007). In going to the lowest luminosities, there is an indication for a leveling-off or even a reduction in E_cyc. Do we possibly see here the transition between the two accretion regimes?

8.2. The long-term decay of E_cyc

With regard to the physical interpretation of the now observed long-term decrease in the cyclotron line energy, we speculate that it could be connected to either unknown effects in the neutron star and its magnetic field, to a geometric displacement of the cyclotron resonant scattering region in the dipole field or to a true physical change in the magnetic field configuration at the polar cap, which evolves due to continued accretion. Apparently, the magnetic field strength at the place of the resonant scattering of photons trying to escape from the accretion mound surface must have changed with time. Putting internal neutron star physics aside, we suggest that it reflects a local phenomenon in the accretion mound: either a geometric displacement of the emission region or a change in the local field configuration, rather than a change in the strength of the underlying global dipole field (here a minimum timescale of a million years is estimated from population studies of rotation-powered pulsars; Bhattacharya et al. 1992; Geppert & Urpin 1994). Our observed timescale, a few tens of years, is extremely short.

The whole issue of accretion onto highly magnetized neutron stars in binary X-ray sources is very complex. Ideas or models with potential relevance to our observations attempting to understand the magnetic field configuration in accreting neutron stars and its evolution over extended periods of continued accretion, can be found in e.g.: Hameury et al. (1983), Konar & Bhattacharya (1997), Brown & Bildsten (1998), Cheng & Zhang (1998), Litwin et al. (2001), Cumming et al. (2001), Melatos & Phinney (2001), Choudhuri & Konar (2002), Payne & Melatos (2004), Payne & Melatos (2007), Wette et al. (2010), Mukherjee & Bhattacharya (2012), Mukherjee et al. (2013a), Mukherjee et al. (2013b). However, as far as we can see, none of them gives the complete picture. Most calculations deal with static solutions that are found under special boundary conditions.

Since the main purpose of this contribution is to report on the discovery of a new observational phenomenon, we refrain from going into any details regarding interpretations of existing models. Instead, we only mention a few areas which we think could have some connection with the observed facts and which may be worthwhile to be explored. Our hope is that the new observational results presented here may boost the motivation for further theoretical studies.

We start by asking whether the observed decrease in E_cyc with time could be a simple movement of the resonant scattering region to a larger distance from the neutron star surface, where the field strength is lower. This would be similar to the decrease in E_cyc during outbursts in high luminosity transients, except that (being in the sub-Eddington regime of accretion) we would not think of the shock region to rise, but rather the total height of the accretion mound may slowly increase with time, such that also the resonant scattering region is displaced to a higher position. For a dipole field with an r^-3 dependence of field strength, the observed ~5 keV reduction in E_cyc from 1992 to 2012 (0.25 keV per year) would correspond to a change in height of ~400 m (starting from the surface itself). The question here is, whether continued accretion really leads to a growth of the accretion mound with time – both in terms of geometrical height and of total mass.

With similar uncertainty, it can be asked whether the accreted material could drag the central field lines radially out, possibly enlarging the total hotspot area and thereby diluting the effective field strength in the region where the resonant scattering takes place (Cheng & Zhang 1998; Zhang & Kojima 2006). Or, whether screening or burial of the magnetic field at the polar caps is possible (Brown & Bildsten 1998; Payne & Melatos 2004; Litwin et al. 2001; Payne & Melatos 2007). It needs to be investigated, how much mass could eventually be stored in the magnetically confined mountains, whether matter is continuously leaking out to larger areas of the neutron star surface (due to plasma pressure exceeding the magnetic pressure) and on what timescales an observational effect can be expected.

Finally, the question of Ohmic dissipation and diffusion of the magnetic field may play a role and physical processes either in the accretion mound or in or below the surface of the neutron star (like hydrodynamic flows) could bury or reduce the surface field (Choudhuri & Konar 2002; Patruno 2012). One would need to investigate whether physical parameters like the characteristic length scale and the relevant conductivity σ for either the crust or the plasma in the accretion mound could be of the right order of magnitude to allow the magnetic diffusion timescale τ = 4πR²σ/c² (Cumming et al. 2001; Ho 2011) to be compatible with the timescale of a few tens of years – as observed for the decrease in the local polar field strength in Her X-1. If magnetic diffusion is indeed relevant, we note that the necessary small length scales and relatively low conductivities would argue for local physics in the hot plasma of the accretion mound, the structure of which is most likely complex because of contributions from higher-order multipoles.

We finally speculate on a possible cyclic behavior of E_cyc on timescales of a few tens to hundreds of years. Could it be that the fast rise of the observed E_cyc values after 1991 (see Fig. 4) represents a special event in which the magnetic field in the accretion mound has rearranged itself as a result of a sudden radial outflow of material? In models by Brown & Bildsten (1998); Litwin et al. (2001); Payne & Melatos (2004); Payne & Melatos (2007); Mukherjee & Bhattacharya (2012); Mukherjee et al. (2013a) the field configuration is shown to change considerably with increased material, leading to a ballooning of the field configuration with diluted field in the symmetry center and increased density of field lines at the circumference of the base of the mound. The estimates of how much mass could be confined by the field vary substantially between the different models. It remains unclear, how important continuous leaking through the outer magnetic boundary may be and what the timescales for semi-catastrophic events might be, in which the field would release (on a short timescale) a substantial fraction of stored material to larger areas of the neutron star surface. For Her X-1 this scenario could mean that we are now in a phase of continuous build-up of the accretion mound with the mass (and the height?) of the mound growing and the observed cyclotron line energy continuously decreasing until another event like the one around 1991 happens again. The mean E_cyc value measured before 1991 of ~35 keV may represent a bottom value. So, when the current decay continues steadily, one may expect another event of a rather fast increase in E_cyc.

In conclusion, we like to urge both observers and model builders to continue to accumulate more observational data as well as more understanding of the physics responsible for the various observed properties of Her X-1 and other objects of similar nature. For model builders a challenge would be to work towards dynamical computations that might eventually lead to self-consistent solutions of the structure and evolution of magnetized accretion mounds of accreting neutron stars with only a few input parameters.

Acknowledgments

This paper is to a large part based on observational data taken by the NASA satellite Rossi X-Ray Timing Explorer (RXTE). We like to acknowledge the dedication of all people who have contributed to the great success of this mission. In the same way, we thank the teams of ESA’s INTErnational Gamma-Ray Astrophysics Laboratory (INTEGRAL), JAXA’s Suzaku and NASA’s Nuclear Spectroscopic Telescope Array (NuSTAR). This work was supported by DFG through grants Sta 173/31-1, 2 and 436 RUS 113/717 and RFBR grants RFFI-NNIO-03-02-04003 and 06-02-16025. The work of K.P. and N.Sh. was also partially supported by RBFR grants 12-02-00186 and 14-02-00657. D.K. is indebted to the Carl Zeiss Stiftung for support. We thankfully acknowledge very useful discussions about the possible physical meaning of the observed effects with D. Bhattacharya, K. Kokkotas, K. Glampedakis and J. Trümper. Finally we thank the anonymous referee for important questions and suggestions.

References

All Tables

All Figures

Fig. 1 Positive correlation between the cyclotron line energy and the maximum X-ray flux of the corresponding 35 day cycle as measured by RXTE/ASM (see Fig. 2 of Staubert et al. 2007) with eight added points: INTEGRAL 2005 (Klochkov et al. 2008a), Suzaku of 2005 and 2006 (Enoto et al. 2008), RXTE 2009, INTEGRAL 2010, Suzaku 2010 and 2012, and NuSTAR 2012. The Suzaku points of 2005/2006 have been corrected upward by 2.8 keV, to account for the difference arising because the Lorentzian profile was used in the analysis by Enoto et al. (2008), while for all others the Gaussian profile was used. The blue rhombs are values observed until 2006, the red dots are from after 2006. The solid red line is a linear fit to data until 2006 with the original slope of 0.66 keV/(ASM cts/s), as found by Staubert et al. (2007). The dotted red line is the best fit to the data after 2006 with the slope fixed to the same value.

In the text

Fig. 2 Mean dependence of cyclotron line energy on pulse phase for the Main-On of 35 d cycle 323, as observed by RXTE/PCA in 2002 November. The solid red line represents a best fit function (a combination of two cosine components). Normalized pulse profiles of the 30–45 keV range are shown for five different 35 d phases: 0.048 (black), 0.116 (red), 0.166 (green), 0.21 (blue) and 0.24 (purple). The main pulse is progressively moving to the right. The right hand scale is normalized flux (0–100) for the pulse profiles.

In the text

Fig. 3 Centroid pulse phase averaged cyclotron line energies at different 35 d phases of Main-On cycle 323. The data points with uncertainties are direct measurements for ten small integration intervals. The blue triangles connected by the solid blue line are values which are calculated by folding observed pulse profiles in the 30–45 keV range with template (Ecyc vs. pulse phase)-profiles for the same 35 d phases (see text for a detailed description).

In the text

Fig. 4 Centroid energy of the phase averaged cyclotron resonance line feature in Her X-1 since its discovery. Data from before 1997 were originally compiled by Gruber et al. (2001), where the original references can be found. The data after 1997 are from observations by RXTE, INTEGRAL (Klochkov et al. 2006, 2008a; Staubert et al. 2007) and Suzaku (Enoto et al. 2008), plus recent values as given in Table 2. Here only values measured at 35d phases <0.20 are shown.

In the text

Fig. 5 Upper panel, left scale: Her X-1 pulse phase averaged cyclotron line energies Ecyc normalized to a reference ASM count rate of 6.8 cts/s using a flux dependence of 0.54 keV/ASM-cts/s. A break in mean Ecyc after MJD 54 000 (>2006) is apparent. Lower panel, right scale: 2–10 keV X-ray flux (of 35 d maximum) from monitoring by RXTE/ASM (from Swift/BAT after MJD 56 000). Blue data points with uncertainties are those fluxes which are used to correlate with Ecyc, the green curve connects measurements of each 35 d cycle.

In the text

	Fig. 6 Her X-1 pulse phase averaged cyclotron line energies Ecyc normalized to a reference ASM count rate of 6.8 cts/s using a flux dependence of 0.44 keV/ASM-cts/s. The data are now consistent with a linear decline of Ecyc with time with a slope of −7.22 × 10-4 keV d-1.
In the text