{'fraction_non_alphanumeric': 0.03603272359591874, 'fraction_numerical': 0.04015074914973803, 'mean_word_length': 4.3339870562855465, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
GRO J1008-57: an (almost) predictable transient X-ray binary

A study of archival RXTE, Swift, and Suzaku pointed observations of the transient high mass X-ray binary GRO J1008-57 is presented. A new orbital ephemeris based on pulse arrival timing shows times of maximum luminosities during outbursts of GRO J1008-57 to be close to periastron at orbital phase -0.03. This makes the source one of a few for which outburst dates can be predicted with very high precision. Spectra of the source in 2005, 2007, and 2011 can be well described by a simple power law with high energy cutoff and an additional black body at lower energies. The photon index of the power law and the black body flux only depend on the 15-50 keV source flux. No apparent hysteresis effects are seen. These correlations allow to predict the evolution of the pulsar's X-ray spectral shape over all outbursts as a function of just one parameter, the source's flux. If modified by an additional soft component, this prediction even holds during GRO J1008-57's 2012 type II outburst.


Introduction
GRO J1008−57 was discovered by CGRO during an X-ray outburst on 1993 July 14 (Wilson et al. 1994;Stollberg et al. 1993). It is a transient neutron star with a Be-star companion of type B0e (Coe et al. 1994) at a distance of 5.8 kpc (Riquelme et al. 2012). In such high-mass X-ray binaries (HMXBs) the compact object is on a wide eccentric orbit around its companion, which itself features a circumstellar disc due to its fast rotation. Periastron passages of the neutron star lead to accretion from the donor's decretion disk and regular X-ray outbursts. Based on RXTE-ASM measurements of the X-ray light curve, Levine & Corbet (2006, see also Levine et al. 2011) derived an orbital period of 248.9(5) d 1 . Since individual X-ray outbursts are relatively short, with durations of around 14 d, determining the parameters of the binary orbit is challenging. Coe et al. (2007) have derived an orbital solution (see Table 2) based on pulse period folding of the 93.6 s X-ray period (Stollberg et al. 1993).
The hard X-ray spectrum of GRO J1008−57 was first studied by Grove et al. (1995) and Shrader et al. (1999). Above 20 keV the spectrum can be well described by an exponentially cutoff power law (Grove et al. 1995;Shrader et al. 1999;Coe et al. 2007). Below 20 keV, Suzaku observations show a com- 1 The error bars are given in units of the last digit shown. plex spectrum with Fe line fluorescent emission and a power law continuum (Naik et al. 2011).
A slight deviation from the continuum at 88 keV in CGRO-OSSE observations was interpreted by Shrader et al. (1999) as a possible cyclotron line feature, i.e., as a feature caused by inelastic scattering of photons off electrons quantized in the strong magnetic field at the neutron star's poles (Schönherr et al. 2007;Caballero & Wilms 2012;Pottschmidt et al. 2012, and references therein). Even if the 88 keV feature is the second harmonic cyclotron line as argued by Shrader et al. (1999), GRO J1008−57 would be one of the most strongly magnetized neutron star in an accreting system to date.
In this paper an analysis of archival observations of GRO J1008−57 with RXTE, Suzaku, and Swift is presented. Section 2 presents the data analysis strategy. In Sect. 3 an updated orbital solution based on data from outbursts in 2005 and 2007 is presented. The X-ray spectrum of GRO J1008−57 is discussed in Sect. 4 and the spectroscopic results are applied and compared in Sect. 5. The paper closes with a summary and conclusions in Sect. 7.

Observations & Data Reduction
Data from five outbursts of GRO J1008−57 are analyzed. The earliest outburst studied here occurred in 2005 February and   Projected orbital plane of the binary. The reference frame is centered on the optical companion, whose size is set to 7 R (Coe et al. 2007). The dashed circle sketches a circumstellar disk of 72 R as proposed by (Coe et al. 2007). Orbital phases of the observations from 2005, 2007, and 2011 April are marked by arrows. The colors are the same as in Fig. 1. Grey arrows represent RXTE observations done in quiescence (Obsids 20132-* and 20123-*). Orbital phases are calculated using the ephemeris given in Table 2.
The binary orbital phases of the 2005, 2007, and 2011 April observations are shown in Fig. 2. Outbursts occurred around the periastron of the orbit, where mass accretion from the circumstellar disk is possible.
In the following the instruments used in the analysis of these outbursts are described. All spectral and timing modeling was performed with the Interactive Spectral Interpretation System (ISIS, Houck & Denicola 2000, version 1.6.2-7). Unless otherwise stated, all error bars are given at the 90% level for each parameter of interest.

RXTE
The Proportional Counter Array (PCA, Jahoda et al. 2006) was one of two pointed instruments onboard the Rossi X-ray Timing Explorer (RXTE) and was sensitive for X-rays between 2 and 90 keV. The PCA consisted of five Proportional Counter Units (PCUs), collimated to a field of view of ∼1 • . Unless stated otherwise, only data from PCU2 were used, which was the best  (Jahoda et al. 2006). Since the source can be weak away from the peak of the outbursts, only data from the top layer of PCU2 were used. Data were reduced using heasoft 6.11 and using standard data reduction pipelines (Wilms et al. 2006). Light curves for the timing analysis were extracted using barycentered GoodXenon data with 1 s time resolution. For spectral analysis, 4.5 keV-50 keV PCA data were used to avoid calibration issues around the Xenon L-edge. Data from 10-20 keV, 20-30 keV, 30-40 keV, and >40 keV were rebinned by a factor of 2, 4, 6, and 8, respectively. In addition, systematic uncertainties of 0.5% were added to the spectra (Jahoda et al. 2006). To compensate for the few percent uncertainty of the PCA's background model, the background was scaled with a multiplicative factor.
Hard X-ray RXTE data were obtained using the High Energy X-ray Timing Experiment (HEXTE, Rothschild et al. 1998). HEXTE consisted of two clusters A and B of four NaI/CsI scintillation detectors, alternating between the source and background positions ("rocking"). In 2006 October, the rocking of cluster A was turned off and the detector was fixed in on-source position due to an anomaly in the rocking mechanism. As no reliable background estimates are available, cluster A was not used for spectral analysis of the post 2006 data. The rocking of cluster B was switched off in 2010 April, thus no HEXTE data was used for analysis during the 2011 April outburst. 20-100 keV HEXTE data were used for spectral analysis and rebinned by a factor of 2, 3, 4, and 10 in the energy intervals 20-30 keV, 30-45 keV, 45-60 keV, and >60 keV, respectively.
The analysis of the GRO J1008−57 data is complicated due to the low Galactic latitude of the source (b = −1. • 827). At such low latitudes, the cumulative effect of the Galactic ridge emission in the ∼1 • radius field of view of the PCA is strong enough to lead to visible features in the X-ray spectrum. When comparing simultaneous PCA measurements with those from imaging instruments, these features mainly show up as excess emission in the Fe Kα band around 6.4 keV and a slight difference in the continuum shape (e.g., Müller et al. 2012). The Galactic ridge emission is believed to originate from many, unresolved X-ray binaries in the field of view (Revnivtsev et al. 2009) and can be empirically modelled by the sum of a bremsstrahlung continuum and iron emission lines from neutral, helium-like, and hydrogenic iron (Ebisawa et al. 2007;Yamauchi et al. 2009). Due to the energy resolution of the PCA, the three lines merge into a blended emission line.
The Galactic ridge emission shows a strong spatial variability. Fortunately, GRO J1008−57 was observed by RXTE in 1996 and 1997 during a quiescent state (proposal IDs 20132 and 20123). The PCA light curves from these observations do not show any evidence for pulsed emission between 91 and 96 s. Figure 3 shows the 181 ks spectrum accumulated from these observations. In the spectral analysis, this component was modeled as the sum of a bremsstrahlung continuum and three Fe K lines and is only present in the PCA data (the ridge emission is too soft to influence HEXTE).

Swift
For the analysis of the Swift observations data from its Xray Telescope's CCD detector, XRT, were used (Gehrels et al. 2004). The 1.5-9 keV spectra were extracted in the "Windowed Timing mode", a fast mode that prevents pileup for sources up to 600 mCrab. Only the observations after the 2011 December outburst and obs. 00031030013 before the giant outburst of 2012 October/November were extracted in the "Photon Counting mode". Source spectra were obtained from a 1 radius circle centered on the source, or an annular circular region if pileup was present. Background spectra were accumulated from two circular regions off the source position with radii up to ∼1. 5. Unless stated otherwise, spectral channels are added up until a minimum signal-to-noise ratio of 20 is achieved. The light curve was extracted with 1 s time resolution.

Suzaku
Similar to Swift-XRT, the X-ray Imaging Spectrometer (XIS) mounted on the Suzaku satellite uses Wolter telescopes to focus the incident X-rays ). After correcting the spacecraft attitude (Nowak et al. 2011) "Normal Clock Mode" 1-9 keV data in 1/4 window from front illuminated CCDs XIS0 and XIS3 and from the back illuminated XIS1 were accumulated from annular circular extraction regions with 1. 5 radius. CCD areas with more than 4% pileup (Nowak et al. 2011) were excluded. Depending on source brightness, the circular exclusion region had a radius of typically 0. 5. The energy ranges 1.72 keV-1.88 keV and 2.19 keV-2.37 keV were ignored due to calibration issues around the Si-and Au-edges (Nowak et al. 2011). For the back illuminated XIS1, the energy range in between the edges (1.88 keV-2.19 keV) was ignored as well due to significant discrepancies with respect to the front illuminated XIS0 and XIS3. The spectral channels of XIS were added such that a combined signal-to-noise ratio of 40 was achieved. For timing analysis, XIS3-light curves with 2 s resolution were used.
Suzaku's collimated hard X-ray detector (HXD; Takahashi et al. 2007) can be used to detect photons from 10 keV up to 600 keV. It consists of two layers of silicon PIN diodes for energies below 50 keV above phoswich counters sensitive >57 keV (GSO). Spectra from the energy ranges of 12-50 keV (PIN) and 57-100 keV (GSO) were used. The PIN energy channels were binned to achieve a signal-to-noise ratio of 20. For the GSO spectrum, a binning factor of 2 was applied to channels between 60 keV and 80 keV and of 4 for higher energies.  Coe et al. (2007) and in this paper. The parameters e, a sin i, and ω remain unchanged in the analysis presented here. Uncertainties of the new values are at the 90% confidence level (χ 2 red = 1.04 for 1024 dof). Coe et al. (2007) Eccentricity e = 0.68(2) Projected semi major axis a sin i = 530 (60)  Note that the pile up in the 2007 XIS data described above was not taken into account in the previous analysis of these data by Naik et al. (2011). If not excising the high count region of the core of the point spread function, photons which arrive in the same or adjacent pixels are misinterpreted by the detection chain as one photon of larger energy. If the joint XIS/HXD spectra are then modeled, this spectral distortion results in the apparent appearance of soft spectral components claimed by Naik et al. (2011). These components are therefore not present in the real source spectrum.

Orbit Determination Using Pulse Arrival Time Measurements
Between the 2007 outburst and the reported periastron passage of the original orbit determination by Coe et al. (2007) using CGRO-BATSE, GRO J1008−57 orbited its companion 21 times. The accumulated uncertainty of the ephemeris of Coe et al. (2007) during this time is about 8.5 d (0.03 P), or almost as long as a typical outburst of the source. Not correcting for this uncertainty in the analysis of the pulsar's spin could introduce false spin-up/spin-down phases. In order to derive the orbit of GRO J1008−57 the method of pulse arrival timing was used (see also Deeter et al. 1981;Boynton et al. 1986, and references therein). Assuming a neutron star with a slowly varying pulse period, the arrival time t(n) of the nth pulse is given using a Taylor expansion of the pulse ephemeris up to the third order, where P 0 is the pulse period in the rest frame of the neutron star at the reference time t 0 , a x sin i is the projected semimajor axis of the neutron star's orbit andṖ andP are the changes in the spin period. The function F describes the time delay due to an  Table 2 are shown in black with (1) and in purple without a spin-up (2). b) Residuals of the model without taking the times of maximum source flux during outbursts into account. c) Residuals of the best-fit model (1). orbit of eccentricity e, longitude of periastron ω, time of periastron passage τ, and mean anomaly θ. The latter is connected to the orbital period P orb via θ = 2π(t − τ)/P orb (see also Kelley et al. 1980;Nagase 1989). Note that for reasons of computational speed and numerical accuracy, t(n) should be evaluated using a Horner (1819) scheme.
To determine the pulse arrival times, a template pulse profile determined by folding part of an observation with a locally determined pulse period (Leahy et al. 1983) was correlated with the observed light curves. Note that an individual profile was created for each instrument used, so that changes in the profile shape caused by time evolution or by different energy ranges and responses of the detectors were taken into account.
Having found the arrival times, the difference between measured and predicted arrival time according to Eq. (1) was minimized by varying the orbital parameters as well as P 0 ,Ṗ, and P using a standard χ 2 minimization algorithm. Since the template pulse profile depends on the orbit correction and pulse ephemeris, i.e., on the fit parameters, the whole process was iterated until convergence was reached.

The Orbit of GRO J1008−57
To determine the orbit of GRO J1008−57 all data available from the 2005 and 2007 outbursts were used. Template pulse pro-  Figure 4 shows the difference between the measured pulse arrival times and a pulse ephemeris which assumes constant pulse periods with values given above. Calculating the arrival time using the ephemeris of Coe et al. (2007, see Table 2) did not result in a good description of the data (Fig. 4, gray lines). In principle, the remaining difference could be explained by a very large spin-up and spin-down. If the difference is explained in this way, the resulting orbital parameters are, however, inconsistent with the phasing of the outbursts seen with all sky monitors (see Sect. 3.3 below). The only orbital parameters that can result in the difference between the gray line in Fig. 4 and the data are the epoch of periastron passage, τ, and the orbital period, P orb . A change in eccentricity, e, ω, or a sin i shifts arrival times in a way that would move the lines in Fig. 4 up and down or stretch the whole figure. Neither change results in a good description of the data. Holding all parameters except P 0 , τ, and P orb fixed, and performing a χ 2 minimization then results in a good description of all arrival times (Fig. 4, black and purple curves).
Close inspection of the best fit revealed, however, a deviation between the best fit and the data from the early measurements during the 2007 outburst, which were taken during the maximum of the outburst. This discrepancy at high luminosities during the 2007 outburst can be explained by assuming that during this phase the pulsar underwent a spin up, i.e.,Ṗ 0 s s −1 and P 0 s s −2 , between MJD 54426 and MJD 54434. This approach is not unreasonable since high luminosity phases are coincident with a high angular momentum transfer onto the neutron star, which has been detected in other transient X-ray binaries as well, e.g., A0535+26 (Camero-Arranz et al. 2012).
Using the resulting orbital parameters, a first inspection of the orbital phases where the outbursts of GRO J1008−57 are detected in RXTE-ASM and Swift-BAT showed that the times of maximum source flux are consistent with a single orbital phase (see Sect. 3.3). To enhance the orbital parameter precision further, especially the orbital period, the times of maximum source flux were fitted simultaneously with the arrival times data. The residuals are very similar to those found previously (compare residual panels b and c of Fig. 4).
The best fit parameters of the final orbital solution and pulse period ephemeris are presented in Table 2. Note that the simultaneous fit of the 2005 and 2007 arrival times data as well as the outburst times in ASM and BAT allow to break the correlation between the orbital period and the time of periastron passage which is generally found when studying pulse arrival times from one outburst only.

ASM-and BAT-analysis
The orbital period of 249.48(4) d found by the combined pulse arrival times and outburst times analysis in the previous section differs by more than 4σ from the value found by Coe et al. (2007, P orb = 247.8(4) d) and by slightly more than 1σ from the orbital period found by Levine & Corbet (2006, P orb = 248.9(5) d) from the RXTE-ASM light curves. Since outbursts occur during each periastron passage, as a consistency check P orb can also be determined by measuring the time of peak flux from the available RXTE-ASM and Swift-BAT data only. The average distance between outburst peaks is 249.7(4) d, which is in agreement with the results of the arrival times analysis of Sect. 3.2 and Levine & Corbet (2006).
Having confirmed the updated orbital period, it is possible to compare the time of maximum source flux with the periastron passages predicted by the revised orbital solution. With the exception of the bright 2007 outburst (see below), the results displayed in Fig. 6 show that outburst maxima are consistent with a mean orbital phase of φ orb = −0.0323 (17) where the uncertainty includes the uncertainties of the orbital parameters. Thus, GRO J1008−57 reaches maximum luminosity during outburst very close to, but significantly before, periastron.
Following the improved orbital parameters and the detected phase shift between the peak flux and the periastron passage, the date of the highest flux of GRO J1008−57 during an outburst can be predicted: where n is the number of orbits since the outburst in 2007. The uncertainty of T max is about 3 d. It is mainly due to the scattering of the outburst times (Fig. 6). A successful prediction of the 2011 April outburst resulted in the RXTE observations during the rise of this outburst, which are analysed in this work. The ASM analysis shows that the 2007 outburst of GRO J1008−57, the brightest outburst seen by RXTE-ASM, was delayed by ∼ 11 days compared to the standard ephemeris. The time of maximum luminosity from this outburst was therefore excluded from the analysis above. However, as the typical duration of one outburst of GRO J1008−57 is ∼14 days (0.056 P orb ), even this outburst is clearly connected with a periastron passage of the neutron star.

Continuum of GRO J1008−57
The continuum emission of accreting neutron stars is produced in the accretion columns over the magnetic poles of the compact object (Becker & Wolff 2007, and references therein). Here, the kinetic energy of the infalling material is thermalized and a hot spot forms on the surface that is visible in the X-ray band as a black body. A fraction of that radiation is then Compton upscattered in the accretion column above the hot spot, forming a hard power law spectrum with an exponential cutoff at around 50-150 keV.
Although theoretical models now start to emerge which allow to model the continuum emission directly (Becker & Wolff 2007;Ferrigno et al. 2009), these models are not yet selfconsistent enough to describe the spectra of all accreting neutron stars. For this reason, empirical models are typically used (see Müller et al. 2012 andDeCesar et al. 2013 for recent summaries).
Due to the mentioned exponential cutoff of the spectrum, data from instruments such as the CGRO-OSSE and -BATSE mainly show the exponential roll over above 20 keV. As a result, the combined CGRO-OSSE and -BATSE spectrum (20-150 keV) during the discovery outburst of GRO J1008−57 in 1993 was initially described by a bremsstrahlung continuum (Shrader et al. 1999). As discussed by Shrader et al. (1999) for the joint 1993 CGRO/ASCA-GIS (0.7-10 keV) spectrum and also by Coe et al. (2007) for INTEGRAL data from the 2004 June outburst, the soft X-ray spectrum cannot be described by a bremsstrahlung continuum. Instead, a power law spectrum with an exponential cutoff of the form was found to be a good description of the spectrum. Here Γ is the photon index and E fold the folding energy.
Applying this model to the 2007 Suzaku data gives a good description of the spectrum above ∼5 keV (Fig. 7a). At softer energies deviations are present. Note that these differences are not the same as those caused by pile up discussed in Sect. 2.3. The broad band deviations are well described by a black body spectrum modified by interstellar absorption. A narrow iron Kα fluorescence line at 6.4 keV is needed to achieve an acceptable fit with a reduced χ 2 red = 1.14 for 1896 degrees of freedom (dof; see  (Coe et al. 2007), results in a poor description below that energy (gray residuals). c) Adding a black body and a narrow iron line, and taking absorption in the interstellar medium into account, improves the fit to χ 2 red ∼ 1.12 (1845 dof).
line. This additional Kβ line is included in all spectral fits in this paper. The χ 2 can be further decreased to χ 2 red = 1.11 (1888 dof) by ignoring the GSO spectrum. The reason are residuals above 75 keV, which might be caused by the putative cyclotron line at 88 keV (Shrader et al. 1999). The signal quality does not, however, allow to constrain the cyclotron line parameters properly. See Sect. 4.3 for further discussion of that feature. Adding the simultaneous RXTE spectrum (93032-03-03-01) to the Suzaku data, and taking the Galactic ridge emission into account, does not appreciably change the best fit parameters. Table 3 2 lists the final best fit parameters. Applying the same model to the quasisimultaneous Swift-XRT, RXTE-PCA and RXTE-HEXTE data (RXTE obsid 93032-03-02-00) from the same 2007 outburst also gives a good description of the spectrum (see Fig. 8 and Table 3 2 ).
In summary, the broad band spectrum of GRO J1008−57 can be described by a model of the form where Fe 6.4 keV is a Gaussian emission line, BBODY describes a black body, and where the Galactic ridge emission is where TBnew is a revised version of the absorption model of Wilms et al. (2000, see also Hanke et al. 2009), using the abundances of Wilms et al. (2000), and where BREMSS is a bremsstrahlung spectrum. The best-fit hydrogen column densities, N H , are consistent between both fits. They are also in  from the interstellar reddening of E is (B − V) = 1.79(5) mag towards GRO J1008−57 (Riquelme et al. 2012, converted to N H as outlined by Nowak et al. 2012).

Combined parameter evolution
A closer inspection of the best fits from both multi-satellite campaigns in Table 3 shows that most of the continuum parameters appear to be constant within their uncertainties. The only significant parameter change between both models is a change in the power law flux by a factor of ∼2 due to the change of flux of the source over the outburst. Fitting all RXTE observations of the 2007 outburst with the basic model of Eq. (4) gives acceptable χ 2 -values for all observations. With the exception of the fluxes of the two continuum components, other spectral parameters appear to be constant within their error bars, although some scatter is still visible. Due to the lower sensitivity of RXTE compared to Suzaku, however, this scatter could be purely statistical in nature.
To illustrate the range of parameter variability, Fig. 9 shows contour maps between several continuum parameters calculated from RXTE data during maximum luminosity in 2007 (93032-03-02-00) and at the end of the 2007 outburst (93032-03-03-04). The contour map on the left demonstrates that the black body flux F BB changes significantly over the outburst, while its temperature kT BB can be described by the same values. Keeping this value fixed during both observations at a consistent value improves the uncertainty on the remaining parameters significantly. In particular, a contour map of Γ and E fold reveals a strong correlation. Similar to the black body temperature and the break energy, the folding energy seems to be consistent with a constant value during the outburst, while the photon index changes.
Inspection of the confidence maps reveals that the only spectral parameters for which significant changes are seen in individual spectral fits appear to be the power law flux, F PL , the photon index Γ, and the black body flux F BB . To constrain these parameters well and to reveal the evolution of the remaining parameters, all datasets available during the 2007 outburst (see Table 1) including the simultaneous Swift and Suzaku are fit simultaneously with the spectral continuum of Eq. (4) (and the Galactic ridge). In these fits the continuum parameters N H , E fold , and kT , and the flux calibration constant C HEXTE are not allowed to vary between the different observations. The free parameters of the simultaneous fit are the fluxes of the spectral components, F PL and F BB , the photon index, Γ, the iron line equivalent width W, and the PCA background scaling factor b PCA . Since the initial fit showed that the equivalent width of the iron line W does not change significantly with source flux, in a final fit this parameter is also assumed to be the same for all observations.
The resulting model describes all 31 spectra of the 2007 outburst with an χ 2 red = 1.12 for 2623 dof. Adding all RXTEdata of the 2005 and 2011 April outbursts to check for possible hysteresis effects and refitting shows that the simplified model still gives a good description of all data if the equivalent width of the iron line is allowed to change between different outbursts (but remaining constant during each outburst). With N H = 1.6 × 10 22 cm −2 , the soft X-ray absorption is too low to be constrained by PCA-data only. For this reason, changes of N H between the outbursts cannot be detected in the PCA data and the column density is held fixed for all outbursts. The best-fit χ 2 for this fit of 68 spectra is χ 2 red = 1.10 for 3639 dof. Finally, to constrain the Galactic ridge emission further, the summed spectrum of the quiescent state in 1996/1997 is added to the simultaneous fit.
The final combined analysis of all 2005, 2007, and 2011 April data (compare Table 1) results in a remarkable fit with an χ 2 red ≈ 1.10 for 3651 dof. The flux independent parameters, N H , E fold , and kT , the outburst dependent parameters (iron line equivalent widths), the Galactic ridge parameters, and the calibration constant are presented in Table 4.
Due to the combined fits of all available data, the uncertainties of all parameters can be reduced significantly by comparing Table 4 with Table 3. This is not only seen for the flux independent parameters, but also for the continuum parameters allowed to vary for each observation, especially for those at low source luminosity, where only four parameters are required to describe the source spectrum. Apart from an instrumental parameter, the PCA background factor, b PCA , the three physical parameters describing the evolution of the spectral continuum are the photon index Γ, the black body flux, F BB , and the power law flux, F PL . As shown in Fig. 10 the three parameters are strongly correlated: the photon index decreases with log F PL , i.e., the source hardens with luminosity, and the black body flux increases linearly with powerlaw flux. Note that the trend from the rise of the 2011 outburst (purple diamonds) is equal to the trends of the declines of both outbursts in 2005 and 2007. Thus, the spectral shape does not depend on the sign of the time derivative of the flux.
The relationships shown in Fig. 10 can be used to reduce the number of parameters needed to describe the spectrum of GRO J1008−57 to one, F PL , if the flux independent parameters are fixed to the values listed in Table 4. A fit to the data shown in Fig. 10 gives the following empirical relationships: and F BB 10 −9 erg s −1 cm −2 = c + d F PL 10 −9 erg s −1 cm −2 with a = 0.834(10), b = −0.184(8) and c = −0.019(7), d = 0.191(9) (χ 2 red = 1.82, 38 dof and χ 2 red = 1.09, 24 dof, respectively). Using the c parameter of the straight line describing the black body flux one finds that the black body starts to contribute significantly once the power law reached a 15-50 keV flux of 10(4) × 10 −11 erg s −1 cm −2 . This corresponds to a total luminosity of 13(5) × 10 35 erg s −1 assuming a distance of 5.8 kpc (Riquelme et al. 2012).
In both fits, the data points resulting from the Suzaku-spectra were ignored, since they differ significantly from the best-fit Table 4. Source flux independent parameters, the iron line equivalent widths, the parameters of the Galactic ridge emission, and the HEXTE calibration constant as determined from the combined spectral analysis (χ 2 red /dof = 1.10/3651). to the correlations. This difference is due to energy calibration differences between RXTE and Suzaku, which lead to a slight change in the fitted powerlaw photon index Γ of ∆Γ ∼ 0.1. Such differences in photon index are typical between missions. See, e.g., Kirsch et al. 2005, who find differences up to ∆Γ = 0.2 between different missions in their analysis of the Crab pulsar and nebula.

Cyclotron Resonant Scattering Feature
As discussed in Sect. 4.1, a possible absorption-like feature above 75 keV is visible in the Suzaku-GSO and RXTE-HEXTE spectra (Figs. 7 and 8, respectively). Initial attempts to fit this feature were unsuccessful due to its rather low significance and strong correlations between the continuum parameters and the line parameters. The flux dependency of the continuum parameters (Eq. 6 and 7) and the large number of flux independent spectral components ( where E cyc is the centroid energy, W the width, and τ the optical depth of the cyclotron line. Performing a fit to the simultaneous Suzaku and RXTE spectrum measured in 2007 results in a satisfactory description of the cyclotron line. The best-fit parameters together with the historic values determined using CGRO-OSSE data are shown in Table 5. The depth is fixed to the value found by Shrader et al. (1999) since the signal-to-noise ratio in the line  core is close to unity, allowing the modelled depth to take any arbitrary value. Although the remaining parameters are consistent with previous results, the data quality at these high energies does, however, not allow constraining any parameters well. In addition, the χ 2 -contours of the line parameters plotted in Fig. 11 reveal strong correlations. The centroid energy of the absorption feature is consistent with nearly any value above 75 keV. Even a combined fit of all RXTE data with the CYCLABS model added does not improve the fits nor does it constrain the parameters better. Despite the large number of observations presented here, it is therefore difficult to distinguish between a cyclotron resonant scattering feature or a modification of the shape of the high energy cutoff away from a pure exponential cutoff at these energies.

The 2011 December Outburst of GRO J1008−57
In 2011 December, GRO J1008−57 underwent an even stronger outburst than the outburst of 2007, which was analyzed in the previous sections. The peak flux in Swift-BAT of ∼320 mCrab was reached right at the expected date around MJD 55913 predicted by the orbital ephemeris from Sect. 3. Swift ToO observations obtained by the authors as a result of this prediction are listed in Table 1. No further X-ray data of sufficient quality of this outburst are available in the archive and, thus, the high energy part of the spectrum above 10 keV is not available for analysis. The measurements made during this outburst, which did not enter the analysis in the previous section, provide a good test for the overall reliability of the simple continuum model behavior found from the earlier outbursts. Fixing the flux independent continuum parameters and N H -value to the values listed in Table 4, and expressing Γ and F BB as a function of power law flux through Eq. 6 and 7, the two free fit parameters are the power law flux and the combined iron line equivalent width.
An initial combined fit of all four Swift-XRT spectra yields an unacceptable fit of χ 2 red = 3.7. This result is due to residuals below 4 keV, which is strong evidence for a change in the absorption column density. Allowing this parameter to vary results in a reasonable fit of χ 2 red = 1.23 for 934 dof. There are, however, slight deviations from a perfect fit left (Fig. 12b). If a partial coverer model of the form where f is the covering factor, is applied to the model, the quality of the fit is slightly increased (χ 2 red = 1.12 for 932 dof). Note, however, that the spectra can also be described by adding a secondary black body component with a temperature of 0.25 keV. This latter approach results in an even better description of the data (χ 2 red = 1.07 for 922 dof). The parameters of both possible models are shown in Table 6 and the spectra with the secondary black body model applied are shown in Fig. 12. Data from the outbursts of GRO J1008−57 in 2012 show that there is indeed an additional soft excess instead of a partial covering material (see Sect. 6 and 7 below).
Twelve days after the last Swift pointing during the outburst, on MJD 55932, Swift started to observe GRO J1008−57 again. Until MJD 55958 six additional observations were made (see Table 1). Even though the Swift-BAT did not detect the source (Fig. 13), Swift-XRT images clearly reveal an X-ray source at the position of GRO J1008−57, indicating that the neutron star was still accreting. Although the extracted Swift-XRT spectra have a few hundred photons only, by performing a combined fit of all six spectra and using the same model as described above (without need of a partial coverer model or secondary black body), the data can be well described with an χ 2 red = 1.05 for 65 dof. The final parameters are listed in Table 7. The hydrogen column density is slightly larger compared to the value during the main outburst in 2011 December. This can be explained by the position of the neutron star on the orbit, which is behind the companion as seen from Earth. Here, the line of sight crosses a much larger region within the system (see Fig. 2), and thus increased absorption caused by the normal orbital modulation of N H is expected (see Hanke et al. 2010 and Hanke 2011 for a discussion of similar behavior in the HMXBs Cygnus X-1 and LMC X-1).
To conclude, the model of GRO J1008−57 found using the earlier data analyzed in the previous sections works for other   (Table 1), scaled to the PCA via the detector calibration constant c XRT (Table 3). Table 7. Spectral parameters of a fit to the Swift-XRT spectra after the 2011 December outburst (χ 2 red = 1.01, 65 dof).
outbursts as well. GRO J1008−57's behavior is therefore independent of the outburst history.

Predicting the Unpredictable: The 2012 November Type II Outburst of GRO J1008−57
Although GRO J1008−57 seems to be predictable in both, the occurrence of outbursts and the spectral shape, an unexpected and unpredictable behavior occurred right after the outburst in 2012 August, which itself occurred at expected time. Instead of fading into quiescence the outburst decay lasted several weeks, during which several flares were detected (see Fig. 13). An epoch folding (Leahy et al. 1983) of the flaring activity in the light curve reveals significant variability on timescales on the order of 10 d and above (see Fig. 13, inset).
Swift observed GRO J1008−57 during the flaring activity on 2012 September 26/27. Due to the low flux (∼1.14 × 10 −9 erg s −1 cm −2 ), the XRT data were grouped to achieve a a minimum signal-to-noise ratio of 5. The same model as used for the previous outbursts was then applied to the data, which re-  August outburst, which occurred as predicted by the orbital ephemeris (green region), and the unexpected flaring activity after the latter outburst (blue region) followed by a giant type II outburst at an orbital phase of nearly 0.3 (red region). Times of observation by Swift (red) and Suzaku (green) are shown as arrows on top (compare Table 1). The inset shows the epoch folding result of the flaring part of the light curve.
sulted in an good fit (χ 2 red = 0.97, 117 dof). The corresponding spectral parameters are listed in Table 8. Although the best-fit iron line equivalent width is much larger than during previous outburst, the lower confidence limit is still in agreement with the value during the 2011 December outburst.
After the flares, on 2012 November 13, the BAT instrument onboard Swift was triggered by a sudden flux increase of GRO J1008−57 (Krimm et al. 2012). The flux reached 1 Crab within a few days, indicating a further outburst of the source. This flux is about three times higher than the strongest outburst detected previously except of the discovery outburst (Shrader Table 8. Spectral parameters of a fit to the Swift-XRT spectrum during the flaring activity after the 2012 August outburst (χ 2 red = 0.97, 117 dof).  (Table 3).
et al. 1999) and nearly one order of magnitude higher than the mean maximal flux over all outbursts detected by ASM and BAT (Fig. 1). According to the orbital parameters listed in Table 2 the unexpected outburst occured close to orbital phase 0.3 (see also Nakajima et al. 2012). This unusual behavior is typical for type II outbursts, which have been seen in many other Be X-ray binaries such as A0535+26  or EXO 2030+375 (Klochkov et al. 2008).
The spectra of GRO J1008−57 as observed by Swift-XRT, Swift-BAT and Suzaku-XIS are shown in Fig. 14. The model applied to the data is the same as used for the previous outbursts. The Swift spectra are in good agreement with this model after excluding calibration features at ∼1.8 keV and ∼2.2 keV, accounting for the calibration of the Si K edge and the Au M edge, respectively (Godet et al. 2007;Hurkett et al. 2008). The BATspectrum was rebinned to a signal-to-noise ratio of 1.5 in the energy range of 14 to 40 keV. Due to the high flux level of the source data from Suzaku-XIS were extracted avoiding regions with more than 2% pile-up fraction.
An initial fit reveals residuals, which are most prominent in Suzaku-XIS and is of similar shape to the residuals seen in the 2011 December outburst (compare Fig. 12). A partial covering model is not able to fit the residual structure. Using a secondary black body component results in a better description of the data. Its temperature of around 0.39 keV is similar to the value found in the 2011 December data. A complex structure below 3.5 keV remains in the residuals of Suzaku-XIS, which is not present in Swift-XRT data. Since the overall continuum shape is modelled successfully, XIS-data below 3.5 keV were removed from the analysis for the time being.
In addition to the continuum model, an unresolved residual structure consistent with the presence of three narrow Fe emission lines at 6.4, 6.67 and 7 keV is visible in the spectra (see Fig. 15), corresponding to neutral Fe Kα fluorescence and recombination lines from H-like and He-like iron. The centroid energies of the fluorescent lines, which have been fixed relative to neutral iron emission at 6.4 keV, are shifted by ∼ +0.1 keV between the front illuminated XIS1 and the back illuminated XIS0 and XIS3. This is probably due to calibration issues, such as a gain shift. To ensure that the line parameters are not affected by this energy shift, data from XIS1 have been excluded in the final fit. With these modifications, a good description of the data is obtained (χ 2 red = 1.08, 2582 dof). The equivalent widths of these lines together with the fluxes, the secondary black body parameters and absorption column density are listed in Table 9.
During the 2012 November outburst the hydrogen column density is comparable to the preceding flaring activity and the 2007 December outburst, which itself is in agreement with the galactic hydrogen column density (Tables 8 and 3). As the observations during and after the 2011 December outburst show, the column density might still change between outbursts (Tables 6  and 7). The extreme flux increase during the 2012 giant outburst Table 9. Spectral parameters of a fit to the five Swift-XRT, Swift-BAT, and Suzaku-XIS spectra of the giant 2012 November outburst (χ 2 red = 1.08, 2582 dof). requires a large amount of material. This is in contrast to the relatively low hydrogen column density. The ionized iron in the system indicates, however, the presence of significant amounts of almost fully ionized plasma, probably due to photoionization caused by the very luminous neutron star, which would not be picked up in absorption by X-ray spectroscopy. This would effectively cause a decrease of the measured neutral hydrogen column density down to the galactic value. Performing a pulse period analysis of the five Swift-XRT lightcurves during the giant outburst, which have been corrected for binary motion using the parameters listed in Table 2, reveals a significant faster rotation period of the neutron star of P = 93.6483(7) s and a spin-up ofṖ = −0.60(4) × 10 −7 s s −1 .

Conclusions
In the previous sections, all available RXTE observations of GRO J1008−57 including three outbursts in 2005, 2007, and 2011 April as well as an observation campaign during source quiescence in 1996/1997 were analyzed. In addition to these data Swift and Suzaku observations during the 2007 outburst were combined with the RXTE analysis, and Swift observations of the 2011 December and "giant" 2012 November outbursts were compared to the spectral results of the previous outbursts.
By performing a detailed pulse arrival times analysis the orbital period of the binary and the time of periastron passage were improved (see Table 2). In addition, a slight non-linear spin-up of the neutron star in the order of 10 −14 s s −2 was detected during high luminosity of the 2007 outburst. This pulse ephemeris implies a spin-up rate of −2 × 10 −8 s s −1 at maximum luminosity Suzaku-XIS spectra of the type II 2012 November outburst of GRO J1008−57. Spectral channels are rebinned for display purposes. b) Combined residuals after applying the same model to the data as used in the previous analysis including the secondary black body (see text). c) Best-fit combined residuals after excluding calibration features between 1.7 and 2.4 keV and adding three narrow Fe lines at 6.4, 6.67 keV, and 7 keV (compare Fig. 15). around MJD 54428. During the giant 2012 outburst, the measured value of −6 × 10 −8 s s −1 is three times higher compared to the 2007 value. This difference is consistent with theory, where the spin-up rateṖ is connected to the luminosity L via (Ghosh & Lamb 1979) −Ṗ ∝ L α with α = 1 for wind-and α = 6/7 for disk-accretion. The maximum luminosity during the giant 2012 outburst was about three times larger than during the 2007 outburst.  Fig. 16. Spin period evolution of GRO J1008−57 since its discovery. Black: historic data taken from Stollberg et al. (1993), Shrader et al. (1999), and Coe et al. (2007). Red: results of this work. The latest period measurement was taken during the "giant" 2012 outburst. To guide the eye, the dashed line shows a spin-down of 2 × 10 −10 s s −1 .
The evolution of the period over the decades as shown in Fig. 16 reveals, however, that the longterm evolution is dominated by a spin-down in the order of 2 × 10 −10 s s −1 . A similar long term behavior was detected in A0535+26 by Camero-Arranz et al. (2012, see also Bildsten et al. 1997. Such long term spin-down might be evidence for the "propeller" regime, where matter near the neutron star gets expelled, removing angular momentum and causing the neutron star to spin-down (Illarionov & Sunyaev 1975;Frank et al. 1992;Bildsten et al. 1997). In the propeller regime, the magnitude of spin down iṡ where B is the magnetic field strength on the polar caps of the neutron star and where all other symbols have their usual meaning. Assuming a mass M = 1.4 M and a radius R = 10 km of the neutron star andṖ as estimated from Fig. 16, the surface magnetic field strength of GRO J1008−57 is B ∼ 2.6×10 12 G. This is a reasonable value for an accreting neutron star in a HMXB and implies a fundamental CRSF at 30 keV in the X-ray spectrum. Note, however, that since the spin-up phases during outbursts lead to a lower observed long-term spin-down, the derived magnetic field strength without taking this spin-up into account is most likely a lower limit. Using the precise orbital solution allows a detailed study of the connection of outbursts with the orbit. In case of GRO J1008−57 all outbursts detected in RXTE-ASM and Swift-BAT until 2012 May are well connected to the periastron passage, with the peak flux occurring close to periastron (Fig. 6). This behavior is in line with results from a recent study of Be outburst behavior by Okazaki et al. (2012). These authors show that in Be systems type I outbursts occur close to periastron, where the neutron star can accrete from a tidally truncated Be-disk. According to Okazaki et al. (2012), type I outbursts occur regularly only in systems with high eccentricity ( 0.6), where the Be disk is truncated close to the critical Roche Lobe radius. The regularity of type I outbursts is increased in systems with large orbital periods, where the disturbance of the Be-disk is minimized. Both conditions are fulfilled in GRO J1008−57, for which Okazaki & Negueruela (2001) calculated the disk of the donor to be truncated at the 7:1 or 8:1 resonance and concluded that this system should show regular type I outbursts. A further prediction is that type II "giant" outbursts can happen if there is a misalignment between the Be-disk and the orbital plane (Okazaki et al. 2012). These outbursts would then be triggered, e.g., by increased activity of the donor star as seen, for example, in A0535+26 (e.g., Yan et al. 2012). There are indications that such a misalignment is also present in GRO J1008−57, since the type I outburst maximum happens slightly before periastron (Fig. 6). Note that another Be source, 2S1845−024, shows similar behavior (Finger et al. 1999), although here the regular type I outbursts occur slightly after periastron.
Between the ordinary type I outburst in August 2012 and the giant type II outburst several flares were detected in Swift-BAT (see Fig. 13). An epoch folding (Leahy et al. 1983) of that time range revealed oscillations in the order of ∼ 10 d (see inset of Fig. 13). The origin of these oscillations is unknown.
As shown in Sect. 4.1, the spectrum of GRO J1008−57 can be well described by a cut-off power law with an additional black body component, with some spectral parameters being independent of flux and outburst ( Table 8). The physical reason might be found in the unique properties of the neutron star in GRO J1008−57, such as its magnetic field strength, mass and radius. Without a working physical model describing the spectra of accreting neutrons stars it is, however, not possible to investigate this aspect any further. Nevertheless, GRO J1008−57 is an ideal candidate to test future physical models, which are still in development (see Becker & Wolff 2007, and references therein). For a distance of 5.8 kpc and a neutron star with 10 km radius, the maximum observed source flux during the 2007 December outburst corresponds to an emission area of around 1.5% of the neutron star's surface. Even during the giant 2012 November outburst the derived area fraction was only 5%. These values are in good agreement with estimates of the hot spot size at the magnetic poles (e.g., Gnedin & Sunyaev 1973;Ostriker & Davidson 1973).
Similar correlations as the one discovered between the photon index Γ and the black body flux F BB with the power law flux F PL (Fig. 10) have been seen in other Be X-ray binaries as well (see, e.g., Reig & Nespoli 2013, for a recent discussion). These authors analyzed the PCA-data of the 2007 outburst of GRO J1008−57 and found a similar correlation between the powerlaw and the source luminosity. In addition, they discovered a correlation between the folding energy and the photon index. Taking the HEXTE-data into account, Fig. 9 shows, however, that the data is consistent with a constant folding energy.
In case of GRO J1008−57, the spectrum hardens with increasing luminosity. Recent theoretical work by Becker et al. (2012) describes several spectral states depending on the source luminosity, where sources with luminosities below the Eddington limit are expected to show this kind of correlation. Assuming a source distance of 5.8 kpc (Riquelme et al. 2012) and using the maximum flux values given in Table 3, the peak luminosity of GRO J1008−57 is around 3 × 10 37 erg s −1 . This value is close to the critical luminosity reported by Becker et al. (2012), which scales like L crit ∼ 1.5 × 10 37 B 16/15 12 . Thus, GRO J1008−57 is likely to be a subcritical accretor during type I outbursts.
The most remarkable result of the analysis of the giant type II outburst of GRO J1008−57 is that the spectral model found for type I outbursts (Eq. 4 with flux independent parameters from Table 3 and flux correlations from Eqs. 6 and 7) is also able to describe the type II hard X-ray (>10 keV) spectrum. The main contribution in that energy range is the powerlaw component, which emerges from inverse Compton effects in the accretion column far above the neutron star's surface. To describe the soft X-ray spectrum, however, an additional soft component is needed (the second black body). This suggests a change in the accretion mechanisms near the neutron star's surface, where Comptonization processes take place and shocks may form as described in Becker et al. (2012). Following the findings of these authors, the derived luminosity of GRO J1008−57 during the giant outburst is close to 10 38 erg s −1 , which indicates that the neutron star is accreting supercritically. Thus, the spectral change in the soft X-rays might be due supercritical accretion. A similar black body component had to be introduced for the 2011 December outburst (see Table 6). Compared to the additional soft component of the 2012 giant outburst, its bolometric flux contributes much less. That indicates that the source was close to or in the transition from the sub-to the supercritial accretion regime during the 2011 December outburst. These findings are in agreement with and confirm the state changes in neutron star Be X-ray binaries as proposed by Reig & Nespoli (2013), where they found a softening of the spectrum close or above the critical luminosity.
Neutral material in the vicinity of the neutron star is believed to be the origin of observed fluorescence lines, e.g., the 6.4 keV iron Kα line. For N H 10 23 cm −2 , self-absorption is not important and the iron line equivalent width and the hydrogen column density scale roughly linearly (Eikmann 2012). The generally low observed equivalent widths indicate that the fluorescence line originates in the absorbing material. In most observations, N H is in agreement with the interstellar value, although some intrinsic absorption cannot be ruled out. During the 2012 August outburst, the best fit Fe Kα equivalent width is much larger (although with large uncertainties). This larger equivalent width could be due to the addition of a line from reflection from, e.g., an enhanced Be disk. If the X-rays are reflected by material, the equivalent width increases significantly up to a factor of 100.
The data analyzed within this paper can not confirm the claimed cyclotron line at 88 keV by Shrader et al. (1999, see Sect. 4.3), who argued that this feature is rather the second harmonic. Although there are slight deviations from the continuum visible in both, RXTE-HEXTE and Suzaku-GSO data, the signal does not allow to constrain any parameters. During the "giant" 2012 November outburst of GRO J1008−57, the signal quality of the Suzaku-PIN and -GSO data could be used to further investigate this claimed high energy cyclotron line. These data were not used here since no background data were available at the time of writing. However, Yamamoto et al. (2013) were able to fit the spectra with preliminary background corrections. These authors claimed a cyclotron feature between 74 and 78 keV, close to the feature at 88 keV reported by Shrader et al. (1999).
In conclusion, the fact that the spectrum of several outbursts of GRO J1008−57 can be modelled with the single simple model of Eqs. (4), (6), and (7) shows that the shape and behavior of the spectrum is well understood. Due to the flux independent parameters and the flux correlations found during the analysis (Table 4 and Fig. 10) GRO J1008−57's spectrum at any time in the subcritial state can be described given only the power law flux. This is a remarkable result which has not been seen in any other transient X-ray binary before. GRO J1008−57 therefore shows a well predictable behavior in both, outburst dates and spectral shape. Therefore this source is an ideal target to clarify more detailed aspects of Be X-ray transients: what drives the peak flux of outbursts? Are there correlations or changes in iron line flux and absorption column density? How can physical models of accretion explain the spectral shape of GRO J1008−57? Are there other sources where similar predictable behavior exists, and are there differences to GRO J1008−57? Does the spectral model of GRO J1008−57 work outside of the flux range covered here?