© ESO, 2016

1. Introduction

X-ray binary pulsars (XRBPs) were discovered more than forty years ago with the pioneering observations of the bright sources Her X-1 (Giacconi et al. 1971) and Cen X-3 (Tananbaum et al. 1972). The origin of the pulsed emission was soon understood to be related to the accretion of the material lost by the companion onto a rotating neutron star (NS). If this is endowed with a sufficiently strong magnetic field (B ~ 10¹¹⁻¹² G), the inflowing material from the donor star is halted at the magnetospheric boundary of the compact object and is funneled toward its magnetic poles, forming one or more accretion columns (Pringle & Rees 1972; Davidson & Ostriker 1973). As the NS magnetic field is known to decay over time, accretion columns are expected to be more extended in relatively young binary systems, as are the high-mass X-ray binaries (HMXB; see e.g. Walter et al. 2015, for a recent review). Some binary systems with intermediate properties between HMXBs and low-mass X-ray binaries (LMXBs) also show evidence of extended accretion columns (IMXBs; see e.g. the case of Her X-1; Becker & Wolff 2007, hereafter BW07). In all these cases, the gravitational potential energy of the accreting material is first converted into kinetic energy within the accretion column and then released in the form of X-rays as the plasma decelerates and settles onto the stellar surface. X-ray pulsations are generated owing to the disalignment between the NS magnetic and rotational axes. A peculiar feature that is observed in the X-ray spectra of many HMXBs and IMXBs is the cyclotron resonant scattering absorption line (CRSF), which provides a direct measurement of the NS magnetic-field strength. This can be estimated by using the relation E_cyc = 11.6 B₁₂ × (1 + z)^-1 keV, where E_cyc is the centroid energy of the fundamental CRSF, B₁₂ is the NS magnetic field strength in units of 10¹² G, and z is the gravitational redshift in the line-forming region (Wasserman & Shapiro 1983). In some cases, higher order harmonics of the fundamental CRSF are also observed at higher energies (see e.g. Walter et al. 2015, and references therein).

The continuum broad-band X-ray emission of accreting X-ray pulsars is usually described by using phenomenological models, including an absorbed power law extending up to ~100 keV with a roll-over at ~30–50 keV or a broken power law (Orlandini 2006). In some cases it has been shown, however, that such simplified models are not able to satisfactorily describe high-quality observations (Ferrigno et al. 2009, hereafter F09). The development of refined spectral models linking the emission from X-ray pulsars to the physics of accretion has been limited so far by the complexity of the radiative and dynamical equations that describe the behaviour of the accreting material under extreme gravity and magnetic field conditions. A major step forward in this respect was made by BW07 by elaborating an analytical approximation of the theoretical model for accretion onto highly magnetized NSs proposed previously by Davidson & Ostriker (1973).

Assuming the case of constant temperature and magnetic field in the emitting region, together with a simplified profile for the velocity of the radiation dominated flow in the accretion column, BW07 found an analytical solution to the Compton reprocessing of seed photons in this region (which is emitted either at the base of the accretion column or along the accretion stream). Through the definition of a proper Green function, these authors managed to approximate the complex magnetic bremsstrahlung emission of the electrons in the accretion flow (Riffert et al. 1999, hereafter R99) as a combination of a line emission at the cyclotron energy and an ordinary bremsstrahlung emission at lower energies. The angular-dependent cross section of the magnetic Compton scattering was then computed by using weighted averages. The BW07 model was first revised by Farinelli et al. (2012a, hereafter F12), relaxing the previous assumption on the velocity profile of the accreting material in the accretion column. These authors considered the case of a free-falling plasma settling onto the NS surface, and carried out the integration of the Green function along the accretion column. They limited their analysis to the case of low luminosity X-ray pulsars (e.g. the supergiant fast X-ray transients in quiescence; Walter et al. 2015), such that the bulk of the seed photons to be Comptonized can be assumed to originate only from the blackbody emission at the base of the column.

In this work, we further extended the revised approach of F12 by including in the Comptonization process the seed photons produced by the bremsstrahlung and cyclotron emissions along the accretion column. We also allow the value of the dipolar magnetic field to vary along the cylindrical column, thus exploring the more realistic situation in which the energy and the intensity of the cyclotron emission are directly linked to the height of the accretion column and the local density of the plasma. This improvement is motivated by the results of F09, who applied the BW07 model to the emission of the Be/X-ray binary 4U 0115+63 and found that the cyclotron emission is characterized by a magnetic field lower than the value estimated from the centroid energy of the fundamental cyclotron resonant scattering feature. They suggested that this discrepancy might have been caused by a difference in the height of the regions producing the continuum cyclotron emission and the scattering feature along the accretion column.

In order to test the presently improved model, we used a number of observations collected with BeppoSAX and the joined Suzaku + NuStar facilities of the three sources Her X-1, Cen X-3, and 4U 0115+63, which are known to possess prominent cyclotron lines in their X-ray spectra. The X-ray spectrum of 4U 0115+63 is also known to be characterized by the presence of a very strong and broad emission-like feature at 10 keV, i.e. right below the energy of the fundamental cyclotron absorption line (F09). This “10-keV feature” (Coburn et al. 2002) has been observed in a number of X-ray pulsars, and it is commonly modelled in the literature by using a broad Gaussian line (F09; Suchy et al. 2008; Vasco et al. 2013) whose origin is not well understood. In our model, we show that the 10-keV feature can be reasonably well explained as the emission from the collisionally excited first Landau level (cyclotron emission) broadened by the Comptonization along the accretion column.

We present an overview of the improved compmag model in Sect. 2, and summarize in Sect. 3 the properties of the X-ray pulsars Her X-1, Cen X-3, and 4U 0115+63 for which we carry out a detailed spectral analysis to test the model. In Sect. 4, we describe our data set and provide some details on the data reduction procedures. All results are summarized in Sect. 5 and discussed in Sect. 6. The theoretical considerations for the seed photon emission terms and the accretion geometry are discussed in Sects. 7 and 8, respectively. We provide our conclusions in Sect. 9.

2. Updated COMPMAG model

The first version of the compmag model was presented by F12, and the main difference with respect to the spectral model developed by BW07 is related to the presence of the second-order bulk-Comptonization term in the energy diffusion space operator, together with the possibility of considering a more general velocity profile for the inflowing matter along the pulsar accretion column. This column was approximated as cylindrical and permeated by a constant magnetic field. The seed photons for the Comptonization were assumed to be uniquely from the blackbody radiation at the bottom of the column and diffused through its height. The algorithm implemented by F12 to solve the Fokker-Planck approximation of the radiative transfer equation (RTE) in the accretion column is based on a general finite-difference iterative procedure for the solution of partial differential equations with two variables.

The first version of the compmag model was applied to the case of the supergiant fast X-ray transients as these sources host pulsars which usually display a relatively low X-ray luminosity (≃10³³⁻³⁴ erg s^-1; Farinelli et al. 2012b). As shown by BW07, in brighter accreting X-ray pulsars the contribution of the Comptonized blackbody is virtually negligible compared to that of the Comptonized bremsstrahlung and the cyclotron emission, which thus needs to be considered if bright X-ray pulsars are studied. This extension is presented in this work and is the key difference between the previous and the updated version of the compmag model.

For electrons embedded in a strong magnetic field, the motion in the direction perpendicular to field lines has discrete energy states (the Landau levels), while along the field lines the velocity is distributed according to a function f(v). The electron distribution function f(v) can be determined from the solution of the integro-differential Boltzmann equation with a particle diffusion term and a Coulomb scattering kernel (Riffert et al. 1999). The proton-electron Coulomb interaction without photon emission tends to bring the electrons to a thermodynamical equilibrium situation, for which f(v) is a Maxwellian at the proton temperature T. If the electron velocity is such that , where E_cyc is the cyclotron energy, the Coulomb scattering cross section σ(v → v′) presents some specific features (Neugebauer et al. 1996). The first consists of two narrow peaks at v′ = ± v, where the one corresponding to forward scattering (v′ = v) is many orders of magnitude higher than that of backward scattering (v′ = −v); in this case there is no photon emission. The second feature, associated with photon emission, is a broad peak around the energies , where and E_el = 1 / 2m_ev² are the electron energies after and before scattering, respectively. In the latter case the intermediate state of the Coulomb scattering brings the electron from the fundamental to higher Landau levels (n ≥ 1) with a decay time which is much shorter than the particle relaxation time. The outcome is a significant loss of the initial electron energy, which causes a depopulation of the tail of the Maxwellian distribution that can be described by using a simple electron-proton Coulomb interaction. The energy threshold of this process implies that the higher the electron temperature kT_e is compared to the magnetic field, the more efficient the depopulation effect is.

Numerical computations performed by R99 showed that the emissivity has a narrow energy peak at the cyclotron energy for photons emitted in the direction perpendicular to the magnetic field, while the peak gets broader for photons progressively emitted in the magnetic field direction. When averaging over all the possible directions, the net result can be qualitatively described as a broad cyclotron emission feature overlapped to a smooth continuum resembling a non-magnetized bremsstrahlung emissivity spectrum. However, a fully self-consistent treatment is particularly difficult and barely feasible in a run-time code (see Neugebauer et al. 1996, for a solution). A suitable approximation can be found by splitting the emissivity into two different seed photon sources (see also BW07), and this is the approach we follow in our treatment. Instead of modelling the cyclotron line emission with a delta function δ(E − E_cyc), we assume a Gaussian profile with a centroid energy of E_cyc and a variable width σ_cyc. The latter quantity parametrizes all the different effects leading to the intrinsic cyclotron line broadening, such as thermal broadening, spatial variations of the magnetic field not related to the dipolar structure, and the release of photons at all emission angles (this is of particular importance in the case of the phase-averaged spectral analyses such as those we carry out in this paper; see Sect. 4).

The RTE in the observer frame for subrelativistic bulk motion and arbitrary geometry was first derived by Blandford & Payne (1981, hereafter BP81). This particular form of the RTE has since been widely used in the literature (e.g. Lyubarskii & Syunyaev 1982; Titarchuk et al. 1997; Farinelli et al. 2008). On the other hand, Psaltis & Lamb (1997, hereafter PL97) outlined how the result derived by BP81 was based on the wrong assumption that the ratio between the zeroth and the second moment of the specific intensity is equal to 1/3, which is strictly true only in the fluid frame where the bulk velocity is zero. The error introduced in the RTE when deriving it in the observer reference frame is at the first order comparable to the adimensional fluid velocity β (see derivation in Appendix B). Under particular assumptions for the matter velocity profile, the RTE can be solved analytically by using the variable separation method for which the solutions are determined by the space boundary conditions for a bounded medium (Sunyaev & Titarchuk 1980; Lyubarskii & Syunyaev 1982; Titarchuk et al. 1997).

It is worth noting that analytical solutions of the RTE with bulk motion, even in the case of ad hoc velocity profiles, can be obtained if the energy diffusion operator depends only on the electron temperature, assumed to be constant through the medium. If second-order bulk Comptonization effects are taken into account, a space-dependent term containing a β² factor appears in the energy operator. In this case, the variable-separation method can no longer be adopted and solutions must be obtained through numerical methods. The contribution of the space dependent term to the solutions of the RTE has been investigated by e.g. Psaltis & Lamb (1997); comparison between analytical and numerical solutions can be found in Titarchuk et al. (1997). As this term provides in all cases more reliable solutions to the RTE, although with an increasing effect depending on the magnitude of β, we included it in our model.

In the presence of a strong magnetic field (B ≳ 10¹² G), two polarization radiation modes arise as a consequence of the vacuum polarization (e.g. Ventura 1979). These two modes are commonly labelled as ordinary (O) and extraordinary (E). The approximate electron-scattering cross section for O-photons is (1)where (2)The approximate cross section of the E-photons is instead (3)where Φ_l is a normalized function of the line profile, and (4)In Eqs. (1) and (3), ϕ is the angle between the photon propagation direction and the magnetic field (Arons et al. 1987). Both O-photons and E-photons interact with the electrons of the plasma with an efficiency dictated by their relative cross sections, but mode conversion also occurs at a given scattering rate (Arons et al. 1987). A system of two equations for both polarization modes containing the transition-rate term should thus be solved. As outlined by BW07, the angular dependence of the cross sections must be replaced by some angle-average values in order to treat the problem in analytical and numerical codes. Additionally, the explicit dependence on energy in Eq. (2) is neglected and instead we use energy-averaged cross sections. On the one hand, E-photons mostly escape from the bounded configuration, except for energies close to the cyclotron value and its harmonics where the scattering cross section becomes highly resonant (see Eq. (3)) and photons are substantially trapped in the line. On the other hand, O-photons are more efficiently scattered, particularly at orthogonal propagation angles with respect to the magnetic field (Eq. (1)). Building on the general approach of BW07, it has been assumed as an underlying assumption in the compmag model that the continuum Comptonization spectrum has to be computed only for O-photons, while the contribution of E-photons is phenomenologically taken into account by using Gaussian absorption features in the X-ray spectra.

The RTE for the zeroth-moment photon occupation number n(x) in the cylindrically symmetric case can be written as

(5)

where x = E/kT_e, Θ = kT_e/m_ec², f_b = 1 + m_ev²/ 3kT_e, β = v/c, σ_∥ = 10^-3σ_T, , and σ_T is the Thomson cross section. Equation (5) can be easily derived from Eq. (18) in BP81 by adopting a few modifications. The first is the addition of the second-order bulk Comptonization term f_b in the energy distribution operator, as previously discussed (see Eq. (33) in PL97 for its derivation). The second modification is the inclusion of the free-free absorption as expressed in the observer reference frame with an accuracy of the order of v/c – higher terms of the series expansion are of the order of v²/c² and can thus be neglected (see Eq. (B4) in PL97). The free-free absorption was not considered in the earlier treatment of F12 or in the analytical treatment by BW07. The absorption coefficient can be written as (6)where (7)is the Gaunt factor and K₀(x/ 2) is the modified Bessel function of the second kind (Titarchuk 1988).

It is important to note that the expression for α_ff used in Eq. (5) is for unmagnetized plasma. In the presence of a strong magnetic field, the behaviour is different and depends on the photon polarization mode: for O-photons α_ff is a smooth function of energy which can be closely approximated by that of unmagnetized case. For E-photons, the Gaunt factor exhibits strong resonance at the cyclotron fundamental energy and its harmonics E_n = nE_cyc (Pavlov & Panov 1976). Since Eq. (5) deals with O-photons, the expression for the absorption coefficient in Eq. (6) is a good approximation, while the free-free absorption resonances of E-photons are not treated in our continuum model. To take into account their contributions in the source spectra, we use phenomenological Gaussian absorption line components. The E-photons mainly interact with the electrons via the cyclotron resonance, but this is effectively an elastic scattering process that preserves the angle-averaged photon distribution (Nagel 1980). Hence, no term describing this process should be included in the RTE.

The function M in the last term of Eq. (5) is the first moment occupation number, defined as (see BP81) (8)The addition of the photon escape term is linked to the solution of the RTE along the vertical direction. Actually, the system also has a finite size along the radial direction, and the rate of photon escape through the walls of the accretion column must be taken into account. In principle, it is necessary to find solutions of the RTE considering both vertical and radial space coordinates (Davidson 1973). However, together with the addition of the energy operator, the problem would become cumbersome and more suitable to be treated by Monte Carlo codes (Odaka et al. 2013, 2014). We parametrize this effect using the characteristic photon time escape defined as t_esc = r₀τ_⊥/c, where τ_⊥ is the optical depth perpendicular to the magnetic field defined assuming a cross section σ_⊥ = σ_T.

As previously mentioned, we consider the bremsstrahlung and cyclotron emissions as separate sources of seed photons in Eq. (5) (see Appendix A). It is worth noting that the emission coefficient expressed in the observer frame for subrelativistic flows comprises a zeroth term plus correction factors of the order of v²/c² (see Eq. (B4) of PL97). As the accuracy of the Eq. (5) is ~v/c, with the exception of the Comptonization operator, we only retained the zeroth-order term for the emission coefficient.

Using the definition dτ = n_eσ_∥dz for the optical depth along the z-axis, performing a logarithmic sampling of the adimensional energy through the change of variable x = e^q, and defining n = J/x³, after some lengthy calculations we obtain the adimensional equation for the zero moment intensity, (9)where , ξ = 15.5r₀/ṁ, and , with H = (σ/σ_∥)Θ.

The parameter r₀ contained in the ξ term, is the adimensional column radius expressed in units of NS Schwarzschild radius through the relation km, while ṁ = Ṁ/Ṁ_Edd is the adimensional accretion rate in Eddington units.

The explicit form of the cyclotron seed photon term as a function of the x variable is (10)with x_cyc = E_cyc/kT_e and B₁₂ = B/ 10¹². We also used the normalized Gaussian function (11)in place of the δ-function corresponding to the monochromatic line emission case. The bremsstrahlung source term is (12)with G(x) given in Eq. (7). The derivation of the source terms in Eqs. (10) and (12) is reported in Appendix A. The functions Φ and Ψ in Eq. (9) are related to the zeroth and first order term of the free-free absorption coefficient and are defined as (13)and (14)where c₁ = 2.1 × 10^-13 and c₂ = 4c₁.

Equation (9) is solved by using the numerical procedure described in F12. We considered as natural boundary conditions that (i) the intensity vanishes at the extremes of the energy domain and at the top of the column; and (ii) the approximation reported in Eq. (31) of F12 holds at the base of the accretion column. We assumed that the NS magnetic field within the accretion column is dipolar, and thus , where is the magnetic field strength at the base of the accretion column (i.e. the NS surface at z₀). In the updated version of the compmag model being developed here, the free parameters regulating the cyclotron emissions to be determined through the fit to the X-ray data are thus and σ_cyc. The additional improvement with respect to the former compmag version that we introduce here is the possibility of treating the height of the accretion column z_max as a free parameter. We also stress that in the new version of the model the mass accretion rate is used as a free parameter in Xspec in place of the vertical optical depth adopted previously (see also Table 2). The relation between these two quantities depends on the form of the velocity profile of the accreting material. Starting from the continuity equation (15)and defining the first adimensional velocity profile as (16)the optical depth measured at the top of the column is (17)where , while β₀ is the matter velocity in units of c at the NS surface, and z₀ is the adimensional NS radius. Given the definition of r₀, we obtain km. Here and throughout the paper we assume m = 1.4 and R_ns = 10 km, implying z₀ = 2.42.

If the inflowing matter is assumed to decelerate towards the NS surface according to the law (see also Eq. (27) of BW07) (18)then, the relations between τ₀ and ṁ is given by (19)It is worth pointing out that the adimensional velocity of Eq. (18) is limited to be less than unity at the top of the column. In the above equation, τ₀ is the approximated vertical optical depth calculated with the reduced Thomson cross section for photons propagating along the lines of the magnetic field: (20)The above reported equations are implicitly derived under the assumption that the accreting column has a pure cylindrical shape with constant radius over height. Actually, for a dipolar magnetic field the column has a cone-like shape of half-angle θ_c given by (Frank et al. 2002) (21)where R_∗ is the NS radius, R_m is the radius at which matter starts to be channelled by the magnetic field lines, and ψ is the inclination angle between the magnetic axis and the NS equatorial rotating plane.

Finally, the last important improvement we introduced in the new version of the compmag model is related to the geometry of the X-ray emission. We now differentiate between the pencil beam case, when the X-ray radiation is emitted upwards with respect to the NS surface, and the fan beam geometry, when the radiation is released from the lateral boundaries of the accretion column. Assuming a cylindrical accretion column, the adimensional photon flux in the pencil beam case is given by (22)where β^t is the adimensional accreting matter velocity at the top of the column.

In order to pass to physical units, it is necessary to consider the relation between the intensity and the occupation number (23)which can be easily derived once hν = xkT_e has been derived (see Eq. (A.3)). The compmag normalization N_comp is defined in this case as the ratio between the energy flux multiplied by the top-column surface , and divided by 4πD², where D is the source distance. Using Eqs. (22) and (23), we obtain (24)where r₀ and d₁₀ are the adimensional column radius and source distance in units of 10 kpc, respectively.

In the fan-beam case, the high optical depth allows the photons to diffuse and escape through the lateral boundaries of the accretion columns. As compmag is a 1D model and the variations of the physical quantities can be computed only along the vertical (z) and not the radial direction, we adopted some simplifications. More specifically, considering the zeroth-moment specific intensity J(x,τ) = x³n(x,τ) at the walls of the column and assuming a uniform brightness B(x,τ) = πJ(x,τ), the vertically integrated adimensional photon flux is given by (25)where dz/ dτ can be easily derived by differentiating either Eq. (17) or Eq. (19), and P_cyl = 2πR₀ is the column perimeter.

From Eqs. (23) and (25) we then obtain (26)We note that we expressed the fluxes in Eqs. (24) and (26) in keV because these are the standard units when writing Xspec models. The conversion to cgs units is performed by the fitting package at the end of the runs.

The compmag normalization factor N_comp, as derived in Eqs. (24) and (26), deserves some consideration. Once the source distance is fixed, the other parameters that determine the observed flux are kT_e and r₀, to be inferred from the fit to the X-ray data. This means that, in principle, N_comp should be kept fixed during the fit to a value of the order of ~, where d₁₀ is the distance of the source in units of 10 kpc. However, in the computation presented above (and also reported in BW07) we considered only the approximated case of a purely cylindrical accretion column and a flat space-time. In a more realistic situation, it is likely that hydrodynamical and General Relativity (GR) effects will affect the spectral energy distribution of the X-ray radiation. In particular, the gravitational-lensing effect – due to the curved space-time close to the NS surface – is expected to enhance the emerging X-ray flux at some angles more than it does in the Newtonian case. For this reason, the normalization N_comp should be proportional to f×, where f> 1 is a numerical factor that parametrizes all the effects that cannot be taken into account (yet) in the present version of the compmag model.

2.1. New convergence procedure of the algorithm

In the first version of the compmag model proposed by F12, the authors developed an empirical criterion to halt the iteration during the resolution of the RTE in the pulsar accretion column. The code first estimated the approximated energy index α from the fit to the source spectrum in the energy range E_min<E<E_max, where E_min> 3kT_bb and E_max ≲ 10 keV. These boundaries were chosen in order to limit the fit to the region where the energy spectral distribution is better approximated with a power law, and the iteration stopped when 1 − α_m/α_{m + 1}<ε. Additionally, the vertical optical depth τ was always divided into N = 10 bins independently of its value in order to achieve a reasonably fast computing time.

This method was shown to be well suited to perform fits to the relatively low-quality data that were available to F12. The data that are used in the present paper to study some of the brightest X-ray pulsars in our Galaxy (see Sect. 4) require, instead, a significant improvement in the convergence procedure for the resolution of the radiative transfer problem in the accretion column. We verified that the presence of cyclotron emission features in the spectra of these sources makes the estimation of the spectral slope through a two-point linear interpolation largely uncertain, resulting in the non-convergence of the algorithm or an unsustainable high number of numerical iterations. We thus changed the initial differential approach of F12 to the convergence problem into an integral one. Defining and as the computed values of the source energy spectral distribution at the mth and m + 1th iterations, with (i,j) identifying the steps in the pre-defined grid of energy and optical depth, the new stopping criterion for the convergence of the algorithm can be expressed as (27)The maximum allowed value of ε is fixed by verifying at each iteration step that the χ² obtained from the fit to the data is monotonically decreasing. We checked a posteriori that ε = 10^-4 achieves the required accuracy in the estimate of the model parameters during all fits performed, maintaining at the same time a reasonably short computational time.

We also checked that for the new data set used in this paper different values of the number N_τ of bins selected for τ resulted in remarkably different values of the model parameters obtained from the fits. After a number of trials, we found that fixing N_τ = 50 provides the best solution, as higher values do not significantly alter the results of the fits and the computational time of each fit is still reasonably short (a few hours on a quad-core Linux machine with a 2.4 GHz processor).

3. Sources

3.1. Cen X-3

The high-mass X-ray binary Cen X-3 was discovered by Uhuru (Giacconi et al. 1971) and is known to host a ~4.8 s spinning NS with an estimated mass of 1.34 ± 0.15 M_⊙ (van der Meer et al. 2007), and an O 6-8 III supergiant star with a mass of 20.5 ± 0.7M_⊙. The orbital period of the system is ~2.1 days, and the NS is on a nearly circular orbit eclipsed by the massive companion for about 20% of the time (Hutchings et al. 1979; Ash et al. 1999). The distance to the source is poorly constrained. A lower limit of 6.2 kpc was suggested by Krzeminski (1974), while Thompson & Rothschild (2009) have more recently derived an estimate of 5.7 ± 1.5 kpc. For the purpose of our work and to ease comparisons with published works, we adopt hereafter a distance of 8 kpc. The system is known to undergo episodes of both disk and wind accretion, as shown by the positive and negative fluctuating derivative of the NS spin period (see e.g. Bildsten et al. 1997).

Prominent fluorescence iron lines for different ionization states are present at 6.4, 6.7, and 6.97 keV in the X-ray spectrum of the source and have been reported at all orbital phases with variable intensities. Observations performed with the gratings on board Chandra also show the presence of Doppler broadened Si xiii, Si xiv, and Fe xxv complexes caused by recombination in a photoionized plasma (Wojdowski et al. 2003; Iaria et al. 2005). It has been argued that the fluorescent lines originate relatively close to the NS (probably in the outer regions of the accretion disk), while the other lines are produced at greater distances from the NS in the wind of the massive companion (Iaria et al. 2005; Naik et al. 2011). The stellar wind velocity measured from these lines is significantly smaller than that expected for an O-type star, providing convincing indications for a strong perturbation of the wind by the NS X-ray radiation.

The source X-ray spectrum is generally well described by a power-law model with exponential cut-off and a Gaussian-like emission feature centred at ~13 keV (Suchy et al. 2008). A cyclotron absorption feature was detected at ~30 keV (Burderi et al. 2000). A partial covering component has also been introduced in some cases to model the variable column density of the source and interpret the hour-long dips in the source lightcurve (see e.g. Naik et al. 2011). The X-ray variability of Cen X-3 ranges from timescales comparable to the orbital period, down to fractions of the NS spin period. The spectral energy distribution also shows a remarkable dependence on the spin and orbital phase (see e.g. Suchy et al. 2008; Raichur & Paul 2008; Naik et al. 2011, and references therein).

3.2. 4U 0115+63

Hard X-ray radiation from 4U 0115+63 was first discovered in 1969 (Johns et al. 1978) during one of the giant-type II outbursts displayed by the source. On that occasion, the peak luminosity of the event reached a few 10³⁷ erg/s, two orders of magnitude above the usual quiescent emission level of the source. A periodicity of 3 yr for the outbursts from 4U 0115+63 was later suggested by Whitlock et al. (1989). The source orbital period was measured by Rappaport et al. (1978) using SAS data at P_orb = 24.3 d. These authors also estimated the eccentricity of the orbit as e = 0.34, and its semimajor axis a_Xsini = 140.1 lt−s. X-ray pulsations with a period of P_S = 3.6 s were discovered by Cominsky et al. (1978).

The source has been widely studied in the optical and IR bands, leading to the identification of its optical counterpart (V 635 Cas; Johns et al. 1978) and the determination of its distance (7–8 kpc; Negueruela & Okazaki 2001). The X-ray spectrum of the source above ~10 keV is usually well described by power-law model with an exponential cut-off (e.g. Coburn et al. 2002). The centroid energy of the fundamental cyclotron line is at ~11 keV, and up to six harmonics have been observed and reported in the literature (Wheaton et al. 1979; White et al. 1983; Heindl et al. 1999; Santangelo et al. 1999; F09). This is a peculiar characteristic of 4U 0115+63 as in all other sources displaying cyclotron lines only the second harmonic is usually detected (the only exception being V 0332+53, which also showed evidence of the third cyclotron line harmonic; Coburn et al. 2002; Orlandini 2004; Tsygankov et al. 2006; Caballero et al. 2007). The variations of the centroid energy of the fundamental cyclotron line with the source luminosity were discussed by Mihara et al. (2004), Nakajima et al. (2006) and Tsygankov et al. (2007), but these results were later criticized by Müller et al. (2013) who showed that slightly different fits to the continuum emission could produce artificial (but significant) variations of the centroid energy. On the other hand, F09 presented the results of the fit to the broad-band X-ray spectrum of 4U 0115+63 with the BW07 model, aiming at constraining the properties of the X-ray emitting region and the high-energy radiation mechanisms. These authors showed that the cyclotron emission has a dominant role in cooling the material that is flowing through the accretion column in this source, and inferred a magnetic field strength from the cyclotron continuum that is approximately a factor of 2 lower than the one estimated from the centroid energy of the fundamental cyclotron line. This suggests that the heights of the continuum and cyclotron line formation along the pulsar accretion column are substantially different.

3.3. Her X-1

Her X-1 hosts a NS with a spin period of 1.24 s and a mass of 1.5 ± 0.3 M_⊙ orbiting around an A/F donor star, HZ Her. The orbital period of the system is 1.7 d and the high inclination angle (~85 deg) with respect to the observer line of sight gives rise to extended X-ray eclipses, lasting about 20% of the binary orbit (Scott et al. 2000). The super-orbital modulation of the source X-ray emission with a period of 34 d is usually ascribed to the presence of a warped accretion disk (Staubert et al. 2009, and references therein) or to the free precession of the NS-disk system (Staubert et al. 2013). During the super-orbital period, Her X-1 displays a “main-on” state lasting 10–11 days where the source achieves the highest X-ray luminosity and a short slightly dimmer state which typically last 5 to 7 days. At other super-orbital phases, the source is still visible in the X-ray domain, although at a much lower luminosity level.

Her X-1 is characterized by low Galactic absorption and is situated at a distance of 6.6 ± 0.4 kpc (Reynolds et al. 1997). The source X-ray spectrum shows emission lines produced from the photoionized accretion disk and on the companion surface, as shown by Chandra and XMM-Newton grating observations (see e.g. Zane et al. 2004; Jimenez-Garate et al. 2005; Ji et al. 2009, and references therein). The fundamental CRSF was initially detected with a centroid energy of 40 keV (Truemper et al. 1977), but it has been shown that this energy has decreased in the past 20 yr by 4.2 keV (Staubert et al. 2014; Klochkov et al. 2015). The centroid energy is also known to positively correlate with the overall source X-ray flux over a wide range of different timescales (from seconds to days; Staubert et al. 2007; Klochkov et al. 2011). This behaviour has been interpreted in terms of the so-called subcritical accretion regime, during which the height of the cyclotron scattering forming region decreases at higher mass accretion rates (Becker et al. 2012). This conclusion was criticized by Mushtukov et al. (2015). At the luminosity of Her X-1 where the cyclotron line is typically monitored, these authors showed that the source should be in the super-critical accretion regime, and thus an opposite trend of the cyclotron centroid energy with the source luminosity should be observed. An alternative possible solution to this problem is that the magnetic field of Her X-1 is strongly non-dipolar, resulting in multiple accretion sites on the NS surface with non-cylindrical sections (Shakura et al. 1991). It is worth noting that the combination of a pencil and a fan-beamed emission from a misaligned dipole would also be able to reproduce the pulse profiles of Her X-1 with a reasonable accuracy (Leahy 2004). This conclusion was also independently confirmed by Blum & Kraus (2000), who used a pulse profile deconvolution method to infer the pattern of the X-ray emission close to the NS surface and the geometrical properties of the corresponding emitting region.

More detailed studies of the source pulse profile were presented by Vasco et al. (2013) and Enoto et al. (2008). The first authors analysed RXTE observations of Her X-1 and found a puzzling double peak in the variation of cyclotron line energy and photon index before and at the same time as the source pulse maximum. Enoto et al. (2008) analysed a 60 ks Suzaku observation of the source and reported on the possible detection of the second cyclotron harmonic close to the maximum of the pulse profile. A extensive analysis of three combined Suzaku/NuSTAR observations has been carried out by Fürst et al. (2013, hereafter F13), who used a phenomenological model to describe the source emission during the main-on state.

4. Observations and data analysis

In order to properly test the improved version of the compmag model presented in this paper, we used a high-quality data set that provides a broad and continuous coverage in the energy range 0.5–100 keV. Our sample exploits archive BeppoSAX (see Boella et al. 1997a, for a description of the mission) data of Cen X-3 and 4U 0115+63, together with a combined NuSTAR (Harrison et al. 2013) and Suzaku (Mitsuda et al. 2007) observation of Her X-1 during one of the source main-on states (see Sect. 3). We summarize all the characteristics of our data set in Table 1, together with the references to previously published papers that made use of these data. We note that the two Cen X-3 observations listed in Table 1 were obtained with the source at different luminosities and we removed the time intervals corresponding to the source eclipse before carrying out any further analysis.

For the BeppoSAX observations, we analysed data from the high-energy narrow field instruments LECS (0.7–4.0 keV, Parmar et al. 1997), MECS (1.5–10.5 keV, Boella et al. 1997b), HPGSPC (7–44 keV, Manzo et al. 1997), and PDS (15–100 keV, Frontera et al. 1997) after performing the full processing through the standard pipeline (saxdas v.2.1, running only on heasoft v. 5.2). Background subtraction was performed by using the Earth occultation for the HPGSPC, the off-source pointing for the PDS, and the standard calibration files for the MECS and the LECS. The LECS and MECS source extraction radii were set to 8′ to optimize the signal-to-noise ratio.

We adopted standard reduction procedures for the Suzaku and NuSTAR data, devoting particular attention to evaluating the effect of pile-up during the Suzaku/XIS observation. In the latter case only the outer portion of the instrument chip was active, and so we made use of data exclusively collected through the optimally calibrated XIS3 (see Fürst et al. 2013, for further details). Considering the scope of this paper, we limited our analysis to the 3–8.5 keV energy range for the Suzaku/XIS, 16–80 keV for the Suzaku/HXD, and 3.5–65 keV energy range for NuSTAR. Any spectral contribution coming from the accretion disk and the surrounding corona was thus excluded from further analysis. We only considered each source spectra averaged over the entire spin period and did not perform a phase-resolved spectral analysis.

All spectral analyses were performed with the Xspec version 12.8 (Arnaud 1996), using the compmag model described in Sect. 2. Inter-calibration constants were included in the fits and found to be compatible with typically expected values¹. The fundamental and all detected harmonics of the CRSF features in each of the considered sources were accounted for in the fits by using the gabs model in Xspec: (28)Here τ_n, E_cn, and σ_n denote the optical depth, energy, and standard deviation of the absorption for the nth harmonic, respectively. The line width is fixed to previously measured values whenever it cannot be robustly determined by the fit. For Cen X-3 and Her X-1, we linked the intensity of the magnetic field in the compmag model used to reproduce the continuum to the centroid energy of the fundamental CRSF; for 4U 0115+63 it was linked to the first harmonic as this is known to be less affected by biases in the continuum modelling of the source (Müller et al. 2013).

5. Results

The updated version of the compmag model that we introduce in the present paper compute the spectra of X-ray pulsars emerging in both a pencil-beam or a fan-beam geometry. The former case can be obtained by using Eq. (16) for the velocity profile of the infalling material in the pulsar accretion column and Eq. (24) to describe the emerging X-ray flux. The parameter β₀ defines in this geometry the terminal velocity of the infalling material at the NS surface. The corresponding equations for the fan-beam geometry are Eqs. (18) and (26).

A preliminary fit to the available spectra considered in this work revealed that reasonably good values of the χ² could be obtained with all the compmag parameters left free to vary, but the uncertainties on these parameters were not very constraining. We verified that this is largely due to the computation in the model of the matter velocity profile inside the accretion column, as the latter determines the optical depth of the column in the vertical direction and, in turn, the total source luminosity. As all sources in this paper are characterized by a luminosity that is far above the value at which the fan-beam emission geometry is expected to set in, we first fixed the normalization of the compmag model (N_comp) and concentrate on the fan-beam case with efficient matter deceleration along the accretion column due to the strong radiation pressure (see Eqs. (18) and (26)). This eases the comparison with both the previous theoretical treatment of the model reported in BW07 and its first application discussed by F09. For completeness, however, we also provide in the next sessions a discussion about the results of the fits to all data with the compmag model and the pencil-beam geometry.

In all cases, we modelled the effect of the interstellar photoelectric absorption on the X-ray emission from each of the considered sources with the phabs spectral component in Xspec, assuming solar abundances.

5.1. Cen X-3

We first fit the low luminosity (LL) and high luminosity (HL) spectra of Cen X-3 obtained from the two available BeppoSAX observations (see Table 1) with the compmag model in the fan-beam geometry case, plus a cyclotron Gaussian absorption line centred at ~30 keV. The results of the fits were not yet acceptable (χ²/d.o.f. = 4049/318 and χ²/d.o.f. = 675/241 for the LL and HL case, respectively) owing to evident residuals around 6.7 keV. The addition of a Gaussian emission line at this energy significantly improved the results of both the LL (χ²/d.o.f. = 509/315) and HL (χ²/d.o.f. = 167/238) observations.

In the case of the LL observation, some residuals were still visible around 5 keV and below 2 keV. The former are most likely due to a known feature of the MECS instrument and were able to be modelled by adding a Gaussian component with a fixed width of σ = 0.1 keV. This improvement led to χ²/d.o.f. = 447/313. The residuals below ~2 keV observed in the LECS data (not available during the HL observation) were interpreted in terms of a second region of absorption partially covering the source. This possibility has been already investigated by Naik et al. (2011) during the phase-resolved spectroscopy of the source performed on some Suzaku data. The authors found a strongly variable local absorption component with N_H ranging from ~2 × 10²² to ~1.3 × 10²³ cm^-2, and a covering fraction in the range 0.3–1. Following the above mentioned work, we added a multiplicative partial covering absorption component (pcfabs) to our Xspec model and obtained a significantly improved result (χ²/d.o.f. = 369/311) with N_H ~ 5 × 10²² cm^-2 and a covering fraction of ~0.4. In the HL spectrum, this component was not necessary to obtain an acceptable fit (see Fig. 1), possibly owing to a higher ionization of the absorbing medium and/or by a weaker constraint of the absorption in the data set due to the lack of LECS data.

Fig. 1 Unfolded spectra of the two BeppoSAX observations of Cen X-3. Also show are the best-fit compmag model and the residuals from the fit in units of σ (bottom panel). All best-fit parameters determined from the fit are given in Table 2.

Although the model discussed above provided a reasonably good fit to both the available observations of Cen X-3, we obtained a similar accretion rate in both observations from the compmag best-fit parameters. This result does not seem to be consistent with the eight times higher X-ray luminosity of the second observation (see Table 2). We verified that this computational issue is due to the degeneracy in the fits between the compmag normalization N_comp and the mass accretion rate ṁ, as they both regulate the source luminosity (see Sect. 2). To avoid this degeneracy, we performed a new fit to the LL data with the compmag normalization fixed to the value obtained from the fit to the HL data (N_comp ~ 60). The fit with this model provided χ²/d.o.f. = 449/312, a slightly worse value than that obtained with a normalization left free to vary, but with the accretion rate parameters ṁ now tracing reasonably well the source luminosity variation between the two observations (see Table 2). No significant changes of the N_H and the covering fraction parameters in the pcfabs components were altered by fixing the compmag normalization. It is noteworthy that the determined value of the N_H is comparable to that measured by Naik et al. (2011) during their phase-averaged spectral analysis (see their Table 1). The most significant change was recorded for the mass accretion rate estimated from the fit, that is now in rough agreement with that expected from the ratio of the X-ray luminosities measured from the LL and HL observations (see Table 2). For completeness, we also tried to fix the normalization of the compmag model in the fit of the HL data to the value obtained from the original fit to the LL data (with N_comp left free to vary). However, in this case, we were not able to obtain an acceptable fit to the HL data. From the values of the different parameters obtained through our best-fit models (Fig. 1), we measured a fraction of about 80% and 30% for the Comptonized cyclotron emission of the LL and the HL observation, respectively (see Fig. 4 and Table 2).

To perform a fit to the HL and LL data with the pencil-beam geometry (assuming a free-fall velocity profile for the material within the accretion column), we have to introduce the flow terminal velocity β₀ at the NS surface as an additional model parameter. We found that the fits to the data from both the HL and LL observations were unable to provide reasonable constraints to all model parameters. The situation did not improve by fixing the compmag normalization as different values of N_comp resulted in virtually equivalent χ². By performing fits with values of N_comp comprised between a few tens and a few thousands, we found that a tight relation exists between this parameter, the accretion rate ṁ, and the radius of the column R₀, i.e. and (as expected according to the parameter definitions; see Sect. 2). We measured the slopes b₁ = 1.12 ± 0.03 and b₂ = 0.44 ± 0.01 from the fits to the LL observation, b₁ = 0.94 ± 0.02 and b₂ = 0.49 ± 0.04 from the HL data. By assuming for the pencil-beam geometry the same value of R₀ measured for the fan-beam geometry, we can use the above relations and coefficients to estimate the mass accretion rate in the former case. We obtained ṁ ~ 1.6 and ṁ ~ 1.9 for the LL and HL observation, respectively. In the pencil-beam geometry, the spectral parameters related to the cyclotron feature and the continuum are consistent (to within the uncertainties) between the LL and HL observations. The electron temperature is slightly higher for both observations than that measured from the fits with the fan-beam geometry, but in both geometries a consistent decrease of kT_e from the LL to the HL observation is measured. From the fits, we found β₀< 0.1 and z_max ≲ 4 km for either the LL and HL data.

5.2. 4U 0115+63

In the case of 4U 0115+63, we first tried to fit the spectra extracted from Obs. I and II (see Table 1) using an absorbed compmag model in the fan-beam geometry, plus two absorption Gaussian lines at energies corresponding to the fundamental and the second harmonic of the CRSF. This gave χ²/d.o.f. = 710/422 and χ²/d.o.f. = 753/526 for Obs. I and II, respectively. In both cases, we noticed a systematic excess at energies ≲2 keV together with an absorption feature around ~35 keV (the latter were more significant in Obs. I). We thus first added a partial covering component to the previous spectral model (pcfabs in Xspec), leading to χ²/d.o.f. = 518/420 and χ²/d.o.f. = 598/524 in Obs. I and II, and then added to the model a further Gaussian line to take into account the presence of the third CRSF harmonic; we obtained χ²/d.o.f. = 444/418 for Obs. I and χ²/d.o.f. = 583/522 for Obs. II, respectuvely. In both cases, the Gaussian widths of the third CRSF could not be constrained during the fit and thus we fixed this parameter to σ = 2 keV in agreement with the measured width of the second CRSF harmonic. In Obs. I, we additionally noticed some residuals from the fit around ~45 keV where the fourth CRSF harmonic is usually detected (see e.g. F09). Adding a fourth Gaussian line to the spectral model used for Obs. I led to a further improvement in the fit with χ²/d.o.f. = 363/415.

No evidence of the Kα iron emission line was found in the spectra of Obs. I and II. We list the best determined parameters of the spectral fits described above in Table 2. The unfolded spectra and the residuals from the fits are shown in Fig. 2. We also show in Fig. 4 the separated contributions of the Comptonized cyclotron and bremsstrahlung emission to the total X-ray spectra of 4U 0115+63; the former component is about 60% and 70% in Obs. I and II, respectively.

Fig. 2 Same as Fig. 1 but for the two BeppoSAX observations of 4U 0115+63.

Concerning the applicability of the compmag model in the pencil-beam geometry to 4U 0115+63 (assuming a free-fall velocity for the material inside the accretion column), the same considerations of Cen X-3 apply (see Sect. 5.1). From the fit to the source spectrum in Obs. I, we determined b₁ = 1.06 ± 0.02 and b₂ = 0.52 ± 0.01 for the two relations and , respectively. The values measured from the fit to the source spectrum of Obs. II are instead b₁ = 0.93 ± 0.02 and b₂ = 0.48 ± 0.01. The corresponding values of the mass accretion rate for a column radius R₀ equal to the value obtained in the fan-beam geometry are ~0.76 and ~0.52. Interestingly, the accretion rate ratio is equal to the luminosity ratio of the two observations (see Table 2). At odds with the case of Cen X-3, the flow terminal velocity and column height determined for the two 4U 0115+63 observations remain roughly constant, with β₀ ~ 0.2 and z_max ~ 3 km.

5.3. Her X-1

A fit with an absorbed compmag model, together with a cyclotron absorption line, to the combined Suzaku and NuStar observation of Her X-1 provided an unacceptable result (χ²/d.o.f. = 7387/3732) due to a strong emission emission line around 6.4–6.7 keV. The addition of an emission line at these energies led instead to a formally acceptable fit (χ²/d.o.f. = 3909/3729; see Table 2). It is worth noting that F13 modelled similar residuals by using two narrow lines with centroid energies of 6.4 and 6.5 keV, respectively (the estimated line widths were σ_G ~ 0.1 keV). One of the two lines was included mainly to smear out the sharp drop of the adopted phenomenological cut-off power-law model, and thus the second line is likely to have a model-dependent rather than a physical origin.

Fig. 3 Same as Fig. 1 but for the combined Suzaku and NuStar observation of Her X-1.

Fig. 4 From top to bottom: comptonized bremsstrahlung and cyclotron spectra of Cen X-3, 4U 0115+63, and Her X-1 obtained from the best-fit models given in Table 2.

Even though the fit was formally acceptable, some residuals were still visible around ~10 keV. These residuals were also observed by F13, who argued that they might have been caused by some instrumental effect and excluded the energy interval 10–12 keV from the spectral analysis. The fact that the same residuals emerged when using a completely different spectral model convinced us that they could have instead a physical origin. We found that the addition of an iron absorption edge at 9.1 keV correctly accounted for the residuals around 10 keV and slightly improved the fit (χ²/d.o.f. = 3769/3727). The estimated contributions of the cyclotron and bremsstrahlung Comptonized component to the total X-ray emission of the source in the 3–100 keV energy range were 60% and 40%, respectively. We report all values of the best-fit model to the data of Her X-1 in Table 2, and show the unfolded spectrum with the residuals from the fit in Fig. 3. The high statistical quality of the Her X-1 observations allowed all the parameters of the compmag model to be constrained reasonably well and provided the most accurate test of the model available so far.

Concerning the details of the CRSF in Her X-1, it is noteworthy that the width of this feature as measured by F13 was found to be dependent on the specific continuum model adopted by the authors and varied between 6 and 8 keV. The width of the CRSF estimated through the use of our compmag model is compatible with the value measured by F13 when a high-energy cut-off (highecut in Xspec) is used to characterize the data. This is also found by F13 to be the best spectral model, statistically speaking. The centroid energy of the CRSF measured by F13 and in the current paper are only marginally different (E_cyc = 37.7−39.1 keV and E_cyc = 37.2 ± 0.2 keV, respectively).

Finally, in the case of Her X-1 it was not possible to find a formally acceptable fit to the data by using the compmag model in the pencil-beam geometry and considering a free-fall velocity profile for the material within the accretion column.

6. Discussion

In this paper, we introduced an updated version of the compmag model, which permits us to compute the X-ray luminosity emerging from a highly magnetized accreting pulsar in the fan-beam and pencil-beam geometry.

Two particularly critical parameters of the model are the normalization N_comp and the mass accretion rate ṁ. In a purely Newtonian case, the mass accretion rate entirely regulates the source X-ray luminosity and it is expected that the normalization of the compmag model can be fixed to a constant proportional to the source distance (see Sect. 5). In a more realistic case, the X-ray emission produced relatively close to the NS surface will be strongly affected by GR effects, and the normalization of the compmag model is no longer trivially connected to the source distance. As in this work the compmag model was tested against observational data in which GR effects cannot be disentangled from the local X-ray emission, we left N_comp free to vary in the fits and determined its value for the different sources from the data. A better definition of N_comp would require the extension of the compmag model to include a full GR treatment of the X-ray emission from the pulsar, but this is outside the scope of the present paper.

The simplified approach adopted here could be successfully applied to the currently best available broad-band X-ray observations of three bright X-ray pulsars: 4U 0115+63, Cen X-3, and Her X-1. We used publicly available BeppoSAX archival data for Cen X-3 and 4U 0115+63, and a recent combined Suzaku +NuStar observation of Her X-1. In all cases, the broad-band coverage of the data ensured a detailed characterization of the source spectrum and provided the possibility of performing tests of the compmag model. In the case of 4U 0115+63, the fit performed to the two BeppoSAX observations with N_comp and ṁ left free to vary provided fully self-consistent physical results. We obtained a ratio of ~1.3 between the mass accretion rate in Obs. I and Obs II, and a ratio of ~1.8 between the corresponding values of N_comp. For comparison, the ratio between the X-ray luminosities recorded during these observations is ~1.5 (0.1–100 keV energy range; see Table 2). On the one hand, the difference in the mass accretion rate evaluated though compmag thus fully reflects the variation in the source luminosity across the two data sets. On the other hand, the difference in the pulse profile of the source that was reported previously among these two observations (Santangelo et al. 1999, F09) also justifies the measured variation in N_comp as some change has certainly occurred in the pattern of the X-ray radiation arising from the pulsar accretion column. A similar conclusion applies to the case of Cen X-3, for which a ratio of ~7.3 was measured from the fits with compmag between the mass accretion rats in the LL and HL observations and the correspondingly known ratio of the X-ray luminosity is ~8 (see Table 2). We note that in this case the value of N_comp of the LL observation had to be fixed to that of the HL observation in order to avoid degeneracies between this parameter and ṁ in the spectral fits (see Sect. 5.1). A similar check was not possible on Her X-1, as in this case only a single observation was available. However, the compmag model could also provide for this source a reasonably good fit to the broad-band X-ray data, giving indications about the structure of the accretion column (see Table 2). Generally, we note that a comparison between the values of ṁ determined from compmag for different systems should be carried out with caution, as uncertainties on the distances to these sources and intrinsic differences in the accretion geometry could lead to systematic biases. Values of ṁ are more easily comparable between different observations of a single source, as in this case the uncertainty on the distance cancels out and variations in these parameters can be directly linked to the physical processes intervening in the release of the X-ray emission.

There are two major limitations in the current version of the compmag model. The first is related to the crude approximation of the scattering cross section, which is described by parallel and perpendicular components to the magnetic field direction (σ_∥ and σ_⊥, respectively). The former appears in the Fokker-Planck treatment of the RTE Eq. (5) through the vertical spatial photon diffusion term, while the latter is related to the photon escape time across the side walls of the accretion column. The cross section , related to the photon diffusion in the energy space and representing the average scattering direction of photons parallel and perpendicular to the magnetic field, is also crudely approximated. These issues were already pointed out by BW07² and described in the first version of the compmag model by Farinelli et al. (2012a). We note them here once more for completeness.

As also briefly mentioned above, the second important limitation of the compmag model is the usage of a full Newtonian treatment to compute the pulsar X-ray emission. Because the emission is emerging close to the NS surface, it is likely that GR will play a key role in shaping the final energy distribution and intensity of this emission before it gets to the observer. The different geometrical parameters inferred from the compmag model, such as the accretion column radius R₀ and its height H, would probably need to be scaled in a non-trivial way to their values in the NS reference frame. This is, however, a limitation that also applies to several spectral models available within Xspec and commonly used in the X-ray astrophysical community. A well-known case is that of the diskbb model (Mitsuda et al. 1984), which assumes a Newtonian geometry and is widely used to characterize the soft spectral component of X-ray binaries hosting NSs and black holes. It is noteworthy that the accretion column radius derived from the application of the compmag model to the cases of 4U 0115+63, Her X-1, and Cen X-3, provides a value of R₀ ~ 1 km that is qualitatively in agreement with that expected when the accreting material is led by a relatively intense magnetic field before arriving close to the NS surface (B ~ 10¹² G; see e.g. Basko & Sunyaev 1976).

Finally, we discuss the use of the RTE (Eqs. (5) and (9)). In the present case, the RTE is written in the static reference frame of the observer with an accuracy of order β, the adimensional bulk velocity of the accreting material. It is important to outline that for a decelerated velocity profile as in Eq. (18), the electron density progressively increases towards the NS surface. As the combined continuum bremsstrahlung plus narrow cyclotron emission is proportional to (see Eqs. (A.1) and (A.2)), most of the seed photons are produced in the region where β → 0. Corrections between the reference frames of the fluid and the observer are thus minimal. The same consideration holds for the escape time of photons and the Comptonization efficiency: both processes are more efficient in a high electron-density environment. This also significantly reduces the uncertainty in Eq. (5) introduced by the use of the relation K/J = 1 / 3 between the second and zeroth moment of the intensity field, which would strictly hold only in the fluid reference frame (see Appendix B).

6.1. Comparison with earlier models

As already reported in Sect. 2 and in greater detail in F12, there are a number of differences between the current implementation of the compmag model and the former analytical treatment presented by BW07. So far, the only application³ of the BW07 model to a broad-band spectrum of an X-ray pulsar was presented by F09 using BeppoSAX data of 4U 0115+63. We thus compare below the results of the present paper for the source 4U 0115+63 with those reported in F09.

One important difference between the two treatments is that the BW07 model was proved unable to provide a satisfactory fit over the entire energy range spanned by the X-ray spectrum of 4U 0115+63. Actually, F09 showed that an additional thermal Comptonization component was needed below ~9 keV to correctly describe the data, while the higher energy part could be interpreted well in terms of Comptonized cyclotron emission. A broad Gaussian line also had to be included in the spectral model to achieve an acceptable result, but the origin of this feature could not be easily explained. They also showed that the magnetic field strength derived from the centroid energy of the cyclotron line was significantly different from that inferred from the fit to the continuum, and the values of ṁ and R₀ had to be fixed in the fit (6 × 10¹⁶ g s^-1 and 600 m, respectively) due to unresolved degeneracies of the different model parameters. As we showed in the previous sections, the compmag model is able to consistently fit the broad-band spectra of the different sources without the need of including additional components (in addition to the cyclotron absorption features).

It should be noted that the lower electron temperature of the bremsstrahlung seed component measured in the compmag and BW07 models (kT_e ~ 8 keV in F09 and kT_e ~ 0.7 keV in the present work) is known to arise from the different approximations of the Compton scattering: BW07 considered only the first-order bulk Comptonization term in the energy-diffusion operator, while we also introduced the second-order term. This component, being space-dependent, cannot be included when performing the analytical treatment of Comptonization through the variable-separation method (as adopted by BW07). It is a general feature of the thermal plus bulk Comptonization models that the inclusion of higher order Compton bulk terms lead to lower electron temperatures for the same spectrum (Titarchuk et al. 1997; Psaltis & Lamb 1997).

In our analysis of the X-ray emission from 4U 0115+63, we did not find any statistically significant evidence of the fifth harmonic of the cyclotron line, which was instead reported by F09. We note, however, that the parameters of the cyclotron lines are known to be strongly affected by the choice of the X-ray continuum and thus the detection of particularly weak wiggles in the high-energy end of an exponentially decaying spectrum can easily turn out to be statistically non-significant in different models.

7. Source term emission

In this section we justify the assumption of a combined thermal bremsstrahlung and cyclotron emission (the latter approximated through a Gaussian function) to generate the seed photons used for the Comptonization in the compmag model. We mentioned in Sect. 1 that this approximation was originally proposed by BW07, given the current impossibility of providing a self-consistent analytical description of the magnetic bremsstrahlung emission produced by the electrons in the accretion flow. While implementing the BW07 treatment of the problem in an Xspec model, F09 adopted an exponentially attenuated δ-function for the cyclotron term of the seed photons. These authors thus assumed that the width of the cyclotron line is intrinsically narrow at the origin and then broadens owing to thermal and bulk Comptonization effects. In the treatment presented in this paper, we used a Gaussian profile instead of a δ-function in order to take into account the natural width of the cyclotron emission feature (the compmag model parameter describing this effect is σ_cyc in Eq. (10)), which is expected to result from both thermal broadening and magnetic field mixing instabilities.

7.1. Bremsstrahlung emission

The fits to the data of the X-ray pulsars considered in this paper with the compmag model provided an electron temperature of the order of a few keV (see Table 2). It is straightforward to verify for each source that the fraction of electrons with velocities is very small (less than 1%) in both the cases of 1D and 3D Maxwellian distributions. The effect of tail depopulation because of the resonant electron-proton interaction with photon emission and the consequent deviation of f(v) from a Maxwellian distribution is thus marginal. Moreover, as previously mentioned in Sect. 2, the angular integration of the bremsstrahlung emission smears out narrow features arising from the effect of the viewing angle, resulting in a spectrum which fairly well agrees with the 3D classical case. We thus concluded that splitting the bremsstrahlung emission process in a strong magnetic field into a smooth continuum similar to the unmagnetized case plus a narrow cyclotron feature is a reasonably good assumption.

It should be noted that, in addition to the proton-electron bremsstrahlung, there is a second channel for the production of cyclotron emission photons due to the photon-electron resonant scattering (RS; see e.g. Canuto et al. 1971; Herold 1979; Daugherty & Harding 1986; Harding & Daugherty 1991). In principle, for a thermal plasma with a temperature of a few keV, only a small fraction of the electrons fulfil the condition E_el>E_cyc. The presence of a significant cyclotron feature in the X-ray data (see Fig. 4) would thus support the idea that the e-γ RS dominates over the e-p RS. In reality, a quantitative estimate of the relative contribution provided by the two processes is very difficult. The bremsstrahlung spectra reported by R99 were computed in the case of an optically thin plasma where the e-p processes dominate over the e-γ interactions. A more realistic solution to the problem would require a self-consistent treatment where both processes are taken into account simultaneously via a non-linear iterative solution scheme. The possibility of exchanging the bremsstrahlung and cyclotron terms for the seed photons in Eq. (5) indicates our inability to develop such a scheme.

7.2. Cyclotron emission

We now consider the broadening of the cyclotron emission feature due to magneto-hydrodynamical effects. First, we emphasize that the cyclotron emission term in the compmag model (Eq. (10)) is self-consistently weighted over the vertical profile of the accretion column because we assume B_z = B₀(z₀/z)³. An intrinsic broadening of the cyclotron emission line energy E_cyc is thus possible owing to geometrical effects (i.e. the distance from the NS) and the assumption on the vertical emissivity profile of the accretion column. The vertical gradient of the magnetic field does not provide, however, a dominant contribution to the broadening of the cyclotron emission feature as the total height of the column where the overall emission process takes place is relatively small (see Table 2) and cannot lead to a value of σ_cyc appreciably larger than zero. A similar conclusion holds for the radial gradient of the magnetic field. The exact expression of a dipolar magnetic field in the vertical direction using cartesian coordinates is (29)where μ is the magnetic moment, and R² = x² + y². The assumption B_z ∝ z^-3 considered so far is thus valid only in the asymptotic limit R ≪ z. However, it is easy to demonstrate that the variations of both B_z and | B | from the centre (R = 0) to the walls of the accretion column (R ≲ 1 km) are negligible and cannot contribute significantly to σ_cyc> 0.

Magneto-hydrodynamical simulations of hyper-critical accretion regimes onto a magnetized star have shown that significant fluctuations of the magnetic field intensity in the vertical direction can be produced as a consequence of field submerging effects by the copious inflowing material in the accretion column (Bernal et al. 2010). Numerical simulations of plasma accreted onto the NS surface also show that interchange instabilities contribute to increasing the gradient of the magnetic field intensity especially at the base of the accretion column (Mukherjee et al. 2013a,b). We thus argue that all these effects can easily give rise to the complex magnetic field structures that are needed to produce the observed broad cyclotron emission feature close to the NS surface.

It is also interesting to discuss the relative contribution of the Comptonized cyclotron and bremsstrahlung components to the total emission of the analysed sources, even though the two terms for the production of the seed photons for the Comptonization should be mutually linked in a more realistic treatment (see Sect. 2 and the discussion earlier in this section). As shown by Neugebauer et al. (1996) through detailed calculations, if the energy of an electron moving along a magnetic field line is higher than the cyclotron energy, the proton-electron Coulomb cross section has four resonances. The first two occur at v′ = ± v, where v and v′ are the velocities of the electron before and after the scattering, respectively. In these cases the p-e scattering process does not lead to the emission of a photon and the forward scattering probability is much higher than the backward owing to the high proton mass (making this a nearly pure elastic scattering). The other two resonances occur at ± v′, such that E′ = (E_el − E_cyc). Here E = 1 / 2m_ev², E′ = 1 / 2m_ev^′2, and the ± sign refers to the forward and backward scattering of the electron, respectively. In this case, the parallel velocity component (and thus the energy) is transferred to the quantized perpendicular electron Landau level, from which the electron decays releasing a cyclotron photon within a characteristic time (30)where B_crit = 4.41 × 10¹³ G. The p-e scattering leading to the emission of a photon is the most interesting process from an observational point of view. In particular, the resonance occurs because the electron in the intermediate state (before the decay) is generally in the excited Landau levels with n = 1. A p-e scattering leading to a photon emission and an electron in the excited Landau level with n = 0 during the intermediate state is also a possibility, but in this case the photons released after the decay would have a lower energy than the cyclotron value. This case is usually considered to be similar to a zero-magnetic-field bremsstrahlung emission. From the represented Comptonized cyclotron and bremsstrahlung emissions in Fig. 4, we conclude that for all sources considered in this work the bulk of the bremsstrahlung emission is produced at energies lower than the cyclotron emission, in line with the above discussion. The contributions of the two components to the total emission of each source are similar, as is expected in the case of p-e Coulomb scatterings where the electrons are populating (for the cyclotron component) or not (for the bremsstrahlung component) the Landau levels n ≥ 1 during their intermediate excited state.

It is also worth noting the strong correlation between the electron temperature and the magnetic field intensity derived from the fits with the compmag model to the data of all sources (Fig. 5). The Spearman correlation coefficient of the data set is R = 0.9, and the relation between the two parameters can be analytically described by the function , with a = 0.81 ± 0.10, b = 0.01 ± 0.008, and c = 4.61 ± 0.70. As the electron temperature in the compmag model determines the exponential roll-over of the spectra, the correlation we found is most likely related to the scaling law between the cut-off energy and the centroid energy of the cyclotron emission lines that was discovered thanks to Ginga and RXTE data through fits with phenomenological models (see e.g. Makishima et al. 1999; Coburn et al. 2002). The advantage of the compmag model in this case is that it can provide a physical interpretation of this correlation. In particular, we note that the electron temperature is about an order of magnitude lower than the cyclotron energy, which provides the maximum temperature at which electrons can be heated. As discussed earlier in this section, this is due to the collisional excitation of the Landaus levels, which limits the electron Maxwellian distribution below the cyclotron energy. As reported by Arons et al. (1987), the cyclotron emission is an effective cooling channel for the material in the accretion column (together with the Compton cooling). This is also the reason why there is an important contribution of the cyclotron emission term in the spectra of all X-ray pulsars analysed in this paper. Arons et al. (1987) also showed that the electron temperature is expected to be generally lower than the average photon energy, because photons have the highest heat capacity. This can be easily verified for 4U 0115+63, Her X-1, and Cen X-3 by using the results listed in Table 2.

While the relation between the magnetic field intensity and the electron temperature is thus well understood in terms of the physical processes regulating the accretion in highly magnetized X-ray pulsars, the possible correlations we identified in Sect. 5 between the different geometrical parameters of the compmag model (i.e. R₀ and z_max) should instead be taken with caution. These relations were derived by assuming a purely Newtonian treatment of the X-ray emission process close to the NS surface, where GR effects (e.g. light bending and lensing) can significantly modify the values of geometrical quantities between the emission and the observer frame. In order to be properly understood and interpreted, the same correlations should be investigated by introducing a ray-tracing code providing solutions for the equations describing the geodesics of the light emerging from the accretion column. This is, however, outside the scope of the present paper.

Fig. 5 Electron temperature as a function of the magnetic field at the base of the column for the sample of accreting X-ray pulsars analysed in this paper. The data values are given in Table 2 and can be fit by a power law (see text for more details).

8. Geometry of the emission region

In modelling the source spectra, we have considered both a pencil-beam geometry with an accelerating velocity profile in the accretion column and a fan-beam geometry where matter decelerates towards the NS surface. The first case is expected to occur in less luminous sources (L_X< 10³⁷ erg s^-1) where the radiation pressure generated by the soft photons from the NS surface is too low to halt the infalling matter in the accretion column and the radiation is preferentially emitted from the top of the column (fan-beam). At higher luminosities (L_X> 10³⁷ erg s^-1), the strong radiation pressure from the NS surface contrasts the pressure of bulk inflow and leads to the formation of a radiative shock above the NS surface. In this case, the energy dependency of the electron-photon scattering cross section leads to substantially different behaviours below and above the cyclotron line energy. Below the cyclotron energy, the optical depth perpendicular to the magnetic field is higher than in the parallel direction. The Comptonized radiation thus emerges preferentially through the lateral walls of the accretion column (fan-beam), while part of the weakly reprocessed radiation can still be emitted along the column. Above the cyclotron energy, the cross sections in both directions become comparable, increasing the emission through the lateral walls of the column. The compmag model is the best choice to fit phase-averaged spectra as in this case the effects described above are averaged out and the treatment of the cross sections in the RTE approximation is justified.

All the sources considered in this paper have X-ray luminosities L_X> 10³⁷ erg s^-1, and thus the usage of the fan-beam geometry in the compmag model was considered more appropriate to fit their X-ray spectra. We have shown in principle that fits with the compmag model and the pencil-beam geometry are possible for Cen X-3 and 4U 0115+63, while no solution could be found in the case of Her X-1. In the pencil-beam case, the empirical dependence discussed in Sects. 5.1 and 5.2 is derived from the definition of the compmag normalization for the top-column emission (see Eq. (22)). Indeed, the N_comp parameter explicitly depends on the emitting area at the top of the accretion column, and on the magnification factor f_geom, which parametrizes all the unknown geometrical and GR effects that are not taken into account in the present version of the compmag model. The inverse proportionality shows that the compmag model correctly takes into account a linear dependence between the source flux and the estimated mass accretion rate. For a given value of the source X-ray flux, a decrease of N_comp has to be compensated by a corresponding linear increase in ṁ. We verified that the changes in ṁ and R₀ derived as a consequence of a variation in the model normalization, N_comp, are sufficient to keep the value of χ²unchanged. As we reported during the analysis of Cen X-3 and 4U 0115+63, both β₀ and z_max are insensitive to variations of N_comp. This can be explained considering that the emerging flux (and thus the source X-ray luminosity) in the Eddington approximation (see Eq. (24)) is mostly regulated by the vertical space gradient of the photon field at the boundary of the accretion column and not by the height of the column z_max. Although the value of ∂n/∂τ at the top of the accretion column is closely related to the overall space photon distribution, the effect of the value of z_max is much smaller in the pencil-beam case than in the fan-beam case. In the latter, this parameter explicitly appears in the integration of the source flux performed along the column height (see Eq. (26)). In this case, z_max has more influence on the vertical gradient of the magnetic field and thus also on the cyclotron emission term.

On the one hand, it is thus evident that we should generally expect some degeneracy in the compmag model parameters and it might not always be possible to verify a priori if the fan-beam or the pencil-beam geometry is better suited to describing the X-ray spectrum of an X-ray pulsar. On the other hand, the fact that the pencil-beam geometry could not be used to satisfactorily fit the X-ray spectrum of Her X-1 gives us confidence that the compmag model might be able to distinguish between the different emission states of these sources, at least when high-quality and wide broad-band data are available. Additional observations of these sources with the next generation of X-ray instruments providing broad-band X-ray coverage (like those on-board Astro-H; see e.g. Takahashi et al. 2010) will give us the possibility of performing additional tests of the compmag model and improving the treatment of all physical processes currently included in the corresponding code.

9. Conclusions

In this work we presented an updated version of the compmag model and applied it to the case of three bright X-ray pulsars for which broad-band high-quality data were available. We found that the model is well suited to describing quantitatively the X-ray spectral energy distribution of these systems and that it provides reasonable values of the physical and geometrical parameters that are used to compute the high-energy emission of the different sources. We discussed in detail the limitations and the difficulties of the computations performed within the model. They are mainly related to the use of vertically and horizontally averaged photon-electron scattering cross sections, together with a cylindrical approximation treatment in which all quantities vary only along the vertical coordinate. For these reasons, the model can be better exploited to carry out the analysis of phase-averaged spectra. Despite its limitations, the compmag model represents an additional significant step forward in the study of highly magnetized accreting X-ray pulsars.

Acknowledgments

The authors thank the anonymous referee who provided many suggestions that significantly improved the first version of the paper. They are also grateful to the ASI data centre for preserving and making available the BeppoSAX data, and to F. Fürst for providing the NuStar and Suzaku spectra of Her X-1 used in this work. This research has made use of data obtained with the NuSTAR mission, a project led by the California Institute of Technology (Caltech), managed by the Jet Propulsion Laboratory (JPL) and funded by NASA, of observations obtained with the Suzaku satellite, a collaborative mission between the space agencies of Japan (JAXA) and the USA (NASA), and with the X-ray astronomy satellite BeppoSAX, a project of the Italian Space Agency (ASI) with participation of the Netherlands Agency for Aerospace Programs (NIVR).

References

Appendix A: Derivation of the bremsstrahlung and cyclotron source terms in the RTE

The cyclotron and bremsstrahlung emissivity in units of erg cm^-3 s^-1 Hz^-1 are (A.1)where and (A.2)with the temperature T_e expressed in Kelvin.

The units of the RTE Eq. (5) are s^-1, while the occupation number is related to the specific intensity by (A.3)where J(E) is given in units of energy cm^-2 s^-1 Hz^-1 ster^-1. We note that in Eq. (A.3) both J(E) and n(E) are averaged over the angle. To obtain the correct units for Eq. (5), the quantity J(E) must be multiplied by the extra factor c.

Using the definitions E = hν and x = E/kT_e, together with the Eq. (A.3), the cyclotron and bremsstrahlung source terms in the forms used within Eq. (5) can be easily obtained as (A.4)We note, however, that the compmag model solves Eq. (9), which is obtained by dividing Eq. (5) by n_eσ_∥cH. By additionally defining Θ = kT_e/m_ec², we finally obtain the expression (A.5)which provides the source term defined in Eqs. (10) and (12) once the proper values of the constants are used. We also note that in order to derive the expression for the cyclotron term, we first used the property of the Dirac δ-function (A.6)and then substituted δ(x − x_cyc) with a normalized Gaussian.

Appendix B: Moments of the intensity radiation field in the observer reference frame

Let us consider a fluid moving with bulk velocity V with respect to a static system of reference, i.e. the observer reference frame. In the diffusion approximation, the specific intensity in the fluid reference frame can be written as (e.g. Pomraning 1973) (B.1)where k_f is the cosine-direction versor, while I_0,f and M_f are the zeroth and first moments of the specific intensity, respectively. The latter are defined as (B.2)and (B.3)respectively. We note that is the spectral flux along the ith_f-coordinate. Together with the zeroth and first moment, the second moment of the radiation field can also be defined as (B.4)For the sake of clarity, we consider here the case of a cartesian coordinate system where only the M^z component is different from zero (slab approximation) so that Eq. (B.1) can be written as (B.5)where is the cosine of the photon direction along the z coordinate. Substituting I_f from Eq. (B.5) into the integrals (B.2) to (B.4), it is straightforward to see that K_f/I_0,f = 1 / 3, where .

We now consider the relation between the zeroth and second moment in the observer reference frame. The transformation law of the photon energy between the fluid and observer reference frame is (B.6)where β = V/c is the fluid velocity, and the relation between the angles u_f and u_l is given by the aberration formula (B.7)Using now the relativistic invariant (B.8)together with the relation (B.6), the specific intensity in the laboratory reference frame can be written as (B.9)where I_f is given in Eq. (B.1). Computing the zeroth, first, and second moment of the specific intensity in the observer reference frame we obtain (B.10)(B.11)and (B.12)We also made use above of Eqs. (B.9) and (B.7). From Eqs. (B.10) and (B.12) we finally obtain (B.13)

All Tables

All Figures

	Fig. 1 Unfolded spectra of the two BeppoSAX observations of Cen X-3. Also show are the best-fit compmag model and the residuals from the fit in units of σ (bottom panel). All best-fit parameters determined from the fit are given in Table 2.
In the text

	Fig. 2 Same as Fig. 1 but for the two BeppoSAX observations of 4U 0115+63.
In the text

	Fig. 3 Same as Fig. 1 but for the combined Suzaku and NuStar observation of Her X-1.
In the text

	Fig. 4 From top to bottom: comptonized bremsstrahlung and cyclotron spectra of Cen X-3, 4U 0115+63, and Her X-1 obtained from the best-fit models given in Table 2.
In the text

	Fig. 5 Electron temperature as a function of the magnetic field at the base of the column for the sample of accreting X-ray pulsars analysed in this paper. The data values are given in Table 2 and can be fit by a power law (see text for more details).
In the text